Correlation Iteration Method of Acoustic Tomography with Incoherent Field Sources

全文:

1. Введение

Пассивная акустическая томография – это метод определения свойств исследуемой среды, не использующий активные источники звука. Он основан на регистрации и последующей обработке акустического излучения от присутствующих в среде источников. Задачи такого рода возникают, например, в медицинской томографии [1–3], где акустическое поле порождается тепловым движением частиц; в геоакустике [4, 5], где его источниками служат микросейсмы; в гидроакустике [6, 7], где оно включает естественные и антропогенные шумы различной природы, в гелиосейсмологии [8, 9]. Отсутствие активных источников звука, во-первых, существенно удешевляет такие исследования, а во-вторых, является более предпочтительным с точки зрения вопросов экологии (в гео- и гидроакустике) или уменьшения воздействия на организм (в медицине). Вместе с тем, методы пассивной томографии оказываются сложнее, чем методы активной томографии. Это связано с необходимостью оценки не только искомых свойств среды, но и описания создающих акустическое поле источников. Такая комбинированная обратная задача излучения-рассеяния представляется сильно недоопределенной, и для ее успешного решения целесообразно наложить те или иные дополнительные ограничения на класс возможных источников поля и допустимые диапазоны свойств среды.

Одним из методов пассивной акустической томографии является метод шумовой интерферометрии. В его основе лежит предположение о том, что акустическое поле создается некогерентными источниками, пространственная плотность мощности которых изотропна по отношению к системе наблюдения. Тогда корреляционная функция сигналов, зарегистрированных в двух точках среды, оказывается пропорциональной разности запаздывающей и опережающей функций Грина этой среды с аргументами в виде координат выбранных двух точек. Корреляционную функцию можно определить на основе экспериментальных данных. Далее она используется для постановки и решения обратной задачи рассеяния. В качестве недостатка метода шумовой интерферометрии можно указать ограничения, накладываемые на область его применимости. Первым ограничивающим фактором является требование инвариантности описывающих среду уравнений по отношению к преобразованию обращения времени, которое нарушается в присутствии течений или поглощения звука. Эти вопросы обсуждались в работах [10–12]. Второй фактор связан с тем, что метод перестает работать, если мощность шума, приходящего с некоторого направления, выше, чем с других. Поскольку такая анизотропия шума встречается часто, это требует осторожности при практическом применении.

В случае, когда анизотропия шумового поля вызвана наличием небольшого числа (по сравнению с числом элементов приемной антенной решетки) дополнительных источников, удается “восстановить” работоспособность метода шумовой интерферометрии [13]. Для этого корреляционной обработке подвергаются все доступные в эксперименте пары сигналов, в результате чего в каждой относительно узкой полосе частот вычисляются их матрицы когерентности. Для каждой из них осуществляется корректировка максимальных собственных чисел так, чтобы минимизировать дисперсию диагональных элементов итоговых матриц после коррекции. Схожая процедура описывалась в [14, 15], где она использовалась для улучшения видимости слабых источников сигнала на фоне более мощных. Физически ее применение означает такую фазировку приемной антенной решетки, при которой сигналы из области пространственной локализации дополнительных источников, создающих анизотропную часть шумового поля, подавляются до среднего уровня шумов. В результате получается набор корреляционных функций, близких по значениям к искомым разностям функций Грина.

В случае, когда число дополнительных источников соизмеримо или превосходит количество элементов приемной антенной решетки, описанная процедура коррекции собственных чисел не приводит к успеху. Тогда предлагается отказаться от использования метода шумовой интерферометрии и рассмотреть задачу акустической томографии в следующей постановке, которая является предметом настоящей статьи. Предполагается, что шумовое поле в среде создается стационарными фоновыми источниками с неизвестной пространственной плотностью мощности $I^{(0)}(r).$ Дополнительно к ним в среде могут появляться независимые контролируемые некогерентные источники с известной пространственной плотностью мощности $I^{(n)}(r).$ Термин “контролируемые” понимается в том смысле, что в эксперименте доступны как сигналы, записанные при наличии лишь фоновых источников $I^{(0)}(r),$ так и сигналы, записанные при наличии источников обоих типов, т.е. с пространственной плотностью мощности $I^{(0)}(r)+I^{(n)}(r).$ Верхний индекс $n,$, который в дальнейшем будет называться ракурсом облучения, при этом означает номер пространственной конфигурации контролируемых некогерентных источников, и при $n=0$ такие источники отсутствуют.

В некоторых случаях описанная постановка задачи может быть формально отнесена к пассивной томографии. Например, в приложениях гидроакустики может рассматриваться акустическое поле в акватории, по которой перемещаются суда, причем их координаты и мощность создаваемого шума известны. Изменение положения судов с течением времени обеспечивает перебор ракурсов облучения; при этом предполагается, что свойства среды за время наблюдений меняются слабо. В других случаях контролируемые источники могут быть созданы специально, что отвечает задаче активной или активно-пассивной томографии. Например, в медицине может идти речь о использовании собственного теплового акустического излучения биологических сред совместно с некогерентным излучением дополнительной “подсветки” [2, 3]. Сделав такое замечание, дальнейшее рассмотрение можно вести независимо от природы контролируемых источников.

2. Корреляционный итерационный метод

Рассматривается рассеяние акустического поля на неоднородности, занимающей внутри всего координатного пространства ${\text{IR}}^d$ размерности $d$ область $\Re$ конечных размеров. Вне $\Re$ находится однородная фоновая среда. Предполагается, что сигналы сосредоточены в узкой полосе частот с центральной частотой ${\omega }_0$, и временнáя зависимость выбирается в виде $~\exp (-i{\omega }_{0}t)$. Координатные зависимости скорости звука $c(r)$ и амплитудного коэффициента поглощения $\alpha (r)$ описываются с помощью комплексного волнового числа $k(r)\equiv \omega /c(r)+i\alpha (r)$. Снаружи неоднородности, т.е. внутри фоновой среды, оно известно: $k(r)\equiv k_0$ при $r\in {\text{IR}}^d\backslash \Re$. Кроме этого, считается известной функция $k_{\text{a}}(r)$, которая равна $k_0$ вне неоднородности, а внутри нее выражает априорную оценку свойств этой неоднородности. Такое предположение не снижает общность рассмотрения, т.к. в отсутствие дополнительной информации можно положить $k_{\text{a}}(r)\equiv k_0$ во всем пространстве.

Потенциал $\phi (r)$ акустического поля удовлетворяет уравнению Гельмгольца:

$\Delta \phi (r)+k_{}^2(r)\phi (r)=F(r),$ $r\in {\text{IR}}^d,$ (1)

где $F(r)$ – источники поля. Его решение, с учетом условия излучения на бесконечности, имеет вид

$\phi (r)={\displaystyle \underset{{\text{IR}}^d}{\int }G(r,r^{\prime })F(r^{\prime })d{r}^{\prime }},$ $r,r^{\prime }\in {\text{IR}}^d,$ (2)

где $G(r,r^{\prime })$ – запаздывающая функция Грина рассматриваемой неоднородной среды. Аналогичное уравнение Гельмгольца можно записать для поля ${\phi }_{\text{a}}(r)$, которое создают те же самые источники в среде с волновым числом $k_{\text{a}}(r)$:

$\Delta {\phi }_{\text{a}}(r)+k_{\text{a}}^2(r){\phi }_{\text{a}}(r)=F(r),$ $r\in {\text{IR}}^d.$ (3)

В этом случае вводится запаздывающая функция Грина $G_{\text{a}}(r,r^{\prime })$, и решение уравнения (3) имеет вид ${\phi }_{\text{a}}(r)={\displaystyle \underset{{\text{IR}}^d}{\int }{G}_{\text{a}}(r,r^{\prime })}F(r^{\prime })d{r}^{\prime }$. Если добавить в обе части уравнения (1) слагаемое $\epsilon (r)\phi (r)$, где $\epsilon (r)\equiv k_{\text{a}}^2(r)-k_{}^2(r)$ – функция рассеивателя относительно среды с волновым числом $k_{\text{a}}(r),$ получится уравнение (3) с источниками поля, равными $F(r)+\epsilon (r)\phi (r)$. Его формальное решение с помощью функции Грина приводит к хорошо известному в квантовой теории поля уравнению Липпмана–Швингера [16]. Оно часто служит основой для решения обратных задач и в акустике [17, 18]:

$\phi (r)={\phi }_{\text{a}}(r)+{\displaystyle \underset{\Re }{\int }{G}_{\text{a}}(r,r^{\prime })\epsilon (r^{\prime })\phi (r^{\prime })d{r}^{\prime }},$   

$r\in {\text{IR}}^d,$ $r^{\prime }\in \Re .$  (4)

Это уравнение отражает тот факт, что полное поле источника $\phi (r)$ складывается из падающего поля ${\phi }_{\text{a}}(r)$, которое он создавал бы в среде с волновым числом $k_{\text{a}}(r)$, и рассеянного поля, выраженного интегралом по области неоднородности $\Re$, где функция $\epsilon (r)$ отлична от нуля. Учитывая, что (4) выполнено при любой конфигурации источников поля, и полагая $F(r)=\delta (r-r^{″})$ при произвольном $r^{″}$, можно записать такое же по форме уравнение относительно функции Грина $G(r,r^{\prime })$:

$G(r,r^{\prime })=G_{\text{a}}(r,r^{\prime })+{\displaystyle \underset{\Re }{\int }{G}_{\text{a}}(r,r^{″})\epsilon (r^{″})G(r^{″},r^{\prime })d{r}^{″}},$

$r,r^{\prime }\in {\text{IR}}^d,$ $r^{″}\in \Re .$ (5)  

Его можно использовать для вычисления функции $G(r,r^{\prime })$, если известна функция $G_{\text{a}}(r,r^{\prime })$. В свою очередь, функцию $G_{\text{a}}(r,r^{\prime })$ можно определить из того же уравнения (5), выполнив замену переменных и рассматривая теперь в качестве априорной оценки однородную среду с волновым числом $k_0$ и функцией Грина $G_{\text{0}}(r,r^{\prime })$:

$G_{\text{a}}(r,r^{\prime })=G_{\text{0}}(r,r^{\prime })+{\displaystyle \underset{\Re }{\int }{G}_{\text{0}}(r,r^{″}){\epsilon }_{\text{a}}(r^{″})G_{\text{a}}(r^{″},r^{\prime })d{r}^{″}},$

$r,r^{\prime }\in {\text{IR}}^d,$ $r^{″}\in \Re .$ (6)

Здесь ${\epsilon }_{\text{a}}(r)\equiv k_{\text{0}}^2-k_{\text{a}}^2(r)$ – функция рассеивателя для априорной оценки $k_{\text{a}}(r)$ относительно фоновой среды. Функция Грина $G_{\text{0}}(r,r^{\prime })$ известна аналитически [2, 17–19]. Например, в двумерном случае ( $d=2$ ), который будет рассматриваться в дальнейшем на этапе численного моделирования, $G_0^{(d=2)}(r,r^{\prime })=-\frac{i}{4}{H}_0^{(1)}\left(k_0\left|r-r^{\prime }\right|\right),$ где $H_0^{(1)}(\cdot )$ – функция Ханкеля первого рода нулевого порядка. В случае, когда априорная оценка $k_{\text{a}}(r)$ близка к истинному значению $k(r)$, $G(r,r^{\prime })\approx G_{\text{a}}(r,r^{\prime })$, и для уравнения (5) справедливо борновское приближение [17, 18, 20, 21]:

$G(r,r^{\prime })\approx G_{\text{a}}(r,r^{\prime })+{\displaystyle \underset{\Re }{\int }{G}_{\text{a}}(r,r^{″})\epsilon (r^{″})G_{\text{a}}(r^{″},r^{\prime })d{r}^{″}},$

$r,r^{\prime }\in {\text{IR}}^d,$ $r^{″}\in \Re .$ (7)

Далее рассматриваются акустические поля, порожденные шумовыми источниками $F^{(n)}(r)$ при каждом $n$ -м ракурсе облучения. Регистрация сигналов осуществляется с помощью $N_R$ приемников. Каждый $j$ -й из них можно задать, определив множество точек $S_j$ его поверхности. Тогда принятые сигналы $U_j^{(n)}$ представляются в виде

$U_j^{(n)}={\displaystyle \underset{S_j}{\int }{\phi }^{(n)}(r)dr},$ (8)

где ${\phi }^{(n)}(r)$, $r\in {\text{IR}}^d$ – поле акустического потенциала при соответствующем ракурсе облучения. На их основе формируется набор матриц когерентности, элементы которых определяются как ${\Gamma }_{ij}^{(n)}=⟨U_i^{(n)}{\left\{U_j^{(n)}\right\}}^{*}⟩$, где “*” означает комплексное сопряжение, а угловые скобки – усреднение по всем реализациям. Последовательно учитывая (8) и (2), а также принимая во внимание, что функции Грина являются детерминированными, и их можно вынести из-под операции усреднения, для матриц когерентности можно получить выражение $\begin{array}{c}{\Gamma }_{ij}^{(n)}=\\ ={\displaystyle \underset{S_i}{\int }d{r}^{\prime }{\displaystyle \underset{S_j}{\int }d{r}^{″}{\displaystyle \underset{{\text{IR}}^d}{\int }d{r}_1{\displaystyle \underset{{\text{IR}}^d}{\int }d{r}_{2}G(r^{\prime },r_1)G^{*}(r^{″},r_2)}{\Gamma }_F^{(n)}(r_1,r_2)}}}.\end{array}$ (9)

Здесь ${\Gamma }_F^{(n)}(r_1,r_2)\equiv ⟨F^{(n)}(r_1){\left\{F^{(n)}(r_2)\right\}}^{*}⟩$ – функция когерентности источников поля, расположенных в двух произвольных точках $r_1,r_2\in {\text{IR}}^d$. Поскольку фоновые источники и контролируемые источники предполагаются независимыми, и, кроме того, контролируемые источники некогерентны, выполняется равенство

${\Gamma }_F^{(n)}(r_1,r_2)={\Gamma }_F^{(0)}(r_1,r_2)+I^{(n)}(r_1)\delta (r_1-r_2),$ $r_1,r_2\in {\text{IR}}^d,$

где ${\Gamma }_F^{(0)}(r_1,r_2)$ и $I^{(n)}(r_1)\delta (r_1-r_2)$ – функции когерентности фоновых и контролируемых источников поля, соответственно. Тогда разность матриц когерентности при $n$ -м и нулевом ракурсах облучения зависит лишь от мощности контролируемых источников $I^{(n)}(r)$: ${\Gamma }_{ij}^{(n)}-{\Gamma }_{ij}^{(0)}={\displaystyle \underset{S_i}{\int }{\displaystyle \underset{S_j}{\int }{\Gamma }_G^{(n)}(r^{\prime },r^{″})}d{r}^{\prime }d{r}^{″}},$ (10)

где

${\Gamma }_G^{(n)}(r^{\prime },r^{″})\equiv {\displaystyle \underset{\Omega }{\int }G(r^{\prime },r)G^{*}(r^{″},r)}{I}^{(n)}(r)dr,$

$r^{\prime },r^{″}\in {\text{IR}}^d.$ (11)

Величина ${\Gamma }_G^{(n)}(r^{\prime },r^{″})$ представляет собой функцию когерентности полей контролируемых источников, которые измеряются в точках $r^{\prime }$ и $r^{″}$. Ее можно представить с помощью уравнения (5) в виде суммы четырех слагаемых:

$\begin{array}{c}{\Gamma }_G^{(n)}(r^{\prime },r^{″})={\displaystyle \underset{\Omega }{\int }dr{I}^{(n)}(r)G_{\text{a}}(r^{\prime },r)}{G}_{\text{a}}^{\text{*}}(r^{″},r)+\\ +{\displaystyle \underset{\Re }{\int }d{r}_1{G}_{\text{a}}^{\text{*}}(r^{″},r_1){\epsilon }^{*}(r_1){\displaystyle \underset{\Omega }{\int }dr{I}^{(n)}(r)}{G}_{\text{a}}(r^{\prime },r)G^{*}(r_1,r)}+\\ +{\displaystyle \underset{\Re }{\int }d{r}_1{G}_{\text{a}}(r^{\prime },r_1)\epsilon (r_1){\displaystyle \underset{\Omega }{\int }dr{I}^{(n)}(r)}G(r_1,r)}{G}_{\text{a}}^{\text{*}}(r^{″},r)+\\ +{\displaystyle \underset{\Re }{\int }d{r}_1{G}_{\text{a}}(r^{\prime },r_1)\epsilon (r_1)}{\displaystyle \underset{\Re }{\int }d{r}_2{G}_{\text{a}}^{\text{*}}(r^{″},r_2){\epsilon }^{*}(r_2)}\times \\ \times {\displaystyle \underset{\Omega }{\int }dr{I}^{(n)}(r)}G(r_1,r)G^{*}(r_2,r),\end{array}$

где $r^{\prime },r^{″}\in {\text{IR}}^d$. Первое слагаемое

${\Gamma }_{G\text{a}}^{(n)}(r^{\prime },r^{″})={\displaystyle \underset{\Omega }{\int }dr{I}^{(n)}(r)G_{\text{a}}(r^{\prime },r)}{G}_{\text{a}}^{\text{*}}(r^{″},r)$ (12)

известно, поскольку оно зависит от известной мощности $I^{(n)}(r)$ и функций Грина $G_{\text{a}}(r,r^{\prime }),$ которые могут быть рассчитаны с помощью (6). Остальные слагаемые выражают нелинейную зависимость от искомой функции $\epsilon (r)$. Такая нелинейность присутствует явно в последнем слагаемом, содержащем произведение двух таких функций, и неявно во втором и третьем слагаемых, поскольку в них входит неизвестная функция Грина исследуемой среды. В этом смысле уравнение (11) схоже с уравнением Липпмана–Швингера (5), которое также нелинейно относительно искомой функции $\epsilon (r).$ Для решения уравнения (11) предлагается применить итерационный подход, который нередко используется [17, 18, 22, 23] при решении уравнения типа (5). Он состоит в последовательном улучшении априорной оценки функции $k(r),$ для чего в уравнениях (5) или (11) используется линеаризующее их борновское приближение (7); после этого функция $G_{\text{a}}(r,r^{\prime })$ пересчитывается согласно (6). В борновском приближении $G(r,r^{\prime })\approx G_{\text{a}}(r,r^{\prime }),$ и поэтому

$\begin{array}{c}{\Gamma }_{ij}^{(n)}-{\Gamma }_{ij}^{(0)}-{\displaystyle \underset{S_i}{\int }d{r}^{\prime }{\displaystyle \underset{S_j}{\int }d{r}^{″}{\Gamma }_{G\text{a}}^{(n)}(r^{\prime },r^{″})}}\approx \\ \approx {\displaystyle \underset{S_i}{\int }d{r}^{\prime }{\displaystyle \underset{S_j}{\int }d{r}^{″}{\displaystyle \underset{\Re }{\int }dr\left[\begin{array}{l}\epsilon (r)G_{\text{a}}(r^{\prime },r){\Gamma }_{G\text{a}}^{(n)}(r,r^{″})+\\ +{\left\{\epsilon (r)G_{\text{a}}^{}(r^{″},r){\Gamma }_{G\text{a}}^{(n)}(r,r^{\prime })\right\}}^{*}\end{array}\right]}}}.\end{array}$(13)

Сходство уравнений (7) и (13) легко прослеживается: роль падающего поля выполняет функция когерентности ${\Gamma }_{G\text{a}}^{(n)}(r^{\prime },r^{″})$; роль полного поля, регистрируемого на антенной решетке, выполняет разность ${\Gamma }_{ij}^{(n)}-{\Gamma }_{ij}^{(0)}$. Для решения (13) надо предварительно вычислить ${\Gamma }_{G\text{a}}^{(n)}(r^{\prime },r^{″})$ с помощью (12), что на практике не очень удобно. Вместо этого вводятся обозначения $Q_{\text{a}}(r^{\prime },r^{″},r)\equiv G_{\text{a}}(r^{\prime },r)G_{\text{a}}^{\text{*}}(r^{″},r)$ и $P_{\text{a}}(r^{\prime },r^{″},r_1,r)\equiv G_{\text{a}}(r^{\prime },r_1)G_{\text{a}}(r_1,r)G_{\text{a}}^{\text{*}}(r^{″},r)$ для произведений двух и трех функций Грина соответственно, и производится перегруппировка множителей:

$\begin{array}{c}{\Gamma }_{ij}^{(n)}-{\Gamma }_{ij}^{(0)}-{\displaystyle \underset{S_i}{\int }d{r}^{\prime }{\displaystyle \underset{S_j}{\int }d{r}^{″}}}{\displaystyle \underset{\Omega }{\int }dr{I}^{(n)}(r)}{Q}_{\text{a}}(r^{\prime },r^{″},r)\approx \\ \approx {\displaystyle \underset{S_i}{\int }d{r}^{\prime }{\displaystyle \underset{S_j}{\int }d{r}^{″}}}{\displaystyle \underset{\Omega }{\int }dr{I}^{(n)}(r)}{\displaystyle \underset{\Re }{\int }d{r}_1\left[\begin{array}{c}\epsilon (r_1)P_{\text{a}}(r^{\prime },r^{″},r_1,r)+\\ +{\epsilon }^{*}(r_1)P_{\text{a}}^{\text{*}}(r^{″},r^{\prime },r_1,r)\end{array}\right].}\end{array}$ (14)

Можно видеть, что уравнение (14) линейно относительно действительной ${\epsilon }^{\prime }(r)\equiv \mathrm{Re}\epsilon (r)$ и мнимой ${\epsilon }^{″}(r)\equiv \mathrm{Im}\epsilon (r)$ частей функции $\epsilon (r)$, которые рассматриваются как независимые переменные. После дискретизации исследуемой среды интегрирование заменяется матричным произведением, и уравнение (14) сводится к простой системе линейных уравнений. Функции $P_{\text{a}}$ в правой части отличаются порядком аргументов и комплексным сопряжением, т.е. соответствующие матрицы – гильбертово сопряженные.

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

1.При нулевом и при каждом $n$ -м ракурсе облучения в эксперименте определяются значения матричных элементов ${\Gamma }_{ij}^{(n)}$, и формируются разности ${\Gamma }_{ij}^{(n)}-{\Gamma }_{ij}^{(0)}$. Первая оценка полагается равной ${\widehat{\epsilon }}_1(r)\equiv {\epsilon }_{\text{a}}(r)=k_{\text{0}}^2-k_{\text{a}}^2(r)$.

2.Организуется бесконечный цикл, в рамках которого на каждой $m$ -й итерации, начиная с $m=1$, делаются следующие действия.

а) Имеющаяся в наличии оценка функции рассеивателя (r) ${\widehat{\epsilon }}_m(r)=k_{\text{0}}^2-{\widehat{k}}_m^2(r)$, где ${\widehat{k}}_m(r)$ – это соответствующая ей оценка волнового числа, подставляется в уравнение (6) вместо величины ${\epsilon }_{\text{a}}(r)$, и рассчитывается функция Грина $G_{\text{a}}(r,r^{\prime })$.

б) С ее помощью определяются входящие в уравнение (14) произведения $Q_{\text{a}}(r^{\prime },r^{″},r)$ и $P_{\text{a}}(r^{\prime },r^{″},r_1,r)$.

в) Уравнение (14) решается относительно функции $\epsilon (r)={\widehat{k}}_m^2(r)-k^2(r)$ . Поскольку выполнено тождество $k_{\text{0}}^2-k^2(r)=\left(k_{\text{0}}^2-{\widehat{k}}_m^2(r)\right)+\left({\widehat{k}}_m^2(r)-k^2(r)\right)={\widehat{\epsilon }}_m(r)+\epsilon (r)$ $k_{\text{0}}^2-k^2(r)=\left(k_{\text{0}}^2-{\widehat{k}}_m^2(r)\right)+\left({\widehat{k}}_m^2(r)-k^2(r)\right)={\widehat{\epsilon }}_m(r)+\epsilon (r)$, следующая оценка функции ${\widehat{\epsilon }}_{m+1}(r)$ полагается равной ${\widehat{\epsilon }}_{m+1}(r)={\widehat{\epsilon }}_m(r)+\epsilon (r)$.

3.Чтобы определить, насколько оценка ${\widehat{\epsilon }}_m(r)$ соответствует входным данным, вводится безразмерная величина ${\delta }_{\Gamma }(\widehat{\epsilon })$ (), равная

Цикл заканчивается, когда она становится ниже заданного порогового значения. При этом относительная невязка

$\delta (\widehat{\epsilon })\equiv \sqrt{{\displaystyle \underset{\Re }{\int }{\left|\widehat{\epsilon }(r)-\epsilon (r)\right|}^{2}dr}/{\displaystyle \underset{\Re }{\int }{\left|\epsilon (r)\right|}^{2}dr}}$ (15)

характеризует близость полученной оценки $\widehat{\epsilon }(r)$ к истинной функции рассеивателя $\epsilon (r)$.

Уравнение (14), лежащее в основе корреляционного итерационного метода, можно в некоторых случаях упростить. Так, если контролируемые источники точечные и при каждом $n$ -м ракурсе облучения расположены в точке $r^{(n)}$, то $I^{(n)}(r)=I^{(n)}\delta \left(r-r^{(n)}\right)$, и приемники тоже точечные и имеют координаты $r_j$, то с учетом этого вместо (14) можно записать

$\begin{array}{c}\frac{{\Gamma }_{ij}^{(n)}-{\Gamma }_{ij}^{(0)}}{I^{(n)}}-Q_{\text{a}}\left(r_i,r_j,r^{(n)}\right)\approx \\ \approx {\displaystyle \underset{\Re }{\int }dr\left[\epsilon (r)P_{\text{a}}\left(r_i,r_j,r,r^{(n)}\right)+{\epsilon }^{*}(r_1)P_{\text{a}}^{\text{*}}\left(r_j,r_i,r,r^{(n)}\right)\right].}\end{array}$(16)

Следует обратить внимание, что итерационные алгоритмы решения обратной задачи рассеяния, как правило, предполагают численное решение прямой задачи рассеяния на каждой итерации [17, 18]. В рамках описанного алгоритма она возникает на втором шаге, когда требуется определить $G_{\text{a}}(r,r^{\prime }).$ Погрешности, которые привносятся при этом, оказывают существенное влияние на точность итоговой оценки. В частности, значительный вклад оказывают процессы многократного рассеяния внутри каждого элемента разрешения. Эти вопросы были затронуты в [24]. Альтернативной возможностью является использование для решения прямой задачи метода конечных элементов. В этом случае удается уменьшить объем вычислений и сократить время расчетов, но требуется аккуратное задание условий на границе области наблюдения поля.

3. Численное моделирование работы корреляционного итерационного метода

С целью проверки работоспособности метода моделировалась задача томографирования двумерной неоднородности $\Re$ квадратной формы, расположенной в начале координат и имеющей размеры $5{\lambda }_0\times 5{\lambda }_0$, где ${\lambda }_0\equiv 2\pi /k_0$ – длина волны в окружающей фоновой среде. Скорость звука $c(x,y)$ и амплитудный коэффициент поглощения $\alpha (x,y)$ внутри неоднородности задавались в виде

$c(x,y)=c_0\left(1+{\beta }_{0}s\left(\frac{x}{{\lambda }_0},\frac{y}{{\lambda }_0}\right)\right)$ 

и

$\alpha (x,y)={\alpha }_0{k}_0\left|s\left(\frac{y}{{\lambda }_0},\frac{x}{{\lambda }_0}\right)\right|,$

соответственно, где $c_0=\omega /k_0$ – скорость звука в фоновой среде. Безразмерные параметры ${\alpha }_0$ и ${\beta }_0$ варьируются, обеспечивая изменение контраста неоднородности. Функция $s(x,y)$ задает пространственное распределение параметров внутри неоднородности. Для удобства данные об исследованных в настоящей статье неоднородностях и значения описанных характеристик объединены в таблицу. $N_R=20$ точечных приемников располагались в точках с полярными координатами $r_j=\left(100{\lambda }_0;2\pi (j-1)/N_R\right);$ $N_F=20$ точечных контролируемых источников – в точках с полярными координатами $r^{(n)}=\left(100{\lambda }_0;2\pi (n-0.5)/N_F\right)$. Шаг дискретизации был выбран равным ${\lambda }_0/8$. При моделировании здесь и далее априорная информация не использовалась, т.е. полагалось $k_{\text{a}}(r)=k_0$.

Таблица. Параметры неоднородностей

Чтобы численно определить сложность восстановления каждой конкретной моделируемой неоднородности, вводится несколько ее характеристик, по-разному связанных со свойствами среды. Первая из них – норма амплитуды рассеяния $f(\theta ,{\theta }^{\prime })$ , равная $‖f(\theta ,{\theta }^{\prime })‖\equiv \sqrt{{\displaystyle \underset{0}{\overset{2\pi }{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}{\left|f(\theta ,{\theta }^{\prime })\right|}^2}d\theta d{\theta }^{\prime }}}$. Амплитуда рассеяния определяется в пространствах разной размерности $d$ аналогично тому, как это сделано в ряде работ [18, 27]; $\theta$ и ${\theta }^{\prime }$ – полярные углы, под которыми распространяются падающая и рассеянная волны. При численном расчете $‖f(\theta ,{\theta }^{\prime })‖$ используется уравнение (6).

Вторая характеристика определяется как максимальный по различным траекториям $\mathcal{l}$ дополнительный набег фазы, который получается при распространении сигнала через неоднородность по сравнению с набегом фазы в фоновой среде: $\Delta \varphi =\underset{\mathcal{l}}{\overline{\max }}{\displaystyle \underset{\mathcal{l}}{\int }\left(k_0-\frac{\omega }{c(r)}\right)d\mathcal{l}}$. Черта над $\overline{\max }$ означает, что эта операция в данном случае нестандартная. Она заключается в том, что вначале определяется наибольшее по абсолютной величине значение аргумента, а затем результату приписывается тот же знак, что и у аргумента. Такое сохранение знака позволяет учесть, что набег фазы в среде может быть как положительным, так и отрицательным. Если для оценки полагать каждую траекторию $\mathcal{l}$ прямой линией, то такая процедура эквивалентна вычислению максимума $\Delta \varphi =\overline{\max }\zeta (a,\theta )$ преобразования Радона $\zeta (a,\theta )=\text{Rad}\left(k_0-\frac{\omega }{c(r)}\right)$ от подынтегрального вы- ражения, которое можно рассчитать быстро.

Третья характеристика равна максимальному модулю отношения $M\equiv \underset{r,r^{\prime }}{\max }\left|\frac{G(r,r^{\prime })-G_{\text{a}}(r,r^{\prime })}{G_{\text{a}}(r,r^{\prime })}\right|$ рассеянного поля $G(r,r^{\prime })-G_{\text{a}}(r,r^{\prime })$ к падающему полю $G_{\text{a}}(r,r^{\prime })$. В зависимости от того, является ли она много меньшей единицы, меньшей единицы или превосходит единицу, неоднородность относят к классу слабых, средних или сильных рассеивателей соответственно [17]. Сильные рассеиватели наиболее сложны для восстановления, и итерационные процедуры, в том числе предлагаемый в настоящей статье метод, не гарантируют успеха. Более того, в двумерном случае, когда объем исходных данных конечен, в классе сильных рассеивателей решение задачи не единственно [18].

На рис. 1а представлено пространственное распределение относительной скорости звука $c(x,y)/c_0$ для неоднородности № 1 с ${\alpha }_0=0.1;$ ${\beta }_0=0.3$. Она описывается функцией

$\begin{array}{c}{s}_1(x,y)\equiv \exp \left(-\frac{{\left(\left|x+iy-0.5\right|-1\right)}^2}{2\times {0.25}^2}\right)-\\ -\frac{1}{2}\exp \left(-\frac{{\left(\left|x+iy+0.5\right|-1\right)}^2}{2\times {0.25}^2}\right),\end{array}$

представляющей собой сумму двух “колец” и имеющей достаточно сложный вид. Форма пространственного распределения коэффициента поглощения такая же: она получается из формы распределения скорости звука симметрией относительно прямой $y=x$. Несмотря на небольшой волновой размер, она относится к классу сильных рассеивателей. Результат ее восстановления с помощью корреляционного итерационного метода приведен на рис. 1б. Здесь представлены выполненные вдоль отрезка $AB$ и нормированные на $k_0^2$ сечения полученной после 500 итераций оценки ${\widehat{\epsilon }}_{500}(x,y)$ и искомой функции рассеивателя $\epsilon (x,y).$ Хорошее визуальное согласование кривых в сочетании с низким достигнутым значением $\delta ({\widehat{\epsilon }}_{500})\approx 0.03$ говорит о работоспособности метода и его применимости для восстановления подобных неоднородностей. На рис. 1в линиями 2 и 4 изображены зависимости от номера итерации $m$ величин $\delta ({\widehat{\epsilon }}_m)$ и ${\delta }_{\Gamma }({\widehat{\epsilon }}_m),$ соответственно. Видно, что обе они достаточно быстро стремятся к нулю.

Для проверки устойчивости метода в его входные данные, т.е. матрицы ${\Gamma }_{ij}^{(n)}$ вносились помехи. С этой целью ко всем принятым сигналам $U_j^{(n)}$ добавлялись независимые комплексные случайные величины, действительные и мнимые части которых были распределены нормально с дисперсией, равной $\frac{{\chi }^2}{2N_F{N}_R}{\displaystyle \sum _{j,n}{\left|U_j^{(n)}\right|}^2}$; $\chi ={10}^{-2}$. Зависимости $\delta ({\widehat{\epsilon }}_m)$ и ${\delta }_{\Gamma }({\widehat{\epsilon }}_m)$ для данных с помехами изображены на рис. 1в линиями 1 и 3. Видно, что обе они лежат выше исходных зависимостей; при этом зависимость $\delta ({\widehat{\epsilon }}_m)$ не стремится к нулю с увеличением номера итерации, а зависимость ${\delta }_{\Gamma }({\widehat{\epsilon }}_m)$ – стремится. Все это говорит о том, что введенные сравнительно небольшие помехи, с одной стороны, замедляют скорость сходимости метода, а с другой стороны, приводит к тому, что результатом восстановления является неоднородность, существенно отличающаяся от искомой. Данное обстоятельство является типичным при решении обратных задач в присутствии сильного рассеивателя.

Рис. 1. (а) – Пространственное распределение относительной скорости звука $c(x,y)/c_0$ для неоднородности № 1. (б) – Действительные (линии 1 и 2) и мнимые (линии 3 и 4) части нормированных на $k_0^2$ оценки ${\widehat{\epsilon }}_{500}(x,y)$ (линии 1 и 3) и искомой функции рассеивателя $\epsilon (x,y)$ (линии 2 и 4) вдоль отрезка $AB.$ (в) – Зависимости величин δ(_m) (линии 1 и 2) и δ_Г(_m) (линии 3 и 4) от номера итерации $m$ при точных входных данных (линии 2 и 4) и при наличии помех (линии 1 и 3).

4. Возможности ускорения итерационной процедуры

При численном моделировании область $\Re$, занимаемая рассматриваемой неоднородностью, разбивается на отдельные элементы с размером, равным выбранному шагу дискретизации и характеризующиеся функцией рассеивателя ${\epsilon }_r$, где $r$ в данном случае – индекс такого элемента. Тогда (14) сводится к системе линейных уравнений стандартного вида A_{nij;r r} = B_nij , где $n=\overline{1;N_F}$ и $i,j=\overline{1;N_R}$. Эта система решается с помощью регуляризации Тихонова: $\widehat{\epsilon }={\left(A^{H}A+\kappa E\right)}^{-1}{A}^{H}B$, где верхний индекс “^H” означает гильбертово сопряжение; $E$ – единичная матрица соответствующего размера, и $\kappa$ – коэффициент регуляризации. На практике выбор значения $\kappa$ может быть затруднен. С одной стороны, если этот коэффициент мал, то решение теряет устойчивость, и относительная невязка δ(_m) при восстановлении неоднородности может быстро возрастать. С другой стороны, большое значение $\kappa$ приводит к уменьшению производимой на каждой итерации добавки к восстанавливаемой функции рассеивателя, а значит, скорость сходимости уменьшается. Подобрать $\kappa$ можно, исходя из вида спектра матрицы $A^{H}A.$ На рис. 2а он представлен для неоднородности № 1 сплошной линией. Можно видеть, что начиная с некоторого номера $n_b$ ее собственные числа $\sigma (n)$ начинают резко уменьшаться вплоть до номера $n_{FR}$. Тогда можно положить $\kappa =\sigma (n_b).$ Однако, такой способ оставляет существенную неопределенность из-за резкого поведения функции $\sigma (n)$ вблизи выбранной точки $n_b$, и это сказывается на работе итерационного метода. Например, на рис. 2б тонкими черными линиями представлены зависимости $\delta ({\widehat{\epsilon }}_m)$ при трех различных $\kappa ,$ которые соответствуют отмеченным римскими цифрами точкам на рис. 2а. Можно видеть, что наилучшая сходимость обеспечивается при $\kappa =\sigma (n_b^{\text{II}})$. Однако это значение может зависеть от заранее неизвестных факторов: контраста неоднородности и возможных помех. Наиболее чувствительны к его выбору первые итерации метода, поскольку в начале вычислений ошибка может привести к получению решения, полностью отличного от искомого (рис. 2б при $\kappa =\sigma (n_b^{\text{III}})$ ). С учетом этого вместо фиксирования постоянного коэффициента $\kappa$ предлагается его экспоненциальное уменьшение с ростом числа итераций: ${\kappa }_{m+1}={\gamma }_{\kappa }{\kappa }_m,$ где на первой итерации выбирается заведомо достаточно большое значение ${\kappa }_1$, а значение параметра ${\gamma }_{\kappa }$ устанавливается немного меньшим единицы. Например, на рис. 2б толстой серой линией представлена зависимость $\delta ({\widehat{\epsilon }}_m)$ при $\kappa =\sigma (n_b^{\text{I}})$ и ${\gamma }_{\kappa }=0.98$. Можно видеть, что описанный подход не только снимает вопрос с точным выбором коэффициента регуляризации, но и ускоряет сходимость метода. Ускорение достигается за счет некоторого увеличения добавок к функции рассеивателя при больших номерах $m$.

Важное наблюдение состоит в том, что значение номера $n_{FR}\approx 2N_F{N}_R.$ Здесь произведение $N_F{N}_R$ соответствует числу трасс источник-приемник, а множитель 2 отражает удвоение числа уравнений системы за счет раздельного приравнивания действительных и мнимых частей. Таким образом, хотя решаемая система состоит из $2N_F{N}_R^2$ уравнений, линейно независимыми из них является гораздо меньшее число. Анализ показывает, что коэффициент корреляции строк матрицы $A$, которые соответствуют разным индексами $n$, близок к нулю, т.е. такие строки можно считать независимыми. Наоборот, коэффициент корреляции строк с равными индексами $n$ достаточно большой. Это наводит на мысль, что можно исключить бóльшую часть уравнений, тем самым существенно уменьшив объем вычислений и потребляемой памяти. Например, можно сформировать укороченную систему уравнений с матрицей $A_{\text{short}}$, которая состоит из строк матрицы $A,$ для которых $i=j$. Число строк в такой матрице равно $2N_F{N}_R$. Уравнения, входящие в такую укороченную систему, имеют простой физический смысл: они соответствуют измерению мощности сигнала каждым приемником. Таким образом, использование укороченной системы позволяет также существенно упростить систему обработки сигналов. С другой стороны, на рис. 2а можно видеть, что спектр матрицы $A_{\text{short}}^H{A}_{\text{short}}$ лежит ниже и левее, чем спектр матрицы $A^{H}A$. Это свидетельствует о потере части информации при переходе от $A$ к $A_{\text{short}}$. В итоге на практике восстановление неоднородностей с помощью укороченной системы уравнений происходит хуже, что особенно сильно проявляется при наличии помех в данных.

Рис. 2. (а) – Упорядоченные по убыванию собственные значения $\sigma (n)$ полной матрицы $A^{H}A$ (линия 1) и укороченной (линия 2) матрицы $A_{\text{short}}^H{A}_{\text{short}}$, нормированные на наибольшее собственное значение $\sigma (1)$. Римскими цифрами обозначены точки возможного выбора номера $n_b$. (б) – Зависимости величин $\delta ({\widehat{\epsilon }}_m)$ от номера итерации $m$ при постоянных коэффициентах регуляризации, соответствующих отмеченным римскими цифрами точкам (тонкие черные линии) и при экспоненциально уменьшающемся коэффициенте регуляризации (толстая серая линия).

Другой подход к улучшению сходимости итераций, применявшийся, например, в [22, 23], состоит в том, чтобы осуществлять фильтрацию пространственного спектра $\tilde{\epsilon }(K)$ функции $\epsilon (r)$ в круге $\left|K\right|\le 2\tau k_0$, где величина $\tau \in (0;1]$. Такая фильтрация связана с тем, что томографическому восстановлению подлежат только компоненты, которые соответствуют векторам пространственных частот $K$, не превосходящих по модулю значения $2k_0$ [17, 18]. Компоненты спектра, возникающие после решения (14), но не удовлетворяющие этому требованию, заведомо не являются физическими, и поэтому они должны быть отброшены. Этому соответствует $\tau =1$. В зависимости от размеров и контраста неоднородности можно дополнительно подвергать фильтрации даже компоненты пространственного спектра, для которых $\left|K\right|\le 2k_0,$ т.е. устанавливать $\tau <1$. В [22] применяется ступенчатое изменение $\tau$ в зависимости от номера итерации $m$. Однако такая зависимость $\tau (m)$ довольно сложная, и задача определения ее параметров для разных неоднородностей представляется избыточной. Вместо этого предлагается использовать простую кусочно-линейную функцию $\tau (m)=\max (1;{\tau }_0+{\gamma }_{\tau }(1-{\tau }_0)m)$, где ${\tau }_0$ определяет начальную область пространственного спектра, подлежащую восстановлению; ${\gamma }_{\tau }^{-1}$ должно быть заведомо меньше предполагаемого числа итераций. Описанная процедура приводит к тому, что на первых итерациях восстанавливаются компоненты пространственного спектра $\tilde{\epsilon }(K)$ с низкими пространственными частотами $\left|K\right|$, т.е. крупномасштабные неоднородности, а затем определяются и более мелкие детали.

Хотя описанная процедура фильтрации может показаться в полной мере обоснованной, она неявно предполагает, что пространственный спектр неоднородности сосредоточен в основном вблизи нуля. В случае если это не так, в ходе первых итераций может получиться оценка функции рассеивателя, состоящая из компонент с низкими пространственными частотами и сильно отличающаяся поэтому от искомой. На фоне этой, изначально неверной, оценки в рамках следующих итераций может быть затруднительно обеспечить сходимость, даже при достижении значения $\tau =1$.

Чтобы проиллюстрировать это обстоятельство, рассматриваются еще две неоднородности. Неоднородность № 2 задается функцией $s_2(x,y)\equiv \exp \left(-\frac{x^2+y^2}{2\times {0.75}^2}\right)$, а неоднородность № 3 – функцией $s_3(x,y)\equiv s_2(x,y)\cos (2\pi x)$. Параметры ${\alpha }_0=0$ и ${\beta }_0=0.4$ полагаются для них равными; система источников и приемников поля совпадает с рассмотренной ранее. Поскольку функция рассеивателя связана со скоростью звука нелинейно, то ширина пространственного спектра неоднородностей может превосходить $2k_0$, особенно это касается неоднородности № 3. Поэтому после вычисления $\epsilon (r)$ для каждой неоднородности соответствующий пространственный спектр $\tilde{\epsilon }(K)$ полагается равным нулю при $\left|K\right|>2k_0$. В итоге, восстановлению подлежат уже отфильтрованные неоднородности. Абсолютные значения их пространственных спектров представлены на рис. 3а, 3б. Исходя из значений характеристик для неоднородностей № 2 и № 3, приведенных в таблице, их можно считать в целом несколько более сильными рассеивателями, чем неоднородность № 1.

На рис. 3в представлены зависимости относительных невязок $\delta ({\widehat{\epsilon }}_m)$ для двух режимов изменения $\tau (m)$. В одном случае этот параметр приравнивается единице, а в другом – изменяется по описанной ранее кусочно-линейной зависимости, где полагается ${\tau }_0=0.1$ и ${\gamma }_{\tau }^{-1}=200$. Можно видеть, что для неоднородности № 2, пространственный спектр которой хорошо локализован вблизи нуля, введение фильтрации с $\tau <1$ позволяет добиться существенного ускорения сходимости итераций. С другой стороны, для неоднородности № 3, значительная часть пространственного спектра которой сосредоточена вблизи точек $(\pm k_0;0)$, ситуация является обратной, и введение фильтрации только замедляет сходимость.

Рис. 3. Нормированные на свой максимум пространственные спектры $\tilde{\epsilon }(K)$ двух неоднородностей скорости звука, заданных функциями (а) – $s_2(x,y)$ и (б) – $s_3(x,y)$ и отфильтрованных внутри изображенных окружностей радиуса $2k_0$. (в) – Зависимости величин $\delta ({\widehat{\epsilon }}_m)$ от номера итерации $m$ при восстановлении неоднородности № 2 (линии 1 и 2) и неоднородности № 3 (линии 3 и 4). Линии 2 и 3 соответствуют итерациям с постоянным $\tau =1$. Линии 1 и 4 соответствует кусочно-линейной зависимости $\tau (m)$.

5. Восстановление рассеивателей разного контраста

Если восстановление неоднородности выполняется разными методами на основе одних и тех же экспериментальных данных, значение относительной невязки $\delta (\widehat{\epsilon })$ является одним из критериев, позволяющих выбрать из них наилучший. Таким образом, можно сравнивать методы между собой или определять множества рассеивателей, для которых тот или иной метод показывает лучшие результаты.

Для корреляционного итерационного метода произвести такое сравнение в общем случае оказывается затруднительным. Это связано с различием во входных данных: корреляционный итерационный метод основан на обработке матриц ${\Gamma }_{ij}^{(n)},$ в то время как входными данными других методов служат либо значения принятых полей, либо вычисляемая на их основе амплитуда рассеяния $f(\theta ,{\theta }^{\prime })$. Однако в наиболее простом случае, когда и контролируемые источники, и приемники являются точечными, а фоновые источники отсутствуют (т.е. $I^{(0)}(r)=0$ ), выражение (9) принимает вид ${\Gamma }_{ij}^{(n)}=I^{(n)}G{\left(r_i,r^{(n)}\right)}^{}{G}^{*}\left(r_j,r^{(n)}\right).$ При этом поле, регистрируемое $j$ -м приемником при $n$ -м ракурсе облучения, согласно (2), равно $\phi {\left(r_j,r^{(n)}\right)}^{}=G{\left(r_j,r^{(n)}\right)}^{}\sqrt{I^{(n)}}$. Следовательно, выполнено равенство ${\Gamma }_{ij}^{(n)}=\phi {\left(r_i,r^{(n)}\right)}^{}{\phi }^{*}\left(r_j,r^{(n)}\right).$ Оно позволяет вычислять данные, необходимые для корреляционного итерационного метода, на основе данных, полученных в обычном активном эксперименте, а значит, производить сравнение методов.

Результаты восстановления неоднородности корреляционным итерационным методом сравнивались с результатами, полученными с помощью функционально-аналитического алгоритма Новикова [25, 26]. Этот метод не использует итераций и в настоящее время является одним из лучших для решения обратных задач рассеяния. Математически показано, что при отсутствии рассеяния назад для неоднородностей с нормой амплитуды рассеяния $‖f(\theta ,{\theta }^{\prime })‖<\frac{1}{3\pi }$ он дает точное решение [27]. Однако на практике работоспособность функционально-аналитических алгоритмов сохраняется в ряде случаев, когда эти требования нарушаются, и величина $‖f(\theta ,{\theta }^{\prime })‖$ оказывается в несколько десятков раз больше [22, 28–30]. При этом отмечается, что в классе сильных рассеивателей неоднородности, внутри которых скорость звука в среднем меньше фоновой, восстанавливаются существенно хуже тех, у которых скорость больше фоновой, при равных значениях $‖f(\theta ,{\theta }^{\prime })‖$ [30]. Входными данными для алгоритма Новикова служит амплитуда рассеяния. Она вычисляется на основе значений $\phi {\left(r_j,r^{(n)}\right)}^{}$, которые интерполируются с помощью преобразования Фурье так, что каждый из углов $\theta$ и ${\theta }^{\prime }$ принимает значения от $0$ до $2\pi$ с шагом $2\pi /1024$.

Численное моделирование проводилось для набора неоднородностей, заданных функцией $s_2(x,y)\equiv \exp \left(-\frac{x^2+y^2}{2\times {0.75}^2}\right)$ (т.е. с той же формой, что и у неоднородности № 2); параметр ${\alpha }_0$ полагался равным нулю, а параметр ${\beta }_0$ изменялся в пределах $[-0.6;1]$. Рассчитанные для каждой из неоднородностей значения $‖f(\theta ,{\theta }^{\prime })‖$ и $\Delta \varphi$ представлены на рис. 4а. Восстановление неоднородностей осуществлялось, в одном случае, на основе точных данных. В другом случае к ним добавлялись помехи, аналогично тому, как это было сделано в пункте 3 настоящей работы; при этом было выбрано значение $\chi ={10}^{-3}$. На рис. 4б приведены значения относительных невязок $\delta ({\widehat{\epsilon }}_{200})$, полученных после 200 итераций при обработке данных корреляционным итерационным методом, а на рис. 4в – значения $\delta ({\widehat{\epsilon }}_{\text{Nov}})$, полученные при обработке алгоритмом Новикова.

Рис. 4. (а) – Зависимости нормы амплитуды рассеяния $‖f(\theta ,{\theta }^{\prime })‖$ (линия 1) и максимального набега фаз $\Delta \varphi$ (линия 2) от параметра неоднородности ${\beta }_0$. Зависимости величин (б) – $\delta ({\widehat{\epsilon }}_{200})$ и (в) – $\delta ({\widehat{\epsilon }}_{\text{Nov}})$ от параметра неоднородности ${\beta }_0$ при точных входных данных (линии 1) и при наличии помех (линии 2); линией 3 изображена зависимость $\delta ({\widehat{\epsilon }}_{\text{Nov}})$ при увеличенном до 32 числе ракурсов излучения и приема.

Во-первых, следует обратить внимание, что при равном объеме входных данных алгоритм Новикова обеспечивает приемлемое качество восстановления неоднородности (невязка не превышает нескольких процентов) в гораздо более широком диапазоне значений ${\beta }_0$, по сравнению с итерационным алгоритмом; это является хорошо известным преимуществом функционально-аналитических методов.

Во-вторых, можно видеть, что норма амплитуды рассеяния не позволяет однозначно определить, удастся ли восстановить соответствующую неоднородность. Так, среди неоднородностей, у которых значение $‖f(\theta ,{\theta }^{\prime })‖\approx \frac{20}{3\pi }$, можно выделить и те, которые восстанавливаются хорошо (например, при ${\beta }_0\approx 0.3$ ), и те, которые не восстанавливаются вовсе (например, при ${\beta }_0\approx -0.5$ ) ни одним из методов. С другой стороны, набег фазы $\Delta \varphi$ оказывается довольно хорошо согласованным с возможностью восстановления неоднородности. Так, для корреляционного итерационного метода неоднородности восстанавливаются при $\left|\Delta \varphi \right|<1.3\pi$, а для алгоритма Новикова – по крайней мере, при $\Delta \varphi >-1.6\pi$ (верхнее значение здесь ограничено диапазоном изменения ${\beta }_0\le 1$ ). Следует отметить, что приведенные конкретные численные значения имеют смысл только для описанных условий моделирования, и при анализе других неоднородностей они будут отличаться.

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

В-четвертых, в диапазоне значений ${\beta }_0\in [-0.1;0.4]$ результаты корреляционного итерационного метода оказываются несколько лучше, чем результаты алгоритма Новикова. Это обстоятельство проявляется при малом количестве ракурсов излучения и приема сигналов. Так, если $N_R=N_F=32$, то значение $\delta ({\widehat{\epsilon }}_{\text{Nov}})$ уменьшается в несколько раз (рис. 4в). Однако в эксперименте объем входных данных часто ограничен, и в таких условиях использование итерационного подхода может оказаться предпочтительным. Кроме того, заслуживает внимания возможность комбинации двух методов, когда оценка, полученная с помощью алгоритма Новикова, используется в качестве априорной оценки для итерационного метода. Такой подход нуждается в отдельном исследовании.

В-пятых, наименьшее значение $\delta (\widehat{\epsilon })$ достигается при ${\beta }_0$, отличном от нуля, как можно было бы изначально предположить. Это объясняется тем, что рассматриваемая величина $\delta (\widehat{\epsilon })$ – относительная невязка, и при малых значениях $\left|{\beta }_0\right|$ значение интеграла ${\displaystyle \underset{\Re }{\int }{\left|\epsilon (r)\right|}^{2}dr}$ в знаменателе выражения (15) тоже мало. Поэтому здесь сказывается вклад ошибок, вызванных численными расчетами.

Заключение

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

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

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

Исследование выполнено за счет гранта Российского научного фонда (проект №19-12-00098).

作者简介

K. Dmitriev

编辑信件的主要联系方式.
Email: presentatio@mail.ru
俄罗斯联邦, Moscow

参考

补充文件

附件文件

动作

1. JATS XML

下载

2. Fig. 1. (a) – Spatial distribution of the relative speed of sound for inhomogeneity No. 1. (b) – Real (lines 1 and 2) and imaginary (lines 3 and 4) parts normalized to estimates (lines 1 and 3) and the desired function scatterer (lines 2 and 4) along the segment (c) – Dependences of the values δ(m) (lines 1 and 2) and δГ(m) (lines 3 and 4) on the iteration number with exact input data (lines 2 and 4 ) and in the pres-ence of interference (lines 1 and 3).

下载 (154KB)

索引源数据

3. Fig. 2. (a) – Eigenvalues of the full matrix (line 1) and shortened (line 2) matrix, ordered in descending order, normalized to the largest eigenvalue. Roman numerals indicate possible number selection points. (b) – Dependences of quantities on the iteration number for constant regularization coefficients corre-sponding to the points marked with Roman numerals (thin black lines) and for an exponentially decreas-ing regularization coefficient (thick gray line).

下载 (155KB)

索引源数据

4. Fig. 3. Spatial spectra of two sound speed inhomogeneities normalized to their maximum, specified by functions (a) – and (b) – and filtered inside the depicted circles of radius . (c) – Dependences of quanti-ties on the iteration number when restoring inhomogeneity No. 2 (lines 1 and 2) and inhomogeneity No. 3 (lines 3 and 4). Lines 2 and 3 correspond to iterations with constant . Lines 1 and 4 correspond to a piecewise linear dependence.

下载 (113KB)

索引源数据

5. Fig. 4. (a) – Dependences of the norm of the scattering amplitude (line 1) and the maximum phase shift (line 2) on the inhomogeneity parameter. Dependences of quantities (b) – and (c) – on the heterogeneity parameter with accurate input data (line 1) and in the presence of interference (line 2); Line 3 shows the dependence with the number of emission and reception angles increased to 32.

下载 (93KB)

索引源数据