Трехмерная модель распространения света в биологических тканях Текст научной статьи по специальности «Математика»

2010 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Сер. 11 Вып. 3

СОВРЕМЕННЫЕ МЕДИЦИНСКИЕ ТЕХНОЛОГИИ

УДК 535.36;535.34;519.245

А. Ю. Сетейкин, И. В. Красников, М. С. Павлов

ТРЕХМЕРНАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ СВЕТА В БИОЛОГИЧЕСКИХ ТКАНЯХ

Амурский Государственный Университет, Благовещенск

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

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

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

— простота в определении, принадлежности текущего положения пакета фотонов, подобласти разбиения;

— простота формы границы области разбиения;

— точность разбиения первоначальной среды.

© А. Ю. Сетейкин, И. В. Красников, М. С. Павлов, 2010

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

Р{8 < Б} = 1 - ехр(-№Б) (1)

£ = Р{8 < Б} (2)

£ = 1 - ехр(-^Б) (3)

Б = - 1п(1 - е)/^, (4)

где е —вероятность того, что расстояние между двумя взаимодействиями пакета фотонов со средой будет меньше ожидаемого, Б — ожидаемое расстояние между двумя взаимодействиями пакета фотонов со средой, Р(в < Б) —вероятность того, что расстояние между двумя взаимодействиями пакета фотонов со средой будет меньше ожидаемого значения; ^ —коэффициент полного ослабления.

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

Рис. 1. Построение трехмерной сетки.

Применение метода Монте-Карло для моделирования распространения фотонов в мутных средах дает довольно гибкий подход к решению. Этот метод моделирует «случайный ход» фотонов в среде, которая обладает поглощением и рассеянием. Метод основан на наборе законов, которые управляют движением фотонов в ткани [3, 4]. Есть три основных величины, определяющих траекторию движения фотона: длина свободного пробега до взаимодействия со средой, угол отклонения и азимутальный угол. На границе раздела двух сред, или двух подобластей с разными оптическими свойствами внутри одной среды фотон может отразиться или пройти через нее. Законы распределения фотонов выражены через распределение вероятности для соответствующих величин в случаях движения, преломления и отражения на границе. Распространение света при моделировании методом Монте-Карло строго, но в достаточной степени наглядно. Этот метод является статистическим и требует вычисления на ЭВМ для подсчета распространения большого числа фотонов.

В данной модели фотон рассматривается как нейтральная частица. Для удобства вычислений можно ввести понятие пакета фотонов, совокупность фотонов, которые движутся вместе. Каждому пакету фотонов присваивается статистический вес, первоначально равный единице. Этот вес показывает, какая часть фотонов содержится в пакете фотонов на конкретный момент времени. При каждом взаимодействии пакета фотонов со средой, часть фотонов из пакета будет поглощаться, соответственно статистический вес пакета будет уменьшаться [7].

Длину свободного пробега фотона удобно задавать в безразмерных величинах длины. Для каждого пробега он задается как 1п(1 — £), где £ — случайная величина. Безразмерная величина длины представляет собой произведение длины на коэффициент полного ослабления ¿в ■ Тогда расстояние, пройденное внутри элемента, определяется, как

(Х2,У2,22)

в = J ^(х,у,г)3,а. (5)

(х1,у1,г1)

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

в = & ■ V7(Х2 - XI)2 + (У2 - У1)2 + (¿2 - zí)2 (6)

а = ^ +2 МХЪУЪ ^ • ^2 - х,)2 + Ы ~ Уг? + - (7)

соответственно.

Когда фотон рассеивается, траектория фотона отклоняется на угол в в интервале [0,п]. Как уже упоминалось, функция плотности распределения вероятности, которая определяет вероятность рассеяния фотонов, это фазовая функция. Фазовая функция обычно определяется эмпирически. Наиболее часто используемой в оптике тканей является функция Хени—Гринштейна. Генерируется случайное число £ из [0, 1]. Для простоты вычислений произведем замену

^ = еов(в). (8)

и

Тогда значение ^ можно выразить как

если g > 0

2Z — 1, если g = 0

и угол отклонения в для пакета определяется как cos 1(^) [3, 4]. Выбор ф также зависит от случайного числа:

ф = 2п£.

(10)

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

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

Внутреннее отражение, Д(а), считается из закона Френели

Закон преломления связывает значение углов падения а и преломления «1 с коэффициентами преломления сред П1 и П2, от которых зависит движение фотона:

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

— если £ < Д(а), то фотон отражен;

— если £ > Д(а), то фотон прошел.

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

После произведения N запусков фотонов для каждого узла сетки мы имеем некоторую выборку (XI, Х2,..., Хп) статистических весов, полученных средой в его окрестности. Для этих величин производится статистическая обработка. Наиболее важными являются математическое ожидание и выборочная дисперсия.

(11)

1, arcsin

ni sin а = ni sin ai.

(12)

0.05 0.04 0.03 0.02 0.01

-0.01 -0.02 --0.03 --0.04-0.05---1-1---L - J- -1-1---1 -1

-0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.05

Рис. 2. Поперечное сечение луча.

Пусть вес всех N фотонов соответствует некоторому значению энергии q. Тогда математическое ожидание полученного веса узлом сетки:

n

X = (l/n)£ Xi, (13)

i=i

где n — число поглощений в окрестности рассматриваемого узла; Xi — это статистический вес, поглощенный в окрестности рассматриваемого i-го узла.

Значение математического ожидания может быть использовано для вычисления плотности энергии следующим образом

g = (Xq)/V, (14)

где V — объем окрестности рассматриваемого узла; q — энергия, которой соответствует полный вес N пакетов фотонов.

Плотность энергии в свою очередь может быть использована для вычисления потока энергии

V = Q/^a, (15)

где ^a — коэффициент поглощения в узле сетки.

Метод Монте-Карло является численным статистическим методом. Причем моделируемый процесс является нормальным, т.е. имеет Гауссово распределение. Для оценки точности воспользуемся правилом «3 сигм».

P{\Q -Х\ < 3(5) = 0,997. (16)

То есть погрешность метода почти наверняка равна 3S. Причем, величина S нам не известна заранее. Оценкой этой величины является выборочная дисперсия.

n

5'2 = (1/«)-^(Хг-Х)2. (17)

i= 1

Тогда можно считать, что

¿«л/52.

(18)

Такая оценка точности полностью учитывает погрешность самого метода.

В качестве тестовой была решена задача распространения лазерного излучения в коже. Расчетная область представляет собой 3 слоя: роговой слой, эпидермис и дерма. Для учета отражения от границы раздела сред: воздух/ роговой слой введен четвертый слой. В таблице приведены оптические параметры для каждого из слоев [5, 6].

Параметры расчетной среды

Название слоя А, длина волны, нм п. показатель преломления 1-1а, коэффициент поглощения, см-1 коэффициент рассеяния, см-1 9> анизотропия 4, размер слоя, см

Роговой слой 400 1,53 230 2000 0,9 0,0019

Эпидермис 1,40 66 800 0,74 0,01

Дерма 1,40 3,5 320 0,74 0,05

<— Дерма

<- Эпидермис <— Роговой слой <- Воздух

Падающее I—I излучение

Рис. 3. Геометрия среды.

О 0.5 1 1.5 2 2.5 3 3.5 4 X. СМ

Рис. 4. Плотность поглощенной энергии.

Дж/см3

О 0.5 1 1.5 г 2.5 3 3.5 4

X, СМ

Рис. 5. Абсолютная погрешность значения плотности поглощенной энергии.

На рисунке 3 приведена геометрия расчетной области. На основе полученных данных вычисляется математическое ожидание и дисперсия в узлах. Полученные результаты для 100 000 пакетов фотонов представлены на рисунках 4 и 5. Точность решения для такого количества пакетов фотонов составляет 10%.

Заключение. В ходе работы произведено построение модели распространения света в биологических тканях, основанное на методе Монте-Карло, созданы алгоритмы решения поставленной задачи. В качестве оценки точности метода использовалась выборочная дисперсия. Решен ряд тестовых задач с произвольной геометрией среды. На основе полученных результатов была проведена оценка эффективности метода. При запуске 100 000 пакетов фотонов погрешность составляет в среднем 10% от величины. Полученные результаты хорошо согласуются с исследованиями других авторов. Предлагаемый метод является весьма гибким и позволяет моделировать процессы для сред, имеющих сложную геометрию.

Литература

1. Seteikin A., Krasnikov I. 2007, Proc. of SPIE. Vol. 6826. P. 1404.

2. Seteikin A., Krasnikov I. Book of Abstracts International Conference Advanced Laser Technologies ALT-07, Levi, Finland. P. 107.

3. Meglinski Igor V. Physics in medicine and biology. 2002. N 47. С. 4271-4285.

4. Meglinski Igor V. Physiological measurement. 2002. N 23. С. 741-753.

5. Jacques S., Wang L. Monte Carlo Modeling of Light Transport in Tissue. In: Welch A. J., Martin J. C. Van Gemert. Optical-thermal response of laser-irradiated tissue Plenum press. New York, 1995. Vol. 12. С. 301-357.

6. Welch A. J. Optical-thermal response of laser-irradiated tissue // A.Welch, M. Vangemert. Plenum Publishing Corporation, 1995. 952 c.

7. Соболь И. М. Метод Монте-Карло. М.: Наука. Гл. ред. физ.-мат. лит-ры, 1985. 80 с.

Статья поступила в редакцию 15 июня 2010 г.