Доклады Российской академии наук. Математика, информатика, процессы управления, 2022, T. 505, № 1, стр. 50-55

Задачи акустической томографии были поставлены довольно давно (см., например, работы [1–4]). Исходной мотивацией для них явилось использование звуковых волн для возможной ранней диагностики рака молочной железы. В последнее время, в связи с распространением инфекции COVID-19, методы акустической томографии предложено использовать также для диагностики заболевания легких. Вычислительное моделирование задач акустической томографии выполнено в работах [5–9] (см. также многочисленные цитирования в них). Ниже мы рассматриваем один из возможных вариантов постановки задачи акустической томографии. В рассматриваемой задаче используются точечные импульсные источники возбуждения акустических сигналов расположенные вне той области, в которой разыскиваются переменные коэффициенты уравнения акустики, приемники акустических сигналов, размещены также вне этой области, а временные импульсы измеряются в небольшой окрестности момента прихода сигнала от источника в соответствующий приемник. Поставленная задача относится к классу лучевых обратных задач, введенных в монографии автора [10]. Характерной особенностью таких задач является распадение исходной задачи об определении нескольких коэффициентов на ряд последовательно решаемых задач об определении одного из искомых коэффициентов. В данном случае такими искомыми коэффициентами являются скорость звука c(x), акустическое поглощение $\sigma (x)$ и плотность среды $\rho (x)$. При этом определение c(x) сводится к решению хорошо известной обратной кинематической задачи, а определение $\sigma (x)$ и $\rho (x)$ приводится к решению некоторой задачи интегральной геометрии на семействе геодезических линий конформной римановой метрики, определяемой коэффициентом c(x). В случае, когда $c(x) \equiv 1$, эта задача является обычной задачей рентгеновской томографии.

Рассмотрим задачу Коши:

в которой $p\, = \,p(x,t)$ – звуковое давление, $c(x)\, > \,0$ – скорость звука, $\sigma (x) \geqslant 0$ – акустическое поглощение, $\rho (x) > 0$ – плотность среды. Уравнение (1) описывает распространение акустических волн в неоднородной поглощающей среде. Это уравнение лежит в основе акустической томографии. Основная проблема акустической томографии состоит в решении обратной задачи для уравнения (1), а именно, построении коэффициентов $c(x)$, $\sigma (x)$ и $\rho (x)$ внутри ограниченной области по заданной вне этой области информации о решениях уравнения (1) на некотором временном интервале [0, T] и для некоторого множества точечных источников $y$. Ниже, для определенности, рассмотрим следующую модель.

Пусть ${{B}_{R}}$ – шар радиуса R с центром в начале координат, ${{S}_{R}}$ – его граница. Примем, что параметры $c(x),\;\sigma (x),\;\rho (x)$ являются непрерывно дифференцируемыми функциями всюду в ${{\mathbb{R}}^{3}}$ (см. ниже дополнительные условия (9)), неизвестны внутри шара ${{B}_{{{{r}_{0}}}}}$, $R > {{r}_{0}} > 0$, а вне его постоянны и заданы,

Обозначим через ${{S}_{R}}(y,{{r}_{0}})$ проекцию шара ${{B}_{{{{r}_{0}}}}}$ на сферу ${{S}_{R}}$ от источника, расположенного в точке $y \in {{S}_{R}}$, т.е. ${{S}_{R}}(y,{{r}_{0}}) = \{ x \in {{S}_{R}}\,{\text{|}}\,{\text{|}}x - y{\text{|}} \geqslant 2\sqrt {{{R}^{2}} - r_{0}^{2}} \} $. Решение задачи (1) обозначим через $p(x,t,y)$, подчеркнув тем самым его зависимость от параметра $y$. Зафиксируем точки $y \in {{S}_{R}}$ и $x \in {{S}_{R}}(y,{{r}_{0}})$. Пусть T₀(x, y) := $\sup \{ \tau > 0\,{\text{|}}\,p(x,t,y) \equiv 0,t \in (0,\tau )\} $.

Рассмотрим следующую постановку задачи акустической томографии (АК).

Задача АК. Пусть $p(x,t,y)$ – решение задачи (1) и известна функция

где $\varepsilon > 0$ фиксировано и произвольно мало. По функции $F(x,t,y)$ требуется найти $c(x)$, $\sigma (x)$, $\rho (x)$ внутри ${{B}_{{{{r}_{0}}}}}$.

Введением новой функции u(x, t, y) = = $p(x,t,y){{\rho }^{{1/2}}}(x)$ задача (1) при $y \in {{S}_{R}}$ преобразуется к виду

в которой коэффициент q(x) вычисляется по формуле

Если коэффициент q(x) известен, то, как следует из (5), функция $m(x) = {{\rho }^{{1/2}}}(x)$ может быть найдена внутри ${{B}_{{{{r}_{0}}}}}$ как решение следующей задачи Дирихле

в которой $r = {\text{|}}x{\text{|}}$. Более того, из условия (2) на функцию $\rho (x)$ и условия ее непрерывной дифференцируемости всюду в ${{\mathbb{R}}^{3}}$, следует, что решение задачи (6) должно удовлетворять требованию $\frac{{\partial m}}{{\partial r}}$ = 0 при $r = {{r}_{0}}$. Это дополнительное условие гарантирует единственность определения функции $m(x)$, следовательно, и $\rho (x)$, по заданной $q(x)$. Поэтому вместо задачи об определении $c(x)$, $\sigma (x)$, $\rho (x)$ в ${{B}_{{{{r}_{0}}}}}$ удобно рассматривать обратную задачу о нахождении в той же области коэффициентов $c(x)$, $\sigma (x)$ и $q(x)$ для уравнения (4) по информации, аналогичной (3):

Здесь $\Phi (x,t,y) = F(x,t,y)\rho _{0}^{{1/2}}$, а $\varepsilon $ фиксировано и произвольно мало.

Подобные задачи рассмотрены к книге [10] для $x \in {{\mathbb{R}}^{2}}$ в следующих случаях:

1) $c(x) \equiv 1$, $\sigma (x)$ и q(x) – неизвестны (раздел 5.2),

2) $\sigma (x) \equiv 0$, $c(x)$ и q(x) – неизвестны (раздел 5.3).

Кроме того, в [10] (в разделе 5.4) рассмотрена двумерная обратная задача электродинамики. В случае, когда свойства среды не зависят от координаты ${{x}_{3}}$, уравнение для компоненты ${{E}_{3}}$ вектора электрической напряженности совпадает с уравнением (1), в котором роль функции $\rho (x)$ играет магнитная проницаемость среды, а коэффициент $\sigma (x)$ определяет электрическую проводимость среды. В этой обратной задаче определяются все три неизвестных коэффициента по динамической информации о решениях трех задач Коши, в которых вместо точечного источника $\delta (x - y)$ используется источник типа плоской волны $\delta (x \cdot {{\nu }^{{(k)}}})$, в котором ${{\nu }^{{(k)}}}$ – единичный вектор, и $k = 1,2,3$.

Рассмотрим риманову метрику, в которой элемент длины $d\tau $ определяется формулой dτ = = ${\text{|}}dx{\text{|/}}c(x)$, ${\text{|}}dx{\text{|}} = {{\left( {\sum\limits_{k = 1}^3 {{{{(d{{x}_{k}})}}^{2}}} } \right)}^{{1/2}}}$. Всюду в дальнейшем мы полагаем, что выполнено следующее

Предположение. Риманова метрика dτ = = ${\text{|}}dx{\text{|/}}c(x)$ является простой в ${{\mathbb{R}}^{3}}$, т.е. любые две точки x и y в пространстве ${{\mathbb{R}}^{3}}$ могут быть соединены единственной геодезической линией $\Gamma (x,y)$.

Заметим, что достаточным условием простоты комформной римановой метрики в ${{\mathbb{R}}^{3}}$ является условие (см. [11]):

в котором $\xi = ({{\xi }_{1}},{{\xi }_{2}},{{\xi }_{3}})$.

Замечание 1. На самом деле, в рамках лучевой постановки обратной задачи, сделанное выше предположение можно ослабить, а именно, требовать простоту римановой метрики нужно только внутри фиксированного шара ${{B}_{{R'}}}$ при $R' > R$. Мы ограничились более жестким предположением, чтобы не делать лишних оговорок в дальнейшем.

Предположим, что существуют конечные числа ${{c}_{0}}$ и ${{c}_{{00}}}$ такие, что

и коэффициенты уравнения (4) обладают следующей гладкостью

Обозначим через $\tau (x,y)$ риманово расстояние между точками $x$ и $y$. С физической точки зрения, $\tau (x,y)$ – время пробега акустического сигнала между точками $x$ и $y$. Известно, что $\tau (x,y)$ является решением следующей задачи Коши для уравнения эйконала:

Чтобы найти $\tau (x,y)$ и уравнения геодезических линий, надо поступить следующим образом. Ввести вектор $\eta = {{\nabla }_{x}}\tau (x,y)$, решить задачу Коши для системы обыкновенных дифференциальных уравнений Эйлера

для всевозможных значений единичного вектора $\nu = (\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta )$. Это можно сделать, и притом однозначно, в силу предположения о простоте римановой метрики и условий, наложенных на $c(x)$. В результате решения задачи (11), (12) определяются уравнение геодезических линий $\xi = \xi (s,\nu ,y)$, а также касательный вектор $\eta = \eta (s,\nu ,y)$ к этой геодезической, и риманово расстояние $\tau = \tau (s,\nu ,y)$ вдоль геодезической. Решая уравнение $x = \xi (s,\nu ,y)$ относительно $s,\;\theta ,\;\varphi $, находим $s = s(x,y)$, $\nu = \nu (x,y)$, т.е. соответствия между парой точек $x$ и $y$ и значениями параметров $s$ и $\nu $. При этом уравнение геодезической $\Gamma (x,y)$ задается равенством $\xi = \xi (s,\nu (x,y),y)$, $s \in [0,s(x,y)]$, а риманово расстояние $\tau (x,y)$, в силу соотношений (11), (12) для этой функции, совпадает с $s(x,y)$. Заметим, что, согласно [10] (см. стр. 35), и сделанному выше предположению о гладкости $c(x)$, функция ${{\tau }^{2}} \in {{C}^{6}}({{\mathbb{R}}^{6}})$.

Для $T > 0$ обозначим через ${{D}_{T}}$ полосу

Теорема 1. Пусть выполнено предположение о простоте римановой метрики и условия (8), (9). Тогда решение задачи (4) представимо в области ${{D}_{T}}$ виде

в котором $H(t)$ – функция Хевисайда: $H(t) = 1$ для $t \geqslant 0$ и $H(t) = 0$ для $t < 0$, функции $\alpha (x,y)$ и $\beta (x,y)$ определены при $x \ne y$ равенствами в которых $\xi \in \Gamma (x,y)$ – промежуточная точка интегрирования, а ds – элемент римановой длины, функция $v(x,t,y)$ непрерывна при $t \geqslant \tau (x,y)$ и $v(x,t,y) \to 0$ при $t \to \tau (x,y) + 0$. Что представляет собой функция $J(x,y)$, будет объяснено ниже.

Чтобы получить формулы (14), (15), надо вначале найти для них дифференциальные уравнения и начальные условия. Эти уравнения находятся с помощью обычного метода подсчета сингулярной и конечной амплитуд решения задачи (4) (см., например, [10]). Для этого надо положить в (13) $v(x,t,y) \equiv 0$, подставить полученное равенство в (4) и приравнять нулю коэффициенты при $\delta '(t - \tau (x,y))$ и $\delta (t - \tau (x,y))$. В результате этого возникают следующие уравнения

Для однородной среды с параметрами $c(x) \equiv c(y)$, $\sigma (x) \equiv 0$, $q(x) \equiv 0$ имеет место формула

Поэтому $\alpha (x,y)$ и $\beta (x,y)$ должны удовлетворять следующим условиям при $x \to y$:

Уравнения (16), (17) являются обыкновенными дифференциальными уравнениями вдоль геодезических $\Gamma (x,y)$. Действительно, согласно первому уравнению (11) и принятому обозначению $\eta = {{\nabla }_{x}}\tau (x,y)$, вдоль геодезической $\Gamma (x,y)$ справедливо равенство

в котором $d\tau = ds$ – элемент римановой длины. Удобно преобразовать выражение для $\Delta \tau (x,y)$, входящее в уравнение (17), используя формулу (2.2.35) из книги [10]. Эта формула справедлива вдоль геодезической $\Gamma (x,y)$ и в данном случае имеет вид: в котором – якобиан перехода от римановых нормальных координат $({{\zeta }_{1}},{{\zeta }_{2}},{{\zeta }_{3}})\, = \,\zeta $ к декартовым $({{x}_{1}},{{x}_{2}},{{x}_{3}})$ = x. Напомним, что римановы нормальные координаты $\zeta $ с центром в точке y могут быть вычислены по формуле $\zeta (x,y) = - \frac{{{{c}^{2}}(y)}}{2}{{\nabla }_{y}}{{\tau }^{2}}(x,y)$. Отсюда, в частности, находим, что $J \in {{C}^{4}}({{\mathbb{R}}^{6}})$ и $J(x,x) = 1$.

Из формулы (21) следует, что

С учетом равенств (20), (22), уравнение (16) преобразуется к виду

в котором $\xi = \xi (s) \in \Gamma (x,y)$, а $ds$ – элемент римановой длины. Так как $J(x,x) = 1$ и τ(x, y) ~ ~ ${\text{|}}x - y{\text{|/}}c(y)$ при $x \to y$, то из равенств (18), (23) следует формула (14).

Уравнение (17) с помощью равенства (16) можно записать в виде

Интегрирование этого уравнения с использованием условия (19) приводит к формуле (15). Очевидно, что $(\alpha \tau ) \in {{C}^{4}}({{\mathbb{R}}^{6}})$, $\beta \in {{C}^{2}}({{\mathbb{R}}^{6}})$.

Функция $v(x,t,y)$ является решением задачи Коши:

в которой оператор L определен формулой (4) и

Из равенств (24), (25) и конечности скорости звука следует, что ${v}(x,t,y) \equiv 0$ для $t < \tau (x,y)$, а значит и справедливость представления (13). Кроме того, фундаментальное решение $G(x,t,\xi )$ для уравнения (24), т.е. решение задачи

имеет вид в котором функция $\alpha (x,y)$ определена формулой (14), а ${{G}_{0}}(x,t,\xi )$ – регулярная непрерывная функция своих аргументов. Для решения задачи (24) справедливо интегральное представление решения в виде

Из формулы (25) следует, что $f(x,t{\kern 1pt} ',\xi ) = 0$ для $t{\kern 1pt} ' \leqslant \tau (x,\xi )$. Поэтому равенство (26) преобразуется к виду

в котором $D(x,y,t)$ – множество представляющее собой внутренность риманова эллипсоида с фокусами в точках $x$ и $y$. При $t \to \tau (x,y)$ этот эллипсоид стягивается к геодезической $\Gamma (x,y)$, а объем множества $D(x,y,t)$ стремится к нулю. Отсюда следует, что $v(x,t,y) \to 0$ при t → $\tau (x,y)$. Непрерывность функции $v(x,t,y)$ при $t \geqslant \tau (x,y)$ следует из принадлежности $\beta (x,y)$ пространству ${{C}^{2}}({{\mathbb{R}}^{6}})$.

Замечание 2. При выполнения условия (2) справедливы равенства

Воспользуемся теоремой 1 для анализа обратной задачи об определении коэффициентов уравнения (4) по данным (7). Прежде всего, заметим, что эти данные однозначно определяют функцию $\tau (x,y)$ для всех $y \in {{S}_{R}}$ и $x \in {{S}_{R}}(y,{{r}_{0}})$. Действительно, в силу представления (13), при фиксированных x и $y$ имеет место равенство

Далее, для тех же $x$ и $y$ по функции $\Phi (x,t,y)$ находятся $\alpha (x,y)$ и $\beta (x,y)$:

Доопределим найденные функции на множество ${{S}_{R}} \times {{S}_{R}}$ с помощью формул

Тогда обратная задача об определении коэффициентов $c(x)$, $\sigma (x)$ и q(x) в ${{B}_{{{{r}_{0}}}}}$ сводится к последовательному решению трех задач: 1) отысканию $c(x)$ по заданной функции $\tau (x,y)$, $(x,y) \in {{S}_{R}} \times {{S}_{R}}$, 2) определению коэффициента $\sigma (x)$ по функции $\alpha (x,y)$, $(x,y) \in {{S}_{R}} \times {{S}_{R}}$, 3) отысканию q(x) по заданной функции $\beta (x,y)$, $(x,y) \in {{S}_{R}} \times {{S}_{R}}$.

Первая из этих задач является обратной кинематической задачей. Оценка устойчивости ее решения найдена в статьях [12, 13]. Эта оценка имеет вид

в котором ${{c}_{1}}(x)$ и ${{c}_{2}}(x)$ – положительные функции, генерирующие простые римановы метрики, ${{\tau }_{1}}(x,y)$ и ${{\tau }_{2}}(x,y)$ – отвечающие им римановы расстояния, а $C = C({{c}_{0}},{{c}_{{00}}})$ – некоторая положительная постоянная, зависящяя от чисел ${{c}_{0}},\;{{c}_{{00}}}$, введенных в (8). Из приведенной выше оценки вытекает однозначность решения обратной кинематической задачи.

Вторая и третья задачи приводят к идентичной задаче интегральной геометрии. В самом деле, если функция $c(x)$ найдена, то геодезические линии $\Gamma (x,y)$, а также функция $J(x,y)$ становятся известными. Тогда из формулы (14) однозначно находятся интегралы

в которых

Задача об определении функции $\hat {\sigma }(x)$ = c²(x)σ(x) из уравнения (31) представляет собой задачу интегральной геометрии. Она исследована в работах [13, 14]. Оценка устойчивости решения этой задачи, найденная в цитированных работах, имеет вид, аналогичный оценке для решения обратной кинематической задачи с естественной заменой c и $\tau $ на $\hat {\sigma }$ и g соответственно. Решение задачи интегральной геометрии на семействе геодезических, определяемых простой римановой метрикой, единственно. Найдя $\hat {\sigma }(x)$, а следовательно и $\sigma (x)$, можно вычислить по формуле (14) функцию $\alpha (x,y)$ для любых $x$ и $y$. Тогда из формулы (15) однозначно вычисляются интегралы

в которых

В результате для определения функции $\hat {q}(x)$ = = c²(x)q(x) из уравнения (32) возникает в точности та же задача интегральной геометрии. Решив ее, находим затем и коэффициент q(x). По нему, решая задачу (6), найдем m(x) и, наконец, определим плотность $\rho (x) = {{m}^{2}}(x)$.

Итогом предыдущих рассмотрений является

Теорема 2. При выполнении условий теоремы 1 и условия (2) обратная задача может иметь только одно решение.

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

ИСТОЧНИК ФИНАНСИРОВАНИЯ

Работа выполнена в рамках государственного задания ИМ СО РАН (проект FWNF-2022-0009).