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

Аннотация

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

Полный текст

ВВЕДЕНИЕ

Пусть имеется однородная непоглощающая среда, в которой скорость звука постоянна и равна c₀. В этой среде, которая будет называться фоновой, имеется область $\mathfrak{R}$ с конечными размерами – так называемая область рассеяния. Область $\mathfrak{R}$ неоднородна по своим акустическим характеристикам – скорости звука $c\left(\overrightarrow{r}\right)$ и амплитудному коэффициенту поглощения $\alpha (\overrightarrow{r},\mathrm{\omega })$. Будет рассматриваться монохроматический случай с частотой ω, однако приводимые ниже соотношения обобщаются на импульсный или многочастотный режим зондирования с помощью фурье-преобразования по частотам. Вне $\mathfrak{R}$ находятся первичные источники $F_0(\overrightarrow{r},\overrightarrow{x})$, которые создают в фоновой среде падающее поле $u_0(\overrightarrow{r},\overrightarrow{x})$; $\overrightarrow{r}$ – произвольная точка пространства. Зависимость от w в аргументах источников и полей будет опускаться в силу монохроматического режима. Радиус-вектор $\overrightarrow{x}$ здесь условно играет роль параметра, задающего падающее поле, которое может быть произвольным. Например, если излучатель квазиточечный, то $\overrightarrow{x}$ – это радиус-вектор центра излучателя. Поле $u_0(\overrightarrow{r},\overrightarrow{x})$ падает на область $\mathfrak{R}$ и рассеивается на неоднородностях скорости звука и коэффициента поглощения. В результате получается поле $u(\overrightarrow{r},\overrightarrow{x})$, удовлетворяющее уравнению [1]:

${\text{∇}}^{2}u(\overrightarrow{r},\overrightarrow{x})+k_0^{2}u(\overrightarrow{r},\overrightarrow{x})=F_0(\overrightarrow{r},\overrightarrow{x})+v\left(\overrightarrow{r}\right)u(\overrightarrow{r},\overrightarrow{x}).$(1)

Здесь $k_0={\omega /c}_0$; $v\left(\overrightarrow{r}\right)$ – функция рассеивателя, которая описывает неоднородности и имеет вид:

$v\left(\overrightarrow{r}\right)={\mathrm{\omega }}^2\left(\frac{1}{c_0^2}-\frac{1}{c^2\left(\overrightarrow{r}\right)}\right)-2i\mathrm{\omega }\frac{\alpha (\overrightarrow{r},\mathrm{\omega })}{c\left(\overrightarrow{r}\right)}$(2)

при временнóй зависимости полей $~\exp (-\mathrm{i\omega t})$. Если в области $\mathfrak{R}$ присутствуют еще неоднородности плотности среды, то в правой части выражения (2) добавляется еще одно слагаемое. Это не изменяет сути последующего изложения. Поскольку вне $\mathfrak{R}$ находится фоновая среда, то $v\left(\overrightarrow{r}\right)\equiv 0$ при $\overrightarrow{r}\notin \mathfrak{R}$.

Рассматривается прямая задача рассеяния, а именно неоднородности в области $\mathfrak{R}$(т.е. рассеиватель), и, следовательно, функция $v\left(\overrightarrow{r}\right)$ считаются известными. Падающее поле $u_0(\overrightarrow{r},\overrightarrow{x})$ тоже известно. Нужно найти поле $u(\overrightarrow{r},\overrightarrow{x})$ и в том числе поле, регистрируемое приемниками. Решение прямой задачи востребовано во многих практических приложениях. Во-первых, это задачи, связанные с отслеживанием структуры поля, распространяющегося через заданную неоднородную среду. Во-вторых, это обратные задачи томографического типа. В них надо восстановить функции, описывающие неизвестные неоднородности, на основе тех или иных характеристик полей, которые измеряются вне области $\mathfrak{R}$ при различных заданных падающих полях u₀ [1]. Чтобы апробировать алгоритмы обработки экспериментальных данных в таких обратных задачах, сначала нужны модельные данные, которые как раз являются решением прямой задачи. Более того, в ряде алгоритмов решения обратной задачи требуется уметь решать прямую задачу в качестве вспомогательной. Однако на практике для области $\mathfrak{R}$ с большим волновым размером при решении прямой задачи могут возникать принципиальные трудности чисто вычислительного характера, если требуется обеспечить очень высокую точность решения. Это приводит к необходимости видоизменения способов решения прямой задачи, и один из них предлагается ниже. В разделе “Заключение” настоящей работы будет дано сравнение предлагаемого способа с решением, которое в настоящее время широко применяется благодаря библиотеке программ “K-Wave”.

РЕШЕНИЕ ПРЯМОЙ ЗАДАЧИ РАЗБИЕНИЕМ НА ПОДОБЛАСТИ

Поле $u(\overrightarrow{r},\overrightarrow{x})$ в произвольной точке $\overrightarrow{r}$ подчиняется уравнению [1]:

$u(\overrightarrow{r},\overrightarrow{x})=u_0(\overrightarrow{r},\overrightarrow{x})+\underset{\Re }{\int }{G}_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})v\left({\overrightarrow{r}}^{\text{'}}\right)u({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}},$(3)

где $G_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})$ – известная функция Грина фоновой среды. При решении прямой задачи совокупность уравнений (3), рассматриваемых для всех точек $\overrightarrow{r}\in \mathfrak{R}$, позволяет при каждом фиксированном падающем поле $u_0(\overrightarrow{r},\overrightarrow{x})$ найти поле $u(\overrightarrow{r},\overrightarrow{x})$ внутри области $\mathfrak{R}$. С другой стороны, прямую задачу можно сначала решать в терминах так называемых вторичных источников $I(\overrightarrow{r},\overrightarrow{x})\equiv v\left(\overrightarrow{r}\right)\text{\hspace{0.33em}}u(\overrightarrow{r},\overrightarrow{x})$, которые возникают в области $\mathfrak{R}$ в ответ на поле в каждой точке $\overrightarrow{r}\in \mathfrak{R}$. Уравнение для них получается умножением соотношения (3) на $v\left(\overrightarrow{r}\right)$:

$I(\overrightarrow{r},\overrightarrow{x})=v\left(\overrightarrow{r}\right)\text{\hspace{0.33em}}{u}_0(\overrightarrow{r},\overrightarrow{x})+v\left(\overrightarrow{r}\right)\underset{\Re }{\int }{G}_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})I({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}},$

$\overrightarrow{r}\in \mathfrak{R}$. (4)

Оба варианта решения – в форме (3) или в форме (4) – равноправны и имеют единственное решение для пассивных сред [1]. Поскольку функция рассеивателя $v\left(\overrightarrow{r}\right)$ известна, задача линейна относительно искомых функций $u(\overrightarrow{r},\overrightarrow{x})$ или $I(\overrightarrow{r},\overrightarrow{x})$, которые в обоих вариантах ищутся при $\overrightarrow{r}\in \mathfrak{R}$.

После нахождения $u(\overrightarrow{r},\overrightarrow{x})$ или $I(\overrightarrow{r},\overrightarrow{x})$ уравнение (3) позволяет непосредственно рассчитать поля $u(\overrightarrow{r},\overrightarrow{x})$ уже в произвольной точке $\overrightarrow{r}$ – как внутри области $\mathfrak{R}$, так и вне нее. В частности, могут быть рассчитаны поля в точках приема, которые характеризуются радиус-вектором $\overrightarrow{y}$:

$u(\overrightarrow{y},\overrightarrow{x})=u_0(\overrightarrow{y},\overrightarrow{x})+\underset{\mathfrak{R}}{\int }{G}_0(\overrightarrow{y},{\overrightarrow{r}}^{\text{'}})I({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}}.$(5)

Для определенности ниже обсуждается процедура нахождения вторичных источников $I(\overrightarrow{r},\overrightarrow{x})$ на основе уравнения (4), которое можно переписать в виде:

$\begin{array}{c}\underset{\mathfrak{R}}{\int }\left\{\delta (\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}})-v\left(\overrightarrow{r}\right)G_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})\right\}\times \\ \times I({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}}=v\left(\overrightarrow{r}\right)\text{\hspace{0.33em}}{u}_0(\overrightarrow{r},\overrightarrow{x}),\end{array}$(6)

где $\delta (\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}})$ – дельта-функция Дирака. Дискретизованная система уравнений (6) решается при каждом фиксированном $\overrightarrow{x}$. Однако в ряде прикладных задач область $\mathfrak{R}$ с линейным размером L имеет большой волновой размер $k_{0}L>>1$, т.е. L составляет десятки длин волн и более. Тогда при решении системы (6) возникают принципиальные трудности вычислительного характера, которые заключаются в следующем. Пусть искомые вторичные источники $I(\overrightarrow{r},\overrightarrow{x})$ при фиксированном $\overrightarrow{x}$ описываются дискретными значениями $I({\overrightarrow{r}}_n,\overrightarrow{x})$, которые соответствуют N дискретным пространственным отсчетам ${\overrightarrow{r}}_n\in \Re$ в пространстве любой размерности. Тогда N искомых неизвестных $I({\overrightarrow{r}}_n,\overrightarrow{x})$ формируются в виде столбца при решении линейной дискретизованной системы (6). В то же время матрица при неизвестных соответствует в (6) интегральному оператору с ядром $\left\{\delta (\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}})-v\left(\overrightarrow{r}\right)G_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})\right\}$ и описывается дискретизованными аргументами $(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})$. Эта матрица будет иметь формат $N\times N$. Например, патологии молочной железы диагностируются в акустических томографах с помощью ультразвука в наиболее часто используемом диапазоне частот 0.5–2 МГц [1—4]. Линейный размер L области $\mathfrak{R}$, соответствующей молочной железе, составляет 0.1–0.2 м, т.е. (50–200)λ₀,  где ${\lambda }_0=2\pi /k{}_0$ – длина волны в фоновой среде. В целях ранней диагностики необходимо учитывать тонкую структуру рассеивателя – детали с линейным размером (0.2–0.5)λ₀, и тогда вспомогательную прямую задачу требуется решать с очень высокой точностью (погрешность менее 1%). Как следствие, шаг пространственной дискретизации не должен превышать (0.1–0.2)λ₀, и количество дискретных отсчетов по одному измерению может достигать, по порядку величины, $\cong {10}^3$. Тогда общее количество двумерных неизвестных составляет $N\cong {10}^6$, и формат матрицы при неизвестных становится порядка ${10}^{\text{6}}\times {10}^6$. Выделение места для такой матрицы вызывает переполнение оперативной памяти даже очень мощной персональной ЭВМ, и решить систему (6) непосредственно становится невозможным. В то же время использование вычислительных кластеров нежелательно, поскольку вычислительные средства должны быть доступны любому медицинскому учреждению.

В связи с обозначенной проблемой предлагается итерационный метод решения прямой задачи рассеяния. Исследуемая область $\mathfrak{R}$ разбивается на совокупность подобластей ${\mathfrak{R}}^m$, решение прямой задачи находится в каждой отдельной подобласти и потом итерационно уточняется. Здесь $m=\mathrm{1,}\mathrm{2,}\mathrm{3,}\mathrm{...},M$ – номер подобласти. Объединение всех подобластей ${\mathfrak{R}}^m$ дает полную область $\mathfrak{R}$. Подобная идея уже использовалась ранее для нахождения термоакустических источников [5, 6], и ниже предлагаемый метод рассматривается подробнее.

Чтобы начать итерационный процесс (j – номер итерации) для решения прямой задачи рассеяния, сначала (j = 0) задается начальное приближение ${\widehat{I}}_{(j=0)}(\overrightarrow{r},\overrightarrow{x})$ для оценки вторичных источников $I(\overrightarrow{r},\overrightarrow{x})$. В простейшем случае можно использовать оценку вторичных источников в приближении Борна, когда неизвестное поле $u\text{(}\overrightarrow{r}\text{,}\overrightarrow{x})$ внутри области $\mathfrak{R}$ заменяются на известное падающее поле $u_0\text{(}\overrightarrow{r}\text{,}\overrightarrow{x})$:${\widehat{I}}_{\left(0\right)}(\overrightarrow{r},\overrightarrow{x})=v\left(\overrightarrow{r}\right)u_0\text{(}\overrightarrow{r}\text{,}\overrightarrow{x}\text{).}$ Однако на практике наиболее часто встречаются рассеиватели, которые сильно искажают падающее на них поле. Для таких рассеивателей начальное приближение в виде борновской оценки не будет гарантировать сходимость итерационной процедуры, и целесообразно использовать более точное начальное приближение. Например, при $j=0$ вторичные источники $(\overrightarrow{r},\overrightarrow{x})$ дискретизуются по $\overrightarrow{r}\in \mathfrak{R}$ с таким крупным (грубым) шагом, чтобы матрица при неизвестных системы (6) не вызывала переполнения оперативной памяти ЭВМ. Тогда решение системы (6) сразу во всей области $\mathfrak{R}$ даст начальную оценку ${\widehat{I}}_{\left(0\right)}(\overrightarrow{r},\overrightarrow{x}).$

Далее задается и фиксируется более мелкий шаг пространственной дискретизации. Оценка ${\widehat{I}}_{\left(0\right)}(\overrightarrow{r},\overrightarrow{x})$ предварительно интерполируется на сетку пространственных отсчетов, которая соответствует выбранному мелкому шагу. При таком шаге прямая задача уже не может быть решена сразу для всей области $\mathfrak{R}$, и рассматриваются подобласти ${\mathfrak{R}}^m$. Начинается итерационный цикл ($j\ge 1$), в котором заданный шаг дискретизации остается неизменным. На каждой итерации с фиксированным номером j поочередно рассматриваются все подобласти ${\mathfrak{R}}^m$, $m=\mathrm{1,}\mathrm{2,}\mathrm{3,}\mathrm{...},M$. Для каждой фиксированной подобласти ${\mathfrak{R}}^m$ интеграл по области $\mathfrak{R}$ в (4) представляется в виде суммы интегралов по выбранной подобласти ${\mathfrak{R}}^m$ и оставшейся части $\mathfrak{R}\backslash {\mathfrak{R}}^m$. Тогда система (4) переписывается в виде:

$\begin{array}{c}I(\overrightarrow{r},\overrightarrow{x})-v\left(\overrightarrow{r}\right)\underset{{\mathfrak{R}}^m}{\int }{G}_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})I({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}}=\\ =v\left(\overrightarrow{r}\right)u_0(\overrightarrow{r},\overrightarrow{x})+v\left(\overrightarrow{r}\right)\underset{\mathfrak{R}\backslash {\mathfrak{R}}^m}{\int }{G}_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})I({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}},\end{array}$

$\overrightarrow{r}\in {\mathfrak{R}}^m$.(7)

Система (7) дискретизуется и решается относительно неизвестных $I(\overrightarrow{r},\overrightarrow{x})$ только в данной подобласти $\overrightarrow{r}\in {\mathfrak{R}}^m$. Эти неизвестные стоят в левой части (7). При этом правая часть системы (7) должна быть известной. С этой целью в правой части (7) вместо неизвестных точных $I({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})$ значений используется их текущая оценка ${\widehat{I}}_{j|m}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})$, которая формируется из уже ранее полученных оценок вторичных источников $\widehat{I}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})$ при ${\overrightarrow{r}}^{\text{'}}\notin {\mathfrak{R}}^m$. Конкретное выражение для ${\widehat{I}}_{j|m}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})$ приводится ниже в (9).

Таким образом, для фиксированной итерации с номером j и текущей подобласти с номером m вместо строгих уравнений (7) используется следующая система уравнений:

$\begin{array}{c}{I}_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})-v\left(\overrightarrow{r}\right)\underset{{\Re }^m}{\int }{G}_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})I_{\left(j\right)}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}}=\\ =v\left(\overrightarrow{r}\right)u_0(\overrightarrow{r},\overrightarrow{x})+v\left(\overrightarrow{r}\right)\underset{\Re \backslash {\Re }^m}{\int }{G}_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}}){\widehat{I}}_{j|m}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}},\\ \overrightarrow{r}\in {\Re }^m.\end{array}$(8)

Здесь оценка ${\widehat{I}}_{j|m}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})$, ${\overrightarrow{r}}^{\text{'}}\notin {\mathfrak{R}}^m$, для фиксированной итерации с номером j и текущей подобласти с номером m строится по принципу:

${\widehat{I}}_{j|m}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})=\left\{\begin{array}{l}{\widehat{I}}_{\left(j\right)}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})\text{при}{\overrightarrow{r}}^{\text{'}}\in {\Re }^k,k=\mathrm{1,}\mathrm{2,}\mathrm{...},(m-1);\\ {\widehat{I}}_{(j-1)}({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})\text{при}{\overrightarrow{r}}^{\text{'}}\in {\Re }^k,k=(m+1),(m+2),\mathrm{...},M.\end{array}\right.$(9)

Система (8) решается при каждом фиксированном $\overrightarrow{x}$ относительно дискретизованных значений $I_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})$ в подобласти $\overrightarrow{r}\in {\mathfrak{R}}^m$. Найденное решение ${\widehat{I}}_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})$ используется далее в качестве уточненной оценки вторичных источников в данной подобласти $\overrightarrow{r}\in {\mathfrak{R}}^m$. При переходе к новой подобласти решается своя система (8), которая формируется для новых точек $\overrightarrow{r}$ при участии уточненной оценки (9). После решения системы (8) для всех подобластей $m=\mathrm{1,}\mathrm{2,}\mathrm{3,}\mathrm{...},M$ получается оценка ${\widehat{I}}_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})$, итоговая для данной итерации j, уже в полной области $\overrightarrow{r}\in \mathfrak{R}$.

Начинается новая итерация $j+1$. Шаг пространственной дискретизации не изменяется. Все подобласти перебираются заново, и предыдущая оценка ${\widehat{I}}_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})$ в каждой подобласти заменяется на уточненную оценку ${\widehat{I}}_{(j+1)}(\overrightarrow{r},\overrightarrow{x}).$ Итерации при фиксированном шаге дискретизации продолжаются до их сходимости, и итерационный цикл завершается.

Если шаг пространственной дискретизации, использованный в итерационном цикле, был крупнее, чем ожидаемая разрешающая способность для итоговых оценок $\widehat{I}(\overrightarrow{r},\overrightarrow{x})$ искомых функций $I(\overrightarrow{r},\overrightarrow{x})$, то шаг дискретизации опять уменьшается. С таким фиксированным шагом выполняется новый итерационный цикл. Начальным приближением для этого цикла служит оценка вторичных источников $\widehat{I}(\overrightarrow{r},\overrightarrow{x})$, полученная в итоге предыдущего итерационного цикла. Как уже упоминалось, прямая задача имеет единственное решение, однако требуется обеспечить сходимость итераций. Поэтому прием постепенного уменьшения шага дискретизации целесообразно использовать как раз для обеспечения сходимости итерационного процесса.

Сходимость итераций в каждом итерационном цикле может контролироваться с помощью следующего критерия. После окончания каждой фиксированной итерации j, т.е. после перебора всех подобластей, считается невязка ${\mathrm{\epsilon }}_j\left(\overrightarrow{x}\right)$ между текущей оценкой ${\widehat{I}}_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})$ и оценкой ${\widehat{I}}_{(j-1)}(\overrightarrow{r},\overrightarrow{x})$, полученной на предыдущей (или еще более ранней) итерации ${\widehat{I}}_{(j-1)}(\overrightarrow{r},\overrightarrow{x})$:

${\epsilon }_j\left(\overrightarrow{x}\right)=\sqrt{\frac{{\displaystyle \underset{\Re }{\int }{\left|{\widehat{I}}_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})-{\widehat{I}}_{(j-1)}(\overrightarrow{r},\overrightarrow{x})\right|}^{2}d\overrightarrow{r}}}{{\displaystyle \underset{\Re }{\int }{\left|{\widehat{I}}_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})\right|}^{2}d\overrightarrow{r}}}}.$(10)

Итерации прекращаются, когда невязка становится меньше наперед заданного числа ε: ${\mathrm{\epsilon }}_j\left(\overrightarrow{x}\right)<\mathrm{\epsilon }$для всех $\overrightarrow{x}$, т.е. для всех падающих полей. Конкретное значение ε определяется точностью, с которой нужно решить прямую задачу при выбранном шаге дискретизации.

Найденная итоговая оценка вторичных источников $I(\overrightarrow{r},\overrightarrow{x})$ позволяет найти поле $u(\overrightarrow{r},\overrightarrow{x})$ в произвольной точке $\overrightarrow{r}$, согласно (3):

$u(\overrightarrow{r},\overrightarrow{x})=u_0(\overrightarrow{r},\overrightarrow{x})+\underset{\Re }{\int }{G}_0(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}})I({\overrightarrow{r}}^{\text{'}},\overrightarrow{x})d{\overrightarrow{r}}^{\text{'}}$

и, в частности, поле на приемниках (5).

АПРОБАЦИЯ ЧИСЛЕННЫМ МОДЕЛИРОВАНИЕМ

Численное моделирование итерационного метода было выполнено с помощью рефракционно-поглощающих неоднородностей, описываемых функцией рассеивателя (2). Для определенности рассматривался двумерный случай, хотя сам алгоритм может применяться в пространстве произвольной размерности. Задаваемые параметры были близки к параметрам ультразвукового томографа [1, 3, 4]. Квазиточечные приемоизлучающие преобразователи располагались равномерно на окружности радиуса 0.1536 м. При частоте 1.25 МГц и с₀ = 1500 м/с длина волны в фоновой среде составляла λ₀ = 1.2 ⋅ 10^–3 м.

Сила рассеивателя оценивается, во-первых, с помощью дополнительного (по сравнению с фоновой средой) набега фазы ΔΨ. Он порождается неоднородностями скорости звука при прохождении волны вдоль какой-либо заданной траектории $l_{\mathfrak{R}}$: $\mathrm{\Delta \psi }\text{\hspace{0.17em}\hspace{0.17em}}\equiv \underset{l_{\mathfrak{R}}}{\int }\left(\frac{\omega }{c_0}-\frac{\omega }{c\left(\overrightarrow{r}\right)}\right)d{l}_{\overrightarrow{r}}=\underset{l_{\mathfrak{R}}}{\int }\frac{\Delta c\left(\overrightarrow{r}\right)/c_0}{1+\Delta c\left(\overrightarrow{r}\right)/c_0}d{l}_{\overrightarrow{r}}$. Здесь $d{l}_{\overrightarrow{r}}$ – длина элемента траектории в окрестности точки $\overrightarrow{r}$; $\mathrm{\Delta }c\left(\overrightarrow{r}\right)=c\left(\overrightarrow{r}\right)-c_0$ – контраст скорости звука. Во-вторых, сила рассеивателя зависит от поглощения в области $\mathfrak{R}$. Так, суммарный коэффициент поглощения вдоль траектории $l_{\mathfrak{R}}$ оценивается как $\underset{l_{\Re }}{\int }\alpha (\overrightarrow{r},\omega )d{l}_{\overrightarrow{r}}$; следовательно, амплитуда волны за счет поглощения уменьшается в $\exp \left(\underset{l_{\Re }}{\int }\alpha (\overrightarrow{r},\omega )d{l}_{\overrightarrow{r}}\right)$ раз.

Для модели, представленной на рис. 1, рассматривалась область $\mathfrak{R}$ с размером $64{\lambda }_0\times 64{\lambda }_0$. Дополнительный набег фазы волны при прохождении вдоль оси OX (сечение $y=0$, рис. 3а) составляет: $\mathrm{\Delta \psi }\approx 1.4\pi$ на участке с положительным контрастом $\mathrm{\Delta }c\left(\overrightarrow{r}\right)>0$, где $\mathrm{Re}v\left(\overrightarrow{r}\right)>0$; $\mathrm{\Delta \psi }\approx -1.1\mathrm{\pi }$ – на участке с отрицательным контрастом $\mathrm{\Delta }c\left(\overrightarrow{r}\right)<0$, где $\mathrm{Re}v\left(\overrightarrow{r}\right)<0$. В то же время при прохождении волны вдоль разных сечений, параллельных оси OY, наибольший (по сечениям) дополнительный набег фазы составляет $\mathrm{\Delta \psi }\approx 2.3\mathrm{\pi }$ на участке с положительным контрастом и $\mathrm{\Delta \psi }\approx -0.9\mathrm{\pi }$ на участке с отрицательным контрастом. Максимальное суммарное поглощение, приводящее к уменьшению амплитуды волны в ≈ 4 раза, соответствует прохождению волны вдоль оси OY (сечение х = 0, рис. 3б). Тем самым рассеиватель достаточно сильно искажает падающие на него поля, т.е. эффекты перерассеяния волн сильны.

Рис. 1. Модель рефракционно-поглощающей неоднородности: общий вид действительной части (а) и мнимой части (б) функции рассеивателя.

Прямая задача решалась для 512 приемоизлучающих преобразователей. Для задания начального приближения ${\widehat{I}}_{\left(0\right)}(\overrightarrow{r},\overrightarrow{x})$ оценки вторичных источников задавался укрупненный шаг дискретизации ${\lambda }_0/2$ вдоль каждой декартовой оси вместо итогового мелкого шага ${\lambda }_0/8$. Тогда количество дискретизованных неизвестных для выбранной модели составляло $N=2^{14}$. Это позволило решить систему (4) сразу во всей области $\mathfrak{R}$. Надо отметить, что переполнение оперативной памяти ЭВМ вызывают матрицы при неизвестных, когда количество неизвестных становится $N>2^{16}$. Абсолютные значения полученной оценки ${\widehat{I}}_{\left(0\right)}(\overrightarrow{r},\overrightarrow{x})$ изображены на рис. 2а для излучателя, положение которого характеризуется углом ${\mathrm{\phi }}_{\overrightarrow{\mathrm{x}}}=1.5\mathrm{\pi }$ в полярной системе координат.

При задании конечного шага дискретизации ${\lambda }_0/8$ вся область $\mathfrak{R}$ разбивалась на 8² подобластей, каждая размером $8{\lambda }_0\times 8{\lambda }_0$. Использовался итерационный алгоритм (8), (9). Относительная итоговая точность для выполнения критерия сходимости (10) задавалась высокой: $\mathrm{\epsilon }={10}^{-4}$. Оценка ${\widehat{I}}_{\left(j\right)}(\overrightarrow{r},\overrightarrow{x})$ после первой итерации (рис. 2б) содержит резкие перепады значений на границах подобластей. В итоговой оценке (после 15 итераций) такие перепады отсутствуют (рис. 2в). Принимаемые поля рассчитывались, согласно (5), для всех 512 угловых положений как излучателей, так и приемников.

Рис. 2. Абсолютные значения оценки вторичных источников: грубая оценка с пространственным разрешением λ₀ / 2 (а) и оценка с высоким разрешением λ₀ / 8 после первой итерации (б) и итоговой пятнадцатой итерации (в).

Правильность решения прямой задачи контролировалась следующим образом. На основе всей совокупности рассчитанных принятых полей решалась уже обратная задача, в которой восстанавливалась оценка $\widehat{v}\left(\overrightarrow{r}\right)$ функции рассеивателя $v\left(\overrightarrow{r}\right)$. Для решения использовался функциональный алгоритм [7–9], который строго учитывает эффекты многократного рассеяния волн на неоднородностях среды. Для уменьшения погрешности восстановления, связанной с эффектами дискретизации, значения полей от 512 приемоизлучающих преобразователей предварительно интерполировались на значения для 2048 преобразователей. Полученная оценка $\widehat{v}\left(\overrightarrow{r}\right)$ практически совпала с истинной функцией $v\left(\overrightarrow{r}\right)$ (рис. 3). Это свидетельствует о работоспособности и высокой точности предложенного итерационного алгоритма решения прямой задачи.

Рис. 3. Контрольное восстановление рассеивателя: сечение y = 0 действительной части (а) и сечение x = 0 мнимой части (б) для истинного рассеивателя v (тонкая линия) и для его восстановленной оценки $\widehat{\mathbf{v}}$ (толстая пунктирная линия).

ЗАКЛЮЧЕНИЕ

В настоящее время широко распространены методы решения прямой задачи рассеяния на основе конечно-разностных схем. Например, один из наиболее популярных пакетов программ “K-Wave” (среда программирования Matlab) использует методы фурье-преобразования по пространственным координатам и введение специального фактора коррекции, что существенно ускоряет вычисления [10–12]. Программы “K-Wave” апробированы в различных прикладных задачах [2, 13]. В то же время при выборе того или иного метода решения прямой задачи надо учитывать прежде всего требуемую точность решения и общее время счета. Так, программы “K-Wave” позволяют учитывать присутствующие неоднородности скорости звука и коэффициента поглощения только в виде значений, определенным образом усредненных. Это приводит к ограничениям на силу неоднородностей и на точность решения [10, 11]. Во многих практических задачах, в том числе в целях терапии или хирургии, точность 5 - 10 % является вполне достаточной [13]. При этом прямую задачу нужно решать для единственного излучателя, который состоит из нескольких сотен элементов, излучающих, однако, одновременно [13]. С другой стороны, в медицинских задачах томографического типа исследуемый объект даже в двумерном случае необходимо облучать поочередно 200 - 300 с направлений как минимум [1–4]. Тогда прямую задачу (которая является вспомогательной при томографии) нужно решить программами “K-Wave” также 200 - 300 раз, что существенно увеличивает общее время счета. В то же время вышеизложенный итерационный метод позволяет получать решения сразу для всех излучателей.

Что касается точности решения, то итерационный метод, в отличие от программ “K-Wave”, остается строгим вне зависимости от контраста и размера неоднородностей. Он позволяет максимально точно учитывать все виды неоднородностей и их тонкую структуру, включая неоднородности показателя степени частотной зависимости коэффициента поглощения [1, 14] (в программах “K-Wave” показатель степени постоянен). Тем самым использование итерационного метода особенно целесообразно, когда требуется получить максимально точное решение. При этом представляется перспективным сочетание программ “K-Wave” с последующим дополнительным уточнением решения итерационным методом за очень малое количество итераций.