CC BY
59
Уфа : УГАТУ, 2013
Ъыьмт QjrAQnQj
Т. 17, № 5 (58). С. 260-269
удк 519.6; 517.958:5
Суперкомпьютерное моделирование в задаче ультразвуковой диагностики
с применением аналитических подходов
1 -) Г. М. Агаян1, С. Ю. Романов2
1 agagal@rambler.ru, 2 romanov60@gmail.com
Московский государственный университет имени М. В. Ломоносова (МГУ) Научно-исследовательский вычислительный центр МГУ (НИВЦ МГУ)
Поступила в редакцию 21.09.2013
Аннотация. Работа посвящена изучению распространения ультразвукового излучения в среде и решению прямой и обратной задачи ультразвуковой томографии. Взаимодействие излучения с неоднород-ностями среды моделируется двумя методами: конечно-разностным и аналитическим. Это позволяет оценить границы применимости рассматриваемых моделей и оценить точность вычислений. Использование двух независимых методов позволяет удостовериться в надежности методов решения прямой и обратной задачи. Методами математического моделирования исследовано влияние плотности вещества на возможность реконструкции. Выбор параметров моделей ориентирован на задачу дифференциальной диагностики заболеваний молочной железы. Использованные алгоритмы базируются на прямом вычислении градиента функционала невязки. Проблема большого объема вычислений при решении обратной задачи преодолевается использованием суперкомпьютеров кластерного типа на основе технологии MPI. Используемые явные конечноразностные схемы идеально подходят для распараллеливания. Приведены результаты модельных расчетов, демонстрирующие эффективность применения предлагаемых подходов.
Ключевые слова: коэффициентные обратные задачи; волновое уравнение; компьютерное моделирование; ультразвуковая томография; параллельные вычисления; суперкомпьютер.
ВВЕДЕНИЕ
В настоящее время в публикациях отмечается большой интерес к задачам ультразвуковой диагностики. Как правило, работы нацелены на решение двумерных задач восстановления двумерных сечений трехмерных диагностируемых объектов. Предлагаемая работа также посвящена изучению распространения ультразвукового излучения в двумерной среде и решению прямой и обратной двумерной задачи ультразвуковой томографии. Взаимодействие излучения с неоднородностями среды моделируется двумя различными подходами. Это позволяет оценить границы применимости рассматриваемых моделей и оценить точность вычислений для различных значений параметров. Использование двух
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, проект № 13-07-00824 А.
независимых методов расчета распространения волн позволяет удостовериться в надежности методов решения прямой и обратной задачи. Одним из параметров, влияющих на распространение ультразвуковых волн, является плотность среды. В настоящей статье методами математического моделирования исследовано влияние плотности на возможность реконструкции неоднородностей для медицинских сред, имеющих вариацию плотности в 10-20 %.
Отметим, что значительный интерес в настоящее время представляет создание ультразвуковых томографов [1-3] для дифференциальной диагностики заболеваний молочной железы. Выбор параметров математических моделей в настоящей работе ориентирован именно на решение этой важной практической задачи.
Для решения прямой и обратной задачи ультразвуковой томографии в общем случае нами использовались конечно-разностные методы. В случае простой неоднородности в виде бесконечного цилиндра, помещенного в однородную
среду, поле прямой задачи рассеяния можно найти аналитически в виде рядов по специальным функциям. Этот подход хорошо исследован в литературе [4-6]. Используя разработанную в этих публикациях технологию, можно рассчитать волновое поле в любой точке пространства для широкого класса форм импульса через ряды. Полученное аналитическое решение использовалось для сравнения с конечно-разностным решением, а также в качестве данных прямой задачи для реконструкции неоднородности. Решив обратную задачу по этим данным, становится возможным сравнить полученный результат с известным точным решением.
Решение обратных задач ультразвуковой томографии как коэффициентных обратных задач с высокой точностью само по себе является сложной задачей. Использованные алгоритмы базируются на прямом вычислении градиента функционала невязки [7-10]. Естественно, что математическая модель должна учитывать явление дифракции, рефракции, переотражения излучения от границ и т. п.
Проблема большого объема вычислений при решении обратной задачи преодолевается использованием суперкомпьютеров кластерного типа на основе технологии MPI. Используемые явные конечноразностные схемы идеально подходят для распараллеливания.
ВЫВОД СКАЛЯРНОГО УРАВНЕНИЯ ЛИНЕЙНОЙ АКУСТИКИ В НЕОДНОРОДНОЙ СРЕДЕ
Уравнение распространения волны в неоднородной среде хорошо известно и может быть получено из следующих соображений [11]. Будем считать, что движение среды описывается адиабатическим движением идеальной жидкости. При адиабатичном движении энтропия, приведенная к единице массы, каждого перемещающегося в пространстве участка остается постоянной по времени, т. е.
^=о,
dt
(1)
где берется полная производная по времени от энтропии.
Будем также считать, что движение среды описывается малыми колебаниями сжимаемой жидкости. Т.е. колебания плотности р(г,£), давления р(г, £) , вектора скорости у(г, £) -малы. Представим функции плотности и давления в виде
p(r, t)=p0(r)+pl(r,t), p(r, t)=Po(r)+Pl(r, t) , (2)
где p0 (r), р0 (r) - постоянные по времени равновесные значения, p (r, t), р (r, t) - малые колебания во времени (p0 >> px, р0 >>рх). pj, Pj, v - малые величины первого порядка. Равновесное значение скорости будем считать равным 0. Заметим что при отсутствии внешнего поля p0 (r) = p0 = const.
Пренебрегая малыми величинами второго порядка, уравнение движения Эйлера запишется в виде
^VM+(v(r, t) • Vv(r, t + =
dty ' p(r)
dv(r,t ) + V( Ро + Pi(r,t )) dt Po(r) + Pi (r, t) dv(r, t) Vp (r, t )
(3)
dt
Po(r )
• = 0.
Колебания в среде также удовлетворяют уравнению непрерывности
dp(r, t ) dt
+p(r,t)div v(r,t) = 0 .
(4)
Таким образом, колебания в среде подчиняются уравнениям (1), (3), (4).
В силу уравнения (1) для уравнения состояния имеем
dp ^ dp 1 dp 1 (Ср dt
dp
dt c2 dt c2 l dt
- + v • Vp
где c =
( dp^
\ r Js
ности приводится к виду
dp оди
1 f dp d
. Тогда уравнение непрерыв-
+ v • Vp I + p(V • v) = 0.
(5)
Подставим выражения (2) в (5), учитывая p0 (r) = p0 = const, имеем
d( po + pi) dt
+ v• V(p0 + pi) 1 + (P0 +Pi)c (V^v) =
= fdp1 + v •Vpi J + P0c 2(V^ v) + Pic 2(V^ v) =
= dpL + P0c 2(V • v) = 0, dt
(6)
где пренебрегли величинами второго порядка малости V • Ур1 и р1с2 (V • у). Далее берем
s
производную по t в (6) и подставляем уравнение Эйлера (3)
1 82 p1 8v
—^тт + Po(V'—) =
с2 8tl
8t
1 82 pi
2 Я2 +Po(V'(—)) =
с2 8t P0
2
1 82 p
T -^L - Po (Vp-1 • Vpi) - pop-1 (V • Vpi) =
г 8t
1 82 p (VPo-VA)
с2 8t2
Po
-Ap = o.
(7)
Таким образом, получили скалярное уравнение линейной акустики в среде с переменной плотностью р0(г) и фазовой скоростью с(г). В случае слабо меняющейся плотности Ур0 (г) мало и, пренебрегая
№ -Ур)
величиной второго порядка малости
Po
получаем стандартное волновое уравнение, которое при наличии точечного источника, генерирующего импульс /(г) , имеет вид (8).
ПОСТАНОВКА И МЕТОДЫ ЧИСЛЕННОГО РЕШЕНИЯ ПРЯМОЙ И ОБРАТНОЙ ЗАДАЧИ
Существуют различные постановки коэффициентных обратных задач для волнового уравнения [12-16]. В настоящей работе обратная задача рассматривается в следующей постановке. Рассмотрим волновое уравнение, которое описывает поле р(г, г) в течение
времени (0, Т) в области Р с Ям (N=2, 3), ограниченной поверхностью 5", с точечным источником, располагающимся в точке г
у(г)р„ (г, г) - Ар(г, г) = 8(г - Го) - /(г), (8) р(г, г = 0) = р (г, г = 0) = о,
д пР =q(г, г). Здесь у 05(г) = с(г) является скоростью волны
в среде, г е Яы - положение точки в пространстве, А - оператор Лапласа по переменной г. Генерируемый источником импульс описывается функцией /(г) , дпр \5Г -производная вдоль нормали к поверхности 5 в области £х(о,Г). Ц_(г,г) - некоторая известная функция, которая может быть найдена по экспериментальным данным на границе 5 из решения внешней для области Р задачи. Будем
предполагать, что неоднородность среды вызвана только изменениями скорости, а вне области неоднородности c(r) = const, где const известна.
Обратная задача состоит в нахождении функции c(r), описывающей неоднородность, по экспериментальным данным измерения волны U(5, t) на границе S области P за время
(o,T) при различных положениях ro источника.
Введем функционал невязки
ф(p(v)) = 1J j(p(s, t) - U(s, t))2dsdt.
o S
Здесь и (я, г) - экспериментальные данные на границе £ области Р за время (0,Т). В рамках настоящей постановки выражения для градиента функционала невязки Ф(р(у)) получены в работах [7-10, 16]. Градиент Ф'(р(с)) имеет вид
ФV (p(v), v) = J Wt (r, t)pt (r, t)dt .
(9)
Здесь р(г,г) - решение основной задачи (8), а ^(г, г) - решение следующей «сопряженной» задачи (10) при заданном с(г)
у(г)щ (г, г) -Ам>(г, г) = 0, (10)
^(г, г = т) = ^ (г, г = т) = 0,
дп™ \зт = Р\зт -и .
Таким образом, для вычисления градиента функционала необходимо решить основную и «сопряженную» задачи. Отметим еще раз, что методы решения прямой и обратной задачи линейной акустики приведены в предположении, что плотность среды, входящая в общее уравнение (7), слабо меняется или постоянна вне и внутри неоднородности. Численные методы решения прямой и обратной задачи описаны в работе [10].
Для решения трехмерных обратных задач широко распространен послойный подход, в котором трехмерный объект разбивается на двумерные слои с некоторым шагом и в пределах слоя решаются независимые двумерные задачи. В дальнейшем будем рассматривать двумерную прямую и обратную задачу в квадратной области Q, в центре которой располагается искомая неоднородность в виде круга, обозначенная цифрой 1 на рис. 1. Значения скорости волны внутри и вне круга будем считать постоянными, но не равными. Неоднородность облу-
o
чается с четырех сторон плоскими зондирующими импульсами, обозначенными стрелками. По периметру квадрата располагаются приемники излучения, обозначенные цифрой 2.
2
Рис. 1. Схема численного эксперимента в двумерном случае
Рассмотрим также задачу распространения и рассеяния волны в трехмерном пространстве для неоднородности в виде цилиндра и плоского падающего перпендикулярно оси цилиндра зондирующего импульса (рис. 2). Если ось цилиндра направлена вдоль оси У, то уравнения задачи не зависят от переменной у, и поэтому задача эквивалентна прямой двумерной задаче распространения волны в рассмотренном выше случае на рис. 1. В трехмерном случае для цилиндра известны аналитические методы решения прямой задачи. Поэтому открывается возможность решения прямой задачи как численно, например, методом конечных разностей, так и независимо аналитически. Это позволяет сравнить результаты решения прямой задачи, а также проверить возможность решения обратной задачи на основе независимых данных, полученных аналитически.
ПОСТАНОВКА И АНАЛИТИЧЕСКИЕ РЕШЕНИЯ ПРЯМОЙ ЗАДАЧИ РАССЕЯНИЯ ДЛЯ ЦИЛИНДРИЧЕСКОЙ НЕОДНОРОДНОСТИ
Рассмотрим поставленную выше прямую задачу в случае неоднородности в виде бесконечного цилиндра. Выпишем аналитическое решение задачи в виде рядов по специальным функциям для этого случая. Для волн, гармонических по времени вида
p(r, t) = Re U(r) • е'ш },
где ю - частота колебаний, волновое уравнение сводится к уравнению Гельмгольца относительно комплекснозначной функции и(г)
Аи + к и = 0, 1 ю
где волновое число к = —. В трехмерном случае
с
решения уравнения Гельмгольца удовлетворяют условию излучения Зоммерфельда, если
im |г| • (—^ + iku) = 0 . ^сю d\r\
(11)
Здесь |r| = t[x2 + y2 + z2 и
предел берется
равномерно по всем направлениям. Если рассмотреть два линейно независимых сферически симметричных решения уравнения
Гельмгольца, u ='
ik|r
,k\r
только u
удовлетворяет условию излучения. Этому решению уравнения Гельмгольца соответствует
Р1(г, £) = Яе^) • }= ^Ю^
имеющий вид
расходящейся волны.
Рассмотрим простейшую задачу рассеяния
на теле Ос Я падающей звуковой волны. Предположим, что внутри тела Ос Я, имеющего границу 5", скорость распространения зву-
ка
равна c(r) = cQ = const
и
плотность
p(r) = pQ = const. Вне тела О скорость звука равна c(r) = c0 = const и плотность равна p(r) = р0 = const. Предположим, что падающее
поле и' (г) ( г i q с r3) представляет собой плоскую волну, распространяющуюся вдоль оси OZ: и' (r) = e~'kz = е cosS . Полное поле вне области Q (r i Q с R3) обозначим через и (г). Оно представляет собой сумму падающей волны и неизвестного рассеянного поля us(r), т. е. u(r) = и' (r) + us (r), где r i Q с R3. Эти функции вне области Q удовлетворяют уравнению
Гельмгольца с волновым числом k = — . Внут-
c0
ри области Q (r е Q с R3) полное поле обозначим через иQ (г). Оно удовлетворяют уравнению Гельмгольца внутри области Q с волновым числом kQ = — . На границе S области Q
cQ
поле и (г) и поле и (г) удовлетворяют условиям сопряжения
r
е
и и2 =
г
r
□
и = и и
ди_ _1_ р0 ду Рп
ди р ду
(12)
д
где — - производная по внешней нормали на
ду
границе 5.
Задача рассеяния на теле □ с Я3, состоит в нахождении рассеянного поля и* (г) вне тела
□ с Я3, удовлетворяющего условиям излучения Зоммерфельда (11) на бесконечности.
Для областей □, имеющих простую форму, решение поставленной задачи может быть выписано в аналитическом виде [4] . Рассмотрим
□ в виде бесконечного цилиндра радиуса Я, см. рис. 2. Введем цилиндрические координаты, т. ч. точка с координатами (х, у, ¿) будет иметь координаты (| г |,у,0).
у*
/ 1,1 иу
о
Рис. 2. Схема численного эксперимента в трехмерном случае
В нашем случае падающая плоская волна не зависит от координаты у и может быть представлена в виде ряда
и' ( г|, 0) = е~*Н С086=Х5И-(~г)т-Тт-Зт (к|г|),
т>0
Г1, т = 0
где Тт = оо8(т0), 5т = \ , J„(x) - функция
[2, т > 0
Бесселя первого рода.
Внутри области □ решение уравнения Гельмгольца, удовлетворяющее условию сопряжения (12), имеет вид
и □ (I г
(| г|, 9) = Х Р -5т -(-'Т-Тт^т (¿Ц) ,
где
^ (кЯ) - Нт (кЯ) - Зт (кЯ) - Нт (кЯ) Зт (к„Я) - Нт (кЯ) - С-3'т (каЯ) - Нт (кЯ)
Решение поставленной задачи рассеяния волна и*(г) (которая является решением уравнения Гельмгольца вне области □, удовлетворяющая условию излучения (11) и такая, что
и' (г) + и* (г) удовлетворяет условию сопряжения (12)) будет иметь вид
и* (|г|, 9) = X Вт -5т -(-Т -Тт -Нт (к|г|) , (13)
т>0
где
В =
С - Зт (кЯ) - Зт (кр Я) - Зт (кЯ) - ^ (кр Я) Зт (кр Я) - Нт (кЯ) - С - Зт (кр Я) - Нт (кЯ)
где
С =
кПр0 _ с0р0
кРп
сррр
Нт (х) - функция Ханкеля второго рода.
Коэффициент С равняется отношению акустических сопротивлений внутри неоднородности и вне ее. Он определяет величину и фазу отраженной волны.
Предположим теперь, что падающее поле представляет собой плоскую волну в виде короткого импульса, распространяющегося вдоль оси 02\
у
р' (г, г) = / (г--),
С0
где функция /(г) отлична от нуля на интервале (г*,г* + х), где х - длительность импульса. Будем рассматривать функцию /(г) в виде
конечного ряда Фурье на отрезке (-Т Т), где Т выбирается из условия полного прохождения импульса интересующей нас области, ( = 0 соответствует началу прохождения импульса через точку х = у = 2 = 0,
N
/(г) = X fn ехр('Юпг),
®п ~ . (14)
Пользуясь линейностью уравнения Гельм-гольца и условий сопряжения, можно решать
а,
задачу рассеяния по формулам (13) (ки = —^)
С0
отдельно для каждой гармоники ши соотношения (14) и получать решения и" (г,шп) . Общим
п=—N
решением поставленной задачи на интервале (-T, T), будет сумма решений, полученных для рассматриваемых гармоник. Заметим, что, так как f (t) действительная функция, проводить расчеты можно только для половины слагаемых. Полученное решение является периодической функцией по времени с периодом 2T. На интервале (-T, T) p(r, t) будет совпадать с искомым решением. Сходимость рядов, представляющих решение задачи рассеяния, исследована, например, в [5].
ЧИСЛЕННЫЕ РАСЧЕТЫ НА СУПЕРКОМПЬЮТЕРЕ
Расчеты прямой задачи на основе аналитических формул (13) проводились в однопроцессорном варианте, поскольку требуется найти волновое поле только в точках границы квадрата и в вычислительном плане эта задача не является объемной. Решение обратных задач проводилось на основе формул (8)-(10) с использованием экспериментальных данных (результатов решения прямой задачи), полученных как ко-нечноразностным методом на основе формул (8), так и независимым методом расчетов на основе аналитических формул (13). Расчеты обратных задач с помощью разработанных ранее алгоритмов и программ [17-19] были проведены на суперкомпьютерах СКЦ МГУ «Чебышёв» и «Ломоносов» с помощью технологии MPI.
Волновое уравнение (8) аппроксимировалось явной трехслойной по времени разностной схемой со вторым порядком точности. Явная схема сравнительно легко распараллеливается на многопроцессорной системе, например, с помощью метода распараллеливания по пространству, поскольку вычисление значения в точке на новом слое по времени зависит только от ближайших соседних точек на текущем и предыдущем слоях. Предлагаемый способ распараллеливания по пространству состоит в том, что общее квадратное поле вычислений разбивается на одинаковые квадратные подобласти, вычисления в которых производятся различными вычислительными ядрами. Было выделено 9 = 3 х 3 подобластей, по границам которых на каждом шаге по времени выполнялся обмен данными.
Следующий уровень распараллеливания связан с тем, что расчеты надо проводить для четырех источников, вычисления для каждого из которых в значительной степени независимы. Таким образом, обратная задача решалась на 36
узлах = 4 источника х 3 х 3 подобластей. Ускорение по сравнению с однопоточным вариантом составило около 25 раз. Время расчета одной обратной задачи - около 20 ч.
Использование численных методов для расчета распространения волн по необходимости предполагает использование ограниченных по пространству областей расчетов, что противоречит реальной постановке задачи. Поэтому приходится ставить искусственные граничные условия по периметру области расчетов. Ранее авторы ставили условия «неотражения» на границе, что выполняется лишь приближенно с некоторой погрешностью. Существуют и другие методы устранения эффектов, связанных с ограничением области расчетов, например, метод PML или его модификация CFS-PML. В настоящей работе ставилась задача выявления возможных погрешностей, вносимых в расчеты только разностными методами, поэтому эффекты, связанные с границей, были полностью исключены простым способом — к исходной области расчетов были добавлены дополнительные буферные полосы с каждой стороны, т. е. границы были отодвинуты. Ширина буферных полос составила 500 точек и выбиралась из условия, чтобы волна, пересекшая наружу границу области, не успела за время T из формулы (10) решения прямой задачи пересечь границу обратно. Это позволило полностью исключить влияние граничных эффектов на внутренних границах, представляющих интерес. Проблема увеличения общего объема расчетов (размер области расчетов составил 2000 х 2000 точек) преодолевалась увеличением числа вычислительных ядер. Такой подход заметно увеличивает объем вычислений и представляется оправданным только в данном случае проведения небольшого числа модельных расчетов.
Одна из проблем, возникших при расчетах — это нехватка памяти для хранения данных, поскольку расчеты одним ядром проводятся на области 667 х 667 точек, а экспериментальные данные требуют хранить 6000 шагов по времени. Для уменьшения объема данных экспериментальные данные были переведены в формат float. Тестовые расчеты не показали заметного ухудшения точности расчетов.
ЧИСЛЕННЫЙ ЭКСПЕРИМЕНТ
Рассмотрим неоднородность, имеющую вид бесконечного цилиндра, с осью вдоль OY и диаметром 12 см. Скорость и плотность внутри ци-
линдра постоянны, са = 1600 м/с. В качестве источника излучения использовалась плоская волна, падающая на неоднородность перпендикулярно оси цилиндра (рис. 2). Выделив приемники, расположенные на сторонах квадрата 20 х 20 см вокруг неоднородности в плоскости, перпендикулярной оси цилиндра, в них были проведены расчеты волнового поля (прямой задачи) по формуле (13). Скорость распространения звука вне неоднородности равнялась 1500 м/с. Размер шага сетки выбирался достаточным для хорошей аппроксимации распространения импульса с длиной волны 7 мм явной разностной схемой, что привело к размеру сетки расчетов 1000 х 1000 точек для рассматриваемого квадрата.
Обратная задача решалась по данным зондирования неоднородности с четырех перпендикулярных направлений плоской волной (четыре источника). Фронт падающей волны параллелен сторонам квадрата. На рис. 3 изображен график короткого импульса в зависимости от времени, который рассматривался в качестве зондирующего.
Рис. 3. Форма короткого зондирующего импульса
Одна из задач, стоявшая в работе, - оценить точность вычислений на основе конечноразно-стных схем. Для этого решалась прямая задача конечноразностным методом в двумерном случае и независимым методом по аналитическим формулам (13) эквивалентная трехмерная задача для цилиндра. На рис. 4, а приведены графики поля на дальней по направлению падающей волны границе, полученные численными и аналитическими методами для случая одинаковой плотности внутри и вне цилиндра. Видно очень хорошее согласие в расчетах.
Другая задача, которая ставилась в работе, проверить надежность методов решения обратной задачи. Для этого обратная задача решалась на основе модельных экспериментальных дан-
ных, полученных независимым методом по аналитическим формулам. На рис. 4, б приведены для сравнения результаты реконструкции по данным численного решения прямой задачи, а на рис. 4, в - по данным аналитического решения прямой задачи. Не удивительно, что результаты хорошо согласуются, так как результаты решения прямой задачи практически совпали.
Как следует из формулы (7), одним из параметров, влияющих на распространение ультразвуковых волн, является плотность среды. Методами математического моделирования исследовано влияние плотности на возможность реконструкции неоднородностей для медицинских сред, имеющих вариацию плотности в 10-20 %. На рис. 5, а приведены графики поля на боковой по направлению падающей волны границе, полученные аналитическими методами для случаев: 1) плотности внутри и вне цилиндра равны, акустическое сопротивление внутри больше еара > с0р0 (зеленый пунктир); 2) плотность внутри меньше, чем вне цилиндра, акустические сопротивления равны еа ра = е0 р0 (синие точки); 3) плотность и акустическое сопротивление внутри меньше, чем вне цилиндра сара < с0р0 (красная сплошная). Видно, что отраженная волна имеет разную фазу, а при равных акустических сопротивлениях отражение исчезает при приближении к нормальному углу падения к поверхности цилиндра. Однако преломленные волны (крайние справа на графике) имеют амплитуду в 10-20 раз большую и практически совпадают. На рис. 5, б и в приведены результаты реконструкции методами (8)-(10) по экспериментальным данным, полученным из аналитических формул для разных плотностей внутри и снаружи цилиндра. На рис. 5, б приведены результаты при равных акустических сопротивлениях, на рис. 5, в — когда акустическое сопротивление внутри меньше, чем вне цилиндра, на рис. 4, в — когда плотности внутри и вне цилиндра равны, акустическое сопротивление внутри больше. Видно, что во всех трех случаях получаются похожие результаты, что связано с тем, что преломленные волны, в которых содержится основная часть энергии, практически совпадают. На рис. 5, в появилась дополнительная граница круга, так как отражение в проти-вофазе, а на рис. 5, б — граница размыта, так как отражение слабое.
а б в
Рис. 4. Решение задачи в предположении, что плотность среды одинакова во всем пространстве: а - графики волны, полученные численным и аналитическим решением прямой задачи;
б - поле скоростей, восстановленное по данным численного решения прямой задачи; в - поле скоростей, восстановленное по данным аналитического решения прямой задачи
а б в
Рис. 5. Решение задачи в предположении, что плотности внутри и вне цилиндра различны: а - вид отраженной и преломленной волн для рр = 1 г/см3 (зеленый пунктир), рр = 0.9375 г/см3 (синие точки), рр = 0.8 г/см3 (красная сплошная); б - восстановленное поле скоростей для ра = 0.9375 г/см3; в - восстановленное поле скоростей для рп = 0.8 г/см3
ВЫВОДЫ
Решение прямой задачи в двумерном случае в аналитической и в конечно-разностной постановках показали хорошее согласие в расчетах волнового поля в неоднородной среде.
В предположении, что плотность среды одинакова для неоднородности и окружающего пространства, решения задачи восстановления поля скоростей для входных данных, полученных аналитически и численно, практически не отличаются. Полученные результаты являются независимым тестом программного обеспечения для решения обратных задач УЗ томографии.
Проведенные модельные расчеты в области параметров, близких к ультразвуковой диагностике заболеваний молочной железы, показали, что восстановление поля скоростей на основе волновой модели (не учитывающей плотность) не сильно чувствительно к наличию разницы в 10-20 % в плотности вещества. Удается с небольшими артефактами восстанавливать как
форму неоднородного объекта, так и абсолютные значения скорости распространения волны. Это оправдывает применение волновой модели в задаче ультразвуковой томографии.
СПИСОК ЛИТЕРАТУРЫ
1. Wiskin J., Borup D.T., Johnson S.A., and Berggren M.
Non-linear inverse scattering: High resolution quantitative breast tissue tomography // J. Acoust. Soc. Am. 2011. Vol. 131, Iss. 5. P. 3802-3813.
2. Glide-Hurst C. K., Duric N., Littrup P. Volumetric breast density evaluation from ultrasound tomography images // Medical Physics. 2008. 35. P. 3988-3997.
3. Jirík R., Peterlík I., Ruiter N., Fousek J., Dapp R., Zapf M., Jan J. Sound-speed image reconstruction in sparse-aperture 3-D ultrasound transmission tomography // IEEE Trans. on Ultrasonics, Ferroelectrics, and Frequency Control. 2012. 59, No. 2. P. 254-264.
4. Varadan V. V., Ma Y., Varadan V. K., Lakhtakia A. Scattering of waves by spheres and cylinders // Field representations and Introduction to Scattering. Amsterdam: North-Holland, 1991. P. 211-324.
5. Korneev V. A., Johnson L. R. Scattering of elastic waves by a spherical inclusion. 1. Theory and numerical results // Geophys. J. Int. 1993. No. 115. P. 230-250.
6. Pike R., Sabatier P. Scattering: scattering and inverse scattering in pure and applied science // San Diego: Academic Press, 2008. p. 2002. ISBN: 0126137625, ISBN13: 9780126137620.
7. Chavent G. Deux resultats sur le probleme inverse dans les equations aux derivees partielles du deuxieme ordre an t et sur I'unicite de la solution du probleme inverse de la diffusion // Paris. C. R. Acad. Sc. 1970. Vol. 270. P. 2528.
8. Natterer F., Wubbeling F. A propagation-backpropagation method for uzltrasound tomography // Inverse Problems. IOP Publishing Ltd, 1995. No. 11. P. 1225-1232.
9. Beilina L., Klibanov M. V. Approximate global convergence and adaptivity for coefficient inverse problems. New York: Springer, 2012.
10. Goncharsky A. V., Romanov S. Y. Supercomputer technologies in inverse problems of ultrasound tomography // Inverse Probl. 2013. Vol. 29, № 7. 075004.
11. Ландау Л. Д., Лифшиц Е. М. Гидродинамика, М.: Наука, 1986.
12. Backushinsky A., Goncharsky A., Romanov S., Seatzu S.
On the identification of velocity in seismics and in acoustic sounding // Pubblicazioni dell'istituto di analisa globale e applicazioni. Serie "Problemi non ben posti ed inversi». Firenze, November 1994. No. 71.
13. Головина С. Г., Романов С. Ю., Степанов В. В. Об одной обратной задаче сейсмики // Вестн. МГУ. Сер. 15. Выч. мат. и киб. 1994. № 4. С. 16-21.
14. Гончарский А. В., Романов С. Ю. Об одной трехмерной задаче диагностики в волновом приближении // ЖВМиМФ. 2000. Т. 40, № 9. С. 13641367.
15. Бакушинский А. Б., Козлов А. И., Кокурин М. Ю. Об
одной обратной задаче для трехмерного волнового уравнения // Журн. вычисл. матем. и матем. физики. 2003. 43, № 8. С. 1201-1209.
16. Гончарский А. В., Романов С. Ю. О двух подходах к решению коэффициентных обратных задач для волновых уравнений // Журн. вычисл. матем. и матем. физики. 2012. 52, № 2. С. 1-7.
17. Гончарский А. В., Романов С. Ю. Суперкомпьютерные технологии в разработке методов решения обратных задач в УЗИ-томографии // Вычислительные методы и программирование. 2012. Т. 13. С. 235-238.
18. Воеводин Вад. В, Овчинников С. Л., Романов С. Ю. Разработка высокоэффективных масштабируемых программ в задаче ультразвуковой томографии // Вычислительные методы и программирование. 2012. Т. 13, № 1. С. 307-315.
19. Романов С. Ю. К вопросу об масштабируемости программы для обратной задачи волновой томографии // Вестник Нижегородского университета им. Н. И. Лобачевского. 2013. № 2-1. С. 160-167.
ОБ АВТОРАХ
АГАЯН Галина Михайловна, науч. сотр. НИВЦ. Дипл.
прикл. математик (МГУ им. М. В. Ломоносова, 1981). Канд.
физ.-мат. наук (там же, 1986). Иссл. в обл. некорректно
поставленных задач.
РОМАНОВ Сергей Юрьевич, вед. науч. сотр. НИВЦ. Дипл.
математик (МГУ им. М. В. Ломоносова, 1982). Канд. физ.-
мат. наук (там же, 1986). Иссл. в обл. обратных задач математической физики.
METADATA
Title: Supercomputer simulations in the problem of ultrasound diagnostics with the use of analytical approaches.
Authors: G. M. Agayan and S. Y. Romanov
Affiliation: Research Computing Center of the M. V. Lomonosov Moscow State University, Russia.
Email: agagal@rambler.ru, romanov60@gmail.com.
Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 17, no. 5 (58), pp. 260-270, 2013. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print).
Abstract: This is a study of ultrasound radiation in the environment and the solution of the direct and inverse problems of ultrasound tomography. The interaction of radiation with inhomogeneities is modeled by two methods: finite-difference and analytic. This allows us to estimate the limits of applicability of the considered models and evaluate the accuracy of the calculations. The use of two independent methods helps to check the reliability of methods for solving direct and inverse problems. Mathematical modeling methods investigated the effects of the density of the substance to the possibility of reconstruction. Choice of model parameters are oriented to the problem of differential diagnosis of breast diseases. The used algorithms are based on the direct calculation of the gradient of the residual functional. The problem of the large amount of computation for solving the inverse problem is overcome by using a supercomputer cluster-type on the basis of MPI technology. Used explicit finite difference scheme is ideally suited for parallelization. The results of model calculations are shown that demonstrate the effectiveness of the proposed approaches.
Key words: coefficient inverse problems; wave equation; computer modeling; ultrasonic tomography; parallel computing; supercomputer.
References (English Transliteration):
1. J. Wiskin, D. T. Borup, S. A. Johnson, and M. Berggren, "Non-linear inverse scattering: High resolution quantitative breast tissue tomography," J. Acoust. Soc. Am. vol. 131, iss. 5, pp. 3802-3813, 2011.
2. C. K. Glide-Hurst, N. Duric, and P. Littrup, "Volumetric breast density evaluation from ultrasound tomography images," Medical Physics, 35, pp. 3988-3997, 2008.
3. R. Jink, I. Peterlik, N. Ruiter, J. Fousek, R. Dapp, M. Zapf, and J. Jan, "Sound-speed image reconstruction in sparse-aperture 3-D ultrasound transmission tomography," IEEE Trans. Ultrasonics, Ferroelectrics, and Frequency Control, 59, no. 2, pp. 254-264, 2012.
4. V. V. Varadan, Y. Ma, V. K. Varadan, and A. Lakhtakia, "Scattering of waves by spheres and cylinders," Field representations and Introduction to Scattering, pp. 211324, Amsterdam: North-Holland, 1991.
5. V. A. Korneev and L. R. Johnson, "Scattering of elastic waves by a spherical inclusion. 1. Theory and numerical results," Geophys. J. Int., no. 115, pp. 230-250, 1993.
6. R. Pike and P. Sabatier, Scattering: scattering and inverse scattering in pure and applied science, p. 2002, San Diego: Academic Press, 2008. ISBN: 0126137625, ISBN13: 9780126137620.
7. G. Chavent, "Deux resultats sur le probleme inverse dans les equations aux derivees partielles du deuxieme ordre an t et sur I'unicite de la solution du probleme inverse de la diffusion," C. R. Acad. Sc., vol. 270, pp. 25-28, Paris, 1970.
8. F. Natterer and F. Wubbeling A propagation-backpropagation method for uzltrasound tomography," Inverse Problems, no. 11, pp. 1225-1232, IOP Publishing Ltd., 1995.
9. L. Beilina and M. V. Klibanov, Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems. New York: Springer, 2012.
10. A. V. Goncharsky and S. Y. Romanov, "Supercomputer technologies in inverse problems of ultrasound tomography," Inverse Probl. vol. 29, no. 7, 075004, 2013.
11. L. D. Landau and E. M. Lifshits, Hydrodynamics, (in Russian). Moscow: Nauka, 1986.
12. A. Backushinsky, A. Goncharsky, S. Romanov, and S. Seatzu, "On the identification of velocity in seismics and in acoustic sounding," Pubblicazioni dell'istituto di analisa globale e applicazioni, Serie "Problemi non ben posti ed inverse", no. 71, Firenze, November 1994.
13. S. G. Golovina, S. Y. Romanov, and V. V. Stepanov, "About an inverse problem of seismology," (in Russian), Vestnik MGU (MSU Bulletin), ser. 15: Calculating Mathematics and Cybernetics, no. 4, pp.16-21, 1994.
14. A. V. Goncharsky and S. Y. Romanov, "On a three-dimensional diagnostic problem in the wave approximation," (in Russian), Computational Mathematics and Mathematical Physics, vol. 40, no. 9. pp. 1364-1367, 2000.
15. A. B. Backushinsky, A. I. Kozlov, and M. Y. Kokurin, "About an inverse problem for the three-dimensional wave equation," (in Russian), Computational Mathematics and Mathematical Physics, vol. 43, no. 8, pp. 1201-1209, 2003.
16. A. V. Goncharsky and S. Y. Romanov, "On two approaches to solving the coefficient inverse problem for the wave equation," (in Russian), Computational Mathematics and Mathematical Physics, vol. 52, no. 2, pp. 1-7, 2012.
17. A. V. Goncharsky and S. Y. Romanov, "Supercomputer technologies in the development of methods for the solution of inverse problems in ultrasound tomography," Numerical Methods and Programming, vol. 13, pp. 235238, 2012.
18. Vad. V. Voevodin, S. L. Ovchinnikov, and S. Y. Romanov, "Development of highly scalable programs in the problem of ultrasound tomography," Numerical Methods and Programming, vol. 13, no. 1, pp. 307-315, 2012.
19. S. Y. Romanov, "On the scalability of the program for the inverse task of wave tomography," Vestnik Nizhegorodskogo Universiteta imeni Lobachevskogo (Bulletin of the Nizhny Novgorod University), no. 2-1, pp. 160-167, 2013.
About authors:
AGAYAN, Galina Michailovna, Candidate of Physical and Mathematical Sciences, Scientific Researcher, Research Computing Center of the M. V. Lomonosov Moscow State University, Russia.
ROMANOV, Sergey Yurievich, Candidate of Physical and Mathematical Sciences, Leading Researcher, Research Computing Center of the M. V. Lomonosov Moscow State University, Russia.
