CC BY
57
Вычислительные технологии
Том 16, № 6, 2011
Моделирование процесса измерения комптоновского рассеяния в позитронной эмиссионной томографии*
И. Г. Казанцев1 , И. П. Яровенко2, И. В. Прохоров2 1 Институт, вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия 2Институт прикладной математики ДВО РАН, Владивосток, Россия e-mail: kazantsev.ivan6@gmail.com, yarovenko@iam.dvo.ru, prh@iam.dvo.ru
Исследуется идеализированная интегральная математическая модель формирования проекционных данных в позитронной эмиссионной томографии на основе рассеянных фотонов. Проведенные численные эксперименты показали хорошее совпадение проекций, полученных аналитически и статистическим методом Монте-Карло.
Ключевые слова: позитронная эмиссионная томография, комптоновское рассеяние, идеализированная модель, статистическое моделирование.
Введение
Идеализированная интегральная модель позитронной эмиссионной томографии (ПЭТ) на первичных (не рассеянных) фотонах хорошо известна [1]. Для распределения внутренних источников активности изотопа f (x,y,z) внутри среды с линейным коэффициентом ослабления ^(x, y, z) и детекторов малых размеров А и B вблизи поверхности просвечиваемого тела модель формирования данных имеет вид
где ¿1' — элементы дайны на линии ЛБ. Функция ^ описывает анатомическую структуру просвечиваемого организма, и ее носитель Б(^) есть ограниченная одноевязная область. Источники с эмиссионной активностью / после их введения в среду Б(^) распределяются в соответствии с особенностями функционирования организма, в силу чего носителем функции / является некоторая многосвязная область Б(/) С Б(^) неизвестной формы. Задача состоит в определении (реконструкции) активности изотопов / и носителя Б(/) по данным РЛБ, регистрируемым большим множеством пар детекторов (Л, Б). В традиционной модели ПЭТ экспоненциальный множитель в (1) считается известным из данных предварительного сканирования (например, методами классической рентгеновской или магнитно-резонансной томографии). Физическая модель ПЭТ основана на использовании большого числа пар гамма-квантов (и, у), разлетающихся
* Работа выполнена при финансовой поддержке РФФИ (гранты № 10-07-00131-аи 11-01-98521), Федеральной целевой программы "Научные и научно-педагогические кадры инновационной России" и гранта 09-П-СО-004 Конкурса интеграционных проектов ДВО РАН с научными учреждениями СО РАН.
B
B
(1)
A
A
приблизительно в противоположных направлениях из точки аннигиляции С (рис, 1, о) в результате столкновения позитрона, излучаемого радионуклидом /, с одним из электронов среды Эти два фотона регистрируются методом совпадений внутри малого
временного промежутка. Данные детекторов РАВ корректируются и затем интерпре-
/
вичные фотоны с энергией Е = 511 кэВ,
Однако на практике в детекторы попадают и фотоны V с энергия ми Е' < Е, претерпевшие комптоновекое рассеяние (рис, 1,6). В зависимости от разрешения детектора по энергии в простейшем случае используется спектральное окно чувствительности Ш = [£, 511] с порогом £ и фотоны с энергией вне окна Ш в алгоритмах ПЭТ не учитываются. Для попавших в окно фотонов методами моделирования оценивается отношение числа первичных фотонов к рассеянным и данные масштабируются в соответствии с этим отношением, что принято называть коррекцией на рассеяние. Выбор порога £ и другие проблемы (одна из них состоит в малой статистике отсчетов) обусловливают широкое привлечение в ПЭТ статистических методов. Модель (1) является идеализированной, однако именно на ее основе получены аналитические формулы
обращения, позволившие ответить на вопросы о принципиальной возможности восста-
/
поетроены операторы прямой задачи, позволяющие применять итерационные методы. Из физических экспериментов известно, что в медицинских сканерах рассеянные фотоны составляют от 30 до 70 % общего числа регистрируемых фотонов, при этом однократно рассеянные (вторичные) фотоны составляют большинство (90-95%) из всех многократно рассеянных фотонов. За последние два десятилетия достигнут прогресс в создании детекторов с достаточно высоким (до 2-3 %) разрешением по энергии. Поэтому актуальной становится задача не только коррекции, или отбрасывания рассеянных фотонов, но и использования их в процессе восстановления активности /(х,у,г) внутренних источников. Для этого необходимо иметь инструмент моделирования регистрации потока рассеянных фотонов, подобный проекционной модели (1) классической томографии. Общие проблемы переноса гамма-излучения и его применения в томографии исследованы в монографиях [2-4],
Рис. 1. Идеализированные модели ПЭТ: а — ПЭТ на первичных фотонах (u, v) с энергией E = 511 кэВ, C € D(f) — точка аннигиляции; б— ПЭТ на однократном комптоновском рассеянии; (u, v) — пара аннигиляционных фотонов, S € D(p) — точка рассея ния, v' — фот он v энергии E', рассеянный с углом 0; точка Р геометрически является потенциальной точкой рассеяния с углом 0, однако находится вне среды у (Р / D(p), ¡л(Р) = 0) и вклада в значения счетчика B не дает; в — поверхность — геометрическое место точек S рассеяния с углом 0 и сферическими координатами (ф, | AS|); точкa Cxy — проекция C на плоскость xAy, ZABS = 0 — ф
В практических задачах оценки однократного рассеяния в ПЭТ в биомедицинских сканерах широко используется аппроксимация SSS (Single Scatter Simulation) [5, 6]. Вывод этой аппроксимации и исследование ее связи с уравнением переноса можно найти в работе [7]. Приближение SSS оценивает количество вторичных комптоновских анни-гиляционных пар фотонов, зарегистрированных парой детекторов A и B как интеграл по всему объему V возможных событий рассеяния S за единицу времени:
ЧАВ - ! dv( aAS(JBS \ У д(Тс\тА | ТВЛ (2)
Ьу J dV \An\AS\2\BS\2 J дп [ +1 1 {2)
V
где
(s в \ S (s в \ B
— I^dl+S v'dl) Г ' — Jvdl+J v'dl) Г
\A S J / fdl, IB = eaae. V S J
IA = tAs 4s e S ' J fdl, IB = tAs tBs e \A S ' Jfdl. (3)
A S
Здесь gaS и aBS — геометрические поперечные сечения детекторов A и B; f — активность источников; a = a(E,S) — линейный коэффициент ослабления, зависящий от
энергии E и положения точки рассеяния 5; ^ - дифференциальное сечение компто-
ÔU
новского рассеяния; tAS и eBS — величины, связанные с эффективностью детекторов A и B, Эффективность е поглощения энергии приемником зависит как от энергии кванта, так и от угла падения фотона на плоскость детектора. Знак штриха означает значение величины после рассеяния. Уравнение (2) симметрично относительно рассеян-
ные фотоны регистрируются обоими детекторами Дифференциальное сечение
комптоновского рассеяния вычисляется по формуле
dac r\( E (E\2 . 2л (E\ Л где re — классический радиус электрона. Полное сечение рассеяния имеет вид
(5)
Модифицируем это уравнение в предположении, что система (A, B) настроена для ре-гиетрирации пар фотонов (u,v), из которых u попадает в детектор A нерассеянным с энергией E = 511 кэВ, а фотон v', отклонившись в точке S на угол 9, зарегистрирован в B с энергией E' (см. рис, 1, б). Связь между энергиями E и E' до и после столкновения
9
E'/E =1/(2 - cos 9). (6)
IB
fs в \ S a ÔOr - f^dl+f M'dl r
Qf^(x,y,z) = JL^eAse'BSe V * Jjfdl, (7)
CA
где (x, y, z) — декартовы координаты точки рассеяния S. Тогда (2) принимает вид
(8)
Это уравнение будет использовано для вывода идеализированного интегрального урав-
9
1. Идеализированная модель
Аппроксимация (2)-(3) апробирована в натурных физических и модельных статистических экспериментах по коррекции на рассеяние и известна [8-10] как достаточно точный инструмент описания процесса регистрации рассеянных фотонов в медицинских сканерах, адаптируемый к объектам V различной формы. Однако для детекторов с идеальным спектральным разрешением и с привлечением геометрических свойств траекторий рассеянных фотонов и поверхностей, на которых сосредоточены точки рассеяния с одинаковым углом в, в уравнении (8) возможно точное описание объема V, и тогда интегрирование в (8) приобретает конкретные пределы, что позволяет вывести интегральное уравнение однократного рассеяния.
Из геометрии однократного рассеяния легко видеть (см, рис, 1,6), что для всех точек рассеяния S с углом рассеяния в справедливо простое геометрическое свойство: ZASB = п — в. Тогда геометрическим местом точек, где происходит рассеяние с определенным углом в, является поверхность Ee тела вращепня Ve, образованного вращением вокруг оси Z дуги ASB (см. рис. 1, в). Используя сферические координаты (ф,^,г) с началом в точке A, поверхно сть Ee и тел о Ve могут быть записаны как
Se = {(ф,р,г)| ф Е [0, 2п), ф Е [0, в], r = |AS|}, (9)
V, = {(ф,ф,г)| ф е [0,2п),ф е [0,9],Г е [0, |АБ|]}. (10)
Легко показать, что радиус окружности, которой принадлежит дуга АББ,
Я- С1 "Л
Используя закон синусов для треугольника ДАББ
_ \ВБ\ _ вт(9 — ф) вт(ф)
получаем выражения для хорд
2R, (12)
И* - ич^Б* - ич£* <">
Б
ми, что 9 < Т для некоторого выбранного угла Т < п/2. Тогда область интегрирования Ут определена формулой (10) и соответствует в (8) носителю V точек рассеяния с уг-9<Т
Идеализированная модель БАВ рассеяния в точках множества Ут будет выведена, используя интеграл (8), вычисленный в точках малых детекторов круглой формы А и В с диаметрами 8 и центрами в А и Б соответственно и усредненный затем по паре кругов (А, В) и их площадям в виде предела
= Ьта jdAídBí** Щ^щ) «'•'<*•»•*>• (14)
7г— а в VT
4
Для малых значений 6 геометрические сечения кругов А и B, перпендикулярные к AS SB
п62 cos p п62 cos^ — p) , vas ~-^-, ctbs ~-^-• (15)
Подставляя (15) в (14) и полагая значения внутренних интегралов J близкими при
VT
варьировании точек A Е А и B Е B, получим
VT
Произведем замену переменных в (16) с декартовых координат (x,y,z) на криволинейные (ф,р,в), где (ф, p) — сферические координаты, а в Е [0,T] (при фиксированном расстоянии |AB|), по правилу
x = X(ф, р,в) = | AS| sin p cos ф,
y = Y(ф, p, в) = | AS| sin p sin ф,
z = Z(ф,р,в) = | AS| cos p, (17)
здесь зависимость |AS| от в ж p выражается уравнением (13), Вычислим элементарный объем
dVT = dxdydz = | J (18)
где J есть якобиан
" дХ/дф дХ/др дХ/дв " J = ду/дф дУ/др дУ/дв , (19)
_ дZ/дф дZ/дp дZ/дв _
| J|
| АЕ> |3 sin2(p) sin2(б1 — р) |Л9|2|£Д|2 1 J| = d6t J =-^-= |АВ| ' (20)
и после замены переменных уравнение (16) принимает вид
t 2п e
Stf = /<«/#/# «/„.«.,Í). (21)
0 0 0
Учитывая (7) и полагая точечные приемники идеально эффективными (е = е' = 1), получим интегральное уравнение рассеяния с углами, меньшими T:
те 2п ( s в \
5- = ¡M¡ j i (Г+í "'"J х
0 0 0
|AS|
х f (ф, p, r)dr. (22)
т
Подынтегральная для J часть
о
в 2п f S B \ |AS|
¿ab [ , cos^cos^-e) Г ^,tp,\AS\)dac -{l^+l^j Г Ь = J J W-^-Же v J J №,<P,r)dr (23)
о о о
может интерпретироваться как мгновенное количество рассеянных с углом в фотонов, при этом
т
SATB = j tfBde. (24)
о
Предварительные результаты, подобные формуле (23), с точностью до скалярного множителя были получены в работе [11] с использованием геометрических вероятностей.
2. Статистическое моделирование
Для проверки адекватности формулы (23) реальному физическому процессу используем имитационное моделирование на основе метода Монте-Карло [12]. Идея метода основана на представлении переноса излучения в среде в виде случайного процесса — движения фотонов через среду, используя известные законы взаимодействия излучения с веществом, В результате такого моделировния накапливается необходимое количество статистической информации для того, чтобы построить оценку величины ^^. Этот подход позволяет имитировать сигнал, который будет наиболее близок к тому, что регистрирует реальный томограф, так как вносит в итоговый сигнал случайные шумы и ошибки, характерные для реального физического процесса.
Алгоритм, используемый в работе, реализуется в следующей последовательности,
1, Розыгрыш точки аннигиляции фотона. Пусть источник с эмиссионной активностью / распределен в некоторой области С0 = ) с О(^). На первом этапе разыгрываем точку г образования позитрона равномерно внутри О0. Будем считать, что позитрон тут же аннигилирует, так что точки образования и аннигиляции позитрона совпадают,
2, Розыгрыш направления движения фотонов, образовавшихся, в результате аннигиляции. Ограничимся рассмотрением наиболее вероятного случая, когда в результате аннигиляции образуются два фотона. Для применяемых в медицине фармпрепаратов отклонение образовавшихся в результате аннигиляции фотонов от антиколлинеарности, как правило, не превышает 10-3 рад, [5], В связи с этим будем считать, что ан-нигиляционные фотоны распространяются строго в диаметрально противоположных направлениях. Вследствие изотропии первоначальных направлений движения фотона равномерно на сфере разыгрываем направление ш = ш3), в котором полетит первый из образовавшихся фотонов. Далее отслеживаем историю фотона, движущегося в направлении ш.
3, Розыгрыш свободного пробега фотона, в веществе. Свободный пробег разыгрывается с учетом экспоненциального закона уменьшения интенсивности потока фотонов с расстоянием. Для нахождения длины свободного пробега в кусочно-неоднородной среде применялся алгоритм, приведенный в [12], Если в результате свободного пробега фотон
покинул томограф, то такая траектория отбрасывается и выполняется переход к п, 1, В противном случае разыгрывается тип взаимодействия фотона с веществом,
4, Розыгрыш взаимодействия фотона с веществом. Для характерного в позитрон-ной эмиссионной томографии диапазона энергий среди всех видов взаимодействия излучения с веществом преобладают фотоэлектрическое поглощение, комптоновское рассеяние и рассеяние по Рэлею [2], Для определения типа взаимодействия излучения с веществом вычисляются вероятности поглощения фотона, рассеяния по Рэлею и по Комп-тону
Ра = Рй = ^(Е)+ ^с(Е)), Рс = 1 - Ръ (25)
где — коэффициент поглощения, ^с _ коэффициенты рассеяния по Рэлею и Комптону соответственно. Здесь и далее символами а, в будем обозначать независимые реализации случайной величины, равномерно распределенной в диапазоне [0,1],
Если а < Ра, то фотон поглотился, и такая траектория дальше не отслеживается. Если а > Ра, то фотон рассеется, и необходимо определить тип рассеяния. При в < Рл фотон рассеется по закону Рэлея, в противном случае фотон рассеется по закону Комп-тона,
5, Рассеяние фотона по Рэлею. Если фотон рассеялся по закону Рэлея, то разыгрываем новое направление движения ш' = (ш', ш2 , ш3) которое связано с первоначальным направлением ш следующими соотношениями:
ш
ш.
[ = Vi — V2 (cos фш3^1 — sin фшъ) j\Jl — Ш2 + V OJi, (26)
2 = л/1 — У2 (cos фшз U)2 — sin фш\) j\j\ — U)\ + УШ2, (27)
ш'3 = —л/l — У2 cos Ф \Jl — Ш3 + VLOs, (28)
где v — косинус угла между направлениями шиш', распределенный с плотностью вероятности
9П = -^(1-и2), (29)
азимутальный угол рассеяния ф = 2пв- Формулы (26)-(28) для нахождения нового направления справедливы при ш3 = ±1, В противном случае используются формулы
= л/ 1 — У2 COS ф, Ш2 = л/ 1 — У2 SÍn0, Шд = и. (30)
После нахождения нового направления движения возвращаемся назад к п, 3,
6, Рассеяние фотона по Комптону. При рассеянии фотона по закону Комптона происходит изменение не только направления, но и энергии фотона. Вероятность того, чтобы рассеянный фотон имел энергию от E до E', в зависимости от значения случайного числа а определяется из уравнения
E' / Emin \ 1
/Ы/Ы =«■ ^
E \ E )
где Em¡n = E/(1 + 2E/E0) — минимальная энергия, которую может приобрести фотон в результате рассеяния, dac/dE — дифференциальное сечение комптоновекого рассеяния
в интервале энергий от Е до Е + йЕ, определяемое комбинацией сечения Кляйна — Нишины — Тамма и функции некогерентного рассеяния [2], позволяющей учитывать влияние связи электронов в атоме вещества на процесс рассеяния.
Решение уравнения (31) — достаточно трудоемкий процесс. Для ускорения в расчетах использовалась заранее подготовленная двумерная таблица, содержащая решения уравнения (31) с заданной точностью на равномерных сетках по Е и а. В промежуточных точках применялась линейная интерполяция. Найденное значение энергии рассеянного фотона позволяет определить косинус угла рассеяния:
где Е0 — энергия покоя электрона. Далее, разыгрывая азимутальный угол рассеяния равновероятно в диапазоне от 0 до 2п и применяя формулы (26)—(30), находим новое направление движения фотона. После этого возвращаемся к п, 3,
7. Регистрация фотона детектором. Если в результате трассировки фотон попал в детектор, запоминаем номер детектора и переходим к отслеживанию фотона в направлении —ш. В случае, когда оба фотона, образовавшиеся в результате аннигиляции, попали соответственно в детекторы А и В и энергия рассеянного фотона находится в заданном диапазоне [а, а + Да], зависящем от в, оценка функции ^^ увеличивается на величину интенсивности активного вещества в точке аннигиляции позитрона. Итоговая
^ АВ
оценка величины имеет вид
где N — общее число моделируемых частиц, а суммирование ведется по всем событиям парного попадания частиц в детекторы.
Данный алгоритм был реализован в виде компьютерной программы, которая использовалась в численных экспериментах. Для генерации случайных чисел, равномерно распределенных в диапазоне [0,1], использовался датчик, описанный в [13], Для проверки формулы (23) были проведены численные эксперименты с фантомом, широко используемым в литературе [14] и позволяющим геометрическое моделирование с контролем границ объекта, вне которого рассеяние считается несущественным. Данный фантом (рис, 2) представляет собой цилиндр радиусом 8 см, заполненный водой.
V = 1 + Ео/Е — Ео/Е',
(32)
(33)
г=1
а
б
А
О
2
Рис. 2. Томографический фантом, сечение уОг (а), общий вид (б); проекции и вычисляются для 100 пар приемников (А, В) с координатами центров = хВ = 0, уА = = 8.0 — (г — 0.5)5, гА = 0, гр = 16, г = 1,..., 100, 5 = 0.16 см, ^ = у! = 0.096 см-1
Рис. 3. Сравнительные результаты моделирования проекций комптоновского рассеяния на томографическом фантоме; пунктирные линии — проекции Мд, рассчитанные методом Монте-Карло, сплошные — проекции проекции от сферы № 1 (а), №2 (б), №3 (в), №4 (г)
Внутрь цилиндра помещены четыре стационарно излучающих сферических источника диаметрами 1.8 см, заполненные фармпрепаратом с единичной активностью: № 1 в центре и №2, 3, 4, смещенные на 4 см от центра. Фотоны регистрируются линейками детекторов A и B, лежащими в плоскости yOz. Считается, что детекторы линейки A регистрируют первичные фотоны, линейки B — рассеянные. Каждая линейка состоит из 100 детекторов сечения 6 = 0.16 см.
Во всех расчетах при моделировании траектории движения фотона отслеживалось до 10 актов взаимодействия со средой. Данные о сечениях взаимодействия излучения с веществом брались из таблиц Хабла [15]. Моделировалось 1011 траекторий (время расчета па компьютере Pentium 3.2 GHz около 5 суток). При расчете профилей сигнала по формуле (23) коэффициент ослабления в воде принимался одинаковым дня всех энергий ^ = ^ = 0.096 см-1. Проекции и Мв масштабировались умножением на скаляр по принципу совмещения точек максимума и приводились к условным единицам в интервале [0,10].
Графики па рис. 3 демонстрируют поведение сравниваемых величин в случае, когда отслеживается рассеяние на угол в = п/3. Наблюдается хорошее соответствие между проекциями, полученными по формуле (23), и статистическим моделированием. Небольшое различие, обусловленное многократным рассеянием фотонов, наблюдается лишь по краям графиков. Аналогичная ситуация наблюдается и дня кривых, соответствующих другим углам рассеяния.
Заключение
В качестве исходной в работе использовалась известная модель регистрации однократно рассеянных фотонов в позитронной эмиссионной томографии, представляющая собой интеграл по области событий рассеяния фотонов. Однако предметом восстановления в эмиссионной томографии являются не точки рассеяния фотонов, а области концентрации активности изотопов. Получено интегральное преобразование, позволяющее для точечных детекторов, обладающих идеальным спектральным разрешением, оценивать комптоновское рассеяние в виде интеграла от распределений коэффициента ослабления и активности радионуклида. Это дает возможность сформулировать прямую задачу в детерминистской постановке в виде оператора именно от тех величин, которые подлежат реконструкции. Для тестового томографического фантома и геометрии сбора данных в виде линеек приемников проведено сравнительное статистическое моделирование регистрации комптоновского рассеяния. Получено хорошее совпадение проекций, генерированных методом Монте-Карло и вычисленных по полученной формуле интегрального преобразования.
Список литературы
[1] Терещенко с.А. Методы вычислительной томографии. М.: Физматлит, 2004. 320 с.
[2] Фано У., Спенсер Л., Бергер М. Перенос гамма-излучения. М.: Госатомиздат, 1963. 284 с.
[3] Аниконов Д.С., Ковтанюк А.Е., Прохоров И.В. Использование уравнения переноса в томографии. М.: Логос, 2000. 224 с.
[4] Горшков В.А., Крёнинг \!.. Аносов Ю.В., Доржгочоо О. Томография на рассеянном излучении. М.: Технополиграфцентр, 2002. 146 с.
[5] Ollinger J.M. Model-based scatter correction for fully 3D PET // Physics in Medicine and Biology. 1996. Vol. 41. P. 153-176.
[61 Watson C.C., Newport D.A., Casey M.E. A single scatter simulation technique for scatter correction in 3D PET // Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine / Eds. P. Grangeat and J.L. Amans. Amsterdam, Netherlands: Kluwer Acad., 1996. P. 255-268.
[7] kosters T. Derivation and Analysis of Scatter Correction Algorithms for Quantitaive Positron Emission Tomography. University of Miinster, 2010. PhD Thesis.
[81 Watson C.C. New, faster, image-based scatter correction for 3D PET // IEEE Trans. Nucl. Sci. 2000. Vol. 47. P. 1587-1594.
[91 Werling A., Bublitz O., Doll J. et al. Fast implementation of the single scatter simulation algorithm and its use in iterative image reconstruction of PET data // Physics in Medicine and Biology. 2002. Vol. 47. P. 2947-2960.
[101 Accorsi R., Adam L.E., Werner M.E, Karp J.S. Optimization of a fully 3D single scatter simulation algorithm for 3D PET // Ibid. 2004. Vol. 49. P. 2577-2598.
[Ill Kazantsev I.G., Matej S., Lewitt R.M. Geometric model of single scatter in PET // IEEE Medical Imaging Conf. San Diego, 2006. P. Mil 313.
[121 Ермаков C.M., Михайлов Г.А. Статистическое моделирование. М.: Наука, 1982. 296 с.
[13] L'Ecuyer P.L., Cote S. Implementing a random number package with splitting facilities // ACM Trans, on Math. Software. 1991. Vol. 17. P. 98-111.
[14] Riauka T., Gortel W. Photon propagation and detection in single-photon emission computed tomography — an analytical approach // Medical Phvs. 1994. Vol. 21. P. 1311-3121.
[15] Hubbell J.H., Seltzer S.M. Tables of X-Rav 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 // Technical Report NISTIR 5632. National Inst, of Standards and Technology, USA, 1995.
Поступила в редакцию 1 июня 2011 г., с доработки — 12 июля 2011 г.
