Содержание к диссертации
Введение
ГЛАВА 1. Терапевтические применения мощного фокусированного ультразвука 16
ГЛАВА 2. Особенности численного моделирования и измерения дифрагирующих акустических волн с ударными фронтами 27
2.1. Введение 27
2.2. Сравнение численных схем при моделировании плоских волн конечной амплитуды 31
2.2.1. Спектральный алгоритм {схема Рунге-Кутта) 34
2.2.2. Временной подход {схема Годунова) 38
2.3. Особенности моделирования нелинейных сильно фокусированных акустических волн с ударным фронтом 45
2.4. Влияние ограниченной полосы пропускания гидрофона на точность измерения параметров акустической волны в режиме развитых разрывов 51
2.4.1. Схема экспериментальной установки 52
2.4.2. Численные и экспериментальные результаты 53
2.5. Выводы главы 2 56
ГЛАВА 3. Фокусировка мощных ультразвуковых пучков в воде и предельные значения параметров разрывных волн 57
3.1. Описание нелинейных фокусированных акустических полей. Обзор литературы 57
3.2. Нелинейное изменение коэффициентов концентрации параметров поля и предельные значения параметров разрывных волн при фокусировке 67
3.3. Сравнение результатов численного моделирования и эксперимента 73
3.4. Пространственное распределение параметров акустического поля 75
3.5. Влияние аподизации поля на излучателе на проявление нелинейных эффектов 82
3.6. Выводы главы 3 86
ГЛАВА 4. Метод определения параметров акустического поля в биологической ткани 88
4.1. Существующие методы расчета акустического поля в биологической ткани. Обзор литературы 88
4.2. Метод определения параметров акустического поля в биологической ткани с учетом нелинейных эффектов 91
4.3. Численная и экспериментальная проверка метода 93
4.3.1. Результаты численного моделирования 93
4.3.2. Сравнение с экспериментальными данными 99
4.3.3. Возможность дирейтинга мощности тепловыделения 103
4.4. Выводы главы 4 107
ГЛАВА 5. Эффект локального сверхбыстрого кипения в ткани при воздействии фокусированным пучком разрывных волн 108
5.1. Введение 108
5.2. Теоретическая модель для описания процесса нагрева биологической ткани 111
5.2.1. Численная нелинейная модель 111
5.2.2. Аналитические оценки нагрева среды и времени до начала кипения 112
5.3. Численная и экспериментальная оценка времени до начала кипения в образце при облучении мощным фокусированным ультразвуком 116
5.3.1. Результаты для гелевого фантома ткани 116
5.3.2. Результаты для биологической ткани (ex vivo печени) 124
5.4. Выводы главы 5 126
Заключение 128
Благодарности 130
Список литературы 131
- Сравнение численных схем при моделировании плоских волн конечной амплитуды
- Влияние ограниченной полосы пропускания гидрофона на точность измерения параметров акустической волны в режиме развитых разрывов
- Пространственное распределение параметров акустического поля
- Метод определения параметров акустического поля в биологической ткани с учетом нелинейных эффектов
Введение к работе
Актуальность темы
Изучение проблем, связанных с распространением высокоинтенсивных фокусированных ультразвуковых пучков, является одним из активно развивающихся современных направлений нелинейной акустики. Возросший в последнее время интерес к этим задачам во многом обусловлен созданием медицинских приборов нового поколения, в которых фокусированный ультразвук высокой интенсивности (или сокращенно HIFU от словосочетания High Intensity Focused Ultrasound) используется для неинвазивного разрушения опухолей, остановки внутренних кровотечений при травмах, ультразвуковой коррекции фигуры, направленной доставки лекарств.
При использовании любого медицинского оборудования необходимо полное понимание всего спектра его возможного влияния на организм человека, как благоприятного, так и, в особенности, неблагоприятного. Разработка подходов для измерений и расчетов параметров волновых полей, создаваемых устройствами для ультразвуковой диагностики и терапии, является необходимым условием их использования. В области HIFU этот вопрос остаётся до сих пор не решенным. Несмотря на то, что HIFU системы уже используются в клиниках, до сих пор не существует общепринятых международных стандартов для описания создаваемых ими полей и их сертификации. Физические механизмы воздействия мощного ультразвука на ткань и соответствующие биоэффекты также требуют дальнейших исследований. Существует проблема выбора наиболее важных параметров акустического поля, которые определяют те или иные биологические эффекты в тканях. Одним из основных осложняющих факторов при описании полей HIFU является сильное проявление нелинейных эффектов, которые приводят к генерации высших гармоник в спектре распространяющейся волны, асимметричному искажению профиля волны, формированию ударных фронтов и дополнительному поглощению энергии волны на разрывах.
На сегодняшний день режимы работы систем ультразвуковой хирургии характеризуются значениями усреднённой по времени и/или по области фокального пятна интенсивности, рассчитанными в приближении линейного распространения гармонической волны. Временной профиль волны при этом остается гармоническим, и величина интенсивности однозначно определяет пиковые значения и пространственные распределения акустического давления. По известному распределению интенсивности линейной волны можно
рассчитать ожидаемый нагрев ткани, а максимальная величина отрицательной фазы давления позволяет оценить вероятность кавитационных эффектов. Однако ясно, что в нелинейном поле знания только интенсивности недостаточно. В нелинейном пучке тепловые эффекты усиливаются за счёт генерации высших гармоник, а при образовании разрывов в профиле волны мощность тепловыделения может возрасти в несколько десятков раз. Максимальное отрицательное давление в асимметрично искаженном нелинейном профиле волны, напротив, всегда меньше, чем в линейном пучке. По мере усиления нелинейных эффектов пространственная структура различных параметров поля, таких, как амплитуда высших гармоник, пиковые положительное и отрицательное давления, полная интенсивность и эффективная мощность тепловыделения волны, будет изменяться различным образом.
Информация о пространственной локализации параметров нелинейных полей, создаваемых медицинскими ультразвуковыми преобразователями, необходима для планирования и предсказания результатов ожидаемого терапевтического воздействия. Измерение нелинейных полей ультразвуковых преобразователей в широком диапазоне используемых мощностей даже в воде является сложной дорогостоящей процедурой, требующей использования миниатюрных широкополосных гидрофонов, и практически невозможно в биологической ткани. Детальное измерение поля в воде обычно проводится лишь при малых интенсивностях, а затем полученные результаты переносятся на большие интенсивности путем линейной экстраполяции. Ясно, что при этом теряются нелинейные эффекты, проявляющиеся различным образом при различных уровнях интенсивности и по-разному изменяющие различные параметры акустического поля. При характерных используемых в HIFU интенсивностях нелинейные эффекты могут даже привести к эффекту насыщения, когда параметры поля в фокусе за счет поглощения энергии волны на разрывах перестают зависеть от мощности излучателя.
Еще более сложной оказывается проблема перенесения полученных в воде результатов на биологическую ткань, что необходимо при планировании хирургической процедуры. В настоящее время это делается путем уменьшения полученного значения интенсивности в воде с учётом поглощения в ткани 0.05-Ю. 15 Нп см" МГц" . Такой упрощенный метод расчета приводит сразу к нескольким ошибкам, связанным с игнорированием или с неправильным учетом нелинейных эффектов, и потому, в любом случае, может дать сильно отличающийся от реальности результат. Дело в том, что даже если нелинейные
поля достаточно точно измерены в воде, эти данные невозможно прямым образом перенести на ткань, в которой поглощение гораздо больше. Соответственно, и накапливающиеся нелинейные эффекты при распространении волны в воде либо ткани будут проявляться различным образом. Чем ближе режим фокусировки в воде к разрывному и, тем более, к насыщению, тем больше будет ошибка при оценке параметров поля в ткани.
Проведение количественных расчётов уровней давления и интенсивности нелинейного фокусированного ультразвука в средах с различным поглощением, подтвержденных данными измерений, является поэтому, несомненно, важным как для определения эффективности воздействия, так и для оценки его безопасности. Исследование влияния нелинейных эффектов на изменение коэффициентов концентрации фокусирующих систем и достижение предельных значений акустических полей при фокусировке в средах с различным поглощением представляет интерес и с фундаментальной физической точки зрения. Количественное определение значений параметров нелинейных фокусированных полей стало возможным с использованием методов численного моделирования, которое может служить мощным дополняющим, а часто и заменяющим реальный эксперимент инструментом для калибровки полей преобразователей, использующихся в ультразвуковой хирургии, перенесения полученных в воде данных на ткань, предсказания возможных биологических эффектов в ткани и разработки стандартов использования HIFU.
До настоящего времени подробное исследование мощных фокусированных ультразвуковых пучков в широком диапазоне мощностей излучателей, применяемых в HIFU-хирургии, провести не удавалось. Большое количество работ посвящено исследованиям акустического поля нефокусирующих или слабо фокусирующих излучателей диагностического ультразвука; часто использовались идеализированные модели гауссовских излучателей, поле которых имеет гладкую структуру, удобную для моделирования, но далекую от реальности; либо исследовались режимы слабой нелинейности в фокусированных пучках. Расчеты и измерения нелинейных сильно фокусированных HIFU полей с учетом образования разрывов проведены лишь для нескольких преобразователей и выбранных конкретных режимов их работы.
В современных приборах, применяемых в ультразвуковой хирургии, уровень интенсивности в фокальной области достигает 25 кВт/см , что приводит к образованию разрывов большой амплитуды (до 80 МПа). При этом
возможен локальный сверхбыстрый, за несколько миллисекунд, нагрев ткани до температур выше 100С и возникновение кипения. До последнего времени эти эффекты сверхбыстрого нагрева до температур кипения в ткани, связанные с образованием разрывов, не исследовались. Однако они являются чрезвычайно важными с точки зрения решения практических медицинских задач, поскольку образование пузырьков пара при кипении в ткани кардинальным образом меняет процесс воздействия ультразвука на ткань.
В данной диссертационной работе численно исследованы нелинейные эффекты при фокусировке мощных ультразвуковых пучков в воде и в биологической ткани в широком диапазоне характерных для HIFU параметров излучателей различной геометрии и излучаемой мощности. Предложен и обоснован новый метод определения параметров нелинейного акустического поля в ткани, основанный на полученных в воде результатах. Последняя часть работы посвящена исследованию нагрева биологической ткани за счет поглощения на разрывах, что необходимо для разработки эффективных и безопасных протоколов неинвазивного хирургического воздействия ультразвука на ткань. Исследованы температурные поля в гелевом фантоме ткани и в реальной биологической ткани, рассчитаны и измерены времена нагрева среды до 100С для различных амплитуд разрывов в фокусе. Исследовано влияние диффузии тепла; показано, для каких интенсивностей она начинает играть существенную роль при расчете времени достижения температуры кипения.
Цели и задачи диссертационной работы
Основной целью настоящей работы стало развитие физических подходов, позволяющих исследовать структуру акустического поля и определить уровни давления в мощных фокусированных пучках, использующихся в современных устройствах для ультразвуковой хирургии, а также получить новые сведения о связи между волновым полем и физическими механизмами возникающих при этом биологических эффектов. В рамках указанной цели решались следующие конкретные задачи:
1. Разработка численного алгоритма, позволяющего моделировать задачи распространения мощных фокусированных нелинейных акустических волн в воде и биологической ткани с учетом образования узких ударных фронтов большой амплитуды и сильной пространственной локализации поля в области фокуса.
Численное исследование влияния нелинейно-дифракционных эффектов и аподизации поля на источнике на изменение пиковых значений и пространственной структуры различных характеристик нелинейных фокусированных полей, создаваемых источниками ультразвуковой хирургии в воде в широком диапазоне параметров таких источников и излучаемой мощности.
Исследование возможности перенесения результатов численного моделирования или измерений нелинейных волновых полей в воде на поглощающие среды. Разработка метода определения параметров нелинейных фокусированных ультразвуковых пучков в биологической ткани, основанного на измерениях или моделировании поля в воде.
Количественное исследование явления усиления теплового воздействия ультразвука на ткань за счет нелинейных эффектов в режимах, характерных для приборов, применяемых в неинвазивной хирургии. Исследование эффекта сверхбыстрого нагрева и возникновения локального кипения в биологической ткани в фокусе излучателя при воздействии разрывной волны.
Научная новизна работы
Впервые задача фокусировки мощных ультразвуковых пучков разрывных волн исследована в широком диапазоне параметров, характерных для устройств ультразвуковой хирургии. Построены калибровочные кривые, позволяющие рассчитать значения пиковых давлений и интенсивности в фокусе для различных режимов работы поршневых преобразователей в воде и поглощающей биологической ткани.
Показаны особенности изменения пространственной структуры различных характеристик ультразвуковых полей, связанные с образованием разрывов. Обнаружен эффект формирования двух разрывов на одном периоде волны в фокальной области пучка за счет интерференции прямой и краевой волн.
Предложен, обоснован и подтвержден в численном и физическом экспериментах новый метод определения параметров нелинейной сильно фокусированной акустической волны в биологической ткани на основе полученных в воде результатов моделирования или измерений.
Показано, что в условиях, характерных для современной клинической практики, образование разрывов и поглощение энергии волны на ударных фронтах могут приводить к локальному нагреву биологической ткани до
температур выше 100С и возникновению кипения в течение нескольких миллисекунд.
Достоверность полученных в работе результатов подтверждается данными экспериментов, проведенных при совместных исследованиях с Центром промышленного и медицинского ультразвука лаборатории прикладной физики университета шт. Вашингтон (Сиэтл, США), а также соответствием результатов теоретическим оценкам и данным численных расчетов, полученным в работах других авторов.
Практическая ценность
Предлагаемая работа, помимо результатов фундаментального характера, представляет собой основу для разработки новых стандартов описания полей, создаваемых устройствами неинвазивной хирургии, методов повышения безопасности и эффективности воздействия ультразвука на ткань, и, таким образом, будет способствовать дальнейшему внедрению этого метода для лечения больных.
Полученные калибровочные кривые по изменению коэффициентов концентрации поля в нелинейном пучке позволяют рассчитать величины пиковых давлений и интенсивности в фокусе произвольного поршневого ультразвукового преобразователя при любом уровне его возбуждения. Данные результаты представляют несомненную практическую важность и могут использоваться для определения величин акустических параметров нелинейного поля в фокусе и выбора оптимальных уровней облучения.
Представленный новый метод определения параметров нелинейного поля в фокусе ультразвукового преобразователя в ткани на основе данных, полученных в воде, чрезвычайно важен для планирования хирургической процедуры. Результаты моделирования могут быть использованы как альтернатива физическим измерениям и совместно с предложенным методом могут применяться для определения значений параметров поля в фокусе с высокой точностью даже при наличии ударных фронтов в профиле волны.
Предсказанный в численном моделировании и получивший подтверждение в эксперименте эффект сверхбыстрого нагревания ткани до температур выше 100С и возникновения кипения за несколько миллисекунд после начала воздействия ультразвуком позволяет использовать возможности диагностического ультразвука для визуализации области воздействия
HIFU при разработке клинических протоколов облучения, прицеливании, мониторинге воздействия в режиме реального времени.
Положения, выносимые на защиту:
Полученные количественные данные для коррекции коэффициентов концентрации и насыщения ультразвукового поля в фокусе нелинейного пучка позволяют рассчитать различные характеристики нелинейно-искаженных профилей волны в широком интервале параметров и мощностей фокусированных поршневых преобразователей, использующихся в устройствах ультразвуковой хирургии.
Интерференция прямой и краевой волн в фокусированных пучках, создаваемых поршневыми источниками, может приводить к формированию двух ударных фронтов на одном периоде волны при её распространении в режиме развитых разрывов.
Разработанный новый алгоритм перенесения результатов измерений акустического поля в воде на поглощающие среды, основанный на масштабировании амплитуды волны на излучателе, позволяет определять значения различных параметров акустического поля, таких как пиковые давления, интенсивность и мощность тепловыделения, в фокальной области излучателей для ультразвуковой хирургии.
Образование ударного фронта в профиле волны в фокальной области пучка в биологической ткани может приводить к локальному сверхбыстрому нагреванию ткани до температуры кипения и формированию паровых пузырьков в течение нескольких миллисекунд. Аналитические оценки тепловыделения, основанные на теории слабых ударных волн, позволяют предсказать время возникновения кипения в ткани с точностью 10%.
При измерении (с 10% точностью) в воде параметров разрывных профилей волн с интенсивностью до 25 кВт/см необходимы следующие требования к ширине полосы пропускания ультразвуковых приемников: 15 гармоник основной частоты для определения пикового отрицательного давления и средней интенсивности, более 40 гармоник - для пикового положительного давления, более 100 гармоник - для мощности тепловыделения.
Апробация работы
Вошедшие в диссертацию результаты докладывались на основных профильных российских и международных конференциях последних лет: международной конференции "Актуальные проблемы нелинейной волновой физики" (Нижний Новгород, 2005), 4-й совместной конференции акустического общества Америки и акустического общества Японии (Гавайи, Гонолулу, 2006), на X и XI Всероссийских школах - семинарах "Волновые явления в неоднородных средах" (Звенигород, 2006 и 2008), на XVIII, XIX и XX сессиях Российского акустического общества (Таганрог, 2006, Нижний Новгород, 2007, и Москва, 2008), на научной школе Therapeutic ultrasound (Франция, Каржез, 2007), на 14-ом «Международном конгрессе по звуку и колебаниям» (Австралия, Керне, 2007), на 7-ом Международном симпозиуме по терапевтическому ультразвуку - ISTU7 (Корея, Сеул, 2007), на 19-ом Международном конгрессе по акустике - ICA2007 (Испания, Мадрид, 2007), на совместной конференции Американского акустического общества и Европейской ассоциации акустиков «Акустика-08» (Франция, Париж, 2008), на 18-ом Международном симпозиуме по нелинейной акустике - ISNA18 (Швеция, Стокгольм, 2008), на 156й, 157й и 158й сессиях Американского акустического общества (Майами, октябрь 2008; Портленд, апрель 2009, и Сан-Антонио, октябрь 2009), на конкурсе студенческих и аспирантских работ при содействии центра энергосберегающих технологий (Москва, 2009), на Международном симпозиуме по акустике IEEE IUS (Италия, Рим, 2009), 9-ом Международном симпозиуме по терапевтическому ультразвуку - ISTU9 (Франция, Экс-ан-Прованс, 2009), на конференции «Фундаментальные и прикладные аспекты инновационных проектов физического факультета МГУ», а также обсуждались на научных семинарах Акустического института им. Н.Н. Андреева, кафедры акустики и лаборатории лазерной оптоакустики кафедры общей физики и волновых процессов физического факультета МГУ.
Публикации
По результатам исследований, представленных в диссертации, опубликовано 23 научные работы, список которых приведен в конце автореферата, в том числе 4 статьи в реферируемых журналах.
Структура и объем диссертации
Сравнение численных схем при моделировании плоских волн конечной амплитуды
Так как решение дифракционной задачи требует значительных затрат машинного времени, то основные особенности нелинейного распространения, формирования ударных фронтов и распространения разрывных волн рассмотрены вначале на примере плоской волны. Распространение плоских волн в среде с квадратичной нелинейностью и поглощением, которое является суммой квадратичного (как в воде) и линейного (как в биологической ткани) по частоте законов поглощения, описывается уравнением типа Бюргерса [1]: которое в безразмерных координатах имеет следующий вид: В качестве граничного условия задана гармоническая волна: P{z = О, в) = -sin в. Здесь р - давление волны, х - координата вдоль распространения волны, r = t-x/cQ — время в сопровождающей системе координат, є — параметр нелинейности среды, р0 - плотность среды, с0 - скорость звука, Ь — диссипативный коэффициент в воде, к — волновой вектор. Безразмерные переменные: Р = р/р0 - акустическое давление, нормированное на исходную амплитуду ро на источнике; z = х/х — координата распространения, нормированная на длину образования разрыва в плоской волне хр = сІр0/єр0й)0, со0 = 2тт/0 - циклическая частота волны; в — со0т — безразмерное время; L(P) — оператор поглощения, имеющий сложный интегральный вид в случае линейного закона поглощения, характерного для биологической ткани; параметр А-хр/хп учитывает относительный вклад нелинейных и диссипативных эффектов. Для сред с квадратичным поглощением, д:п =2p0cl/ba 2 и, таким образом, Ав =Ьа 0/єр0 характеризует поглощение в воде на основной частоте. Для сред с линейным по частоте законом поглощения = 71 (сет - коэффициент поглощения на частоте излучателя), следовательно, Ат = сІр0ат/єр0о)0 определяет поглощение в ткани. При моделировании распространения волны в воде коэффициент поглощения в ткани полагался равным нулю Ат = 0.
При моделировании распространения волны в ткани частотно зависимое поглощение в ткани a{f) в численной модели являлось комбинацией линейного и слабого квадратичного закона поглощения, как в воде: Дисперсия скорости звука рассчитывалась для линейного по частоте закона поглощения с помощью локальных дисперсионных соотношений типа Крамерса-Кронига [85, 103]: Остановимся здесь кратко на основных деталях используемых алгоритмов. Уравнение (2.2) решалось на каждом шаге сетки по z методом расщепления по физическим факторам (поглощения и нелинейности): Оператор поглощения Ьпогп удобнее рассчитывать в спектральном представлении. Для этого решение для акустического давления P{z, в) представляется в виде разложения Фурье: где С„ - комплексная амплитуда п-\л гармоники в спектре волны. Оператор поглощения для каждой из гармоник Сп имеет следующий вид: где К(п) = k(na 0)F = К + ІК" -комплексное волновое число в движущейся системе координат; К и ІТ"- безразмерные действительная и мнимая части К(п):Частотная зависимость поглощения К" в (2.8) моделируется соотношением (2.3), а дисперсия скорости звука К рассчитывается по известной зависимости (2.4). Для уменьшения числа гармоник в моделировании и обеспечения при этом устойчивости численной схемы на стадии развитых разрывов в схему вводится дополнительная искусственная вязкость, изменяющаяся по следующему закону: K"(n)1KK=T-AB-n2exp[-((z-7r/2)/l)2] при2 тг/2и К"(п)1кк =Т-Ап-п2 при z тг/2. Как видно, вязкость, пропорциональная квадрату частоты, возрастает при приближении волны к точке, в которой амплитуда разрыва максимальна. Параметр / характеризует длину области, в которой амплитуда разрыва возрастает от нуля до максимального значения. В приведенных расчетах для плоской волны параметр / = 0.8 был подобран эмпирически таким образом, чтобы кривые, построенные в случае внесения в схему искусственного поглощения и без него, имели минимальное отличие. Параметр Т определяет величину введенного дополнительного поглощения, которое приводит к сглаживанию ударных фронтов и обеспечению устойчивости схемы при уменьшении числа гармоник. Например, если планируется уменьшить количество гармоник в расчете в два раза по сравнению с устойчивым решением, полученным без искусственного поглощения, во столько же раз надо увеличить параметр Т. Отметим, что для схем Годунова введение искусственного поглощения не является необходимым на этапе решения задачи о распространении плоской волны. Однако, как будет показано ниже, это становится необходимым в дальнейшем, при рассмотрении фокусировки разрывных волн. Нелинейный оператор L„e.1UH рассчитывался двумя разными способами: в спектральном представлении (с помощью схемы Рунге-Кутта) и во временном представлении (с использованием схем типа Годунова первого и второго порядков точности). Подробнее об этом будет написано в пунктах 2.2.1 и 2.2.2 данной главы. Результатом численного моделирования являлись амплитуды гармоник Сп и профили волны в каждом узле координатной сетки. Анализируя профиль волны, определялись величины пиковых давлений. Полная интенсивность волны вычислялась как сумма интенсивностей всех гармоник Тп (г): Мощность тепловыделения в среде за счет поглощения энергии волны вычислялась на каждом шаге сетки по z как разность интенсивностей до и после расчета нелинейного и диссипативного операторов. Это позволяет учесть поглощение на ударном фронте, вносимое через собственную вязкость схемы Годунова. В случае использования спектрального подхода мощность тепловых источников также рассчитывалась как сумма тепловыделения за счет поглощения каждой из гармоник:
Численное моделирование распространения плоской волны проводилось для излучателя с частотой = 2 МГц и начальной амплитудой давления ро = 50 МПа. В этом случае расстояние образования разрыва равно 1.38 мм, а на расстоянии z = яг/2 образуется пилообразная волна с амплитудой разрыва 100 МПа. Эта величина соответствует максимальной амплитуде ударного фронта, измеренного в фокусе преобразователя в эксперименте, о котором будет сказано в 2.4. При моделировании в воде использовались следующие физические параметры: плотность /?о=1000 кг/м , скорость звука ?о=1486 м/с; нелинейность среды = 3.5; коэффициент поглощения Ъ = 4.3-10-3 КГ-C -M"1. В этом случае безразмерный параметр поглощения в воде равен Ав = 1.6 10"4. При моделировании распространения плоской волны в биологической ткани были выбраны два характерных значения коэффициента поглощения аг =0.1 Нп/см и сст = 0.05 Нп/см на 1 МГц, характерных для поглощения в мягких тканях [12]. Остальные физические параметры среды (плотность, скорость звука и нелинейность) предполагались такими же, как и в воде. В этом случае Ат = 3.1 10-2 и А,. = 1.5 10"2. Традиционно в спектральном представлении решается система связанных нелинейных уравнений для амплитуд первых Nmax гармоник, а амплитуды остальных гармоник с большими номерами п Nmax полагаются равными нулю. Решение ищется в виде разложения Фурье (2.6), которое при расчете нелинейного оператора дает связанную систему для Фурье-компонент: Эта система интегрируется численно для конечного числа гармоник Nmax (2Nmax числа точек на профиле волны) методом Рунге-Кутта. Данный метод хорош тем, что позволяет строить схемы различного порядка точности. Его основная идея состоит в построении специального алгоритма — такого, чтобы приращение функции на шаге с1С„ Сп+\ — Сп совпадало с приращением dCn, которое определяется из ряда
Тейлора с учетом возможно большего числа членов. При этом вторые и следующие производные определяются не дифференцированием, а путем многократного вычисления функции /(С, z) в некоторых промежуточных точках между z и z+hz. В настоящее время наиболее распространена схема четвертого порядка точности [104], которая и используется в данной работе. Таким методом удобно описывать распространение слабо искаженных сигналов с достаточно узким спектром. Однако взаимное действие диссипативных и нелинейных эффектов приводит к искажению профиля волны и образованию ударного фронта конечной ширины, внутренняя структура которого определяется коэффициентом поглощения в воде Ак: Р = th(0/2AB), —л 6 п [1]. При образовании в профиле волны участков с большим значением производной (узких ударных фронтов) спектр сигнала расширяется. Следовательно, необходимо учитывать все большее и большее количество гармоник, что непосредственно отражается на устойчивости схемы, времени и точности расчета. Ситуация осложняется еще и тем, что все гармоники эффективно взаимодействуют, и по мере распространения волны ошибка, возникающая при определении высокочастотной части, сказывается на всем спектре. На рис. 2.1 показаны профили волны в точке z = я/2, рассчитанные с помощью различного числа гармоник в воде (/1,,= 1.6-10"4, Ат = 0). Как известно, на данном расстоянии амплитуда сформировавшегося разрыва максимальна. Шаг hz вдоль координаты распространения изменялся от 1.25-10"5 для Nmax — 4000 до 5-Ю"6 для Nmwc= 10000. Таким образом, зависимость между шагом по времени и по координате распространения имеет следующий вид: hz = 0.016-Ив. С увеличением числа гармоник растет количество точек на профиле от 1 для Nmax = 4000 до 7 для Nmax= 10000. Можно видеть, что решение становится устойчивым только при Nmax 10000. Во всех остальных случаях, в восстановленном профиле присутствуют осцилляции Гиббса, обусловленные пренебрежением в расчетах высокочастотной частью спектра. Чтобы избежать этого явления, то есть уменьшить амплитуду высокочастотных гармоник, в решение можно вводить искусственную вязкость.
Влияние ограниченной полосы пропускания гидрофона на точность измерения параметров акустической волны в режиме развитых разрывов
Результаты численного моделирования фокусировки пучка использовались для подготовки эксперимента по измерению профилей волны в фокусе излучателя, подтверждения предсказанных и объяснения наблюдаемых эффектов, а также для проверки численного решения: насколько оно соответствует физическому эксперименту. Кроме того, результаты численных расчетов можно использовать для понимания ограничений, которые вносят измерения. Если в численном моделировании ограничение на число гармоник определяется только мощностью компьютера, то в эксперименте точность измерений связана с конечностью частотной полосы пропускания приемников ультразвука [76,77]. Только для современных тонких оптоволоконных гидрофонов, доступных для HIFU экспериментов, полоса пропускания составляет 100 МГц, а для стандартных ПВДФ гидрофонов ограничена 30 МГц. Учитывая, что излучатели мощного ультразвука обычно работают в диапазоне частот 1-3 МГц, в эксперименте удается измерить не более 100 гармоник. Как было показано в 2.2-2.3, при формировании узкого ударного фронта в профиле волны, её спектр может содержать несколько тысяч гармоник. Естественно возникает вопрос, с какой точностью производятся измерения профилей волны в эксперименте и насколько верными оказываются полученные значения. Поэтому цель данного параграфа состоит в том, чтобы определить, сколько гармоник необходимо для восстановления различных параметров волны, таких как пиковые давления, интенсивность и мощность тепловыделения, с точностью 10%. Уравнение ХЗК (2.23) моделировалось в воде в спектральном представлении (так как в этом случае необходимо пользоваться формулой (2.11)), учитывая 500 гармоник в спектре волны. Затем рассчитывались различные параметры волны (пиковые давления, интенсивность, мощность тепловых источников), учитывая различное число гармоник численного решения.
Результаты моделирования сравнивались с экспериментальными данными, полученными на основе измеренной формы волны в воде. Эксперименты проводились в Центре медицинского и промышленного ультразвука лаборатории прикладной физики университета шт. Вашингтон в Сиэтле. Автор диссертационной работы участвовал в экспериментальной части работы на этапе подготовки эксперимента, при анализе и обсуждении полученных данных. Схема установки, использующейся для проведения измерений мощных акустических полей, показана на рис. 2.14 [79]. Источник ультразвука представляет собой одноэлементный сферически сфокусированный пьезокерамический излучатель (PZ 26, Ferroperm Piezoceramics, Denmark) с резонансной частотой f0 = 2.158 МГц, радиусом зо = 21 мм и фокусным расстоянием F = 44.4 мм, что соответствует линейному коэффициенту фокусировки G = 48. Для измерений использовался импульсный режим работы; чтобы избежать развития кавитации на излучатель подавалась последовательность коротких импульсов (30 - 40 периодов волны), которые генерировались с помощью функционального преобразователя (Agilent 33250А, Agilent, Palo Alto, CA) и линейного усилителя высокой частоты (ENI А-300, ENI, Rochester, NY). Для измерения профилей волны в фокусе излучателя в воде и в гелевом полиакриламидном фантоме ткани использовались оптоволоконные гидрофоны (RP Acoustics, Germany) с активным диаметром 100 мкм и с полосой пропускания 30 МГц (гидрофон FOPH500) и 100 МГц (FOPH2000). Таким образом, для излучателя частотой 2.158 МГц первым гидрофоном могло быть измерено только 15 гармоник, вторым - 50 гармоник. Сигнал гидрофона подавался на преобразователь, выходным параметром которого являлось напряжение, и записывался с помощью цифрового осциллографа (Lecroy LT344, Chestnut Ridge, NY). Источник ультразвука был прикреплен к трехмерной позиционирующей системе (Velmex Inc., Bloomfield, NY), контролирующей его положение относительно гидрофона. Компьютерная программа LABVIEW (National Instruments, Austin, TX) использовалась для контроля генератора, осциллографа и позиционирующей системы во время измерений.
Все эксперименты были проведены в большом резервуаре с водой комнатной температуры (20С). Вода предварительно очищалась и дегазировалась во избежание возникновения кавитации. Для сравнения данных эксперимента с численными расчетами моделировалось уравнение ХЗК (2.21) с модифицированным граничным условием [79]: Размер излучателя (2«о = 44.8 мм) и пространственная аподизация амплитуды давления на нем в граничном условии (2.26) модифицированы таким образом, чтобы численное решение в фокальной области пучка наиболее точно соответствовало данным измерений при линейной фокусировке волны. Рисунок 2.15 показывает частотный отклик гидрофонов (а), сравнение профилей волны давления (б) и их спектров (в) в фокусе, промоделированных (сплошная линия) и измеренных двумя разными гидрофонами: FOPH 500 (пунктирная линия, полоса пропускания 30 МГц) и FOPH 2000 (штрих-пунктирная линия, полоса пропускания 100 МГц). Амплитуда давления на источнике была ро - 0.4 МПа (/о= 5 Вт/см ), что соответствует случаю, когда ударный фронт формируется близко к фокусу [108]. Результаты эксперимента для FOPH 2000 показывают, что пиковое положительное давление более чем на 15 МПа ниже, чем полученное в моделировании: р+ = 68 МПа в численном расчете ир+ = 52 МПа в эксперименте. Для FOPH 500, пиковое положительное давление даёт ещё более низкие результаты по сравнению с моделированием: 43 МПа.
С другой стороны, пиковые отрицательные давления практически совпадают между собой на всех трех профилях. Эти графики демонстрируют, что полоса пропускания гидрофона является важной характеристикой для точных измерений ударного фронта и пикового положительного давления. Также можно видеть, что, несмотря на увеличенную полосу пропускания современного гидрофона FOPH 2000, ударный фронт в фокусе HIFU излучателя всё же не может быть измерен с достаточной точностью. Зная спектр волны р,ъ можно восстановить параметры волны, учитывая различное число гармоник в решении. Таким образом, можно рассчитать профиль, интенсивность и мощность тепловыделения волны как функции различного числа гармоник N: где N меняется от 1 до 500 в моделировании или 1 до 50 в эксперименте для гидрофона FOPH 2000. Чем больше гармоник использовать в решении, тем ближе будет результат к тому значению, которое получено при прямых измерениях или расчетах. Рисунок 2.16 иллюстрирует эти зависимости для пиковых давлений (а), интенсивности волны (б) и мощность тепловых источников (в). Горизонтальные штриховые линии на графиках соответствуют уровням на 10% ниже, чем при расчете параметров с максимальным количеством гармоник (500). Видно, что для восстановления интенсивности и пикового отрицательного давления с 10%-й точностью необходимо не более 10 гармоник. Пиковое положительное давление более чувствительно к полосе пропускания гидрофона, поэтому более чем 40 гармоник необходимо, чтобы определить его значения с той же точностью. Для мощности тепловыделения наблюдается большая разница между результатами, полученными в эксперименте и в численном моделировании. Как уже было показано в 2.3-2.4, это самый чувствительный параметр к условиям измерения или моделирования, и для его точного предсказания необходимо знать более сотни гармоник, что в настоящее время невозможно обеспечить в измерениях. Следует отметить, что помимо полосы пропускания гидрофонов существует ещё ряд факторов, ограничивающих точность получаемых результатов в эксперименте. Во-первых, точность калибровки современных гидрофонов составляет порядка 10%. Кроме того, в случае использования оптоволоконного гидрофона, результаты будут зависеть от наклона среза оптоволокна по отношению к фронту волны, от развития кавитации на нем. Использование результатов численного моделирования позволяет избежать всех этих проблем и определять параметры акустического поля с высокой точностью.
Пространственное распределение параметров акустического поля
Различные параметры акустического поля ответственны за различные физические механизмы воздействия, поэтому безусловный интерес представляет задача исследования пространственной локализации тех или иных акустических величин при различных условиях фокусировки [120]. Характерный вид распределений различных безразмерных параметров фокусированной волны вдоль оси пучка показан на рис. 3.10 для случая линейной фокусировки (/V = 0); для случая, когда ударный фронт образуется в профиле волны вблизи фокуса (N = 0.33) и до фокуса (N=1.17). Графики представлены для слабофокусирующего излучателя с линейным коэффициентом фокусировки G— 10. В этом случае, при изменении параметров поля будут наблюдаться те же особенности, что и для излучателей с большим коэффициентом фокусировки, однако графики выглядят более наглядно. Здесь рассматривается распространение волны в непоглощающей среде (воде), поэтому заметный переход энергии волны в тепло будет происходить лишь после образования ударного фронта. Быстрый рост кривой тепловыделения на рис. 3.10 (г) соответствует расстоянию, на котором формируется ударный фронт. Из рисунка видно, что по мере усиления нелинейных эффектов (увеличение N) значение акустического поля в фокальной области меняется по-разному для различных параметров волны, о чем подробно говорилось в 3.2. Фокальные значения пикового положительного давления Р+ и интенсивности волны / (рис. 3.10 (а, в)) вначале возрастают при увеличении N до образования разрыва (N = 0.33), а затем уменьшаются за счет поглощения энергии на этих разрывах (ЛГ= 1.17). При этом пиковое отрицательное давление Р- в фокусе монотонно убывает по мере увеличения N, а мощность тепловых источников О, Из рис. 3.10 также видно, что положения пространственных максимумов для различных параметров волны, совпадающие при N-0, смещаются вдоль оси пучка с увеличением JV различным образом, причем в некоторых случаях немонотонно. Максимум мощности тепловыделения (рис. 3.10 (г)) с ростом N монотонно смещается в сторону излучателя. Как уже говорилось в предыдущем параграфе, значительный рост мощности тепловыделения начинается только после формирования разрыва в профиле волны.
При увеличении амплитуды волны (то есть параметра N) ударный фронт начинает формироваться всё ближе к излучателю, поэтому по положению максимума мощности тепловыделения на оси пучка можно судить о том, на каком расстоянии от излучателя в профиле волны уже образовался разрыв. Максимум положительного пикового давления и интенсивности вначале, наоборот, смещается в сторону геометрического фокуса, переходя даже в каком-то интервале значений N за фокус, а затем - в сторону излучателя. Более подробно за таким смещением максимума параметра Р+ можно проследить на рис. 3.11. Оно объясняется влиянием эффекта самодефокусировки пучка, происходящей за счет несимметричного искажения профиля и увеличения скорости распространения участков фазы сжатия волны вблизи оси распространения. При слабом искажении профиля (N=0.25) эффект самодефокусировки мал и приводит к смещению максимума Р+ от излучателя. В случае максимальной асимметрии (JV=0.33) эффект дефокусировки усиливается и максимум Р+ смещается в сторону излучателя. При дальнейшем области приводит к дополнительному Давления + смещению максимума в сторону излучателя. Скорость распространения участков фазы разрежения волны вблизи оси пучка с ростом N, напротив, уменьшается, что приводит к усилению самофокусировки и монотонному смещению максимума Р- к излучателю (рис. 3.10 (б)). На рис. 3.10 и 3.11 также можно видеть характерный «излом» вблизи максимума кривой давления Р+. После стремительного нарастания давления до максимума происходит его быстрое уменьшение, особенно хорошо заметное при больших нелинейностях (АГ» 1.17 на рис. 3.10). Такой эффект появляется после того, как в профиле волны сформируется ударный фронт, на котором начинается эффективное поглощение энергии. Стоит отметить, что он наблюдается для сильно [121] и слабо фокусированных излучателей [122], а также в случае плоских преобразователей [123]. Ещё одна интересная особенность наблюдается в профиле волны при моделировании на высоких уровнях нелинейности (в два раза больших тех, при которых разрыв образуется в профиле волны вблизи фокуса). На рис. 3.12 представлены распределение пиковых давлений вдоль оси излучателя и профили волны вблизи геометрического фокуса излучателя для сильно фокусированного преобразователя G = 40 и параметра нелинейности N=0.5. На рис. 3.12 (б) ясно видно, что профиль волны содержит два узких ударных фронта, которые при распространении постепенно сближаются и около фокуса сливаются в один разрыв. Ранее существование такого явления уже было предсказано и подтверждено численно и экспериментально для ближнего поля плоского поршневого излучателя [28, 123]. Как известно, поглощение энергии на ударном фронте пропорционально амплитуде разрыва в третьей степени.
Когда на периоде волны наблюдается два разрыва, поглощение энергии пропорционально сумме кубов двух амплитуд, АЪ5Л + Аъвл. После слияния разрывов поглощение пропорционально уже кубу суммы двух амплитуд: (Asl+As2)3, которое, очевидно, значительно больше, чем сумма кубов амплитуд разрыва А\л + A3S2. Поэтому за счет нелинейных эффектов и поглощения энергии волны на образовавшемся разрыве, амплитуда волны резко уменьшается после слияния двух ударных фронтов в один. Так как данный эффект наблюдается до геометрического фокуса, то далее происходит медленное увеличение давления за счет дифракционного усиления вблизи фокуса. разрывов на одном периоде волны вызвано влиянием дифракции, в частности, приходом краевых волн. На рис. 3.13 представлены распределения профилей волны (а) и производной от профиля волны по времени dP{z,6,R)/dd (б) в пространстве, то есть на различных расстояниях R от оси излучателя. Видно, что образование двух разрывов на одном периоде волны происходит не только вдоль оси излучателя, но и на некотором расстоянии от неё. «Ступеньке» на рис. 3.12 (б) соответствует треугольная голубая область на рис. 3.13 (а) или светлая область на рис. 3.13 (б). Стоит отмстить, что явление образования двух разрывов в профиле волны наблюдается как в слабопоглощающей среде (воде), так и в биологической ткани, как для слабо фокусированных излучателей, так и для преобразователей с большим коэффициентом усиления. По мере усиления нелинейных эффектов происходит не только изменение коэффициентов концентрации и смещение пространственных максимумов для различных акустических параметров пучка, но меняется также их пространственная структура [124]. На рис. 3.14 показаны пространственные распределения положительного Р+ и отрицательного Р_ пиковых значений давления, интенсивности / и мощности тепловыделения О в плоскости, проходящей через ось пучка, в случае линейной (а, б) и нелинейной (в-е) фокусировки (G=40). Значение параметра нелинейности JV=0.25 соответствует случаю максимального усиления коэффициента концентрации для положительного пикового давления (рис. 3.4). Графики представлены в линейном масштабе с восемью равными уровнями, изменяющимися от нуля до максимума амплитуды соответствующей величины.
Для линейного пучка (колонка слева) показаны только распределения для амплитуды давления Р± и интенсивности / , поскольку в этом случае распределения пикового положительного и отрицательного давлений совпадают, а распределения интенсивности волны и мощности тепловыделения пропорциональны друг другу, т.е. также имеют одинаковую пространственную структуру. Как видно из рисунков, в нелинейном поле фокальная область положительного пикового давления (в) и особенно тепловых источников (е) существенно более локализована в пространстве по сравнению с линейным полем. Напротив, фокальная область отрицательного пикового давления (д) заметно смещена в сторону излучателя и гораздо менее локализована в пространстве, особенно в поперечном направлении. Распределение интенсивности в нелинейном пучке (г), на удивление, мало отличается от линейного случая, хотя в фокальном профиле волны уже образовался ударный фронт.
Метод определения параметров акустического поля в биологической ткани с учетом нелинейных эффектов
Предложенный метод основан на том, что преобразователи, использующиеся в HEFU терапии, обычно имеют коэффициенты усиления по давлению 20-50; поэтому амплитуда поля в фокальной области существенно выше, чем на пути к фокусу, а длина фокальной области - намного короче, чем фокусное расстояние излучателя [12, 79]. При таких условиях разумно предположить, что нелинейные эффекты в фокальной области будут преобладать над эффектами нелинейного распространения до фокуса. Тогда степень нелинейного искажения формы волны в фокусе будет определяться в основном уровнем давления в фокальной области. Следуя данным предположениям, рассмотрим вначале фокусировку линейного пучка в воде и в ткани при таких давлениях на источнике, когда амплитуда давления pF в геометрическом фокусе преобразователя будет одинаковой в обеих средах. Предполагается, что все акустические параметры воды и ткани, за исключением поглощения, одинаковы. На рис. 4.1 показаны распределения амплитуд давления вдоль оси пучка для одноэлементного поршневого источника с линейным-коэффициентом концентрации G=pp/pQ = 40. Решения получены численно с помощью линеаризованного уравнения ХЗК (2.23). Амплитуда акустического давления р нормирована на амплитуду давления на излучателе ро в воде: Р =р/ро. Сплошная кривая соответствует распространению волны в ткани, пунктирная — в воде. Если безразмерный коэффициент поглощения в ткани на частоте излучателя равен Ат, то амплитуда давления на источнике в случае распространения волны в ткани должна быть в ехр(Ат) раз больше, чем в случае распространения в воде (на рис. 4.1 параметр Ат = 1 и ехр(/4т)= 2.7). В этом случае давления в фокусе будут одинаковыми в обеих средах.
В работе проверяется следующая гипотеза: если увеличить амплитуду давления на излучателе в нелинейном пучке таким же образом, т.е. то профили волны в фокусе также будут близки как по форме, так и по пиковым значениям давления. Если эта гипотеза верна, то профили волны в фокусе, пиковые давления и интенсивность нелинейных ультразвуковых полей в биологической ткани могут быть определены из измерений или моделирования акустического поля в воде, используя лишь линейный коэффициент поглощения Ат. В данной главе вычисления проводились на основе решения уравнения (2.23) в спектральном представлении. Так как метод дирейтинга предложен не только для пиковых давлений и интенсивности, но и для мощности тепловыделения (раздел 4.4.3) с использованием формулы (2.11), то для её правильного вычисления необходимо рассчитывать нелинейный оператор в спектральном представлении. Большинство результатов приведено здесь для параметра поглощения Ат = 1 (см. 2.3), но также проводится сравнение с численными результатами для параметров Ат = 0.25,0.5 и 0.75. На рис. 4-.2 показаны зависимости пикового положительного и отрицательного давлений (рис. 4.2 а) и интенсивности волны (рис. 4.2 б), рассчитанных в ткани (сплошная кривая) и в воде (пунктир) вдоль оси пучка. На малых графиках показаны профили волны (а) и амплитуды первых шести гармоник (б), рассчитанные в воде и в ткани в фокусе излучателя, z= 1. Расчёты в воде проводились для граничного условия (4.5), линейного коэффициента усиления С? = 40 и нелинейного параметра JV=0.25.
При этих условиях ударный фронт формируется вблизи фокуса, а нелинейное увеличение коэффициента усиления для пикового положительного давления максимально (рис. 3.4). Расчеты в поглощающей ткани (АТ=1) проводились для того же значения коэффициента G, но давление на источнике, и соответственно, нелинейный параметр N, были в ехр(1) = 2.7 раз выше, в соответствии с величиной параметра поглощения Ах=\. Давления и интенсивность на рисунках нормированы на соответствующие значения на излучателе, используемые при моделировании в воде, поэтому безразмерная амплитуда давления при z = 0 в ткани равна 2.7, а интенсивности - 7.4. В случае линейного распространения волны, такое увеличение исходного давления при распространении в ткани точно бы компенсировало потери энергии на пути к фокусу, и значения давления и интенсивности в фокусе преобразователя в воде и & ткани были бы одинаковыми (рис. 4.1). Результаты, представленные на рис. 4.2, соответствуют нелинейному распространению волны и образованию ударных- фронтов в профиле волны. Однако видно, что даже в этом случае величины давления и интенсивности в фокусе излучателя в воде и в ткани близки друг к другу. Данные расчетов в ткани несколько ниже, чем в воде, за счет более сильного поглощения в ткани гармоник, генерируемых на пути к фокусу. Однако эти отличия практически незаметны для интенсивности и пикового отрицательного давления, а для пикового положительного давления разница составляет меньше 10%. Профили волны в фокусе (рис. 4.2 а) и спектры (рис. 4.2 б) также очень близки. Отметим, что хотя пиковое положительное давление в профиле волны, рассчитанной в воде (пунктир), на 10% выше, чем в ткани- (сплошная кривая), амплитуды разрыва в обоих профилях практически
Похожие диссертации на Нелинейные эффекты в мощных фокусированных ультразвуковых пучках: моделирование и применение в неинвазивной хирургии
