Журнал вычислительной математики и математической физики, 2021, T. 61, № 12, стр. 2095-2108

ВВЕДЕНИЕ

Постановки обратных задач для уравнений переноса излучения весьма разнообразны и, не смотря на достаточно давнюю историю, интерес со стороны специалистов к обратным задачам теории переноса остается стабильно высоким (см. [1]–[20]). В большинстве задач требуется найти коэффициент ослабления по входящему и выходящему из среды потокам излучения. Такая постановка естественна и традициона для томографии и соответствует модели, когда можно пренебречь рассеянием и наличием внутренних источников в среде. В этом случае коэффициент ослабления в уравнении остается единственной неизвестной величиной, определяющей свойства среды. Особенно ценны для практики постановки обратных задач, в которых определяется хотя бы один из коэффициентов уравнения и не требуется знание остальных. Например, в [5]–[14] изучались задачи определения коэффициента ослабления с внешними источниками специального типа, подавляющими влияние рассеяния и внутренних источников в облучаемой среде, а также задачи нахождения поверхностей разрывов коэффициентов уравнения по информации только о выходящем излучении. В частности, в [7], [8] для определения коэффициента ослабления в стационарном уравнении переноса было предложено использовать специальный источник излучения, имеющий разрыв первого рода по угловой переменной на некотором подмножестве единичной сферы, например, в горизонтальном сечении. Доказано, что решение представимо в виде суммы баллистической и рассеянной компоненты, причем баллистическая составляющая потока излучения имеет разрыв первого рода, а рассеянная компонента является непрерывной функцией. В недавних работах (см. [19], [20]) доказаны подобные утверждения в случае, когда плотность потока входящего излучения имеет разрывы по пространственным переменным на некоторой кривой, принадлежащей границе просвечиваемой области.

Используя указанные разложения, преобразование Радона коэффициента ослабления выражается через величины разрывов плотности потока входящего и выходящего из среды излучения. Таким образом, обратная задача для интегро-дифференциального уравнения сводится к традиционной задаче интегральной геометрии – обращению преобразования Радона в сечении исследуемого трехмерного объекта.

В настоящей статье мы обобщим результаты исследований (см. [7], [8]), полученные для стационарного уравнения переноса излучения, на нестационарный случай. Мы исследуем структуру нестационарного решения прямой задачи и получим формулу для преобразования Радона коэффициента ослабления. Справедливости ради заметим, что для нестационарной модели переноса излучения можно добиться дискриминации рассеивающей составляющей в измеряемом сигнале путем уменьшения длительности зондирующего импульса. Действительно, импульсные источники излучения давно и успешно используются в томографии. Более широкое распространение такие источники получили в оптической томографии, так как в рентгеновском диапазоне создание ультракоротких импульсов сопряжено с более жесткими физическими и технологическими ограничениями. Еще одно важное достоинство импульсных источников в сравнении с традиционными непрерывно излучающими генераторами рентгеновских лучей связано со снижением дозы лучевой нагрузки при томографии биологических объектов.

В работе проведено численное тестирование алгоритма решения обратной задач на известном фантоме (см. [21]) и показано, что с уменьшением длительности импульса качество восстановленных изображений повышается. Отметим, что принцип фильтрации рассеянного поля основан на различной гладкости баллистической и рассеянной компонент излучения, поэтому в предложенном методе томографирования среды нет острой необходимости в использовании ультракоротких импульсов (см. [22], [23]). Для получения синтетических данных о выходящем излучении программно реализован метод Монте-Карло решения нестационарного уравнения переноса с источниками излучения специального типа.

1. ПОСТАНОВКА ОБРАТНОЙ ЗАДАЧИ

Рассмотрим интегро-дифференциальное уравнение следующего вида (см. [17], [22], [24]):

Уравнение (1) описывает нестационарный процесс взаимодействия излучения с веществом, а функция $I(r,\omega ,t)$ интерпретируется как плотность потока частиц в момент времени $t \in [0,T]$ в точке $r \in {{\mathbb{R}}^{3}}$, движущихся со скоростью ${v}$ в направлении единичного вектора $\omega \in \Omega = \{ \omega \in {{\mathbb{R}}^{3}}:\left| \omega \right| = 1\} $. Функции $\mu $ и $\sigma $ имеют смысл коэффициентов ослабления и рассеяния, а через $p$ и $J$ обозначены индикатриса рассеяния и плотность внутренних источников.

Процесс переноса излучения происходит в некоторой многокомпонентной системе $G$, состоящей из объединения конечного числа ограниченных и попарно непересекающихся подобластей ${{G}_{1}},\; \ldots ,\;{{G}_{m}}$, причем замыкание $\bar {G}$ является выпуклым множеством в ${{\mathbb{R}}^{3}}$. Поверхность $\partial{ \bar {G}}$ будем называть внешней границей множества $G$, а $\gamma = \partial G{{\backslash }}\partial{ \bar {G}}$ – внутренней границей множества $G$.

Относительно множества $G$ дополнительно предполагаем следующее: любая прямая, имеющая общую точку с $G$, пересекает границу $\partial G$ в конечном числе точек. Это условие, называемое условием обобщенной выпуклости, типично в теории переноса и необременительно с прикладной точки зрения. Возникающая, казалось бы, несогласованность, связанная с наличием линейных участков на поверхности $\partial G$, во многих случаях может быть преодолена продолжением этих участков с возможным увеличением числа областей ${{G}_{i}}$.

Символом ${{\operatorname{mes} }_{m}}(X)$ в дальнейшем будем обозначать меру Лебега множества $X$ в ${{\mathbb{R}}^{m}}$. Из условия обобщенной выпуклости нетрудно получить, что $\mathop {\operatorname{mes} }\nolimits_3 (\partial G) = 0$ (см. [11]).

Обозначим через ${{L}_{{r,\omega }}}$ луч, исходящий из точки $r \in {{\mathbb{R}}^{3}}$ в направлении $\omega $, ${{L}_{{r,\omega }}} = {\text{\{ }}r + \omega t:t > 0{\text{\} }}$, а   через $d(r,\omega )$ – расстояние от точки $r \in \bar {G}$ до границы $\partial{ \bar {G}}$ в направлении $\omega $, т.е. $d(r,\omega ) = {\text{me}}{{{\text{s}}}_{1}}({{L}_{{r,\omega }}} \cap \bar {G})$. Пусть множество ${{\Gamma }_{{ \pm \omega }}}$ образовано точками $z \in \partial{ \bar {G}}$, которые представимы в виде

Посредством ${{\Gamma }_{\omega }}$ построим множества ${{\Gamma }^{ \pm }} = {{\Gamma }_{{ \pm \omega }}} \times \Omega $, $\Gamma = {{\Gamma }^{ + }} \cup {{\Gamma }^{ - }}$ и для краткости изложения введем ряд обозначений:

Присоединим к уравнению (1) начальные и граничные условия

Для удобства введем в рассмотрение функцию и запишем начальное условие (2) и граничное условие (3) в виде

Задача 1. В прямой задаче требуется определить функцию $I$ из уравнения (1) и условия (4) при заданных $\mu $, $\sigma $, $J$, $p$, ${v}$, $h$.

Исследованию прямых задач для стационарных и нестационарных уравнений переноса излучения подобного рода посвящено не мало работ, отметим лишь несколько широко известных монографий по этой тематике, а именно, [24]–[29].

Задача 2. Под обратной задачей будем понимать задачу определения функции $\mu $ из соотношений (1), (4) и дополнительного условия

при заданных $h$, $H$, ${v}$.

Для стационарного уравнения переноса излучения задача 2 при специальных условиях на внешний источник излучения рассматривалась в [5], [7], [8], [11], [14], [19], [20]. В настоящей работе результаты из [7], [8], [11] не только обобщаются на нестационарный случай, но и проводится численный анализ качества томограмм коэффициента ослабления при уменьшении длительности импульса внешнего источника излучения.

2. ФУНКЦИОНАЛЬНЫЕ ПРОСТРАНСТВА. ОСНОВНЫЕ ОГРАНИЧЕНИЯ

Обозначим через ${{C}_{b}}(X)$ пространство непрерывных и ограниченных на $X$ функций. Относительно неотрицательных функций $\mu (r)$, $\sigma (r)$, $J(r,\omega ,t)$, $p(r,\omega \cdot \omega {\text{'}})$ будем предполагать, что $\mu \geqslant {\text{const}} > 0$, $\mu ,\sigma \in {{C}_{b}}(G)$, $J \in {{C}_{b}}(X)$, а функция $p(r,\omega \cdot \omega {\text{'}}) \in {{C}_{b}}(G \times [ - 1,\;1])$ и для всех $r$ удовлетворяет условию нормировки

Под дифференциальным выражением в левой части уравнения (1) будем понимать производную в точке $(r,t)$ в направлении вектора $({{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}},1{\text{/}}{v})$:

и через ${{d}^{ \pm }}(r,\omega ,t)$ обозначим

Будем говорить, что функция $f(r,\omega ,t)$ принадлежит $D(X)$, если выполнены следующие условия:

1) при всех $(r,\omega ,t) \in X$ функция $f(r + \tau \omega ,\omega ,t + \tau {\text{/}}{v})$ абсолютно непрерывна по $\tau $, $\tau \in [ - {{d}^{ - }}(r,\omega ,t),{{d}^{ + }}(r,\omega ,t)]$,

2) $f(r - {{d}^{ - }}(r,\omega ,t)\omega ,\omega ,t - {{d}^{ - }}(r,\omega ,t){\text{/}}{v}) = 0$,

3) $\left( {\tfrac{1}{{v}}\tfrac{\partial }{{\partial t}} + \omega \cdot {{\nabla }_{r}} + \mu } \right)f \in {{C}_{b}}(X)$.

Определим оператор $\mathcal{L}:D(X) \to {{C}_{b}}(X)$ следующим равенством:

Обозначим через ${{\Omega }_{0}}$ некоторое подмножество нулевой меры на $\Omega $ и рассмотрим оператор $\mathcal{S}:{{C}_{b}}(G \times (\Omega {{\backslash }}{{\Omega }_{0}}) \times [0,T]) \to {{C}_{b}}(X)$, заданный соотношением Так как $\sigma (r) \in {{C}_{b}}(G)$, $p(r,\omega \cdot \omega {\text{'}}) \in {{C}_{b}}(G \times [ - 1,\;1])$, то при $f \in {{C}_{b}}(G \times \Omega {{\backslash }}{{\Omega }_{0}} \times [0,T])$ функция $\mathcal{S}f$ действительно принадлежит пространству ${{C}_{b}}(G \times \Omega \times [0,T])$ (см. [11]), следовательно, определение оператора $\mathcal{S}$ корректно.

В следующем разделе рассматривается задача нахождения функции $f \in D(X)$, удовлетворяющей уравнению

на множестве $X$. Из определения функционального пространства $D(X)$ вытекает, что функция $f$ является решением прямой задачи (1), (4) с однородными граничными и начальным условиями ($h = 0$), а при переходе через границу раздела сред ($r \in \partial G{{\backslash }}\partial{ \bar {G}}$) решение остается непрерывно вдоль направления $\omega $.

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

Предварительно докажем несколько вспомогательных лемм. Рассмотрим некоторое множество $G{\text{'}} \times \Omega $, состоящее из таких точек $(r,\omega ) \in \bar {G} \times \Omega $, что прямая ${\text{\{ }}r + \tau \omega ,\; - {\kern 1pt} \infty < \tau < + \infty {\text{\} }}$ имеет непустое пересечение с $G$. В силу условия обобщенной выпуклости множество $G$ является подмножеством $G{\text{'}}$.

Лемма 1. Пусть функция $\mu (r) \in {{C}_{b}}(G)$ и $\phi (r,\omega ,t) \in {{C}_{b}}(X)$, тогда функция

принадлежит пространству ${{C}_{b}}(G{\text{'}} \times \Omega \times [0,d(r,\omega )])$ и функцияпринадлежит пространству ${{C}_{b}}(X{\kern 1pt} {\text{'}})$, где $X{\kern 1pt} {\text{'}} = G{\kern 1pt} {\text{'}} \times \Omega \times [0,T]$.

Доказательство. Докажем первое утверждение леммы. Для начала заметим, поскольку $\tau \leqslant d(r,\omega ) \leqslant d$, где $d$ – диаметр множества $\bar {G}$, то из ограниченности функции $\mu $ на множестве $G$ следует ограниченность функции $a$. Докажем непрерывность функции $a(r,\omega ,\tau )$ на $G{\text{'}} \times \Omega \times [0,d(r,\omega )]$.

Фиксируем произвольную точку $(r,\omega ,\tau ) \in G{\text{'}} \times \Omega \times [0,d(r,\omega )]$ и для любой последовательности точек $({{r}_{n}},{{\omega }_{n}},{{\tau }_{n}}) \in G{\text{'}} \times \Omega \times [0,d(r,\omega )]$, сходящихся к точке $(r,\omega ,\tau )$ при $n \to \infty $, рассмотрим выражения

где ${{\psi }_{n}}(\tau {\text{'}}) = {{\tau }_{n}}\mu ({{r}_{n}} + \tau {\text{'}}{{\tau }_{n}}{{\omega }_{n}})$ и $\psi (\tau {\text{'}}) = \tau \mu \left( {r + \tau {\text{'}}\tau \omega } \right)$.

Заметим, что $\left| {{{\psi }_{n}}(t{\kern 1pt} {\text{'}})} \right| \leqslant {\text{const}}$ для любого $n$. Докажем, что при $n \to \infty $ последовательность ${{\psi }_{n}}(t{\kern 1pt} {\text{'}})$ стремится к $\psi (t{\kern 1pt} {\text{'}})$ почти всюду на $[0,\;1]$. В силу условия обобщенной выпуклости при любом $n$ луч ${{L}_{{{{r}_{n}},{{\omega }_{n}}}}} = {\text{\{ }}{{r}_{n}} + \tau {\text{'}}{{\omega }_{n}},\tau {\text{'}} \geqslant 0{\text{\} }}$ пересекает границу $\partial G$ в конечном числе точек. Пусть ${{r}_{n}} + {{\tau }_{{i,n}}}{{\tau }_{n}}{{\omega }_{n}}$ ($i = 1,\;2,\; \ldots ,\;{{q}_{n}}$, ${{q}_{n}} < \infty $) – точки пересечения отрезка ${\text{\{ }}{{r}_{n}} + \tau {\text{'}}{{\tau }_{n}}{{\omega }_{n}},\;0 \leqslant \tau {\text{'}} \leqslant 1{\text{\} }}$ с границей $\partial G$, тогда множество

является счетным подмножеством отрезка $[0,\;1]$, поскольку представляет собой счетное объединение конечных множеств.

Таким образом, поскольку $\mu (r) \in {{C}_{b}}(G)$, то все функции ${{\psi }_{j}}(\tau {\text{'}})$ из последовательности ${\text{\{ }}{{\psi }_{n}}(\tau {\text{'}}){\text{\} }}$, $n = 1,\;2,\; \ldots $, непрерывны по $\tau {\text{'}}$, $\tau {\text{'}} \in [0,\;1]{{\backslash }}\Pi $. Следовательно, ${{\psi }_{n}}(\tau {\text{'}}) \to \psi (\tau {\text{'}})$ почти всюду на $[0,\;1]$. Тогда по теореме Лебега о предельном переходе под знаком интеграла получаем

В итоге, мы доказали непрерывность функции $a(r,\omega ,\tau )$ на множестве $G{\text{'}} \times \Omega \times [0,d(r,\omega )]$.

Докажем второе утверждение леммы. Покажем, что функция ${{d}^{ \pm }}(r,\omega ,t)$ принадлежит ${{C}_{b}}(X{\kern 1pt} {\text{'}})$. Так как множество $\bar {G}$ выпуклое, то функцию $d(r,\omega )$ на этом множестве можно определить следующим образом:

где $\chi (r)$ – характеристическая функция множества $\bar {G}$, диаметр которого равен $d$. Поскольку $\chi (r) \in {{C}_{b}}(G)$, то, согласно доказанному утверждению в первой части леммы, функция $d(r,\omega ) \in {{C}_{b}}(G{\text{'}} \times \Omega ) \subset {{C}_{b}}(G \times \Omega )$.

По определению ${{d}^{ + }}(r,\omega ,t) = min{\text{\{ }}d(r,\omega ),{v}(T - t){\text{\} }}$ и ${{d}^{ - }}(r,\omega ,t) = min{\text{\{ }}d(r, - \omega ),{v}t{\text{\} }}$, следовательно, функции ${{d}^{ \pm }}(r,\omega ,t)$ принадлежат пространству ${{C}_{b}}(X{\kern 1pt} {\text{'}})$.

Покажем непрерывность функции ${{\Phi }^{ + }}$, непрерывность ${{\Phi }^{ - }}$ доказывается аналогично. Ограниченность функции ${{\Phi }^{ + }}$ следует из ограниченности функций $a(r,\omega ,t)$, $\phi (r,\omega ,t)$, ${{d}^{ + }}(r,\omega ,t)$. Для доказательства непрерывности функции ${{\Phi }^{ + }}$ на множестве $X{\kern 1pt} {\text{'}}$ зафиксируем произвольную точку $(r,\omega ,t) \in X{\kern 1pt} {\text{'}}$ и для любой последовательности точек $({{r}_{n}},{{\omega }_{n}},{{t}_{n}}) \in X{\kern 1pt} {\text{'}}$, $({{r}_{n}},{{\omega }_{n}},{{t}_{n}}) \to (r,\omega ,t)$ при $n \to \infty $, рассмотрим выражения

где и Так как ${{d}^{ + }}(r,\omega ,t) \in {{C}_{b}}(X{\text{'}})$, то аналогично доказанному в первой части утверждения показывается, что ${{\psi }_{n}}(t{\kern 1pt} {\text{'}}) \to \psi (t{\kern 1pt} {\text{'}})$ при $n \to \infty $ почти всюду на $[0,\;1]$. Из теоремы Лебега заключаем, что функция ${{\Phi }^{ + }}(r,\omega ,t)$ непрерывна на множестве $X{\text{'}}$. Утверждения леммы доказаны.

Лемма 2. При любой функции $F \in {{C}_{b}}(X)$ существует единственное решение $f \in D(X)$ уравнения

которое при всех $(r,\omega ,t) \in X$ задается формулойи удовлетворяет условиюгде $\bar {\mu } = \mathop {inf}\limits_{r \in G} \left| {\mu (r)} \right|$.

Доказательство. Непосредственно проверяется, что функция $f$ из (12) удовлетворяет уравнению (11), условиям 1) и 2) из определения пространства $D$. По лемме 1 функция $a(r,\omega ,\tau )$, определенная соотношением (9), принадлежит пространству ${{C}_{b}}(G{\text{'}} \times \Omega \times [0,d(r,\omega )])$, поэтому $a(r, - \omega ,\tau )$ принадлежит ${{C}_{b}}(G{\text{'}} \times \Omega \times [0,d(r, - \omega )])$. Функция $F(r,\omega ,t)$ принадлежит ${{C}_{b}}(X)$ по условию, следовательно, согласно лемме 1, функция $f$ в (12) принадлежит пространству ${{C}_{b}}(X)$. Функция $\left( {\tfrac{1}{{v}}\tfrac{\partial }{{\partial t}} + \omega \cdot {{\nabla }_{r}} + \mu } \right)f$ принадлежит ${{C}_{b}}(X)$, поскольку $f$ удовлетворяет уравнению (11), в котором правая часть $F \in {{C}_{b}}(X)$.

Для доказательства оценки (13) воспользуемся представлением (12):

Лемма доказана.

Из утверждений леммы 2 вытекает, что оператор ${{\mathcal{L}}^{{ - 1}}}:{{C}_{b}}(X) \to D(X)$, определенный формулой

существует и ограничен в ${{C}_{b}}(X)$.

В пространстве $D(X)$ введем норму

Из (13) вытекает очевидное неравенство на основании которого заключаем, что из сходимости в пространстве $D(X)$ следует сходимость в ${{C}_{b}}(X)$. Нетрудно показать, что линейное множество функций $D$, наделенное нормой (16), образует банахово пространство.

Из леммы 2 вытекает, что решение уравнения (8) в $D$ эквивалентно решению уравнения

Теорема 1. Уравнение $(8)$ в пространстве $D$ однозначно разрешимо.

Доказательство. Поскольку функция $p$ неотрицательна и удовлетворяет условию нормировки и, кроме того, справедливо неравенство (17), то для ${{\left\| {{{\mathcal{L}}^{{ - 1}}}\mathcal{S}f} \right\|}_{D}}$ получаем

где $\overline \lambda = \left\| {\tfrac{\sigma }{\mu }} \right\|$. Так как $\mu (r) \geqslant \sigma (r)$, то $\overline \lambda \leqslant 1$ и из (19) вытекает, что Так как норма оператора ${{\mathcal{L}}^{{ - 1}}}\mathcal{S}$, действующего в банаховом пространстве, меньше единицы, то уравнение (18) однозначно разрешимо и решение может быть найдено методом последовательных приближений где ${{f}_{0}} = {{\mathcal{L}}^{{ - 1}}}J$. Теорема 1 доказана.

Отметим, что условие $\overline \lambda \leqslant 1$ излишне для существования единственного ограниченного решения на конечном промежутке времени $t \in [0,\;T]$. Действительно, если $\sigma > \mu $, то в уравнении (8) можно сделать замену $f = {{f}_{\lambda }}{{e}^{{\lambda t}}}$, где $\lambda $ – произвольное число, удовлетворяющее условиям $\mu (r) + \lambda {\text{/}}{v} > 0$, $\tfrac{{\sigma (r)}}{{\mu (r) + \lambda {\text{/}}{v}}} \leqslant 1$ для всех $r \in G$. Тогда из теоремы 1 при соответствующем выборе нормы в пространстве $D$, например,

вытекает существование и единственность решения ${{f}_{\lambda }}$, а следовательно, и решения $f$ уравнения (8). Тем не менее ограничение $\overline \lambda \leqslant 1$ для неограниченного временного интервала $t \in [0,\infty )$ гарантирует стабилизацию нестационарного решения $f$ при $t \to \infty $.

4. ИССЛЕДОВАНИЕ ОБРАТНОЙ ЗАДАЧИ ДЛЯ УРАВНЕНИЯ ПЕРЕНОСА ИЗЛУЧЕНИЯ

В предыдущем разд. 3 мы рассмотрели начально-краевую задачу для уравнения переноса излучения с однородными начальными и граничными условиями. Решение неоднородной начально-краевой задачи (1), (4) можно представить в виде $I = f + {{I}_{0}}$, где

и функция $f$ принадлежит $D(X)$ и удовлетворяет уравнению Так как оператор $\mathcal{S}:{{C}_{b}}(G \times (\Omega {{\backslash }}{{\Omega }_{0}}) \times [0,T]) \to {{C}_{b}}(X)$ обладает сглаживающими свойствами и $J \in {{C}_{b}}(X)$, то функция $\mathcal{S}{{I}_{0}} + J$ принадлежит ${{C}_{b}}(X)$, даже если функция $h$ имеет разрывы по переменной $\omega $. Следовательно, согласно результатам предыдущего раздела, уравнение (21) однозначно разрешимо в $D(X)$. Таким образом, функция $I$ может быть представлена в виде суммы двух функций – $f$ и ${{I}_{0}}$. Функция $f$, имеющая смысл рассеянного поля, непрерывна на множестве $X{\kern 1pt} {\text{'}} \subset \bar {G} \times \Omega \times [0,T]$, а нерассеянное излучение ${{I}_{0}}$ может содержать разрывы по переменной $\omega $ на множестве $X{\kern 1pt} {\text{'}}$. Такая особенность структуры решения уравнения переноса излучения лежит в основе метода решения обратной задачи.

Обозначим величину разрыва функции $h$ по переменной $\omega $ при $\omega \to {{\omega }^{0}} = ({{\omega }_{1}},{{\omega }_{2}},0)$ через

и наложим ограничение на функцию $h$: Тогда для граничных точек $(\eta ,{{\omega }^{0}},t) \in {{Y}^{ + }}$, удовлетворяющих ограничению $d(\eta , - {{\omega }^{0}}) \leqslant {v}t \leqslant {v}T - d(\eta , - {{\omega }^{0}})$, получаем Так как $f \in {{C}_{b}}(X{\text{'}})$, то $[f](\eta ,{{\omega }^{0}},t) = 0$, и тогда из соотношения (23) при выполнении условия (22) вытекает равенство В левой части равенства (24) стоит преобразование Радона функции $\mu $ в любой горизонтальной плоскости ${{r}_{3}} = {\text{const}}$. После замены переменных соотношение (24) можно переписать в виде где $r$ – любая точка области $G$, а переменная $t$ принадлежит временному промежутку $[d(r, - {{\omega }^{0}}){\text{/}}{v},\;T - d(r,{{\omega }^{0}}){\text{/}}{v}]$ при достаточно большом $T$.

По аналогии со стационарным случаем (см. [8], [11]) доказывается следующая теорема единственности решения обратной задачи.

Теорема 2. Пусть при некотором $t \in (0,T)$ справедливо

где $I{\text{'}}$, $I{\text{''}}$ – решения начально-краевых задач $(1)$–$(3)$ для двух совокупностей коэффициентов $\left\{ {\mu {\text{'}}(r),\sigma {\text{'}}(r),p{\text{'}}(\omega \cdot \omega {\text{'}}),J{\text{'}}(r,\omega ,t)} \right\}$ и $\left\{ {\mu {\text{''}}(r),\sigma {\text{''}}(r),p{\text{''}}(\omega \cdot \omega {\text{'}}),J{\text{''}}(r,\omega ,t)} \right\}$ с одной и той же функцией $h(\xi ,\omega ,t)$, удовлетворяющей условию $(22)$, тогда справедливо соотношение $(25)$ и $\mu {\text{'}}(r) = \mu {\text{''}}(r)$ почти всюду в $G$.

Особенность теоремы заключается в том, что из двух наборов коэффициентов уравнения (1):

утверждается совпадение только $\mu {\text{'}}$ и $\mu {\text{''}}$, а об остальных ничего не говорится. Из формулы (25) следует, что неизвестные функции $\sigma $, $p$, $J$ не влияют на процедуру определения функции $\mu $. Привлекая физическую терминологию, можно сказать, что влияние рассеяния и наличие радиоактивных источников в среде подавляются за счет выбора внешнего источника. Так как функция $h$ имеет разрыв в точке ${{\omega }^{0}} = ({{\omega }_{1}},{{\omega }_{2}},0)$, то функция $\mu (r)$ восстанавливается “послойно” в горизонтальных плоскостях ${{r}_{3}} = {\text{const}}$, что является традиционным в рентгеновской томографии (см. [30]).

Время $t$ в соотношениях (25) и (26) может быть любым, лишь бы в этой точке выполнялось условие (22), гарантирующее наличие ненулевого разрыва функции $h$. Однако при численном решении обратной задачи (и при практической реализации метода) желательно, чтобы скачки функций $h$ и $H$ были максимальными в моменты времени $t$ и $t + d(\eta , - {{\omega }^{0}}){\text{/}}{v}$ соответственно. Например, если среда зондируется импульсом, имеющим максимум интенсивности баллистической компоненты в момент времени $t = {{t}_{0}}$, то в точке $(\eta ,{{\omega }^{0}}) \in {{\Gamma }^{ + }}$ выходящий из среды сигнал предпочтительнее измерять в момент времени $t = {{t}_{0}} + d(\eta , - {{\omega }^{0}}){\text{/}}{v}$. В этот момент времени отношение компоненты нерассеянного излучения к компоненте рассеянного поля будет максимальным в точке $(\eta ,{{\omega }^{0}})$.

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

Для проверки алгоритма решения обратной задачи 2 необходимо знание функции $H$, которая в физическом эксперименте является измеряемой величиной. При математическом моделировании, чтобы найти $H$, нужно решить прямую задачу 1 при всех заданных коэффициентах уравнения (1) и функции $h$. Среди большого разнообразия численных методов решения уравнений переноса излучения в многомерном случае методам Монте-Карло альтернативы практически не существует. Мы будем использовать одну из модификаций метода Монте-Карло, называемую методом “сопряженных блужданий”, применение которой оправдано при нахождении плотности потока излучения в фиксированной точке фазового пространства $X$ (см. [31], [32]).

Обозначим через $\mathcal{P}$ оператор $\mathcal{P}:{{C}_{b}}({{Y}^{ - }}) \to D(X)$, определенный формулой

тогда приближенное решение задачи (1)–(3) с неоднородным граничным условием можно записать в виде усеченного ряда Неймана Функция ${{I}_{N}}$ является аппроксимацией решения начально-краевой задачи $I$, и представление (28) служит основой для построения алгоритма Монте-Карло. Каждый элемент последовательности ${{I}_{n}}$, $n = 1,\;2, \ldots ,\;N$, определяет вклад излучения, испытавшего от $1$ до $n$ актов рассеяния, а слагаемое ${{I}_{0}}$ учитывает вклад нерассеянного излучения.

Для вычисления усеченного ряда Неймана (28), представляющего собой сумму многомерных интегралов, воспользуемся методом Монте-Карло. Построим траекторию частицы

где точки марковской цепи (29) определяются по следующему правилу: В (30) случайная величина ${{\tau }_{{i + 1}}}$ распределена на промежутке $[0,{{d}^{ - }}({{r}_{i}},{{\omega }_{i}},{{t}_{i}})]$ с плотностью вероятности а случайный вектор ${{\omega }_{{i + 1}}}$ распределен на единичной сфере $\Omega $ c переходной плотностью вероятности $p({{r}_{{i + 1}}},{{\omega }_{i}} \cdot {{\omega }_{{i + 1}}})$.

Пусть $\mathbb{E}[\Theta ]$ – математическое ожидание случайной величины $\Theta $, тогда, согласно методу Монте-Карло, для значения функции ${{I}_{N}}$ в точке $(r,\omega ,t)$ можно записать выражение

где ${{\theta }_{0}} = {{I}_{0}}({{r}_{0}},{{\omega }_{0}},{{t}_{0}})$ и случайные величины ${{\theta }_{n}}$ при $n > 0$ определяются рекуррентным образом: Повторяя описанную процедуру $M$ раз, получим выборку объемом $M$ для случайной величины ${{\Theta }_{N}}$. Среднее значение выборки дает оценку математического ожидания случайной величины ${{\Theta }_{N}}$ и, следовательно, приближенное значение ${{I}_{N}}$ в точке $(r,\omega ,t)$. Простой анализ рекуррентных соотношений (33) показывает, что реализованный весовый метод Монте-Карло учитывает поглощение и “вылет” частиц из области. Учет этих эффектов приводит к относительно небольшому увеличению трудоемкости метода, но в значительной степени способствует уменьшению дисперсии оценки математического ожидания случайной величины ${{\Theta }_{N}}$ и повышению точности решения уравнения переноса излучения (см. [32]).

6. ВЫЧИСЛИТЕЛЬНЫЕ ЭКСПЕРИМЕНТЫ

В данном разделе приведены результаты численных экспериментов по восстановлению внутренней структуры среды. Тестирование алгоритма проводилось в два этапа. На первом этапе для заданных параметров среды с помощью метода Монте-Карло, описанного в предыдущем разделе, вычислялось выходящее излучение $H$. При решении прямой задачи предполагалось, что нестационарный источник облучения $h$ исследуемой среды импульсного типа (различной длительности) имел разрыв первого рода по угловой переменной в плоскости ${{\omega }_{3}} = 0$. На втором этапе решалась обратная задача с использованием формулы (25), осуществляющей фильтрацию рассеянного излучения, либо с помощью формулы классической томографии

Формула (34) отличается от (25) тем, что под знаком логарифма стоят значения функций $h$, $H$ в точках границы области, а не величины скачков этих функции. Для нахождения функции $\mu $ из уравнений (25), (34) можно использовать различные алгоритмы обращения преобразования Радона. Например, широко распространенный в современной литературе алгоритм свертки и обратной проекции (см. [30]).

Тестирование алгоритма решения обратной задачи проводилось на предложенном в [21] фантоме, который специально сконструирован для формальной процедуры оценки качества различных алгоритмов компьютерной томографии. Фантом представляет собой цилиндр высотой и диаметром в $10$ см. Внутренний объем цилиндра поделен на пять цилиндрических слоев, каждый из которых играет свою роль при тестировании тех или иных качеств алгоритма реконструкции.

При проведении экспериментов мы ограничились первой секцией фантома, предназначенной для тестирования линейности (см. [33]) восстановленных значений коэффициента ослабления. Данная характеристика является одной из наиболее важных при диагностике и планировании лечения (см. [34]), и сущность данного показателя заключается в следующем. Традиционно в компьютерной томографии значения коэффициента ослабления для различных материалов приводятся в единицах Хаунсфилда ($HU$). Они связаны со значениями коэффициента ослабления вещества $\mu $ по линейному закону

где $\mu $ – коэффициент ослабления вещества, а ${{\mu }_{{{\text{water}}}}}$, ${{\mu }_{{{\text{air}}}}}$ – коэффициенты ослабления для воды и воздуха соответственно. Применение шкалы Хаунсфилда хорошо себя зарекомендовало в классической компьютерной томографии при диагностировании мягких тканей. В случае восстановления структуры более плотных тканей, таких как кости, относительно малая энергия применяемого излучения и увеличение доли рассеянного излучения приводят к тому, что включения, имеющие одинаковую плотность, дают различные значения коэффициентов ослабления в шкале Хаунсфилда в зависимости от того, где расположено включение. В результате линейная связь (35) между единицами Хаунсфилда и коэффициентом ослабления среды нарушается. Особенно остро эта проблема стоит в конусно-лучевой компьютерной томографии (см. [35]).

Схематичное изображение используемого фантома изображено на фиг. 1. Данная секция представляет собой цилиндр диаметром $100$ мм и высотой $17.5$ мм, выполненный из пластика с коэффициентом ослабления, эквивалентным воде ($HU = 0$). Секция содержит еще четыре цилиндрических включения диаметром $13$ мм и высотой $17.5$ мм с различными значениями коэффициента ослабления ($ - 1000,\;30,\; - 30,\;100\;HU$). Следуя методике, предложенной в [21], для проверки линейности вычислялись средние значения коэффициента ослабления для каждого из материалов, входящих в фантом, после этого полученные значения сравнивались с эталонными значениями в единицах Хаунсфилда. Источником излучения выбирался Гауссиан по времени, имеющий вид

где ${{I}_{0}}$ – амплитуда импульса, ${{t}_{0}}$ – момент времени, соответствующий максимальной мощности сигнала, величина ${{t}_{p}}$ определяет длительность импульса на половине амплитуды, и В экспериментах были выбраны следующие значения указанных величин: ${{I}_{0}} = 32\,000$, ${{t}_{0}} = 3$ нс, а величина ${{t}_{p}}$ варьировалась в экспериментах, принимая значения $3$, $30$, $300$, $3000$ пс.

Результаты реконструкции фантома при различной длительности зондирующего импульса приведены на фиг. 2: слева приведены результаты восстановления источника с помощью обращения преобразования Радона, используя формулу (34), а справа – после фильтрации рассеяния за счет применения источника специального типа, основанные на формуле (25). Средние значения восстановленных коэффициентов ослабления в единицах Хаунсфилда приведены в табл. 1. При этом слева указано среднее значение, вычисленное по формуле (34) без предварительной фильтрации рассеяния, а справа приведено значение, полученное при расчетах по формуле (25).

Как видно из табл. 1 и приведенных томограмм, облучение короткими импульсами позволяет фильтровать рассеянную составляющую в выходящем излучении. Поэтому при облучении фантома короткими импульсами получается реконструкция приличного качества даже без применения дополнительной фильтрации рассеяния по формуле (25). В целом погрешность восстановления коэффициента ослабления – на уровне, полученном в [21]. Отметим, что даже в этом “хорошем” случае дополнительная обработка сигнала с помощью формулы (25) дает лучшее соответствие между эталонными и восстановленными значениями коэффициента ослабления. При увеличении длительности импульса качество реконструкции без предварительной обработки сигнала стремительно падает. При длительности импульса $3000$ пс не только наблюдается огромное различие в численных значениях восстановленных коэффициентов, но и само наличие трех включений из четырех визуально определяется с трудом. Применение источника специального типа с последующей обработкой, напротив, дает хорошие результаты. Как видно на фиг. 2 справа, внутренняя структура фантома восстановилась с приемлемым качеством при любой длительности импульса. Однако по мере роста длительности импульса качество реконструкции незначительно, но все равно падает. Таким образом, для достижения наилучшего качества реконструкции рекомендуется комбинировать импульсный источник с источником специального типа.

ЗАКЛЮЧЕНИЕ

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