Содержание к диссертации
Введение
ГЛАВА I. Различные схемы и модификации метода монте-карло, используемые в лучевой терапии 10
1.1 Схемы метода монте-карло, используемые для переноса заряженных частиц
1.2 Существующие реализации метода монте-карло, разработанные для использования в программах планирования дистанционной лучевой терапии 16
1.2.1 VMC 17
1.2.2 VMC++ 18
1.2.3 MMC 18
1.2.4 PEREGRINE 19
1.2.5 SMC 20
1.2.6 MCDOSE 21
1.2.7 DPM 21
1.3 Существующие методы повышения эффективности расчетов при расчете поля фотонов в лучевой терапии 22
1.3.1 Изменение длины свободного пробега 22
1.3.2 Предварительное вычисление источника первых взаимодействий 24
1.3.3 Русская рулетка и расщепление 24
1.3.4 Экспоненциальное преобразование 25
1.3.5 Увеличение числа вторичных частиц 26
1.3.6 Разделение задач, использование заранее рассчитанных результатов 27
1.4 Методы увеличения эффективности расчетов, используемые в лучевой терапии для электронов 28
1.4.1 Поглощение электронов внутри области интереса 28
1.4.2 Исключение электронов по пробегу 28
1.5 Макроскопические теории учета гетерогенностеи при вычислении дозы 29
1.5.1 Теория Брэгга-Грея 29
1.5.2 Теория Спенсера-Аттикса 30
1.5.3 Учет больших полостей в поле фотонного излучения 31
1.5.4 Теория Берлина для фотонного излучения 32
1.6 Выводы 33
ГЛАВА II. Метод pl-оценок потока в точке 35
2.1 Построение дерева траекторий в однородной бесконечной водной среде 38
2.2 PL обработка дерева траектории 39
2.2.1 Построение обратной подобной траектории 40
2.2.2 Вычисление вклада в дозу в детекторной точке от текущего звена траектории 41
2.3 Методы ускорения расчетов при использовании PL-оценок 43
2.3.1 Керма-приближение 43
2.3.2 Сглаживание электронных треков 51
2.4 Учет дивергенции первичного излучения от немононаправленных источников 54
2.4.1 Комптоновское (некогерентное) рассеяние 55
2.4.2 Зависимость углов сдвига от местоположения точки первого взаимодействия 57
2.4.3 Зависимость весового коэффициента от углов рассеяния и сдвига.58
2.4.4 Зависимость весового коэффициента от энергии начальных фотонов 61
2.5 Учет гетерогенностей в методе pl оценок 63
2.5.1 Гетерогенности 1-го рода 63
2.5.2 Гетерогенности П-го рода .- 64
2.6 Выводы 69
ГЛАВА III. Моделирование источников излучения радиотерапевтических установок 71
3.1 Реконструкция эффективного спектра тормозного излучения на основании обработки пространственных дозных распределений 72
3.2 Обработка файлов фазового пространства радиотерапевтических установок 90
3.3 Моделирование источников в методе PL-оценок 111
3.3.1 Моделирование энергетического спектра начальных частиц 112
3.3.2 Моделирование интенсивности 114
3.3.3 Моделирование расходящегося источника 117
3.4 Выводы 119
Заключение 121
Приложение. Назначение и общая структура программного модуля MCPL 123
Ввод исходных данных 124
Ввод и обработка дерева траекторий 124
Описание геометрии 125
Литература
- Существующие реализации метода монте-карло, разработанные для использования в программах планирования дистанционной лучевой терапии
- Предварительное вычисление источника первых взаимодействий
- Построение обратной подобной траектории
- Моделирование источников в методе PL-оценок
Введение к работе
Актуальность работы.
В последние два десятилетия в техническом обеспечении лучевой терапии произошел качественный скачок, связанный с широким распространением медицинских ускорителей, сложнейших систем коллимирования излучения и высокоточной аппаратуры для экспериментальных измерений. Этот скачок естественно сопровождался быстрым развитием новых методов и алгоритмов расчета дозовых распределений, обеспечивающих жесткие требования к точности расчета. Особенно сильное влияние на разработку новых современных методов расчета дозовых распределений оказало развитие и внедрение в клиническую практику техники поперечной модуляции интенсивности пучков, систем обратного планирования и оптимизации облучения. В результате в большинстве клиник стали использоваться системы 3-х мерного дозиметрического планирования с соответствующими алгоритмами 3-х мерного расчета дозы.
В настоящее время в таких передовых системах применяются в основном два метода: 1) метод свертки/суперпозиции, использующий дозовые ядра дифференциального тонкого луча; 2) метод тонкого луча в различных модификациях. Эти методы пришли на смену эмпирическим и полуэмпирическим методам. Они обладают достаточной быстротой и обеспечивают в случае однородности среды высокую точность расчета. Однако при наличии гетерогенности в области расчета, а также вблизи поверхности тела пациента погрешность расчета с помощью этих методов сильно возрастает. Принято считать, что единственным методом, который в условиях сложной геометрии расчета может обеспечить рекомендуемую международными комиссиями точность, является вероятностный метод Монте-Карло. Важно, что метод не требует записи соответствующего уравнения переноса, фактически здесь нет проблем с формулировкой дискретной модели. В обоих случаях реализация вычислительного алгоритма не требует упрощений. Метод максимально при-
5 способлен для использования на ЭВМ. Единственным ограничением являются затраты машинного времени для получения результата с необходимой точностью. В настоящее время эти проблемы не являются решенными. Но в последнее время и в этом направлении ведутся целенаправленные работы и имеются определенные достижения.
Таким образом, метод Монте-Карло является наиболее точным методом расчета доз в лучевой терапии, и поэтому, ускорение расчетов путем модификации метода Монте-Карло является актуальной научно-технической задачей.
Целью работы является усовершенствование одной из модификаций метода Монте-Карло, а именно, метода PL-оценок потока в точке [9], а также разработка новых моделей источников для использования их в методе Монте-Карло.
При этом решаются следующие основные задачи:
1.Поиск возможных способов ускорения расчетов поглощенной дозы с использованием метода PL-оценок. На основе исследований предлагаются методы ускорения расчетов.
Метод PL-оценок потока в точке обобщается на случай расчета доз в неоднородной среде.
Разработка корректных математических моделей источников излучения радиотерапевтических ускорителей и облучателей с б0Со с целью использования их в методе PL-оценок.
Разработка моделей фотонных источников излучения радиотерапевтических установок для использования их в расчетах различными алгоритмами метода Монте-Карло в лучевой терапии.
Теоретическая база исследования. Для решения поставленных в работе задач использовалась теория переноса излучения, теория вероятностей и математическая статистика.
Научная новизна настоящей работы заключается в следующем:
1. Впервые предложены методы ускорения метода PL-оценок потока в
точке.
2.Впервые показана возможность использования метода PL-оценок потока в точке при расчетах доз в химически неоднородных средах.
3.У совершенствованы модели источников используемых в методе PL-оценок потока в точке.
Предложен оригинальный метод формирования функций распределения вероятности, пригодных для выборки при моделировании фотонных источников для расчетов поглощенных доз методом Монте-Карло.
Предложен метод восстановления энергетического спектра дистанционного радиотерапевтического облучателя на основе пространственного распределения энергопоглощений.
Практическая значимость результатов работы состоит прежде всего в том, предложенные в работе быстрые модификации метода PL-оценок с корректным учетом неоднородностей, а также новые математические модели источников дистанционной лучевой терапии, позволили создать на основе метода Монте-Карло методику расчета доз в лучевой терапии за приемлемое для клиник расчетное время.
Достоверность полученных результатов работы определяется использованием корректных теоретических методов, строгостью применяемого математического аппарата, а также хорошим соответствием с экспериментом и результатами расчетов других авторов.
Основные положения, выносимые на защиту:
Способы ускорения расчетов при использовании метода PL-оценок.
Методика учета гетерогенностей при расчете доз в лучевой терапии при использовании метода PL-оценок.
Новые модели фотонных источников излучения для использования их в методе PL-оценок.
Методика реконструкции эффективного спектра тормозного излучения на основании обработки пространственных дозных распределений, измеренных при вводе в эксплуатацию радиотерапевтических ЛУЭ и кобальтовых
7 облучателей.
Апробация работы. Основные материалы и положения работы докладывались и обсуждались на следующих научных конференциях и семинарах:
Научная конференция в Московском инженерно-физическом институте «Научная сессия МИФИ 2004,2005». - Москва, 2004г, 2005 г.
II Евразийский конгресс по медицинской физике и инженерии «Медицинская физика 2005», Москва, 2005 г.
Публикации.
Материалы диссертации опубликованы в 5 печатных работах: 1 статье в научных журналах, 4 научных трудах и тезисах конференций. 1 статья принята к изданию в журнале «Атомная энергия».
Структура и объем диссертации. Диссертация состоит из введения, трех глав, заключения, списка литературы и приложения. Общий объем составляет 131 страницу печатного текста, включая 50 рисунков. Список литературы включает 42 наименования, из них 30 иностранных.
Краткая характеристика работы.
В первой главе проведен анализ литературы по обозначенной тематике. На основе анализа были обобщены и систематизированы существующие на данный момент различные реализации и схемы метода Монте-Карло. Даны разъяснения по каждому типу схем метода Монте-Карло, такие как область применимости и реализованные в них методики ускорения. Описаны существующие программы, реализующие метод Монте-Карло для лучевой терапии. Для каждой программы показаны её преимущества перед стандартной реализацией метода Монте-Карло. Приведен список основных производителей систем дозиметрического планирования лучевой терапии, у которых запланирована (или уже введена) реализация метода Монте-Карло.
Существующие реализации метода монте-карло, разработанные для использования в программах планирования дистанционной лучевой терапии
Большое количество машинного времени требуемое основными реализациями метода Монте-Карло такими как EGS4[23], ITS [24], MCNP [25] или PENELOPE[26], делает их неприменимыми для рутинных расчетов дозовых распределений, проводимых в онкологических центрах. Для адаптации метода Монте-Карло к клиническим расчетам, некоторые М.-К. программы специально адаптируются под условия, которые являются наиболее типичными в лучевой терапии. Например, используются только источники внешнего излучения ограниченного диапазона энергий, а в качестве сред в объеме используются либо определенный набор материалов, либо материалы с маленьким атомным номером. Естественно, что каждое ограничение позволяет рассматривать более узкую область применения алгоритма, что, в свою очередь, позволяет применять различные способы, ведущие к увеличению эффективности расчетов. В последнее время разработаны специальные программные реализации метода М.-К., которые предназначены для использова ния в онкологических учреждениях. В качестве примера можно привести такие программы как ММС [27,28], VMC и XVMC [29-32], VMC++ [37], SMC[33], PEREGRINE [34], MCDOSE[35] и DPM [36].
Несомненно, что ускорение расчетов с использованием метода Монте-Карло вызывает интерес производителей систем планирования лучевой терапии. На данном этапе пока используются системы с упрощениями процесса переноса частиц, но с развитием вычислительной техники таких упрощений будет все меньше и меньше. Уже сейчас некоторые производители включают различные модификации метода М.-К. в свои коммерческие продукты: CMS, Elekta, BrainLAB (XVMC) North American Scientific (PEREGRINE) Nucletron (VMC++) Varian (MMC, VMC++) Эта программа предназначена для расчета дозового распределения от электронных источников. Алгоритм использует 2-х или 3-х мерную геометрию, которая может быть получена с КТ изображений. В отличие от традиционного метода Монте-Карло , этот алгоритм использует различные приближения и упрощения в описании процесса взаимодействия электронов с веществом. Однако использование этих упрощений уменьшает время расчета примерно в 35 раз без существенной потери точности. Для тестирования алгоритма были рассчитаны дозовые распределения в стандартных фантомах. Эти распределения сравнивались с результатами расчетов в аналогичных условиях по программе EGS4 и алгоритму Хогстрома. Во всех случаях наблюдались хорошее соответствие между EGS4 и VMC особенно вблизи границ раздела сред, где алгоритм Хогстрома дает несоответствие дозы до 40%. 1.2.2 VMC++ VMC++ является программным модулем, в котором реализована схема укрупненных соударений при моделировании прохождения фотон-электронного каскада через трехмерную геометрию. Эта методика частично основывается на идеях и алгоритмах разработанных еще в VMC и XVMC методах, однако туда также включены и различные оригинальные разработки моделирования физических процессов и схемы укрупненных столкновений. VMC++ разработан на языке программирования C++ с использованием объектно-ориентированного подхода. В частности, в взаимодействия, относящиеся к моделированию взаимодействия фотонов или электронов, реализованы как различные объекты. Это делает программирование более простым, а также предоставляет более широкие возможности в доработке и модификации программного продукта. Взаимодействие между различными частями программы, такими как моделирование взаимодействий, обработка геометрии и модуль расчета дозы, реализованы через специальные абстрактные методы, что делает код максимально гибким.
Предварительное вычисление источника первых взаимодействий
В случае если фотон покидает область интереса, то к нему можно применить технику называемую русской рулеткой. Эта техника заключается в следующем: разыгрывается случайное число от 0 до 1. Если случайное число выше некоего порога а, то фотон исключается из дальнейшего рассмотрения. Если же разыгранное число ниже порога, то фотон считается «выжившим», а его вес увеличивается в V раз.
Расщепление является процессом обратным русской рулетке. Если фотон находится в области интереса, то генерируется N фотонов с такими же фазовыми координатами (за исключением орбитального угла, который разыгрывается из равномерного распределения) и вес каждого фотона изменя 1 ется в Т7 раз. Метод расщепления подходит, например, для более детального изучения вклада тормозных фотонов в общую дозу.
Экспоненциальное преобразование
Эта техника применяется для повышения эффективности в задачах, в которых важно изучение проникновение излучения в глубокие слои (расчет защиты) или для изучения эффектов на поверхности. Часто этот метод используется для нейтронов, но может использоваться и для фотонов.
Суть метода заключается в следующем. Для простоты будем рассматривать прямоугольный фантом с одинаковой плотностью вещества во всех точках. Одна из плоскостей фантома параллельна оси OZ. Масштабированное расстояние до взаимодействия фотона определяется из формулы A = A(l-C{i), где X расстояние до взаимодействия в единицах длины свободного пробега, X масштабированное расстояние, ju косинус угла между направлением движения фотона и осью OZ, С - параметр масштабирования. Вероятность взаимодействия в этом случае будет где множитель (1 - Cjil) введен для того, чтобы обеспечить нормировку СО плотности вероятности на 1, т. е. )Р [Лрл - 1. о Для Cfi=0 имеет место несмещенная плотность вероятности. Для 1, среднее расстояние до взаимодействия увеличивается, и для С/и 0 расстояние уменьшается. Для того, чтобы обеспечить несмещенность оценок необходимо ввести весовой коэффициент исходя из следующей пропорции Решая это уравнение относительно со , получим Расстояние до следующего взаимодействия с учетом взаимодействия разыгрывается из формулы
В некоторых приложениях желательно изучить поведение вторичных частиц в режимах, где они появляются с малой вероятностью. Например, рентгеновское излучение в диагностических аппаратах есть не что иное, как тормозное излучения. Однако, в диагностических режимах энергия электронов около 70 кэВ. При такой энергии вероятность образования тормозного излучения во много раз меньше вероятности образования дельта-электронов. Следовательно, изучение процессов образования тормозного излучения займет большое количество времени. Другим примером может служить процесс образования пар для материалов с низким Z или при энергии первичных фотонов близкой к пороговой. При таких условиях основным процессом является комптоновское рассеяние. Для того, чтобы изучить влияние позитронов на общий вклад в дозу потребуется очень большое количество машинного времени.
Чтобы уменьшить время расчета и увеличить вклад изучаемых частиц предложено использовать технику, называемую «увеличение числа вторичных частиц». Эта методика аналогична методу расщепления. В точке образо 27 вания вторичной частицы вместо одной рождается несколько частиц. Импульсные и энергетические характеристики рассчитываются из соответствующих распределений образования вторичных частиц. Для того чтобы обеспечить «честность» розыгрыша, вес каждой вторичной частицы уменьшается в ТУраз, где N- число рожденных частиц. При этом подходе нарушается закон сохранения энергии для траектории, т.к. рождается целый веер вторичных частиц, в то время как поглощенная энергия считается как разница между первичной и одной из вторичных частиц. Однако при большом числе историй и множестве взаимодействий, закон сохранения энергии в целом выполняется.
Построение обратной подобной траектории
Продолжительность расчета в методе PL-оценок прямо пропорционально числа ветвей в траектории. В текущей реализации метода есть два способа расчета. В первом берутся во внимании только электронные звенья траектории, во втором только фотонные. Очевидно, что при использовании электронных звеньев уменьшается погрешность вычислений, но из-за их большого количества это ведет к увеличению продолжительности расчетов. Второй способ работает значительно быстрее, но точность результатов ниже. Кроме того, этот способ можно применять только в тех областях, где есть электронное равновесие.
Как видно из табл. 2.1, скорость расчета с использованием только фотонных звеньев выше в 1,5 - 2,5 раза в зависимости от энергии. Уменьшение относительного числа фотонных звеньев с ростом энергии связано с увеличением энергии вторичных электронов, что в свою очередь влияет на пробег электронов в среде, и соответственно на количество электронных треков.
Для того чтобы наиболее эффективно применять метод PL-оценок, в на стоящей работе предложено использовать способ расчета, который является комбинацией электронного и фотонного способов.
Введем понятие уровня взаимодействия для вторичной частицы. Уровень взаимодействия частицы - это количество фотонов, образовавшихся на траектории до рождения текущей частицы. На примере это выглядит так. Допустим, что начальный фотон испытал первое взаимодействие со средой, и это взаимодействие является некогерентным рассеянием. Для понятия уровней взаимодействия процесс рассматривается как уничтожение родительской частицы и образования нескольких вторичных. У вторичных частиц будет первый уровень взаимодействия, т.к. первичный фотон испытал одно столкновение. Далее, допустим, что вторичный фотон образовал пару электрон-позитрон. Тогда у родившихся частиц будет второй уровень взаимодействия и т. д. Заметим, что при взаимодействии электрона с веществом не происходит изменения уровня продуктов взаимодействия.
Электрон достаточно быстро теряет свою энергию с ростом количества взаимодействий. Зависимость вклада в дозу от электронов определенного уровня взаимодействия показана на рис. 2.1
Как видно из рис. 2.1, вклад в дозу от электронов падает примерно в 10 раз за 3 - 4 взаимодействия. Поэтому при низком уровне взаимодействия можно рассчитывать дозу от электронов, а при высоком уровне взаимодействия - от фотонов.
Для дальнейших рассуждений понадобится понятие «граничный уровень взаимодействия». Граничный уровень взаимодействия - уровень взаимодействия вторичной частицы, ниже которого ведется расчет по электронной составляющей траектории, а выше - по фотонной составляющей. Вводя в расчет граничный уровень взаимодействия, можно изменять его, добиваясь тем самым оптимального соотношения скорость-точность для расчета дозы.
Предположим, что в выборке траекторий максимальный уровень взаимодействия есть N Для звеньев с уровнем взаимодействия I К вклад в дозу рассчитывается по электронной составляющей траектории, а для остальных - по фотонной. Ниже приведена зависимость коэффициента а от выбора граничного уровня взаимодействия К. Рис. 2.2, Рис. 2.3.
Из рис. 2.2 и 2.3 видно, что предложенная методика расчета дает максимальный эффект при граничном уровне взаимодействия равным 1. С увеличением граничного уровня время расчета резко возрастает и становится примерно равным времени расчета по стандартной электронной схеме.
Такая зависимость обусловлена изменением количества электронов и фотонов различных уровней взаимодействия. Для того чтобы оценить количество электронов и фотонов в выборке траекторий, были проведены расчеты. Результаты показаны на рис. 2.4 и 2.5. Количество частиц каждого уровня взаимодействия отнормировано на общее число частиц. Рассматривались моноэнергетические источники с энергией 1,25 и 5 МэВ в полубесконечной гомогенной водной среде, рис 2.4, рис.2.5.
Если сравнивать рис. 2.4 и 2.5, то можно заметить, что количество электронов для источника с энергией 5 МэВ выше при малых уровнях взаимодействия, чем для источника с энергией 1,25 МэВ. Это связано с тем, что при энергии 5 МэВ эффект образования пар играет более весомую роль. Электроны, образованные взаимодействием 5 МэВ фотона имеют более высокую энергию, а значит и более длинный пробег, что также влияет на увеличение их количества. Более резкий спад в количестве фотонов и электронов при более высоких уровнях взаимодействия для энергии 1,25 МэВ, очевидно, связан с меньшей энергией вторичных частиц.
Моделирование источников в методе PL-оценок
Процессом некогерентного (комптоновского) рассеяния называется процесс рассеяния гамма-кванта на электронных подоболочках атома, сопровождающийся выбиванием одного из атомных электронов. Дифференциальное сечение некогерентного рассеяния гамма-кванта часто принято заменять на дифференциальное сечение комптоновского рассеяния на свободном покоящемся электроне (формула Клейна-Нишины-Тамма [7]):
Поскольку реально гамма-квант рассеивается на связанных электронах оболочек атома, то для учета связанности электрона в атоме сечение dcsKNT/dQ умножается на функцию некогерентного рассеяния S(q,Z), где q - передаваемый импульс, Z - порядковый номер химического элемента. Но в нашем случае, это неважно, т.к. используется отношение сечений комптоновского рассеяния для различных углов в.
Первоначальный набор траекторий, используемый в методе PL-оценок для расчета дозы, насчитывается в гомогенной водной среде от мононаправ ленного источника. При расчете дозы в детекторной точке, каждое звено траектории протягивается до точки первого взаимодействия. Очевидно, что состояние любой частицы в траектории зависит от состояния «родительских» частиц. Допустим, что после первого взаимодействия в наборе траекторий фотон рассеялся на угол в. Для расходящихся источников, угол между начальным фотоном и фотоном, после первого взаимодействия будет составлять в . Это, в свою очередь, ведет к изменению энергии рассеянного фотона, и изменению состояния всех последующих частиц.
Учет расхождения лучей. Сплошная линия -траетория из заранее насчитанного набора, пунктирная линия - прямая между точкой первого взаимодействия и точечным источником, в - угол рассеяния фотона для мононаправленного источника, в -угол рассеяния фотона при расходящемся источнике
Для того чтобы компенсировать изменения в различиях между углами рассеяния, предлагается ввести весовой коэффициент, который в общем виде выражается как комптоновское дифференциальное сечение рассеяния на угол -комптоновское сечение рассеяния на угол в, т.е. угол рассея ния, взятый из исходной траектории. Далее, вклад в дозу от текущего звена траектории рассчитывается с учетом весового коэффициента со. D co-D , (2.28) где D- - вклад в дозу от текущего звена траектории, без учета расхождения лучей.
Допустим, что после совмещения текущего звена траектории с детекторной точкой, точка первого взаимодействия попала в область первых столкновений и имеет координаты (xo,yo,zo)- Размер поля равен 2Dx2D Источник имеет координаты (0,0, -R) . Тогда косинус угла сдвига вычисляется как
Условие попадания в область первых столкновений cos(A#) cos(0), где cos(0) максимальный угол раствора поля, определяемый из соотноше ния Ниже приведены максимальные углы раствора, для различных полей и SSD. Как видно из таблицы углы раствора не превышают 10 градусов. Однако поля 40x40 и 20x20 см, являются очень большими и достаточно редко применяются. Поэтому максимальные углы раствора обычно имеют значения до 3 градусов.
Точка первого взаимодействия находится внутри области первых взаимодействий. Следовательно, угол сдвига будет не больше угла максимального раствора поля и будет линейно зависеть от поперечного расстояния от точки первого взаимодействия до оси источника.
Для определенности будем считать, что энергия начального фотона равна EQ. Как известно, вероятность рассеяния на определенный угол при ком-птоновском эффекте задается формулой (2.26) , рис.2.11 демонстрирует эту зависимость для энергий 1.25 МэВ и 5 МэВ.
Очевидно, что изменение угла падения начального фотона, можно интерпретировать как сдвиг функции распределения вдоль оси X на величину Ав. Тогда весовой коэффициент вычисляется, согласно определению (2.27).
