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

1. ВВЕДЕНИЕ

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

В рентгеновской томографии существует множество алгоритмов подавления деструктивного влияния рассеяния. Основная часть из них основана на компенсационном учете рассеянного сигнала и коррекции проекционных данных [1–6]. Другая часть алгоритмов базируется на конструировании специальных источников излучения и последующей обработке прошедшего через среду излучения [7–13]. Обе группы алгоритмов являются приближенными и, либо требуют знания значительной априорной информации о характеристиках среды и больших вычислительных ресурсов, либо предъявляют высокие аппаратные требования к конструкции источников и детекторов излучения.

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

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

Идея использования ультракоротких импульсов для повышения качества томографических изображений не нова и подобные методы успешно используются в диффузионной оптической томографии [15]. Применение коротких импульсов либо серий коротких импульсов разных длительностей позволяет добиться увеличения пространственного разрешения томографических методов диффузионной томографии, что является критической проблемой для медицинского применения данных методов. В случае рентгеновской томографии проблема пространственного разрешения не столь критична и основное внимание направлено на повышение точности восстановления численных значений коэффициента ослабления. Существенный вклад в проекционные данные рассеянного излучения от неколлимированного источника приводит к тому, что реконструкция идентичных по химическому составу включений в зависимости от их пространственного расположения дает различные значения коэффициентов ослабления по шкале Хаунсфилда. Особенно остро проблема паразитного влияния рассеянного излучения стоит в конусно-лучевой компьютерной томографии, где эффект коллимации зондирующего излучения незначителен [16]. Преимуществом предлагаемого алгоритма над существующими является возможность повышения качества томографических изображений, оставаясь в рамках имеющейся на сегодняшний день приборной базы и вычислительных мощностей и, несмотря на необходимость проведения повторных просвечиваний, применение импульсов малой длительности понижает общую лучевую нагрузку в сравнении непрерывным режимом облучения.

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

Рассмотрим интегро-дифференциальное уравнение, описывающее нестационарный процесс взаимодействия излучения с веществом, вида [17]

Функция $I(r,\omega ,t)$ интерпретируется как плотность потока частиц в момент времени $t \in [0,T]$ в точке $r \in {{R}^{3}}$, движущихся со скоростью с в направлении единичного вектора $\omega \in \Omega = \{ \omega \in {{R}^{3}}$ : |ω| = 1}. Функции $\mu = \mu (r)$ и $\sigma = \sigma (r)$ имеют смысл коэффициентов ослабления и рассеяния, а через $p = p(r,\omega \cdot \omega {\kern 1pt} ')$ обозначена индикатриса рассеяния. Исследуемый объект содержится в шаре G с центром в начале координат и диаметром d, G = = $\{ r \in {{R}^{3}}:\left| r \right| = d{\text{/}}2\} $. Обозначим через ${{\Pi }_{\omega }}$ плоскость, касательную к границе области G и перпендикулярную направлению ω, Π_ω = $\{ r\, \in \,{{R}^{3}}:r \cdot \omega $ = = d/2}. Очевидно, что плоскости ${{\Pi }_{{ - \omega }}}$ и ${{\Pi }_{\omega }}$ параллельны и находятся на расстоянии d друг от друга. Сканирование объекта осуществляется путем синхронного поворота плоскостей с источниками излучения ${{\Pi }_{{ - \omega {\kern 1pt} *}}}$ и детекторами ${{\Pi }_{{ + \omega *}}}$ в горизонтальных сечениях $\omega {\kern 1pt} * \in \Omega {\kern 1pt} *\, = \,\{ \omega \, = \,({{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}})\, \in \,\Omega :{{\omega }_{3}}$ = 0}.

В задачах томографии обычно требуется определить хотя бы один из коэффициентов уравнения (1), если для каждого $\omega {\kern 1pt} * \in \Omega {\kern 1pt} *$ его решение ${{I}_{{\omega {\kern 1pt} *}}}$ известно на множестве ${{\Pi }_{{ - \omega {\kern 1pt} *}}} \times \Omega {\kern 1pt} *\; \times \left[ {0,T} \right]$ и ${{\Pi }_{{ + \omega {\kern 1pt} *}}} \times \Omega {\kern 1pt} *$ × [0, T].

Рассматривая начально-краевую задачу для уравнения (1), мы будем опускать параметрическую зависимость решения уравнения ${{I}_{{\omega {\kern 1pt} *}}}(r,\omega ,t)$ от направления $\omega {\kern 1pt} *$. Для краткости изложения введем ряд обозначений: $X = G \times \Omega \times \left[ {0,T} \right]$, ${{X}_{0}} = G \times \Omega $ × {t = 0}, ${{\Omega }_{{ - \omega {\kern 1pt} *}}} = \{ \omega \in \Omega : - \omega {\kern 1pt} *\; \cdot \omega > 0\} $, ${{Y}^{ - }} = {{\Pi }_{{ - \omega {\kern 1pt} *}}} \times {{\Omega }_{{ - \omega {\kern 1pt} *}}} \times \left[ {0,T} \right]$, ${{X}^{ - }} = {{Y}^{ - }} \cup {{X}_{0}}$. Присоединим к уравнению (1) начальные и граничные условия

где ${{h}_{0}} = {{h}_{0}}(r,\omega )$ характеризует состояние процесса в начальный момент времени, а функция h_ext = = ${{h}_{{ext}}}(z,\omega ,t)$ описывает плотность потока излучения, падающего на объект. На множестве ${{X}_{0}} \times {{Y}^{ - }}$ построим функцию и объединим условия (2) и (3) в одно начально-краевое условие следующего вида

В прямой задаче требуется определить функцию I из уравнения (1) при условии (4) при заданных функциях μ, σ, p, c, h. Исследованию прямых задач для стационарных и нестационарных уравнений переноса излучения подобного рода посвящено немало работ, отметим лишь несколько широко известных монографий по этой тематике [17–19].

Переходя к постановке обратной задачи для простоты изложения. предположим, что серийное облучение среды, зависящее от направления $\omega {\kern 1pt} *$, осуществляется прямоугольными импульсами длительностью δ вида

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

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

в которых величины c, d, δ и функции $h(\xi ,\omega {\kern 1pt} *,t)$, $H(\eta ,\omega {\kern 1pt} *)$ при $(\xi ,\omega {\kern 1pt} *,t) \in {{\Pi }_{{\omega {\kern 1pt} *}}} \times \Omega {\kern 1pt} *\; \times [0,T],$ (η, ω*) ∈ ∈ ${{\Pi }_{{\omega {\kern 1pt} *}}} \times \Omega {\kern 1pt} *$ заданы. В условии переопределения (5) мы требуем знания лишь осредненных значений плотности потока по промежутку, равному ширине импульса, смещенному на время пробега баллистической оставляющей зондируемого сигнала от источника к приемнику, что несколько снижает требования к временному разрешению детекторов.

3. ОЦЕНКА ВКЛАДА РАССЕЯННОГО ИЗЛУЧЕНИЯ. ПРОЦЕДУРА ЭКСТРАПОЛЯЦИИ

Обозначим через ${{L}_{{r,\omega }}}$ луч, исходящий из точки r ∈ ${{\mathbb{R}}^{3}}$ в направлении ω, ${{L}_{{r,\omega }}} = \{ r + \omega \tau :\tau > 0\} $, а через $d(r, - \omega )$ расстояние от точки $r \in G$ до плоскости Π_–ω* в направлении $ - \omega $, и пусть d(r, –ω, t) = = $\min \{ d(r, - \omega ),ct\} $.

Для решения начально-краевой задачи (1), (4) справедливо представление в виде равномерно сходящегося ряд Неймана [17–19]

где операторы $A$ и $S$ определяются выражениями

В представлении (6) функция

имеет смысл интенсивности нерассеянного поля, а функция ${{I}_{n}} = {{\left( {AS} \right)}^{n}}{{I}_{0}}$ при n = 1, 2, … описывает n-кратно рассеянное поле.

Можно показать, что для всех $\eta \in {{\Pi }_{{\omega {\kern 1pt} *}}}$ для функции I₁ справедлива оценка

при $d{\text{/}}c \leqslant t \leqslant d{\text{/}}c + \delta $, где $\overline \sigma = \mathop {\sup }\limits_{r \in G} \sigma $. Из (10) находим где $\varepsilon (\delta ) = c\delta {\text{/}}d$ и $C \leqslant \overline \sigma d{\text{/}}4$. При осреднении функции I₀ по промежутку времени $[d{\text{/}}c,d{\text{/}}c + \delta ]$ для всех $\eta \in {{\Pi }_{{\omega {\kern 1pt} *}}}$ получаем

Из неравенства (11) видно, что при уменьшении длительности зондирующего импульса осредненные значения функции I₁ стремятся к нулю. Можно показать, что осредненные по промежутку времени значения функций I₂, I₃, … также сходятся к нулю при δ → 0. Ограничимся приближением однократного рассеяния, оставляя в ряде Неймана только два первых слагаемых, и при проведении численных экспериментов мы покажем, что при зондировании импульсами небольшой длительности даже такое грубое приближение позволяет построить хорошую аппроксимацию решения уравнения переноса излучения в многократно рассеивающей среде.

Таким образом, из (10), (11), (13) в приближении однократного рассеяния получаем

Константу $C$ можно определить, например, если облучить среду двумя импульсами различной длительности δ₁, δ₂, 0 < δ₂ < δ₁, ε_I = ε(δ_i), а затем выразить интеграл от функции μ

Левая часть соотношения (15) является лучевым преобразованием Радона искомой функции μ(r). Таким образом, решение обратной задачи свелось к обращению преобразования Радона и линейных интегралов на множестве ${{\Pi }_{\omega }}_{{{\kern 1pt} *}} \times \Omega {\kern 1pt} *$ достаточно для реализации алгоритма обращения [20].

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

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

Для описания серийного облучения среды использовался импульсный источник, зависящий от направления зондирования ω*, с длительностью импульсов 100 и 300 пикосекунд (δ₁ = 300, δ₂ = 100).

Алгоритм тестировался в два этапа. На первом этапе задавались значения коэффициентов, описывающих среду и с помощью метода Монте-Карло [21] вычислялись значения функций $\Psi ({{\varepsilon }_{i}}) = H(\eta ,\omega {\kern 1pt} *,{{\varepsilon }_{i}})$. В экспериментах моделировался томограф, содержащий 101 детектор, а дискретизация по углу составляла 200 направлений. Для каждого детектора имитировалось 10⁵ траекторий фотонов, при этом отслеживалось до 10 актов взаимодействия фотона с веществом, что соответствует n = 10 в ряде Неймана (6) для решения уравнения переноса излучения. Выбранные параметры позволяют найти значения $H(\eta ,\omega {\kern 1pt} *,{{\varepsilon }_{i}})$ с относительной погрешностью, не превышающей 0.2%.

На втором этапе выполнялась экстраполяция данных о выходящем излучении и согласно формуле (15) находилось преобразование Радона функции μ. Также вычислялись значения преобразования Радона функции μ без экстраполяции

что соответствует формуле классической томографии без применения методов подавления рассеяния. После чего, применяя для обращения преобразования Радона алгоритм свертки [20], решалась задача восстановления коэффициента ослабления.

Традиционно значения коэффициента ослабления приводятся в единицах Хаунсфилда (HU), связанных со значениями коэффициента ослабления вещества μ следующим образом HU = = $\frac{{\mu - {{\mu }_{{water}}}}}{{{{\mu }_{{water}}} - {{\mu }_{{air}}}}} \times 1000$, где μ – коэффициент ослабления вещества, а ${{\mu }_{{water}}}$, ${{\mu }_{{air}}}$ – коэффициенты ослабления для воды и воздуха соответственно.

Фантом, выбранный для тестирования алгоритма, был предложен в работе [22] и представлял собой цилиндр высотой и диаметром в 10 см. Внутренний объем цилиндра поделен на пять цилиндрических слоев, каждый из которых играет свою роль при тестировании различных качеств алгоритмов реконструкции. При тестировании было выбрано сечение, приведенное на рис. 1a, которое проверяет возможности предлагаемого нами подхода для восстановления слабоконтрастных включений в случае, когда поблизости находятся сильно-рассеивающие объекты.

Среда включает в себя цилиндр, заполненный полимером с коэффициентом ослабления около 35 единиц Хаунсфилда. В свою очередь, этот цилиндр содержал еще три цилиндрических включения диаметром 9 мм со значениями коэффициента ослабления в –80, 0, и 80 HU. Цилиндрические включения имеют коэффициенты ослабления, близкие к фоновому значению, и представляют собой пример слабоконтрастных включений. Кроме того, фантом содержит 9 цилиндрических полостей диаметром 3 мм, расположенных в виде человеческой челюсти. Данные полости дают возможность устанавливать в них металлические штифты, играющие роль сильных рассеивателей, создавая тем самым плохие условия для реконструкции слабоконтрастных включений. В численных экспериментах мы ограничились случаем одного металлического включения, в то время как остальные полости оставались незаполненными и содержали воздух. Отметим, что подобный фантом достаточно тяжелый для реконструкции слабоконтрастных включений. Все дело в том, что он содержит в себе как материал с наименьшим значением в шкале Хаунсфилда (HU = = –1000), так и с максимальным (HU = 1000), что не позволяет применять для выделения слабоконтрастных включений, например методы сжатия динамического диапазона.

Обработка выходящего излучения с помощью формулы (15) проводилась для каждой проекции независимо. Пример такой обработки для одной проекции приведен на рис. 2. Анализ графиков показывает, что отношение баллистической составляющей к суммарному сигналу составляет порядка 50% для импульса длительности 100 пикосекунд, а для импульса длительностью 300 пикосекунд – менее 35%. Так как с увеличением длительности зондирующего импульса графики соответствующих проекций смещены вверх по оси ординат, то восстановленные по таким проекциям значения коэффициента ослабления будут заниженными. В свою очередь применение формулы (15) дает график проекции, проходящий существенно ближе к “идеальной” проекции, нежели графики проекций без обработки, и находится ниже кривой, соответствующей проекции без рассеяния. Для сравнения на рисунке приведена кривая, полученная путем линейной экстраполяции функции $\Psi (\varepsilon )$, т.е. в разложении (14) функция $\Phi (\varepsilon )$ заменена на ε. Линейная экстраполяция дает несколько худшие результаты и ее применение оправдано при малых значениях параметра ε, поэтому далее мы ее не использовали.

Результаты реконструкции фантома приведены на рис. 1б, в, г. Рисунок 1в содержит результаты восстановления фантома по данным о выходящем излучении, которые получены при зондирующих импульсах 100 и 300 пс. Увеличение длительности зондирующего импульса приводит к значительному снижению качества реконструкции. Даже визуально видно, что коэффициент поглощения восстановился с большой погрешностью. Более того центральное слабоконтрастное включение, имеющее значение коэффициента поглощения HU = 0, наиболее близкое к значению основной среды, стало практически не различимым. Рисунок 1г содержит результаты реконструкции, соответствующие данным, обработанным при помощи формулы (15).

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

Как видно из томограмм, влияние сильно-рассеивающего включения особенно сильно вблизи него и приводит к появлению артефактов, которые вызываются наличием сильно-рассеивающего включения и недостаточным уровнем дискретизации по угловой переменной. Для анализа данной ситуации при подсчете средних значений восстановленного коэффициента ослабления нами дополнительно выделялись два региона, отмеченные красными прямоугольниками на рис. 1a. Первый из них, обозначенный 1, находился в непосредственной близости от включений, вызывающих появление артефактов. Второй регион, обозначенный на рис. 1a как 2, напротив, находился далеко от рассеивающих включений.

Получившиеся в результате экспериментов значения приведены в табл. 1. В первой строке содержатся эталонные значения коэффициентов. В последующих строках приведены средние значения коэффициентов ослабления, полученные в результате реконструкции фантома с помощью формулы (16) при δ₁ = 300 и δ₂ = 100 и при реконструкции фантома по данным, полученным после обработки выходящего излучения по формуле (15). Из табл. 1 видно, что фильтрация рассеяния существенно улучшает качество реконструкции значений коэффициента поглощения и, как и ожидалось, при анализе экстраполяции одиночных проекций, приведенных на рис. 2, восстановление фантома по экстраполированным данным приводит к завышению реальных значений коэффициента ослабления. Однако для включения, содержащего металлический штифт (1000 HU), данные значения оказались почти в два раза меньше оригинальных значений. Мы связываем данное обстоятельство с малыми размерами самого включения и недостаточной дискретизацией по угловой переменной.

5. ЗАКЛЮЧЕНИЕ

В работе предложен приближенный экстраполяционный алгоритм нахождения коэффициента ослабления уравнения переноса излучения и была рассмотрена его простейшая реализация по результатам двух измерений. Увеличение количества измерений различных по длительности импульсов, прошедших через среду, повысит качество экстраполяции и приведет к повышению точности решения обратной задачи. В то же время проведение дополнительных измерений увеличивает радиационную нагрузку на изучаемый объект, усложняет устройство томографа и делает его более дорогостоящим. Тем не менее теоретические исследования в этой области весьма перспективны, и мы связываем с данной тематикой наши дальнейшие планы. Интересен вопрос выбора устойчивого к ошибкам измерений и оптимального с точки зрения точности многопараметрического разложения функции $\Psi (\varepsilon )$ в выражении (14).