Известия РАН. Серия физическая, 2021, T. 85, № 6, стр. 875-882

ВВЕДЕНИЕ

Последние два десятилетия наблюдается быстрое развитие методов неинвазивной ультразвуковой хирургии с помощью высокоинтенсивного сфокусированного ультразвука (HIFU) [1]. При этом терапевтический эффект осуществляется либо за счет нагрева ткани и ее тепловой денатурации, либо посредством механического разрушения (“гистотрипсии”) ткани на фрагменты субклеточных размеров [2–4]. Область применения HIFU методов гистотрипсии очень обширна и включает, например, лечение опухолей печени и почек, а также ликвификацию абдоминальных гематом [5, 6].

Общий подход к осуществлению HIFU состоит в фокусировке ультразвукового пучка в заданную область тела через согласующую среду и лежащие на пути ткани. К преимуществам метода относятся неинвазивность, отсутствие требования стерильной операционной, уменьшение риска заражения во время операции, возможное ускоренное восстановление пациента после манипуляций [1, 5]. В то же время, одним из существенных ограничений метода HIFU в настоящее время является отсутствие надежных методов предварительного планирования оперативного вмешательства для каждого конкретного случая. Поэтому одной из важных задач в неинвазивной ультразвуковой хирургии является разработка методов планирования безопасных и эффективных протоколов терапевтического воздействия.

Для решения этой задачи необходимо иметь количественную информацию о параметрах ультразвукового поля в теле человека. Экспериментальное определение таких параметров в клинических условиях не представляется возможным. Альтернативным способом является использование методов численного эксперимента [7–9]. Задача о прохождении пучка в теле человека осложняется тем, что различие в акустических свойствах разных биологических тканей и органов делает среду распространения неоднородной, что может приводить к искажению структуры фокального пятна, а также к смещению положения фокуса за счет эффектов рефракции [10, 11]. Отдельной задачей является определение акустических свойств различных тканей пациента тем или иным методом.

В наших предыдущих работах рассматривалась задача облучения почки человека [12]. Была промоделирована фокусировка ультразвукового пучка с использованием акустической модели, созданной на основании данных компьютерной томографии. Результаты моделирования показали, что за счет рефракции на неоднородностях ткани фокальный максимум смещается, и его форма искажается относительно фокального максимума в однородной среде. При этом эффективность и точность терапевтического воздействия снижается. В связи с этим возникает задача разработки методов компенсации или уменьшения искажений, вносимых неоднородной средой.

В работе решалась задача по исследованию влияния неоднородностей брюшной стенки на фокусировку ультразвука в зависимости от положения и ориентации излучателя относительно тела человека [13, 14]. Задача рассматривалась на примере использования HIFU облучения для ликвификации абдоминальных гематом (рис. 1). В клинических условиях лечение предполагается проводить следующим способом: мощный ультразвуковой пучок фокусируется в некоторую область гематомы за брюшной стенкой (рис. 1). Под его воздействием в режиме гистотрипсии происходит разжижение гематомы. Для отведения образовавшихся продуктов разрушения в область облучения вставляется тонкая игла и производится аспирация содержимого [6]. Мы рассматривали случай, когда гематома располагается за брюшной стенкой недалеко (5–6 см) от поверхности кожи. На первом этапе исследования задача решалась в линейном приближении.

АКУСТИЧЕСКАЯ МОДЕЛЬ ТЕЛА ЧЕЛОВЕКА

Важным шагом в работе является создание и подготовка реалистичной модели акустически неоднородного торса человека с гематомой, так как для построения численного решения соответствующих модельных уравнений в некоторой области необходимо знание всех акустических параметров ткани как функций координат. Такими параметрами в случае линейной модели распространения волны служат скорость звука, плотность и коэффициент поглощения. Биологическая среда является неоднородной по отношению к данным характеристикам. Внутри мягких тканей скорость звука и плотность связаны с относительными пропорциями воды, белка (коллаген и гемоглобин) и липидов, и, следовательно, по своей природе неоднородны [10, 11, 14]. Высокое содержание белка в ткани увеличивает скорость звука и плотность по сравнению с этими же параметрами в воде или крови, в то время как при высоком содержании липидов они становятся ниже, например, в жире.

Пространственные распределения акустических параметров были получены из данных компьютерной томографии торса тела человека. Компьютерная томография (КТ) является одним из наиболее успешных методов неразрушающего послойного исследования внутренней структуры тела человека. Метод основан на измерении и компьютерной обработке разности ослабления рентгеновского излучения различными по плотности и составу тканями [15]. Из КТ изображений, яркость которых описывается по шкале Хаунсфилда, можно восстановить плотность тканей в каждой точке тела пациента [15].

Анализ плотности и скорости звука различных тканей показал, что эти два параметра сильно коррелируют друг с другом и между ними можно установить линейную зависимость [16]. Пользуясь установленными соотношениями, можно восстановить распределение скорости звука в каждой точке акустической модели. С другой стороны, было выявлено, что коэффициент поглощения плохо коррелирует с плотностью, поэтому восстановление этого параметра приходится проводить другим методом. С этой целью выполнялась сегментация акустической модели, то есть выделение различных типов тканей и органов с присвоением им определенного номера. Из изображений КТ, полученных на основе выбора пороговых значений скорости звука, были выделены по типам ткани: кожа, жир, мышцы, гематома, кость [12]. При построении трехмерных матриц акустических параметров отсегментированным областям были присвоены известные из литературы значения коэффициента поглощения [17]. Далее работа велась только с отсегментированными данными и двумя матрицами, содержащими распределения скорости звука и плотности в среде.

Данные акустической модели неоднородной среды были получены с шагами дискретизации изображений КТ в декартовой системе координат по направлениям x и y с шагами Δx = Δy = 0.77 мм (поперек торса) и по направлению z с шагом Δz = = 2.5 мм (вдоль торса). Такой шаг матрицы неоднородностей является достаточно крупным и сравнимым с размерами фокальной области ультразвукового пучка. Поэтому, если непосредственно использовать такую матрицу в акустическом моделировании, то анализ изменения формы фокуса, образования побочных максимумов и других важных характеристик становится невозможным. Для обеспечения точности расчетов акустического поля было проведено уменьшение шагов дискретизации до значений Δx = Δy = Δz = = 0.5 мм с использованием метода трилинейной интерполяции, реализованного в среде MathLab. Далее распределения акустических параметров среды переносились на численную сетку в декартовой системе координат, связанной с излучателем, в которой ось z направлена вдоль оси ультразвукового пучка (рис. 1).

При определении наиболее выигрышного направления облучения, важно избегать попадания ультразвукового пучка на кости и воздухсодержащие органы пациента. Кости сильно поглощают и отражают ультразвук, что может привести к их перегреву и ожогам. Наличие пузырьков воздуха на пути ультразвукового пучка вызывает его отражение и рассеяние. Расположение излучателя таким образом, чтобы создаваемый им пучок проходил самый короткий путь через мягкие ткани до точки фокуса, представляется наиболее выигрышным решением. Это позволяет уменьшить эффекты поглощения в ткани и эффекты рефракции, приводящие к смещению фокуса. В качестве объекта терапевтического воздействия HIFU рассматривается гематома, располагающаяся за брюшной стенкой человека, которая является “простым” объектом облучения, потому что находится непосредственно около кожи, далеко от ребер и тазовых костей, а также далеко от воздушных полостей, образуемых в кишечнике и других газосодержащих органах. Основным параметром выбора оптимального расположения излучателя стал угол, при котором ультразвуковой пучок проникает в мягкие ткани практически перпендикулярно к поверхности кожи.

После выбора области и направления облучения 3D матрица акустических параметров, полученных из КТ изображения, поворачивалась так, чтобы ось z новой системы координат совпадала с осью ультразвукового пучка. С точки зрения математического моделирования осуществляется поворот не излучателя относительно модели, а торса человека относительно главной оси излучателя. Для этого необходимо имеющиеся 3D матрицы, содержащие в себе информацию о значениях скорости звука, плотности и коэффициента поглощения в каждом объеме 0.5 мм × 0.5 мм × × 0.5 мм среды облучения в имеющейся системе координат, перевести в новую систему координат. Как правило, для этого используют матрицу поворота – матрицу, которая используется для выполнения собственного ортогонального преобразования в евклидовом пространстве и содержащую в себе информацию о величине угла поворота и оси, относительно которой совершается вращение. Соответственно, любое вращение в трехмерном пространстве может быть представлено как композиция поворотов вокруг трех ортогональных осей. Композиция представляет собой матрицу, являющуюся результатом произведения соответствующих трех матриц поворота для каждой оси. Для удобства, матрицу вращения выражают через углы Эйлера – углы, описывающие поворот абсолютно твердого тела в трехмерном евклидовом пространстве. Это углы, которые соответствуют вращению относительно одной из осей, а повороты вокруг каждой из осей осуществляются последовательно. Сначала поворот осуществлялся относительно оси Oz на угол $\alpha ,$ который называют углом прецессии, затем – относительно “новой” оси Ox' (образованной после первого поворота) на угол $\beta $, называемый углом нутации, а последний поворот – относительно “новейшей” оси Oz' (образованной после первого и второго поворотов) на угол $\gamma ,$ называемый углом собственного вращения. Для каждого этапа вращения можно написать свою матрицу поворота:

Где результирующую матрицу поворота можно получить путем перемножения этих трех матриц:

Результирующая матрица (2) и была использована для пересчета координат из имеющейся системы в выбранную как “наиболее выигрышную”. Выбор угла основывался на использовании одного из нескольких наиболее вероятных предполагаемых расположений, где соблюдалась максимальная перпендикулярность пучка к коже. Результат поворота показан на рис. 2, где рис. 2а соответствует положению излучателя до поворота, а рис. 2б – после поворота. Белым цветом выделена гематома, а черным – излучатель. Все остальные типы биологических тканей и согласующая среда (вода) показаны в различных оттенках серого.

ТЕОРЕТИЧЕСКАЯ МОДЕЛЬ И ПАРАМЕТРЫ ФОКУСИРОВКИ

Моделирование распространения ультразвукового пучка с частотой f₀ = 1 МГц в неоднородной биологической ткани проводилось с использованием программного обеспечения k-Wave: свободно доступного инструментария для среды MathLab, предназначенного для описания акустических волновых полей [18]. Численные методы, заложенные в k-Wave, позволяют выполнить моделирование акустических волн в неоднородных средах со степенным законом поглощения [19]. В общем случае для расчета поля высокоинтенсивной ультразвуковой волны в теле человека, определяющие уравнения должны учитывать комбинированные эффекты дифракции, нелинейности, поглощения, дисперсии и неоднородности акустической среды, т.е. различие в акустических параметрах тканей и органов, таких как скорость звука, плотность, поглощение и параметр нелинейности. Соответствующие уравнения могут быть записаны в виде системы дифференциальных уравнений в частных производных первого порядка. В нашей работе используется линеаризованный вид данных уравнений:

Здесь $\vec {u}$ – вектор колебательной скорости в среде, $\vec {d}$ – вектор смещения частиц среды, p – акустическое давление, ρ – флуктуации плотности, ${{\rho }_{0}}(\vec {r}),$ ${{c}_{0}}(\vec {r})$ и $\delta (\vec {r})$ – равновесные плотность, скорость звука и диссипативный коэффициент неоднородной среды, соответственно, $\vec {r}$ – радиус-вектор. Диссипативный коэффициент $\delta $ полагался равным нулю в воде, а для сегментов неоднородной брюшной стенки использовалась однозначная связь $\delta = {{2\alpha {{c}_{0}}} \mathord{\left/ {\vphantom {{2\alpha {{c}_{0}}} {(2\pi {{f}_{0}})}}} \right. \kern-0em} {(2\pi {{f}_{0}})}}$ между диссипативным коэффициентом $\delta $ и коэффициентом поглощения на частоте f₀ = 1 МГц, полученным из литературы [17]. Значения коэффициентов поглощения, диапазоны скоростей звука и плотностей различных мягких тканей, использовавшихся в модели тела человека, представлены в табл. 1.

Систему (3) решали с использованием псевдоспектрального метода, где пространственные градиенты акустического поля рассчитываются в k‑пространстве на основе преобразования Фурье, а двухслойная схема интегрирования по времени включает поправку для компенсации численной дисперсии. При моделировании рассматривали излучатель с рабочей частотой 1 МГц, фокусным расстоянием 85 мм, диаметром центрального отверстия 30 мм и апертурой 85 мм. Граничное условие в виде начальной колебательной скорости на плоскости z = 0 рассчитывали с помощью интеграла Рэлея для шага сетки 0.5 мм в предположении равномерности амплитуды нормальной составляющей скорости колебаний поверхности излучателя. Фокусировка пучка осуществлялась в область гематомы на глубину 64 мм в положении до поворота (рис. 2а) и на 55 мм в положении после поворота (рис. 2б). Изменяя направление прохождения волны через границу воды и кожи, можно было определить оптимальное расположение излучателя, при котором поле акустического давления было бы по амплитуде и форме наиболее близким к амплитуде и форме поля при фокусировке в однородной среде (воде). Расчеты проводили в системе координат, связанной с излучателем. Пространственные распределения акустических параметров поворачивали на заданные углы вокруг точки геометрического фокуса излучателя (рис. 2), после чего определяли параметры среды в узлах численной сетки, связанной с излучателем. Результаты расчетов акустического поля для каждого угла сравнивали между собой.

РЕЗУЛЬТАТЫ

Результаты моделирования представлены на рис. 3, где показаны пространственные распределения интенсивности, нормированной на интенсивность на излучателе, вдоль осей z (рис. 3а), x (рис. 3б), y (рис. 3в), и на рис. 4, где показаны пространственные распределения амплитуды давления поля, нормированной на давление на излучателе, в осевых плоскостях излучателя xz и yz, и в фокальной плоскости xy для случаев различного расположения излучателя относительно тела, а также в однородной среде (воде). Для проверки правильности и точности численных расчетов показано также аналитическое решение на оси пучка, полученное с помощью интеграла Релея в воде (рис. 3). Данный график практически полностью совпадает с графиком решения для поля интенсивностей в воде, полученного методом k-Wave, что показывает корректность работы программы и позволяет использовать ее и для более сложных случаев, где аналитическое решение рассчитано быть не может.

На рис. 4а–4в показаны пространственные распределения амплитуды давления поля, нормированной на давление на излучателе, в однородной среде без поглощения (воде), которые в данной работе считаются “эталонными”. При положении излучателя как на рис. 2а (до поворота), акустическое поле в гематоме представлено на рис. 4г–4е, а при положении как на рис. 2б (после поворота) – на рис. 4ж–4и. Видно, что в присутствии неоднородностей происходит искажение фокального максимума – он изменяет форму и смещается из исходного положения на оси излучателя (рис. 4г–4е и 4ж–4и) в сравнении со случаем в воде (рис. 4а–4в). Когда облучение идет под наклоном к поверхности кожи (рис. 2а), фокальное пятно искажается из-за преломления лучей на границах слоев ткани, а также смещается от оси пучка по направлениям x и y на 0.5 мм (рис. 3б, 3в), а вдоль z – на 3.5 мм ближе к излучателю (рис. 3а). Кроме того, рефракция на неоднородностях ткани совместно с поглощением приводят к уменьшению амплитуды давления в фокусе, что затрудняет реализацию терапевтического воздействия ультразвука на ткань. В положении, когда ось пучка перпендикулярна касательной к поверхности кожи (рис. 2б) – фокальное пятно практически не меняет форму и смещается вдоль оси пучка z на 2 мм ближе к излучателю (рис. 3а). При этом величина смещения по оси y составляет порядка 0.5 мм, в то время как по оси x смещение фокуса практически отсутствует (рис. 3б, 3в).

При фокусировке под углом интенсивность в фокусе уменьшается на 83% (рис. 2а) относительно фокусировки в воде, в то время как при фокусировке перпендикулярно поверхности кожи (положение излучателя как на рис. 2б) – на 51%, что можно заметить на всех графиках рис. 3. Соответственно, после подбора более выигрышного положения интенсивность в фокусе выросла больше чем на 30% по сравнению с интенсивностью в фокусе при положении до поворота.

Была также проведена оценка вклада поглощения в общие потери интенсивности в фокусе. Для этого был промоделирован случай фокусировки в неоднородной среде с нулевым поглощением при положении излучателя в позиции до поворота (рис. 3). Было показано, что вклад в потери от поглощения составляет 63% от всех потерь по интенсивности в фокусе. Однако в других случаях неоднородности брюшной стенки или при более глубоком залегании фокуса вклады поглощения и рефракции в потери могут меняться. В будущих работах планируется изучение более сложных случаев расположения гематом в теле человека, однако уже на рассмотренном случае, который можно классифицировать как “простой”, подтверждается важность и эффективность предварительного планирования оперативного вмешательства с подбором оптимального пути облучения.

ЗАКЛЮЧЕНИЕ

На основе численного решения уравнений линейной акустики в неоднородной среде с частотно-зависимым поглощением с использованием пакета k-Wave была построена теоретическая модель для описания фокусировки ультразвуковых пучков в неоднородных мягких тканях тела человека. Проведено моделирование ультразвукового пучка, создаваемого терапевтическим излучателем, при фокусировке в область абдоминальной гематомы на глубину 55 и 64 мм в зависимости от положения излучателя. Показано, что уменьшение интенсивности в точке фокуса может достигать около 80% за счет поглощения ткани и эффектов аберраций пучка в неоднородной среде по сравнению с фокусировкой в воде. Примерно половину от этих потерь составляет непосредственно поглощение, которое можно компенсировать увеличением мощности на излучателе. Преломление пучка на поверхности тела и неоднородности скорости звука при распространении волны приводят к смещению точки фокуса в продольном и поперечном направлениях на расстояние до 1 мм, а также к искажению формы фокального пятна. Влияние данных эффектов зависит от расположения излучателя, что необходимо учитывать при разработке протоколов ультразвукового облучения.

Работа выполнена при поддержке РФФИ (проект № 20-02-00210а).