Содержание к диссертации
Введение
Глава 1. Состояние вопроса (обзор) 16
1.1. Математические м одели радиацион ных пол ей 16
1.2. Численные методы расчета радиационных полей 26
1.3. Математическая модель тепловых полей 32
1.4. Численные методы расчета тепловых полей 36
1.5. Программные комплексы для расчета радиационных и тепловых полей 41
Глава 2. Метод характеристических функций в оценивании математического ожидания случайных величин с бесконечной дисперсией 45
2.1. Асимптотики характеристической функции 47
2.2. Оценки математического ожидания 53
2.3. Смещение 56
2.4. Дисперсия 59
2.5. Погрешности 62
2.6. Выбросы 65
2.7. Численное исследование оценок 67
Глава 3. Моделирование радиационных и тепловых полей 79
3.1. Математическая модель радиационных полей 80
3.2. Численный метод расчета радиационных полей 84
3.3. Математическая модель переноса тепла 92
3.4. Решение биотеплового уравнения 95
Глава 4. Программный комплекс для моделирования радиационных и тепловых полей. Тестирование и приложения 111
4.1. Общие характеристики комплекса 111
4.2. Тестирование программного комплекса 113
4.3. Численное исследование локальных оценок 121
4.4. Практическое моделирование лазерных процедур 126
Заключение 135
Список используемой литературы 137
- Численные методы расчета радиационных полей
- Программные комплексы для расчета радиационных и тепловых полей
- Численный метод расчета радиационных полей
- Тестирование программного комплекса
Введение к работе
Настоящая работа посвящена моделированию процессов переноса излучения и тепла в биологических тканях, подвергающихся интенсивному лазерному облучению.
Такие исследования актуальны для лазерной терапии и хирургии, биооптики, теплофизики. В медицине лазерные источники излучения и оптоволоконная техника позволили создать целый ряд новых эффективных медицинских технологий, существенно уменьшающих сроки лечения, его стоимость, риск и тяжесть осложнений. К таким технологиям, в частности, относятся термотерапия (лазерная гипертермия), каналирование и перфорация мягких и костных тканей, фото динамическая терапия, технологии, основанные на термокоагуляции и термоабляции тканей. Все они нашли широкое применение в онкологии, отоларингологии, дерматологии, урологии, нейрохирургии, кардиохирургии, косметологии и других отраслях медицины.
Главными физическими агентами воздействия во всех перечисленных процедурах, вызывающими конечный медицинский эффект, являются излучение и тепло. Следовательно, для исследования воздействия лазерного излучения на биоткани необходимо, прежде всего, иметь информацию о радиационных и тепловых полях, возбуждаемых в процессе облучения.
Одним из способов получения такой информации являются приборные методы. Большинство используемых методов являются проникающими: регистрирующие датчики помещаются внутри ткани. Преимуществом таких методов является их точность и надежность. Но они имеют существенные недостатки:
введение датчиков создает дополнительную травму ткани, которая
часто недопустима по медицинским показаниям, т.к. ведет к
повреждению жизненно важных объектов (сосуды, нервные узлы);
получаемая информация, относится лишь к малому числу точек.
Существующие непроникающие приборные методы
(дистанционные и контактно-поверхностные) позволяют получать лишь, поверхностное распределение нужных характеристик. Ведутся разработки и для непроникающих измерений радиационных и температурных полей внутри тканей, но пока они далеки от внедрения в клиническую практику.
Альтернативным методом получения информации о радиационных и тепловых полях является метод математического компьютерного моделирования. Он, в принципе, позволяет получать пространственные распределения характеристик; проводить исследования «без повреждения» исследуемого объекта и существенных затрат, которых требует эксперимент.
Разработка средств компьютерного моделирования включает 3 задачи:
создание математической модели (исходные уравнения, граничные условия, параметры);
разработка алгоритмов реализации модели (аналитические, численные, статистические методы);
создание компьютерной программы, реализующей эти алгоритмы.
К настоящему времени в литературе предложен ряд моделей, численных алгоритмов и компьютерных программ для моделирования радиационных и тепловых полей. Проведенный анализ литературы (Глава 1) позволяет сделать следующие выводы.
1. Наиболее адекватной математической моделью переноса лазерного излучения в мутных средах является кинетическая модель, т.е.
кинетическое уравнение переноса излучения в сочетании с законами Френеля, описывающими граничные условия. Именно эта модель принята в настоящей работе.
Наиболее адекватные математические модели переноса тепла в биотканях основаны на нестационарном уравнении теплопроводности и отличаются разными способами учета конвективного переноса тепла за счет кровотока.
Учет теплового эффекта кровотока в больших сосудах (вены, артерии) осуществляется явным образом. Из-за сложной геометрии кровеносных систем такие модели не являются универсальными, они разрабатываются отдельно для конкретных биотканей.
Учет капиллярного кровотока осуществляется двумя методами: «кен »-метод и «bio heat transfer»-MeTOfl.
а) Идея <е#»-метода заключается в учете теплового эффекта
капиллярного кровотока в коэффициенте теплопроводности. Такой
способ не требует изменения вида исходного уравнения (уравнения
теплопроводности) и введения дополнительных параметров, но на наш
взгляд, является не вполне обоснованным физически.
б) «Bio heat transfer»-MeTOfl - физически более обоснован. Учет
капиллярного кровотока осуществляется здесь введением в уравнение
теплопроводности дополнительного конвективного слагаемого,
пропорционального разности температуры среды и температуры крови.
Такой метод содержит дополнительные параметры и требует решения
уже более сложного уравнения, называемого в мировой литературе «bio
heat transfer equation». Поскольку в отечественной литературе перевод
этого термина пока не закрепился, мы примем свой: «биотепловое»
уравнение. В стационарном случае оно совпадает с уравнением
Гельм гольца.
Сопоставления, проведенные в ряде работ, позволяют заключить, что, несмотря на простоту, первая модель может применяться в случае
достаточно мелких капилляров. Вторая модель позволяет учитывать капилляры практически любых размеров.
В настоящей работе в качестве тепловой модели выбрана комбинация этих моделей, включающая в качестве частных случаев keff и «bio heat transfer» модели, а также уравнение теплопроводности без конвекции.
2. Наиболее точным методом реализации кинетической модели, позволяющей учитывать сложную геометрию среды и источников излучения, является метод Монте-Карло. Во всех работах по биооптике реализован, так называемый, нелокальный метод Монте-Карло, что является существенным недостатком при расчете радиационных полей. Такие методы дают значения радиационных характеристик, усредненные по некоторым областям фазового пространства, что может приводить к большим систематическим погрешностям в случаях больших градиентов радиационных полей, весьма характерных для практических лазерных процедур. В то же время в теории методов Монте-Карло для расчета полей ионизирующего излучения разработаны локальные алгоритмы, позволяющие рассчитывать нужные характеристики радиационного поля в заданных точках пространства. Наиболее известной и часто применяемой на практике является локальная оценка Калоса (см. Глава 1 «Математические модели радиационных полей»), позволяющая рассчитывать плотность потока частиц (гамма-квантов, нейтронов) в фиксированной точке. Платой за это преимущество является бесконечность дисперсии этих оценок и, как следствие, пониженная скорость сходимости, плохое качество сходимости (велика вероятность больших отклонений от оцениваемой величины), что существенно затрудняет практическое использование этих оценок. Несмотря на то, что для преодоления этого недостатка предложен ряд, так называемых, ускоренных локальных методов, проблема далека от окончательного решения.
Литературный анализ показал, что численные методы для принятой тепловой модели не достаточно развиты. Во всех работах используются громоздкие конечно-разностные методы. В то же время эта модель может быть реализована любым численным методом, разработанным для уравнения теплопроводности, решению которого посвящена обширная литература.
Для решения стационарного уравнения теплопроводности наибольшее применение получили: метод конечных элементов; конечно-разностный метод; метод взвешенных невязок (решение параметризуется некоторой функцией; ее параметры находятся из условия ортогональности невязки к заданным весовым функциям). Для различных приложений, включая биомедицинские, имеющих дело со сложной геометрией и гетерогенными средами, хорошо зарекомендовал себя метод конечных элементов.
Для решения нестационарного уравнения теплопроводности, в основном, применяются методы, основанные на его сведении к системе более простых уравнений. Такое сведение осуществляется разными способами. Наиболее распространенным является сведение к системе обыкновенных дифференциальных уравнений (по времени) применением какого-нибудь численного алгоритма к пространственной части уравнения, т.е. к стационарному уравнению теплопроводности, в предположении известной производной по времени. Для решения стационарного уравнения используются те же методы: метод конечных элементов, конечно-разностный метод. Для решения системы обыкновенных дифференциальных уравнений используются в основном методы численного интегрирования Рунге-Кутта, конечно-разностный метод, метод взвешенных невязок.
Недостатком этого способа является громоздкость: приходится численно решать две системы. Реализация этого способа особенно
трудна для задач, в которых параметры уравнения зависят от температуры.
Другим способом решения нестационарного уравнения является сведение к системе дифференциальных уравнений применением какого-нибудь конечно-разностного алгоритма к временной части. Наиболее интересными, в практическом смысле, являются схемы, приводящие к пошаговому решению уравнения по времени; температурное поле в следующий момент времени находится из температурного поля в предыдущий момент времени. В таких схемах на каждом временном шаге приходится решать одно стационарное уравнение. Часто — это уравнение типа Гельмгольца. Этот подход проще и естественнее предыдущего, не возникает проблем, связанных с зависимостью параметров от температуры.
Для численной реализации биотеплового уравнения, на наш взгляд, последний способ является более естественным: биотепловое уравнение и в стационарном случае является уравнением типа Гельмгольца.
3. В настоящее время в мире предприняты попытки создания компьютерных программ для моделирования радиационных и тепловых полей. Наиболее полными являются программные комплексы: LATIS (London R.A. et al., Lawerence Livermore National Laboratory, California) и LITCIT (Roggan Л. et al., Laser-Medizin-Zentrum gGmbH, Berlin).
Но они являются чисто лабораторным продуктом и не доступны для сторонних пользователей. Достоинством этих комплексов является использование наиболее адекватных математических моделей: кинетическая модель переноса излучения (т.е. кинетическое уравнение переноса излучения в сочетании френелевскими граничными условиями) для моделирования радиационных полей; и биотепловое уравнение (нестационарное уравнение теплопроводности с конвективным слагаемым для учета переноса тепла за счет капиллярного кровотока) для моделирования тепловых полей.
Но эти программы имеют существенные недостатки по части принятых численных методов. Основным недостатком является использование только нелокального метода Монте-Карло (его недостатки указаны выше) для решения уравнения переноса излучения. Кроме того, для решения биотеплового уравнения используются громоздкие конечно-разностные методы,
В результате проведенного анализа в качестве исходных в данной работе приняты следующие модели и подходы к их реализации.
Для радиационных полей:
модель: квазистационарное уравнение переноса излучения с
френелевскими граничными условиями;
подход: метод Монте-Карло с использованием модифицированных
локальных оценок с конечной дисперсией.
Для тепловых полей:
модель: биотепловое уравнение с граничными условиями,
учитывающими взаимодействие биоткани с окружающей средой, и источниками тепла, генерируемыми лазерным излучением;
подход: метод взвешенных невязок для сведения биотеплового
уравнения к стационарному уравнению Гельмгольца и метод конечных элементов для решения этого уравнения.
ЦЕЛЬЮ данной работы является создание единого программного комплекса для моделирования радиационных и тепловых полей в гетерогенных биологических тканях, подвергаемых лазерному облучению.
Для достижения этой цели предполагается решить следующие ЗАДАЧИ:
1. Для математического ожидания случайных величин с бесконечной дисперсией предложить и исследовать статистические
оценки, имеющие лучшие характеристики сходимости и устойчивости по сравнению с выборочным средним.
На основе этих оценок разработать ускоренный локальный алгоритм метода Монте-Карло для решения уравнения переноса излучения. Конкретизировать его для моделирования полей лазерного излучения в мутных средах.
Разработать алгоритм решения биотеплового уравнения на основе методов конечных элементов и взвешенных невязок. Конкретизировать его для моделирования нестационарных тепловых полей в биологических тканях.
Реализовать алгоритмы в виде единого программного комплекса.
Провести тестирование комплекса и промоделировать ряд практических лазерных процедур,
НАУЧНАЯ НОВИЗНА
Разработаны новые оценки математического ожидания случайных величин с бесконечной дисперсией, имеющие большую скорость сходимости и большую устойчивость по отношению к выбросам по сравнению с выборочным средним.
На основе этих оценок разработан новый ускоренный локальный алгоритм метода Монте-Карло для расчета локальных радиационных характеристик полей лазерного излучения.
Предложена новая схема решения биотеплового уравнения, сочетающая метод взвешенных невязок и метод конечных элементов.
На основе принятых моделей и численных методов создан единый программный комплекс для расчета радиационных и тепловых полей лазерного излучения.
Впервые проведено моделирование радиационных и тепловых полей в тканях головного мозга, облучаемых лазерным излучением с
длиной волны 1064 нм и 805 нм; печени крысы, облучаемых лазерным излучением с длиной волны 1064 нм.
АПРОБАЦИЯ
Основные результаты работы докладывались и обсуждались на Всероссийской конференции молодых ученых «Физика в биологии и медицине» (Екатеринбург, 1999, 2001) (доклады были отмечены дипломами конференции); международной конференции по биомедицинской оптике «International Simposium in Biomedical Optics "BiOS'2000"» (Сан-Хосе, Калифорния, 2000), на семинаре Уральского научно-практического центра радиационной медицины (Челябинск, 2004), а также неоднократно на семинарах по биомедицинской оптике в Челябинском государственном университете (1997-2004 г.г.).
ПРАКТИЧЕСКАЯ ЗНАЧИМОСТЬ
Разработанный ускоренный локальный метод может быть применен для практических расчетов локальных радиационных характеристик полей как лазерного излучения, так и ионизирующего излучения.
Созданный программный комплекс предназначен для моделирования нестационарных радиационных и тепловых полей в гетерогенных биологических тканях, описываемых в виде двумерных осесимметричных многозонных сред. Он может быть использован для разработки и оптимизации лазерных хирургических операций, методов лазерной термотерапии и биостимуляции, для сравнительного исследования воздействия лазерного излучения с различными длинами волн; может служить основой математического обеспечения экспериментального определения оптических и теплофизических параметров биотканей.
3. Созданный программный комплекс может быть использован во всех организациях, занимающихся разработкой и оптимизацией лазерных медицинских технологий. В частности, он был использован в Челябинском государственном институте лазерной хирургии для оптимизации лазерных процедур на головном мозге; в МУЗ ГКБ №1 для оптимизации процедур лазерной термотерапии и исследования полей лазерного излучения в присутствии фотосенсибилизатора в процедурах фотодинамической терапии; в Челябинском государственном университете для математического обеспечения экспериментального определения оптических параметров биотканей, оптимизации метода экспериментального измерения температуры ткани.
ОСНОВНЫЕ ПОЛОЖЕНИЯ И РЕЗУЛЬТАТЫ, ВЫНОСИМЫЕ НА ЗАЩИТУ
Предложенные оценки математического ожидания случайных величин с бесконечной дисперсией обладают существенными преимуществами по сравнению с выборочным средним. Они имеют большую скорость сходимости и большую устойчивость по отношению к выбросам, к ним применим стандартный способ одновременного оценивания статистической погрешности (через оценивание дисперсии). Установленный критерий для выбора свободного параметра оценок обеспечивает близкое к оптимальному сочетание устойчивости и точности оценок.
Предложенная схема решения биотеплового уравнения абсолютно устойчива и имеет второй порядок сходимости как по времени, так и по пространству.
Созданный программный комплекс позволяет рассчитывать основные характеристики нестационарных радиационных и тепловых полей в биологических тканях, подвергаемых электромагнитному
облучению в оптическом диапазоне, в приближении двумерных гетерогенных многозонных сред.
4. Выполненные расчеты полей представляют практический интерес для лазерной хирургии. С их помощью разработаны и оптимизированы лазерные операции на головном мозге, щитовидной железе, печени и др.
ПУБЛИКАЦИИ
Результаты работы изложены в 9 статьях и 4 тезисах докладов на научных конференциях.
Решению поставленных задач посвящены 4 главы диссертации.
ПЕРВАЯ ГЛАВА содержит обзор известных моделей, методов и программ для моделирования радиационных и тепловых полей в биотканях, подвергающихся лазерному облучению.
ВТОРАЯ ГЛАВА посвящена решению задачи №1: разработке и исследованию статистических оценок математического ожидания одномерных неотрицательных случайных величин с бесконечной
дисперсией, функция распределения F(;e) = P{<*} которых удовлетворяет либо условию
l-F(x) = 0(x-a),x-»<»;a>l, (1)
либо его частному случаю
l~F(x) = cx'a +о(х~р},лг-*а=;1<ог<2, а*(3;с>0. (2)
Общая идея предлагаемого в настоящей работе метода оценивания математического ожидания заключается в использовании априорной информации о случайной величине , а именно, условия (1) или (2).
Построение новой ускоренной локальной оценки основано на статистическом оценивании характеристической функции g случайной
величины Щ и использовании асимптотических связей между
характеристической функцией и математическим ожиданием // = М^.
Эти связи установлены и сформулированы в виде теоремы для произвольного распределения, имеющего асимптотическое разложение степенного типа, включая частные случаи (1), (2), Этот результат используется для построения двух оценок математического ожидания: первая оценка применима для оценивания математического ожидания в случае (1) (и, разумеется, в его частном случае (2)), вторая - в случае (2). Обе оценки зависят от некоторого неотрицательного свободного параметра t > О.
Для этих оценок: получены асимптотические оценки скорости изменения смещения и дисперсии при f-»0, построен критерий устойчивости, предложены практические способы оценивания смещения и статистической погрешности.
Оценки численно исследованы на задаче о вычислении математического ожидания /г неотрицательной случайной величины с функцией распределения F вида:
П, х*\
1-F(x)~.
-3/2 л (3)
В этом случае Щ имеет бесконечную дисперсию и математическое ожидание /и = 3.
ТРЕТЬЯ ГЛАВА посвящена решению задач №2, 3: разработке численных алгоритмов расчета характеристик радиационных и тепловых полей лазерного излучения в гетерогенных биологических тканях, описываемых в виде двумерных осесимметричных многозонных сред.
Для моделирования радиационных полей приняты:
математическая модель: квазистационарная кинетическая модель переноса излучения с френелевскими граничными условиями;
численный подход: нелокальный и ускоренный локальный метод Монте-Карло, который включает оценку Калоса и полученные в Главе 2 оценки математического ожидания.
Для моделирования тепловых полей приняты:
математическая модель: модель биотеплового уравнения с граничными условиями, учитывающими взаимодействие биоткани с окружающей средой, и источниками тепла, генерируемыми лазерным излучением. Эта модель включает в качестве частных случаев keff и «bio heat transfer» модели, а также уравнение теплопроводности без конвекции;
численный подход: метод взвешенных невязок для сведения биотеплового уравнения к решению стационарного уравнения (типа Гельмгольца) на каждом временном шаге и метод конечных элементов для решения этого уравнения.
ЧЕТВЕРТАЯ ГЛАВА посвящена решению задач № 4,5: разработке и тестированию программного комплекса и моделированию с его помощью ряда практических лазерных процедур. Настоящая глава содержит:
общую характеристику программного комплекса (область применимости, входные параметры, рассчитываемые характеристики полей излучения, возможности комплекса);
результаты тестирования программного комплекса, которое включало сопоставление с известными аналитическими решениями, результатами других программ, а также результатами эксперимента;
исследование ускоренных локальных оценок для задач лазерной медицины;
результаты практического моделирования лазерных процедур: радиационных и тепловых полей в ткани мозга человека, подвергающейся облучению двух видов лазеров; в коже человека; в ткани печени крысы.
Численные методы расчета радиационных полей
Основное применение эта модель получила в работах по экспериментальному определению оптических параметров биотканей [65, 83], точнее параметра 2. Экспериментальные значения интенсивности на оси коллимированного пучка аппроксимируются экспоненциальной зависимостью, показатель которой и является искомым параметром 2.
Следует также отметить, что данная модель используется в нескольких работах и для расчета радиационных полей [154, 68, 116]. В частности, в работе [54] - для моделирования процесса фотокоагуляции сетчатки глаза и в [82, 101] - лазерной тепловой абляции (испарения твердых структур вещества) тканей зуба. Была принята плоская однородная геометрия среды, поверхность которой подвергается облучению тонким лазерным пучком. Оптические свойства среды описывались одним параметром 2, значение которого предполагалось постоянным и одинаковым для всей ткани. Полученные распределения I{z) использовались для расчета генерируемых лазерным излучением источников тепла, необходимых для моделирования тепловых полей (см. математические модели тепловых полей) в случае [54], и для расчета скорости абляционного процесса [82, 101].
Но широкого применения для биомедицинских задач эта модель не получила. Это связано, как было уже отмечено выше, с очень узкой областью применения этой модели: она не описывает большинство биотканей, являющихся сильнорассеивающими неоднородными средами; формула (1.1) позволяет рассчитывать только одномерные распределения интенсивности (в зависимости от глубины z), в то время как поля лазерного излучения сильно анизотропны и, Существуют различные модификации закона (1.1), позволяющие, в некотором смысле, расширить область применимости модели. В частности, в работе [53] предложена модификация (1.1) для учета рассеяния излучения при моделировании лазерной коагуляции. Принята двумерная осесимметричная геометрия среды, поверхность которой облучается лазерным пучком гауссова профиля. Предполагается, что основная характеристика поля излучения, интенсивность /(г), удовлетворяет уравнению вида: Здесь R френелевский коэффициент отражения; /0 — интенсивность падающего излучения; р - радиус пучка (расстояние, на котором интенсивность падает в е раз); л - сечение рассеяния. Эта модель, в принципе, позволяет рассчитывать уже двумерные радиационные поля, но она справедлива лишь для однородных сред и, следовательно, не является адекватной для большинства биотканей. Таким образом, несмотря на серьезные преимущества (аналитический вид, простота в реализации) модель Бугера-Ламберта-Бэра не является адекватной для задач биомедицинской оптики. Диффузионная модель Основной характеристикой поля излучения также является интенсивность /(г). Предполагается, что она удовлетворяет уравнению диффузии с поглощением и источниками (т.е. неоднородному уравнению Гельмгольца): Здесь S(г) — плотность мощности источника излучения в точке г, равная количеству энергии, испускаемой единицей объема среды с центром в точке г в единицу времени. Оптические свойства среды определяются двумя параметрами: D - коэффициент диффузии иссечение поглощения.
Данная модель, в принципе, может применяться к любой геометрии среды. Область ее применимости ограничивается слабоанизотропными полями и источниками излучения в среде. В отсутствие источников в среде (что характерно для лазерной медицины) слабоизотропные поля возникают в слабонеоднородных участках среды вдали от границ.
Эта модель имеет большую область применимости, чем предыдущая; содержит два параметра (D и 2а); основана на достаточно простом уравнении, для решения которого хорошо развиты численные методы. Но для задач лазерной медицины эта модель также не адекватна. Это связано, прежде всего, с сильной анизотропией лазерных пучков, что противоречит основным приближениям модели. К существенным недостаткам этой модели следует также отнести неопределенность граничных условий. Фактически диффузионное приближение, само по себе справедливо лишь для области, достаточно удаленной от границ и источников.
Диффузионная модель в биомедицинской оптике применяется в основном для расчета тепловых источников (плотность поглощенной мощности), используемых при моделировании тепловых полей. При этом используются следующие приближения.
Программные комплексы для расчета радиационных и тепловых полей
Если в качестве функции f принят полином вида (1.25) и весовыми функциями являются сами базисные функции Nx ..t,NK, получаем «метод Бубнова-Галеркина» [6]. В настоящее время предприняты попытки создания программных комплексов для расчета радиационных и тепловых полей биологических тканей, облучаемых лазером. Наиболее известными являются программные комплексы: LAT1S, MCML, LITCIT, HELIOS.
Программа MCML, разработанная авторами Lihong Wang (Optical Imaging Laboratory, Texas A&M University) и Steven L, Jacques (Oregon Medical Laser Center) [145-148], предназначена для расчета характеристик радиационных полей лазерного излучения в гетерогенной двумерной осесимметричной среде. Распространение излучения в среде рассматривается в рамках кинетической теории. Уравнение переноса излучения решается аналоговым методом Монте-Карло. Расчет тепловых полей в рамках данного комплекса не предусмотрен.
Программа HELIOS, разработанная Веассо СМ., Mordon S.R., Brunetaud J.M. (INSERM, France) [53], предназначена для моделирования процессов тепловой коагуляции в тканях, облучаемых лазером. Для расчета радиационных полей использован «расширенный» закон Бэра (1.2) для расчета пространственного распределения интенсивности лазерного излучения в цилиндрической системе координат (зависимость от угла не учитывается) с учетом рассеяния. Но такая модель (как было упомянуто выше) не является вполне обоснованной и справедлива только в приближении сильного поглощения и тонкого лазерного пучка. Расчет тепловых полей с учетом влияния кровотока основан на решении биотеплового уравнения конечно-разностным методом.
Наиболее полными являются программные комплексы: LATIS, LITCIT. Программные комплексы LATIS (LAser TISsue), разработанный коллективом авторов London R.A., GHnsky М.Е., Bailey D.S., Zimmerman G.B., Eder D.C. (Lawerence Livermore National Laboratory, California) [97] и LITCIT (Laser Induced Temperature Calculation In Tissue), разработанный коллективом авторов Roggan A., Muller G. (Laser-Medizin-Zentrum gGmbH, Berlin) [120] предназначены для расчета характеристик радиационных и тепловых полей лазерного излучения.
В данных комплексах принята двумерная гетерогенная осесимметричная модель среды. Моделирование радиационных полей основано на решении уравнения переноса излучения аналоговым методом Монте-Карло, моделирование тепловых полей - на решении биотеплового уравнения конечно-разностным методом. При расчете радиационных и тепловых полей LATIS учитывает следующие физические процессы: — «материальный»: основан на решении уравнений состояния, определяющего внутреннюю энергию и давление; фазовые переходы, химические процессы: протеиновая денатурация, коагуляция ткани; — «гидродинамический»: основан на учете упругих и пластических деформаций, эффектов кавитации и абляции. Учет всех перечисленных эффектов осуществляется изменением значений оптических и теплофизических параметров биотканей. Среди отечественных известна лишь одна программа для расчета радиационных и тепловых полей лазерного излучения. Это программа разработана коллективом Саратовского государственного университета под руководством проф. Тучина В.В. [39]. Принята двумерная гетерогенная осесимметричная модель среды. Моделирование радиационных полей основано на решении уравнения переноса излучения аналоговым методом Монте-Карло, моделирование тепловых полей - на решении стационарного уравнения теплопроводности методом конечных элементов. Учет кровотока не производится. Перечисленные комплексы являются наиболее полными на данный момент, но они не являются программным продуктом и предназначены только для использования самими авторами. К существенным недостаткам всех программ следует отнести применение только нелокальных методов расчета радиационных полей. Кроме того, программы не используют метод конечных элементов, который на наш взгляд, лучше для биомедицинских задач. В результате проведенного анализа в качестве исходных в данной работе приняты следующие модели и подходы к их реализации. Для радиационных полей: модель: квазистационарное уравнение переноса излучения с френелевскими граничными условиями; подход: метод Монте-Карло с использованием модифицированных локальных оценок с конечной дисперсией. Для тепловых полей: модель: биотешювое уравнение с граничными условиями, учитывающими взаимодействие биоткани с окружающей средой, и источниками тепла, генерируемыми лазерным излучением; подход: метод взвешенных невязок для сведения биотеплового уравнения к стационарному уравнению Гельмгольца и метод конечных элементов для решения этого уравнения.
Численный метод расчета радиационных полей
Заметим, что результат (2.52) сильнее простой асимптотической несмещенности оценок ц]у т.е. равенства Mr}} y(t) = Ь} (t) + о(1), /- 0, которое выполняется для любой оценки с математическим ожиданием, стремящимся к нулю при t -» 0. Можно сказать, что г) удовлетворяет условию асимптотической относительной несмещенности. Только оно гарантирует возможность уменьшением t безгранично уменьшать относительную систематическую погрешность оценки смещения. Оценки смещения г)п содержат помимо / новый свободный параметр я 1. Согласно (2.54) параметр а не следует выбирать близким к 1. С другой стороны при а »1 произведение at может стать слишком большим, не обеспечивающим близости математического ожидания г} к оцениваемому смещению. Целесообразно выбирать значение я в области нескольких единиц, скажем, я = 2. Применение оценок г} требует знания и параметра распределения у. В случае JMJ он определяется параметром порядка а (у=а-1), который часто, как отмечалось, может быть определен из общих представлений о конструкции случайной величины Щ. В случае fa значение параметра у, как правило, неизвестно. Но даже в отсутствие точного знания / применение оценок и имеет смысл. Используя в (2.51) вместо точного у какую-нибудь нижнюю оценку у (всегда можно положить у = а-1 т.к. уЄ(а-\,ї\), приходим к оценке f]jr, являющейся оценкой смещения сверху. Она является достаточно хорошей, поскольку отношение i?ur Mj,r во всех случаях не превосходит величины rf}a_l/rjil=(a \)/(aa 1-1), которая достаточно близка к I при рекомендуемых значениях параметра а и значениях параметра а не слишком близких к 1. Например, ;; a-i/ ;1 «2.5 при Нетрудно построить оценки смещения, и не предполагая известным параметр у. Эти оценки содержат комбинацию уже из трех оценок /л , вычисленных при разных значениях параметра t. Оценки более громоздки и имеют существенно большую дисперсию. На наш взгляд они практически менее интересны, и мы их рассматривать не будем. В случае, когда не выполняется ни одно из условий (2.39)-(2.41) или нельзя проверить их справедливость оценки rjn имеет смысл использовать для грубого оценивания смещения оценок /Uj. Известно, что характер разброса выборочного среднего % на разных выборках существенно различается в случаях D оо и Dlf = «=. В отличие от первого во втором случае Dif = оо и распределение if не нормально на сколь угодно больших выборках: оно имеет «тяжелый хвост», приводящий к неустойчивости по отношению к выбросам, т.е. большой (по сравнению с первым случаем) вероятности появления выборок, на которых далеко отклоняется от характерных значений Mlf. Для рассматриваемых случайных величин (2.1), (2.2) такая неустойчивость усугубляется с уменьшением параметра а, уже для значений а — 3/2 оценка становится из-за этого мало приемлемой для практики. Оценки jU; при любом D oo имеют конечную дисперсию и их распределения асимптотически нормальны. Поэтому они устойчивы по отношению к выбросам при достаточно больших объемах выборок п. Однако, как отмечалось, ц}- % при f- 0 на любой фиксированной выборке. Следовательно, в случае D =« уменьшение t должно приводить к ухудшению устойчивости іі} в отношении выбросов. В терминах распределений это объясняется следующим. При достаточно большом п распределение каждой оценки [i} близко к нормальному.
При уменьшении t (при фиксированном п) оно отклоняется от нормального и приближается к предельному распределению , т.е. к некоторому устойчивому закону, «хвост» которого намного «тяжелее» нормального. Качественными рассуждениями можно построить некий простой критерий для выбора значений t, обеспечивающих устойчивость fi} по выбросам. Рассмотрим случай а не слишком близких к 2, скажем а з/2, когда проблема выбросов особенно актуальна. В этом случае выбросы , как правило, связаны с выбросом одного из элементов выборки (1 ..,,1 ), соизмеримого с суммой всех остальных элементов [16]. Пусть для определенности t [л=М (к = 1,...,п -1) и Щп Хпц. Тогда (1 + X)fi, т.е. выброс - (л составляет около Я 100%. На той же выборке при tfi к 1 из (2.31), (2.32) получаем: следует ожидать устойчивости оценок jUj в отношение выбросов для распределений с параметром а, не слишком близким к 2. Для а близких к 2 этот критерий, по-видимому, менее точен, но и проблема выбросов здесь менее актуальна. Правое неравенство критерия обеспечивает близость смещенной оценки /и к несмещенной Щ на регулярных выборках (не дающих выбросов), что гарантирует не слишком большие, а часто достаточно малые смещения ц . Во всяком случае поиск оптимальных t следует производить в диапазоне значений, удовлетворяющих (2.55). В качестве неизвестного математического ожидания (л в (2.55) можно принять его оценку, pi или , вычисленную на малой (пробной) выборке.
Тестирование программного комплекса
Принятые в данной работе математические модели и численные методы были реализованы в программном комплексе, предназначенном для расчета радиационных и тепловых полей лазерного излучения в биотканях. Основные результаты этой главы были опубликованы в [1, 4], а также в материалах международной конференции по биомедицинской оптике «International Simposium in Biomedical Optics "BiOS 2000"» (Сан-Хосе, Калифорния, 2000) [93] и научно-информационном сборнике Лазерной Ассоциации РФ [3].
Цель; расчет радиационных и тепловых полей лазерного излучения в биологических тканях. Область применения1: лазерная медицина, дозиметрия, биооптика, теплофизика. Методы расчетов: радиационных полей: квазистационарное уравнение переноса излучения с френелевскими граничными условиями; метод Монте-Карло с использованием модифицированных локальных оценок с конечной дисперсией. тепловых полей: биотепловое уравнение с граничными условиями, учитывающими взаимодействие биоткани с окружающей средой, и источниками тепла, генерируемыми лазерным излучением; метод взвешенных невязок для сведения биотеплового уравнения к стационарному уравнению Гельмгольца и метод конечных элементов для решения этого уравнения. Типы излучения (частиц); квазифотоны. Типы взаимодействия: упругое рассеяние, поглощение. Источники излучения; произвольное пространственное и угловое распределение. Геометрия облучаемых объектов: двумерная осесимметричная модель ткани, представляющая собой набор однородных кольцевых зон прямоугольного сечения (рис. 4.1, совпадает с рис. 3.1). Каждая зона описывается набором оптических и теплофизических параметров: коэффициент рассеяния 2Д, коэффициент поглощения Иа, средний косинус угла рассеяния g, показатель преломления /г, толщина слоя /, удельная теплоемкость с, плотность р, коэффициент теплопроводности к. Рассчитываемые характеристики: пространственные распределения плотности поглощенной мощности, отраженного и прошедшего излучения; пространственное и временное распределение температуры внутри ткани. Сервис: интерактивный ввод данных, включая выбор задачи, ввод начальных данных, оптических и теплофизических параметров ткани; прерывание расчетов и продолжение ранее прерванных расчетов. Разработанный программный комплекс прошел ряд тестовых расчетов, включающих сравнение с известными аналитическими решениями, результатами других программ, а также результатами эксперимента. На рис. 4.2-4.7 приводятся отдельные результаты тестирования, подтверждающие надежность разработанного программного комплекса. Моделирование радиационных полей Для проверки работоспособности радиационной части программного комплекса были реализованы следующие задачи. 1. Сравнение с аналитическим решением. Рассчитывалась интенсивность излучения на оси тонкого цилиндрического пучка в бесконечной однородной слаборассеивающей среде. Согласно закону Бугера-Ламберта-Бэра в таком случае интенсивность нерассеянного излучения на оси пучка /(г) = /(z) с увеличением глубины z убывает по экспоненциальному закону с показателем 2 = 2 + 2а. Лазерный пучок: диаметр d = 0.01 см, мощность 1 Вт, интенсивность Параметры среды: 2в = 100.0см"1, 2. Сравнение с другими программами. Рассчитывалась поле плотности поглощенной мощности Qir z) в бесконечной однородной среде, облучаемой цилиндрическим лазерным пучком. Лазерный пучок: диаметр /=0.1 см, мощность 1 Вт. Параметры среды: 2Д= 10.0см"1, 2S = 1.0см"1, # = 0.9, « = 1.5. На рис. 4.3 представлено распределение Q(r,z) на различной глубине на оси пучка. Для сравнения использовались результаты расчетов MCML кода [147]. Моделирование тепловых полей 1. Сравнение с аналитическим решением. Рассчитывалось тепловое поле внутри однородного бесконечного цилиндра радиусом 1 см, помещенного в среду с фиксированной температурой О С. Боковая поверхность контактирует с окружающей средой, коэффициент теплообмена равен 10 Вт/(см2-К). Внешние источники тепла отсутствуют. Параметры среды: р = 1 г/см3, k = 1 Вт/(см К), с = 1 Дж/(г-К). Начальная температура равна 100 С Рис. 4.4 демонстрирует полное согласие результатов расчетов нашего комплекса с аналитическими значениями температуры для простой нестационарной тепловой задачи. 2. Сравнение с другими программами, В этом примере реализована более сложная задача, включающая расчет радиационных полей и тепловых полей с источниками тепла (плотность поглощенной мощности), вызванными лазерным облучением. Рассчитывалось пространственное распределение температуры в различные моменты времени внутри ограниченного однородного цилиндра на глубине 2 = 0.55 см, облучающегося цилиндрическим мононаправленным лазерным пучком. Верхний торец цилиндра (z = 0) теплоизолирован (нет теплообмена), остальные поверхности контактируют со средой постоянной температуры 20 С; коэффициент теплообмена равен 1 Вт/(см2 К). В начальный момент времени t = 0 температура цилиндра постоянна и равна 20 С. Лазерный пучок: цилиндрический профиль, диаметр d = 1см, мощность 5 Вт. Параметры среды: д=0.2см , 2, =50.0см-1, g = 0.95,n = L4, к = 4.8 мВт/(см К), /7 = 1.075 г/см3, с = 34882 Дж/(г К). Полученные результаты были сравнены с данными моделирования LITT кода [120], предоставленные авторами для сравнения (рис. 4.5).
Сравнение с экспериментом. На рис. 4.6 демонстрируются результаты сравнения наших численных результатов, аналитических и экспериментальных значений температуры для задачи о нагреве бесконечной однородной среды точечным источником тепла.
Измерения были выполнены с помощью собственной оригинальной установки для определения теплофизических параметров ткани: удельной теплоемкости с и коэффициента теплопроводности к ткани. В качестве объекта для исследования был выбран картофель размером около 10 см. Электрический тепловой источник размером около I мм был помещен в центр картофеля, температура измерялась с помощью термопар размером 0.3 мм [95].
Похожие диссертации на Математическое моделирование радиационных и тепловых полей в биологических тканях, подвергаемых лазерному облучению
