Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org



         Исследование томографических схем низкочастотной
         ультразвуковой диагностики на суперкомпьютерах*
                                                         С.Ю. Романов
                 Московский государственный университет им. М.В. Ломоносова

            Статья посвящена сравнению различных схем ультразвуковой (УЗ) томографии. Об-
            ратные задачи УЗ томографии рассматриваются как коэффициентные обратные зада-
            чи для волнового уравнения. Алгоритмы решения задачи УЗ томографии основаны
            на прямом вычислении градиента функционала невязки через решение сопряжённой
            задачи. Проведён сравнительный анализ возможностей различных томографических
            схем на прохождение, на отражение в трёхмерном случае. Полученные результаты
            сравниваются с результатами моделирования обычных УЗИ аппаратов, в которых
            изображение получается с одного ракурса. Для расчета прямых и обратных задач УЗ
            томографии использовался суперкомпьютер на GPU.


1. Введение
    Целью работы является исследование различных схем волновой томографии, используе-
мых для реконструкции внутренней структуры объектов. Результаты работы могут найти при-
менение в таких областях как ультразвуковая (УЗ) диагностика в медицине [1-4], неразрушаю-
щий контроль в промышленности, инженерной сейсмике, сейсморазведке, электромагнитном
зондировании приповерхностных слоев Земли и т.д. [5-11].
    Настоящая работа посвящена в первую очередь проблемам УЗ томографии рака молочной
железы. Популярные в медицине УЗИ аппараты работают по схеме на отражение, когда источ-
ники и приёмники располагаются с одной стороны исследуемого объекта и улавливается отра-
жённое от объекта излучение. Важнейшим вопросом проектирования УЗ томографов является
выбор наиболее эффективной схемы томографических исследований. В рентгеновских томо-
графических комплексах используются данные, полученные с разных ракурсов. По аналогии с
рентгеновской томографией представляется естественным ожидать, что использование много-
ракурсных УЗ томографических схем позволит более детально восстанавливать структуру не-
однородного объекта по сравнению с обычными одноракурсными УЗИ аппаратами. С помо-
щью математического моделирования исследованы возможности различных схем УЗ томогра-
фии — полноракурсных, на прохождение, на отражение, с использованием нескольких поло-
жений источников, окружающих объект с разных сторон.

2. Постановка прямой и обратной задачи ультразвуковой томографии
      В работе обратная задача рассматривается в скалярном волновом приближении. Давление
акустического поля u(r,t) в области   R3 , создаваемое точечным источником, находящимся
в точке r0 и генерирующим короткий широкополосный импульс, который описывается
функцией f(t), удовлетворяет волновому уравнению:

         c(r )utt (r, t )  u(r, t )   (r  r0 )  f (t ) ,                  (1)

       u(r, t  0)  ut (r, t  0)  0 ,  n u | ΓT  p( r, t ) ,               (2)

где t – время, 0<t<T, c(r)=1/v2(r) , v(r) - скорость звука в среде,  – оператор Лапласа по
r  R3 , Γ – граница области  ,  n u |ΓT – производная вдоль нормали к поверхности Γ в об-


*
    Исследования поддержаны РФФИ (грант 14-07-00078 а).


                                                                 199
   Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org


ласти Γ  0, T  , p(r,t) – некоторая известная функция. Скалярная волновая модель позволяет
описать дифракцию, рефракцию, переотражение волн в среде.
    Рассмотрим обратную по отношению к задаче (1–2) задачу, которая состоит в нахождении
неизвестной скорости звука v(r) по экспериментальным данным измерения акустического дав-
ления U  , t  на границе Γ области  за время (0,T) при различных положениях источника r0.
Обратная задача ставится как задача минимизации по c(r) функционала невязки:


                              u(  , t )  U (  , t ) ddt .
                        T
                  1
      ( u( с ))  
                                                       2
                                                                             (3)
                  20         Γ


Здесь u( r, t ) – решение прямой задачи (1–2) при заданной c(r)=1/v2(r).
    Обратные задачи томографии в волновой модели являются нелинейными, некорректно по-
ставленными задачами [12]. Минимизация функционала невязки в этом случае является слож-
ной проблемой [13]. Прорывные результаты в области решения задач волновой томографии
связаны с возможностью прямого вычисления градиента функционала uc  по формуле (4)

                        T
      C ( u( c ))   wt ( r , t )ut ( r , t )dt ,                         (4)
                        0


где u(r,t) – решение задачи (1–3), а w(r,t) – решение «сопряженной» задачи. Зная C из (4),
можно построить различные итеративные схемы для минимизации функционала невязки (3).
Для решения трёхмерной обратной задачи была использована явная разностная схема для вол-
нового уравнения второго порядок точности. Для данной постановки задачи эти результаты
опубликованы в [4, 14]. В похожих постановках выражение для градиента получено в [15, 16].

3. Возможности GPU-кластеров для решения задач УЗ томографии
    Обратные задачи УЗ томографии являются очень сложными с вычислительной точки зре-
ния. Для получения достаточно высокого разрешения необходимо решать 3D задачу восста-
новления трехмерной функции с(x,y,z) на сетке размерностью до 5003 точек. Таким образом,
необходимо решать нелинейную обратную задачу с общим числом неизвестных более 100 млн.
    Предлагаемый в работе численный метод с явной разностной схемой относится к SIMD ал-
горитмам (Single instruction – multiple data) и имеет высокую степень параллелизма. Такие зада-
чи могут быть эффективно реализованы на вычислительных системах с параллельной архитек-
турой. В работах [14, 17, 18] приведены результаты эффективного распараллеливания 2D зада-
чи на кластерной архитектуре с использованием процессоров общего назначения, при исполь-
зовании до десятков тысяч процессорных ядер. Графические процессоры (GPU) также позво-
ляют получить очень высокую производительность на задачах такого типа. Вычислительная
конфигурация с использованием GPU необходимая для решения задач УЗ томографии может
быть реализована в рамках одной стойки, что делает возможным установку небольших, но вы-
сокопроизводительных суперкомпьютеров непосредственно в медицинских учреждениях и ре-
шение задач УЗ томографии на практике за приемлемое для медицинских исследований время.
В настоящее время GPU успешно используются в задачах УЗ томографии [19].
    За счёт наличия высокопроизводительной памяти современные GPU позволяют получить
существенный выигрыш на data-intensive задачах, которые содержат относительно немного
арифметических операций и много операций обращения в память. Именно к таким задачам от-
носится явная разностная схема, используемая при решении обратной задачи УЗ томографии.
    Несмотря на огромные размеры массивов, заданных на сетках размерностью ~n4 для 3D
задачи, разработанный алгоритм решения таков, что хранить в памяти необходимо только дан-
ные для двух шагов по времени явной разностной схемы, размер которых ~ n3 (порядка 0.1% от
числа вычисляемых значений). Это позволяет разместить их полностью в высокопроизводи-


                                                           200
   Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org


тельной памяти GPU. Во внешней памяти хранятся данные только для границ. Ввиду малого
объема граничных данных, время обмена данными по медленному каналу с внешней памятью
составляет около 2% по сравнению со временем расчетов нового слоя по времени. В результате
применение GPU в такой задаче только за счет более скоростной памяти даёт ускорение расчё-
тов по сравнению с однопоточными расчетами на обычном процессоре до 30 раз.
    В таблице 1 приведены времена расчетов на 1 GPU NVIDIA Tesla X2070 суперкомпьютера
«Ломоносов» СКЦ МГУ для 4 источников 20 итераций для размеров сетки 500×500 и
1000×1000 точек в 2D задаче УЗ томографии. Ниже приведены для сравнения результаты ана-
логичных расчетов на процессорах общего назначения суперкомпьютера «Ломоносов» для раз-
личных разбиений области расчетов на подобласти (число источников × разбиение по X × раз-
биение по Y). Как видно, время работы одного процессора GPU совпадает со временем расчета
на ~ 25-35 ядрах для размера сетки 5002 точек и на ~ 13-20 ядрах для размера сетки 10002 точек.
                      Таблица 1. Время выполнения 20 итераций программы
                      Размер сетки расчетов                   500×500     1000×1000
                         Время, сек. (GPU)                       98          695
                                     Число процессов
                                        16(=4×2×2)              214
            Время, сек. (CPU)          144(=4×6×6)               17            94
                                      400(=4×10×10)                            22
                                      784(=4×14×14)                            14
    Распараллеливание вычислений выполняется следующим образом. При решении обратной
задачи УЗ томографии используются данные, полученные при различных положениях источ-
ников. Суммарный градиент функционала невязки метода градиентного спуска на каждой ите-
рации равен сумме градиентов grad_i, полученных при решении подзадачи для каждого источ-
ника i. Заметим, что основной объем вычислений суммарного градиента выполняется незави-
симо в каждой подзадаче по вычислению grad_i. Это говорит о том, что структура рассматри-
ваемого алгоритма такова, что имеется возможность эффективного распараллеливания вычис-
лений по источникам. Предлагаемая структура вычислений предполагает, что каждая графиче-
ская карта проводит расчеты для одного источника (одной подзадачи). Обмены между картами
происходят только после каждой итерации для суммирования градиентов каждой подзадачи.
Это позволяет распараллелить задачу на число узлов (GPU), равное числу положений источни-
ков, практически без потерь эффективности.
    Количество положений источников может достигать нескольких сотен в реальных задачах
томографии в медицине. В этом случае теоретически можно использовать несколько сотен карт
GPU с высокой эффективностью. При недостаточном количестве GPU в системе на каждом
GPU могут выполняться последовательно подзадачи для нескольких источников.
    При проведении вычислений в 3D по явной разностной схеме трёхмерный объём данных
для текущего момента времени обрабатывается последовательно по оси Z, это так называемый
“z-marching” метод. При проведении расчетов вся область размера порядка 400×400 в плоско-
сти ХУ разбивается на прямоугольные блоки размера 32×16 точек сетки. Выбор размера блока
определяется используемым типом GPU. Каждый мультипроцессор GPU обрабатывает отдель-
ный блок 32×16×Zn точек сетки последовательно по слоям, где Zn – размер расчётной области
вдоль оси Z. При последовательных послойных расчетах автоматически активно используются
регистры и кэш-память GPU.

4. Модельные расчёты трёхмерных задач ультразвуковой томографии
    Вопросам моделирования в обратных задачах УЗ томографии посвящены публикации [20-
22]. Проведенные в настоящей работе модельные расчёты позволяют оценить возможности
различных томографических схем в 3D варианте. Для сокращения времени расчётов размер
исследуемого объекта был уменьшен до ≈60 мм. Расчётная область представляет собой куб
размером H×H×H, где H=100 мм. При шаге разностной сетки 0,3 мм количество точек сетки
составляет 320×320×320. Схема трёхмерного вычислительного эксперимента показана на


                                               201
   Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org


Рис. 1 слева. В расчётах использовалось 24 положения источника S в двух плоскостях
z=h0=0.33H и z=h1=0.66H. Длина волны зондирующего импульса λ равна 5 мм, диапазон значе-
ний скорости звука внутри фантома 1,43-1,6 км/c, скорость в окружающей среде — 1,5 км/c.
     В томографической схеме с полными данными приёмники располагались на всех гранях
куба с шагом 2 мм. Каждая грань куба, таким образом, является массивом из 54×54 детекторов.
В модельной задаче на прохождение, когда регистрируется только прошедшее через объект
излучение, приёмники R располагались только на противоположной от источников S грани ку-
ба, как показано на Рис. 1 в центре. Для каждой из боковых граней куба данные с детекторов
снимаются для 6 положений источников на противоположной грани. Эксперимент повторяется
для всех 4 боковых граней куба, таким образом, всего используется 24 положения источника. В
модельной задаче на отражение, когда регистрируется только отраженное от объекта излуче-
ние, приёмники R располагаются на той же грани куба, что и источник S, как показано на
Рис. 1 справа. Подобная схема используется в существующих УЗИ аппаратах.




                                 Рис. 1. Схемы 3D эксперимента
    Сечение восстановленного трёхмерного распределения скорости звука в плоскости z=const
показаны на Рис. 2. Использование суперкомпьютера позволило решать обратную задачу за
время ~ 2 часов. Количество задействованных GPU-устройств соответствовало количеству ис-
точников ультразвука. На Рис. 2 слева результаты без внесения шума в данные, на Рис. 2 в цен-
тре при уровне шума 2.5%, на Рис. 2 справа при уровне шума 10%. Видно, что в схеме с пол-
ными данными хорошо восстанавливаются даже самые мелкие включения диаметра ~ 1 мм.




            Рис. 2. Сечения восстановленного 3D изображения по схеме с полными данными
    Восстановленные по схеме на прохождение изображения в плоскости z=const показаны на
Рис. 3. На Рис. 3 слева без внесения шума в данные, на Рис. 3 в центре при уровне шума 2.5%,
на Рис. 3 справа при уровне шума 10%.




                                               202
    Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org




             Рис. 3. Сечения восстановленного 3D изображения по схеме на прохождение
    Как видно из сравнения Рис. 3 и 4, качество восстановленных изображений в схеме на про-
хождение хуже, чем в схеме с полным диапазоном углов. Мелкие детали изображения особенно
для зашумленных данных уже не различаются.
    Восстановленное по схеме на отражение изображение без шума в данных показано на
Рис. 4 слева. Модельные расчёты показывают, что возможности томографических исследова-
ний 3D объектов в рамках схемы на отражение весьма ограничены.




                 Рис. 4 Сечения восстановленного 3D изображения по схеме на отражение
    Таким образом, как видно из проведённых модельных расчётов, томографические схемы с
полным диапазоном углов обладают безусловными преимуществами по сравнению со схемами
на прохождение или на отражение. В схеме с полным диапазоном данных хорошо восстанавли-
вается не только форма неоднородных включений, но и пространственное распределение ско-
рости звука. В схеме на отражение можно восстанавливать только границы неоднородностей,
которые, как правило, сдвинуты вследствие неоднородной скорости. Восстановление абсолют-
ного значения скорости невозможно. В задаче диагностики рака молочной железы это означает,
что видимые на отражение неоднородности крайне трудно характеризовать, т.е. выяснить яв-
ляются ли они злокачественной опухолью.
    Проблема восстановления изображений на отражение является характерной не только для
медицинских УЗ исследований. Аналогичная проблема возникает, например, при электромаг-
нитном зондировании приповерхностных слоёв Земли, при поиске мин, в сейсморазведке, ин-
женерной сейсмике. Во всех этих задачах используется диагностика на отражение.
    Ситуация становится ещё хуже, если скорость звука внутри объекта варьируется в широких
пределах. При увеличении диапазона изменения скорости звука в 2 раза в схеме на отражение
получается изображение, показанное на Рис. 4 справа. Как видно, увеличение вариации скоро-
сти внутри объекта приводит к появлению артефактов, таких как удвоение границ и т.п.

Литература
1. Wiskin J., Borup D.T., Johnson S.A., Berggren M. Non-linear inverse scattering: high resolution
   quantitative breast tissue tomography // J. Acoust. Soc. Am. 2012. Vol. 131, No. 5. P. 3802–3813.
2. Pratt R.G., Huang L., Duric N., Littrup P. Sound-speed and attenuation imaging of breast tissue
   using waveform tomography of transmission ultrasound data // Medical Imaging 2007: Physics of
   Medical Imaging, Proceedings of the SPIE. 2007. No. 6510. P. 65104S.



                                                203
    Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org


3. Буров В.А., Зотов Д.И., Румянцева О.Д. Восстановление пространственных распределений
   скорости звука и поглощения в фантомах мягких биотканей по экспериментальным данным
   ультразвукового томографирования // Акуст. журн. 2015. Т. 61, № 2. С. 254–273.
4. Goncharsky A.V., Romanov S.Y. Supercomputer technologies in inverse problems of ultrasound
   tomography // Inverse Problems. 2013. Vol. 29, No. 7. P. 075004.
5. Батраков Д.О., Батракова А.Г., Головин Д.В., Кравченко О.В., Почанин Г.П. Определение
   толщин слоев дорожной одежды методом георадиолокационного зондирования // Физиче-
   ские основы приборостроения. 2014. Т. 3, № 2(11). С. 46–57.
6. Гайкович П.К., Хилько А.И., Гайкович К.П. Метод многочастотной ближнепольной аку-
   стической томографии объёмных неоднородностей морского дна // Известия высших учеб-
   ных заведений. Радиофизика. 2011. Т. 54, № 6. С. 431–443.
7. Кругликов А.А., Лазоренко Г.И., Шаповалов В.Л., Хакиев З.Б., Явна В.А. Компьютерное
   моделирование распространения акустических колебаний в насыпях железных дорог //
   Вестник Ростовского гос. университета путей сообщения. 2012. № 3. С. 135–140.
8. Данилин А.Н., Ерохин Г.Н., Кремлев А.Н., Пестов Л.Н., Саломонс Б., Круд Ф. Численное
   решение задачи определения сверхслабых дифракторов в сложной акустической среде //
   Вестник Балтийского федерального университета им. И. Канта. 2014. № 10. С. 115–119.
9. Головина С.Г., Романов С.Ю., Степанов В.В. Об одной обратной задаче сейсмики // Вест-
   ник Московского университета. Сер. 15: Вычисл. матем. и киберн. 1994. № 4. С. 16–21.
10. Базулин Е.Г. О возможности использования в ультразвуковом неразрушающем контроле
    метода максимальной энтропии для получения изображения рассеивателей по набору эхо-
    сигналов // Акустический журнал. 2013. Т. 59, № 2. С. 235–254.
11. Пикалов В.В., Баландин А.Л., Родионов Д.Г., Власенко М.Г., Князев Б.А. Терагерцовая то-
    мография низкоконтрастных объектов: алгоритмы и экспериментальные измерения // Вест-
    ник Новосибирского государственного университ. Серия: Физика. 2010. Т. 5, № 4. С. 91-97.
12. Бакушинский А.Б., Козлов А.И., Кокурин М.Ю. Об одной обратной задаче для трехмерного
    волнового уравнения // Журн. вычисл. матем. и матем. физики. 2003. Т. 43, № 8. С. 1201-
    1209.
13. Сулимов В.Д., Шкапов П.М. Гибридные алгоритмы вычислительной диагностики гидроме-
    ханических систем // Вестник Московского государственного технического университета
    им. Н.Э. Баумана. Серия: естественные науки. 2014. № 4 (55). С. 47-63.
14. Goncharsky A.V., Romanov S.Y. Inverse problems of ultrasound tomography in models with at-
    tenuation // Physics in medicine and biology. 2014. Vol. 59, No. 8. P. 1979–2004.
15. Natterer F. Possibilities and limitations of time domain wave equation imaging // Contemporary
    Mathematics. 2011. Vol. 559. P. 151–162.
16. Beilina L., Klibanov M.V., Kokurin M.Yu. Adaptivity with relaxation for ill-posed problems and
    global convergence for a coefficient inverse problem // Journal of Mathematical Sciences. 2010.
    Vol. 167, No. 3. P. 279-325.
17. Воеводин Вад.В., Овчинников С.Л., Романов С.Ю. Разработка высокоэффективных мас-
    штабируемых программ в задаче ультразвуковой томографии // Вычислительные методы и
    программирование: новые вычислительные технологии. 2012. Т. 13, № 1. С. 307–315.
18. Романов С.Ю. К вопросу об масштабируемости программы для обратной задачи волновой
    томографии // Вестник Нижегородского университета им. Н.И. Лобачевского. 2013. № 2-1.
    С. 160-167.
19. Goncharsky A.V., Romanov S.Y., Seryozhnikov S.Y. Inverse problems of 3D ultrasonic tomog-
    raphy with complete and incomplete range data // Wave Motion. 2014.Vol. 51, No. 3. P. 389–404.




                                                204
   Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org


20. Агаян Г.М., Романов С.Ю. Суперкомпьютерное моделирование в задаче ультразвуковой
    диагностики с применением аналитических подходов // Вестник Уфимского государствен-
    ного авиационного технического университета. 2013. Т. 17, № 5 (58). С. 260–269.
21. Агаян Г.М., Воеводин Вад.В., Романов С.Ю. О применимости послойных моделей в реше-
    нии трехмерных задач ультразвуковой томографии // Вычислительные методы и програм-
    мирование: новые вычислительные технологии. 2013. Т. 14, № 1. С. 533–542.
22. Романов С.Ю. Моделирование на суперкомпьютере в коэффициентных обратных задачах
    волновой томографии с учетом поглощения. // Вестник Уфимского государственного авиа-
    ционного технического университета. 2015. Т. 19, № 1(67). С. 224–231.




                                               205
   Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org



Research of tomographic schemes of low-frequency ultrasonic
diagnostics on supercomputers
Sergey Romanov
Keywords: Ultrasonics, 3D coefficient inverse problems, supercomputer, transmission and
reflection tomography
Various ultrasonic tomography schemes are compared. Inverse problems of ultrasonic
tomography are addressed as coefficient inverse problems for wave equation. The algorithms
for solving inverse problems of ultrasonic tomography are based on direct computation of the
gradient of the residual functional by solving the conjugate problem for the wave equation.
The potentialities of different reflection and transmission tomographic schemes are compared.
The results of the computation of model problems are analyzed in the 3D case and compared
to those of modeling of common ultrasonic diagnostic schemes, where the image is taken
from a single vantage point. Direct and inverse ultrasonic tomography problems are solved
using a GPU-based supercomputer.