Моделирование решения акустической обратной задачи рассеяния для трехмерной нестационарной среды

1. ВВЕДЕНИЕ

Будем предполагать, что на множестве ℝ³×[0, T] определена функция c(x, t) > 0, зависящая от координат x = (x, y, z) и времени t$\in$[0, T]. Эта функция моделирует фазовую скорость звука в пространстве ℝ³, которое считается заполненным некоторой средой с постоянной плотностью. Более того, мы будем считать, что c(x, t) = c₀ = const > 0 вне заданной области X$\subset$ℝ³ для всех t$\in$[0, T]. Таким образом, под термином “нестационарная среда” мы будем иметь в виду среду, в которой имеются нестационарные акустические неоднородности только в области X, т.е. скорость звука в ней зависит от пространственной переменной и времени. Далее, мы считаем, что в ℝ³ при t$\in$[0, T] определено скалярное акустическое волновое поле p(x, t), создаваемое источниками F₀(x, t), которые локализованы в известной ограниченной области S, X∩S = $\varnothing$. Это поле регистрируется в области Y (X$\subset$ℝ³, Y∩X = $\varnothing$, Y∩S = $\varnothing$). Рассматриваемая нами обратная задача в общей форме заключается в нахождении неизвестной функции c(x, t) для x$\in$ X, t$\in$[0, T] по данным p(x, t), x$\in$Y, t$\in$[0, T].

Для задач такого типа часто используется приближение линейной акустики. Тогда математической моделью, связывающей p(x, t) и c(x, t), будет система равенств следующего вида:

$\left\{\begin{array}{l}\begin{array}{l}{\Delta }_{x}p\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-\frac{1}{c^2\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}\frac{{\partial }^{2}p\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial t^2}\mathrm{<mo>=</mo>}\\ =F_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>,}x\in {ℝ}^3\mathrm{,}t\in \mathrm{<mo\; stretchy="false">(</mo>0,}T\mathrm{<mo\; stretchy="false">]</mo>}\end{array}\hfill \\ p\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,0<mo\; stretchy="false">)</mo><mo>=</mo>}\phi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>,}\frac{\partial p\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,0<mo\; stretchy="false">)</mo>}}{\partial t}\mathrm{<mo>=</mo>}\psi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>,}x\in {ℝ}^3\mathrm{,}\hfill \end{array}\right.$ (1)

т.е. прямой задачей (задачей нахождения p(x, t) по заданным c(x, t), $\phi$(x), ψ(x), F₀(x, t)) в данном случае будет задача Коши (1), а обратной задачей – обратная коэффициентная задача для волнового уравнения из (1): найти c(x, t) по известным $\phi$(x), ψ(x), F₀(x, t) и дополнительным данным p(x, t), x$\in$Y, t$\in$[0, T]. Функциональные пространства, в которых будут решаться эти задачи, а также требования на величины c(x, t), $\phi$(x), ψ(x), F₀(x, t) будут указаны ниже.

Модель (1) хорошо известна, и прямая задача для нее подробно изучена теоретически в различных аспектах и в различных функциональных пространствах (см., например, [1, 2] и др.). Обратная задача также рассматривалась многими авторами. Однако в этом случае на зависимость функции c(x, t) от переменных x и t накладывались серьезные ограничения. Так, в немногих работах считалось, что c(x, t) = c(t), и для таких обратных задач устанавливались теоремы существования и единственности решения (см., например, [3]). В большинстве работ по этой тематике предполагается, что c(x, t) = c(x). Обратные задачи с такой функцией c(x) исследованы весьма подробно. В частности, выделены многие функциональные классы, на которых задача имеет единственное решение (см., например, [3–7]). Предложены и обоснованы численные алгоритмы решения обратных задач такого рода с различными вариантами данных ([8–15] и др.).

Уравнение (1) для случая c(x, t) = c(x) рассматривалось также и в спектральной форме. Тогда для гармонического источника вида f(x, ω)e^iωt с известной частотой ω поле p(x, t) может быть найдено в виде p(x, t) = u(x, ω)e^iωt, где комплексная амплитуда u(x, ω) удовлетворяет уравнению

$\begin{array}{c}\Delta u(x,\omega )+k_0^{2}u(x,\omega )=\\ =f(x,\omega )+{\omega }^2\xi (x)u(x,\omega ),x\in {ℝ}^3.\end{array}$(2)

Здесь $k_0=\frac{\omega }{c_0},$ $\xi (x)=c_0^{-2}-c^{-2}(x)$. Нахождение такой функции u(x, ω) по заданному коэффициенту c(x) и заданному f(x, ω) составляет прямую задачу для (2). Обратная задача для этого уравнения обычно ставится так: зная для некоторого набора частот ω комплексную амплитуду поля u(x, ω) в области Y, найти коэффициент ξ (x), а значит, и функцию c(x), определяющую акустические неоднородности в области X.

Обратная задача для (2) также достаточно хорошо исследована теоретически: изучены в определенных аспектах вопросы существования и единственности ее решения для различных видов областей X, Y, вопросы ее устойчивости по отношению к возмущениям данных и т.д. (см. [7, 12, 15] и др.) Например, известно, что для единственности решения такой обратной задачи в общем случае необходимо использовать по меньшей мере счетное число частот. Однако на практике реализовать это невозможно, и приходится ограничиваться сравнительно небольшим их числом. Анализ этой проблемы дан, например, в [14]. Тем не менее, такое изучение акустических рассеивателей на сравнительно небольшом числе частот часто встречается на практике. Для получения в этом случае определенного решения обратной задачи необходимо использовать дополнительные априорные ограничения, сужающие множество возможных решений.

Во многих работах рассмотрены разнообразные численные методы решения обратных задач для (2) в двумерной и трехмерной постановках. В связи с этим отметим работы [13, 15–22]. Обзор этих и подобных подходов можно найти, например, в [14]. Отметим, что все эти методы требуют значительного времени вычислений при решении трехмерной обратной задачи.

В этой работе мы рассматриваем обратную задачу для (1) в предположении, что области X и Y имеют вид плоских слоев, в случае специального источника F₀ и для функции c = c(x, t) общего вида, о которой ниже будут сделаны некоторые дополнительные предположения. Наша цель – предложить и апробировать численный алгоритм, позволяющий решать такую обратную задачу на персональном компьютере (ПК) за минуты. Рассматриваемая обратная задача является, вообще говоря, некорректно поставленной. Она может не иметь решений в рассматриваемом функциональном классе, решение в случае его существования может быть не единственным и неустойчивым по отношению к возмущениям данных случайной ошибкой. Поэтому решить нашу обратную задачу без дополнительных предположений о решении нельзя. Такими предположениями могут быть условия гипотетической теоремы существования и единственности решения обратной задачи (если такие теоремы для случая c = c(x, t) имеются) или другие допущения. Мы не занимаемся теоретическим исследованием свойств существования и единственности решений рассматриваемой обратной задачи. Более того, из-за отсутствия значимых теоретических результатов в этом направлении для обратной задачи c = c(x, t) мы будем предполагать существование ее решения в рассматриваемом ниже классе функций. Из-за возможной неединственности решения (квазирешения) мы будем пытаться выделить ее решение с некоторыми оптимальными свойствами. Нас устраивает и альтернатива – предположение о том, что для имеющихся данных обратная задача имеет единственное решение в рассматриваемом функциональном классе. Рассмотренные нами модельные обратные задачи подобного рода показывают, что такие допущения вполне возможны для небольших по размеру локальных нестационарных неоднородностей среды. В рамках этих допущений предлагаемый нами алгоритм дает достаточно точные приближенные решения при “малых” погрешностях данных. Эти приближения численно устойчивы в том смысле, что при стремлении погрешностей данных к нулю приближения сходятся к решению, полученному для точных данных.

2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ЗОНДИРОВАНИЯ СРЕДЫ

Примем следующую упрощенную модель зондирования среды с переменной во времени скоростью звука. Будем считать, что зондирование производится в моменты времени $t_k\mathrm{<mo>=</mo>}k\Delta t\mathrm{,<mi>\; </mi>}k\mathrm{<mo>=</mo>0,1,...,}n\mathrm{,<mi>\; </mi>}\Delta t\mathrm{<mo>=</mo>}\frac{T}{n}$ импульсным источником вида

$F_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^n}g\mathrm{<mo\; stretchy="false">(</mo>}t-t_k\mathrm{<mo\; stretchy="false">)</mo>.}$

Здесь финитная функция f(x) с носителем S описывает пространственную локализацию источника, а функция g(t) определяет его временной профиль. Предполагается, что g(t) финитна и “быстро убывает” на своем носителе [0, h_t], 2h_t << Δ_t. Таким образом, во временной области источник представляет собой последовательность коротких импульсов, следующих друг за другом с интервалами Δt – h_t (Δt – h_t >> h_t). На этих временных интервалах происходит сбор и начало обработки информации. Пример характерной зависимости F₀(x, t) при фиксированном x приведен на рис. 1.

Рис. 1. Пример функции источника F₀(x, t) при фиксированном x.

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

Предположения 1. а) На каждом временном отрезке [t_k–1, t_k] пространственное распределение скоростей звука в среде практически не меняется, т.е. можно считать, что c_k(x, t) = c_k(x) для t$\in$[t_k–1, t_k], где c_k(x) – некоторое усредненное на рассматриваемом временном отрезке распределение скорости. б) На каждом временном отрезке [t_k–1, t_k] за время Δt – h_t после импульса вызванные им звуковые колебания среды в областях X и Y практически затухают, и среда в этих областях возвращается в исходное состояние. С точки зрения математической модели эти предположения означают, что на каждом интервале зондирования (t_k–1, t_k), k = 1, 2, …, n можно считать выполненными аналоги соотношений (1):

$\left\{\begin{array}{l}\begin{array}{l}{\Delta }_{x}p(x,t)-\frac{1}{c_k^2(x)}\frac{{\partial }^{2}p(x,t)}{\partial t^2}=\\ =f(x)g(t-t_{k-1}),x\in {ℝ}^3,t\in (t_{k-1},t_k],\end{array}\hfill \\ p(x,t_{k-1})=\mathrm{0,}\frac{\partial p(x,t_{k-1})}{\partial t}=\mathrm{0,}x\in {ℝ}^3.\hfill \end{array}\right.$(3)

Относительно коэффициентов этой задачи Коши сделаем следующие

Предположения 2. а) c_k(x)$\in$C(ℝ³); б) 0 < c_min ≤ c_k(x) ≤ c_max, $\forall$x$\in$ℝ³, где c_min, c_max – известные константы; c_k(x) = c₀ при x$\in$ℝ³\X; в) f(x)$\in$C(S), g(t) равна нулю вне отрезка [0, Δt], g(t)$\in$C[0, Δt], причем |g(t)| – быстро убывающая на [0, Δt] функция; в качестве g(t) на [0, Δt] можно, например, взять g(t) = e^–μt, μ$\in$C с достаточно большим значением Reμ > 0.

Известно, что задача (3) имеет единственное классическое решение: p_k(x, t)$\in$C²(ℝ³×[t_k–1, t_k]) (см., например, [23]).

Обсудим смысл Предположения 1б. Формально в случае постоянной скорости c_k(x) = c₀ и точечного источника, т.е. при f(x) = δ(x – x₀), x₀$\in$S, эта задача имеет решение

$p\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{g\left(t-\frac{r}{c_0}\right)}{4\pi c_0^{2}r}\mathrm{,}r\mathrm{<mo>=</mo><mo>|</mo>}x-x_0\mathrm{<mo>|</mo>,}$

(см., например, [24, с. 514]), и поэтому из-за Предположения 2в о виде функции g(t) оно практически равно нулю в каждой фиксированной точке области X и Y при «достаточно большом времени» t. Таким образом, в этом случае выполнено Предположение 1б, если взять достаточно большой интервал времени (t_k–1, t_k), т.е. “большое” Δt. Само же Предположение 1б можно рассматривать как неявное ограничение на класс функций c_k(x) и функцию источника g(t).

Пусть p_k(x, t), k = 0, 1 …, n, – решение задач (3). Из Предположений 1 и 2 следует, что существуют преобразования Фурье F[p_k(x, t)](ω) этих функций по переменной t, и можно считать выполненными формулы:

$\begin{array}{c}{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}F\mathrm{<mo\; stretchy="false">[</mo>}{p}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo><mo\; stretchy="false">(</mo>}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\\ ={\displaystyle {\int }_{t_{k-1}}^{t_k}{p}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}{e}^{-i\omega \mathrm{<mo\; stretchy="false">(</mo>}t-t_{k-1})}dt}\mathrm{,}\\ u_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}\in C_x^2\mathrm{<mo\; stretchy="false">(</mo>}{ℝ}^3\mathrm{<mo\; stretchy="false">)</mo>.}\end{array}$

Из (3) тогда получим

$\begin{array}{c}\Delta u_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}+k_0^2{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\\ =f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\tilde{g}\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{<mo\; stretchy="false">)</mo>}+{\omega }^2{\xi }_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>,}x\in {ℝ}^3\mathrm{,}\end{array}$ (4)

в котором $k_0\mathrm{<mo>=</mo>}\frac{\omega }{c_0}$, g̃(ω) = F[g](ω) и ${\xi }_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{c}_0^{-2}-c_k^{-2}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>.}$ ${\xi }_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{c}_0^{-2}-c_k^{-2}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>.}$

Обозначим множество функций c_k(x), подчиненных условиям из Предположений 1 и 2, как C.

Пусть теперь нам известны величины u_k(x, ω), x$\in$Y. Тогда можно для каждого k поставить обратную коэффициентную задачу для уравнения (4): по данным u_k(x, ω), x$\in$Y найти функцию ξ_k(x). Эту обратную задачу мы будем называть вспомогательной обратной задачей. Зная ее решение ξ_k(x) для каждого k, можно найти c_k(x)$\in$C и построить решение обратной задачи зондирования в форме

$c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^n}{c}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\chi \mathrm{<mo\; stretchy="false">(</mo>}{t}_{k-1}\mathrm{,}{t}_k)\mathrm{.}$ (5)

Здесь χ(∙) – характеристическая функция множества. Итак, для решения обратной задачи зондирования нестационарной среды (в нашей постановке) при t$\in$[0, T] необходимо решать вспомогательную обратную задачу на каждом временном интервале работы источника. Для такой многократной процедуры требуется достаточно “быстрый” алгоритм решения вспомогательной задачи.

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

3. РЕДУКЦИЯ ВСПОМОГАТЕЛЬНОЙ ОБРАТНОЙ ЗАДАЧИ К ОДНОМЕРНЫМ ИНТЕГРАЛЬНЫМ УРАВНЕНИЯМ

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

$\begin{array}{c}G\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}=G\mathrm{<mo\; stretchy="false">(</mo><mo>|</mo>}x-x\mathrm{<mo>\text{'}</mo><mo>|</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}=\\ =-\frac{e^{-i{k}_{0}R}}{4\pi R}\mathrm{,}R=|x-x\mathrm{<mo>\text{'}</mo><mo>|</mo>}\end{array}$

для уравнения Гельмгольца (4), можно свести обратную коэффициентную задачу для (4) к нелинейной относительно неизвестных u_k(x', ω), ξ_k(x'), x'$\in$X системе интегральных уравнений:

$\begin{array}{c}{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}+\\ +{\omega }^2{\displaystyle {\int }_X}G\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}{\xi }_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo><mo\; stretchy="false">)</mo>}{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}dx\mathrm{<mo>\text{'}</mo>,}x\in X\mathrm{,}\\ {\omega }^2{\displaystyle {\int }_X}G\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}{\xi }_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo><mo\; stretchy="false">)</mo>}{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}dx\mathrm{<mo>\text{'}</mo>}=\\ \mathrm{<mo>=</mo>}{w}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>,}x\in Y,\end{array}$ (6)

(см., например, [7, 12–14] и др.). Здесь u_k(x, ω)$\in$C² (ℝ³), ξ_k(x)$\in$C(ℝ³). Входящие в (6) функции

$\begin{array}{c}{u}_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\tilde{g}\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle {\int }_S}G\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo><mo\; stretchy="false">)</mo>}dx\mathrm{<mo>\text{'}</mo>,}\\ w_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}-u_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>,}x\in Y,\end{array}$

– это известные (вычислимые) функции, такие, что u₀(x, ω)$\in$C² (ℝ³), w_k(x, ω)$\in$C² (ℝ³).

Во всех этих уравнениях частота ω играет роль параметра. Поэтому мы придерживаемся следующей схемы решения нелинейной системы (6) как обратной задачи для каждой используемой частоты ω:

1) решаем второе уравнение (линейное интегральное уравнение Фредгольма первого рода), записанное в виде

${\omega }^2{\displaystyle {\int }_X}G\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}{v}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}dx\mathrm{<mo>\text{'}</mo>}=w_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>,}x\in Y$ (7)

относительно функции

$v_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\xi }_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo><mo\; stretchy="false">)</mo>}{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>,}x\mathrm{<mo>\text{'}</mo>}\in X$;

2) по найденной v_k(x', ω) вычисляем функцию u_k(x, ω), x$\in$X из первого равенства системы (6), записанного в форме

$\begin{array}{c}{u}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}+\\ +{\omega }^2{\displaystyle {\int }_X}G\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}{v}_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}dx\mathrm{<mo>\text{'}</mo>,}x\in X\mathrm{<mo>;</mo>}\end{array}$ (8)

3) находим функцию ξ_k(x) из уравнения v_k(x, ω) = ξ_k(x)u_k(x, ω), x$\in$X, используя найденные функции v_k(x, ω), u_k(x, ω).

Обратная задача (6) рассматривалась, например, в [11, 13–19], и там разработаны алгоритмы, позволяющие решать ее за значительное время на суперкомпьютере и другой сложной вычислительной технике. С другой стороны, в работах [25, 26] был предложен экономичный метод, позволяющий решать ее в реальном времени на (ПК) средней производительности. Этот метод и будет применен нами ниже. Он основан на сведении обратной задачи (6) к одномерным интегральным уравнениям в предположении, что области X, Y имеют специальный вид плоских слоев: $X\mathrm{<mo>=</mo>}{ℝ}_{xy}^2\times \mathrm{<mo\; stretchy="false">[</mo>}{z}_1\mathrm{,}{z}_2\mathrm{<mo\; stretchy="false">]</mo>,}Y\mathrm{<mo>=</mo>}{ℝ}_{xy}^2\times \mathrm{<mo\; stretchy="false">[</mo>}{z}_3\mathrm{,}{z}_4\mathrm{<mo\; stretchy="false">]</mo>}$. На рис. 2 представлена схема расположения этих областей. Там же условно показано и возможное положение источников.

Рис. 2. Геометрическая схема регистрации данных обратной задачи: X – область акустических неоднородностей, Y – область регистрации данных u_k(x, ω), звездочки – возможные положения источников поля.

Будем далее для краткости опускать индекс k в уравнениях (7), (8). Введем двумерное преобразование Фурье F[∙] по переменным (x, y), которое для произвольной функции $a\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}z\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}\in L_2({ℝ}_{xy}^2)$ задается в виде

$\begin{array}{c}\tilde{A}\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}F\left[a\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}z\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}\right]\left(\Omega \right)\mathrm{<mo>=</mo>}\\ ={\displaystyle {\int }_{{ℝ}_{xy}^2}}a\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}z\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}{e}^{i\mathrm{<mo\; stretchy="false">(</mo>}{\Omega }_{1}x+{\Omega }_{2}y)}dxdy\mathrm{,}\Omega \mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\Omega }_1\mathrm{,}{\Omega }_2\mathrm{<mo\; stretchy="false">)</mo>,}\end{array}$

а также соответствующее обратное преобразование

$F_{\Omega }^{-1}\left[\tilde{A}\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo>}\right]\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo>.}$

Используя преобразования Фурье Ũ, Ũ₀, G̃, Ṽ, W̃ функций u, u₀, G, v, w, выразим через них по теореме о свертке интегралы

$\begin{array}{c}{\displaystyle {\int }_X}G\mathrm{<mo\; stretchy="false">(</mo><mo>|</mo>}x-x\mathrm{<mo>\text{'}</mo><mo>|</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}dx\mathrm{<mo>\text{'}</mo>,}x\in X\mathrm{<mo>;</mo>}\\ {\displaystyle {\int }_X}G\mathrm{<mo\; stretchy="false">(</mo><mo>|</mo>}x-x\mathrm{<mo>\text{'}</mo><mo>|</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>\text{'}</mo>,}\omega \mathrm{<mo\; stretchy="false">)</mo>}dx\mathrm{<mo>\text{'}</mo>,}x\in Y\mathrm{,}\end{array}$

входящие в формулы (7), (8). Тогда из этих формул получим совокупность равенств:

$\begin{array}{c}{\omega }^2{\displaystyle {\int }_{z_1}^{z_2}\tilde{G}\mathrm{<mo\; stretchy="false">(</mo>}z-z\text{'}\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo>}\tilde{V}\mathrm{<mo\; stretchy="false">(</mo>}z\text{'}\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo>}dz\text{'}}\mathrm{<mo>=</mo>}\\ =\tilde{W}\left(z\mathrm{,}\omega \mathrm{,}\Omega \right)\mathrm{,}z\in \mathrm{<mo\; stretchy="false">[</mo>}{z}_3\mathrm{,}{z}_4\mathrm{<mo\; stretchy="false">]</mo>}\end{array},$ (9)

$\begin{array}{c}\tilde{U}\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\tilde{U}}_0\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo>}+\\ +{\omega }^2{\displaystyle {\int }_{z_1}^{z_2}\tilde{G}\mathrm{<mo\; stretchy="false">(</mo>}z-z\text{'}\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo>}\tilde{V}\mathrm{<mo\; stretchy="false">(</mo>}z\text{'}\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo>}dz\text{'}},z\in \mathrm{<mo\; stretchy="false">[</mo>}{z}_1\mathrm{,}{z}_2\mathrm{<mo\; stretchy="false">]</mo>,}\end{array}$ (10)

а также формулу

$\begin{array}{c}\tilde{V}\mathrm{<mo\; stretchy="false">(</mo>}z\text{'}\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\\ =F[\xi \mathrm{<mo\; stretchy="false">(</mo>}x,y\mathrm{,}z\text{'}\mathrm{<mo\; stretchy="false">)</mo>}{F}_{\Omega }^{-1}[\tilde{U}\mathrm{<mo\; stretchy="false">(</mo>}z\text{'}\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo>}]\mathrm{<mo\; stretchy="false">(</mo>}x,y\mathrm{<mo\; stretchy="false">)</mo>}](\Omega )\mathrm{,}\\ z\text{'}\in \mathrm{<mo\; stretchy="false">[</mo>}{z}_1\mathrm{,}{z}_2\mathrm{<mo\; stretchy="false">]</mo>.}\end{array}$ (11)

Отметим, что равенства (9)‒(11) связывают величины Ũ(z, ω, Ω), G̃(z – z', ω, Ω), Ṽ(z', ω, Ω), как функции аргументов z, z', а числа ω, Ω играют роль параметров.

Как уже отмечалось, уравнения (7) для каждого фиксированного ω (или для конечного набора таких частот) могут иметь не единственное решение. Однако каждое такое уравнение имеет единственное нормальное решение (см., например, [27]), т.е. то решение v(x, ω), которое из всех возможных, составляющих множество V, имеет минимальную норму:

${‖\overline{v}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}‖}_{L_2\mathrm{<mo\; stretchy="false">(</mo>}X\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\min \{{‖v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}‖}_{L_2\mathrm{<mo\; stretchy="false">(</mo>}X\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>:</mo><mi>\; </mi>}v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}\in \mathcal{V}\}\mathrm{.}$

Аналогичные положения верны и для одномерных уравнений (9): для каждого из них существует единственное нормальное решение V(z', ω, Ω) с минимальной нормой в пространстве L₂[z₁, z₂] из всех возможных решений. Можно установить связь нормальных решений уравнений (9) и нормального решения уравнения (7):

$\overline{v}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{F}_{\Omega }^{-1}\left[\overline{V}\mathrm{<mo\; stretchy="false">(</mo>}z\text{'}\mathrm{,}\omega \mathrm{,}\Omega \mathrm{<mo\; stretchy="false">)</mo>}\right]\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo>.}$

Для поиска нормальных решений линейных интегральных уравнений имеются разработанные и многократно апробированные численные методы, устойчивые по отношению к возмущениям данных (регуляризующие алгоритмы, РА) (см., например, [2, 3, 27, 28] и др.).

4. АЛГОРИТМ РЕШЕНИЯ ВСПОМОГАТЕЛЬНОЙ ОБРАТНОЙ ЗАДАЧИ

Равенства (9)‒(11) используются для решения вспомогательной обратной задачи (6), причем описанная выше схема ее решения 1)‒3) переходит в следующую процедуру.

Алгоритм. Для каждой рассматриваемой частоты ω выполняем следующие шаги.

Шаг 1. По известным величинам u₀, w, G вычисляем их преобразования Фурье Ũ₀(z, ω, Ω), W̃₀(z, ω, Ω), G̃(z – z', ω, Ω).

Шаг 2. Решаем одномерное линейное интегральное уравнение Фредгольма первого рода (9) относительно функции Ṽ₀(z, ω, Ω), z$\in$[z₁, z₂], для каждого Ω.

Шаг 3. Вычисляем функцию Ũ(z, ω, Ω), z$\in$[z₁, z₂], из равенства (10), используя найденную Ṽ(z, ω, Ω).

Шаг 4. С помощью обратного преобразования Фурье по переменной Ω находим по Ṽ(z, ω, Ω) и Ũ(z, ω, Ω) функции v(x, ω), u_k(x, ω).

Шаг 5. По ним вычисляем функцию ξ_k(x) из уравнения v(x, ω) = ξ_k(x)u_k(x, ω), x$\in$X.

Шаг 6. Вычисляем функцию c_k(x) из равенства ${\xi }_k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{c}_0^{-2}-c_k^{-2}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>.}$

Обсудим конкретную реализацию шагов алгоритма.

А) При выполнении шага 1 проводится дискретизация задачи. Для этого вводится параллелепипед П = [a_x, b_x]×[a_y, b_y]×[a_z, b_z], содержащий области X, Y, S. В нем области X и Y заменяются слоями П_X = [a_x, b_x]×[a_y, b_y]×[z₁, z₂], и П_Y = [a_x, b_x]×[a_y, b_y]×[z₃, z₄] соответственно. Заданные в П_Y функции u₀(x, ω), w (x, ω) аппроксимируются на равномерной трехмерной сетке в этой области трехмерными массивами чисел. Тогда двумерные преобразования Фурье по переменным x, y функций u₀(x, ω), w (x, ω) соответствуют вычислению дискретных преобразований Фурье по этим переменным в прямоугольнике [a_x, b_x]×[a_y, b_y], которое выполняется с помощью алгоритма БПФ (быстрого преобразования Фурье). В итоге получаются параметрические векторы ${\widehat{u}}_0\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}[{\widehat{U}}_0\mathrm{<mo\; stretchy="false">(</mo>}{z}_p\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>}]\mathrm{,}\widehat{w}\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}[\widehat{W}(z_l\mathrm{,}\omega \mathrm{,}\widehat{\Omega })],$ ${\widehat{u}}_0\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}[{\widehat{U}}_0\mathrm{<mo\; stretchy="false">(</mo>}{z}_p\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>}]\mathrm{,}\widehat{w}\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}[\widehat{W}(z_l\mathrm{,}\omega \mathrm{,}\widehat{\Omega })],$ где ${\mathrm{<mo>\{</mo>}{z}_p\mathrm{<mo>\}</mo>}}_{p\mathrm{<mo>=</mo>1}}^{N_1}\mathrm{,<mi>\; </mi><mo>\{</mo>}{z}_l{\mathrm{<mo>\}</mo>}}_{l\mathrm{<mo>=</mo>1}}^{N_2}$ – одномерные сетки в [z₁, z₂] и в [z₃, z₄] соответственно, а Ω есть переменная двумерного дискретного преобразования Фурье. Аналогичная процедура проводится и для функции G(|x – x'|, ω), так что в итоге получается параметрическая матрица {G̃(z_l – z_j, ω, Ω)}. В результате уравнения (9), (10) переходят в системы линейных уравнений (СЛАУ)

$\begin{array}{c}{\omega }^2{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^{N_1}}{\nu }_{lj}\tilde{G}\mathrm{<mo\; stretchy="false">(</mo>}{z}_l-z_j\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>}\widehat{V}\mathrm{<mo\; stretchy="false">(</mo>}{z}_j\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\\ =\widehat{W}\left(z_l\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\right)\mathrm{,}l\mathrm{<mo>=</mo>1,...,}{N}_2\mathrm{,}\end{array}$ (12)

$\begin{array}{c}\widehat{U}\mathrm{<mo\; stretchy="false">(</mo>}{z}_p\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\widehat{U}}_0\mathrm{<mo\; stretchy="false">(</mo>}{z}_p\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>}+\\ +{\omega }^2{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^{N_1}}{\nu }_{pj}\tilde{G}\mathrm{<mo\; stretchy="false">(</mo>}{z}_p-z_j\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>}\widehat{V}\mathrm{<mo\; stretchy="false">(</mo>}{z}_j\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>,}\\ p\mathrm{<mo>=</mo>1,...,}{N}_1\mathrm{,}\end{array}$ (13)

в которых конечномерными аналогами неизвестных Ṽ(z', ω, Ω), Ũ(z, ω, Ω), являются векторы ${[\widehat{V}\mathrm{<mo\; stretchy="false">(</mo>}{z}_j\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>}]}_{j\mathrm{<mo>=</mo>1}}^{N_1}\mathrm{,}{[\widehat{U}\mathrm{<mo\; stretchy="false">(</mo>}{z}_p\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>}]}_{p\mathrm{<mo>=</mo>1}}^{N_2}$, а ν_lj, ν_pj – квадратурные коэффициенты аппроксимации интегралов. Детали описанной процедуры дискретизации можно найти в [25, 26].

Б) Шаг 2 алгоритма связан с решением линейного интегрального уравнения Фредгольма первого рода (9) в пространстве L₂. Это уравнение и получающаяся в результате его дискретизации система линейных уравнений (12) наследуют типичное свойство исходного уравнения (7) – некорректную постановку этой обратной задачи. Уравнения (7), (9), (12) могут не иметь решений для используемых экспериментальных данных w просто потому, что математическая модель не может их адекватно описать. В случае разрешимости этих обратных задач, возможна неоднозначность решения и неустойчивость решений по отношению к малым вариациям данных. Из теории некорректно поставленных задач известно, что эти трудности преодолеваются с помощью применения специальных РА, позволяющих устойчиво найти приближенные нормальные решения обратных задач ([2, 3, 27, 28] и др.).

По этой причине мы применяем для реализации шага 2 регуляризующий алгоритм, ориентированный на поиск нормальных решений (тихоновская регуляризация, метод TSVD). В результате мы получаем приближенное решение обратной задачи, которое имеет минимальную норму “вторичных источников” v(x, y, z, ω) = u(x, y, z, ω)ξ(x, y, z) (см. [14]), и это позволяет выделить соответствующее единственное решение обратной задачи. В случае малых линейных размеров нестационарных акустических неоднородностей такое нормальное решение оказывается близким к истинным вторичным источникам. В общем же случае для получения на шаге 2 решений, близких к “истинным”, следует учитывать всю доступную дополнительную информацию о вторичных источниках.

Указанные методы регуляризации обоснованы и апробированы в ряде работ (например, в [12, 13] при решении другой обратной задачи). Отметим, что шаг 2 алгоритма является наиболее трудоемким пунктом при использовании алгоритма. Мы решали СЛАУ (12) с помощью различных вариантов метода регуляризации А.Н. Тихонова [16, 17, 27, 28] и с помощью метода TSVD [18, 28]. Обоснование применения этих методов дано с общих позиций теории регуляризации некорректных задач в [12, 13] для аналогичного интегрального уравнения, и здесь мы не повторяем это обоснование. Наилучшие результаты в расчетах получились для метода TSVD с “обрезанием” (см. [18, 28]) по уровню δ₀^2/3 сингулярных чисел матриц $A\mathrm{<mo>=</mo>}[A_{lj}]\mathrm{<mo>=</mo>}[{\nu }_{lj}\tilde{G}\mathrm{<mo\; stretchy="false">(</mo>}{z}_l-z_j\mathrm{,}\omega \mathrm{,}\widehat{\Omega }\mathrm{<mo\; stretchy="false">)</mo>}]$. Здесь число δ₀ – уровень среднеквадратичной ошибки данных, который будет точно определен ниже. Эти результаты мы и будем приводить в дальнейшем.

В) Шаги 3, 4 не вызывают вопросов. Что касается шага 5, то на нем уравнения v(x, ω) = ξ_k(x)u_k(x, ω), x$\in$X можно было бы решать отдельно для каждой фиксированной частоты ω. Но тогда вычисляемая величина ξ_k(x) будет, вообще говоря, зависеть от частоты. Чтобы избежать этого, следует решать эту совокупность зависящих от ω уравнений как систему с применением МНК.

Г) Реализация шага 6 очевидна.

5. ПРИЛОЖЕНИЕ ПРЕДЛАГАЕМОГО МЕТОДА И РЕШЕНИЕ МОДЕЛЬНЫХ ЗАДАЧ

В этом разделе мы проиллюстрируем решение основной обратной задачи в форме (5) с помощью многократного применения алгоритма решения вспомогательной задачи.

Первая модельная обратная задача заключается в нахождении по данным зондирования u_k(x, ω), k = 1, 2, …, n, измеряемым в области Y, траектории движущего в области X малого (по сравнению с размером области) тела. Опишем более детально параметры модели. Считаем, что все уравнения модели (т.е. (1), (3), (4), (6) и их следствия (9), (12)) записаны в безразмерной форме, так что c₀ = 1 и k₀ = ω. Области П, П_X и П_Y имеют вид:

$\begin{array}{l}\begin{array}{c}\Pi \mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}-\mathrm{10,10<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>}-\mathrm{10,10<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>}-\mathrm{0.6,6.5<mo\; stretchy="false">]</mo>,}\\ {\Pi }_X\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}-\mathrm{10,10<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>}-\mathrm{10,10<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>}-\mathrm{0.6,1.5<mo\; stretchy="false">]</mo>,}\end{array}\hfill \\ {\Pi }_Y\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}-\mathrm{10,10<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>}-\mathrm{10,10<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>6.01,6.5<mo\; stretchy="false">]</mo>.}\hfill \end{array}$

Пространственное распределение источника f(x) представлено дискретной аппроксимацией δ-функции: f(x) = δ(x – L), где L – отрезок прямой x = 0, y$\in$[–5,5], z = 6, т.е. источник локализован вдоль этого отрезка. Временная зависимость источника определяется формулой g(t) = e^–20tsin(20πt), t$\in$[0, Δt], Δt = 1. Рассматривается 20 импульсов такого источника (n = 20). Для того чтобы выяснить скорость работы алгоритма решения обратной задачи на одной частоте, “зондирование” проводится при единственном ω = 3. Модельное решение обратной задачи задается в форме “узкого” трехмерного сферически симметричного гауссиана δc(x, t), наложенного на фон:

$\begin{array}{c}c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{c}_0+\delta c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>;</mo>}\delta c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\\ =A\exp \left(-\frac{\mathrm{<mo>|</mo>}x-h\mathrm{<mo\; stretchy="false">(</mo>}t{\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}^2}{{\sigma }^2}\right)\mathrm{,}A\mathrm{<mo>=</mo>0.4,}\sigma \mathrm{<mo>=</mo>0.3,}\end{array}$

и движущегося по спирали h(t) = (η_x, η_y, η_z):

$\begin{array}{c}{\eta }_x\mathrm{<mo>=</mo>}-4+\frac{7}{20}t\mathrm{,}{\eta }_y\mathrm{<mo>=</mo>2.5}\cos \left(\frac{\pi t}{10}\right)\mathrm{,}\\ {\eta }_z\mathrm{<mo>=</mo>0.5}+0.65\sin \left(\frac{\pi t}{10}\right)\mathrm{,}t\in \mathrm{<mo\; stretchy="false">[</mo>0,20<mo\; stretchy="false">]</mo>.}\end{array}$

Геометрическая схема зондирования вместе с модельным решением показана на рис. 3.

Рис. 3. Модельная геометрическая схема зондирования. Движущаяся неоднородность представлена в дискретные моменты времени как поверхности уровня функции δc(x, t) (на уровне 0.5 от ее максимального значения). Звездочками условно показано расположение источника.

Данные модельного зондирования u(x, ω), x$\in$П_Y для решения обратной задачи были получены путем решения прямой задачи типа (2) с указанными параметрами в области П×[0, 20] по методу конечных элементов. Для их обработки при решении обратной задачи применялся предложенный в разделе 4 алгоритм. В алгоритме использовались равномерные сетки в областях П_X и П_Y размеров N×N×M и N×N×M' соответственно с N = 160, M = 61, M' = 41.

В некоторых численных экспериментах данные для решения обратной задачи задавались с модельной ошибкой измерений. В наших расчетах это делалось наложением на функцию u(x, ω), x$\in$П_Y аддитивной нормально распределенной псевдослучайной помехи с нулевым среднем так, что получаемая в итоге приближенная функция u_δ(x, ω) удовлетворяла бы условию

${‖u_{\delta }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}-u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}‖}_{L_2({\Pi }_Y)}\le \delta {‖u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\omega \mathrm{<mo\; stretchy="false">)</mo>}‖}_{L_2({\Pi }_Y)}\stackrel{\text{def}}{\mathrm{<mo>=</mo>}}{\delta }_0\mathrm{.}$

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

На рис. 4а показаны в сравнении точное решение обратной задачи и приближенное ее решение c_calc(x, t) для точных данных (δ = 0), полученное с помощью предлагаемого алгоритма.

Далее исследовалось влияние ошибки данных на качество получаемого решения. Результаты решения обратной задачи для разных уровней ошибки δ показаны на рис. 4б. Сравнение полученных приближенных решений показывает их значительную чувствительность к уровню возмущения данных.

Рис. 4. (а) – Сравнение точного решения и приближенного решения для точных данных (δ = 0). (б) – Приближенные решения для возмущенных данных с δ = 10^–6 и δ = 10^–4.

Представляет интерес нахождение траектории движущего объекта как линии перемещения его “центра масс”, положение которого моделируется для каждого фиксированного t значением x, усредненным по вычисленному относительному возмущению поля скоростей $\delta c_{\text{calc}}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}=\frac{c_{\text{calc}}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-c_0}{{‖c_{\text{calc}}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-c_0‖}_{L_1({\Pi }_X)}}$, т.е. вектором

$(x_M(t),y_M(t),z_M(t))={\displaystyle \underset{{\Pi }_X}{\int }(x,y,z)}\delta c_{\text{calc}}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}dx.$

Соответствующие результаты показаны на рис. 5. Ошибка в нахождении траекторий вычислялась по правилу:

${\Delta }_{tr}\mathrm{<mo>=</mo>}\frac{{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^n}\left|\eta \mathrm{<mo\; stretchy="false">(</mo>}{t}_k\mathrm{<mo\; stretchy="false">)</mo>}-{\eta }_{\text{appr}}\mathrm{<mo\; stretchy="false">(</mo>}{t}_k\mathrm{<mo\; stretchy="false">)</mo>}\right|}{L\mathrm{<mo\; stretchy="false">(</mo>}\eta \mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}$

где η(t) – точная траектория, η_appr(t) – приближенно найденная траектория, а L(η) – длина точной траектории. Отметим, что даже в случае δ = 0 ошибка данных содержит часть, связанную с конечномерной аппроксимацией решаемой обратной задачи. Именно по этой причине величина Δ_tr оказывается не нулевой при δ = 0.

Рис. 5. (а) – Сравнение точной траектории (сплошная линия) и найденных приближенных траекторий (линия с кружками): (а) – для невозмущенных данных (δ = 0), (б) – для возмущенных данных с δ = 10^–6.

Вторая модельная задача, которую мы представим, заключается в нахождении по данным зондирования движущегося прямолинейно локального распределения скоростей в форме гауссиана, который распределен по тору. Все параметры зондирования сохраняются без изменений. Модельное решение в данном случае имеет вид c(x, t) = c₀ + + 0.4exp(–5|R(x – η)|), где функция

$\begin{array}{c}R\left(\mathbf{x}\mathbf{\right)}=({(\mathrm{x}-0.2)}^2+{(y+0.3)}^2+z^2+\\ +R_0^2-r_0^2{)}^2-4R_0^2({(x-0.2)}^2+z^2)\end{array}$

определяет в уравнении R(x) = 0 поверхность тора с параметрами R₀ = 0.3, r₀ = 0.1, а вектор-функция

$\begin{array}{l}\eta \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left({\eta }_x\mathrm{,}{\eta }_y\mathrm{,}{\eta }_z\right)\mathrm{,<mi>\; </mi>}{\eta }_x\mathrm{<mo>=</mo>}-4+7t\mathrm{,}\\ {\eta }_y\mathrm{<mo>=</mo>2.5<mo\; stretchy="false">(</mo>1}+t\mathrm{<mo\; stretchy="false">)</mo>,}{\eta }_z\mathrm{<mo>=</mo>0,<mi>\; </mi>}t\in [\mathrm{0,20<mo\; stretchy="false">]</mo>}\end{array}$

задает прямую, по которой происходит сдвиг неоднородности.

На рис. 6 дается сравнение точного решения обратной задачи и ее приближенного решения для точного задания данных. Видно, что приближенное решение достаточно хорошо с качественной точки зрения отражает особенности точного решения. Для этой модельной задачи также проведено исследование влияния различного уровня возмущений данных на точность получаемого решения. Результаты показаны на рис. 7. Так же, как и для первой модельной задачи, можно видеть значительную чувствительность приближенных решений к уровню помех данных. Поэтому в практических приложениях предлагаемой методики следует вводить в постановку обратной задачи адекватные априорные ограничения на искомое решение и соответственно модифицировать алгоритм из раздела 4. Однако это требует дополнительных исследований, выходящих за рамки данной статьи.

Рис. 6. Сравнение (а) – точного и (б) – приближенного решения для точных данных.

Рис. 7. Приближенные решения для разных уровней ошибок данных: (а) – δ = 10^–6, (б) – δ = 10^–4.

6. ОБСУЖДЕНИЕ СВОЙСТВ ПРЕДЛАГАЕМОГО АЛГОРИТМА РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ И ВЫВОДЫ

Все вычисления проводились в системе МАТЛАБ на ПК с процессором Intel (R) Core (TM) i7-7700 CPU 3.60 GHz, ОЗУ 16Гб. Алгоритм решения трехмерной обратной задачи оказался достаточно быстрым: нахождение трехмерной функции c(x, t) вида (5) для одного значения t (т.е. на интервале типа (t_k–1, t_k)) требует около 10–15 с для пространственных сеток указанного выше размера.

Решаемая обратная задача весьма чувствительна к ошибкам входных данных, как можно видеть при рассмотрении рис. 5, 8. Это связано с весьма быстрым убыванием сингулярных чисел матриц [A_lj] = [ν_ljG(z_l – z_j, ω, Ω)] СЛАУ (12), решаемых на шаге 2 алгоритма. Такое убывание является специфической особенностью решаемой обратной задачи, а точнее, ядра интегрального уравнения (9), т.е. функции Грина прямой задачи. Аналогичное свойство обратной коэффициентной задачи для волнового уравнения отмечалось ранее в работах [12, 13]. Подтверждающие это соответствующие теоретические оценки устойчивости при различных априорных предположениях на точное решение можно найти в [2, 3].

Тем не менее, при достаточно малых возмущениях данных алгоритм для рассматриваемой трехмерной обратной задачи позволяет достаточно точно решить ее за разумное время на ПК средней производительности.

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

1. Рассматриваемая трехмерная обратная задача акустического зондирования нестационарной среды (в рамках использованной модели) может быть решена для достаточно мелких сеток за сравнительно малое время с помощью предложенного алгоритма. Алгоритм позволяет достаточно надежно определять положения и форму траектории движения небольших неоднородностей акустической среды при данных с малыми ошибками. Алгоритм может быть модифицирован для задач зондирования нестационарной среды в цилиндрической или сферической области по аналогии с работой [26].

2. Исходная обратная задача, основанная на уравнениях (4) и рассматриваемая для конечного набора частот, имеет, вообще говоря, не единственное решение. Для выделения определенного решения необходимо использовать дополнительную информацию об искомой функции c(x, t) или о связанных с ней величинах. При отсутствии значимой информации такого рода в нашем алгоритме используется предположение (5) о ее виде и поиск нормального решения интегрального уравнения (7), т.е. минимального значения нормы величины v(x', ω) = ξ_k(x')u_k(x', ω).

В случае, когда уравнение (7) имеет единственное решение при рассматриваемом наборе частот, это решение совпадает с найденным нормальным.

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

Работа выполнена при финансовой поддержке РНФ (проект № 22-71-10070) для первого автора и Программы повышения конкурентоспособности Национального исследовательского ядерного университета МИФИ (проект 02.a03.21.0005 от 27.08.2013) для второго автора. Основные результаты разделов 1–3 получены за счет гранта РНФ.