Численный метод решения обратной задачи для волнового уравнения в среде с локальной неоднородностью Текст научной статьи по специальности «Математика»

УДК 519.632

С. Г. Головин^, Е. В. Захаров2

ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ ДЛЯ ВОЛНОВОГО УРАВНЕНИЯ В СРЕДЕ С ЛОКАЛЬНОЙ НЕОДНОРОДНОСТЬЮ*

Для волнового уравнения рассматривается обратная задача определения границы локальной неоднородности по измерениям поля в ограниченной области расположения приемников в трехмерной среде. Задача сведена к системе интегральных уравнений. В работе предлагается итерационный метод решения обратной задачи и приводятся результаты вычислительных экспериментов.

Ключевые слова: волновое уравнение, уравнение Гельмгольца, неоднородная среда, интегральные уравнения, численный метод.

1. Введение. Граничная задача для волнового уравнения является классической моделью при изучении рассеяния акустических волн локальным включением. Подобные задачи получили широкое распространение в инженерной геофизике, акустике океана, сейсморазведке, дефектоскопии, физике атмосферы, медицине и многих других областях.

В работе исследуется среда, в которой распространение акустических колебаний описывается источниками с временной зависимостью exp(—iut) и рассеяние волн описывается уравнением Гельмгольца. Решение прямой задачи для уравнения Гельмгольца методом интегральных уравнений было предложено в [1-4]. В данной работе рассмотрена задача распространения акустических волн в трехмерной однородной среде, содержащей локальную неоднородность с гладкой поверхностью. Обратная задача состоит в определении поверхности, являющейся границей неоднородности по измерениям отраженного от неоднородности волнового поля, возбуждаемого точечными источниками. Такая постановка задачи распространена в задачах акустического зондирования.

2. Постановка задачи. Рассмотрим распространение акустических волн в трехмерной безграничной однородной среде fio, содержащей односвязную ограниченную область !1[ £ й3 с достаточно гладкой границей díl\ — локальную неоднородность. Однородная среда fio имеет постоянное значение скорости распространения волн со, а в области fii определена кусочно-гладкая функция с(М), где М = (x,y,z). Неоднородность облучается различными точечными источниками Mf = (xf,yf,Zf), расположенными в ограниченной области пространства fi/. Измерение поля, рассеянного неоднородностью, может проводиться лишь в пределах области приема fip. Границы всех областей не имеют общих точек.

Распространение акустических волн малой амплитуды можно описать с помощью волнового уравнения [5, 6]:

AF(M, t) - = F(t)S(M - Mf),

где А — оператор Лапласа, F(M, t) — поле в среде, зависящее от пространственной переменной М(ж, у, z) € i?3, t ^ 0, с(М) — скорость распространения волны в среде, S(-) — дельта-функция Дирака, F(t)S(M — Mf) — функция источника.

В предположении, что распространение акустических волн является установившимся (с временной зависимостью exp(^¿o;í)), поставим следующую граничную задачу для уравнения Гельмгольца:

A v{M) + h¡lv{M) = f{u)S{M-Mf), Me fi0,

Av(M) + kjv(M) = 0, Mefi i, ^

1 Факультет BMK МГУ, ст. преп., к.ф.-м.н., e-mail: sgolovina-msuQmail.ru

2 Факультет BMK МГУ, проф., д.ф.-м.н., e-mail: zspecQcs.msu.ru

* Работа выполнена при финансовой поддержке РФФИ, проект № 17-01-00525.

с условиями сопряжения на границе 1

'ду(М)

НМ)}д о, =0

дп

= О

ап.

и условиями излучения Зоммерфельда на бесконечности (при г = у/./-1' + //- + г1' оо)

= „ = о а

аг \г / \г

где у(М) = у(М,ш) — поле давления, [-]ап! — скачок значения функции на границе раздела сред,

п — внешняя нормаль к границе, ко = —, кЛМ) = ——-. Функцию источника можно аппрок-

с0 с(М)

симировать дельта-функцией /(ш)6(М — М/), где /(ш) — амплитуда поля, Mf — координаты источников поля, ш — частота. Эта аппроксимация основывается на предположении о малости линейных размеров источников в сравнении с масштабом решаемой задачи.

Введем новую функцию с(М) = с^2 — с~2(М), тогда с(М) = со^— с(М)сд и запишем

уравнение (1) в виде

АУ(М) + к$у(М) = ¡(ш)8(М - Мг) + ш2у(М)с(М). (2)

В этом уравнении правая часть состоит из источников первичного поля и источников вторичного поля ш2'о{М)с{М). Полное поле в уравнении (2) можно представить как сумму первичного и вторичного полей: и(М) = г>о(М) + у\(М). Функции г>о(М) и (М) удовлетворяют уравнениям Гельмгольца с разными правыми частями:

АУО(М) + к?ММ) = ¡(ш)5(М - МД

(3)

Д«1 (М) + ф! (М) = ш2у(М)с(М).

Заметим, что первичное и вторичное поля удовлетворяют условиям излучения Зоммерфельда на бесконечности. Используя функцию Грина

0(М, Р) = —I-ехр \i-RMP

4тгЯМр \ с0

где i?(M, Р) = у/{х — хо)2 + (у — уо)2 + (г — го)2 — расстояние между точками М(х,у,г) и Р(х, у, г), от дифференциального уравнения (2) перейдем к интегральному уравнению для поля:

у{М) = «о(М) + ш2 J 0(М, Р)с(Р)у(Р) йР, (4)

Пг

где г>о(М) = ] £?(М, Р)/(ш)5(Р — М/) йР имеет физический смысл волны, излучаемой источником

П;

в однородной среде.

Запишем уравнение (4) для случаев, когда М принадлежит локальной неоднородности П1 и области расположения приемников Пр (см. рис. 1):

у{М) = «о(М) + ш2 J 0(М, Р)С(Р)У(Р) йР, М е Пь

П1

(5)

V!(М) = иг I 0(М, Р)с(Р)у(Р)ёР, М е пр.

Интегрирование проводится по локальной области П1 С Д3, где функция с(М) отлична от нуля. Система (5) является нелинейной относительно функций с(М) и и(М) при М ё 01.

При решении прямой задачи, когда граница неоднородности П1 известна и соответственно известна функция с(М), М е П1, необходимо решить интегральное уравнение Фредгольма второго рода и найти функцию г>(М), М € О1, далее вычислить значение функции (М), М € Ор,

12 ВМУ, вычислительная математика и кибернетика, № 4

Рис. 1. Расположение исследуемых областей: 1 область расположения источников; 2 область расположения приемников; 3 область с локальным включением

используя систему (5) (при фиксированном источнике). Заметим, что при решении обратной задачи использование в вычислениях нескольких частот ш или нескольких источников существенно для единственности решения обратной задачи. В работе [7] приводится пример, когда решение системы (5) не единственно, если измерения рассеянного поля проводятся на одной частоте и при одном источнике.

Для проверки предложенной модели распространения акустических волн в однородной среде. содержащей сферическую локальную неоднородность, в работе [8] были проведены модельные расчеты. Наличие сферической симметрии позволяет применить для решения прямой задачи метод разделения переменных и использовать представление искомого поля в виде разложения в ряд по собственным функциям волнового оператора, построить решение в аналитической форме и сравнить полученные результаты с модельными расчетами. Выполненные расчеты показали достаточно хорошее совпадение решений прямой задачи аналитическим методом и методом интегральных уравнений, погрешность отклонения составила не более 2%.

Разработка численных методов решения обратной задачи, когда неизвестными являются функции с(М) и г>(М), М € Оь встречает ряд трудностей: данные и решение лежат в разных областях, различные размерности множества входных данных и решения, не выполнены обычные условия регулярности [9, 10].

В работе [11] рассмотрен абстрактный аналог системы (5) в виде операторного уравнения -Р(г) = 0, где оператор Р действует из гильбертова пространства И\ в другое гильбертово

пространство Н2, г = (у(М), с(М)) вектор неизвестных. Для решения использован метод Ньютона Гаусса [12]:

2„+1 = ащ'ттЦ^,,) + Р'п{% - гн)||2 = гп - (Р'*К',)-1

г£Я

где Р'п = Р'(гп). В итеративно-регуляризованном варианте метода Ньютона Гаусса на каждом шаге минимизируется по г функционал

Ф{а.п,гп,г) = + - гн)||2 + ап\\г - £„||2,

где ап параметр регуляризации, £„ некоторый элемент Н\. Тогда итерационный процесс можно записать в следующем виде:

• 1 — ( ¡'Ц РЦ "Ь " /< О ( РЦ Р„МСп), ) + М„( ч/( )) ■

При решении поставленной задачи предложенным методом возникают значительные трудности, связанные с размерностью массивов данных.

3. Построение итерационного процесса и результаты численных расчетов. Рассмотрим оптимизированный метод решения системы (5) с использованием метода Ньютона Гаусса. Запишем систему (5) в операторном виде:

Кг(у,с) = г'о(М) + ш2 [ С(М, Р)с(Р)у(Р) ЛР - г>(М) = О, М е Оь

П1

К2(у,с) = со2 I 0(М, Р)с(Р)у(Р) ЛР - гч(М) = О,

п,

М € О

•рч

где V = г>(М), с = с(М).

Для решения этой задачи будем использовать следующий итерационный процесс:

1) зададим начальное приближение с0 = 0 для функции с;

2) решив уравнение К1 (г'°, с0) = 0 стандартным методом, найдем функцию г>°;

3) решив уравнение Е^Ас1) = 0 регуляризованным методом Ньютона Гаусса, найдем с1;

4) решив уравнение А'^Дс1) = 0, найдем функцию V1 и т.д.

На каждом шаге предложенного итерационного процесса приходится численно решать интегральное уравнение, что замедляет скорость вычисления итераций, но позволяет на порядок снизить размерность обрабатываемых массивов. Отметим, что первый шаг процесса совпадает с решением обратной задачи в борновском приближении, когда отраженное от неоднородности поле заменяется первичным. Нахождение функции с(Р) сводится к решению уравнения Фредгольма первого рода:

гч (М) =ш2 I С(М, Р)с(Р)г>0(Р) ЛР.

Приведем результаты модельных расчетов, где в качестве входных данных использовалось решение прямой задачи с внесенной погрешностью 5%. Два источника располагались в области Q,f в точках с координатами (0, 0, 271) и (200, 200, 271), измерения проводились на сетке приемников 25 х 25, расположенных в области

Ор = {г = 275, 0 < ж < 800, 0 < у < 800}.

Исследуемая область О0 имела форму куба размера 200 х 200 х 200:

О0 = {70 < г < 270, 0 < я < 200, 0 < у < 200},

в ней имелась неоднородность О1 размера 90 х 20 х 90. Параметр регуляризации и частота выбирались следующим образом: а.п = 0.01, ш = 400.

На рис. 2 изображено точное решение обратной задачи (заданная неоднородность) и начальное приближение для итерационного процесса, в качестве которого бралась нулевая функция.

Рис. 2. Исследуемая неоднородность: а точное: решение; б начальное: приближение в итерационном

процессе

На рис. 3, 4 приводятся результаты вычислений предложенным итерационным методом. 13 ВМУ, вычислительная математика и кибернетика, № 4

Рис. 3. Результаты работы итерационного процесса: а поело первой итерации; б поело второй итерации; в поело третьей итерации; я поело четвертой итерации

в

Рис. 4. Результаты работы итерационного процесса: а, поело пятой итерации; б после шестой итерации; в поелс: седьмой итерации; я поелс: десятой итерации

Рисунки показывают, что уже начиная с пятой итерации форма неоднородности восстанавливается с хорошей точностью. Для получения искомого решения с погрешностью не более 2%

достаточно было выполнить 10 итераций. Результаты численного моделирования показали, что для одночастотного случая количество приемников должно быть достаточно большое и они должны располагаться в области большей, чем характерные размеры искомых границ неоднородностей.

СПИСОК ЛИТЕРАТУРЫ

1. Kolton D., Kress R. Inverse Acoustic and Electromagnetic Scattering Theory. Third edition. Vol. 93. Springer-Verlag, 2013.

2. Колтон Д., Кресс Р. Метод интегральных уравнений в теории рассеяния. М.: Мир, 1987.

3.Дмитриев В. И., Захаров Е. В. Метод интегральных уравнений в вычислительной электродинамике. М.: МАКС Пресс, 2008.

4. Тихонов А.Н., Самарский A.A. Уравнения математической физики. М.: Наука, 1972.

5. Горюнов A.A., Сасковец A.B. Обратные задачи рассеяния в акустике. М.: Изд-во МГУ, 1989.

6. Жданов М. С. Теория обратных задач и регуляризации в геофизике. М.: Научный мир, 2007.

7. Davaney A. J. Nonuniqueness in the inverse scattering problem //J. Math. Phys. 1978. 19. N 7. P. 15261531.

8. Головина С. Г., Никитина Е. В. Численный анализ методов определения спектральной амплитуды акустического поля // Вестн. Моск. ун-та. Сер. 15. Вычиел. матем. и киберн. 1997. № 3. С. 20-23. (Golovina S.G., Nikitina Е. Numerical analysis of methods for determining the spectrial amplitude of the acoustic field // Moscow Univ. Comput. Math, and Cybern. 1997. 15. N 3. P. 20-23.)

9. Тихонов A.H., Ареенин В.Я. Методы решения некорректных задач. М.: Наука, 1986.

10. Рамм А. Г. Многомерные обратные задачи рассеяния. М.: Мир, 1994.

11. Бакушинский А.Б., Левитан С.Ю. Некоторые модели и численные методы нелинейной вычислительной диагностики // Сборник трудов ВНИИСИ АН СССР. Вып. 13. М.: Изд-во ВНИИСИ, 1991. С. 3-25.

12. Бакушинский А.Б., Гончарский A.B. Итеративные методы решения некорректных задач. М.: Наука, 1989.

Поступила в редакцию 05.06.17