CC BY
67
Вычислительные технологии
Том 17, № 6, 2012
Метод решения задачи томографии, основанный на специфике комптоновского рассеяния*
И. П. Яровенко Институт прикладной математики ДВО РАН, Дальневосточный федеральный университет, Владивосток, Россия e-mail: [email protected]
Предлагается метод решения задачи томографии, основанный на специфике комптоновского рассеяния. Постановка задачи формализуется с помощью методов теории переноса излучения. Особое внимание уделяется обсуждению практической применимости предлагаемого метода. Проводится ряд численных экспериментов, показывающих влияние несовершенства детекторов на точность работы метода.
Ключевые слова: теория переноса излучения, томография, комптоновское рассеяние, преобразование Радона.
Введение
Точность и эффективность алгоритмов томографического исследования среды напрямую зависит от того, какие эффекты взаимодействия излучения со средой учитываются в математической модели, используемой при построении алгоритма. Для большинства энергий, характерных для искусственных и природных источников рентгеновского излучения, преобладают три основных вида взаимодействия излучения с веществом — фотоэлектрическое поглощение, релеевское рассеяние и рассеяние по Комптону [1]. При этом если поглощение хорошо учитывается в рамках классической томографии, то рассеяние до сих пор остается основной проблемой, влияющей на качество реконструкции [2].
Классический метод решения данной проблемы заключается в коллимации исходного и принимаемого сигналов. Коллиматоры приводят к уменьшению случайных составляющих в детектируемом сигнале [2], однако при этом снижается чувствительность томографа и на порядок увеличивается время накопления необходимого объёма информации. Другой подход состоит в применении методов, позволяющих отфильтровывать рассеяние. Так, например, в монографии [3] при использовании в качестве источника излучения разрывной функции удалось восстановить неизвестный коэффициент ослабления в среде, где возможные типы взаимодействия излучения с веществом ограничиваются поглощением и релеевским рассеянием. Следует упомянуть о работах, использующих специфику комптоновского рассеяния. Как правило, это различные приближения однократного рассеяния, с помощью которых удаётся свести задачу томографии к восстановлению распределения плотности электронов по семейству интегралов по окружностям [4-6].
* Работа выполнена при финансовой поддержке конкурсного проекта ДВО РАН (проект 12-Ш-В-01М-004).
Настоящая работа посвящена исследованию метода решения задачи томографии, основанного на специфике комптоновского рассеяния. Идея решения задачи томографии заключается в том, что при комптоновском рассеянии энергия фотонов может только уменьшиться. Таким образом, при значении энергии, совпадающем с максимальной энергией источника излучения, интенсивность принимаемого сигнала должна описываться простым экспоненциальным законом (все рассеянные фотоны перейдут на более низкие уровни энергии). По сути данный подход аналогичен фильтрации рассеяния, используемой в [3], за исключением того, что не требует построения специальных источников излучения. В работе это формализуется на основе уравнения переноса излучения, а точнее его варианта, соответствующего случаю, когда в среде преобладает комптоновское рассеяние [7].
Особое внимание в представленном исследовании уделяется обсуждению практической применимости предлагаемого метода. Основная проблема здесь состоит в том, что для эффективной работы метода требуются детекторы с идеальным разрешением по энергии. Однако используемые в современных томографах сцинтилляционные детекторы в зависимости от конструкции детектора и энергии измеряемого сигнала позволяют определять энергию фотона с погрешностью от 2 до 20 % [2]. В статье приведены результаты численных экспериментов, показывающие влияние несовершенства детекторов на точность работы метода. Для этой цели предложенная формула реконструкции применяется к данным, полученным с помощью имитационного моделирования работы реального томографа. Для описания сигнала, принимаемого сцинтилляционными детекторами, использован метод Монте-Карло. Такой подход позволил имитировать сигнал, который будет наиболее близок к регистрируемому реальным томографом. Кроме того, метод Монте-Карло вносит в итоговый сигнал случайные шумы и ошибки, характерные для реального процесса, что более естественно при тестировании метода решения обратной задачи, чем использование каких-либо модельных шумов.
1. Прямая задача для уравнения переноса с чисто комптоновским рассеянием
Пусть процесс переноса излучения рассматривается внутри некоторой выпуклой ограниченной области С в трёхмерном евклидовом пространстве Е3. Введём следующие обозначения: г Е С — точка в Е3, ш', ш — направление движения фотона до и после рассеяния, ш', ш € П = {ш € Е3 : |ш| = 1}, а', а — энергия фотона до и после рассеяния, а', а Е I = [а, а]; 0 < а < а < то. Будем считать, что в процессе взаимодействия излучения с веществом фотоны могут рассеиваться только по закону Комптона. Это предположение приводит к тому, что при переходе в результате рассеяния фотона с характеристиками (ш, а) в фотон с характеристиками (ш', а') данные переменные связаны соотношением Комптона
а' = д(ш,ш',а), д(ш,ш',а)
1 + а(ш ■ ш' — 1)
где ш ■ ш' означает скалярное произведение векторов ш и ш', а переменная ш' принадлежит подмножеству единичной сферы Пш,а = {ш' : ш' € П, ш ■ ш' > 1 — 1/а + 1/а} и верны неравенства а < д(ш,ш', а) < а.
В качестве математической модели указанного процесса выбирается стационарное уравнение переноса [7], имеющее вид
ш • Уг f (т,ш,а) + г, а)/(т,ш,а) = к(т,ш,ш',а)/(т,ш',д(ш,ш',а))йш' + 3(т,ш,а). (2)
Здесь f (т,ш,а) — плотность излучения в точке т Е О, распространяющегося в направлении ш Е П и имеющего энергию а Е I; ^(т, а) — коэффициент полного взаимодействия излучения со средой в точке т при энергии а; к(т,ш,ш',а) — индикатриса рассеяния; 3(т,ш,а) — плотность внутренних источников излучения.
Для дальнейшего использования при формулировке рассматриваемых задач обозначим через = {т + шЬ,Ь > 0} луч, исходящий из точки т в направлении ш, ¿(т,ш) = тев\{Ьг,ш П О} — расстояние от точки т до границы области О в направлении ш. Здесь символом тев\ обозначена мера Лебега на прямой. Введём в рассмотрение множества Г± = {т Е дО : П О = 0}. В точках множества Г- (Г+) излучение в направлении ш входит в О (выходит из О). Обозначим также Г± = {(т,ш,Е) Е дО х П х I : т Е Г±}. Добавим к уравнению (2) следующее граничное условие:
f(£,ш,а) = к(£,ш,а), (£,ш,а) Е Г-, (3)
где функция к интерпретируется как плотность потока излучения, входящего в О.
Задача. Задачу определения функции f из уравнения (2) и условия (3) при известных ^,к,3,к будем называть прямой задачей (2), (3).
Введём обозначения: X = О х П х I, Сь(Х) — банахово пространство функций, определённых и ограниченных на Хо = X и непрерывных на Хо, с нормой
1Мк(Хо) = Эйр Мх)1
Сформулируем ограничения на коэффициенты уравнения (2). Будем считать, что ^ Е СЬ(О0 х1),к Е Сь(О0 хПхПх1). Относительно функции к предположим к Е Сь(Г-).
Введём также обозначение ^(т, ш,а) = ш • Уг f (т,ш, а) + ^(т, а)f (т, ш, а), где под выражением ш •У f (т, ш, а) будем понимать производную по пространственной переменной в направлении ш
дf (т + шЬ,ш, а)
ш • Уг f (т, ш, а) =
дЬ
4=0
Обозначим в виде
(Б<)(т,ш,а) = J к(т,ш,ш',а) <(т,ш',д(ш,ш',а))3,ш' (4)
линейный оператор, соответствующий интегралу столкновений, Б : Сь(Х0) ^ Сь(Х0). Определение. Функцию f (т,ш,а) назовём решением задачи (2), (3) если:
1) ^ Е Сь(Х0);
2) для всех (т,ш,а) Е Х0 она удовлетворяет уравнению ^ = Sf + 3;
3) выполняется граничное условие f (т — ¿(т, -ш)ш,ш,а) = к(т — ¿(т, -ш)ш,ш,а). При сделанных предположениях выполняются все ограничения работы [7] и справедлива следующая теорема.
п
Теорема 1. Решение краевой задачи (2), (3) существует и единственно.
Следует отметить, что теорема 1 не содержит традиционного неравенства для коэффициентов уравнения, характерного для многих подобных утверждений в теории переноса [3]. Наличие этого математического эффекта является следствием закона компто-новского рассеяния, устанавливающего функциональную зависимость между потерей энергии фотона при рассеянии и изменением направления его движения. Указанный математический эффект отсутствует в моноэнергетическом случае, где игнорируется потеря энергии при взаимодействии частиц со средой.
2. Постановка и решение обратной задачи томографии
Пусть дополнительно к краевому условию (3) задано условие
f (п,ш,а) = Н (п,ш,а), (п,ш,а)Г+. (5)
Функция Н(п, ш, а), определённая на Г+, задает выходящий из среды поток излучения. Сформулируем задачу томографии, которую мы будем изучать.
Задача томографии. Пусть задан некоторый уровень энергии а0. Определить функцию т,а0) из уравнения (2) и краевых условий (3), (5), если известны только функции Н(£,ш,а0), Н(£,ш,а0).
Для решения поставленной задачи введём дополнительные условия на функцию Н. Пусть функция Н, описывающая входящее излучение, отлична от нуля только для а < ао, т.е. в определении множества (1) будем иметь а = а0 В данном случае при приближении энергетической переменной а к величине а мера множества Пш,а, определяемого формулой (1), будет уменьшаться и при достижении величины а оно выродится в множество нулевой меры. В результате при а = а в уравнении (2) интегральное слагаемое будет отсутствовать.
Таким образом, учитывая, что в силу определения
Н(г + С(т, ш)ш — С(т + С(т, ш)ш, —ш)ш, ш, а) = Н(г — С(т, —ш)ш, ш, а)
и
(1(т+(1(г,ш)ш,—ш) с1(т,ш)
У ^(т + <С(т,ш)ш — шЬ,а)с<Ь = J ^(т + шЬ,а)<И
о -й(т-ш)
при а = а имеем
(й{т,ш)
— J ^(т + шЬ,а)с<Ь
-(1(т,—ш)
Следовательно, справедлива формула
с1(т,ш)
[ . ч, л Н (т — С(т, —ш)ш,ш,а)
ц(т + шг,а)М = Ь—^-,/ , ' / . (7)
] Н(т + <С(т, ш)ш, ш, а) 4 7
-<1(т,—ш)
В итоге рассматриваемая задача сведена к задаче обращения преобразования Радона от функции Как известно, эта задача имеет единственное решение в широком
классе функций [8]. Из единственности обращения преобразования Радона вытекает утверждение, обеспечивающее единственность решения задачи томографии.
Теорема 2. Пусть имеются две совокупности коэффициентов уравнения (2): {ц', к',7'} и {ц", к'', J"}, а /' и /'' — соответствующие решения прямых задач (2), (3) с одной и той же функцией Н. Пусть функция Н удовлетворяет условиям теоремы 1 и при а = а0 справедливо равенство /' = /'' на Г = Г- и Г+. Тогда ц' (г, а0) = ц'' (г, а0) почти всюду в О.
Доказательство данного утверждения почти полностью совпадает с соответствующим доказательством аналогичного утверждения, приведённого в [3].
Обсудим практическую применимость предложенного метода восстановления. Для работы формулы необходим источник излучения, скачком обращающийся в ноль при достижении некоторого значения энергии а. Математический пример функции Н, описывающей источник данного типа, привести несложно, однако на практике такой источник получить практически невозможно. В реально действующих рентгеновских трубках излучение генерируется в результате взаимодействия ускоренных электронов с атомами материала анода [2]. В процессе такого взаимодействия могут испускаться два типа рентгеновских квантов, один из которых — кванты характеристического излучения с энергиями, определяемыми энергетическими уровнями атомов анода. Они излучаются в результате взаимодействия ускоренных электронов с электронами атомной оболочки. Другой тип — это кванты тормозного излучения, которые возникают в результате взаимодействия падающих электронов с потенциалом ядра атома. Для тормозного излучения характерен непрерывный спектр. Суммарный спектр, испускаемый трубкой, является суммой характеристического и тормозного излучений. Качественное поведение спектра излучения рентгеновской трубки зависит преимущественно от ускоряющего анодного напряжения. На рис. 1 приведены характерные спектры от рентгеновской трубки с анодом из вольфрама (Ш) при различных ускоряющих напряжениях. Видно, что с увеличением ускоряющего напряжения возрастают интенсивность спектра и энергия рентгеновских квантов. Максимум спектра сдвигается в область высоких энергий. При напряжениях на трубке, превышающих порог возбуждения характеристических линий вольфрама (65, 67 кэВ), на спектре появляются пики характеристического излучения. Из приведённых на рисунке графиков следует, что невозможно выбрать подходящую энергию, при которой интенсивность скачкообразно обращается в ноль. Учиты-
а б в
Рис. 1. Характерные спектры от рентгеновской трубки с анодом из вольфрама при ускоряющих напряжениях 50 (а), 90 (б), 130 (в) кВ
вая, что разница между интенсивностями характеристического и тормозного спектров может достигать нескольких порядков, в качестве приближения такой энергии можно попытаться выбрать энергию, соответствующую характеристическому излучению.
Другим слабым звеном предлагаемого метода является предположение о том, что детекторы излучения имеют идеальное разрешение по энергии. Как уже указывалось, применяемые на практике сцинтилляционные детекторы в зависимости от конструкции детектора и энергии измеряемого сигнала позволяют определять энергию фотона с погрешностью от 2 до 20 %. В результате детектор будет регистрировать не только баллистические (прямолетящие) фотоны, но и фотоны, рассеявшиеся на небольшие углы, что также внесёт дополнительные искажения в томограмму [2].
Чтобы оценить влияние указанных факторов на качество реконструкции, в работе проведен ряд численных экспериментов по использованию предлагаемой формулы реконструкции к данным, полученным с помощью имитационного моделирования работы реального томографа.
3. Имитационная модель томографа
Будем считать, что множество детекторов образует кольцо с диаметром О и шириной Ь. Приведём описание алгоритма расчёта зарегистрированного томографом сигнала.
1. Розыгрыш начального направления и энергии фотона. На первом этапе необходимо определить начальное направление движения фотона, который излучает рентгеновская трубка. Направление движения разыгрывается равномерно в конусе с вершиной в точке расположения источника излучения и осью — прямой, соединяющей эту точку с центром кольца детектора. Раствор конуса задаётся параметрами томографа и, как правило, составляет угол порядка 50°. Первоначальная энергия фотона определяется как независимая реализация случайной величины, заданной с плотностью распределения, пропорциональной энергетическому спектру рентгеновской трубки. Спектры рентгеновских трубок рассчитывались с помощью алгоритмов, приведённых в [9].
2. Розыгрыш свободного пробега фотона в веществе. Свободный пробег разыгрывается с учётом экспоненциального закона уменьшения интенсивности потока фотонов с расстоянием. Для нахождения длины свободного пробега в кусочно-неоднородной среде применялся алгоритм, описанный в [10]. Если в результате свободного пробега фотон покинул рабочее пространство томографа, то такая траектория не отслеживается и выполняется переход к п. 1. В противном случае разыгрывается тип взаимодействия фотона с веществом.
3. Розыгрыш типа взаимодействия фотона с веществом. В данной работе ограничимся диапазоном энергий, при котором среди всех видов взаимодействия излучения с веществом преобладают фотоэлектрическое поглощение, комптоновское рассеяние и рассеяние по Релею [1]. Для определения типа взаимодействия излучения с веществом вычисляются вероятности поглощения фотона, рассеяния по Комптону и рассеяния по Релею:
Ра = 11а(Е )/ц(Е), Рс = 1 — Рк, Рк = цп(Е )/Ы(Е) + 1ю(Е)),
где ца — коэффициент поглощения, цс, — коэффициенты рассеяния по Комптону и Релею соответственно.
Далее символами в,1 будем обозначать независимые реализации случайной величины, равномерно распределенной в диапазоне [0,1].
Если в < Pa, то фотон поглотится, и такая траектория далее не отслеживается. Если в > Pa, то фотон рассеется и необходимо определить тип рассеяния. При y < Pr фотон рассеется по закону Релея, в противном случае — по закону Комптона.
4. Рассеяние фотона по Релею. Если фотон рассеялся по Релею, то разыгрываем новое направление движения ш' = , ш'2, ш'3), которое связано с первоначальным направлением ш следующими соотношениями [10]:
ш' = л/ 1 — V2 (cos фш3ш\ — sin фш2) !\Jl — ш"2 + РШ\, (8)
ш2 = //1 — V2 (cos фш3ш2 — sin фш\) ! -\Jl — ш"2 + vш2, (9)
ш3 = —/l — V2 cos ф\! 1 — ш2 + vш3, (10)
где V — косинус угла между направлениями ш и ш', распределённый с плотностью вероятности
3
g(V )=Ш (1 + V2),
азимутальный угол рассеяния ф = 2nY. Формулы (8)-(10) для нахождения нового направления движения справедливы при ш3 = ±1. В противном случае используются выражения
ш' = л/1 — V2 cos ф, ш'2 = л/1 — V2 sin ф, ш'3 = V. (11)
После нахождения нового направления возвращаемся к п. 2.
5. Рассеяние фотона по Комптону. В этом случае происходит изменение не только направления, но и энергии фотона. Вероятность того, что рассеянный фотон будет иметь энергию от а до а' в зависимости от значения случайного числа y, определяется из уравнения [10]
a í amin \
it da ( = Y. (12)
a \ a /
где amjn = а/(1 + 2а) — минимальная энергия, которую может приобрести фотон в результате рассеяния, dac/da — дифференциальное сечение комптоновского рассеяния в интервале энергий от а до а + da, определяемое комбинацией сечения Кляйна — Нишины — Тамма и функции некогерентного рассеяния [11], позволяющей учитывать влияние связи электронов в атоме вещества на процесс рассеяния.
В расчётах для розыгрыша энергии рассеянного фотона использовалась заранее подготовленная двумерная таблица, содержащая решения уравнения (12) с заданной точностью на равномерных сетках по а и y .В промежуточных точках применялась линейная интерполяция. Найденное значение энергии рассеянного фотона позволяет определить косинус угла рассеяния:
V
1 + 1/а — 1/а'.
Затем, разыгрывая равновероятно в диапазоне от 0 до 2п азимутальный угол рассеяния и применяя формулы (8)—(10), находим новое направление движения фотона. После этого возвращаемся к п. 2.
6. Регистрация фотона детектором. Будем считать, что в томографе применяется сцинтилляционный детектор, преобразующий энергию регистрируемой частицы
в световые вспышки, интенсивность которых пропорциональна энергии, потерянной частицей в сцинтилляторе. Если в результате трассировки фотон попадает в кольцо детекторов, мы продолжаем отслеживать его дальнейшую историю до тех пор, пока фотон либо не поглотится, либо не выйдет за пределы детектора. Пусть в результате такой трассировки фотона получены набор (xi,yi,zi), i = 1,..,N, координат точек, в которых фотон взаимодействовал с веществом детектора, и набор соответствующих значений энергии , потерянной фотоном в результате каждого взаимодействия. Номер детектора, который принял сигнал, определялся по координатам энергетического центра тяжести [12]:
N N N N N N
x = ^ ai, У = ^ ai, z = ^(аiZi)/^ а.
i=l i=l i=l i=l i=l i=l
Так как интенсивность световых вспышек в детекторе пропорциональна энергии, потерянной частицей в сцинтилляторе, сцинтилляционный детектор можно применять для определения энергии регистрируемой частицы. В данном исследовании регистрируемая детектором энергия моделировалась нормально распределенной случайной величиной с математическим ожиданием, равным энергии, которую потеряла частица, и дисперсией, определяемой разрешением детектора по энергии. Энергетическое разрешение детектора считалось пропорциональным квадратному корню из энергии [12], а коэффициент пропорциональности принимался одинаковым для всех уровней энергий и определялся из чувствительности детектора, характерной для конкретного томографа.
7. Оценка интенсивности зарегистрированного излучения. Пусть в результате трассировки p фотонов с энергией в интервале (а, а+Да], летящих в направлении, лежащем в конусе Дш, попали в детектор, расположенный в точке г и занимающий объём Дг. Тогда оценка интенсивности, зарегистрированной данным детектором, будет иметь вид
и, \ 1 ^ h(ai,Ui)
H(г,ш,а) > -:—т—т—,
v' ' ; N^ ДаДгДш
i=l
где ai,ui — начальные энергия и направление движения фотона (до взаимодействия со средой).
Алгоритм был реализован в виде компьютерной программы, которая использовалась в дальнейших экспериментах. При проведении экспериментов моделировался сканер, содержащий 39 детекторных колец диаметром 842 мм. Каждое кольцо содержало 642 детектора из ультрабыстрой керамики (UFC) размером 4 х 4 х 20 мм. Раствор конуса лучей рентгеновской трубки составлял 50°. Моделировалась рентгеновская трубка с анодом из вольфрама, оснащённая алюминиевым фильтром для подавления влияния фотонов с низкой энергией [9]. Данные параметры соответствуют томографу Biograph™ TruePoint PET • CT фирмы Сименс.
Во всех расчётах при моделировании траектории фотона отслеживалось до 10 актов взаимодействия со средой. Данные о сечениях взаимодействия излучения с веществом брались из таблиц Хабла — Зельтцера [13]. Для генерации случайных чисел, равномерно распределённых в диапазоне [0,1], использовался датчик, описанный в [14].
4. Численные эксперименты по решению обратной задачи
Проведённые в работе эксперименты можно разделить на две группы. Цель первой группы — показать влияние несовершенства источника и детектора излучения на отклонение сигнала от идеального, рассчитанного по формуле (6). Для этой цели используется достаточно простой фантом, представляющий собой цилиндр диаметром 20 мм и высотой 300 мм, заполненный водой, ось которого совпадает с осью кольца детектора. Излучение, рассчитанное с помощью метода Монте-Карло, сравнивалось с излучением, рассчитанным при помощи простого экспоненциального закона ослабления интенсивности.
При проведении этой группы экспериментов интенсивность зарегистрированного томографом сигнала вычислялась при двух значениях ускоряющего напряжения на рентгеновской трубке. В первом случае ускоряющее напряжение было сравнительно невелико (70 кВ) и спектр излучения рентгеновской трубки состоял только из фотонов тормозного излучения. Как следует из рис. 2, а, график генерируемого источником спектра достаточно быстро убывает при увеличении энергии, и данный случай можно рассматривать как приближение скачкообразного источника. Такое поведение спектра вызвано наличием алюминиевого фильтра [9]. В качестве уровня энергии, на котором регистрируется сигнал, было выбрано 78 кэВ. На рис. 2, б, в приведено сравнение поведения графиков проекций, полученных методом Монте-Карло и рассчитанных по закону экспоненциального ослабления (6). Видно, что даже при погрешности детектора в 8% качественное поведение наблюдается достаточно хорошо, хотя и имеется достаточно большой разброс данных, вызванный регистрацией рассеянных фотонов. При уменьшении погрешности детектора до 2 % расчётный и моделируемый сигналы фактически совпадают (см. рис. 2, в). В данном случае можно ожидать хороших результатов восстановления исследуемой среды.
При увеличении ускоряющего напряжения до 120 кВ в спектре излучения появляются пики, соответствующие характеристическому излучению (рис. 3, а). Сигнал регистрировался на уровне энергии 67 кэВ, что соответствует второму пику характеристического излучения. Данная серия тестов была проведена для того, чтобы продемон-
Рис. 2. Результаты моделирования работы томографа при ускоряющем напряжении 70 кВ: а — график спектра, создаваемого рентгеновской трубкой; б, в — сравнительные графики томографических проекций при чувствительности детектора 8 и 2 % соответственно. Точки — проекции, полученные методом Монте-Карло, линии — проекции, рассчитанные по закону экспоненциального ослабления (6)
стрировать, к чему приведёт выбор в качестве энергии просвечивания одного из пиков характеристического излучения.
На рис. 3, б представлены результаты сравнительного моделирования в случае погрешности детектора, составляющей 8%. Хотя разница между интенсивностями характеристического и тормозного спектров достаточно велика, наличие "хвоста" тормозного излучения приводит к тому, что на выбранный энергетический уровень попадает достаточно большое количество фотонов, изначально имеющих другую энергию. В результате наблюдается большое различие между графиком, построенным с помощью метода Монте-Карло, и теоретической кривой, полученной по закону экспоненциального ослабления, что в итоге снижает качество реконструкции.
Цель второй группы экспериментов — выяснить влияние несовершенства детектора и источника излучения на пространственное разрешение метода. Для этого используется широко известный фантом Дерензо. Эксперименты данной группы состояли из двух этапов. На первом методом Монте-Карло имитировался сигнал, регистрируемый детекторами томографа, на втором — методом свёртки и обратной проекции [8] восстанавливалась функция / и строилась томограмма.
Рис. 3. Результаты моделирования работы томографа при ускоряющем напряжении 120 кВ: а — график спектра, создаваемого рентгеновской трубкой; б — сравнительный график томографических проекций при чувствительности детектора 8%. Точки — проекции, полученные методом Монте-Карло, линии — рассчитанные по закону экспоненциального ослабления
аб
Рис. 4. Результаты восстановления структуры фантома Дерензо при ускоряющем напряжении на рентгеновской трубке 70 кВ: а — схематичное изображение структуры фантома, б — результаты восстановления коэффициента полного ослабления на уровне энергии 78 кэВ
Результаты моделирования фантома Дерензо приведены на рис. 4. Рис. 4, а иллюстрирует структуру фантома, который представляет собой цилиндр из оргстекла диаметром 20 см, условно разбитый на шесть секторов, в каждом из которых имеются отверстия определённых диаметров: 6, 5, 4, 3.5, 3 и 2.5 мм. Отверстия заполнены водой. Данный фантом, как правило, применяется для определения степени разрешения алгоритма реконструкции. Реконструкция фантома Дерензо, полученная по результатам просвечивания источником излучения с ускоряющим напряжением 70 кВ, приведена на рис. 4, б. Как видно, включения диаметром более 3 мм достаточно хорошо различимы. Это свидетельствует о том, что предложенный метод реконструкции имеет пространственное разрешение не хуже, чем применяемые в томографах алгоритмы [2].
Список литературы
[1] Фано У., Спенсер Л., Бергер М. Перенос гамма-излучения. М.: Госатомиздат, 1963.
[2] Cierniak R. X-Ray Computed Tomography in Biomedical Engineering. London: Springer, 2011. 319 p.
[3] Аниконов Д.С., Прохоров И.В., КовтАнюк А.Е. Использование уравнения переноса в томографии. М.: Логос, 2000. 224 с.
[4] Truong T.T., Nguyen M.K., Zaidi H. The mathematical foundations of 3d compton scatter emission imaging // Intern. J. of Biomedical Imaging. 2007. Article ID 92780.
[5] Cruvinel P.E., Balogun F.A. Compton scattering tomography for agricultural measurements // Eng. Agric. 2006. Vol. 26, No. 1. P. 151-160.
[6] Palamodov V.P. An analytic reconstruction for the Compton scattering tomography in a plane // Inverse Problems. 2011. Vol. 27, No. 12. P. 125004.
[7] Аниконов Д.С., Коновалова Д.С. Краевая задача для уравнения переноса с чисто комптоновским рассеянием // Сибирский матем. журнал. 2005. Т. 46, № 1. C. 3-16.
[8] Наттерер Ф. Математические аспекты компьютерной томографии. Пер. с англ. М.: Мир, 1990. 288 с.
[9] Boone J.M., Seibert J.A. An accurate method for computer-generating tungsten anode x-ray spectra from 30 to 140 kV // Medical Phys. 1997. Vol. 24, No. 11. P. 1661-1670.
[10] Михайлов Г.А. Весовые методы Монте-Карло. Новосибирск: Изд-во СО РАН, 2000.
[11] Hubbell J.H., Veigele W.J., Briggs E.A. et al. Atomic form factors, incoherent scattering functions, and photon scattering cross sections //J. Phys. Chem. Ref. Data. 1975. Vol. 4. P. 471-538.
[12] Manno I. Event Reconstruction Strategies. ([Online] Available: http://www.rmki. kfki.hu/^ manno/Borexino.html)
[13] Hubbell J.H., Seltzer S.M. Tables of X-Ray Mass Attenuation Coefficients and Mass Energy-Absorption Coefficients 1 keV to 20 MeV for Elements Z = 1 to 92 and 48 Additional Substances of Dosimetric Interest. Gaithersburg, 1995 (Preprint National Institute of Standard and Technology).
[14] L'Ecuyer P., Cote S. Implementing a random number package with splitting facilities // ACM Trans. on Math. Software. 1991. Vol. 17, No. 1. P. 98-111.
Поступила в 'редакцию 20 августа 2012 г.
