Статья посвящена сравнению различных схем ультразвуковой (УЗ) томографии. Обратные задачи УЗ томографии рассматриваются как коэффициентные обратные задачи для волнового уравнения. Алгоритмы решения задачи УЗ томографии основаны на прямом вычислении градиента функционала невязки через решение сопряжённой задачи. Проведён сравнительный анализ возможностей различных томографических схем на прохождение, на отражение в трёхмерном случае. Полученные результаты сравниваются с результатами моделирования обычных УЗИ аппаратов, в которых изображение получается с одного ракурса. Для расчета прямых и обратных задач УЗ томографии использовался суперкомпьютер на GPU.
r R3 , Γ – граница области , nu |ΓT – производная вдоль нормали к поверхности Γ в области Γ 0,T , p(r,t) – некоторая известная функция. Скалярная волновая модель позволяет описать дифракцию, рефракцию, переотражение волн в среде.
Рассмотрим обратную по отношению к задаче (
0 Γ 2 1 T
u( ,t ) U ( ,t )2ddt .
Здесь u(r,t) – решение прямой задачи (
Обратные задачи томографии в волновой модели являются нелинейными, некорректно
поставленными задачами [12]. Минимизация функционала невязки в этом случае является
сложной проблемой [13]. Прорывные результаты в области решения задач волновой томографии
связаны с возможностью прямого вычисления градиента функционала uc по формуле (
T C ( u( c )) wt ( r ,t )ut ( r ,t )dt ,
0
где u(r,t) – решение задачи (
Предлагаемый в работе численный метод с явной разностной схемой относится к SIMD алгоритмам (Single instruction – multiple data) и имеет высокую степень параллелизма. Такие задачи могут быть эффективно реализованы на вычислительных системах с параллельной архитектурой. В работах [14, 17, 18] приведены результаты эффективного распараллеливания 2D задачи на кластерной архитектуре с использованием процессоров общего назначения, при использовании до десятков тысяч процессорных ядер. Графические процессоры (GPU) также позволяют получить очень высокую производительность на задачах такого типа. Вычислительная конфигурация с использованием GPU необходимая для решения задач УЗ томографии может быть реализована в рамках одной стойки, что делает возможным установку небольших, но высокопроизводительных суперкомпьютеров непосредственно в медицинских учреждениях и решение задач УЗ томографии на практике за приемлемое для медицинских исследований время. В настоящее время GPU успешно используются в задачах УЗ томографии [19].
За счёт наличия высокопроизводительной памяти современные GPU позволяют получить существенный выигрыш на data-intensive задачах, которые содержат относительно немного арифметических операций и много операций обращения в память. Именно к таким задачам относится явная разностная схема, используемая при решении обратной задачи УЗ томографии.
Несмотря на огромные размеры массивов, заданных на сетках размерностью ~n4 для 3D задачи, разработанный алгоритм решения таков, что хранить в памяти необходимо только данные для двух шагов по времени явной разностной схемы, размер которых ~ n3 (порядка 0.1% от числа вычисляемых значений). Это позволяет разместить их полностью в высокопроизводительной памяти 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 итераций программы Время, сек. (CPU) Размер сетки расчетов Время, сек. (GPU) Число процессов 16(=4×2×2) 144(=4×6×6) 400(=4×10×10) 784(=4×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. Модельные расчёты трёхмерных задач ультразвуковой томографии Вопросам моделирования в обратных задачах УЗ томографии посвящены публикации [2022]. Проведенные в настоящей работе модельные расчёты позволяют оценить возможности различных томографических схем в 3D варианте. Для сокращения времени расчётов размер исследуемого объекта был уменьшен до ≈60 мм. Расчётная область представляет собой куб размером H×H×H, где H=100 мм. При шаге разностной сетки 0,3 мм количество точек сетки составляет 320×320×320. Схема трёхмерного вычислительного эксперимента показана на Рис. 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%. Рис. 3. Сечения восстановленного 3D изображения по схеме на прохождение Как видно из сравнения Рис. 3 и 4, качество восстановленных изображений в схеме на прохождение хуже, чем в схеме с полным диапазоном углов. Мелкие детали изображения особенно для зашумленных данных уже не различаются.
Восстановленное по схеме на отражение изображение без шума в данных показано на Рис. 4 слева. Модельные расчёты показывают, что возможности томографических исследований 3D объектов в рамках схемы на отражение весьма ограничены.
Рис. 4 Сечения восстановленного 3D изображения по схеме на отражение Таким образом, как видно из проведённых модельных расчётов, томографические схемы с полным диапазоном углов обладают безусловными преимуществами по сравнению со схемами на прохождение или на отражение. В схеме с полным диапазоном данных хорошо восстанавливается не только форма неоднородных включений, но и пространственное распределение скорости звука. В схеме на отражение можно восстанавливать только границы неоднородностей, которые, как правило, сдвинуты вследствие неоднородной скорости. Восстановление абсолютного значения скорости невозможно. В задаче диагностики рака молочной железы это означает, что видимые на отражение неоднородности крайне трудно характеризовать, т.е. выяснить являются ли они злокачественной опухолью.
Проблема восстановления изображений на отражение является характерной не только для медицинских УЗ исследований. Аналогичная проблема возникает, например, при электромагнитном зондировании приповерхностных слоёв Земли, при поиске мин, в сейсморазведке, инженерной сейсмике. Во всех этих задачах используется диагностика на отражение.
Ситуация становится ещё хуже, если скорость звука внутри объекта варьируется в широких пределах. При увеличении диапазона изменения скорости звука в 2 раза в схеме на отражение получается изображение, показанное на Рис. 4 справа. Как видно, увеличение вариации скорости внутри объекта приводит к появлению артефактов, таких как удвоение границ и т.п. Литература 3. Буров В.А., Зотов Д.И., Румянцева О.Д. Восстановление пространственных распределений скорости звука и поглощения в фантомах мягких биотканей по экспериментальным данным ультразвукового томографирования // Акуст. журн. 2015. Т. 61, № 2. С. 254–273. 20. Агаян Г.М., Романов С.Ю. Суперкомпьютерное моделирование в задаче ультразвуковой диагностики с применением аналитических подходов // Вестник Уфимского государственного авиационного технического университета. 2013. Т. 17, № 5 (58). С. 260–269. 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.