АКУСТИЧЕСКИЙ ЖУРНАЛ, 2008, том 54, № 3, с. 469-482
ОБРАБОТКА АКУСТИЧЕСКИХ СИГНАЛОВ ^^^^^^^^ И КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ
УДК 534.2:517.9
РЕШЕНИЕ ТРЕХМЕРНОЙ ОБРАТНОЙ ЗАДАЧИ АКУСТИЧЕСКОГО РАССЕЯНИЯ. МОДИФИЦИРОВАННЫЙ АЛГОРИТМ НОВИКОВА
© 2008 г. Н. В. Алексеенко, В. А. Буров, О. Д. Румянцева
Московский государственный университет им. М.В. Ломоносова, физический факультет
119992, ГСП-2, Москва, Ленинские горы Тел.: (495) 939-3081; Факс: (495) 932-8820 E-mail: burov@phys.msu.ru Поступила в редакцию 24.07.07 г.
Впервые осуществлено численное восстановление модельных трехмерных рассеивателей различной силы и размера на основе монохроматического функционально-аналитического алгоритма Новикова. Алгоритм учитывает процессы многократного рассеяния и не имеет жестких ограничений на силу рассеивателей. Получаемая оценка рассеивателя приближается к истинному значению после ограничения ширины пространственного спектра рассеивателя областью с радиусом около 2k0. Помехоустойчивость алгоритма к случайным ошибкам в экспериментальных данных достаточно высока для практических целей диагностики. Однако требуемое количество вычислительных операций велико.
PACS: 43.60.Pt , 43.35.Wa
1. ВВЕДЕНИЕ
Трехмерные системы акустического томогра-фирования в настоящее время находятся на начальном этапе развития как в техническом, так и в алгоритмическом плане [1]. Разрабатываемые практические системы - это объединение большого числа двумерных систем, так как обмен между слоями предполагается слабым и почти не учитывается [2]. В настоящей работе обсуждаются впервые полученные результаты численного моделирования решения обратной трехмерной задачи рассеяния акустических (скалярных) волн с помощью монохроматического функционально-аналитического алгоритма Новикова [3], учитывающего эффекты многократного рассеяния при постановке задачи, близкой к строгой. Достаточно большой список публикаций, относящихся к функционально-аналитическим методам, приведен в [4]. Поэтому далее приводятся только ссылки, непосредственно относящиеся к рассматриваемому алгоритму или к способу получения данных рассеяния для него. С точки зрения используемого подхода, алгоритм Новикова представляет собой модификацию предшествующего ему трехмерного алгоритма Новикова-Хен-кина [4].
Функция рассеивателя (рассеивающий потенциал) v(r, ю), восстанавливаемая при решении монохроматической обратной задачи рассеяния, характеризует неоднородность в виде отклонения скорости звука c(r) и амплитудного коэффициента поглощения или усиления a(r, ю) внутри обла-
сти рассеяния от их значений в фоновой среде. В случае однородного непоглощающего фона, в котором скорость звука равна с0 и волновое число
к0: у(г, ю) = к2 - к2с (г, ю). Здесь кс(г, ю) = к0п(г) + + гк0д(г, ю) - локальное комплексное волновое число, действительная часть которого определяется показателем преломления п(г) = с0/с(г), а мнимая часть а(г, ю) = к0д(г, ю) представляет собой амплитудный коэффициент поглощения (при а(г, ю) > 0) или усиления (при а(г, ю) < 0); временная зависимость поля полагается ~ ехр(-г'ю£). Тогда в терминах п и а функция рассеивателя приобретает вид:
у(г, ю) = к0{ 1 - п2(г) + q2(г, ю) - 2'п(г)ю)} = 2(1 1 ^ 2, ч 0. а(г, ю) (1)
---2—)+а( г'ю)-2 ' ю -Ътт-
с0 с (гу С(г)
Из (1) следует, что присутствие неоднородности скорости звука описывается за счет Яе V, но влияет также и на 1т V. В то же время, присутствие поглощения описывается, в первую очередь, членом 1т V, однако влияет и на Яе V в виде квадратичной поправки а2. За счет этого даже в случае отсутствия рефракции (п = 1) чисто поглощающий или усиливающий рассеиватель описывается функцией у(г, ю) = а2(г, ю) - 2г'юа(г, ю)/с0 с ненулевой действительной частью. С другой стороны, для рассеивателя, обладающего поглощением или усилением и, одновременно, фокусирующим показателем преломления п(г) > 1, возможна ситуация
/ 0 (рЦ / k-/ /k+ kR-P/2"--
Л® (p) ^^^
^^ k«
k2 = l2 = k0
при k, l e С3.
(2)
Если величина |к71 не бесконечно мала, то из
условий (2) следует, что к ± кд, кД - к] = к0; 1/ ±
± 1д, \Д - 12 = к0. В трехмерном случае дополнительно полагается, что к = 17, и тогда:
kR — k2 = k0, k7 ^ kR, k7 ^ lR,
iR =
R
k / = l /.
(3)
Рис. 1. Взаимная ориентация векторов, используемых в трехмерном алгоритме. Действительные компоненты волновых векторов кд и 1д образуют вектор р = кд - 1д; в плоскости, перпендикулярной к р, расположены альтернативные мнимые компоненты
волновых векторов к/ = к± и взаимно ортогональные орты 0(р), ю(р).
д2(г, ю) = я2(г) - 1, т.е. а2(г, ю) = ю2(с~2(г) - с- ), приводящая к Яе у(г) = 0. Однако часто в акустических задачах поправка а2 пренебрежимо мала, и тогда (1) переходит в обычно используемое выражение у(г, ю) = ю2(с-2 - с~2(г)) - 2г'юа(г, ю)/с(г).
Экспериментальными данными для трехмерного алгоритма Новикова является классическая амплитуда рассеяния /(к, 1), пропорциональная асимптотике рассеянного поля в дальней зоне [4];
к, 1 е К3 - волновые векторы падающей плоской монохроматической волны и рассеянной волны в
дальней зоне, соответственно; к2 = 12 = к2. Эти данные могут как непосредственно измеряться в эксперименте, так и синтезироваться из данных рассеяния, зарегистрированных в ближней зоне. Например, одна из возможных томографических схем сбора информации о трехмерном объекте предлагается в [1].
Функционально-аналитические методы решения обратной задачи рассеяния основаны на формальном распространении волновых векторов на область комплексных значений к, 1 е С : к = кд + г'к/, 1 = 1д + И1, где кд = Яе к, к = 1т к, 1д = = Яе 1, 1/ = 1т 1. Направление мнимой части волнового вектора к/ = |к/ |у можно характеризовать
единичным вектором у е К3, |у | = 1. Одновременно требуется сохранение условия
Таким образом, мнимые компоненты 1/ = к/ Ф 0 ортогональны к плоскости, в которой лежат векторы кд и 1д, образующие действительный вектор р = к - 1 = кд - 1д и, следовательно, они ортогональны вектору р (рис.1). Поскольку к/ = 1/ и к/ ± р, то проекция любого вектора к, образующего фиксированный вектор р, на направление вектора р всегда действительна и составляет р/2. Следовательно, соотношения
22 k = ko,
2 kp = p2
(4)
эквивалентны требованиям (3). Проекция вектора к на плоскость, перпендикулярную р, составляет к - р/2, и квадрат проекции равен (к - р/2)2 =
= к2 - р2/4 в силу (4). Возможны две альтернативные, противоположные ориентации направляющего вектора у = у± мнимой компоненты к/ = к± (величина которой не бесконечно мала), лежа-
р х(к - р/2) _
щей в этой плоскости: Y±(k, p) = ± -
Pl • k - p/2|
= ±
P х(kr - p/2)
. Геометрическая ситуация по-
IpI • 4kl-IpI2/4 яснена на рис. 1.
Следует подчеркнуть, что бесконечно малая, но ориентированная вдоль направления у, мнимая добавка k —► 0 не обязана быть перпендикулярной к kR. В этом случае обобщенная амплитуда рассеяния h(k, l) переходит в предельные функции hY(k, l) = hY(kR, lR) = lim h (kR + i\kj|у, lR), зави-
N ^ o
сящие от ориентации у; \lR \ = \kR \ = k0; k, l e R". Предельные значения hY(k, l) однозначно находятся из экспериментальных данных f (k, l) с помощью уравнений Фаддеева [5]:
hY(k, l) = f (k, l) + 2ni J hY(k, m)x
При комплексных волновых векторах классическое волновое поле переходит в обобщенное поле у(г, к), амплитуда рассеяния/(к, 1) (к, 1 е К ) переходит в обобщенную амплитуду рассеяния й(к, 1),
(к, 1 е С3).
х 0[(m - k, у)]8(m2 - k2)f (m, l)dm,
ni г
т.е. hY(k, l) -— I hY(k, m)х
k 0 J
I = kn
х 0[(m - k, у)] f (m, l)dm = f (k, l);
k, l e R3.
3
m e
m
Функция Хевисайда 6(5) = {0 при 5 < 0; 1 при 5 > 0} имеет в качестве аргумента скалярное произведение векторов т - к и у. Надо обратить внимание, что вектор у один и тот же в левой и правой частях уравнений Фаддеева (5). Поэтому подынтегральные значения Лу(к, т) зависят от взаимной ориентации как векторов к и у, так и векторов т и у. Тогда при численном решении значения Лу(к, т) находятся из системы уравнений, являющейся дискретным аналогом уравнений (5), в которой фиксируются векторы к и у, а вектор 1 принимает тот же набор дискретизованных значений, что и вектор т под интегралом.
Для удобства осуществления математических операций для комплексных к, 1 вводится функция Н(к, р) = Л(к, 1 = к - р), удовлетворяющая так называемому д -уравнению [6]:
= -2п | С,Я(к, Ч)х
Ze
(6)
х H(k + Z, p + Z)S(Z+2kZ)dZ, k e M,
Z2 + 2kZ = 0, т.е. (k + Z)2 =
(7)
v (X = -P) =
lim H(k, p),
Ik k e M
(8)
где Н(к, р) = Л(к, 1 = к - р); к, 1 е С3.
Для слабых, борновских, рассеивателей равенство (8) справедливо уже при |к71 —- 0. По мере увеличения силы рассеивателя становится необ-
ходимым учитывать эффекты перерассеяния внутри неоднородности. Именно для "ослабления" эффектов перерассеяния нужна мнимая добавка к/ Ф 0 к волновым векторам, причем чем сильнее рассеиватель, тем большие значения |к71 требуются для достижения асимптотического равенства V (-р) = Н(к, р). Соответственно увеличивается и модуль полного вектора |к |. С целью оценки минимальной величины |к | = |кт1П |, при которой асиптотика уже справедлива, использовалось уравнение типа Липпана-Швингера для обобщенного поля ц(г, к) = ехр(г'кг)у(г, к) и соответствующей ему обобщенной функции Грина 0(г, к)
[7]: ц(г, к) = 1 + ^ О (г - г', кМг')ц(г', к)й?г'.
Здесь ^ - область рассеяния; значение ц(г, к) = 1 отвечает падающему полю. С учетом выражения для пространственного спектра функции О(г, к), это уравнение принимает вид:
|(r, k) = 1 - -Ц f dr' f dX iO-n-V J J
(2 n)
r' e ^
X e[
exp [ i X ( r - r ' ) ]
x2 + 2 k X
где ]' = 1, 2, 3 - индекс декартовых компонент вектора; М - множество всех векторов к е С , которые при любом фиксированном векторе р е К одновременно удовлетворяют двум условиям (4). Черта над к ■ означает дифференцирование по комплексно сопряженной переменной. Функция 5^2 + 2kZ) выделяет те векторы Z, для которых
Следовательно, при фиксированном к е С3 действительный вектор Z принадлежит окружности радиуса |кд | с центром в точке -кд, и эта окружность лежит в плоскости, перпендикулярной вектору к7, так как из (7) следует
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.