CC BY
22
Программное обеспечение вычислительных, телекоммуникационных и управляющих систем
DOI: 10.18721/JCSTCS.11404
УДК 519.63-688; 534-8, 534.222
РЕШЕНИЕ ЗАДАЧИ ПОСТРОЕНИЯ УЛЬТРАЗВУКОВЫХ ИЗОБРАЖЕНИЙ
МЕТОДОМ ВОЛНОВОЙ ИНВЕРСИИ С ИСПОЛЬЗОВАНИЕМ АЛГОРИТМОВ РАСПАРАЛЛЕЛИВАНИЯ ВЫЧИСЛЕНИЙ
К.А. Верещагин, И.Н. Белых
Санкт-Петербургский политехнический университет Петра Великого,
Cанкт-Петербург, Российская Федерация
Активно развивающееся направление ультразвуковой томографии требует разработки эффективных подходов решения прямых и обратных задач для волнового уравнения как на теоретическом, так и на прикладном уровне. В статье изучены методы, применяемые для решения задач реконструкции изображений из волнового поля, и предложен оригинальный алгоритм построения двухмерных срезов ультразвуковых изображений. Для акустически неоднородной модели упругой среды проведено решение прямой задачи для волнового уравнения в двухмерной постановке с учетом поглощения с целью получения полного поля отраженных волн. Представлен оптимизационный алгоритм 2D-волновой инверсии ультразвукового поля для решения обратной задачи во временной области. Приведены основные архитектурные элементы программного решения, рассмотрены способы распараллеливания вычислений для повышения производительности процесса моделирования. Обсуждены результаты тестирования и проверки на корректность различных этапов алгоритма. Достоинствами предложенного подхода являются точность и устойчивость решения при повышенной вычислительной производительности.
Ключевые слова: волновое уравнение, отраженные волны, волновая инверсия, оптимизация, параллельные вычисления, ультразвуковая томография.
Ссылка при цитировании: Верещагин К.А., Белых И.Н. Решение задачи построения ультразвуковых изображений методом волновой инверсии с использованием алгоритмов распараллеливания вычислений // Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление. 2018. Т. 11. № 4. С. 49-62. DOI: 10.18721/JCSTCS.11404
ULTRASOUND IMAGE RECONSTRUCTION BASED ON WAVEFORM INVERSION METHOD USING PARALLEL CALCULATION ALGORTIHMS
K.A. Vereshchagin, I.N. Belykh
Peter the Great St. Petersburg Polytechnic University, St. Petersburg, Russian Federation
Ultrasound tomography is actively developing area that requires effective research and development efforts in forward and inverse problems solution for wave equation. In this work different methods of image reconstruction from elastic wave field are
investigated and original algorithm for two-dimensional ultrasound imaging is suggested. Full reflection wave field including absorption was modeled by means of forward problem for 2-D acoustic wave equation solution for inhomogeneous elastic media. Invers problem is solved by optimization algorithm for two-dimentional full waveform inversion of reflection wave field. The main features of the software application design are described, using parallel calculation techniques for the modeling process efficiency improvement. The testing results are discussed and algorithm steps are validated for correctness. The advantages of the proposed approach are accuracy and stability of the solution with improved computational efficiency.
Keywords: wave equation, reflection waves, waveform inversion, optimization, parallel calculations, ultrasound tomography.
Citation: Vereshchagin K.A., Belykh I.N. Ultrasound image reconstruction based on waveform inversion method using parallel calculation algortihms. St. Petersburg State Polytechnical University Journal. Computer Science. Telecommunications and Control Systems, 2018, Vol. 11, No. 4, Pp. 49-62. DOI: 10.18721/JCSTCS.11404
Введение
Ультразвуковые исследования (УЗИ) широко используются в медицине и неразру-шающем контроле для дистанционного определения внутренней структуры изучаемых объектов. В настоящее время активно развиваются различные подходы к ультразвуковой томографии (УЗТ), являющейся, в отличие от компьютерной (КТ) и магниторе-зонансной (МРТ) томографии, не инвазив-ным методом, позволяющим обследовать целевые объекты в трехмерном пространстве. Подходы УЗТ различаются по типам используемых волн — отраженных или проходящих, что достигается различной геометрией расположения источников и приемников ультразвуковых (УЗ) колебаний, а также методами реконструкции пространственных изображений из зарегистрированного волнового поля. В одном из развитых подходов трансмиссионной УЗТ на проходящих волнах [4] исследуемый объект располагается между источниками и приемниками, лежащими в одной плоскости, которая пошагово перемещается в перпендикулярном направлении. Такая геометрия расстановки датчиков накладывает ограничения на размер объекта.
В настоящей статье рассматривается алгоритм получения изображения трехмерной модельной среды, содержащей несколько объектов с различными акустическими свойствами, путем регистрации на ее поверхности отраженного волнового поля линейкой трансдьюсеров и получения
наборов двухмерных срезов при ее перемещении в направлении перпендикулярном линии датчиков.
При распространении волн в гетерогенных упругих средах происходят различные эффекты, искажающие волновое поле. Это многократные отражения волн от границ раздела сред, рассеивание и затухание энергии в различных материалах, а также эффект дифракции волн на неоднородно-стях. В результате решения обратной задачи УЗТ эти эффекты минимизируются с целью получения итогового пространственного изображения внутренней структуры исследуемого объекта.
Большой практический интерес представляет эффективное решение обратной задачи в двухмерной постановке для отраженных волн при расположении линейки трансдьюсеров на поверхности исследуемой среды с возможностью ее перемещения. Такая схема наблюдения не накладывает ограничений на размер объекта, являясь промежуточным этапом УЗТ, путем построения высокоразрешающих объемных УЗ-изображений набором двухмерных срезов.
В теоретических работах по сейсморазведке развито множество методов решения обратной задачи, идея большинства из которых заключалась в упрощении физической картины рассматриваемого явления распространения волн и использования [8] алгоритмов корректировки возникающих эффектов реверберации, дифракции и поглощения.
Устранение артефактов дифракции позволяет более точно определять форму объектов и идентифицировать мелкие объекты на ультразвуковых изображениях. Долгое время применялись, в основном, аппаратные методы устранения дифракционных эффектов. Теоретические подходы для разработки программных методов, основанные на расчете дифракционных искажений и их подавлении, исследованы в [14]. В основу этих методов легли идеи, использовавшиеся в задачах обработки сейсмических данных. В этой области наибольшее распространение получили различные методы сейсмической миграции, такие как миграция Кирхгофа, конечно-разностная миграция и миграция в Фурье-области [15]. Эти методы были популярны и активно развивались в то время, когда мощностей вычислительных систем не хватало для полного численного решения волнового уравнения.
Наиболее перспективный метод решения обратной задачи — метод волновой инверсии, основанный на решении прямой и обратной задач построения УЗ-изображений методом численного решения волнового уравнения, наиболее точно описывающего физическую картину явления. Эти методы ранее нашли применение в области инверсии сейсмических данных при решении задач сейсморазведки, являясь, однако, достаточно ресурсоемкими вследствие большого объема данных. Актуально исследование применимости этих методов в области построения УЗИ, где в силу повышенных требований к разрешающей способности объем данных соизмерим или может даже превосходить объемы трехмерных сейсмических данных.
В последние годы, в связи с быстрым развитием вычислительной техники, в частности, массивно-параллельных систем, таких как ОРи-кластеры и суперкомпьютеры, все больше внимания стало уделяться возможностям применения сложных волновых методов как в сейсмике, так и в области ультразвуковых исследований.
Из аналитических трудов по сейсморазведке можно отметить [1, 9, 13], где прора-
ботана базовая теория, которая в дальнейшем легла в основу многочисленных исследований.
За последнее десятилетие существенного теоретического и практического продвижения в решении обратных задач с применением для УЗ медицинских приложений удалось достичь лаборатории Los Alamos в содружестве с коллективами из нескольких университетов США [7] и группе авторов из Московского государственного университета имени М.В. Ломоносова [16]. Большинство этих исследований посвящено тороидальной системе трансдьюсеров, позволяющей регистрировать полное волновое поле как проникающих, так и отраженных волн, но накладывающей ограничения на размер исследуемого объекта. В этих работах рассматривались варианты решения волнового уравнения разностными схемами во временной области [16] либо волновой инверсии в частотной области [7].
Из разнообразия рассматриваемых систем наблюдения можно отметить [4], где применялась модель двух параллельных плоскостей, реализованная пошаговым сканированием объекта двумя линейками трансдьюсеров, расположенных по разные стороны от него. Эта система принимает как проникающие, так и отраженные волны, но с потерями части волновой информации. В работе подробно рассмотрено влияние разных видов регуляризации на качество результирующих изображений.
В настоящей статье приводятся решения прямой и обратной задач построения ультразвуковых изображений методом инверсии поля отраженных волн во временной области для двухмерного случая.
Решение прямой задачи
Определим формально прямую задачу построения ультразвуковых изображений. Входные данные задачи:
акустические параметры гетерогенной среды известны;
линейка из N трансдьюсеров расположена на поверхности исследуемой области; форма исходного импульса задана.
Для каждого трансдьюсера я моделирование волнового процесса происходит следующим образом: трансдьюсер я, работая в режиме источника сигнала, излучает начальный импульс, после чего все N тран-сдьюсеров, включая я, работают в режиме приемника в течение некоторого времени Т, принимая отраженные от среды волны. В результате каждого излучения получается ансамбль из N реализаций волнового процесса, в котором каждая /-я реализация представляет собой сигнал длительностью Т, зарегистрированный одним приемником. Решением прямой задачи является набор из N таких ансамблей, рассчитанных для каждого из N трансдьюсеров.
В методах волновой инверсии решение прямой задачи осуществляется при помощи численного решения волнового дифференциального уравнения, имеющего вид:
1 д2 р( х, ) 2
c 2( х)
ды2
-V2p(х, t; х ) = s(х, t; х ), (1)
где хя — положение источника импульса; р (х, t; х5) — амплитуда волны в точке х в момент времени с(х) — скорость звука в точке х; s(х, t;х5) — первоначальный сигнал.
Существуют различные подходы к численному решению такого уравнения, рассмотрим самые распространенные из них.
Метод конечных разностей. Это дифференциальное уравнение относится к гиперболическому типу, поэтому для его численного решения на некоторой сетке удобно применять метод конечно-разностной аппроксимации с крестообразным или Т-образным шаблоном для аппроксимации производных второго порядка, используя известные соотношения для производных:
д2 р _ рТ^+рТ
dt2
- + 0(т2),
д2Р = PU - 2Рk + PU
+ 0(h2),
(2)
(3)
дх2 к2
где (/, ]) — узел сетки; т и к — шаги сетки по оси времени и по пространственной оси.
Эта схема является явной, т. к. значения в узлах шаблона на нижних временных слоях должны быть известны, поэтому она более простая с вычислительной точки зрения.
Однако на практике, при реализации метода на тестовой 2Б-модели были выявлены определенные недостатки: условия на устойчивость и сходимость метода требуют выбора очень маленького шага дискретизации по пространственной оси, около 10 точек на длину волны сигнала [11]. В результате задача моделирования распространения высокочастотных волн, таких как ультразвук, становится очень затратной с вычислительной точки зрения, что критично для реализации следующего важного этапа алгоритма обращения — решения обратной задачи, при котором моделирование прямой задачи происходит многократно и итеративно.
Псевдоспектральный метод. Рассмотрим альтернативные методы решения волнового уравнения во временной области, обладающие меньшей вычислительной сложностью, чем метод конечных разностей. Одним из самых популярных методов, используемых для решения подобной задачи на практике, является псевдоспектральный метод, активно применявшийся для моделирования акустических волн в задачах сейсморазведки [10].
Основная идея метода состоит в том, что вместо аппроксимации пространственных производных конечными разностями применяется метод Фурье, основанный на численном нахождении производных в спектральной области и использовании обратного преобразования Фурье. При этом производные по времени аппроксимируются конечной разностью. Для уменьшения количества точек по времени применяется метод ^-пространств, поглощающие граничные условия реализуются методом PML (Perfect Matched Layer) [11]. Данный метод обладает заметным преимуществом по скорости работы.
Недостаток заключается в том, что метод Фурье, несмотря на свою скорость, допус-
кает заметные погрешности, выражающиеся в виде появления мелких осцилляций в пространстве (эффект Гиббса) в том случае, если в спектре исходного сигнала есть резкие перепады частот, превышающие частоту, поддерживаемую дискретной временной шкалой по теореме Котельникова— Найквиста [11]. Это накладывает определенные ограничения на среду моделирования и требует повышенного внимания к спектральной картине, а в некоторых ситуациях может быть необходимо применение оконных сглаживающих функций. В данной работе рассматриваемые сигналы почти полностью лежат в ограниченной полосе спектра, что избавляет от этих проблем и позволяет применять метод на практике.
Для использования метода в прикладных задачах в разных языках программирования существуют специальные библиотеки численного моделирования, например, библиотека k-Wave Toolbox для языка MATLAB.
Решение обратной задачи
Исходные данные обратной задачи построения ультразвуковых изображений — данные о волновом поле в результате решения прямой задачи: N ансамблей, каждый из которых состоит из N реализаций длительностью T каждая. Решением обратной задачи является восстановленное по исходному волновому полю пространственное распределение параметров неоднородной среды, через которую прошел излученный сигнал, в данном случае — значения скорости звука в каждой точке пространства.
Метод волновой инверсии во временной области. Основная идея метода волновой инверсии при решении обратной задачи состоит в использовании знаний о математической модели распространения волнового поля для численного моделирования физического процесса и нахождения тех параметров уравнения, которые позволяют наиболее точно смоделировать данные, полученные в результате физического эксперимента.
Задача волновой инверсии обычно формулируется в виде задачи минимизации некоторой целевой функции, характеризующей отклонение смоделированных данных от экспериментальных. В качестве целевой функции часто берут расстояние между данными, порожденное Х2-нормой [4]:
E(c) = ^ZZit5^' t; xrfdt, (4)
2 s r 0
Sp(Xr, t; Xs) = p(xr, t; Xs) - pobs(X-, t; Xs), (5)
где s — номер источника сигнала; r — номер приемника сигнала; c — текущая модель скоростей; р(х„ t; xs) и pobs(xr, t; xs) — модельные и реальные данные, полученные приемником r в момент времени t в результате эксперимента с источником s.
Задача для волнового уравнения решается численно. Для определения направления минимизации на каждой итерации вычисляется градиент функции E в аппроксимации ¿к) « VE(c(k)), после чего происходит сдвиг текущей модели скоростей по следующей формуле [3]:
c( k+1) = c(k) +а( k) p( k), (6)
где c(k) — модель скоростей, построенная на к-й итерации; а(к) — шаг минимизации; p(k) — направление минимизации, зависящее от оптимизационной стратегии,
например, в методе градиентного спуска
р(к) = —gg, (7)
а в методе сопряженных градиентов
p(k) =-g (k) +ß(k) p(k-1), (8)
где ß(k) — параметр, который может быть выбран по-разному, в зависимости от конкретизации метода.
Оба этих метода относятся к методам локальной оптимизации, т. е. в зависимости от выбора начального приближения модели скоростей с(0) алгоритм может сойтись к экстремуму, не являющемуся глобальным минимумом целевой функции, что особенно проявляется при решении нелинейных и невыпуклых задач, к кото-
рым относится и задача инверсии данных при помощи волнового уравнения.
Выполнение каждой итерации алгоритмов подразумевает вычисление двух основных параметров: вектора градиента и скалярного значения шага оптимизации.
Вычисление градиента целевой функции. Обычно для вычисления вектора частных производных используют конечно-разностные соотношения вида:
д (Xl,..., xn)
ду
/ (XI,..., х{ + к,..., Хп)-/(XI,..., Хп )
к '
(9)
где к/ — шаг дискретизации для переменной х/.
Однако для данной задачи такая схема является крайне неэффективной, т. к. для вычисления всего вектора градиента потребуется N полных решений прямой задачи, где N — количество переменных в модели скоростей.
Вместо этого в задачах такого масштаба для вычисления приближенного значения градиента может использоваться метод сопряженного состояния, в котором значение частной производной вычисляется по следующей формуле [1]:
дЕ
дс( хк) с( хк )3
хУр2 р( хк, t; х)
(10)
s 0
дt2
р (хк, Т -1; х8 )Л,
1 д2 р (хк ,Т - г; х8)
с2( х)
ы2
-V2 р (хк, Т -1; х8) _
_^8р(х^, V; х8)8(х - х^),
(11)
новной прямой задачи, вычисляется разность 8р(хр ?; хя), после чего происходит ее обращение во временной оси. Полученные данные трактуются как синхронный импульс для каждого трансдьюсера, и для решения сопряженного волнового уравнения моделирование процесса производится с использованием уже готовой реализации.
Данный метод обладает заметным преимуществом по вычислительной производительности, т. к. требует только одного дополнительного решения прямой задачи для вычисления градиента. Недостатком, как видно из формулы, является необходимость моделирования и запоминания не только сенсорных данных, но и полного волнового поля решений основного и сопряженного уравнения, что существенно увеличивает объем сохраняемых в памяти данных.
Для упрощения дальнейшего вычисления шага минимизации можно считать градиент нормированным; направление от этого не меняется, а алгоритм вычисления шага становится более универсальным [2]. Формулы для вычисления направления на к-й итерации принимают следующий вид:
р(к) _-
у(к)
(к)
- + Р(к) р(к-1),
(13)
где ^ — градиент, В(к) — коэффициент, который может принимать разные значения, например [3]
где р"(хк, ?; хя) — решение прямой задачи для сопряженного уравнения:
В(к) _. PFR -
(к)
(к-1)
(14)
Г ,,
в методе Флетчера-Ривза, или
(8(к))Т(8(к) -8(к-1)
В®
)
(к-1)
(15)
8р(х8, V; х8) _ р(х8, ?; Xs) - роЬ8(х8, V; х8). (12)
Моделирование решения сопряженного уравнения происходит следующим образом: для каждого источника я на основе уже имеющихся данных от каждого приемника g, посчитанных при решении ос-
в методе Полака-Рибьера.
Вычисление шага минимизации. После вычисления направления минимизации необходимо оценить размеры шага а, который нужно сделать в этом направлении для получения новой модели скоростей. Оценка вычисляемого шага обычно произ-
2
8
водится при помощи проверки выполнения определенных условий, гарантирующих, что шаг не будет слишком большим или маленьким. На практике часто используют условия Агтцо-Оо№1ет [2]:
E (c + ap) < E (c) + Pap,
(16)
E(c + ap) <E(c) + c1aVE p,
|VE (c + ap)T
< c
VETp,
(17)
(18)
Таким образом, регуляризованная функция ошибки имеет следующий общий вид:
Er (c) = E (c) + XR(c),
(19)
где р е (0,1), обычно берут значения порядка 10—4—10—3.
Алгоритм вычисления шага выглядит следующим образом: вначале берется некоторый шаг атах, после чего происходит проверка условий, в случае их невыполнения берется новый шаг та, где т — коэффициент от 0 до 1, после чего алгоритм повторяется до тех пор, пока условия не будут выполнены.
Выполнение этого шага обычно требует дополнительного решения 1-2 прямых задач. Это увеличивает вычислительную сложность каждой итерации, но улучшает их точность и уменьшает их количество. Соблюдение этого баланса в целом — главная проблема при решении оптимизационных задач.
Данный алгоритм гарантирует только то, что шаг не будет слишком большим. Для более точных оценок при оптимизации часто используют более жесткие условия Вольфе [2]:
где X — коэффициент регуляризации, определяющий степень влияния регуляризиру-ющей функции на результат. Значение коэффициента определяется эмпирически для конкретной задачи. В качестве регуляризи-рующей функции часто берут, например, L1 или L2 норму (регуляризация Тихонова). Для обращения волновых данных важно, чтобы регуляризирующая функция, помимо прочего, сохраняла края объектов на изображении для большей разрешающей способности. Поэтому при построении томографических изображений часто используют TV-регуляризацию (total variation), функция которой имеет следующий вид [5]:
Лту(с(х)) = ^/е2 + ||Ус( х )||2, (20)
х
где х — пространственная координата; с(х) — распределение скорости; е — малое слагаемое, обеспечивающее дифференци-руемость при значении нормы, равной нулю во всей области.
Выражение для вычисления градиента этой функции выглядит следующим образом [6]:
VRTV(c( х)) = V-
Vc
Vе2 + II Vc( х )||2
(21)
где с е (0,1), С2 е (сь 1).
Однако выполнение этих условий требует вычисления градиента на каждом этапе, что очень сильно усложняет процесс вычислений, поэтому в данной работе было решено использовать условия первого типа.
ТУ-регуляризация. Одним из способов преодоления неустойчивости оптимизационных задач является регуляризация — добавление некоторого дополнительного слагаемого к условию функции ошибки. Чаще всего регуляризация представляет собой некоторую функцию, зависящую от каких-то априорных данных о модели, в данном случае — изначальном распределении скоростей.
Метод упорядоченных подмножеств. Как
упоминалось ранее, вычисление градиента методом сопряженного состояния (adjoint-state method) осложнено необходимостью считывать и обрабатывать на каждой итерации огромные объемы данных, что при ограничениях на систему ввода-вывода существенно увеличивает время на т. н. «накладные расходы». Для ускорения процесса обращения данных в задачах томографии можно использовать метод упорядоченных подмножеств (ordered-subsets method) [5]. Метод состоит в том, что данные сенсоров делятся на несколько наборов, каждый из которых содержит опреде-
ленное число срезов (под срезом здесь подразумевается набор данных, полученных от одного источника). Для каждого такого набора производится шаг итерации, минимизирующий некоторую локальную функцию ошибки, определенную отдельно для каждого набора к. Для того чтобы минимизирующий алгоритм работал корректно для небольших подмножеств, функцию ошибки было решено нормализовать относительно реальных данных срезов данного подмножества. Таким образом, из формулы (4) для определенного далее нормализующего коэффициента получаем формулу:
Е{с)=2Кт\ ЕЕЯ^,';х*)]2(22)
2К (с) *еИк г 0
8р( хг, г; х*) = р( хг, г; х*) - роЫ (хг, г; х*), (23) К (с) = Е Е (Роь* (хг, г; х5 ))2 йг, (24)
г о
где индексы суммирования я и г обозначают, соответственно, номера источников
и приемников сигнала; Мк — индексное множество всех источников, принадлежащих к-му набору; с — текущая модель скоростей; р(Хг, Г; х) и роЬАх„ Г; х,) — модельные и реальные данные, полученные приемником г в момент времени / в результате эксперимента с источником К(с) — функция нормализации ошибки.
Таким образом, каждая итерация требует работы со значительно меньшим объемом данных, это увеличивает количество шагов, сделанных при анализе полного набора, что существенно ускоряет процесс без особого ухудшения качества итогового изображения [5].
Программная реализация
Общая архитектура программы. Для программной реализации рассмотренных ранее методов был выбран язык МЛТЬЛБ. На рис. 1 представлена иМЬ-диаграмма основных классов, реализующих необходимые структуры данных и алгоритмы решения прямой и обратной задач.
Рис. 1. Диаграмма взаимодействия основных классов проекта Fig. 1. Diagram of project main classes interaction
Представленные на рисунке основные классы взаимодействуют следующим образом:
• KWaveFPSolver — реализует алгоритм решения прямой задачи для псевдоспектрального метода. Класс реализует интерфейс ForwardProblemSolver, содержащий метод solve, принимающий объект класса Model и возвращающий решение прямой задачи;
• InverseProblemCUDASolver — класс, реализующий весь необходимый функционал для решения обратной задачи: вычисление градиента, регуляризация, вычисление шага минимизации и т. д. При создании объекта в конструктор передается объект класса KWaveFPSolver, используемый далее в методах для решения прямой задачи, а также параметры регуляризации Lambda и Epsilon;
• Model — класс, агрегирующий в себе объекты, описывающие все основные параметры моделируемого эксперимента. Объекты этого класса передаются в качестве параметра всем методам, работающим с экспериментальной моделью. Полями объекта класса Model являются объекты классов Medium, Sensor, Sources, Grid. Это позволяет легко модифицировать модель при реализации разных алгоритмов;
• Medium — класс параметров среды, содержит поля значений скорости звука и плотности в каждой точке;
• Sensor — класс, задающий расположение приемников сигнала относительно исследуемой области;
• Sources — класс параметров источника сигнала, задает расположение источников и форму начального сигнала для каждого источника;
• Grid — класс параметров дискретизации модели, задает шаг дискретизации и количество отсчетов по осям X, Y и T.
Для моделирования волнового процесса при решении прямой задачи использовалась библиотека k-Wave Toolbox, неоднократно применявшаяся в научных исследованиях [10].
Распараллеливание вычислений на CPU и GPU. Последовательные реализации алгоритмов решения прямой и обратной задачи являются вычислительно затратными: этапы нахождения направления и шага оптимизации суммарно требуют дополнительного решения как минимум 5-6 прямых задач. Для моделируемой области средних размеров (около 200x300 дискретных значений) с более чем 100 трансдьюсерами, решение одной прямой задачи на типичном персональном компьютере (4 ядра, 2.6 ГГц, 8 Гб RAM, 500 Гб HDD), использованном для тестирования, занимает более получаса.
Кроме того, на этапе вычисления градиента методом сопряженного состояния требуется сохранять значения амплитуд не только в местах расположения трансдьюсе-ров, но и во всей остальной области в каждый момент времени. Объем обрабатываемых на каждой итерации данных для упомянутой модели составляет примерно 120 Гб. Данные такого объема невозможно хранить в оперативной памяти, поэтому при выполнении необходимых вычислений происходит нагруженный дисковый обмен. Суммарное время операций чтения, обработки и записи для одной итерации для HDD дисков составляет примерно 1,5 ч.
Вследствие высокой ресурсоемкости процесса моделирования становится актуальной задача оптимизации процесса решения прямой задачи с использованием технологий распараллеливания вычислений на CPU и GPU. Библиотека k-Wave, используемая для численного решения волнового уравнения при моделировании прямой задачи, содержит отдельные модули, предназначенные для более эффективного моделирования волнового процесса в 3D-областях большого масштаба. Эти модули написаны на языке программирования C++ с применением технологий OpenMP и CUDA и позволяют эффективно использовать вычислительные ресурсы компьютера.
Для возможности использования этих модулей в решении прямой задачи построения ультразвуковых 2D-изображений, двухмерная задача была представлена в ви-
де трехмерной с пренебрежительно малым шагом в третьем измерении. Результат моделирования для такого представления оказался абсолютно соответствующим обычной 2Б-реализации, что подтверждает корректность такого использования модулей. Таким образом, в работе заложена возможность расширения предлагаемого двухмерного подхода до трех пространственных переменных в продолжение работы [8] для решения задачи получения томографических ультразвуковых изображений. Также для распараллеливания процесса решения уравнений для разных источников и чтения многочисленных файлов сенсорных данных использовался MATLAB Parallel Toolbox, позволяющий распределить нагрузку по вычислению данных каждого источника между несколькими worker-процессами при помощи операций выделения пула и распараллеливания циклов for.
В результате такой оптимизации время последовательного вычисления одной полной итерации для рассмотренной далее тестовой модели уменьшилось с 3,5 до 1,5 ч при использовании пула worker-процессов и до 27 мин при подключении библиотеки, использующей GPU. Распараллеливание проводилось с применением CPU Intel Core i5-6200U (4 ядра, 2.6 ГГц) и GPU NVIDIA GeForce GTX 950M (640 CUDA-ядер, 914 МГц).
Тестирование программы и анализ результатов
Тестирование минимизации функционала.
Проверим, что оптимизирующий алгоритм работает корректно, т. е. значение функции ошибки убывает с каждой итерацией.
Многократное осуществление итераций занимает очень много времени для крупных моделей, поэтому тестирование проводилось на небольшой, но структурно представительной модели, приведенной на рис. 2. Количество дискретных значений в пространстве N = 256, N = 150; шаг дискретизации в пространстве (х = ¿у = 1 • 10—4 м; дискретизация по времени N = 1160 с, Ж = 3 • 10—8 с; скорость звука в среде составляет 1480 м/с, в центр помещен объект плотностью 1500 м/с. В объект добавлены вкрапления, скорость звука внутри которых равна, соответственно, 1570 и 1640 м/с. Трансдьюсеры расположены вдоль оси X с шагом равным четырем.
Возьмем в качестве начальной модели однородную среду плотностью 1480 м/с и запустим 45 итераций оптимизирующего алгоритма.
График зависимости функции ошибки от числа итераций изображен на рис. 3. По графику видно, что зависимость функции ошибки от количества итераций убывающая, что подтверждает корректность оптимизирующего алгоритма.
Рис. 2. Структура модели среды в виде пространственного распределения скорости звука в м/с, отображаемой в полутоновой шкале
Fig. 2. Model structure as the sound speed dimensional distribution represented in grayscale
ylQ-7
Рис. 3. Зависимость значения функции ошибки от числа итераций
Fig. 3. Error function dependency on number of iterations
Тестирование ТУ-регуляризации. Для
проверки корректности ТУ-регуляризации и реализации ее градиента был запущен отдельно градиентный спуск для функции ЯТУ. Значение функции может вести себя по-разному в зависимости от параметра е, поэтому было проведено два эксперимента с разными значениями е.
Из рисунка 4 видно, что при очень малых значениях параметра е градиент становится слишком чувствительным к величине шага, и процесс минимизации происходит некорректно. При достаточно большом е алгоритм работает корректно.
а)
а)
10 20 30 40 50
итерации
Рис. 4. Зависимость значения функции регуляризации от числа итераций для малого (а) и большого (б) параметра £
Fig. 4. Regularization function dependency on number of iterations for small (a) and large (б) value of £
Рис. 5. Восстановленное пространственное распределение скоростей звука для исходной модели (а) и со сдвинутыми вправо по оси Х внутренними объектами (б)
Fig. 5. Restored sound speed spatial distribution for initial model (а) and with internal objects shifted to the right along the Х axis (б)
3,2
3
2.8
2.6
2.4
2.2
а)
б)
Рис. 6. Восстановленное пространственное распределение скоростей звука для исходной модели с поглощением (а) и со сдвинутыми вправо по оси Х внутренними объектами (б)
Fig. 6. Restored sound speed spatial distribution for initial model with acoustic absorption (a) and with internal objects shifted to the right along the Х axis (б)
Тестирование разрешающей способности алгоритма. Разрешающая способность алгоритма в первую очередь зависит от выбираемого им направления оптимизации.
Некоторые детали могут быть видны после первой итерации, а какие-то так и не проявляются. Для изучения этого вопроса рассмотрим несколько небольших представительных моделей и построим изображения восстановленного пространственного распределения скоростей звука после выполнения нескольких итераций.
Для начала возьмем модель, рассмотренную ранее. Из рисунка 5 видно, что четко выделяются передняя и задняя границы объекта среды и мелких вкраплений, при этом плохо прорисованы боковые участки объекта.
Это может быть связано с тем, что волны, отражающиеся от боковых субвертикальных участков границ, слабо улавливаются сенсорами, т. к. объект расположен строго в центре. Рассмотрим случай объекта, сдвинутого относительно линии сенсоров, переместив его на 60 единиц вправо по оси X.
На рис. 5 б видно, что небольшой сдвиг никак не повлиял на степень детализации боковых участков, из чего можно сделать вывод, что длина линии сенсоров должна значительно превосходить размеры предмета для его качественного отображения.
Теперь изучим влияние параметра поглощения среды на качество изображения
объекта. Будем считать, что поглощение происходит по экспоненциальному закону, зависящему от частоты сигнала [10].
На рис. 6 видно, что восстановленное изображение получается более размытым, менее четко прорисованы задние участки границ, т. к. часть волновой информации поглощается. Кроме того, сильнее проявляются поперечные артефакты в виде полос, особенно это видно по изображению со сдвинутым объектом. Такие артефакты могут быть подавлены дальнейшей обработкой полученных изображений и последующим применением алгоритмов сегментации.
Заключение
Предложен оригинальный алгоритм решения прямой и обратной задач построения плоских ультразвуковых изображений методом волновой инверсии отраженного поля во временной области. Алгоритм обладает рядом преимуществ по сравнению с другими алгоритмами, упомянутыми выше. Во-первых, предложенный алгоритм вычислительно более устойчив при наличии шумов, чем алгоритмы, использующие конечно-разностные схемы. Во-вторых, алгоритм работает только с отраженной частью волнового поля, регистрируемой на поверхности среды, что снимает ограничения на размеры изучаемых объектов. Это преимущество позволяет расширить область применения предлагаемого метода не
только для медицинских приложений, но и для задач неразрушающего контроля. Объемное изображение может быть представлено набором двухмерных срезов. Недостатком подхода на отраженных волнах может быть отсутствие отражений при больших углах падения волн, что налагает определенные требования к размерам и геометрии системы наблюдения. Разработанный алгоритм программно реализован с распараллеливанием вычислений, что примерно на порядок повысило производительность процесса моделирования.
Тестирование алгоритма проведено на упрощенных, но акустически представительных моделях. Упрощение пространственной модели среды, состоящей из не-
скольких объектов относительно простой формы, принципиально не влияет на решение задачи и было обусловлено только вычислительными мощностями тестового компьютера.
Дальнейший интерес представляет решение задачи в трехмерной постановке для получения полномасштабных высокоразрешающих томографических ультразвуковых изображений сложно устроенных модельных сред. Прямые и обратные задачи в такой постановке необходимо алгоритмизировать и решать на высокопроизводительных параллельных системах, таких как суперкомпьютер или ОРи-кластер, вследствие существенного роста объема данных и вычислений.
СПИСОК ЛИТЕРАТУРЫ
1. Boonyasiriwat, P. Valasek, P. Routh, W., et al.
An efficient multiscale method for time-domain waveform tomography // Geophysics. 2009. Pp. WCC59-WCC68.
2. Nocedal J., Wright S. Numerical optimization. Springer, 2000. Pp. 30-63.
3. Liu Y., Teng J., Xu T., Badal J., Liu Q., Zhou B. Effects of conjugate gradient methods and step-length formulas on the multiscale full waveform inversion in time domain: numerical experiments // Pure and Applied Geophysics. 2017. Pp. 1983-2006.
4. Lin Y., Huang L., Zhang Z. Ultrasound waveform tomography with the total-variation regulariza-tion for detection of small breast tumors // SPIE Medical Imaging. 2012. Pp. 832002-1-832002-9.
5. Hudson H.M., Larkin R.S. Accelerated image reconstruction using ordered subsets of projection data // IEEE Trans. Med. Imaging. 1994. Pp. 100-108.
6. Peyre G. The numerical tours of signal processing - advanced computational signal and image processing // IEEE Computing in Science and Engineering. 2011. Pp. 94-97.
7. Sandhu G.Y., Li C., Roy O., Schmidt S., Duric N. Frequency domain ultrasound waveform tomography: breast imaging using a ring transducer // Physics in Medicine and Biology. 2015. Vol. 60. No. 14. Pp. 5381-5398.
8. Belykh I., Markov D. Volumetric ultrasound imaging: modeling of waveform inversion // Computer Science Research Notes 2703. 2017. Pp. 77-82.
9. Virieux J., Operto S. An overview of full waveform inversion in exploration geophysics // Geophysics. 2009. Pp. WCC1—WCC26.
10. Treeby B.E., Jaros J., Rendell A.P., Cox B.T. Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method // The Journal of the Acoustical Society of America. 2012. Pp. 4324—4336.
11. Langtangen P.H., Linge S. Finite difference computing with PDEs // Texts in Computational Science and Engineering. 2017. Pp. 93—205.
12. Li C., Duric N., Huang L. Clinical breast imaging using sound-speed reconstructions of ultrasound tomography data // Medical Imaging Int. Society for Optics and Photonics. 2008. Pp. 692009-1—692009-9.
13. Pratt R.G. Seismic waveform inversion in the frequency domain. Part 1: Theory and verification in a physical scale model // Geophysics. 1999. Vol. 64. Pp. 888—901.
14. Huang L., Williamson M. Detection of breast microcalcifications with super-resolution ultrasound imaging: A clinical study // SPIE Medical Imaging. 2013. Pp. 867510-1—867510-10.
15. Li J. Kirchhoff common offset migration velocity analysis via differential semblance // Master's Thesis, Rice University. 2007. Pp. 1—61.
16. Гончарский А.В., Романов С.Ю., Серёж-ников С.Ю. Обратные задачи послойной ультразвуковой томографии с данными на цилиндрической поверхности // Вычислительные методы и программирование. 2017. Т. 18. С. 267—276.
Статья поступила в редакцию 09.07.2018.
REFERENCES
1. Boonyasiriwat, P. Valasek, P. Routh, W., et al.
An efficient multiscale method for timedomain waveform tomography. Geophysics, 2009, Pp. WCC59-WCC68.
2. Nocedal J., Wright S. Numerical optimization. Springer, 2000, Pp. 30-63.
3. Liu Y., Teng J., Xu T., Badal J., Liu Q., Zhou B. Effects of conjugate gradient methods and step-length formulas on the multiscale full waveform inversion in time domain: numerical experiments. Pure and Applied Geophysics, 2017, Pp. 1983—2006.
4. Lin Y., Huang L., Zhang Z. Ultrasound waveform tomography with the total-variation regulariza-tion for detection of small breast tumors. SPIE Medical Imaging, 2012, Pp. 832002-1—832002-9.
5. Hudson H.M., Larkin R.S. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imaging, 1994, Pp. 100—108.
6. Peyre G. The numerical tours of signal processing — advanced computational signal and image processing. IEEE Computing in Science and Engineering, 2011, Pp. 94—97.
7. Sandhu G.Y., Li C., Roy O., Schmidt S., Duric N. Frequency domain ultrasound waveform tomography: breast imaging using a ring transducer. Physics in Medicine and Biology, 2015, Vol. 60, No. 14, Pp. 5381—5398.
8. Belykh I., Markov D. Volumetric ultrasound imaging: modeling of waveform inversion. Computer Science Research Notes 2703, 2017, Pp. 77—82.
9. Virieux J., Operto S. An overview of full waveform inversion in exploration geophysics. Geophysics, 2009, Pp. WCC1—WCC26.
Received 09.07.2018.
10. Treeby B.E., Jaros J., Rendell A.P., Cox B.T.
Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method. The Journal of the Acoustical Society of America, 2012, Pp. 4324—4336.
11. Langtangen P.H., Linge S. Finite difference computing with PDEs. Texts in Computational Science and Engineering, 2017, Pp. 93—205.
12. Li C., Duric N., Huang L. Clinical breast imaging using sound-speed reconstructions of ultrasound tomography data. Medical Imaging Int. Society for Optics and Photonics, 2008, Pp. 692009-1—692009-9.
13. Pratt R.G. Seismic waveform inversion in the frequency domain. Part 1: Theory and verification in a physical scale model. Geophysics, 1999, Vol. 64, Pp. 888—901.
14. Huang L., Williamson M. Detection of breast microcalcifications with super-resolution ultrasound imaging: A clinical study. SPIE Medical Imaging, 2013, Pp. 867510-1—867510-10.
15. Li J. Kirchhoff common offset migration velocity analysis via differential semblance. Master's Thesis, Rice University, 2007, Pp. 1—61.
16. Goncharskiy A.V., Romanov S.Y., Seryozh-nikov S.Y. Obratnyye zadachi posloynoy ultrazvu-kovoy tomografii s dannymi na tsilindricheskoy poverkhnosti [Inverse problems of layer-by-layer ultrasonic tomography with the data measured on a cylindrical surface]. Vychislitelnyye metody i pro-grammirovaniye [Numerical methods and programming], 2017, Vol. 18, Pp. 267—276. (rus)
СВЕДЕНИЯ ОБ АВТОРАХ / THE AUTHORS
ВЕРЕЩАГИН Константин Алексеевич VERESHCHAGIN Konstantin A.
E-mail: kostya.veresh@mail.ru
БЕЛЫХ Игорь Николаевич BELYKH Igor N.
E-mail: ibelykh.spb@gmail.com
© Санкт-Петербургский политехнический университет Петра Великого, 2018
