Журнал вычислительной математики и математической физики, 2020, T. 60, № 6, стр. 1013-1026

1. ВВЕДЕНИЕ

Пусть скалярная функция $p({\mathbf{x}},t)$ определяет акустическое волновое поле, зависящее от координат ${\mathbf{x}} \in {{\mathbb{R}}^{3}}$ и времени $t \geqslant 0$. Поле создается в бесконечной среде источниками, локализованными в известной области $S$. Среда характеризуется локальной фазовой скоростью звука $c({\mathbf{x}})$ и имеет постоянную плотность. При этом известно, что $c({\mathbf{x}}) = {{c}_{0}} = {\text{const}}$ вне заданной области $X$, удовлетворяющей условию $X \cap S = \emptyset $. В самой области $X$ функция $c({\mathbf{x}})$ может быть переменной, и это интерпретируется как нахождение там акустических неоднородностей. В этом случае для гармонического источника вида $f({\mathbf{x}},\omega ){{e}^{{i\omega t}}}$ с известной частотой $\omega $ поле $p({\mathbf{x}},t)$ в рамках приближения линейной акустики может быть найдено в форме $p({\mathbf{x}},t) = u({\mathbf{x}},\omega ){{e}^{{i\omega t}}}$, где комплексная амплитуда $u({\mathbf{x}},\omega )$ удовлетворяет уравнению

и условию излучения (см., например, [1]). Здесь ${{k}_{0}} = \tfrac{\omega }{{{{c}_{0}}}}$, а $\xi ({\mathbf{x}}) = c_{0}^{{ - 2}} - {{c}^{{ - 2}}}({\mathbf{x}})$.

Нас интересует следующая обратная задача для уравнения (1): зная комплексную амплитуду поля $u({\mathbf{x}},\omega )$ в области $Y,Y \cap X = \emptyset $, $Y \cap S = \emptyset $, для некоторого набора частот $\omega $ найти коэффициент $\xi (x)$, т.е. функцию $c({\mathbf{x}})$, определяющую акустические неоднородности в области $X$. Вводя функцию Грина для уравнения Гельмгольца (1) в ${{\mathbb{R}}^{3}}$: $G(\rho ,\omega ) = - \tfrac{{exp\left( {i\omega \rho {\text{/}}{{c}_{0}}} \right)}}{{4\pi \rho }}$, можно при определенных предположениях о гладкости функций $u({\mathbf{x}},\omega )$, $f({\mathbf{x}},\omega )$, $c({\mathbf{x}})$ (см., например, [1]) свести обратную задачу к нелинейной относительно неизвестных $u({\mathbf{x}}{\text{'}},\omega )$, $\xi ({\mathbf{x}}{\text{'}})$, ${\mathbf{x}}{\text{'}} \in X$, системе интегральных уравнений:

Входящие в (2) функции – это известные (вычислимые) функции, а величины $u({\mathbf{x}}{\text{'}},\omega )$, $\xi ({\mathbf{x}}{\text{'}})$, ${\mathbf{x}}{\text{'}} \in X$, подлежат определению.

Задача (2) достаточно хорошо исследована теоретически: изучены вопросы существования и единственности ее решения (см., например, [2]–[4] и др.).

Во многих работах рассмотрены различные численные методы решения задачи (2) в двумерной и трехмерной постановках. Так, в [1] система (2) сводится к нелинейному операторному уравнению, которое затем решается специальным итерационным методом. При этом явно не учитывается, что это уравнение является некорректно поставленной задачей и итерационный метод не регуляризируется. Тем не менее метод работает для модельных “рассеивателей средней силы” [1, с. 101–105]. В работе [4] для решения системы уравнений (2) в трехмерном аксиально-симметричном случае использован регуляризованный метод Гаусса–Ньютона. В работе [5] для аналогичной задачи были применены специальный градиентный метод и метод Флетчера–Ривса. Другие градиентные методы были использованы в [6], [7].

Отметим, что имеются и альтернативные подходы в решении обратной задачи акустического зондирования, не связанные с системой (2). В частности, оригинальный метод граничного управления, позволяющий решать трехмерные задачи типа (2), был предложен и развит в работах [8], [9]. В [10] для решения двумерной обратной задачи рассеяния был применен известный метод Р.Г. Новикова [11], а в последующей работе [12] проведен сравнительный анализ разновидности этого метода и некоторых других функционально-аналитических методов решения двумерных обратных задач акустического рассеяния. Весьма многообещающими при обработке реальных экспериментальных данных оказались методы М.В. Клибанова, суммированные в монографии [13], а также методы из монографии [14]. Большой интерес представляют также недавние работы [15], [16].

Все упомянутые методы решения обратной задачи акустического зондирования требуют в трехмерной постановке значительных вычислительных ресурсов. В связи с этим отметим статьи [17], [18], в которых исходная обратная коэффициентная задача для волнового уравнения, из которой собственно и получается уравнение (1), сводится к трехмерному интегральному уравнению Фредгольма I рода с правой частью, содержащей специальные интегралы от регистрируемого поля. Предложенный метод решения этого уравнения оказался очень эффективным численно, и позволяет решать трехмерные обратные задачи для достаточно мелких сеток на персональном компьютере (ПК) даже без распараллеливания. Этот метод оказывается также очень полезным и для задачи (2), и это будет продемонстрировано ниже.

В данной статье мы придерживаемся следующей схемы решения нелинейной системы (2), сводящей последнюю к решению линейных интегральных уравнений:

1) решаем второе уравнение, записанное в виде

относительно функции ${v}({\mathbf{x}}{\text{'}},\omega ) = \xi ({\mathbf{x}}{\text{'}})u({\mathbf{x}}{\text{'}},\omega )$, ${\mathbf{x}}{\text{'}} \in X$;

2) вычисляем функцию $u({\mathbf{x}},\omega )$, ${\mathbf{x}} \in X$, из первого равенства системы (2), записанного в форме

3) находим функцию $\xi ({\mathbf{x}})$ из уравнения ${v}({\mathbf{x}},\omega ) = \xi ({\mathbf{x}})u({\mathbf{x}},\omega )$, ${\mathbf{x}} \in X$.

Близкая схема использовалась, например, в работе [19]. Однако там для решения некорректной задачи – трехмерного уравнения, аналогичного (3), не использовались методы регуляризации, а решение обратной задачи искалось на достаточно узком классе функций с “кусочно-постоянным током”. В итоге для решения такой трехмерной обратной задачи на сравнительно мелких сетках также требуются значительные вычислительные ресурсы.

Ниже мы покажем, что при некоторых, не очень обременительных, специальных предположениях о множестве регистрации данных $Y$ и о множестве $X$, где ищется решение обратной задачи, схема 1 )–3) эффективно реализуется численно и позволяет решать соответствующую трехмерную обратную задачу на достаточно мелких сетках даже на персональном компьютере средней производительности за несколько минут (без распараллеливания). Предлагаемый алгоритм решения обратной задачи и его численное исследование являются основными результатами этой работы.

2. СХЕМА РЕГИСТРАЦИИ ДАННЫХ ОБРАТНОЙ ЗАДАЧИ И РЕДУКЦИЯ ЗАДАЧИ

Будем рассматривать конкретную схему регистрации данных, т.е. функции $u({\mathbf{x}},\omega )$, “в плоском слое” $Y$ для декартовых координат ${\mathbf{x}} = (x,y,z)$. Эта схема была ранее использована в [17], [18] для решения обратной задачи с другими данными. На фиг. 1 схематично изображена соответствующая геометрическая схема: область неоднородностей поля скоростей $X$, область $S$, в которой расположены источники, и область регистрации рассеянного поля $Y$. Эти области имеют вид плоских слоев, по переменным $x$, $y$, перпендикулярных оси $Oz$. Звездочками показано возможное положение источников зондирующего поля. Таким образом,

Будем использовать двумерное преобразования Фурье ${{F}_{r}}[{\kern 1pt} \cdot {\kern 1pt} ]$ по переменным $r = (x,y)$, которое для функции $a(r,z,\omega ) = a(x,y,z,\omega )$ задается в виде

а также обратное преобразование $F_{\Omega }^{{ - 1}}[\tilde {A}(z,\omega ,\Omega )](r)$. Предполагая существование соответствующих образов Фурье $\tilde {U}$, ${{\tilde {U}}_{0}}$, $\tilde {G}$, $\tilde {V}$, $\tilde {W}$ функций $u$, ${{u}_{0}}$, $G$, ${v}$, $W$, формально преобразуем интегралы входящие в формулы (3), (4). Тогда из этих формул получим совокупность равенств: Эти равенства используются в дальнейшем для решения прямой и обратной задачи. Формулы (5) верны, если функции $u$, ${{u}_{0}}$, $G$, $W$ принадлежат классу медленно растущих функций [20], а функция $\xi $ финитна. Именно это будет предполагаться в дальнейшем. Особо отметим, что равенства (5) связывают величины $\tilde {U}(z,\omega ,\Omega )$, $\tilde {G}(z - z{\text{'}},\omega ,\Omega )$, $\tilde {V}(z{\text{'}},\omega ,\Omega )$ как функции аргументов $z$, $z{\text{'}}$, а числа $\omega $, $\Omega $ играют роль параметров.

3. АЛГОРИТМ РЕШЕНИЯ ПРЯМОЙ ЗАДАЧИ

Рассматриваемая нами прямая задача заключается в нахождении из системы (5) функции $\tilde {W}(z,\omega ,\Omega )$, $z \in [{{z}_{3}},{{z}_{4}}]$, по известной финитной функции $\xi (r,z{\text{'}})$ и заданной функции источников ${{\tilde {U}}_{0}}(z,a,\Omega )$. Для решения этой задачи используем следующий

Алгоритм 1

Шаг 1. Вычисляем преобразования Фурье  ${{\tilde {U}}_{0}}(z,\omega ,\Omega )$, $\tilde {G}(z - z{\text{'}},\omega ,\Omega )$ (алгоритмически – применяя быстрое дискретное преобразование Фурье (БПФ)).

Шаг 2. Для каждого из параметров $\omega $, $\Omega $ (матрично) реализуем итерационный процесс решения первых двух уравнений из (5):

начинающийся с функции ${{\tilde {U}}_{0}}(z,\omega ,\Omega )$.

Шаг 3. Останавливаем процесс по некоторому правилу на итерации номер $\nu $ и получаем приближенное решение ${{\tilde {U}}_{\nu }}(z,\omega ,\Omega )$.

Шаг 4. Последовательно вычисляем величины

и принимаем функцию ${{\tilde {W}}_{\nu }}(z,\omega ,\Omega )$ (или ее обратное преобразование Фурье) как приближенное решение прямой задачи.

Мы не будем здесь проводить теоретический анализ сходимости алгоритма 1, так как формально он не используется нами при решении обратной задачи. Он нужен лишь для генерации ее модельных данных. Отметим лишь, что, как следует из общей теории решения интегральных уравнений II рода (см., например, [21]), алгоритм будет быстро сходиться, по крайней мере, для “малых” $\omega $. Ниже будут продемонстрированы численные примеры, подтверждающие это положение.

4. ОБРАТНАЯ ЗАДАЧА

Она состоит в нахождении из системы (5) по заданной величине $\tilde {W}(z,\omega ,\Omega )$, $z \in [{{z}_{3}},{{z}_{4}}]$, функции $\xi (r,z{\text{'}})$, $z{\text{'}} \in [{{z}_{1}},{{z}_{2}}]$. Для этого мы применим

Алгоритм 2

Шаг 1. При заданной функции $\tilde {W}(z,\omega ,\Omega )$ решаем для всех рассматриваемых параметров $\omega $, $\Omega $ одномерные уравнения I рода

используя подходящий метод регуляризации (регуляризующий алгоритм, РА) этих некорректно поставленных задач.

Шаг 2. По найденной функции $\tilde {V}(z{\text{'}},\omega ,\Omega )$ вычисляем величину

и затем функции

Шаг 3. Находим $\xi (r,z{\text{'}})$ из уравнения $u(r,z{\text{'}},\omega )\xi (r,z{\text{'}}) = V(r,z{\text{'}},\omega )$. Это можно сделать, используя метод наименьших квадратов (МНК по величине $\omega $) или просто вычисляя решение в виде $\xi (r,z{\text{'}}) = V(r,z{\text{'}},\omega ){\text{/}}u(r,z{\text{'}},\omega )$. В последнем случае результат будет, вообще говоря, зависеть от $\omega $.

Сделаем некоторые замечания об алгоритме 2.

1. Для реализации ш.1 необходимо детализировать предположения о функциях $V$, $W$. Мы используем включения $V({\mathbf{x}},\omega ) \in {{L}_{2}}(\mathbb{R}_{{xyz}}^{3})$, $W({\mathbf{x}},\omega ) \in {{L}_{2}}(\mathbb{R}_{{xyz}}^{3})$ для каждой рассматриваемой частоты $\omega $. Первое из них вытекает из финитности функции $\xi $, а второе постулируется. В этом случае для уравнений (8) применимы многие методы решения линейных некорректно поставленных задач (например, [4], [5], [22], [23], [24] и др.).

2. Уравнение (3) (и порождаемые им уравнения (8)) могут иметь не единственное решение для используемого конечного набора частот $\omega $. Поэтому для реализации ш.1 алгоритма мы применяем РА, ориентированные на поиск нормальных решений (тихоновская регуляризация, метод TSVD). В случае единственности решения для (3) оно будет совпадать с вычисленным единственным нормальным решением [17], [18].

3. Уравнения (8) можно решать отдельно для каждой фиксированной частоты $\omega $, а можно решать их и как систему уравнений в совокупности для всех используемых частот. Соответственно в последнем случае ш.3 алгоритма реализуется по МНК.

5. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Везде далее полагается, что уравнения (5) записаны в безразмерном виде с ${{c}_{0}} = 1$, так что ${{k}_{0}} = \omega $. Мы будем также считать, что модельные источники задаются в форме

где ${{{\mathbf{x}}}_{m}}$ – координаты точечных $\delta $-образных источников. Тогда и фурье-образ этой функции ${{\tilde {U}}_{0}}(z,\omega ,\Omega )$ можно выписать аналитически.

Для численного исследования предлагаемых алгоритмов рассмотрены прямые и обратные модельные задачи с функциями $\xi (x,y,z)$ вида

где Эти функции моделируют небольшие локальные неоднородности среды, положение которых и соответствующие им распределения скоростей нужно найти. Именно на поиск таких неоднородностей “настраивается” алгоритм 2. Геометрически эти неоднородности имеют форму эллипсоидов.

Величина ${{A}_{0}}$ определяет “контраст”

искомого решения. В расчетах использовалось значение ${{A}_{0}} = 0.3$, которое соответствует контрасту 25.2. По классификации из [1, с. 33] такой рассеиватель можно считать “сильным”, учитывая его характерные размеры $l \sim 0.4$ и значения ${{c}_{0}} = 1$, $\omega = 2$: $\tfrac{{\Delta c}}{{{{c}_{0}}}} \gg \tfrac{{{{c}_{0}}}}{{l\omega }}$. Подобные рассеиватели достаточно часто встречаются на практике.

Исходные уравнения (3), (4) аппроксимировались конечно-разностным методом на равномерных сетках в областях

Размерности сеток в этих областях выбирались как $N \times N \times M$ и $N \times N \times {{M}_{1}}$ с числами $N$, $M$, ${{M}_{1}}$, которые будут указаны ниже. Число $\varepsilon > 0$, характеризующее “толщину слоя”, в котором измеряются данные для решения обратной задачи, выбиралось как $\varepsilon = 0.5$ (толстый слой) или $\varepsilon = 0.02$ (тонкий слой). Дискретные аналоги функций ${{\tilde {U}}_{0}}(z,\omega ,\Omega )$, $\tilde {G}(z - z{\text{'}},\omega ,\Omega )$ вычислялись для задаваемых частот $\omega $ по известным величинам ${{u}_{0}}(x,\omega )$ и $G(\rho ,\omega )$ с помощью быстрого преобразования Фурье. Затем они использовались в алгоритмах 1, 2.

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

при ${{A}_{m}} = 1$ и с аналогичными источниками в точках результаты получились весьма близкими. Поэтому мы приводим ниже результаты численных экспериментов только для источников первого вида.

6. НЕКОТОРЫЕ СВОЙСТВА АЛГОРИТМА 2

Все вычисления проводились в системе МАТЛАБ на ПК с процессором Intel (R) Core (TM) i7-7700 CPU 3.60 GHz, ОЗУ 16Гб (без распараллеливания). Алгоритм 2 решения обратной задачи оказался достаточно “быстрым”. Приведем данные соответствующих численных экспериментов для решения обратной задачи на одной частоте $\omega = 2$. Величины $M$, ${{M}_{1}}$ определяют скорость решения одномерного интегрального уравнения (8) при фиксированном $\Omega $. Соответствующее время решения ${{t}_{0}}(M,{{M}_{1}})$ практически мало меняется при переходе от одного уравнения к другому. Это время управляется требуемой разрешающей способностью алгоритма по переменной $z$. Для полного времени решения обратной задачи на указанных сетках фактически справедлива оценка $t(N,M,{{M}_{1}}) \approx {{t}_{0}}(M,{{M}_{1}}){{N}^{2}}$, и число $N$ управляется здесь требуемой разрешающей способностью по осям $x$, $y$. Приведем на фиг. 10 эту зависимость $t(N) = t(N,M,{{M}_{1}})$, полученную для фиксированных $M = 51$, ${{M}_{1}} = 51$.

Отметим, что решаемая обратная задача весьма чувствительна к ошибкам входных данных. При ее решении с двойной точностью, внесение в правую часть уравнения (8) случайных ошибок с амплитудой порядка ${{10}^{{ - 5}}}$ приводит к чрезмерному искажению приближенного решения при использовании и метода TSVD, и метода регуляризации. Это связано с весьма “быстрым” убыванием сингулярных чисел матриц ${{A}^{{(m)}}}$ для СЛАУ, решаемых на ш. 1 алгоритма 2, и это является специфической особенностью решаемой обратной задачи. Аналогичное свойство обратной коэффициентной задачи для волнового уравнения отмечалось ранее в работах [12], [13]. Соответствующие оценки погрешности при различных априорных предположениях на точное решение можно найти в [2], [3].

7. ВЫВОДЫ

Из численных экспериментов, проведенных в работе, можно сделать следующие выводы.

1. Трехмерная обратная коэффициентная задача для уравнения колебаний, возникающая при моделировании акустического зондирования, может быть решена численно с помощью предлагаемого алгоритма для достаточно мелких сеток за время порядка нескольких минут на ПК (даже без распараллеливания). Для этого следует, например, использовать схему регистрации данных обратной задачи “в плоском слое”.

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

3. Предлагаемый алгоритм позволяет достаточно надежно определять положения небольших локальных неоднородностей акустической среды при данных с малыми ошибками. Данные могут задаваться даже в “тонком слое”.

4. Алгоритм 2 при отдельном его применении для разных частот дает различную точность приближенных решений. Поэтому точность совместного решения обратной задачи на нескольких частотах сразу (по МНК) может оказаться несколько хуже, чем на некоторых конкретных частотах из используемого набора.

Авторы благодарны О.Д. Румянцевой и А.С. Шурупу за полезные замечания.