CC BY
38
Низкочастотная ультразвуковая томография: математические методы и эксперимент
А. В. Гончарский, С. Ю. Романов," С. Ю. Серёжников
Московский государственный университет имени М.В. Ломоносова, Научно-исследовательский вычислительный центр. Россия, 119234, Москва, Ленинские горы, д. 1, стр. 4.
Поступила в редакцию 18.05.2018, после доработки 01.10.2018, принята к публикации 29.10.2018.
Статья посвящена исследованию возможностей низкочастотной ультразвуковой томографии в медицине. Основным приложением является разработка томографических методов дифференциальной диагностики рака молочной железы. Для решения обратной задачи использована математическая модель, описывающая как явления дифракции и рефракции, так и поглощение ультразвука в неоднородной среде. Разработанные алгоритмы реконструкции скоростного разреза тестировались на стенде для экспериментальных исследований. Для зондирования использовались широкополосные ультразвуковые источники в диапазоне рабочих частот от 50 кГц до 600 кГц. В экспериментах на стенде использовался гидрофон ТеЫупе Яе80п ТС4038 с чувствительностью ~ 4 мкВ/Па, рабочей полосой частот 10800 кГц. Для усиления сигналов использовался предусилитель ТеЫупе Яе80п УР1000. Исследования проводились на фантомах с акустическими параметрами близкими к значениям параметров в мягких тканях человека. Для реконструкции 3Б скоростного разреза использовалась послойная модель. Реконструированные по экспериментальным данным сечения скоростного разреза имеют разрешение ~ 2 мм при центральной длине волны ~ 4 мм. Важным результатом работы является экспериментальное подтверждение адекватности модели физическим процессам.
Ключевые слова: ультразвуковая томография, волновое уравнение, обратные задачи, гидрофон, суперкомпьютерные технологии.
УДК: 534.2, 519.6, 517.958. PACS: 02.30.Zz, 02.60.Cb, 43.35.Wa.
ВВЕДЕНИЕ
Первоочередной задачей в медицине является разработка ультразвуковых томографов для дифференциальной диагностики рака молочной железы. Смертность среди женской половины человечества от этого заболевания стоит на первом месте среди всех онкологических заболеваний. Прогресс в лечении рака молочной железы напрямую связан с ранней диагностикой, которую можно обеспечить с помощью регулярных обследований. Для дифференциальной диагностики рака молочной железы используются различные технологии, в том числе и томографические. Рентгеновские томографы, несмотря на все свои достоинства, нельзя использовать для регулярных обследований из-за высокой лучевой нагрузки. Магнитно-резонансная томография и позитронно-эмиссионная томография наиболее эффективны при использовании контрастных агентов, которые также не рекомендованы при регулярных обследованиях. Существующие в медицине диагностические приборы ультразвуковых исследований (УЗИ) регистрируют сигнал на отражение и могут определять лишь границы неоднородностей. Получить данные о величине скорости звука внутри объекта, необходимые для характеризации тканей, в схемах на отражение невозможно [1]. Проектируемые в настоящее время ультразвуковые томографы, использующие как отраженое, так и проходящее излучения, позволяют осуществлять характеризацию мягких тканей. Разработки ультразвуковых томографов ведутся в США, Европе, Японии, России уже более 10 лет. На настоящий момент промышленно выпускаемых ультразвуковых томографов не существует. Разработки находятся на стадии макетов, стендов и опытных образцов [2-7]. В большинстве разработок диапазон частот ультразвуковых источников превышает 1 МГц. В настоящей работе рассматривается низкочастотная
а Е-шаИ: [email protected]
томография в диапазоне частот до 600 кГц. Кроме частот, отличаются методы и алгоритмы реконструкции скоростного разреза. Разработка ультразвуковых томографов сопряжена с целым рядом трудностей. По сравнению с рентгеновской томографией обратные задачи ультразвуковой томографии являются намного более сложными. Математическая модель в ультразвуковой томографии должна учитывать такие эффекты, как дифракцию, рефракцию, переотражение волн и т. п. Обратные задачи ультразвуковой томографии в таких моделях являются трехмерными и нелинейными. В работах [8-15] разработаны итерационные методы решения обратных задач волновой томографии, которые апробированы на многочисленных модельных расчетах. В работе [16] предложено направление низкочастотной томографии и показана его эффективность на модельных примерах. Однако модельные задачи не могут дать ответ об адекватности используемой модели реальным физическим процессам. Ответ на этот вопрос может дать только эксперимент. В настоящей работе обсуждаются результаты томографических исследований, полученных на разработанном и изготовленном стенде для экспериментальных исследований по низкочастотной ультразвуковой томографии. Алгоритмы решения обратных задач ультразвуковой томографии апробированы на экспериментальных данных, полученных на стенде, в послойной модели томографии. Результаты реконструкции скоростного разреза по экспериментальным данным показывают, что методы ультразвуковой низкочастотной томографии обеспечивают разрешение ~ 2 мм при центральной длине волны ~ 4 мм. Сравнимое разрешение томографии по проходящему излучению продемонстрировано в работе [17], где использованы более высокие частоты, а обратная задача решается в упрощенной параболической модели. Важным результатом работы является экспериментальное подтверждение адекватности модели физическим процессам.
1. ПОСТАНОВКА ОБРАТНОЙ ЗАДАЧИ ПОСЛОЙНОЙ УЛЬТРАЗВУКОВОЙ ТОМОГРАФИИ И МЕТОДЫ ЕЕ РЕШЕНИЯ
Основным объектом исследования в настоящей работе являются мягкие ткани человека, в которых поперечные волны практически отсутствуют. Поэтому в работе используется скалярная волновая модель, которая описывает волновые эффекты, такие как дифракция, рефракция, переотражение и поглощение ультразвука. В скалярной волновой модели акустическое давление и(г, £) описывается дифференциальным уравнением гиперболического типа. Скалярная модель позволяет рассчитывать акустическое давление и(г, £) г е (Ж = 2, 3) по заданным начальным данным, пользуясь уравнениями
е(г)ий(г, £) + а(г)и4(г, £) — Аи(г, £) =
= 6(г — го) ■ д(£) и(г, £)^=о = 0, и (г,£)|<=0 = 0. (2)
(1)
Здесь c(r) = 1/v2(r), v(r) — скорость волны в среде, а(г) — коэффициент поглощения в среде, А — оператор Лапласа по переменной r. Генерируемый источником в точке г0 импульс описывается функцией g(t). Расчет волнового поля при заданном c(r), a(r) из уравнений (1), (2) представляет собой прямую задачу.
Приведем постановку обратной задачи ультразвуковой томографии. Рассматриваемая обратная задача является коэффициентной обратной задачей, в которой необходимо по данным измерений волнового поля на поверхности окружающей объект исследования восстановить коэффициенты c(r) и a(r) в уравнении (1). В статье рассматривается послойная схема 3D ультразвуковой томографии. Схема эксперимента в отдельном слое приведена на рис. 1. Зондирующие акустические импульсы поочередно излучаются источниками ультразвукового излучения, обозначенными цифрой 1 на рис. 1. Цифрой 2 на рис. 1 обозначены приемники, регистрирующие акустическое давление распространяющегося от источников волнового поля как функцию от времени u(t). Такая схема эксперимента позволяет измерять волновое поле как прошедшее через объект исследования G, так и отраженное от объекта. Детальное измерение волнового поля осуществляется приемниками, расположенными на окружности Г с шагом, не превышающим половину центральной длины волны зондирующего импульса. Источники и приемники располагаются в плоскости z = const в водоподобной среде L, окружающей объект G. Скорость распространения волны v(r) и коэффициент поглощения a(r)
1
ю,
■ ..... L ,
♦ \
• ш *
в области G неизвестны. В среде L скорость волны постоянна и известна v(r) = const = v0, а коэффициент поглощения равен 0.
Таким образом, обратная задача ультразвуковых томографии заключается в восстановлении неизвестной функции скорости звука v(r) и коэффициента поглощения a(r) из уравнения (1) в исследуемой области G, по данным измерения зондирующей ультразвуковой волны на границе Г. В используемой в настоящей статье послойной схеме 3D ультразвуковой томографии обратная задача решается независимо в каждой из плоскостей z = const.
Математически обратная задача сводится к следующему. Запишем функционал невязки Ф(с, a) от аргументов c(r), a(r) между экспериментальными данными и рассчитанным полем на границе Г
M T Г
Ф(с,а) = £ J J(uj (y, t; c, а) - Uj (7,t))2 d7dt. (3)
j=1 o г
Здесь Uj (Y,t) — экспериментальные данные на границе Г за время (0, T) от j-го источника (j = 1,..., M), M — число источников в эксперименте. В формуле (3) uj (y, t; c, а) — волновое поле, полученное при решении прямой задачи (1)-(2) при заданных коэффициентах c(r), a(r). Скорость распространения волны v(r) находится из соотношения c(r) = 1/v2(r). Заметим, что функционал невязки представляет собой сумму по j = 1,..., M значений невязки, полученных для каждого источника. При каждом фиксированном источнике j интеграл суммируется по времени (0,T) и по окружности Г, где расположены все приемники, принимающие сигналы для выбранного источника j. Обратная задача ставится как задача поиска функций cg (r), ag (r), минимизирующих функционал невязки min Ф(и(с, a)) = Ф(и^, ag)). Функции cg(r), ag(r)
c(r),a(r)
принимаются за приближенное решение обратной задачи. Согласно теории решения обратных задач cg(r), ag (r) стремятся к точному решению обратной задачи при 6 ^ 0.
Эффективными методами минимизации функционала невязки Ф(^ a) являются градиентные методы. Выражение для градиента функционала невязки в строгой математической постановке получено в работах [10, 18]. Имея явное выражение для градиента, можно построить различные итеративные алгоритмы минимизации функционала невязки. Градиент функционала Ф(^ a) имеет вид
M
ФСМ)^ /
j=1 о
' (r, t)uj (r,t) dt,
M
(4)
] - (г, ¿К(г,*)ё*,
^ = ! 0
где и — решение «основной» задачи (1), (2), а — — решение «сопряженной» задачи (5), (6) для у — источника и заданных с(г), а(г).
«Сопряженная» задача имеет вид
Рис. 1. Схема эксперимента послойной томографии. Источники обозначены цифрой 1, приемники — 2
c(r)wtt(r,t) — a(r)wt(r,t) — Aw(r,t) = uj (r,t)|rer — Uj (y, t); w(r,t = T ) = 0, wt(r,t = T ) = 0.
(5)
w
На границе расчетной области ставится условие неотражения [19, 20]. Для расчетов использовались явные разностные схемы четвертого порядка точности [21].
Для решения «основной» задачи (1), (2) необходимо задавать генерируемый источником в точке г0 импульс д (£). Однако в настоящей работе использовался следующий подход. Экспериментальные данные и? (7, £) измерялись на всей границе Г с шагом менее половины центральной длины волны для каждого источника ^ в однородной среде L без объекта в. Далее решалось уравнение (1) распространения волны в обратном времени в однородной среде без объекта, начиная от конечного момента времени Т, с нулевыми конечными условиями, нулевой правой частью уравнения и граничными условиями и? (г, £)|геГ = и? (7,£) [22]. Полученные значения и(г, £ = т) и (г, £ = т) для некоторого малого момента времени т > 0 принимались за начальные значения для «основной» задачи.
При таком подходе для решения «основной» задачи (1), (2) вместо импульса д (£) задаются ненулевые начальные условия. Это позволяет более точно смоделировать зондирующий акустический импульс, неявно учесть сложные процессы, происходящие при преобразовании источником электрического импульса в акустический, учесть эффекты дисперсии, диаграммы направленности. Итерационные алгоритмы реконструкции скоростного разреза апробированы на суперЭВМ на решении большого количества модельных задач. Целью настоящей работы является апробация алгоритмов на экспериментальных данных.
2. СТЕНД ДЛЯ УЛЬТРАЗВУКОВЫХ ТОМОГРАФИЧЕСКИХ ИССЛЕДОВАНИЙ
Экспериментальные данные для ультразвуковых томографических исследований были получены на разработанном и изготовленном томографическом стенде. Экспериментальный стенд предназначен для проведения исследований по томографическому восстановлению акустических свойств объектов с акустическими параметрами, близкими к параметрам мягких тканей. Стенд обеспечивает прецизионное перемещение излучателя и приемника ультразвукового импульса вокруг исследуемого объекта и позволяет собирать данные, необходимые для восстановления томографических изображений как в послойном варианте, так и в трехмерной томографической схеме с данными на цилиндрической поверхности. Полученные на стенде экспериментальные данные используются для апробации методов акустической томографии.
Экспериментальный стенд включает цилиндрическую ванну, ультразвуковые излучатели, гидрофон, двухкоординатные (ротационный и линейный) приводы для излучателя и гидрофона, блок управления механическим узлом, генератор ультразвуковых импульсов, модуль АЦП. Конструкция экспериментального стенда изображена на рис. 2.
Внутри цилиндрической ванны 1, заполненной водой, закрепляется исследуемый объект (фантом) 2. Вокруг объекта перемещаются ультразвуковой излучатель 3 и гидрофон 4. Излучатель с помощью держателя крепится к линейному приводу 5. Гидрофон 4 с помощью держателя крепится к линейному приводу 6. В качестве линейных приводов были использованы
а С,
Рис. 2. Экспериментальный стенд для УЗ томографических исследований: а — 3Б-модель; б — общий вид стенда
линейные модули KK40 фирмы Hiwin с шаговыми двигателями. Линейные приводы обеспечивают перемещение излучателя и гидрофона по вертикали в пределах 14 см. Шаг перемещения линейных приводов составляет 5 мкм. Скорость перемещения линейных приводов в процессе сбора данных составляет 3 мм/c. Положение линейных приводов контролируется с помощью оптических концевых выключателей, входящих в комплект поставки линейных модулей KK40, с точностью 25 мкм.
Вращение излучателя и гидрофона вместе с линейными приводами вокруг исследуемого объекта осуществляется с помощью ротационных приводов 7, 8. Угол поворота контролируется с помощью прецизионных оптических датчиков RESA фирмы Renishaw 9, 10. Оптические датчики обеспечивают измерение угла поворота с точностью 5". Для управления механическими приводами экспериментального стенда с помощью персонального компьютера был разработан и изготовлен электронный блок управления на программируемых микроконтроллерах семейства AVR фирмы Atmel, подключаемый к персональному компьютеру через интерфейс USB.
Общий вид экспериментального стенда показан на рис. 2, б. Ванна 1 имеет форму цилиндра диаметром 500 мм и высотой 400 мм. По периметру ванны расположен поглощающий материал, предотвращающий переотражение ультразвука от металлических стенок ванны. Управление линейными приводами излучателя 3 и приемника 4 и ротационными приводами 7, 8 осуществляется с помощью блока управления 11. Акустический импульс для излучателя 3 формируется цифровым генератором 12 и усиливается амплитудным усилителем 13. Волновое поле на цилиндрической поверхности при фиксированном положении источника измеряется с помощью гидрофона 4. Измеренный сигнал усиливается предварительным усилителем 14, оцифровывается модулем АЦП 15 и поступает в память управляющего компьютера.
В экспериментах на стенде использовался гидрофон Teledyne Reson TC4038 (слева вверху) c предусилите-лем Teledyne Reson VP1000 (внизу) на рис. 3. Данный гидрофон имеет чувствительность ^15 мкВ/Па, круговую диаграмму направленности, рабочую полосу частот 10-800 кГц. Предусилитель Teledyne Reson VP1000 обеспечивает коэффициент усиления до 32 дБ
Рис. 3. Фотография приемника, источника и предусилителя
и полосу пропускания 1 МГц. Сферический пьезоке-рамический ультразвуковой излучатель (справа вверху на рис. 3), рассчитан на применение в диапазоне частот 50-600 кГц и обеспечивает акустическое давление не менее 3 Па-м/В.
Для формирования зондирующих импульсов используется генератор сигналов специальной формы АКТАКОМ AWG-4110 (12). Генератор сигналов обеспечивает с высокой точностью и стабильностью формирование сигнала произвольной формы в полосе частот до 10 МГц методом цифроаналогового преобразования. Сигнал, формируемый генератором, загружается в цифровом виде с компьютера. Электрический импульс, сформированный генератором, усиливается амплитудным усилителем АКТАКОМ АУА-1745 (13), который обеспечивает выходное напряжение ±85 В при токе до 500 мА для питания ультразвуковых излучателей. Полоса пропускания усилителя составляет 1 МГц.
Для оцифровки принятых гидрофоном акустических данных используется модуль АЦП Е20-10 (15) производства ЬСагё, обеспечивающий оцифровку сигналов
с частотой выборки до 10 МГц. Разрядность АЦП составляет 14 бит, пределы измерения — ±300 мВ, ±1 В или ±3 В, точность оцифровки в диапазоне ±300 мВ составляет ^50 мкВ.
Функциональная схема стенда приведена на рис. 4. В акустическом тракте экспериментального стенда электрический импульс формируется цифровым генератором по сигналу, полученному от блока управления. Электрический импульс усиливается амплитудным усилителем и поступает на пьезоэлектрический ультразвуковой излучатель. Прошедший через объект акустический импульс принимается гидрофоном, усиливается предусилителем и поступает на модуль АЦП для оцифровки. Оцифрованные данные поступают в память управляющего компьютера. Объем данных для расчета скоростного разреза в одном слое составляет ^100 Мбайт. Для реконструкции скоростного разреза использовался графический процессор NVidia Tesla K40s суперкомпьютера «Ломоносов-2» в режиме удаленного доступа.
Для автоматизации экспериментов по ультразвуковой томографии было разработано программное обеспечение, обеспечивающее работу экспериментального стенда под управлением персонального компьютера. Для функционирования программы необходима операционная система Windows XP или выше, библиотека lusbapi.dll поддержки модулей АЦП производства L-Card и драйвер связи с микроконтроллерным блоком ftd2xx фирмы Future Technology Devices. Программа реализована на языке С++. Микропрограммы для блока управления реализованы на AVR Assembler.
Программа позволяет управлять движением линейных и ротационных приводов стенда посредством блока управления, получать информацию о точном положении приводов, запускать внешний генератор сигналов в заданные моменты времени, контролировать текущее состояние экспериментального стенда с помощью графического интерфейса и управлять процессом сбора данных.
Как показало тестирование, наиболее оптимальным является режим непрерывного сбора данных. В этом режиме акустические сигналы регистрируются в точках поверхности цилиндра, координаты которых снимаются с высокой точностью с датчиков углового положения ротационных приводов. В процессе сбора данных гидрофон непрерывно поворачивается вокруг исследуемого объекта.
Управляющий компьютер
Суперкомпьютер ■= >
«Ломоносов-2»
Генератор импульсов
>
Излучатель Объект Гидрофон
Амплитудный усилитель
А)
АЦП
Предусилитель
Двухкоординатные приводы излучателя и гидрофона
I
I
Блок управления
<
Рис. 4. Функциональная схема стенда
u(t), мВ 200
100
0
-100 -200
A(F) 1
u(t), мВ 600
400 200 0
-200
A(F)
0
70 80 90 100 мкс 0 200 400 600 800 кГц по 120 130 140 мкс "0 200 400 600 кГц Рис. 5. Акустический импульс и его спектр для разных задающих импульсов
За один оборот гидрофона снимается 550 пакетов данных, соответствующих положениям ротационного привода гидрофона с шагом 36'. Каждый пакет данных содержит 4 серии из 2048 16-битных отсчетов. Частота выборки отсчетов составляет 10 МГц. Серии данных снимаются с интервалом 2 мс и усредняются для повышения точности измерений. Интервал между сериями выбран исходя из того, чтобы переотражения ультразвука от элементов конструкции успели рассеяться. Для сбора данных на одном сечении требуется ~10 мин.
3. ПРОБЛЕМА ФОРМИРОВАНИЯ ЗОНДИРУЮЩИХ ИМПУЛЬСОВ В УЛЬТРАЗВУКОВОЙ ТОМОГРАФИИ
Одна из основных проблем решения обратных задач ультразвуковой томографии — нелинейность обратной задачи реконструкции скоростного разреза. Эффективность алгоритмов решения обратной задачи во многом зависит от формы зондирующего импульса. Показано, что оптимальным вариантом является короткий широкополосный импульс с минимальным количеством переколебаний. На рис. 5, а, б приведена форма акустического импульса и его спектр, если для возбуждения пьезокерамического излучателя использован короткий импульс длительностью 0.2 мкс. Видно, что акустический импульс имеет большое количество переколебаний, связанных с собственными резонансными частотами пьезокерамического элемента. Для того чтобы скорректировать форму импульса, необходимо решать обратную задачу. Пусть f (t) — возбуждающий электрический импульс, а u(t) — акустический импульс на выходе излучателя. В линейной модели связь между u(t) и f (t) описывается интегральным уравнением типа свертки [23]. Ядро уравнения измеряется экспериментально, если на вход акустического тракта подать возбуждающий импульс короче 1 мкс. С формальной точки зрения задача корректировки акустического импульса является некорректно поставленной. Методы решения таких задач разработаны в конце прошлого столетия. При использовании одной из стандартных программ монографии [24] удается скорректировать форму акустического импульса.
На рис. 5, в, г приведена форма скорректированного акустического импульса u(t) и его спектр. Как видно из рис. 5, в, короткий широкополосный сигнал u(t) содержит малое количество переколебаний. Сформированный таким образом зондирующий импульс является стабильным. Многочисленные эксперименты показали, что акустические импульсы стабильно воспроизводятся с точностью не хуже 1% от максимальной амплитуды сигнала. Именно такие импульсы были использованы в эксперименте. С физической точки зрения использование методов решения обратных задач позволило рассчитать такой электрический импульс f (t),
который обеспечивает подъемы на низких и высоких частотах, чтобы скомпенсировать ослабление эффективности ультразвукового излучателя на этих частотах.
4. РЕКОНСТРУКЦИЯ СКОРОСТНОГО РАЗРЕЗА ПО ЭКСПЕРИМЕНТАЛЬНЫМ ДАННЫМ
В экспериментах по реконструкции ультразвуковых томографических изображений использовался фантом, помещенный в воду. Фантом имеет акустические параметры близкие к параметрам мягких тканей человека. Скорость распространения звука в мягких тканях варьируется от 1400 м/с до 1550 м/с. Наименьшее значение — в жировой ткани (1400-1450 м/с), в железистой ткани — 1500-1550 м/с. Большая часть раковых опухолей характеризуется повышенной скоростью распространения ультразвука — 1600 м/с и выше [25].
На рис. 6 приведена фотография и 3Э-модель фантома, который представляет собой силиконовый цилиндр диаметром 62 мм и длиной 140 мм, содержащий три цилиндрические полости. Первая полость диаметром 10 мм расположена наклонно и имеет высоту 140 мм. Вторая полость начинается у основания фантома и имеет диаметр 15 мм и высоту 45 мм. Третья полость диаметром 15 мм и высотой 45 мм начинается от верха фантома. В нее помещен однородный пластиковый стержень диаметром 7 мм. Фантом прикрепляется к основанию с помощью трех металлических игл диаметром 1 мм. Скорость звука в силиконе составляет 1440 м/с; в окружающей воде, заполняющей полости, — 1500 м/с, в пластиковом стержне ^ 1800 м/с.
Для реконструкции скоростного разреза использовались данные при 24 положениях излучателя и 500 положениях приемника в каждом сечении. На рис. 7 приведены восстановленные по экспериментальным данным скоростные разрезы фантома в трех различных
Рис. 6. Фантом: а — 3Б-модель; б — фотография
20
40
60
80 мм 20
20 40 60 80
40
60 80
мм 20
40
60 80
мм 20
40 60
1.8 1.7 1.6 — 1.5 1.4
80 v, км/с
Рис. 7. Реконструированные сечения скоростного разреза
— 1.60
1.58
1.56
1.54
1.52
- 1.50
J 1.48
я 1.46
и 1.44
и 1.42
® 1.40
— 1.60
1.58
1.56
1.54
1 1.52
■ 1.50
U 1.48
11.46
1 1.44
11.42
™ 1.40
0.050 0.045 0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005 0.000
Рис. 8. Реконструированные сечения: а-б скоростного разреза; в — коэффициента поглощения
горизонтальных сечениях: г1 = 25 мм, г2 = 70 мм, г3 = 115 мм. На рис. 7, а представлено сечение при г1 = 25 мм, которое пересекает нижнюю и наклонную полости. Видно, что восстановленное значение скорости ультразвука в полостях совпадает со скоростью в воде. На рис. 7, б представлено среднее сечение при г2 = 70 мм, которое содержит только наклонную полость и не содержит нижнюю и верхнюю полости 3 и 4, поскольку проходит между ними. На рис. 7, в представлено сечение при ¿3 = 115 мм, которое пересекает наклонную и верхнюю полости. В верхней полости фантома видно изображение пластиковый стержень с более высокой скоростью волны. Из рис. 7, в видно, что хорошо восстанавливаются не только форма неоднородных областей, но и скорость звука внутри них. Реконструированное значение скорости в стержне составило 1800 м/с. Видно, что наклонная полость меняет свое положение на сечениях с разными г. На сечениях хорошо видны металлические иглы крепления диаметром ~1 мм.
Максимальное различие скорости звука в неоднород-ностях фантома составляет ±15%. Тем не менее, как видно из рис. 7, разработанные алгоритмы обеспечивают высокое разрешение по скорости звука, несмотря на низкий контраст.
На рис. 8, а, б приведены восстановленные по экспериментальным данным скоростные разрезы фантома в двух горизонтальных сечениях. Фантом представляет собой однородный силиконовый цилиндр, на наружной стороне которого сделана канавка вдоль образующей цилиндра шириной 2 мм и глубиной 1 мм (в сечении — половина окружности радиуса 1 мм). Для крепления цилиндра использовалась одна металлическая спица. Видно, что, несмотря на малую разницу (меньше 5%) скоростей волн в воде и в цилиндре, надежно восстанавливается форма неоднородности. Рис. 8 демонстрирует, что пространственное разрешение не хуже 2 мм при длине волны ~4 мм.
Реконструкция скоростного разреза проводилась в модели с поглощением (1). На рис. 8, в приведено распределение коэффициента поглощения a(r) в сечении на рис. 8, а. Более светлые участки соответствуют более сильному поглощению. Коэффициент поглощения восстанавливается хуже, чем скоростной разрез [13]. Аномальные значения скорости и поглощения в металлических иглах крепления связаны с ограничениями используемых методов, поскольку акустические параметры в металле очень сильно отличаются от акустических параметров в мягких тканях.
Объем собираемых данных для одного сечения в проведенных экспериментах составил ~100 МБайт. Длительность регистрации сигнала — 200 мкс. Шаг расчетной сетки по пространственным координатам — 0.33 мм, по времени — 0.15 мкс. Количество итераций градиентного спуска — 30. Для реконструкции томографических изображений использовался графический процессор NVidia Tesla K40s в составе суперкомпьютера «Ломоносов-2» МГУ [26]. Использование одной карты NVidia Tesla K40s позволяет в 30 раз ускорить реализацию разработанных авторами алгоритмов [27]. Время реконструкции одного слоя скоростного разреза составило порядка 5 мин.
ЗАКЛЮЧЕНИЕ
1. Разработанные в работах [27-31] алгоритмы решения обратных задач ультразвуковой томографии прошли апробацию на экспериментальных данных. Для тестирования алгоритмов на экспериментальных данных был разработан и изготовлен стенд. Для зондирования использовались низкочастотные широкополосные источники с рабочими частотами менее 600 кГц.
2. Реконструкция скоростного разреза осуществлялась в рамках послойной модели. Показано, что пространственная разрешающая способность в горизонтальной плоскости составляет ~2 мм. Время расчетов
составило ^S мин при использовании одного графического процессора NVidia Tesla K40s для расчета одного слоя. Важным результатом работы является экспериментальное подтверждение адекватности модели физическим процессам.
3. Проведенные эксперименты показали, что послойная модель может быть использована для реконструкции 3D скоростного разреза. Однако более адекватной моделью является ßD-модель, в которой обратная задача реконструкции ßD-изображения решается не в послойном варианте, а сразу по всем данным, собранным на цилиндрической поверхности. В этом случае в обратной задаче неизвестным является функция c(r) во всех точках трехмерной сетки. Эта задача является намного более сложной. В работах [32] разработаны алгоритмы реконструкции 3D скоростного разреза в такой постановке. Алгоритмы апробированы на модельных расчетах. Конструкция стенда позволяет сформировать данные для реконструкции 3D скоростного разреза. Экспериментальная проверка на стенде 3D-алгоритмов является предметом дальнейших исследований авторов статьи.
Исследование выполнено за счет гранта Российского научного фонда (проект № 17—11—0106s). Исследование выполнено в Московском государственном университете им. М. В. Ломоносова. Работа выполнена с использованием ресурсов суперкомпьютерного комплекса МГУ им. М. В. Ломоносова.
СПИСОК ЛИТЕРАТУРЫ
1. Goncharsky A., Romanov S., Seryozhnikov S. // Ultrasonics.
2016. 67. P. 136.
2. Буров В. А., Зотов Д. И., Румянцева О. Д. // Акустический журнал. 2014. 60, №4. С. 443.
3. Буров В. А., Зотов Д. И., Румянцева О. Д. // Акустический журнал. 201S. 61, №2. С. 2S4.
4. Jifik R., Peterlik I., Ruiter N. et al. // IEEE Trans. on Ultrasonics, Ferroelectrics, and Frequency Control. 2012. 59, N 2. P. 2S4.
5. Ruiter N. V., Zapf M., Hopp T. et al. // In: Proceedings, Medical Imaging 2017: Ultrasonic Imaging and Tomography.
2017. 10139. 101391N.
6. Schmidt S, Duric N, Li C. et al. // Medical Physics. 2011. 38, N 2. P. 998.
7. Wiskin J. W., Borup D. T., Iuanow E. et al. // EEE Trans. Ultrason. Ferroelectr. Freq. Control. 2017. 64, N 8. P. 1161.
8. Natterer F. // In: AMS: Tomography and Inverse Transport Theory. 2011. 559. P. 1S1.
9. Natterer F. // In: Handbook of Mathematical Methods in Imaging. Springer Nature. 2014. P. 1.
10. Beilina L., Klibanov M.V., Kokurin M.Yu. // J. of Math. Sci. 2010. 167, N 3. P. 279.
11. Kuzhuget A., Beilina L., Klibanov M. et al. // Inverse Problems. 2012. 28, N 9. 095007.
12. Goncharsky A. V., Romanov S. Y. // Inverse Problems. 2013. 29, N 7. P. 075004.
13. Goncharsky A. V., Romanov S. Y. // Physics in Medicine and Biology. 2014. 59, N 8. P. 1979.
14. Boehm C., Martiartu N. K, Vinard N. et al. // In: Proc. SPIE. Medical Imaging 2018: Ultrasonic Imaging and Tomography. 2018. 10580. 105800H.
15. Yong P., Liao W., Huang J., Li Z. // Inverse Problems. 2018. 34, N 4. 045006.
16. Goncharsky A. V., Romanov S. Y., Seryozhnikov S. Y. // Dok-lady Physics. 2016. 61, N 5. P. 211.
17. Lenox M.W., Wiskin J., Lewis M.A. et al. // Int. J. of Biomed. Imaging. 2015. 2015. ID 454028.
18. Goncharsky A. V., Romanov S. Y. // Inverse Problems. 2017. 33, N 2. 025003.
19. Engquist B., Majda A. // Mathematics of Computation. 1977. 31, N 139. P. 629.
20. GivoliD., Keller J. B. //Wave Motion. 1990. 12. 261279.
21. Bilbao S. // Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics. Wiley, Chichester. Wiley, Chichester, 2009.
22. FinkM. // Contemporary Physics. 1996. 37, N 2. P. 95.
23. Гончарский А. В., Романов С.Ю., Серёжников С.Ю. // Вычислительные методы и программирование: Новые вычислительные технологии. 2018. 19. С. 150.
24. Tikhonov A. N., Goncharsky A. V., Stepanov V. V., Yago-la A. G. // Numerical Methods for the Solution of Ill-Posed Problems. Springer Netherlands, 1995.
25. Hendee W. R., Ritenour E. R. // Medical Imaging Physics. Wiley-Blackwell, 2002.
26. Sadovnichy V., Tikhonravov A., Voevodin V., Opanasenko V. // In: Contemporary High Performance Computing: From Peta-scale toward Exascale. Chapman & Hall/CRC. Computational Science, Boca Raton, USA, CRC Press 2013. 283307.
27. Goncharsky A., Seryozhnikov S. // Supercomputing. RuSCDays 2017. Communications in Computer and Information Science. Springer, Cham. 2017. 793. P. 363.
28. Romanov S. // Supercomputing. RuSCDays 2017. Communications in Computer and Information Science. Springer, Cham. 2017. 793. P. 67.
29. Гончарский А. В., Романов С. Ю., Серёжников С. Ю. // Вычислительные методы и программирование: Новые вычислительные технологии. 2017. 18. С. 267.
30. Гончарский А. В., Романов С. Ю., Серёжников С. Ю. // Вычислительные методы и программирование: Новые вычислительные технологии. 2017. 18. С. 312.
31. Bazulin E., Goncharsky A., Romanov S., Seryozhnikov S. // Lobachevskii Journal of Mathematics. 2018. 39, N 4. P. 486.
32. Goncharsky A. V., Romanov S. Y., Seryozhnikov S. Y. // Wave Motion. 2014. 51, N 3. P. 389.
Low Frequency Ultrasonic Tomography: Mathematical Methods and Experimental Results A.V. Goncharsky, S.Yu. Romanova, S.Yu. Seryozhnikov
Research Computing Center, Lomonosov Moscow State University, Moscow 119234, Russia E-mail: a [email protected].
This paper describes investigation of the possibilities of low-frequency ultrasonic tomography in medicine. The main application is the development of tomographic methods for differential diagnosis of breast cancer. To solve the inverse problem, a mathematical model is used that describes both the phenomena of diffraction and refraction and the absorption of ultrasound in an inhomogeneous medium. The algorithms developed for reconstructing the velocity structure were evaluated using the test bench for experimental studies. Broadband ultrasonic sources in the 50-600 kHz range were used for sounding. A Teledyne Reson TC4038 hydrophone with a sensitivity of approximately 4 ^V/Pa and a frequency range of 10-800 kHz was used as a receiver. A Teledyne Reson VP1000 preamplifier was used to amplify the signals.
The studies were carried out on phantoms with acoustic parameters close to the properties of human soft tissues. A layerwise model was used to reconstruct the three-dimensional velocity structure. The sound speed cross-sections reconstructed from the experimental data have a resolution of approximately 2 mm, while the central wavelength is approximately 4 mm. An important result of the work is an experimental confirmation of the correspondence of the model to physical processes.
Keywords: ultrasonic tomography, wave equation, inverse problems, hydrophone, supercomputer technologies. PACS: 02.30.Zz, 02.60.Cb, 43.35.Wa. Received 18 May 2018.
English version: Moscow University Physics Bulletin. 2019. 74, No. 1. Pp. 43-51.
Сведения об авторах
1. Гончарский Александр Владимирович — доктор физ.-мат. наук, профессор, зав. лабораторией; тел.: (495) 939-27-59, e-mail: [email protected].
2. Романов Сергей Юрьевич — доктор физ.-мат. наук, вед. науч. сотрудник; тел.: (495) 939-27-59, e-mail: [email protected].
3. Серёжников Сергей Юрьевич — канд. физ.-мат. наук, сотрудник; тел.: (495) 939-27-59, e-mail: [email protected].
