Климанов Радиобиологическое и дозиметрическое планирование 2011
.pdfrp |
|
|
|
E0 Eth f |
|
z |
] |
|
|
|
0 |
[0,042 |
|
|
|
|
|
||||
|
|
|
|
|||||||
|
|
|
|
Mc2 |
|
|
RCSDA |
|
|
(8.108) |
|
[1 erf((RCSDA z zshift (E0 )) / total )] |
1 |
, |
|||||||
|
2 |
|||||||||
|
|
|
|
|
|
|
|
|
|
где total 2mono inelastic2 .
Если протоны являются полиэнергетическими, то во всех трех формулах (8.106) –(8.108) ζmono следует заменить на ζp. Значения
zshift для разных начальных энергий протонов приводятся в табл.
8.8.
Рис. 8.30. Уменьшение флюенса первичных протонов разных энергий с учетом ядерных взаимодействий и флуктуации потерь энергии [44]
В работе [44] Улмер получил также формулы для тормозной способности протонов с учетом ядерных взаимодействий и флуктуаций. Однако эти выражения довольно громоздки, Учитывая, что в литературе имеются подробные данные для тормозной способности протонов (например, [26]), эти выкладки здесь не рассматриваются. Приведем в заключение новые формулы, предложенные Улмером для параметра p в формуле (8.93):
• нерелятивистский случай
p 5 10 15 E06 5 10 12 E05 2 10 |
9 E04 |
(8.109) |
|||
4 10 7 |
E03 5 10 5 |
E02 0,0003 E0 1,6577; |
|||
|
81
• релятивистский случай
p 4 10 15 E06 4 10 12 E05 2 10 9 E04 |
(8.110) |
|||||||
4 10 7 E03 5 10 5 |
E02 |
0,0027 E0 1,6576. |
||||||
|
|
|||||||
|
|
|
|
|
|
Таблица 8.8 |
||
Значение параметра zshift в формулах (8.107) и (8.108) для протонов |
||||||||
разных энергий [44] |
|
|
|
|||||
|
|
|
|
|
|
|
|
|
Энергия (расчетная), МэВ |
97,66 |
|
160,15 |
178,09 |
141,31 |
230,54 |
|
|
Энергия (номинальная), МэВ |
97,00 |
|
160,00 |
177,00 |
141,11 |
229,92 |
|
|
zshift , см |
0,17 |
|
0,23 |
0,43 |
0,31 |
0,26 |
|
11. Расчет дозового распределения на основе измерения флюенса протонов в воздухе
Большинство алгоритмов расчета дозы для пучков протонов и систем дозиметрического планирования создавались для конкретных ускорителей и способов облучения. В них широко используются результаты измерения дозовых распределений в водном фантоме для различных комбинаций входных параметров конкретной машины. Более универсальным является подход, в котором методика расчета дозы не зависит от особенностей машины, а специфика конструкции ускорителя учитывается через конфигурационные измерения основных характеристик линии пучка и их обработку. Такой подход был развит Б. Шаффнер в работе [45]. Остановимся на ней подробнее в части, относящейся к машинам с фольговой системой расширения пучка.
11.1. Общая характеристика метода расчета
Алгоритм расчета протонной дозы, базирующийся на определе-
нии распределения флюенса в воздухе, включает три основных ша-
га [45]:
1. Расчет полного водоэквивалентного расстояния для каждой точки расчетной сетки.
82
2.Расчет распределения флюенса протонов в воздухе для каждой энергетической интервала в спектре пучка (энергетического слоя) при заданной конфигурации линии пучка.
3.Свертка/суперпозиция флюенса для каждого энергетического слоя с дозовым распределением узкого мононаправленного пучка
(в англоязычной литературе часто используется краткий термин бимлет (англ. beamlet)), масштабированным вдоль оси пучка, ис-
пользуя полное водоэквивалентное расстояние для каждой расчетной точки.
Первые два шага являются независимыми друг от друга. На третьем шаге соединяются данные пациента и поперечное распределение протонов, производимое линией пучка, с физической моделью бимлета. Физическую модель бимлета протонов в воде Шаффнер взяла из работы Улмера [44], рассмотренной в предыдущем разделе.
11.2. Определение конфигурационных параметров из измерения флюнса в воздухе
Устройства и их расположение на линии пучка в настоящее время сильно отличаются между собой у разных производителей и для различных способов облучения. Поэтому для систем планирования облучения является обязательным такое построение, при котором они могли бы поддерживать различные машины, моделируя линию пучка с помощью небольшого количества характерных параметров с возможность их определения на базе конфигурационных измерений. Ключевыми величинами ко всем характеристическим параметрам в работе [45] служат номинальная энергия на входе носика и толщина водного эквивалента от входа в носик (NeT). Результаты аппроксимации представлялись графически как функции от этих двух переменных и как функции одного комбинированного параметра, которым являлся остаточный пробег:
Rres (E, NeT ) RCSDA NeT , |
(8.111) |
где RCSDA – пробег протонов с энергией E в приближении непре-
рывного замедления, который рассчитывался по методу Улмера
[44].
В число параметров, используемых для расчета флюенса в воздухе, входят виртуальное и эффективное SAD и размер эффективного
83
источника для однородного сканирования или применения на линии пучка техники рассеивающих фольг.
11.2.1. Определение виртуального SAD
Под виртуальным SAD, согласно [46], понимается расстояние от изоцентра до положения виртуального источника. SADvirt определяется в работе [45] с помощью линейной аппроксимации зависимости 50 % 50 % ширины поля от расстояния до изоцентра. Пози-
ция виртуального источника находится в точке, где линейная аппроксимация пересекает ось абсцисс (рис. 8.31,а). Процесс включает два шага:
1. Определение 50 % 50 % ширины поля на разных расстояни-
ях с помощью аппроксимации функцией: |
|
|
|
|
|
|
|
|
|
||||||||||
|
N |
FS |
|
x |
|
x FS |
|
|
|
||||||||||
f (x) |
|
erf |
|
|
|
|
|
erf |
|
|
|
|
|
, |
(8.112) |
||||
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
2 |
|
2 |
|
|
|
|
2 |
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
где N – нормировочный множитель; FS+ – полуширина поля в положительном направлении оси x; FS_– полуширина поля в отрицательном направлении оси x; ζ+ – параметр спада в положительном направлении оси x (мм); ζ_– параметр спада в отрицательном направлении оси x (мм). Параметры ζ+ и ζ_ находятся через процедуру подгонки.
2. Линейная аппроксимация размера поля ( FS FS FS ) в
зависимости от расстояния до изоцентра (z) для нахождения пересечения с осью абсцисс.
Положение SADvirt в системах с рассеивающими фольгами, в общем случае, зависит как от номинальной энергии, так и от эквивалентной толщины линии пука в носике. Однако последние конструкции коммерческих ускорителей отличаются тем, что с помощью комбинации материалов на линии пучка в носике достигается минимальное изменение виртуального SADvirt. Имеется еще зависимость SADvirt от энергетического слоя, но для убыстрения расчетов в работе [45] используется постоянное значение, которое берется для самого проникающего энергетического интервала в спектре. Погрешность в расчете дозы от этого приближения считается пренебрежимо малой.
84
Рис. 8.31. Схематическое изображение измерений, требуемых для конфигурации систем рассеивающих фольг или однородного сканирования: а – измерение поперечных профилей широкого поля на разных расстояниях для определении SADvirt; б – измерение флюенса в воздухе для определения SADeff; в – измерения
пенумбры для определения эффективного размера источника (адаптировано из
[45])
85
11.2.2. Определение эффективного SAD
Эффективное SAD определяется, согласно [46], как положение источника, при котором выполняется закон обратных квадратов для ослабления пучка в воздухе (рис. 8.31,б). SADeff в общем случае зависит от энергии протонов на входе в носик и количества поглощающего материала на линии пучка. Для его определения проводится измерение флюенса в зависимости от положения детектора на оси пучка. Затем в работе [45] выполнятся подгонка к результатам измерений функции
SAD 2
f (z) N eff , (8.113)
(SADeff z)2
где N – нормировочный множитель; SADeff – расстояние между изоцентром и эффективным источником, являющееся подгоночным параметром.
Измерения рекомендуется проводить на центральной оси пучка в пределах ± 20 см от изоцентра, т.е. в области, где располагается пациент.
11.2.3. Определение эффективного размера источника
Концепция эффективного размера источника и его использование для планирования протонной терапии впервые были предложены в работе [36]. В работе [45] эффективный размер находится в номинальном положении источника. Входными данными для определения системой планирования размеров эффективного источника являются измерения пенумбры в воздухе при частичном перекрытии пучка половинным блоком, располагаемым на разных расстояниях вдоль оси (рис. 8.31,в). Процедура состоит из следующих шагов:
1. Определение параметра пенумбры ζp с помощью подгонки к экспериментальным данным функции:
86
|
|
|
|
|
|
x0 |
|
|
||
N |
|
|
|
x |
|
|
||||
|
|
|
erf |
|
|
|
|
|
1 , |
|
2 |
|
|
p |
|||||||
2 |
||||||||||
|
|
|
|
|
|
|
|
|
|
|
f (x) |
|
|
|
|
|
|
|
|
|
|
N |
|
|
|
|
||||||
|
|
|
x0 x |
|
||||||
|
|
|
|
erf |
|
|
|
|
|
1 , |
|
|
|
|
|
|
|
||||
|
2 |
|
|
|
|
2 p |
|
|
||
|
|
|
|
|
|
|
|
|
|
блокируются отрицательные x;
блокируются положитель ные x,
(8.114)
где x0 – положение проекции блока по отношению к измеряемой величине (т.е. 50 % флюенса); ζp – параметр спада пенумбры (мм).
2. Линейная аппроксимация зависимости параметра ζp от расстояния до изоцентра (z). Эффективный размер источника находится как значение аппроксимационной прямой в точке расположения номинального источника.
11.3. Использование конфигурационных параметров для расчета дозы
Параметры, полученные в результате обработки конфигурационных измерений в воздухе, табулируются как функции ключевых переменных для каждого энергетического слоя. Из них определяются распределения флюенса и дозы.
11.3.1. Расчет распределения флюенса
Расчет флюенса в воздухе при отсутствии компенсатора выполняется в [45] при следующих упрощающих допущениях:
распределение флюенса внутри поля является однородным;
блок является непроницаемым, бесконечно тонким барьером для протонов, т.е. рассеяние от блока отсутствует и, следовательно, нет геометрических эффектов, вызываемых конечной толщиной апертуры.
Распределение флюенса в воздухе в этом случае описывается следующим выражением:
87
0 |
|
SADeff ,l z 2 |
|
|
x (z) |
|
|
x (z) |
|
|||||||||||||||
l |
N |
|
|
|
|
|
|
erf |
|
|
|
|
|
|
|
erf |
|
|
|
|
||||
SAD |
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
|
|
|
|
|
|
|
|
|
|
|||||||||||||||
|
|
|
|
|
|
2 p,l |
(z) |
|
|
2 p,l (z) |
|
|
||||||||||||
|
|
|
|
eff ,l |
|
|
|
|
|
|||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
y |
(z) |
|
|
|
|
y |
(z) |
|
|
|
|
|
|
|
|||||||
|
erf |
|
|
|
|
|
|
|
erf |
|
|
|
|
|
|
|
, |
(8.115) |
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||
|
|
|
|
2 p,l (z) |
|
|
|
|
2 p,l (z) |
|
|
|
|
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
где l0 – флюенс для энергетического слоя l в точке (x,y,z) в отсут-
ствие компенсатора; N – нормировочный фактор, зависящий от веса слоя l и от условий нормировки; ось z имеет начало в изоцентре и положительное направление вперед к источнику пучка;
x , y – расстояния от точки расчета до проекции блока на ось z,
при этом используется виртуальное SAD; ζp,l – параметр пенумбры для слоя l в положении z, которое определяется из эффективного размера источника для слоя l и корректируется на позицию блока, если она отличается от калибровочной. Расчет ζp,l проводится по формуле:
p,l |
|
IBDcurrent z |
|
SADnom IBDcurrent |
p,l , (8.116) |
SADnom IBDcurrent |
|
||||
|
|
|
SADnom IBDref |
где IBDcurrent, IBDref – текущее и референсное расстояния изоцентрблок. Первый член в (8.116) преобразует эффективный размер ис-
точника в параметр пенумбры, второй член конвертирует эффективный размер источника, полученный из конфигурационных измерений при референсных условиях в эффективный размер источника для текущей позиции блока.
Расчет флюенса в открытой от блока части выполняется заменой произведения функций ошибок в уравнении (8.115) на сумму функций ошибок. Поправки на неоднородность флюенса в открытом поле и конечную толщину апертуры могут быть, в принципе, включены в (8.116) как добавочные факторы, используя рекомендации работы [47].
Большинство облучений в протонной терапии выполняется с индивидуальными для каждого пациента компенсаторами. В работе [45] рассеяние в компенсаторе моделируется как возмущение флюенса для каждого слоя. Расчет ведется по формуле:
88
|
|
0 |
|
|
|
0,5xz |
|
xi |
|
|
|
|
0,5xz |
|
xi |
|
|
||||||||
c |
|
l |
|
|
|
|
|
|
|
||||||||||||||||
l |
|
|
erf |
|
|
|
|
|
|
|
erf |
|
|
|
|
|
|
|
|||||||
4 |
|
|
|
|
c |
|
|
|
c |
||||||||||||||||
|
|
|
|
|
|
|
|||||||||||||||||||
|
|
|
i, j |
|
|
2 i, j |
|
|
|
|
|
2 i, j |
|
|
|||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(8.117) |
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
0, yz |
|
yi |
0,5 yz |
|
yi |
|
|
|||||||||||||||
|
|
|
|
|
|
|
|
||||||||||||||||||
|
|
erf |
|
|
|
|
|
|
|
erf |
|
|
|
|
|
|
|
, |
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
2 ic, j |
|
|
|
2 ic, j |
|
|
|
|
|||||||||||
|
|
|
|
|
|
|
|
||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
где функции ошибок являются интегралами от гауссовского поперечного распределения протонов, вызванного рассеянием в i,j- пикселе компенсатора (суммирование ведется по всем пикселям компенсатора); l0 – флюенс для слоя l в точке (x,y,z) без компен-
сатора (уравнение (8.116)); xz и yz – размеры пикселей в дивергентной расчетной сетке в позиции z (дивергентная (веерная) сетка характеризуется SADvirt и размером пикселя в изоцентре);
xi , yi – расстояние от точки расчета до проекции ближних пикселей компенсатора, который представляется в той же дивергентной сетке; i, j (z) – стандартное отклонение гауссовского распре-
деления в позиции z протонного пучка, рассеянного в i,j-пикселе компенсатора. Величина стандартного отклонения в работе [45] рассчитывалась по методике работы [37].
11.3.2. Расчет распределения дозы
Методика расчета дозового распределения из распределения флюенса протонов в воздухе, примененная Шаффнер, основывается на использовании аналитической модели Улмера (см. предыдущий раздел). Для каждого энергетического слоя рассчитывается трехмерное дозовое распределение, создаваемое в воде соответствующим бимлетом (бесконечно тонким мононаправленным моноэнергетическим источником или тонким лучом). При необходимости определения дозы в другом материале проводится масштабирование этого распределения в соответствии с интегральным водяным эквивалентом вдоль центральной оси и умножение на двухмерную ступенчатую функцию, равную единице внутри пикселя расчетной сетки и нулю вне пикселя. Полная доза для каждого энергетического слоя равна сумме вкладов от всех бимлетов, умноженных на величину флюенса в позиции бимлета. Таким обра-
89
зом, дозовое распределение рассчитывается как свертка распределения флюеса в воздухе и дозового распределения бимлетов.
Результаты расчета распределений поглощенной дозы, выполненные в работе [45], вполне удовлетворительно согласуются с экспериментальными данными для разных машин.
12. Применение метода Монте-Карло для расчета доз от протонов
Метод тонкого луча, который используется в большинстве коммерческих систем дозиметрического планирования, при всех его достоинствах (высокая скорость расчета, приспособленность к лучевой терапии с модуляцией интенсивности пучка и др.) имеет ограниченную точность в задачах с существенными негомогенностями. Этот недостаток связан с одномерным масштабированием данных по ТЛ протонов в воде на другие среды. Для преодоления
этого недостатка разрабатываются более сложные методы масштабирования [40,48,49] (см. также раздел 9 настоящей главы), однако
вблизи границ раздела разных сред они не обеспечивают требуемую точность. Дальнейший прогресс возможен только на основе использования метода Монте-Карло.
Другая серьезная проблема в протонной лучевой терапии – погрешности, связанные с укладкой пациента и движением органов. Эти неопределенности имеют серьезное влияние на процесс облучения, потому что протонные пучки можно направить с высокой точностью к специфической области внутри пациента. Вместе с тем, небольшое смещение пика Брэгга может привести к ―недодозированию‖ опухоли и, наоборот, к ―передозированию‖ критических органов. В настоящее время эта проблема преодолевается, главным образом, путем включения в CTV дополнительного объема. Но в этом случае частично теряется высокая точность протонных дозовых распределений. Крайне желательно поэтому уменьшать искусственное увеличение CTV, учитывая движение органов при планировании облучения. Однако это затруднительно сделать в алгоритме ТЛ, потому что требуется проводить лучевой анализ изменяющегося во времени распределения плотности. В то же время для алгоритма метода Монте-Карло такой проблемы не существует.
90