CC BY
110
Наука и Образование. МГТУ им. Н.Э. Баумана. ВсШИС Электрон, журн. 2014. № 12. С. 350-364."
Б01: 10.7463/1214.0738862
Представлена в редакцию: 15.11.2014
Сетевое научное издание
© МГТУ им. Н.Э. Баумана
УДК 535.2
Статистическое моделирование переноса излучения и переходные характеристики многослойной биоткани
Макаров С. Ю.1'* "тьи-роа^таЛл
волгоградский государственный университет , Волгоград, Россия
В статье описан метод Монте-Карло для моделирования распространения монохроматического излучения в многослойных биологических тканях с оптическими параметрами, близкими к коже и подкожным тканям. Статистический подход позволил вычислить стационарное распределение светового поля в такой биоткани. Процесс вычисления занимает примерно полчаса для современных ПК. При учете конечного значения скорости света в исходном алгоритме стало возможным получить зависящие от времени коэффициенты поглощения слоев биологической ткани. Эти зависящие от времени функции имеют вид переходных характеристик и могут помочь оценить возможность использования стационарной модели переноса излучения в многослойных биологических тканях при воздействии коротких световых импульсов.
Ключевые слова: Метод Монте-Карло, освещенность, коэффициент поглощения, функция Хеньи-Гринштейна
Введение
Методы статистического моделирования получили широкое распространение для решения различных научных и прикладных задач [1][2]. Особенно эффективно их применение в ситуациях, когда получение решения другими методами затруднительно, например, в силу сложного характера аналитических зависимостей в исходных уравнениях или граничных условиях. Хотя метод влечет за собой необходимость в большом числе испытаний (для набора статистики), а значит и расчетов на ЭВМ, постоянный рост производительности компьютеров делает этот недостаток все менее существенным при решении вопроса о применимости метода в конкретной ситуации. Весьма эффективным оказывается применение метода Монте-Карло для нахождения характеристик распределения светового излучения в мутных (рассеивающих) средах. К задачам такого рода приводит, например, применение лазеров в биомедицинской диагностике [3]. Более того, в задачах распространения излучения в рассеивающих и поглощающих средах метод Монте-Карло утвердился как своеобразный эталонный метод, по которому оценивают адекватность тех или иных приближенных методов решения уравнения переноса излучения в таких средах [4]. Это
Наука йОбразо
МГТУ им. Н.Э. Баумана
объясняется его универсальностью (метод работает при произвольных соотношениях оптических коэффициентов, любых кратностях рассеяния и конфигурациях граничных условий). В данной статье метод Монте-Карло используется для моделирования распространения излучения в многослойной биоткани. Достаточно подробно описывается технология такого моделирования. Далее алгоритм модифицируется для учета конечной скорости распространения фотона в слоях. На этой основе автором вводятся световые переходные характеристики для такой биоткани, дающие критерий применимости стационарных моделей переноса излучения при использовании коротких лазерных импульсов в процессах, использующих поглощение энергии лазерного излучения в глубине биоткани. Кроме того, результаты моделирования используются для получения распределения энергии излучения в многослойной биоткани от пучка с гауссовым профилем интенсивности.
1. Анализ проблемы
При распространении света через среду со случайными рассеивателями приходится учитывать два эффекта - поглощение и собственно рассеяние. Оба они приводят к тому, что интенсивность строго коллимированного пучка, попадающего в полубесконечную среду, будет уменьшаться по мере продвижения вглубь среды по закону (Бугера) I(z)=I(0) exp(-ft z), fit - коэффициент экстинкции (ослабления) для данной среды, ft=fa+fs , где fia -коэффициент поглощения (вероятность поглощения фотона на единицу длины пути), fs -коэффициент рассеяния (вероятность рассеяния фотона на единицу длины пути), соответственно вероятность "выживания" фотона W0 = fs / ft , а вероятность "гибели" (т.е. поглощения) соответственно fa\/ft • В общем случае требуется определить результирующее световое поле в заданных точках после рассеяния, поглощения и отражения от границ. Так как в рассеивающей среде присутствует некогерентная составляющая излучения (превалирующая над когерентной), то в точке r для заданного направления (луча) s нужно записывать среднюю плотность потока энергии в единичном телесном угле - лучевую интенсивность I [5, с. 165]:
dP
I(r, s) V/ /гл , Вт м-2 ср-1
da dQ
Введя интегральную интенсивность (освещённость)
U(r) = JI(r,s) dQ
4л
можно выразить поглощаемую единицей объёма среды мощность в виде
A(r) =juaU (r) (1)
Идея метода Монте-Карло для исследования стационарного распределения светового поля в рассеивающей и поглощающей среде состоит в нахождении статистическим образом распределения плотности вероятности поглощения фотона в объёме среды, а затем использовании формулы (1) для восстановления искомого распределения.
2. Модель многослойной биоткани
В качестве многослойной биоткани будем рассматривать кожу человека. С оптической точки зрения кожа представляет собой поглощающую среду с ярко выраженными рассеивающими свойствами. Основными компонентами, определяющими поглощающие свойства кожи в видимом и ближнем ИК спектральном диапазоне являются меланин, содержащийся в пигментированном эпидермисе, а также гемоглобин крови. Поглощение жировой ткани определяется поглощением липидов, воды и пигмента в-каротина. Рассеяние определяется, в основном, фиброзной структурой дермы кожи, хотя определенный вклад в светорассеяние вносит и эпидермис кожи, основными рассеивателями в котором являются митохондрии клеток. Для дермы реальное рассеяние света происходит на колла-геновых волокнах, образующих ее структуру, и узлах образованных сплетением отдельных волокон [6, с. 10].
Для исследования распределений света представим кожу как плоскую многослойную рассеивающую и поглощающую свет среду, на поверхность которой нормально падает лазерный пучок. При этом предполагается, что каждый слой характеризуется своим набором параметров: коэффициентами /ла и Цц, фактором анизотропии g, толщиной ё и показателем преломления п. Оптические характеристики слоев зависят от длины волны излучения. Для определенности используем модель из 6 слоёв, оптические свойства которых для длины волны 633 нм приведены в таблице 1 [7, с. 521] [8, с. 871].
Таблица 1. Оптические параметры модели кожи, Х=633 нм
№ Слой Ца см-1 Ц8 см-1 £ п Б мкм
1 Эпидермис 4,3 107 0,79 1,5 100
2 Папилярная дерма 2,7 187 0,82 1,4 200
3 Поверхностное сосудистое сплетение 3, 3 192 0,82 1,4 200
4 Ретикулярная дерма 2,7 187 0,82 1,4 900
5 Глубокое сосудистое сплетение 3,4 194 0,82 1,4 600
6 Гиподерма (жировая ткань) 1,1 130 0,82 1,4 2000
3. Моделирование методом Монте-Карло распространения излучения в
многослойных средах
Будем рассматривать бесконечно тонкий световой пучок, перпендикулярно падающий на многослойную среду. Кроме задания параметров каждого слоя (таблица 1), необходимо задать показатель преломления окружающей среды сверху и снизу всех слоёв для правильного учета отражения от границ. Статистический метод моделирования в данном случае предполагает многократный запуск фотона, представляемого в виде объекта-
пакета, которому приписывается исходный вес Ж=1, в биоткань [9][10]. Далее рассчитывается распространение фотонов в 3-х измерениях и записываются собранные вследствие поглощения части веса фотонного пакета в каждой ячейке пространственного массива А (х, у, г), после чего рассчитывается освещенность по формуле (1). Плотность распределения вероятности отражения и прохождения (см'2 ср'1) моделируется путем подсчета веса покинувших поверхности фотонов. Задача считается цилиндрически симметричной, так что пространственный массив А (х, у, г) представлен двумерной матрицей А (г, г), г - радиальная координата цилиндрической системы координат (ц. с. к), от точки падения луча, 2- глубина, отсчитываемая от поверхности 1-го слоя. Декартова система координат (д. с. к.) используется для трассировки распространения фотона, г координата ц. с. к. используется также для расчета плотности распределения диффузного отражения Я (г, а) и пропускания Т (г, а), где а - это угол между направлением «убегания» фотона от поверхности и внешней нормалью к поверхности (-2 для отражения и +2 для пропускания). Подвижная сферическая с. к., ось г которой динамически подстраивается к текущему направлению распространения фотона, используется для расчета изменения направления движения фотона из-за рассеяния. В этой с. к. сначала рассчитываются угол отклонения в и азимутальный угол ф, полученные из-за рассеяния, а затем направление распространения фотона обновляется в неподвижной декартовой с. к. с использованием направляющих косинусов.
Для реализации метода Монте-Карло прежде всего необходимо смоделировать случайную величину х с произвольным распределением вероятности р(х) на отрезке [а, Ь] при помощи случайной величины £ с равномерным распределением на отрезке [0, 1], генерируемой стандартными средствами. Для этого предполагается существование взаимнооднозначного отображения [0, 1] и [а, Ь], например, неубывающей функцией, так что интегральные функции распределения обеих случайных величин равны:
X
4 = ^ (£) = ^ (х) = | р(Х¥Х
а
Это равенство затем решается относительно х (аналитически или численно). Например, чтобы получить представление для случайной величины длины свободного пробега фотона в, можно использовать классическую форму интегральной функции вероятности для длины свободного пробега: Р{в<в'} =1- вхр(-^1 в') и, соответственно, для плотности распределения вероятностей: р(в) = ехр(-р.1 в). В итоге получим программное выражение для случайной длины свободного пробега фотона в виде
„ - 1п(1 -4) Иг .
Учитывая равномерный характер распределения на отрезке [0,1], это выражение эквивалентно следующему
- 1п4
^ = -
Иг
В частности, отсюда легко получить среднюю длину свободного пробега: <s> = <-ln / |it = 1/ |it , что проверяется непосредственно по распределению вероятностей. Свободный пробег фотона описывается приращениями:
х = х + juxs
У = У + Uys (2)
z = z + Uzs
Как только фотон совершает свой шаг, некоторая часть его веса уменьшается вследствие поглощения при взаимодействии со средой:
AW = W — Ut
Собранный в элементе массива A(r,z) вес будет увеличен на это же количество: A(r, z) =A(r, z) + AW, а вес самого фотона уменьшен: W =W - AW. После пробега и уменьшения веса вследствие поглощения фотон может испытать рассеяние. Угол отклонения ве [0, п], и азимутальный угол у е[0, 2 п) моделируются также статистически. А именно, используется распределение вероятности для cos в:
1-g2
pw (cos())) = -ту,-2—;-, (Henyey and Greenstein, 1941).
2 (1 + g - 2g cos в)
Таким образом
+1
< cos в > = Jpw(cos в) d cos в = g.
Типичное значение для биотканей g ~ 0.8 в видимом спектре. Соответственно, получаем статистические представления:
cos0 = — 2 g
1 + g2 -
1 - g2
1 - g + 2 g£
, g > 0
cos0 = 1, g = 0
Представление для равномерно распределенного азимутального угла рассеяния проще: у=2п£. Сразу же после выбора cosQ и у направляющие косинусы рассеянного фотона принимают новые значения:
sin#
ßx = I _ (ßxßz coW-ßy sinv)+ßx cos°
V1 -ßz
ßv = ~F=== (ßeßz coW + ßx sinv)+ßy
л/1-ßz
2
ßz =- sin^cos^^/1 - ßz 2 + ßz cosd
Однако, если исходное направление фотона очень близко к нормали поверхности ткани, (например, /л2 >1-10-6) то используются менее "погрешные" формулы: /л'х = sind
-1
<
cosy, ¡л'у = sinQ siny, ¡¡'z = SIGN(¡z) cosQ. Этап моделирования рассеяния завершается обновлением направляющих косинусов: ¡¡х = ¡¡'х, ¡¡у = ¡¡'у, ¡¡z = ¡¡'z.
В процессе движения согласно полученной длине шага фотон может пересечь границу раздела «ткань - окружающая среда». Тогда для него есть две возможности:
1) пройти границу и наблюдаться как часть диффузного отражения или прохождения (если речь идёт о нижней границе);
2) испытать внутреннее отражение и вернуться в ткань.
Для решения возникающей проблемы рассчитывается так называемая сокращенная длина свободного пробега: s1 = (z0 - z) /¡¡z при ¡¡z <0, s1 = (z1 - z) / ¡¡z - при ¡¡z>0, где z0 и z1 это z-координаты верхней и нижней границ текущего слоя, s1 представляет собой расстояние от текущего положения фотона до границы по направлению движения фотона. Если фотон испытает полное внутреннее отражение, то он пролетает остаток пути s=s-s1. Если нет, то рассчитывается вероятность для данного фотона испытать внутреннее отражение, которая зависит от угла падения на границу ai = arccos \¡z\. Закон Снеллиуса даёт связь между углами падения при прохождении ni sin ai = nt sin at , а сам коэффициент отражения рассчитывается с использованием формул Френеля
1
) = 2
2 2
sin (a¡ -at) tan (ai -at)
2 2
sin (ai +at) tan (a¡ +at)
как среднее для двух ортогональных направлений поляризации. Далее определяется, будет ли данный фотон испытывать внутреннее отражение, путем генерации случайного числа и сравнения его значения с рассчитанным коэффициентом отражения: если £ < R(а.i).¡ то фотон испытает внутреннее отражение, если £ > R(аi), то фотон покинет ткань. Если реализовалось внутреннее отражение, то фотон остаётся на поверхности, до которой он долетел по сокращенной длине, а его направляющие косинусы обновляются: (^х, /лу, /л2) = (^х, ИУ, ~Иг). В этой точке проверяется остаток длины свободного пробега Если она достаточно большая, чтобы опять достигнуть границы раздела «ткань-окружающая среда», то описанный выше процесс повторяется. Если остаток длины умещается в пределах одного слоя, то фотон продвигается на этот остаток длины, и затем происходят описанные выше этапы поглощения и рассеяния. В случае, если остаток длины «упирается» в границу между слоями, то эта ситуация обрабатывается подобным образом, а именно: сначала фотон продвигается по укороченному пути в1 до границы слоев, с параметрами, например, слой 1: /иа1, цв1, п1, слой 2: /ла2, /л&, п2. Далее обновляется шаг фотона: в = в - в1. Затем по Френелю определяется, будет ли этот фотон отражен, или пройдет во второй слой. Первый вариант рассматривается аналогично границе с внешней средой. В случае же прохождения, шаг должен быть скорректирован так, чтобы оптическая длина пути при переходе границы раздела не изменилась: в = в • ци/ря. Описанный процесс повторяется, пока на очередном этапе путь фотона не уместиться в пределах одного слоя, после чего происходят поглощение и рассеяние в обычном порядке. Если же реализовалось покидание фотоном ткани, то инкрементируется накопленное отражение или прохождение в элементе [г,а] соответствующих массивов. Поскольку такой фотон покинул ткань, его трассировка пре-
кращается и в биоткань запускается новый фотон. Если же фотон подолгу остаётся в ткани, не вылетая наружу, то его вес, очевидно, стремиться к 0, и когда достигнет порога (например, W=0,000001), то его трассировку необходимо оборвать, т.к. такой фотон вносит очень малый вклад в рассчитываемые величины, а затрачиваемое на него машинное время может быть значительным. Однако просто прервать процесс нельзя, т. к. число моделируемых фотонов как правило, весьма значительно (~ 10) и отбрасывание может привести к нарушению баланса энергии. Для корректного прерывания трассировки пороговых фотонов в программах моделирования по методу Монте-Карло обычно используется приём под названием «русская рулетка»: генерируется случайное число и подпороговому фотону даётся один шанс из m (m=10, например), выйти из рулетки с надпороговым весом (mW): W = mW при £ < 1/m, W = 0 при £>1/m. Такой способ является корректным с точки зрения сохранения энергии при больших числах моделируемых фотонов.
Собранный вес от поглощения фотонов хранится в массиве Arz[r, z], где r и z представляют сетку в r- и z- направлениях с шагом Ar и Az соответственно. Данные массива представляют собой полный собранный вес фотонов в каждом элементе двумерной сеточной системы. Чтобы получить распределение по какому-либо одному измерению, нужно просуммировать (свернуть) двумерный массив по другому измерению:
Nr-1
Az[z] = £ Arz[r, z]
r=0
Полный фотонный вес, собранный в слое (layer) A1 и полный фотонный вес, поглощенный всей тканью A затем получаются из массива Az[z]:
Al[layer] = £ Az[z]
z in layer
Nz -1
A = £ Az[z]
z=0
После соответствующего масштабирования получаются искомые величины для вероятностей. Для плотностей вероятности поглощения:
Arz[r,z] -з
Arzir, zl =-, см
1 ' J ASAzN '
Лф1 -i
Azlzi =-, см
L J AzN '
где N - общее число трассированных фотонов, ASAz - элемент объёма в ц. с. к.
Соответственно для вероятностей поглощения:
Al[layer] = Al[layer], A = A
L 7 J N N
Для получения разрешенных по времени вероятностей поглощения фотона, стандартный программный код метода Монте-Карло был модифицирован следующим образом: введена сетка по времени c шагом At = 1-10-12 с, константа скорости света c = 3-1010
см/с, а для фотонного пакета добавлен счетчик времени, инициализируемый с очередным запускаемым в биоткань фотоном. Соответственно, к формулам (2) добавляется приращение счетчика времени: I = I + в п / с, где п - показатель преломления того слоя, в котором находится фотон. Далее, для каждого слоя вводится массив PAt[] для отслеживания во времени поглощения света в этом слое и получения разрешенной по времени вероятности поглощения. Для найденного показания счетчика времени фотона массив инкрементиру-ется: PAt[t] = PAt[t] + А Ж (соответствующее этому t поглощение в слое). Введём световую переходную характеристику как функцию времени, получаемую в результате суммирования найденной путём моделирования разрешенной по времени вероятности поглощения фотона в слое (после нормирования на полное число фотонов):
Р[ь] =±РАфк]
к=0
На рис. 1 представлен вид рассчитанных таким образом характеристик для некоторых слоёв применяемой модели кожи. Видно, что для глубоких слоёв дополнительно наблюдается некоторая задержка в отклике, обусловленная, очевидно, необходимым временем пролета фотонов от поверхности биоткани до этих слоев. Из рисунка также видно, что стационарная модель переноса излучения для данной структуры биоткани применима, если временной параметр импульса, определяющий рассматриваемый эффект воздействия лазерного излучения на биоткань при поглощении, составляет не менее 60 пс.
Итоговые результаты представлены в таблице 2. Число трассируемых фотонов в вычислительном эксперименте равнялось 10 , шаг сетки 0,01 мм. Затраченное время - около 30 мин при расчете на процессоре типа 17.
ретикулярная дерма
0,2
0,15
0,1
0,05
40 60
пикосекунды
100
Рис. 1. Световые переходные характеристики подкожной ткани.
0
Таблица 2. Результаты моделирования
№ Слой л1 стационарное поглощение, отн. ед. Ь время установления, пс
1 Эпидермис 0,1382 ~30
2 Папилярная дерма 0,1482 ~40
3 Поверхностное сосудистое сплетение 0,1402 ~40
4 Ретикулярная дерма 0,2267 ~60
5 Глубокое сосудистое сплетение 0,0478 ~60
6 Гиподерма (жировая ткань) 0,0128 ~60
Рассчитанный коэффициент френелевского отражения составил 0,04, диффузного -около 0,25. Наибольшее значение для данной структуры биоткани имел коэффициент поглощения (биоткани в целом) - более 0,71.
В заключении необходимо отметить, что, хотя результаты, полученные методом Монте - Карло описанным выше образом, относятся к коллимированному пучку бесконечно малого диаметра, применение аппарата функций Грина позволяет получить решение для пучков конечных размеров с произвольными профилями интенсивности, без не-
обходимости моделировать эти случаи отдельно. А именно, решение получается путем свертки найденного решения для бесконечно тонкого пучка и функции профиля интенсивности пучка конечных размеров:
+œ +œ
C(х,y,z) = J JG(x', y', z) S(x-x',y-y') dx'dy
—œ —œ
где G(x,y,z) - функция Грина для среды, S(x,y) - профиль интенсивности пучка конечных размеров. В случае цилиндрической симметрии получим
œ 2л ~ ~
C(r,z) = Jdr'G(r', z) r' J S(r2 + r'2 - 2rr'cos)) d)
0 0
Для пучка с гауссовым профилем интенсивности с эффективным радиусом пучка R
S (r) = S exp(—2(r / R)2)
œ 2л
C(r, z) = S (r) J dr'G(r', z) exp (-2 (r'/R)2) r ' J exp(4rr'cœ)/R2) d) 0 0
Второй интеграл выражается аналитически через модифицированную функцию Бесселя
1 2л
I0 ( x) = — J exp( x cos))d) 2л 0
и в итоге
œ
C (r, z) = S (r) J dr'G(r', z)exp(- 2 (r'/R)2 ) I0(4rr'/R2 ) 2лг'
0
0
Таким образом, остаётся только один интеграл, который можно рассчитать численными методами, т.к. цифровое представление функции Грина уже получено как решение для бесконечно узкого пучка.
На рис. 2 представлена рассчитанная таким образом освещенность в глубине используемой модели многослойной биоткани для гауссова пучка с нормированной мощностью 1 Вт с эффективным радиусом Я=1 мм. Из рис. 2 следует, что при заданных параметрах модели освещенность представляет собой непрерывную, достаточно гладкую функцию координат г и г. Это обусловлено равенством показателей преломления соседних слоев (кроме эпидермиса) и гауссовым характером распределения интенсивности в падающем пучке. Поглощенная мощность уже не будет непрерывной функцией из-за резкого различия в коэффициентах поглощения слоёв модельной биоткани.
Рис. 2. Распределение освещенности в модельной биоткани от гауссова пучка
Выводы
В силу статистического характера процессов поглощения и рассеяния света метод Монте-Карло является удобным и универсальным средством математического моделирования в задачах распространения излучения в биологических тканях. Универсальность метода позволяет использовать его в качестве эталонного при тестировании приближенных методов решения уравнения переноса излучения. Кроме того, неуклонный рост производительности ЭВМ позволяет все чаще использовать метод для решения прикладных задач биомедицинской диагностики. При этом дополнительным плюсом метода является его применимость для областей произвольной формы.
В представленной статье метод Монте-Карло использовался для моделирования распространения монохроматического излучения в многослойной биоткани, со слоями, близкими по оптическим характеристикам к коже и подкожным тканям. По результатам моделирования можно сделать следующие выводы.
(1) Метод Монте-Карло позволил решить задачу нахождения распределения излучения в многослойной рассеивающей и поглощающей среде на современных процессорах за приемлемое время. В результате удалось рассчитать коэффициенты френелевского и диффузного отражения от слоистой биоткани, а также коэффициент пропускания и поглощения такой средой.
(2) Модификация стандартного алгоритма метода для учета времени свободного пробега фотонов позволила рассчитать разрешенные по времени вероятности поглощения фотона в слоях - световые переходные характеристики многослойной биоткани. Эти характеристики могут служить критерием применимости стационарной модели переноса излучения в многослойной биоткани.
(3) Использование математического аппарата функций Грина позволяет быстро получить решение для лазерного пучка с произвольным профилем интенсивности на основе однократно проведенного моделирования методом Монте-Карло для пучка бесконечно малой толщины.
(4) Для пучка с гауссовым профилем интенсивности освещенность внутри многослойной биоткани оказывается непрерывной функцией координат (при одинаковых коэффициентах преломления соседних слоев). Соответственно, поглощенная мощность будет разрывной функцией на границах раздела слоев с резко различающимися коэффициентами поглощения.
Список литературы
1. Соболь И М. Метод Монте-Карло. М.: Наука, 1978. 64 с.
2. Хан Г., Шапиро С. Статистические модели в инженерных задачах: пер. с англ. М.: Мир, 1969. 395 с.
3. Быков А.В., Кириллин М.Ю., Приезжев А.В. Моделирование методом Монте-Карло сигнала оптического когерентного доплеровского томографа: влияние концентрации частиц в потоке на восстановленный профиль скоростей // Квантовая электроника. 2005. Т. 35, № 2. С. 135-139.
4. Братченко И. А., Захаров В. П. Приближенный метод расчета распределения энергии оптического излучения в многократно рассеивающих средах // Компьютерная оптика. 2008. Т. 32, № 4. С. 370-374.
5. Исимару А. Распространение и рассеяние волн в случайно неоднородных средах. В 2 т. Т. 1. Однократное рассеяние и теория переноса: пер. с англ. М.: Мир, 1981. 280 с.
6. Башкатов А.Н., Генина Э.А., Тучин В.В. Исследование оптических и диффузионных явлений в биотканях при воздействии осмотически активных иммерсионных жидкостей: Общий биофизический практикум. Саратов: Саратовский гос. ун-т, 2005. 71 с.
7. Тучин В.В. Исследование биотканей методом светорассеяния // Успехи физических наук. 1997. Т. 167, № 5. С. 518-525.
8. Башкатов А. Н., Генина Э. А., Кочубей В. И. Тучин В. В. Оптические свойства подкожной жировой ткани в спектральном диапазоне 400-2500 нм // Оптика и спектроскопия. 2005. Т. 99, № 5. С. 868-874.
9. Prahl S.A. Light transport in tissue: PhD dissertation. University of Texas at Austin, 1988.
10. Wang L.-H., Jacques S.L. Monte Carlo Modeling of Light Transport in Multi-layered Tissues in Standard C. University of Texas, M. D. Anderson Cancer Center, 1992. 173 p.
Science and Education of the Bauman MSTU, 2014, no. 12, pp. 350-364.
DOI: 10.7463/1214.0738862
Received:
15.11.2014
Science ^Education
of the Bauman MSTU
ISSN 1994-0448 © Bauman Moscow State Technical Unversity
Statistical Modeling of Radiative Transfer and Transient Characteristics for Multilayer Biological Tissue
S.Yu. Makarov
1,*
msu-po£t@mailju Volgograd State University, Volgograd, Russia
Keywords: Monte Carlo method, illumination, absorption coefficient, Henyey-Greenstein function
The Monte-Carlo method [1] already long ago proved itself as a powerful and universal tool for mathematical modelling in various areas of science and engineering. Researchers often choose this method when it is difficult to find a solution by other ways (or impossible at all), e.g. because of sophisticated analytical dependences, area of modelling or boundary conditions. Certainly, this necessarily statistical and flexible method requires significant computation time, but a continuously increasing computation capability makes it more and more attractive for a choice in specific situation.
One of the promising areas to use the method of statistical modelling is description of light propagation in the turbid (scattering) media. A high motivation for development of this approach is widely used lasers in biomedicine [3]. Besides, owing to its flexibility, the Monte-Carlo method is also of importance in theoretical researches, in particular, to estimate a degree of adequacy of the offered approximation methods for solving a radiative transfer equation [4].
It is known that key parameters of turbid media are an absorption coefficient (characterizes absorption probability of a photon per unit of path length) and a scattering coefficient (characterizes scattering probability of a photon per unit of path length). The ratio of each of the coefficients to their sum (extinction) defines a probability of "death" or "survival" of a photon, respectively, in interaction with lenses. Generally, in the scattering medium there is a non-coherent radiation component, which in turbid media such as biological tissues, already at the insignificant depth becomes prevailing over the coherent one (residual of the incident laser beam) [5].
The author used the Monte-Carlo method to simulate optical radiation propagation in the multilayer biological tissues with their optical characteristics corresponding to the skin and subcutaneous tissues. Such a biological tissue is the absorbing and scattering medium. The main absorbers in visible and near-IR spectral range in this case are melanin (contains in pigmented epidermis) and blood hemoglobin. Scattering is defined both by the fibrous structure of skin, and by the cell organelles. The optical characteristics, absorbing and scattering properties of skin, and subcutaneous tissues are rather well studied [6][7][8].
The aim of modelling was to have both the stationary distribution of light radiation in multilayer biological tissue and the transitional characteristics corresponding to such biological tissue structure. These characteristics may appear when a biological tissue is subjected to short laser pulses. Discovering of such characteristics and their analysis allows us to formulate a criterion for the validity of stationary models of radiation transfer in multilayer biological tissue. As to transitional characteristics (i.e. time-allowed absorption coefficients in layers), the standard algorithm of the Monte-Carlo method [9][10] was modified to take into consideration a final speed value of light distribution in layers. In such a way, in particular, it was found that time about 60 Ps (for radiation with the wavelength of 633 nm) is necessary to obtain the steady-state (stationary) values of the absorbed power in skin and subcutaneous tissues. The entire modelling process (based on the i7 processor) took about 30 minutes, thus the number of the traceable photons made 100 million.
The results obtained for a beam of infinitesimal diameter allowed us to find the solution on distribution of light energy in biological tissue for a laser beam of a final width with a Gaussian profile of intensity. For this purpose, was used the mathematical apparatus of Green function. Its digital representation in this case is a modelling result for a beam of infinitesimal diameter. Since indexes of refraction of the adjacent layers in model were specified to be equal, the received function of integrated intensity for a Gaussian beam in the depth of biological tissue became a continuous and smooth function. In turn, it means that the function of absorbed power for the considered model of biological tissue will not be any more a continuous function, as the absorption coefficients of the adjacent layers are different.
Thus, the Monte-Carlo method applied to simulate a distribution of the optical radiation in multilayer biological tissue allowed us to have a required distribution of the resultant light field in a stationary case for the acceptable time, as well as to obtain time-allowed absorption coefficients in layers (transitional characteristics of light), which can be used to estimate the validity of stationary models of radiation transfer when biological tissue is subjected to short light pulses.
References
1. Sobol' I.M. MetodMonte-Karlo [Monte Carlo Method]. Moscow, Nauka Publ., 1978. 64 p. (in Russian).
2. Hahn G.J., Shapiro S.S. Statistical models in engineering. New York, Wiley 1967. (Russ. ed.: Hahn G.J., Shapiro S.S. Statisticheskie modeli v inzhenernykh zadachakh. Moscow, Mir Publ., 1969. 395 p. ).
3. Bykov A.V., Kirillin M.Yu., Priezzhev A.V. Monte Carlo simulation of an optical coherence Doppler tomograph signal: the effect of the concentration of particles in a flow on the reconstructed velocity profile. Kvantovaya elektronika, 2005, vol. 35, no. 2, pp. 135-139. (English translation: Quantum Electronics, 2005, vol. 35, no. 2, pp. 135-139. DOI: 10.1070/QE2005v035n02ABEH003337 ).
4. Bratchenko I. A., Zakharov V. P. Approximate method of optical energy distribution calculation in multiple scattering mediums. Komp'yuternaya optika = Computer Optics, 2008, vol. 32, no. 4, pp. 370-374. (in Russian).
5. Ishimaru A. Wave Propagation and Scattering in Random Media. Vol. 1. Single Scattering and Transport Theory. N.Y. Academic Press, 1978. (Russ. ed.: Ishimaru A. Rasprostranenie i rasseyanie voln v sluchayno neodnorodnykh sredakh. V 2 t. T. 1. Odnokratnoe rasseyanie i teoriyaperenosa. Moscow, Mir Publ., 1981. 280 p.).
6. Bashkatov A.N., Genina E.A., Tuchin V.V. Issledovanie opticheskikh i diffuzionnykh yavleniy v biotkanyakh pri vozdeystvii osmoticheski aktivnykh immersionnykh zhidkostey [Investigation of optical and diffusion phenomena in biological tissues under the influence of osmotically active immersion liquids]. Saratov, SSU Publ., 2005. 71 p. (in Russian).
7. Tuchin V.V. Light scattering study of tissues. Uspekhi fizicheskikh nauk, 1997, vol. 167, no. 5, pp. 518-525. (English translation: Physics-Uspekhi, 1997, vol. 40, no. 5, pp. 495-515. DOI: 10.1070/PU1997v040n05ABEH000236 ).
8. Bashkatov A. N., Genina E. A., Kochubey V. I. Tuchin V. V. Optical properties of the subcutaneous adipose tissue in the spectral range 400-2500 nm. Optika i spektroskopiya, 2005, vol. 99, no. 5, pp. 868-874. (English translation: Optics and Spectroscopy, 2005, vol. 99, no. 5, pp. 836-842. DOI: 10.1134/1.2135863 ).
9. Prahl S.A. Light transport in tissue: PhD dissertation. University of Texas at Austin, 1988. 10. Wang L.-H., Jacques S.L. Monte Carlo Modeling of Light Transport in Multi-layered Tissues
in Standard C. University of Texas, M. D. Anderson Cancer Center, 1992. 173 p.
