автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Математическое моделирование радиационных и тепловых полей в биологических тканях, подвергаемых лазерному облучению

На правах рукописи

Аникина Алла Степановна

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАДИАЦИОННЫХ И ТЕПЛОВЫХ ПОЛЕЙ В БИОЛОГИЧЕСКИХ ТКАНЯХ, ПОДВЕРГАЕМЫХ ЛАЗЕРНОМУ ОБЛУЧЕНИЮ

05.13.18 - математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Челябинск -2004

Работа выполнена в ГОУВПО Челябинском государственном университете

НАУЧНЫЙ РУКОВОДИТЕЛЬ:

доктор физико-математических наук, профессор

Лаппа Александр Владимирович

ОФИЦИАЛЬНЫЕ ОППОНЕНТЫ:

доктор физико-математических наук, профессор

Учайкин Владимир Васильевич кандидат физико-математических наук, профессор

Заляпин Владимир Ильич

ВЕДУЩАЯ ОРГАНИЗАЦИЯ:

Государственное образовательное учреждение высшего профессионального образования «Уральский государственный университет»

Защита диссертации состоится « 29 » декабря 2004 г.

в 11.00 часов на заседании Диссертационного Совета Д 212.296.02 при Государственном образовательном учреждении высшего профессионального образования «Челябинский государственный университет» по адресу: 454021, г. Челябинск, ул. Бр. Кашириных, 129.

С диссертацией можно ознакомиться в библиотеке Челябинского государственного университета Автореферат разослан «д[в» ноября 2004 г.

диссертационного совета:

Ученый секретарь

АКТУАЛЬНОСТЬ ТЕМЫ

Настоящая работа посвящена моделированию процессов переноса излучения и тепла в биологических тканях, подвергающихся интенсивному лазерному облучению. Такие исследования актуальны для лазерной терапии и хирургии, биооптики, теплофизики. В медицине лазерные источники излучения и оптоволоконная техника позволили создать целый ряд новых эффективных медицинских технологий, существенно уменьшающих сроки лечения, его стоимость, риск и тяжесть осложнений. К таким технологиям, в частности, относятся термотерапия (лазерная гипертермия), каналирование и перфорация мягких и костных тканей, фотодинамическая терапия, технологии, основанные на термокоагуляции и термоабляции тканей.

Во всех перечисленных процедурах главными физическими агентами воздействия, вызывающими конечный медицинский эффект, являются излучение и тепло. Следовательно, для исследования воздействия лазерного излучения на биоткани необходимо, прежде всего, иметь информацию о радиационных и тепловых полях, возбуждаемых в процессе облучения.

Метод математического компьютерного моделирования получения такой информации естественным образом делится на два этапа: первый представляет собой расчет характеристик радиационного поля, второй — расчет характеристик тепловых полей с вычисленными на предыдущем этапе источниками тепла.

В настоящее время, предприняты попытки создания компьютерных программ для решения этой задачи. Наиболее полными являются программные комплексы LATIS [London R.A. et al. // Appl Optics, 1997, vol. 36] и LITCIT [Roggan A., Muller G. // Proceedings of SPIE, 1994, vol. 2327]. Но они являются чисто лабораторным продуктом и не доступны для сторонних пользователей. Достоинством этих комплексов является использование наиболее адекватных математических моделей: кинетическая модель переноса излучения (т.е. кинетическое уравнение переноса излучения в сочетании френелевскими граничными условиями) для моделирования радиационных полей; и биотепловое уравнение (нестационарное уравнение теплопроводности с конвективным слагаемым для учета переноса тепла за счет капиллярного кровотока) для моделирования тепловых полей.

Но эти программы имеют существенные недостатки по части принятых численных методов. Для моделирования радиационных полей используется нелокальный метод Монте-Карло, который дает значения радиационных характеристик, усредненные по некоторым областям фазового пространства, что может приводить к большим систематическим погрешностям в случаях больших градиентов радиационных полей, весьма характерных для практических лазерных процедур. В то же время в теории методов Монте-Карло для расчета полей ионизирующего излучения разработаны локальные алгоритмы, позволяющие рассчитывать нужные характеристики радиационного поля в заданных точках пространства. Наиболее известной и часто применяемой на практике является локальная оценка Калоса [Kalos М.Н. // NucI Sci and Engng, 1963, vol. 16], позволяющая рассчитывать плотность потока частиц (гамма-квантов, нейтронов) в фиксированной точке.

Расчет плотности потока сводится к вычислению математического ожидания случайной величины с бесконечной дисперсией. Стандартное оценивание математического ожидания предполагает использование «универсальной» оценки (не использующей никакой информации о случайной величине): выборочного среднего. Для случайных величин с бесконечной дисперсией выборочное среднее обладает рядом недостатков, существенно затрудняющих его практическое использование: пониженная скорость сходимости; плохое качество сходимости (велика вероятность больших отклонений (выбросов) от оцениваемой величины). Перечисленные негативные свойства выборочного среднего являются «продолжением» его уникального достоинства: применимости к любым случайным величинам с конечным математическим ожиданием. Для построения более качественных оценок естественно отказаться от такой универсальности, попытаться учесть, по возможности более полно, имеющуюся априорную информацию о распределении осредняемой случайной величины. В настоящее время предложен ряд, так называемых, ускоренных локальных методов, но, несмотря на это, проблема далека от окончательного решения.

Для моделирования тепловых полей в указанных программах используются громоздкие конечно-разностные методы, реализация которых

особенно трудна для задач, в которых параметры уравнения зависят от температуры.

Более интересными, в практическом смысле, являются схемы, приводящие к пошаговому решению уравнения по времени: температурное поле в следующий момент времени находится из температурного поля в предыдущий момент времени. В частности, метод взвешенных невязок для биотеплового уравнения приводит к решению на каждом временном шаге уравнения, совпадающего с уравнением типа Гельмгольца. Для различных приложений, включая биомедицинские, имеющих дело со сложной геометрией и гетерогенными средами, для решения этого уравнения хорошо зарекомендовал себя метод конечных элементов.

В настоящее время не существует полного алгоритма решения биотеплового уравнения, сочетающего преимущества этих методов: они реализованы только в виде отдельных программ.

ЦЕЛЬЮ данной работы является создание единого программного комплекса для моделирования радиационных и тепловых полей в гетерогенных биологических тканях, подвергаемых лазерному облучению.

Для достижения этой цели предполагается решить следующие ЗАДАЧИ:

1. Для математического ожидания случайных величин с бесконечной дисперсией предложить и исследовать статистические оценки, имеющие лучшие характеристики сходимости и устойчивости по сравнению с выборочным средним.

2. На основе этих оценок разработать ускоренный локальный алгоритм метода Монте-Карло для решения уравнения переноса излучения. Конкретизировать его для моделирования полей лазерного излучения в мутных средах.

3. Разработать алгоритм решения биотеплового уравнения на основе методов конечных элементов и взвешенных невязок. Конкретизировать его для моделирования нестационарных тепловых полей в биологических тканях.

4. Реализовать алгоритмы в виде единого программного комплекса.

5. Провести тестирование комплекса и промоделировать ряд практических лазерных процедур.

НАУЧНАЯ НОВИЗНА

1. Разработаны новые оценки математического ожидания случайных величин с бесконечной дисперсией, имеющие большую скорость сходимости и большую устойчивость по отношению к выбросам по сравнению с выборочным средним.

2. На основе этих оценок разработан новый ускоренный локальный алгоритм метода Монте-Карло для расчета локальных радиационных характеристик полей лазерного излучения.

3. Предложена новая схема решения биотеплового уравнения, сочетающая метод взвешенных невязок и метод конечных элементов.

4. На основе принятых моделей и численных методов создан единый программный комплекс для расчета радиационных и тепловых полей лазерного излучения.

5. Впервые проведено моделирование радиационных и тепловых полей в тканях головного мозга, облучаемых лазерным излучением с длиной волны 1064 нм и 805 нм; печени крысы, облучаемых лазерным излучением с длиной волны 1064 нм.

АПРОБАЦИЯ

Основные результаты работы докладывались и обсуждались на Всероссийской конференции молодых ученых «Физика в биологии и медицине» (Екатеринбург, 1999, 2001) (доклады были отмечены дипломами конференции); международной конференции по биомедицинской оптике «International Simposium in Biomedical Optics "BiOS'2000"» (Сан-Хосе, Калифорния, 2000), на семинаре Уральского научно-практического центра радиационной медицины (Челябинск, 2004), а также неоднократно на семинарах по биомедицинской оптике в Челябинском государственном университете (1997-2004 г.г.).

ПРАКТИЧЕСКАЯ ЗНАЧИМОСТЬ

1. Разработанный ускоренный локальный метод может быть применен для практических расчетов локальных радиационных характеристик полей как лазерного излучения, так и ионизирующего излучения.

2. Созданный программный комплекс предназначен для моделирования нестационарных радиационных и тепловых полей в гетерогенных биологических тканях, описываемых в виде двумерных осесимметричных многозонных сред. Он может быть использован для разработки и оптимизации лазерных хирургических операций, методов лазерной термотерапии и биостимуляции, для сравнительного исследования воздействия лазерного излучения с различными длинами волн; может служить основой математического обеспечения экспериментального определения оптических и теплофизических параметров биотканей.

3. Созданный программный комплекс может быть использован во всех организациях, занимающихся разработкой и оптимизацией лазерных медицинских технологий. В частности, он был использован в Челябинском государственном институте лазерной хирургии, для оптимизации лазерных процедур на головном мозге; в МУЗ ГКБ №1 для оптимизации процедур лазерной термотерапии и исследования полей лазерного излучения в присутствии фотосенсибилизатора в процедурах фотодинамической терапии; в Челябинском государственном университете для математического обеспечения экспериментального определения оптических параметров биотканей, оптимизации метода экспериментального измерения температуры ткани.

ОСНОВНЫЕ ПОЛОЖЕНИЯ И РЕЗУЛЬТАТЫ, ВЫНОСИМЫЕ НА ЗАЩИТУ

1. Предложенные оценки математического ожидания случайных величин с бесконечной дисперсией обладают существенными преимуществами по сравнению с выборочным средним. Они имеют большую скорость сходимости и большую устойчивость по отношению к выбросам, к ним применим стандартный способ одновременного оценивания статистической погрешности (через оценивание дисперсии). Установленный критерий для выбора свободного

параметра оценок обеспечивает близкое к оптимальному сочетание устойчивости и точности оценок.

2. Предложенная схема решения биотеплового уравнения абсолютно устойчива и имеет второй порядок сходимости как по времени, так и по пространству.

3. Созданный программный комплекс позволяет рассчитывать основные характеристики нестационарных радиационных и тепловых полей в биологических тканях, подвергаемых электромагнитному облучению в оптическом диапазоне в приближении двумерных гетерогенных многозонных сред.

4. Выполненные расчеты полей представляют практический интерес для лазерной хирургии. С их помощью разработаны и оптимизированы лазерные операции на головном мозге, щитовидной железе, печени и др.

ПУБЛИКАЦИИ

Результаты работы изложены в 9 статьях и 4 тезисах докладов на научных конференциях.

СТРУКТУРА И ОБЪЕМ ДИССЕРТАЦИИ

Диссертация изложена на 152 страницах и состоит из введения, четырех глав, заключения и списка цитированной литературы. Диссертация содержит 30 рисунков и 4 таблицы.

СОДЕРЖАНИЕ РАБОТЫ

Во введении обоснована актуальность темы и сформулированы основные задачи работы. Изложено краткое содержание диссертации по главам.

Первая глава содержит обзор известных моделей, методов и программ для моделирования радиационных и тепловых полей в биотканях, подвергающихся лазерному облучению. В результате проведенного анализа в качестве исходных в данной работе приняты следующие модели и подходы к их реализации.

Для радиационных полей: • модель: квазистационарное уравнение переноса излучения с френелевскими граничными условиями;

• подход: метод Монте-Карло с использованием модифицированных локальных оценок с конечной дисперсией.

Для тепловых полей:

• модель: биотепловое уравнение с граничными условиями, учитывающими взаимодействие биоткани с окружающей средой, и источниками тепла, генерируемыми лазерным излучением;

• подход: метод взвешенных невязок для сведения нестационарного биотеплового уравнения к уравнению типа Гельмгольца на каждом временном' шаге и метод конечных элементов для решения этого уравнения.

Во Второй главе предложены и исследованы статистические оценки математического ожидания одномерных неотрицательных случайных величин с бесконечной дисперсией, функция распределения которых

удовлетворяет либо условию:

Условие (1), в принципе, допускает и распределения с конечной дисперсией , но этот случай неинтересен, и в настоящей работе всегда считается = оо.

К этому типу случайных величин принадлежит и локальная оценка Калоса: она является неотрицательной случайной величиной типа (2) с а — 3/2.

Общая идея предлагаемого в настоящей работе метода оценивания математического ожидания заключается в использовании априорной информации о случайной величине а именно, условия (1) или (2). Построение новой ускоренной локальной оценки основано на статистическом оценивании характеристической функции случайной величины и использовании асимптотических связей между характеристической функцией и математическим ожиданием , установленных в следующей теореме.

(1)

либо его частному случаю:

\-F{x) = cx'a+ о{х'1),х^юо;\<а<2, а</3;с>0. (2)

Теорема 1. Пусть функция распределения F неотрицательной случайной величины £ имеет вид:

N

1 - F(x) = Y,cix~"' + fM> x>0;N = 0,1,...; 1 < or, <... <aN < 2\cj * 0, (3) y-i

о

где г(*) = фГвА,),*-»°°; a0e(l,2), £ =0.

Тогда характеристическая функция этой случайной величины представима в виде: g{t) = 1 + it/л + ¿сД {t)\t\"> + cp(t), (4)

где *e(0 = Re4, +i\mAasign{t), Аа =i\z'a{e,z-\)dz, ^(0 = o(|f),f-»0. (5)

0

В двух частных случаях, когда l)s(x) = 0(x'fi), х —>к>, f}>aN и 2) = о{х'р), х -> оо, р > aN имеем:

О(И'), /?<2,

1) ?(/) = ■ 0(/21п|/|), /? = 2, /—>0; 2) <p(t) = Kt2+o(t2),0> 2,

ее

где К = - ^e(x)xdx, |А'|<со.

4Г)> р<2>

Kt2 +o(t2),0>2,

Непосредственно при построении оценок используются два частных результата этой теоремы: первая оценка использует результат теоремы при N = 0, вторая-при N = 1.

Обе оценки могут быть записаны одной формулой:

Mj* MjW--£zj,(t), 7 = 1.2;

/ ч _ Im(exp) -1) - S,Ba Re(exp) -1) XjA')- t '

(6)

(7)

где (#,,...,£,) - независимая выборка реализаций />0 - свободный параметр,

Оценка ^ применима для оценивания в случае (1) (и, разумеется, в его частном случае (2)); оценка цг применима в случае (2).

Эти оценки обладают целым рядом достоинств по сравнению со стандартной оценкой математического ожидания — выборочным средним. В отличие от последнего они имеют конечную дисперсию при и, как следствие, они

обладают большей скоростью убывания статистической погрешности, к ним применим стандартный способ ее оценивания, их сходимость имеет более регулярный характер. Платой за эти преимущества является смещенность и даже несостоятельность оценок при любом значении параметра 1>0. Но эти недостатки исчезают в асимптотике при что делает оценки интересными с

практической точки зрения.

Теоретическое исследование позволило установить следующие факты, существенные для практического применения оценок

• При 1-^0 смещение о ц е нр -г^)аниченно убывает, а их дисперсия неограниченно растет.

• Более точно асимптотическое поведение смещения и дисперсии оценок при

сформулировано в следующих теоремах для неотрицательной случайной величины £с функцией распределения ¥, удовлетворяющей условию (1) (/' = 1) и условию (2) (/' = 2). Теорема 2.

• При D£ = oo оценки fij, в отличие от выборочного среднего, устойчивы по

отношению к выбросам при не слишком малых Практический критерием

устойчивости является неравенство:

1/и/*<«/«1/// (8)

В качестве неизвестного математического ожидания ц можно принять его

оценку: или £ , вычисленную на малой (пробной) выборке.

• Оценка предпочтительней в условиях применимости обеих оценок. Преимущества оценки //, по сравнению с 4 проявляются только при = оо.

Для построения корректных алгоритмов вычисления смещений оценок fit нужны более точные результаты об их асимптотическом поведении. Теорема 2 устанавливает такой результат только в одном случае: для при Для

расширения этого результата необходимо усилить требования к функции распределения F случайной величины Введем в рассмотрение усиленный вариант условия (1):

1 - F(x) = сх'а + о(х~" ),jc-»co;l<a<2;c>0, (9)

Введем также два усиленных варианта условия (2):

1 - F(x) = сх'а + с^'" + о(х'р ),х->°о-,1<а<р<2;с,с{*0, (10)

или

F(x) = cx'a +о{х'"),х-^<х>-,\<а<2<р-, с>0. (Ц)

Каждое из введенных условий гарантирует бесконечность

В условиях (9)-(11) асимптотики смещений являются точно степенными, что позволяет построить (неединственным образом) статистические оценки смещения на тех же выборках, что и сами оценки Построенная в настоящей работе

оценка смещения имеет вид:

Г)J At) = , г > о, д > \,у е [а-1,1] (12)

Здесь у = 1, / = я -1 вслучае(9), j = 2,у = @-\ вслучае(Ю), j = 2, у = \ вслучае (11).

Оценка смещения 7/,г(') позволяет оценивать систематическую погрешность вычисления математического ожидания одновременно с оцениванием самого Ц. Она содержит помимо параметра / новый свободный параметр а. Асимптотические свойства оценки смещения 7у,Д') сформулированы в следующей теореме.

Теорема 4

Пусть выполняется любое из условий (9)-(11). Тогда для оценки (12) справедливо:

М77Л,(/) = (М^(0- + о(1)), а/ О,

Численное исследование оценок ц^ было выполнено для задачи о бросании наудачу точки в шар единичного радиуса. Здесь £ = р'г (р — расстояние от точки до центра шара) имеет функцию распределения вида:

[1, х<\,

(13)

В этом случае £ имеет дисперсию 0£ = со и математическое о ж и д а н^ЗЭ т а задача отражает характерные особенности локальной оценки Калоса для расчета характеристик радиационных полей: функция распределения (13) имеет такую же степенную асимптотику

На рис. 1 приведены значения относительного смещения (значения и относительного стандартного отклонения оценок (значения

!¡л). Представлены: точные значения; оцененные статистически (со

статистикой и = 106); асимптотические значения. Видно, что смещения (по

модулю) обеих оценок монотонно уменьшаются, а стандартные отклонения

увеличиваются с уменьшением параметра в области Область действия

асимптотик обеих оценок примерно одинакова. Смещения и стандартные

отклонения, оцененные статистически, практически совпадают с точными при 1-2

в области меньших значений начинает сказываться статистическая погрешность.

Рис 1. Относительные смещения (а) и стандартные отклонения (б) оценок ¡1Х (1) и цг (2) как функции параметра I. Открытые маркеры - точные значения; сплошные - статистическая оценка, линии - асимптотики.

Рис. 2 иллюстрирует значительно более регулярный характер сходимости предлагаемых оценок к математическому ожиданию в области выполнения критерия (8) при увеличении объема выборки п в сравнении с выборочным средним 4 . Приведены относительные отклонения (/^ —//)/// для двух значений

параметров / = Ю"5 и / = 10~2 и (£-//)/>. Видно, что выборочное среднее в данном расчете имеет два выброса: малый, около 10%, при и«104 и большой, до 30%, при ия5-105, влияние которого распространяется дальше и = 106.

Обе предлагаемые оценки при «средних» значениях параметра / = 10~2, удовлетворяют критерию (8) и имеют совершенно устойчивый характер сходимости, без выбросов. Хорошо видно смещение оценки Его значение

(8%) достигается быстро, практически к , и далее значение оценки уже

существенно не меняется. Смещение оценки цг (0.5%) не заметно, ее отклонение

от математического ожидания уменьшаются до 1% к , а после

оно становится устойчиво меньше.

При очень малых значениях свободного параметра смещения

оценок существенно меньше: оно незаметно даже для оценки но характер сходимости в обоих случаях значительно хуже. Это соответствует нарушению критерия устойчивости (8): значения оценок с уменьшением t приближаются к выборочному среднему, и начинают проявляться все его негативные свойства.

Рис. 2. Относительное отклонение от математического ожидания выборочного среднего £ и оценок //,, /Лг для / = 10~5 и 10~2 в зависимости от объема выборки п. Использовалась одна последовательность реализаций £,,£2,„.дм всех оценок Пунктирными линиями обозначено точное относительное смещение оценок.

Третья глава посвящена разработке алгоритмов расчета характеристик радиационных и тепловых полей лазерного излучения на основе принятых математических моделей и подходов к их реализации (см. Глава 1).

В качестве модели гетерогенной среды принята двумерная осесимметричная модель, представляющая собой набор однородных изотропных кольцевых зон прямоугольного сечения со своими оптическими и теплофизическими параметрами (рис. 3).

Рис.3 Модель среды.

Основной рассчитываемой характеристикой радиационных полей в данной работе является плотность поглощенной мощности которая представляет

собой тепловой источник, генерируемый лазерным излучением, и используется при моделировании тепловых полей. Эта величина является 'локальной характеристикой, и для ее расчета в работе реализованы нелокальный и

ускоренный локальный методы Монте-Карло. Вычисление 0(г*) в обоих методах сводится к вычислению математического ожидания некоторого случайного функционала <Ц на траекториях частиц, распространение которых рассматривается как марковский процесс.

В рамках нелокального метода функционал имеет вид:

^ ~(~Дг|{о> г, й Дг'

и его математическое ожидание вычисляется как выборочное среднее:

(14)

где - число построенных траекторий частиц; - координата точки поглощения частицы на траектории, Дг, |Дг| - область, содержащая г* и ее объем, Р - мощность падающего излучения.

Ускоренный локальный метод сочетает локальную оценку Калоса и построенные в Главе 2 оценки математического ожидания (_/ = 1,2). Случайная

величина (оценка Калоса) определяется на аналоговых траекториях частиц следующим образом:

где - координата столкновения и направление движения частицы

перед ним; N - число столкновений частицы на траектории; £1* = (г* — Г()/|г* —

- оптическое расстояние между точками

Г = + коэффициент (макроскопическое сечение) поглощения и

рассеяния.

Перенос тепла за счет теплопроводности с учетом капиллярного кровотока в среде О описывается биотепловым уравнением с начальным и граничными условиями, учитывающими перенос тепла на внутренней Г' и внешней Г границах, и источниками тепла, генерируемыми лазерным излучением:

Здесь Т = 7*(г,/) —7^; Т"(г,/) -температура среды в точке г в момент времени V, Ть- температура крови в сосудах; к —коэффициент теплопроводности среды;

- удельная теплоемкость при постоянном давлении и плотность среды; О" = всьиь; Сь,1)ь — удельная теплоемкость при постоянном давлении и «объемный расход» кровотока в сосудах; в - параметр, характеризующий теплообмен между сосудом и тканью; - плотность поглощенной мощности лазерного

излучения (вычислена на предыдущем этапе); ТА,кл, Тв,кв - температура и коэффициент теплопроводности в смежных зонах А и В (с общей границей Г'); П, ПАВ - внешняя нормаль к Г и нормаль к Г' (направление от зоны А к зоне В); т — признак граничного условия, заданного на некотором подмножестве границы Г. Надлежащим выбором параметра т возможен учет следующих физических условий на внешней границе Г: - фиксированная температура

тепловой поток т — 3 - теплообмен по закону Ньютона: к — коэффициент теплоотдачи, Т3 — температура окружающей среды вблизи границы Г. Условия (19) обеспечивают физические требования на внутренних границах: непрерывность температуры и плотности тока тепла.

Численный подход для реализации тепловой модели сочетает идеи метода взвешенных невязок и метода конечных элементов.

Метод, основанный на идее метода взвешенных невязок, используется для сведения нестационарной задачи к решению стационарной на каждом временном шаге: решение параметризуется линейной функцией; ее единственный параметр находится из условия ортогональности невязки к весовой функции В

результате на каждом временном шаге решается стационарное

неоднородное уравнение Гельмгольца:

V(*V74r))-a7Xr) + £(r) = 0. (20)

Здесь Т{г) = Т(г,т2) - неизвестное температурное поле в момент времени г2, e(r) = r[e(r,rJ) + e(r,r1)]/2 + a7'(r,r1) + v(wr(r,rl)), T(r,zt) - уже известное

температурное поле в предыдущий момент времени г, (на первом шаге оно совпадает с начальным условием), к = тк/2, a.—cp + a'cll, a — cp — axll.

Граничные условия на каждом временном шаге имеют одинаковый вид и совпадают с условиями (18)-(19).

Полученная в результате схема решения нестационарного уравнения относится к неявным с параметром 1/2. В [Калиткин Н.Н. Численные методы. 1978] установлено, что данная схема имеет второй порядок точности по времени и абсолютно устойчива.

Для решения уравнения (20) используется метод конечных элементов. Идея метода заключается в нахождении функции Т из множества допустимых функций U, удовлетворяющей условию: %{T) = mfz(l>). В качестве U принято

множество непрерывных представителей удовлетворяющих условию

(у-Г,)|г=0,т.е. условию (18) при m = 1.

В этом случае удалось построить функционал, минимизирующая функция которого автоматически удовлетворяет исходному уравнению (20) и всем граничным условиям. Этот функционал имеет вид:

Для численного решения полученной вариационной задачи каждая зона среды дополнительно разбивалась на треугольные элементы с неизвестными значениями искомой функции Г в вершинах треугольника. В качестве аппроксимирующей функции внутри элемента был принят линейный полином.

Полученная в результате конечно-элементной аппроксимации численная схема относится к конформным схемам, подробно изученным в [Сьярле Ф. Метод конечных элементов для эллиптических задач, 1980]. Согласно результатам этой работы данная схема имеет второй порядок точности по пространственным координатам.

Четвертая глава посвящена разработке программного комплекса на основе полученных в Главе 3 алгоритмов расчета радиационных и тепловых полей лазерного излучения и моделированию с его помощью ряда практических лазерных процедур.

Программный комплекс предназначен для моделирования радиационных и тепловых полей излучения оптического диапазона в мутных средах и может быть использован для разработки и оптимизации лазерных хирургических операций, методов лазерной термотерапии и биостимуляции, может служить основой математического обеспечения экспериментального определения оптических и теплофизических параметров биотканей.

В качестве модели среды принята двумерная осесимметричная модель, представляющая собой набор однородных кольцевых зон прямоугольного сечения (рис. 3) со своими оптическими и теплофизическими параметрами. Каждая зона описывается набором оптических и теплофизических параметров: коэффициент рассеяния коэффициент поглощения Е0, средний косинус угла рассеяния (параметр индикатрисы рассеяния, в качестве которой принята традиционная для биооптики функция Хени-Гринштейна), показатель преломления удельная теплоемкость плотность коэффициент

теплопроводности

Моделирование осуществляется в два этапа. Первый этап представляет собой расчет характеристик радиационных полей. Основными являются следующие характеристики: пространственные распределения плотности поглощенной мощности, плотности тока отраженного и прошедшего излучения. Второй этап -моделирование нестационарных тепловых полей с известными источниками тепла, вычисленными на предыдущем этапе.

Программный комплекс прошел полномасштабные тестовые расчеты, включающие сопоставление с известными аналитическими решениями, результатами других подобных программ (рис. 4), а также результатами эксперимента. Полученные результаты показывают хорошую работоспособность комплекса. Программный комплекс был использован для моделирования

радиационных и тепловых полей в ткани мозга человека; в коже человека; в ткани

печени крысы и др.

равна 20 °С Параметры среды <: = 4 8мВт/(см К), /> = 11 г/см', с = 3 49 Дж/(гК)

расчеты, маркеры - результаты ЦТСГГ кода

Рис 4 Температура в различные моменты времени внутри ограниченного

однородного цилиндра на глубине 2 = 0 55 см, облучающегося цилиндрическим мононаправленным лазерным пучком Верхний торец цилиндра (г=0) теплоизолирован (нет теплообмена), остальные поверхности контактируют со средой температура, которой постоянна и равна 20 °С, коэффициент теплообмена равен 1 Вт/(см2 К) В начальный момент времени температура цилиндра постоянна и £„ = 02 см"1, 2, = 50 0 см"1, £ = 0 95, л = 14, Сплошная линия - собственные

На рис б приведены результаты применения ускоренных локальных оценок для типичной задачи лазерной медицины. Рассчитывалась интенсивность рассеянного лазерного излучения от единичного мононаправленного

цилиндрического лазерного пучка равномерной интенсивности (рис. 5).

Оптические параметры среды =1.29 см'1, 5^=63.4 СМ §;=0 87. В качестве точного значения

математического ожидания принято 4, рассчитанное на пробной выборке (я = 105):

Видно, что выборочное среднее £ в данном расчете имеет большой выброс до 100% при я»104, влияние которого распространяется дальше п = 10®. Обе оценки при «средних» значениях параметра / = 1(Г2, удовлетворяют критерию устойчивости оценок ^ при и = 104 и имеют совершенно устойчивый характер сходимости, без выбросов. Смещение обеих оценок практически одинаково, его значение (5%) достигается к и = 103, и далее уже существенно не меняется. При меньших значениях параметра смещения оценок имеют практически

такое же значение (5%), но характер сходимости в обоих случаях значительно

Рис 5 Схема облучения

хуже: оценки полностью повторяют выброс. Этот случай соответствует нарушению левой части критерия (8).

Рис 6 Относительное отклонение от математического ожидания выборочного среднего £ и оценок р2 для {= 10~3 и 1(Г2 в зависимости от объема выборки п Использовалась одна последовательность реализаций ...для всех оценок

На рис. 8-9 приведено распределение плотности поглощенной мощности и стационарные изотермы в полубесконечной среде белого вещества головного мозга человека, облучаемой лазерными пучками различной длины волны (805 нм и 1064 нм). Образец ткани представлял собой однородный цилиндр (рис 7) Теплообмен на границах не учитывался (все границы - теплоизолированы) Оптические и теплофизические параметры ткани для соответствующих длин волн приведены в таблице [Тучин В.В. Лазеры и волоконная оптика в биомедицинских исследованиях. Саратов, 1998].

Начальная температура для всей ткани равна 37 °С. Лазерный пучок в обоих случаях имел равномерное распределение плотности мощности, диаметр 0.04 см, мощность 1 Вт. Режим облучения — непрерывный.

-Ч

0.04 см

,2 см

/У

Рис 7 Схема облучения

Длина волны Ъа, 1/см 1/см г п Аг, 103 Вт/см К р, г/см3 с, Дж/(г К)

1064 нм 32 469 0 0 87 1 45 6 48 1 4 18

805 нм 02 400 0 09 1 38

яям рив нгяк мкрямя я«чш ,

Рис. 8. Пространственное распределение плотности поглощенной мощности в белом веществе мозга человека, облучаемого лазерным пучком (слева - ?.=805 нм, справа -Х=1064 нм).

2.9 •

2

, в-

0.9 -0.0-

лагерный н5чок ¡¡у^^о-Г .иирный щчок пучка, хЮ"1 см

Рис. 9. Стационарные изотермы в белом веществе головного мозга человека, облучаемого лазерным пучком (слева- Х=805 нм, справа- ^.=1064 нм).

ОСНОВНЫЕ РЕЗУЛЬТАТЫ И ВЫВОДЫ

Поставленные в работе задачи (см. стр. 5) выполнены. Цель, сформулированная на стр. 5, достигнута: создан единый действующий программный комплекс для моделирования радиационных и нестационарных тепловых полей лазерного излучения в мутных средах, подвергаемых облучению оптического диапазона. С помощью выполненных расчетов полей лазерного излучения разработаны и оптимизированы лазерные операции на головном мозге, щитовидной железе, печени и др.

Основные результаты диссертации изложены в следующих работах:

1. Аникина А.С., Лаппа А.В. Моделирование температурных полей в биологических тканях, облучаемых лазером // Лазерные технологии в медицине: сборник научных работ. - Челябинск, 1998. - Вып. 1. - С. 125-127.

1 I з

расстояние огося

2. Аникина А. С. Моделирование радиационно-температурных полей в биологических тканях, облучаемых лазером // Физика в биологии и медицине: сборник научных работ. - Екатеринбург, 1999. - С. 26-27.

3. Аникина А.С. Моделирование температурных полей в биологических тканях при взаимодействии с лазерным излучением // Студент и научно-технический прогресс: тезисы научных докладов. - Челябинск, 1999. - С. 10-11.

4. Аникина А.С, Камалов В.А., Лаппа А.В. Программный комплекс для моделирования радиационно-тепловых полей в биологических тканях, облучаемых лазером // Лазерные технологии в медицине: сборник научных работ. - Челябинск, 1999. - Вып. 2. - С. 187-198.

5. Lappa A.V., Anikina A.S., Kamalov V.A. A new computer code for calculation of radiation and heat fields in laser-irradiated tissues. // Laser-Tissue Interaction XI: Photochemical, Photothermal, and Photomechanical: Proceedings of SPIE; D.D. Duncan, J.O. Hollinger, S.L. Jacques; Eds. - 2000. - V. 3914. - P. 502-511.

6. Аникина А.С, Лаппа А.В. Компьютерное моделирование радиационно-тепловых полей в биологических тканях, облучаемых лазером // Использование лазеров для диагностики и лечения заболеваний: научно-информационный сборник (приложение к бюллетеню «Лазер-Информ» Лазерной Ассоциации). - Москва, 2001. - Вып. 3. - С. 46-52.

7. Аникина А.С Исследование радиационно-тепловых полей биологических тканей, облучаемых лазером, методами компьютерного моделирования // Физика в биологии и медицине: сборник научных работ. - Екатеринбург, 2 0 01.-С 4-6.

8. Лаппа А.В., Рогальский Ю.К., Аникина А.С Возмущение спектра диффузного отражения биологической ткани поглощающими экзогенными агентами. Там же, С. 29-30.

9. Аникина А.С, Лаппа А.В. Исследование радиационно-тепловых полей в биологических тканях, облучаемых лазером, методами компьютерного моделирования // Лазерные технологии в медицине: сборник научных работ. -Челябинск, 2001.-Вып. 3.-С. 176-181.

10. Privalov V.A., Lappa A.V., Anikina A.S. et al.Clinical Trials of a New Chlorin Photosensitizer for Photodynamic Therapy of Malignant Tumors // Optical Methods

»24156

for Tumor Treatment and Detection. Mechanisms and Techniques in Photodynamic Therapy XI: Proceedings of SPIE; T.J. Dougherty; Ed. - 2002. - Vol.4612. -P. 178-190.

П.Привалов В.А., Лаппа А.В., Аникина А.С. и др. Клинические результаты лазерной фотодинамический терапии злокачественных новообразований с использованием нового сенсибилизатора Радахлорин // Новые технологии и фундаментальные исследования в медицине: сборник научных трудов. -Челябинск, 2002. - С. 101-104.

12. Лаппа А.В., Бахвалов Е.В., Аникина А.С. Метод характеристических функций в оценивании математического ожидания случайных величин с бесконечной дисперсией// Известия Челябинского научного центра. — Челябинск, 2004. -Вып. 2.-С. 1-6.

13. Лаппа А.В., Бахвалов Е.В., Аникина А.С. Оценки математического ожидания случайных величин с бесконечной дисперсией в «методе характеристических функций». Свойства, погрешность, устойчивость // Известия Челябинского научного центра. - Челябинск, 2004. - Вып. 3. - С. 6-11.

Подписано в печать 25.11.04 Формат 60x84 1/16. Бумага офсетная. Печать офсетная. Усл. печ.л. 1,4. Уч.-издл. 1,0 Тираж 100 экз. Заказ 357 Бесплатно

ГОУВПО «Челябинский государственный университет»

454021 Челябинск, ул. Бр. Кашириных, 129 Полиграфический участок Издательского центра ЧелГУ 454021 Челябиинск, ул. Молодогвардейцев, 576

Введение.

Глава 1. Состояние вопроса (обзор).

1.1. Математические модели радиационных полей.

1.2. Численные методы расчета радиационных полей.

1.3. Математическая модель тепловых полей.

1.4. Численные методы расчета тепловых полей.

1.5. Программные комплексы для расчета радиационных и тепловых полей.

Глава 2. Метод характеристических функций в оценивании математического ожидания случайных величин с бесконечной дисперсией.

2.1. Асимптотики характеристической функции.

2.2. Оценки математического ожидания.

2.3. Смещение.

2.4. Дисперсия.

2.5. Погрешности.

2.6. Выбросы.

2.7. Численное исследование оценок.

Глава 3. Моделирование радиационных и тепловых полей.

3.1. Математическая модель радиационных полей.

3.2. Численный метод расчета радиационных полей.

3.3. Математическая модель переноса тепла.

3.4. Решение биотеплового уравнения.

Глава 4. Программный комплекс для моделирования радиационных и тепловых полей. Тестирование и приложения

4.1. Общие характеристики комплекса.

4.2. Тестирование программного комплекса.

4.3. Численное исследование локальных оценок.

4.4. Практическое моделирование лазерных процедур.

Настоящая работа посвящена моделированию процессов переноса излучения и тепла в биологических тканях, подвергающихся интенсивному лазерному облучению.

Такие исследования актуальны для лазерной терапии и хирургии, биооптики, теплофизики. В медицине лазерные источники излучения и оптоволоконная техника позволили создать целый ряд новых эффективных медицинских технологий, существенно уменьшающих сроки лечения, его стоимость, риск и тяжесть осложнений. К таким технологиям, в частности, относятся термотерапия (лазерная гипертермия), каналирование и перфорация мягких и костных тканей, фотодинамическая терапия, технологии, основанные на термокоагуляции и термоабляции тканей. Все они нашли широкое применение в онкологии, отоларингологии, дерматологии, урологии, нейрохирургии, кардиохирургии, косметологии и других отраслях медицины.

Главными физическими агентами воздействия во всех перечисленных процедурах, вызывающими конечный медицинский эффект, являются излучение и тепло. Следовательно, для исследования воздействия лазерного1 излучения на биоткани необходимо, прежде всего, иметь информацию о радиационных и тепловых полях, возбуждаемых в процессе облучения.

Одним из способов получения такой информации являются приборные методы. Большинство используемых методов являются проникающими: регистрирующие датчики помещаются внутри ткани. Преимуществом таких методов является их точность и надежность. Но они имеют существенные недостатки:

• введение датчиков создает дополнительную травму ткани, которая часто недопустима по медицинским показаниям, т.к. ведет к повреждению жизненно важных объектов (сосуды, нервные узлы);

• получаемая информация, относится лишь к малому числу точек.

Существующие непроникающие приборные методы дистанционные и контактно-поверхностные) позволяют получать лишь, поверхностное распределение нужных характеристик. Ведутся разработки и для непроникающих измерений радиационных и температурных полей внутри тканей, но пока они далеки от внедрения в клиническую практику.

Альтернативным методом получения информации о радиационных и тепловых полях является метод математического компьютерного моделирования. Он, в принципе, позволяет получать пространственные распределения характеристик; проводить исследования «без повреждения» исследуемого объекта и существенных затрат, которых требует эксперимент.

Разработка средств компьютерного моделирования включает 3 задачи:

1. создание математической модели (исходные уравнения, граничные условия, параметры);

2. разработка алгоритмов реализации модели (аналитические, численные, статистические методы);

3. создание компьютерной программы, реализующей эти алгоритмы.

К настоящему времени в литературе предложен ряд моделей, численных алгоритмов и компьютерных программ для моделирования радиационных и тепловых полей. Проведенный анализ литературы (Глава 1) позволяет сделать следующие выводы.

1. Наиболее адекватной математической моделью переноса лазерного излучения в мутных средах является кинетическая модель, т.е. кинетическое уравнение переноса излучения в сочетании с законами Френеля, описывающими граничные условия. Именно эта модель принята в настоящей работе.

Наиболее адекватные математические модели переноса тепла в биотканях основаны на нестационарном уравнении теплопроводности и отличаются разными способами учета конвективного переноса тепла за счет кровотока.

Учет теплового эффекта кровотока в больших сосудах (вены, артерии) осуществляется явным образом. Из-за сложной геометрии кровеносных систем такие модели не являются универсальными, они разрабатываются отдельно для конкретных биотканей.

Учет капиллярного кровотока осуществляется двумя методами: «keff »-метод и «bio heat transfer»-MeTOfl. а) Идея «keff »-метода заключается в учете теплового эффекта капиллярного кровотока в коэффициенте теплопроводности. Такой способ не требует изменения вида исходного уравнения (уравнения теплопроводности) и введения дополнительных параметров, но на наш взгляд, является не вполне обоснованным физически. б) «Bio heat 1гапзГег»-метод - физически более обоснован. Учет капиллярного кровотока осуществляется здесь введением в уравнение теплопроводности дополнительного конвективного слагаемого, пропорционального разности температуры среды и температуры крови. Такой метод содержит дополнительные параметры и требует решения уже более сложного уравнения, называемого в мировой литературе «bio heat transfer equation». Поскольку в отечественной литературе перевод этого термина пока не закрепился, мы примем свой: «биотепловое» уравнение. В стационарном случае оно совпадает с уравнением Гельмгольца.

Сопоставления, проведенные в ряде работ, позволяют заключить, что, несхмотря на простоту, первая модель может применяться в случае достаточно мелких капилляров. Вторая модель позволяет учитывать капилляры практически любых размеров.

В настоящей работе в качестве тепловой модели выбрана комбинация этих моделей, включающая в качестве частных случаев keff и «bio heat transfer» модели, а также уравнение теплопроводности без конвекции.

2. Наиболее точным методом реализации кинетической модели, позволяющей учитывать сложную геометрию среды и источников излучения, является метод Монте-Карло. Во всех работах по биооптике реализован, так называемый, нелокальный метод Монте-Карло, что является существенным недостатком при расчете радиационных полей. Такие методы дают значения радиационных характеристик, усредненные по некоторым областям фазового пространства, что может приводить к большим систематическим погрешностям в случаях больших градиентов радиационных полей, весьма характерных для практических лазерных процедур. В то же время в теории методов Монте-Карло для расчета полей ионизирующего излучения разработаны локальные алгоритмы, позволяющие рассчитывать нужные характеристики радиационного поля в заданных точках пространства. Наиболее известной и часто применяемой на практике является локальная оценка Калоса (см. Глава 1 «Математические модели радиационных полей»), позволяющая рассчитывать плотность потока частиц (гамма-квантов, нейтронов) в фиксированной точке. Платой за это преимущество является бесконечность дисперсии этих оценок и, как следствие, пониженная скорость сходимости, плохое качество сходимости (велика вероятность больших отклонений от оцениваемой величины), что существенно затрудняет практическое использование этих оценок. Несмотря на то, что для преодоления этого недостатка предложен ряд, так называемых, ускоренных локальных методов, проблема далека от окончательного решения.

Литературный анализ показал, что численные методы для принятой тепловой модели не достаточно развиты. Во всех работах используются громоздкие конечно-разностные методы. В то же время эта модель может быть реализована любым численным методом, разработанным для уравнения теплопроводности, решению которого посвящена обширная литература.

Для решения стационарного уравнения теплопроводности наибольшее применение получили: метод конечных элементов; конечно-разностный метод; метод взвешенных невязок (решение параметризуется некоторой функцией; ее параметры находятся из условия ортогональности невязки к заданным весовым функциям). Для различных приложений, включая биомедицинские, имеющих дело со сложной геометрией и гетерогенными средами, хорошо зарекомендовал себя метод конечных элементов.

Для решения нестационарного уравнения теплопроводности, в основном, применяются методы, основанные на его сведении к системе более простых уравнений. Такое сведение осуществляется разными способами. Наиболее распространенным является сведение к системе обыкновенных дифференциальных уравнений (по времени) применением какого-нибудь численного алгоритма к пространственной части уравнения, т.е. к стационарному уравнению теплопроводности, в предположении известной производной по времени. Для решения стационарного уравнения используются те же методы: метод конечных элементов, конечно-разностный метод. Для решения системы обыкновенных дифференциальных уравнений используются в основном методы численного интегрирования Рунге-Кутта, конечно-разностный метод, метод взвешенных невязок.

Недостатком этого способа является громоздкость: приходится численно решать две системы. Реализация этого способа особенно трудна для задач, в которых параметры уравнения зависят от температуры.

Другим способом решения нестационарного уравнения является сведение к системе дифференциальных уравнений применением какого-нибудь конечно-разностного алгоритма к временной части. Наиболее интересными, в практическом смысле, являются схемы, приводящие к пошаговому решению уравнения по времени: температурное поле в следующий момент времени находится из температурного поля в предыдущий момент времени. В таких схемах на каждом временном шаге приходится решать одно стационарное уравнение. Часто - это уравнение типа Гельмгольца. Этот подход проще и естественнее предыдущего, не возникает проблем, связанных с зависимостью параметров от температуры.

Для численной реализации биотеплового уравнения, на наш взгляд, последний способ является более естественным: биотепловое уравнение и в стационарном случае является уравнением типа Гельмгольца.

3. В настоящее время в мире предприняты попытки создания компьютерных программ для моделирования радиационных и тепловых полей. Наиболее полными являются программные комплексы: LATIS (London R.A. et al., Lawerence Livermore National Laboratory, California) и LITCIT (Roggan A. et al., Laser-Medizin-Zentrum gGmbH, Berlin).

Но они являются чисто лабораторным продуктом и не доступны для сторонних пользователей. Достоинством этих комплексов является использование наиболее адекватных математических моделей: кинетическая модель переноса излучения (т.е. кинетическое уравнение переноса излучения в сочетании френелевскими граничными условиями) для моделирования радиационных полей; и биотепловое уравнение (нестационарное уравнение теплопроводности с конвективным слагаемым для учета переноса тепла за счет капиллярного кровотока) для моделирования тепловых полей.

Но эти программы имеют существенные недостатки по части принятых численных методов. Основным недостатком является использование только нелокального метода Монте-Карло (его недостатки указаны выше) для решения уравнения переноса излучения. Кроме того, для решения биотеплового уравнения используются громоздкие конечно-разностные методы.

В результате проведенного анализа в качестве исходных в данной работе приняты следующие модели и подходы к их реализации. Для радиационных полей:

• модель: квазистационарное уравнение переноса излучения с френелевскими граничными условиями;

• подход: метод Монте-Карло с использованием модифицированных локальных оценок с конечной дисперсией.

Для тепловых полей:

• модель: биотепловое уравнение с граничными условиями, учитывающими взаимодействие биоткани с окружающей средой, и источниками тепла, генерируемыми лазерным излучением;

• подход: метод взвешенных невязок для сведения биотеплового уравнения к стационарному уравнению Гельмгольца и метод конечных элементов для решения этого уравнения.

ЦЕЛЬЮ данной работы является создание единого программного комплекса для моделирования радиационных и тепловых полей в гетерогенных биологических тканях, подвергаемых лазерному облучению.

Для достижения этой цели предполагается решить следующие ЗАДАЧИ:

1. Для математического ожидания случайных величин с бесконечной дисперсией предложить и исследовать статистические оценки, имеющие лучшие характеристики сходимости и устойчивости по сравнению с выборочным средним.

2. На основе этих оценок разработать ускоренный локальный алгоритм метода Монте-Карло для решения уравнения переноса излучения. Конкретизировать его для моделирования полей лазерного излучения в мутных средах.

3. Разработать алгоритм решения биотеплового уравнения на основе методов конечных элементов и взвешенных невязок. Конкретизировать его для моделирования нестационарных тепловых полей в биологических тканях.

4. Реализовать алгоритмы в виде единого программного комплекса.

5. Провести тестирование комплекса и промоделировать ряд практических лазерных процедур.

НАУЧНАЯ НОВИЗНА

1. Разработаны новые оценки математического ожидания случайных величин с бесконечной дисперсией, имеющие большую скорость сходимости и большую устойчивость по отношению к выбросам по сравнению с выборочным средним.

2. На основе этих оценок разработан новый ускоренный локальный алгоритм метода Монте-Карло для расчета локальных радиационных характеристик полей лазерного излучения.

3. Предложена новая схема решения биотеплового уравнения, сочетающая метод взвешенных невязок и метод конечных элементов.

4. На основе принятых моделей и численных методов создан единый программный комплекс для расчета радиационных и тепловых полей лазерного излучения.

5. Впервые проведено моделирование радиационных и тепловых полей в тканях головного мозга, облучаемых лазерным излучением с длиной волны 1064 нм и 805 нм; печени крысы, облучаемых лазерным излучением с длиной волны 1064 нм.

АПРОБАЦИЯ

Основные результаты работы докладывались и обсуждались на Всероссийской конференции молодых ученых «Физика в биологии и медицине» (Екатеринбург, 1999, 2001) (доклады были отмечены дипломами конференции); международной конференции по биомедицинской оптике «International Simposium in Biomedical Optics "BiOS'2000"» (Сан-Хосе, Калифорния, 2000), на семинаре Уральского научно-практического центра радиационной медицины (Челябинск, 2004), а также неоднократно на семинарах по биомедицинской оптике в Челябинском государственном университете (1997-2004 г.г.).

ПРАКТИЧЕСКАЯ ЗНАЧИМОСТЬ

1. Разработанный ускоренный локальный метод может быть применен для практических расчетов локальных радиационных характеристик полей как лазерного излучения, так и ионизирующего излучения.

2. Созданный программный комплекс предназначен для моделирования нестационарных радиационных и тепловых полей в гетерогенных биологических тканях, описываемых в виде двумерных осесимметричных многозонных сред. Он может быть использован для разработки и оптимизации лазерных хирургических операций, методов лазерной термотерапии и биостимуляции, для сравнительного исследования воздействия лазерного излучения с различными длинами волн; может служить основой математического обеспечения экспериментального определения оптических и теплофизических параметров биотканей.

3. Созданный программный комплекс может быть использован во всех организациях, занимающихся разработкой и оптимизацией лазерных медицинских технологий. В частности, он был использован в Челябинском государственном институте лазерной хирургии для оптимизации лазерных процедур на головном мозге; в МУЗ ГКБ №1 для оптимизации процедур лазерной термотерапии и исследования полей лазерного излучения в присутствии фотосенсибилизатора в процедурах фотодинамической терапии; в Челябинском государственном университете для математического обеспечения экспериментального определения оптических параметров биотканей, оптимизации метода экспериментального измерения температуры ткани.

ОСНОВНЫЕ ПОЛОЖЕНИЯ И РЕЗУЛЬТАТЫ, ВЫНОСИМЫЕ НА ЗАЩИТУ

1. Предложенные оценки математического ожидания случайных величин с бесконечной дисперсией обладают существенными преимуществами по сравнению с выборочным средним. Они имеют большую скорость сходимости и большую устойчивость по отношению к выбросам, к ним применим стандартный способ одновременного оценивания статистической погрешности (через оценивание дисперсии). Установленный критерий для выбора свободного параметра оценок обеспечивает близкое к оптимальному сочетание устойчивости и точности оценок.

2. Предложенная схема решения биотеплового уравнения абсолютно устойчива и имеет второй порядок сходимости как по времени, так и по пространству.

3. Созданный программный комплекс позволяет рассчитывать основные характеристики нестационарных радиационных и тепловых полей в биологических тканях, подвергаемых электромагнитному облучению в оптическом диапазоне, в приближении двумерных гетерогенных многозонных сред.

4. Выполненные расчеты полей представляют практический интерес для лазерной хирургии. С их помощью разработаны и оптимизированы лазерные операции на головном мозге, щитовидной железе, печени и др.

ПУБЛИКАЦИИ

Результаты работы изложены в 9 статьях и 4 тезисах докладов на научных конференциях.

Решению поставленных задач посвящены 4 главы диссертации.

ПЕРВАЯ ГЛАВА содержит обзор известных моделей, методов и программ для моделирования радиационных и тепловых полей в биотканях, подвергающихся лазерному облучению.

ВТОРАЯ ГЛАВА посвящена решению задачи №1: разработке и исследованию статистических оценок математического ожидания одномерных неотрицательных случайных величин £ с бесконечной дисперсией, функция распределения ^(л;) = Р{^<л;} которых удовлетворяет либо условию l-F(x) = о{х~а}, д:оо; а > 1, (1) либо его частному случаю l-F{x) = cx-a +о(х'р), х-*°с\1<а<2, а*/3;с> 0. (2)

Общая идея предлагаемого в настоящей работе метода оценивания математического ожидания заключается в использовании априорной информации о случайной величине а именно, условия (1) или (2). Построение новой ускоренной локальной оценки основано на статистическом оценивании характеристической функции g случайной величины £ и использовании асимптотических связей между характеристической функцией и математическим ожиданием (Л = М£.

Эти связи установлены и сформулированы в виде теоремы для произвольного распределения, имеющего асимптотическое разложение степенного типа, включая частные случаи (1), (2). Этот результат используется для построения двух оценок математического ожидания: первая оценка применима для оценивания математического ожидания в случае (1) (и, разумеется, в его частном случае (2)), вторая - в случае (2). Обе оценки зависят от некоторого неотрицательного свободного параметра t > 0.

Для этих оценок: получены асимптотические оценки скорости изменения смещения и дисперсии при / —>0, построен критерий устойчивости, предложены практические способы оценивания смещения и статистической погрешности.

Оценки численно исследованы на задаче о вычислении математического ожидания [г неотрицательной случайной величины £ с функцией распределения F вида:

1,

-3/2 - (3)

X , х>1.

В этом случае £ имеет бесконечную дисперсию и математическое ожидание /г = 3.

ТРЕТЬЯ ГЛАВА посвящена решению задач №2, 3: разработке численных алгоритмов расчета характеристик радиационных и тепловых полей лазерного излучения в гетерогенных биологических тканях, описываемых в виде двумерных осесимметричных многозонных сред.

Для моделирования радиационных полей приняты:

• математическая модель: квазистационарная кинетическая модель переноса излучения с френелевскими граничными условиями;

• численный подход: нелокальный и ускоренный локальный метод Монте-Карло, который включает оценку Калоса и полученные в Главе 2 оценки математического ожидания.

Для моделирования тепловых полей приняты:

• математическая модель: модель биотеплового уравнения с граничными условиями, учитывающими взаимодействие биоткани с окружающей средой, и источниками тепла, генерируемыми лазерным излучением. Эта модель включает в качестве частных случаев keff и «bio heat transfer» модели, а также уравнение теплопроводности без конвекции;

• численный подход: метод взвешенных невязок для сведения биотеплового уравнения к решению стационарного уравнения (типа Гельмгольца) на каждом временном шаге и метод конечных элементов для решения этого уравнения.

ЧЕТВЕРТАЯ ГЛАВА посвящена решению задач №4,5: разработке и тестированию программного комплекса и моделированию с его помощью ряда практических лазерных процедур. Настоящая глава содержит:

• общую характеристику программного комплекса (область применимости, входные параметры, рассчитываемые характеристики полей излучения, возможности комплекса);

• результаты тестирования программного комплекса, которое включало сопоставление с известными аналитическими решениями, результатами других программ, а также результатами эксперимента;

• исследование ускоренных локальных оценок для задач лазерной медицины;

• результаты практического моделирования лазерных процедур: радиационных и тепловых полей в ткани мозга человека, подвергающейся облучению двух видов лазеров; в коже человека; в ткани печени крысы.

ЗАКЛЮЧЕНИЕ

1. Разработаны и реализованы статистические оценки математического ожидания случайных величин с бесконечной дисперсией. Эти оценки обладают существенными преимуществами по сравнению с выборочным средним: они имеют большую скорость сходимости и большую устойчивость по отношению к выбросам, к ним применим стандартный способ одновременного оценивания статистической погрешности (через оценивание дисперсии).

Для этих оценок:

• установлен критерий для выбора свободного параметра, который обеспечивает близкое к оптимальному сочетание устойчивости и точности оценок;

• построены оценки смещения, которые позволяют оценивать систематическую погрешность вычисления математического ожидания (л одновременно с оцениванием самого [л.

2. На основе этих оценок разработан ускоренный локальный алгоритм метода Монте-Карло для решения уравнения переноса лазерного излучения в мутных средах.

3. Разработан алгоритм решения биотеплового уравнения, сочетающий метод взвешенных невязок и метод конечных элементов. Схема абсолютно устойчива и имеет второй порядок сходимости как по времени, так и по пространству.

4. Разработанные алгоритмы реализованы в виде программного комплекса для расчета характеристик нестационарных радиационных и тепловых полей в биологических тканях, подвергаемых электромагнитному облучению в оптическом диапазоне в приближении двумерных гетерогенных многозонных сред. Проведенные тестовые исследования позволяют судить о хорошей работоспособности комплекса.

5. С помощью комплекса проведено моделирование радиационных и тепловых полей в тканях головного мозга, облучаемых лазерным излучением с длиной волны 1064 нм и 805 нм; печени крысы, облучаемых лазерным излучением с длиной волны 1064 нм. Выполненные расчеты полей представляют практический интерес для лазерной хирургии. С их помощью разработаны и оптимизированы лазерные операции на головном мозге, щитовидной железе, печени и др.

1. Аникина А.С. Моделирование радиационно-температурных полей в биологических тканях, облучаемых лазером // Физика в биологии и медицине: сборник научных работ Уральской конференции молодых ученых. Екатеринбург, 1999. - С. 26-27.

2. Антюфеев B.C., Михайлов Г.А., Назаралиев М.А. Модификации локальных оценок с учетом осевой симетрии в задачах атмосферной оптики. // Вероятностные методы решения задач математической физики. Новосибирск: Изд. ВЦ СО АН СССР, 1971, С. 26-43.

3. Беляев Н.М., Рядно А.А. Методы теории теплопроводности. М.: Высшая школа, 1982. Т. 1,2. - 327 е., 304 с.

4. Бондаревский И .Я., Бордуновский В.Н., Астахова Л.В. Опыт применения высокоинтенсивного лазерного излучения при операциях на печени (экспериментальное исследование) // Лазерные технологии в медицине: Сборник научных работ. -Челябинск, 1999.-С. 114-121.

5. Гаплагер Р. Метод конечных элементов. Основы. М.: Мир, 1984. -428 с.

6. Деклу Ж. Метод конечных элементов. Пер. с франц. Б.И.Квасова; Под ред. Н.Н. Яненко. М.: Мир, 1976. 95 с.

7. Дэйвид Г. Порядковые статистики. М.: Наука, 1979. 336 с.

8. Ермаков С.М. Метод Монте-Карло и смежные вопросы. М.: Наука., 1975.-470 с.

9. Зарубин B.C. Инженерные методы решения задач теплопроводности. М.: Энергоатомиздат, 1983. 328 с.

10. Зенкевич О., Морган К. Конечные элементы и аппроксимация. М.: Мир, 1986.-318 с.

11. Зенкевич О.С., Тейлор P.JL, Ту Дж.М. Метод конечных элементов в технике. М.: Мир. 1975. - 543 с.

12. Золотарев В.М. Одномерные устойчивые распределения. М.: Наука, 1983.-304 с.

13. Ибрагимов И.А., Линник Ю.В. Независимые и стационарно связанные величины. М.: Наука, 1965. 524 с.

14. Исимару А. Распространение и рассеяние волн в случайно-неоднородных средах. М.: Мир, 1981. -Т. 1,2.- 280 е., 317 с.

15. Капиткин Н.Н. Численные методы. М.: Наука, 1978. 511 с.

16. Кольчужкин A.M., Учайкин В.В. Введение в теорию прохождения частиц через вещество. М.: Атомиздат, 1978. — 256 с.

17. Ландсберг Г.С. Оптика. М.: Высшая школа, 1976. 926 с.

18. Лаппа А.В., Бахвалов Е.В., Аникина А.С. Метод характеристических функций в оценивании математического ожидания случайных величин с бесконечной дисперсией // Известия Челябинского научного центра. Челябинск, 2004. Вып. 2. С. 1-6.

19. Лыков А.В. Теория теплопроводности. М.: Высшая школа, 1967.600 с.24,25,26