CC BY
54
Наука и Образование
МГТУ им. Н.Э. Баумана
Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 12. С. 87-109.
Б01: 10.7463/1215.0826701
Представлена в редакцию: Исправлена:
© МГТУ им. Н.Э. Баумана
25.10.2015 10.11.2015
УДК 004.9+535.2
Моделирование распространения ненаправленного излучения и неинвазивное определение параметров рассеяния и поглощения излучения в биоткани
Макаров С. Ю.1'* 'т&ц-роы:@таЦ:ш
волгоградский государственный университет,
Волгоград, Россия
В статье исследуется распространение ненаправленного излучения от ламбертова пучка света с естественной поляризацией, освещающего поверхность биоткани. Для моделирования методом Монте-Карло такого излучения используется найденная функция статистического представления для угла падения фотонов, верифицированная аналитическим выражением для коэффициента прохождения мощности. Описан метод фиксации мощности, передаваемой ненаправленным излучением в биоткань, и применение диффузионного приближения для этого случая. Выявлены особенности поведения решения на поверхности биоткани внутри и вне падающего пучка, которые предложено использовать для неинвазивного определения параметров диффузионного приближения расчетным методом. Также предложен принцип измерения подповерхностной освещенности, необходимый для такого расчета, и описан принцип измерения полного коэффициента отражения, знание которого при известных (или измеренных) параметрах диффузионного приближения позволяет определить дополнительно ещё два параметра, характеризующих биоткань как рассеивающую среду - коэффициент рассеяния и параметр анизотропии.
Ключевые слова: перенос излучения, диффузионное приближение, неинвазивные методы
Введение
Целью исследования в статье является моделирование распространения рассеянного (ненаправленного) излучения от ламбертова пучка света с естественной поляризацией, освещающего поверхность биоткани, а также возможность использования особенностей распространения такого излучения для неинвазивного определения параметров распространения и поглощения монохроматического излучения внутри биоткани. Такими параметрами являются параметры диффузионного приближения [1], а также коэффициент рассеяния и параметр анизотропии. Соответственно, постановка задачи может быть сформулирована следующим образом: 1) создание математической и численной модели переноса излучения в диффузионном приближении для пучка рассеянного света, падающего на поверхность биоткани, 2) создание на основе построенной модели метода неинвазивного оп-
ределения параметров рассеяния и поглощения излучения внутри биоткани. При этом для верификации модели и метода используется статистическое моделирование (метод Монте-Карло).
Для моделирования методом Монте-Карло [2] ненаправленного излучения используется полученная автором функция статистического представления угла падения дискретных фотонов [3] в таком излучении. Для удобства применения диффузионного приближения, в том числе для целей измерения параметров, используется принцип фиксации мощности, передаваемой в биоткань ненаправленным излучением [3], при этом статистическое моделирование применяется параллельно с диффузионным приближением на каждом шаге с целью верификации результатов.
Для однозначного восстановления двух неизвестных параметров диффузионного приближения предлагается использовать результаты измерений двух величин: значение выходящей из биоткани вследствие внутреннего рассеяния мощности излучения и значение подповерхностной освещенности, принцип измерения которой, совместимый с граничными условиями, описан в статье. Также изложен принцип измерения полного коэффициента отражения, знание которого позволяет восстановить значения ещё двух параметров, характеризующих биоткань как рассеивающую среду - коэффициента рассеяния и параметра анизотропии.
По отношению к результатам исследований других авторов в данной области предлагаемый метод может быть отнесен к классу непрямых двух- или трех-параметрических (если дополнительно восстанавливается коэффициент рассеяния) CW (т.е. стационарных) методов [4][5], но выгодно отличается от методов с пространственно-разрешенным диффузным отражением своей простотой (измеряется только интегральное диффузное отражение, без пространственной зависимости). Использование при моделировании рассеянного излучения найденной функции статистического представления угла падения фотонов, предложенный принцип измерения подповерхностной освещенности, выявленные особенности поведения решения на поверхности вне падающего пучка в условиях фиксации передаваемой от излучения мощности, а также использование этих особенностей для неинвазивного нахождения параметров переноса излучения в биоткани, составляют научную новизну публикуемых результатов.
Следует отметить, что в отличие от случая измерения оптических характеристик образцов биологических тканей in vitro, для которого разработаны наиболее широко используемые и точные методы на основе прямых измерений ослабления в тонких срезах и использования интегрирующих сфер [6], определение оптических параметров для живых тканей in vivo, как правило, сопряжено с более значительными аппаратными затратами и сильнее зависит от используемой математической модели распространения света в рассеивающей среде. Тем не менее, развитию таких неинвазивных методов уделяется большое внимание, стимулируемое применением лазеров в медицинской диагностике и терапии. Так, помимо использования оптоволоконных датчиков для доставки и приема излучения [7], было улучшено математическое обеспечение для получения значений оптических параметров из результатов измерений, в частности, с использованием ускоренных вариантов метода Монте-Карло [8] и алгоритмов нейронных сетей [9]. Интересным явля-
ется также использование возможностей конфокальной микроскопии для одновременного измерения коэффициента рассеяния и параметра анизотропии [10]. Однако для разработанных методов "точечный" характер измерений, даваемый, например, волоконно-оптическими датчиками, может, что называется, выйти боком при работе на реальных живых тканях из-за возможного наличия внутренних или поверхностных неоднородностей, что потребует проведения нескольких процедур измерения с последующим усреднением результата. В этих условиях предлагаемый автором метод, использующий интегральные по поверхности измерения, может иметь определенные преимущества.
Знание параметров распространения и поглощения монохроматического излучения для конкретной биоткани позволяет, во-первых, определить распределение излучения в объеме ткани, и во-вторых, найти мощность возникающих тепловых источников от поглощения энергии лазерного излучения, следовательно, дает возможность предсказать локальное повышение температуры в различных точках биоткани при заданной мощности источника излучения. С учетом уже имеющихся в теплокровных организмах внутренних (биологических) источников тепла [11], поддерживающих внутреннюю температуру, действие возникающих тепловых источников от лазерного нагрева может привести к термическому повреждению живой ткани. Поскольку разница между стационарной внутренней температурой (~37°С -38°С) и нижней границей термического повреждения клеточных мембран (~45°С) составляет всего несколько градусов, то возможность неинвазивного определения вышеупомянутых параметров для живых тканей приобретает особую актуальность.
Везде далее под биотканью может подразумеваться и произвольная мутная, т.е. рассеивающая и поглощающая (в оптическом или ближнем ИК диапазоне) среда, с оптическими параметрами, близкими к таковым для непрозрачных мягких биотканей человеческого организма, в т. ч. и фантомы соответствующих биотканей [12, с. 273].
1. Применение метода Монте-Карло для моделирования распространения ненаправленного излучения
Как известно, падающая на единицу поверхности мощность (освещенность поверхности) при освещении ламбертовым пучком, т.е. ненаправленным излучением с независящей от угла падения лучевой интенсивностью, будет равна, как и в законе Ламберта для излучения, величине
ж/2
^ = I2ж |соэв- Этвйв = пI,
0
где I - лучевая интенсивность падающего излучения. При вычислении прошедшей мощности нужно учитывать наличие критического угла (считаем первой средой воздух), кроме того, нужно, используя отношение векторов Пойтинга и дифференциалов телесных углов для обеих сред, выразить лучевую интенсивность прошедшего излучения через интенсивность падающего:
I
t
«3С0%| ^ i.
cose
Для естественного света (хаотическая поляризация), как известно
/ х2 /■ ч 2
2сов#
+
И2 = 1 [W+(tJ ]=1
IcosOj cos^j + Л cos^j
У
п соб^ + соб^
Используя эти формулы, автором было получено аналитическое выражение коэффициента прохождения мощности Тр = р / Р^ , имеющее достаточно громоздкий вид (при п>1):
^ . » 1, 1, 1 (п - 1)(3п2 +1) п2
Тр(п) = — (п + 1)--(п — 1) ,----—---- ,----- -
р 4 2 12 (п, 1)2 (п2 , 1)2
' (п + 1)2(п — 1)3 (п2 — 1)2 (п2 — 1)(3п4 + 2п2 + 3) (п + 1) + --^----Ь ----— 1п [п (п + 1)] — --—-:-- +
(л2 + i)2
2
2л
л3 — л 2 + 3л + 1 2л (л 2 — 1)
(л2 + 1)
(л2 + 1)
4л (л + 1)
+ 2
(л 2 — 1)2
1Л
1 +
(Л2 — 1) л (Л2 + 1)
(Л — 1)2
(Л2 + 1)
1
С
>
расчет по которому дает для л=1,4 значение Тр = 0,9232 (соответственно для коэффициента отражения 0,0768).
При моделировании методом Монте-Карло, при условии, что модельные фотоны падают на освещаемую поверхность под произвольными углами в из диапазона [0, п/2] с равной вероятностью, для л=1,4 был получен коэффициент прохождения по мощности, равный (с точностью до 4-го знака) 0,8693 (и соответственно коэффициент отражения 0,1307). Расхождение с результатом аналитического расчета вызвано тем, что при моделировании методом Монте-Карло ненаправленного излучения вместо равномерного распределения по углам падения для дискретных фотонов необходимо использовать функцию статистического представления в виде [3]:
в = arcsm(V^) (1)
где £ — случайное число из диапазона [0, 1]. Формул а (1) вытекает из требования соблюдения условия Ii = const для представления ламбертова пучка в модели дискретных фотонов с учетом выражения для дифференциала телесного угла. Статистическое моделирования с учетом (1) показало полное согласие с аналитическим расчетом для величины Тр . Таким образом, при моделировании распространения ненаправленного излучения методом Монте-Карло, для получения правильных значений энергетических величин в рассеивающей и поглощающей среде следует использовать статистическое представление (1), что в дальнейшем будет подразумеваться по умолчанию.
2. Сведения из теории переноса излучения
Везде далее подразумевается, что временные параметры воздействия излучения на биоткань таковы, что допустимо использовать стационарные модели переноса излучения. Для кожи и подкожных тканей это справедливо, например, при наносекундных и более длительных импульсах [2]. В этих условиях можно использовать стационарное уравнение теории переноса излучения [13, с. 49]:
МрЙ = -^/(г,8) + / IРаО,^Дг^')(2)
оя Ах
где / - коэффициент экстинкции (ослабления), / - коэффициент рассеяния, арп(8,8') задает распределение вероятности рассеяния по телесным углам с направления 8' в направление 8. Это распределение (в литературе иногда называемое фазовой функцией), определяется дифференциальным сечением рассеяния:
p0(s,s') , a, = \addd
s - J ^ d1 as 4л
Для осесимметричного рассеяния фазовая функция определяется только углом (или косинусом угла) между направлением падающего и рассеянного фотона, а так как dd = 2 л sin ede, то можно перейти к распределению по косинусу угла рассеяния или к
индикатрисе рассеяния %(&) = 4л pd (в), при этом
л 1 л
J pn dd. = J pn (в)2л sin в de = - J х(в) sin в de =
4л 0 2 0
- +1 +1
— J^(cose) dcose = Jpcos j(cose)dcose = 1
2 -1 -1
Соответственно, средний косинус угла рассеяния (параметр анизотропии), вычисляется стандартным образом:
л 1 л
g=< cose >= J cose pn (в)2лsinв de = 1J cose^(e)sine de = 0 2 0 1 +1 +1 (3) — J cose^(cose) d cose = J cose pcos jcose) d cose.
2 -1 -1
Разлагая индикатрису рассеяния в ряд по полиномам Лежандра
2j +1 1
г ""j о J ' j=0 2 -1
и, принимая во внимание, что P0(x) = 1, P1(x)=x, найдем с учетом (3)
z(e) -Z0 +х ^(cose) = 1+з g cose.
Таким образом, параметр анизотропии g=<cos$> определяет "вытянутость" индикатрисы рассеяния по сравнению с изотропным случаем, т.е. большие значения g означают преимущественное рассеяние на малые углы. Значение g ~ 1 для реальных индикатрис означает, что рассеяние происходит в основном вперед, g ~ -1 - в основном назад, а g = 0
z(cose) =Y.XjPj (cose), Xj = Jx(cose) Pj (cose)d (cose),
соответствует изотропному рассеянию. Большие (близкие к 1) значения фактора анизотропии позволяют падающему излучению глубже проникать в биоткань, т.к. при этом увеличивается эффективная длина пробега фотона, при котором он "забывает" свое первоначальное направление движения, равная 1/(^(1-g)). Так, например, для дермы кожи
значение g=0,8 для желтого света увеличивает среднюю длину пробега в 5 раз, в среднем с 50 мкм до 250 мкм.
Поскольку аналитическое решение для уравнения переноса (2) не может быть получено в силу его незамкнутости, то обычно его преобразуют так, чтобы преобразованные уравнения включали в качестве искомых интегральные, измеримые на опыте величины, например, полную (интегральную) освещенность и плотность потока энергии
U(r) = JI(r,s) dQ,
4л
F(r) = JI(r,s) s dQ.
4л
При этом, если падающее на биоткань излучение представляет собой упорядоченное (например, однонаправленное) излучение, как в коллимированном лазерном пучке, то целесообразно в решении выделить компоненту U , оставшуюся от падающего излучения, прошедшего границу раздела и толщу среды, ослабленную рассеянием и поглощением [1]: U(r) = U(r) + Ц) (r), т.к. эта компонента сразу находится интегрированием по поверхности, на которую падает излучение:
UT(r) =ff1 (r0,S) cosв0da0, r0 еГ, z = J^ds, s = ^^, cose0 = (n,—s) Г | r — r0 |2 r0 | r — r01
Остающаяся компонента U^(r) представляет собой диффузное излучение, обусловленное рассеянием в среде.
3. Моделирование в диффузионном приближении для пучка рассеянного
излучения, освещающего биоткань
Диффузионное приближение для нахождения U^r) выводится в предположении, что
изотропная часть в интенсивности существенно больше неизотропной, порождающей поток энергии:
1 3
ID (r,s) = — UD (r) + — Fd (r) s, UD (r) >> \Fd (r)| (4)
4л 4л
В силу того, что диффузная компонента генерируется только внутри среды из падающего излучения при рассеянии, при таком подходе появляется необходимое граничное условие, состоящее в том, что при выходе на границу среды в направлении внутрь среды должна остаться только интенсивность падающего излучения:
ID (r,s)r = 0 Vs: (s • n i) > 0, (5)
которое в диффузионном приближении строго не может быть выполнено в силу условия в (4). Поэтому в уравнениях диффузионного приближения вместо (5) используются при-
ближенные граничные условия [14], например равенство нулю потока энергии от этой компоненты внутрь среды. Кроме того, в уравнениях для Ц/У) и ^(г) необходимо дополнительно учитывать эффективные источники, описывающие преобразование энергии падающего излучения, прошедшего в биоткань, в энергию диффузного излучения, и соответственно, сами уравнения усложняются [1, с. 198].
С учетом этих соображений представляется целесообразным использовать ненаправленное, изначально рассеянное излучение, с помощью конструкции, обеспечивающей фиксацию мощности, передаваемой в биоткань излучением (рис. 1). Кроме того, в этом случае отпадает необходимость разделения искомого решения и(г) на отдельные компоненты и, соответственно, уравнения диффузионного приближения принимают простую форму, а граничные условия - стандартный вид. Устройство для фиксации передаваемой мощности представляет собой камеру с хорошо отражающими излучение стенками, внутрь которой вводится лазерное излучение с мощностью Р^.
Рис. 1. Принцип фиксации мощности, передаваемой в биоткань излучением. 1 - биоткань, 2 - рассеиватель,
3 - камера с хорошо отражающими стенками
Внутри конструкции излучение предполагается изотропным, для чего в камере можно установить рассеиватели, а также использовать диффузно отражающие стенки, покрытые соответствующим материалом, например, на основе сульфата бария. £ - площадь поверхности ткани, на которую падает пучок рассеянного излучения, £ - площадь стенок, отражающих излучение, Ят - полный коэффициент отражения по мощности от биоткани,
включающий диффузное (из глубины) и френелевское (от границы раздела) отражение, Я - коэффициент отражения от стенок (близкий к единице). В этом случае из рисунка видно, что Рг = Р - Р= Рь. Таким образом, в стационарном состоянии мощность, передаваемая от
устройства в биоткань, оказывается не зависящей от характеристик отражения излучения биотканью (показателя преломления и других оптических параметров, степени изотропности излучения, шероховатости поверхности и т.п.), а определяется только мощностью источника Р .
В дальнейшем будем предполагать наличие осевой симметрии (r - расстояние от оси симметрии). Более подробно баланс энергии описывается следующим образом. Если F -плотность потока падающего на поверхность биоткани излучения, F - плотность потока
прошедшего границу раздела с биотканью падающего излучения и вошедшего в биоткань, FOMÎ(r) - плотность потока диффузного отражения из глубины биоткани через поверхность
биоткани наружу, Fz(r) - результирующая плотность потока внутрь биоткани, F - плотность потока френелевского отражения от границы раздела, P с индексами - соответствующие мощности, то
Fin = Fi - Fsp = Fi ■ TF (n)
FZ (r ) = Fin - Fout(r X
Pi = PL + Pr = PL + Psp + Pout,
P = Rt ■ P, , P = PL
i
(1-Rt )
рг = рь.
В общем случае некоторая часть излучения неизбежно поглощается стенками. При учете коэффициента отражения от стенок, не равного единице, можно получить следующие выражения:
Р
P =
- L
S '
(1-Rt ) + (1-R)S
ST
р=_р_
(1 — )
Таким образом, если А и А близки по порядку величины, то устройство позволяет фиксировать передаваемую в биоткань мощность излучения мощностью источника Рь, вне зависимости от коэффициентов френелевского и диффузного отражения от ткани, при условии (1 — Я) /(1 — Я) «1, которое для биоткани несложно выполнить, поскольку (1 — Я ) составляет величину порядка единицы [2].
Теперь, поскольку падающее излучение по предположению является ненаправленным, то в данном случае Ц^г) = и(г), соответственно в основном уравнении диффузионного приближения [1, с. 198] нет источников (ни естественных, ни эффективных), и оно принимает простой вид:
V 2и (г)—(г) = 0,
2 , (6)
/е^Г = 3/а /г ,
/г =/а (1 — ё) = /а +/1,
где ц - транспортный коэффициент, - приведенный коэффициент рассеяния.
Уравнение для потока, ввиду отсутствия источников, также сильно упрощается:
Ж(г) =--!- VU (г).
Как известно, для однозначного решения требуется задать граничные условия. Везде далее передаваемая излучением в биоткань мощность считается фиксированной (рис.1). Для конкретности, с целью сравнения решения в диффузионном приближении с результатами статистического моделирования, рассмотрим расчетную область, изображенную на рис. 2. Граничные условия для такой области будут иметь вид:
u(г = r2,z) = 0, z е [0,zmax], u(г, zтах) = 0, г е[0, д2], дп
— (г,z = 0) = 0, г е [r1,я,], (7)
дz
^(г,z = 0) = -3^г(г), г <rl. дz
Р1
зеркальная поверхность
Р2=0 ■ ^2 = 0
и=0 ^таж и=0
и=0
Рис. 2. Расчетная область для сравнения результатов моделирования
В первых двух уравнениях (7) учтено, что вся передаваемая биоткани мощность поглощается в рассматриваемой области, если R2 и гтах выбраны достаточно большими, что проверяется моделировнием методом Монте-Карло для заданных значений коэффициентов поглощения и рассеяния данной биоткани. Третье уравнение обусловлено условием зеркального отражения. В последнем уравнении фигурирует плотность потока (результирующего) входящей в биоткань энергии. Как уже отмечалось, полная входящая в биоткань мощность является фиксированной величиной:
Я^ (г) dSl = Рь,
2
Sl
^ = пЯ1
Моделирование методом Монте-Карло показало, что даже если падающее на поверхность биоткани излучение является ламбертовым с независящей от г плотностью падающего потока, в плотности потока входящей (передаваемой в биоткань) мощности появляется радиальная зависимость в силу наличия такой зависимости у диффузно-отраженной мощности, выходящей из глубины биоткани. При этом результирующая входящая мощность в центре пучка будет меньше, чем на краю, т.к. там больше диффузное отражение. Для моделирования такой зависимости можно использовать, например, степенную функцию
(г) = Рь - а
1 + с
V
г \ г
V Я1)
(8)
В (8) достаточно выбрать два параметра - Ь и с, поскольку значение параметра а определится из условия нормировки на полную входящую мощность:
1 Ь + 2 1
а =-т-= — Р-
п Я? Ь + 2 + 2 с ^
Легко видеть, что значение освещенности на поверхности в центре пучка (г=0) определяется параметром Р, при заданном значении которого будет постоянным и амплитудный коэффициент к , через который выражается параметр с:
а ка =ЬР),
а 2 р
С = (Ь + 2)-ка.
При больших Ь модельное решение не меняется при одновременном изменении Ь и с, если их отношение сохраняется. На рис.3 для сравнения представлены для заданной расчетной области (Я1=3мм, Я2=8мм, гтах=10мм) результаты моделирования для освещенности на поверхности и(г) методом Монте-Карло и моделирования в диффузионном приближении для трех режимов: один режим соответствовал с=0 (равномерное распределение входящей мощности по всей площади пучка), два других рассчитывались для ка=0,1333
(выбранного из условия близости модельного решения на поверхности в центре пучка к результату, даваемого моделированием Монте-Карло) с двумя выборами значения искомой функции на границе пучка (терминаторе).
Для режима, названного «модель 1», граничное значение выбиралось из условия непрерывности производной (условие гладкости), а для режима, названного «модель 2» граничное значение выбиралось из условия и(г=Я1, г=0) = 0,8и(0,0). Значения оптических параметров принимались равными характерным для подкожных тканей в ближнем ИК-
диапазоне: показатель преломления п=1,4, коэффициент поглощения А =1,5 см-1, коэффициент рассеяния л =150 см-1, фактор анизотропии £=0,8.
Рис. 3. Сравнение результатов моделирования
Как можно видеть, диффузионное приближение (с=0) с равномерным распределением ¥=еот1 дает большую ошибку в определении подповерхностной освещенности в центральной области, тогда как учет неравномерного по радиусу распределения входящей мощности в модельном решении по формуле (8) позволяет получить верное значение в центральной области вдали от края-терминатора. Однако при этом ухудшается аппроксимация для области вблизи терминатора, что обусловлено выбором модельной функции (8). Однако каков бы ни был выбор этой функции, обеспечить точное соответствие на всей границе невозможно, т.к. сами граничные условия имеют разрывный характер. С другой стороны, неопределенность в точках разрыва граничных условий можно использовать для построения численного решения (будем называть его граничным решением), лучшим образом аппроксимирующего значения искомой функции в тех или иных областях непрерывности граничных условий или даже в одной точке (где производятся измерения). Это проиллюстрировано на рис. 3, где выбор граничного значения и(г=Я1, 2=0) = 0,8Ц(0,0) позволил получить модельному решению «модель 2» лучшее соответствие для освещенности в большей части центральной области. Но тот же рисунок иллюстрирует, что вдали от пучка (область "чистой диффузии") модельное решение практически совпадает с диффузной аппроксимацией при ¥=еот1. Более того, моделирование Монте-Карло для иммерсионного случая показало, что в этой области нет зависимости и от показателя преломления биоткани. Это означает, что решение в этой области практически не зависит от конкрет-
ного распределения входящей мощности по поверхности и от показателя преломления биоткани по отношению к внешней среде, а определяется только величиной полной входящей в биоткань мощности (фиксированной и равной Р ) , а также внутренними параметрами поглощения и рассеяния излучения в биоткани. В п. 4 этот факт будет использован для решения обратной задачи - восстановления неизвестных параметров поглощения и рассеяния, задающих распределение световой энергии в биоткани.
4. Неинвазивное определение параметров распространения и поглощения монохроматического излучения в биоткани
Рассмотренный в п. 3 принцип фиксации мощности, передаваемой в биоткань излучением, позволяет расчетным путем (используя методы математического моделирования) осуществить неинвазивное определение двух важных параметров, описывающих распространение и поглощение излучения в диффузионном приближении - коэффициента поглощения ра и транспортного коэффициента р а значит, и коэффициента nf = 3 ¡ла p,tr в диффузионном уравнении. Для решения задачи неинвазивного определения параметров ра и р автором предложено использовать отмеченную в п. 3 независимость решения вне
пучка падающего излучения (вдали от него) от фактического распределения входящей мощности по поверхности (т.е. по радиусу пучка в условиях осевой симметрии), при фиксации мощности, передаваемой излучением в биоткань. Для нахождения граничного решения, как и в п. 3, используется численное моделирование переноса излучения в диффузионном приближении. Неизвестные параметры определяются итерационно, путем сравнения рассчитанных значений подповерхностной освещенности и мощности выходящего из биоткани излучения, с данными измерений соответствующих величин. При этом используется расчетная область, изображенная на рис. 4, с модификацией граничных условий на верхней границе: 1) в центре имеется зеркальная область радиусом R0<Rb 2) зеркальная кольцевая область R1<r<R2 заменена абсолютно поглощающей (показана черным цветом на рис. 4) - для моделирования свободного выхода излучения и/или поглощения его датчиком, измеряющим полный поток излучения (т.е. выходящую мощность Pout) из этой области. Таким образом, излучение входит в кольцевую область R0<r<R1 на поверхности биоткани, причем в математическом модели принято равномерное распределение входящей мощности по этой кольцевой области. Граничное условие для зеркальной области r < R1, как и ранее, задается равенством F=0, а для поглощающей области для нахождения граничных значений £/r(r) использовалось условие потоковой изотропности в условиях диффузионного приближения: F = F+r (с поправочным коэффициентом, близким к 1).
P,
P.
out
Rn
зерк.
<->
погл.
U=0
z
max
Fz=0
U(0,0) Fz(r)
Ri
R
U=0
Ur(r)
U=0
Рис. 4. Расчетная область и граничные условия для неинвазивного определения параметров распространения и поглощения излучения в диффузионном приближении
Измерение полной выходящей мощности Pout производится соответствующим фотоэлектронным прибором с учетом требования, что должно быть собрано всё излучение, выходящее из области R\<r<R2. Принцип измерения подповерхностной освещенности в центральной точке U(0,0) будет изложен ниже.
Так же, как это было сделано ранее, сравним результаты моделирования для заданной расчетной области и для биоткани с теми же оптическими параметрами, выбрав решение с равномерным распределение входящей мощности по кольцевой области Ro<r<Ri и с заданными условием потоковой изотропности граничными значениями Ur(r) (рис. 5).
Так же, как и в п. 3, равномерное распределение входящей мощности по площади пучка порождает существенную ошибку в определении подповерхностной освещенности в центральной части области входа излучения. Однако для значений освещенности в области "свободной диффузии" r<R0 это обстоятельство, как видно из рис. 5, роли не играет - в данном случае диффузионное приближение с "зеркальными" граничными условиями в области r<R0 демонстрирует хорошее совпадение с результатами статистического моделирования.
Граничное решение задавалось условием U(r=R1,0) =0,25 Umax, где Umax представляет собой максимальное значение U в кольцевой области R0<r<R1. При таком выборе граничного решения значение интеграла от выходящего потока по области R1<r<R2 оказывается близким к таковому для моделирования методом Монте-Карло. Последнее же представля-
ет собой измеряемую на опыте полную выходящую мощность из области Я1<г<Я2. При этом не требуется, чтобы было строгое соответствие значений и(г) в этой области наблюдаемым на опыте, т.к. задание полного потока на границе 2=0 не эквивалентно заданию распределения плотности потока (задача Неймана), которое неизвестно как в области пучка (там реальное распределение заменено в модели равномерным), так и вне его (там принимается потоковая изотропность). Но для определения значений параметров достаточно использовать какую-либо однозначно определяемую граничными условиями модельную функцию и(г, 2), при заданных параметрах биоткани дающую те же (или близкие) результаты измерений. Тем не менее, как видно из рис. 5, при фиксированной входящей мощности и использованных в модели граничных условиях, модельное решение оказывается достаточно близким к наблюдаемой на опыте функции в областях вне пучка.
Вт/см2
■ Монте-Карло -модельное решение — г = Р0 ......г =
Рис. 5. Модельное решение для определения параметров
Процедуру восстановления неизвестных параметров с использованием модельного решения наглядно демонстрирует рис. 6, где построены две линии, которые соответствуют удовлетворяющим экперименту (результатам измерений) зависимостям Значения параметров находятся по точке пересечения графиков. Эта процедура вполне аналогична используемой в предложенном автором методе неинвазивного измерения перфузии крови и коэффициента теплопроводности для живых тканей [15]. При построении соответствующих графиков, линия, отвечающая измерению и(0,0), имеет малый наклон, а ли-
ния, соответствующая измерению Р01а должна представлять собой гиперболу, так как, согласно диффузионному уравнению (6)
2
сот1
Ма =■
Метод проверялся при значениях £=0,8; 0,85; 0,9. Было выяснено, что наибольшую точность (погрешность метода на уровне 5% или менее) удается получить при £ = 0,8 (типичное значение для мягких биотканей в ближней ИК - области). При £=0,85 погрешность восстановления параметров составила ~ 10%, а при £=0,9 - около 15%. Если имеется дополнительная информация о величине £, то настройкой граничного решения в модели можно снизить погрешность метода в соответствующем диапазоне значений £.
Рис. 6. Графическое определение неизвестных параметров ^ и р.
Для £=0,8 и типичных для мягких биотканей значений и р радиус центральной
(зеркальной) области Р0=1мм оказывается вполне достаточным, что говорит о хороших условиях применимости диффузионного приближения. Однако для некоторых биологических сред параметр анизотропии может достигать больших значений (в видимом диапазоне). Для приемлемой погрешности при £>0,8 радиус Я0 следует увеличить до 2 мм, чтобы на этом расстоянии "диффузия успела отработать" (дать для £7(0,0) правильное значение), как это показано на рис. 7. При значениях £>0,9, по-видимому, диффузионное приближение становится неприемлемым, что и определяет границы применимости метода.
Рис.7. Выбор радиуса Я0 для биоткани с параметром £=0,9: (а) - неправильный; (б) - правильный
Для измерения значения подповерхностной освещенности в точке (0,0) можно предложить следующий метод, согласующийся с граничными условиями для рассматриваемой расчетной области, и поясняемый рис. 8.
Рис. 8. Принцип измерения подповерхностной освещенности или полного коэффициента отражения. 1 - биоткань с излучением от стороннего источника (при измерении освещенности) либо от одного источника PL0 ( при измерении коэффициента отражения), 2 - рассеиватель, 3 - камера с отражающими стенками, 4 - отражающая поверхность, ST - площадь освещаемой источником PL0 поверхности биоткани, SM - площадь отверстия для вывода измеряемой мощности
Это же устройство может использоваться для измерения полного коэффициента отражения от биоткани при освещении ненаправленным излучением, что, как будет показано ниже, позволяет получить дополнительно значения ещё двух параметров (коэффициента рассеяния и параметра анизотропии). Кроме того, устройство можно использовать и
для определения мощности самого лазерного источника (оптической мощности), если она изначально неизвестна.
В методе измерения освещенности или полного коэффициента отражения, как и в методе фиксации передаваемой мощности, используется камера с хорошо отражающими стенками, содержащая рассеянное лазерное излучение от дополнительного источника Рьо и от самой биоткани, но теперь излучение может выходить не только в биоткань, но и на измеритель мощности через дополнительное отверстие с площадью Зы. При измерении подповерхностной освещенности мощность этого лазерного излучения меняют так, чтобы добиться равенства РЬ0 = Ры. Это означает, что поток энергии из биоткани становится
равным нулю, т.е. реализуется "зеркальное" граничное условие как на рис. 4 для центральной области г < Я0. В этих условиях подповерхностная освещенность найдется по формуле, получаемой из фотонной модели и связи освещенности с плотностью энергии:
Ц. =
4 Р,
ы
Г х \
Б
п
V Т1/п )
4 Ры 2 4 Ры
п2 =.1,96 (для п = 1,4), (9)
Бы Бы
ы
где Т - значение коэффициента прохождения мощности для ламбертова пучка излучения из воздуха в среду с показателем преломления п (см. п.1.4), а Т17 - в обратном направлении. В формуле (9) учтено, что:
т - 1т
Тп ^ п2 ±п
I
сг )
Более точными для измерений представляются иммерсионные системы, для которых п2 =1, либо использующие гель с соответствующим показателем преломления. Таким гелем можно покрыть только поверхность биоткани для устранения имеющихся неровностей (как на поверхности кожи), и тем самым фиксировать величину (Тп /Т1/п ).
Измерение коэффициента отражения от биоткани для пучка ненаправленного излучения производится с помощью того же устройства с использованием формулы
Ят =1
г Р \
1ь о ,1
Р V ы
Бм „ ч
. (10)
Наконец, саму мощность лазерного источника также можно измерить на том же устройстве, закрыв зеркальной стенкой выход излучения на биоткань, и тогда РЬ0 = Ры.
Получив значения ра и р их можно использовать для нахождения параметров р и £.
Для этого используется измеренное значение коэффициента отражения (полного или диффузного - для иммерсионных систем). В самом деле, обозначив уже известное значение (р1г - ра) = Ар, можно построить зависимость р (£):
Ам
М< = (Т-М) ■ о»
Далее, для трех значений £ = {0,8; 0,85; 0,9} и соответствующим им по формуле (11) значений р при известном с помощью моделирования Монте-Карло находим соответствующие значения коэффициента отражения Ят и по трем точкам строим график Ру^).
Найдя пересечение этого графика с горизонтальной прямой, задаваемой измеренным значением Ят, восстановим значение
Альтернативный, более быстрый способ получения значения коэффициента рассеяния, заключается в использовании формулы (11), в которую нужно подставить значение £ для заданной длины волны из зависимости, полученной эмпирическим путем для данной биоткани, если она известна. Например, для кожных тканей эта зависимость имеет вид [16]:
£ (Я) = 0,7645 + 0,2355
1 - ехр
Г Я-500 нм^ ч 729,1 нм у
Таким образом, все параметры, характеризующие рассеяние и поглощение монохроматического излучения в среде биоткани, оказываются измеренными. Знание этих параметров позволяет решать прямую задачу нахождения в биоткани светового поля, например, методом Монте-Карло, уже для любых интересующих условий.
Заключение
Особенности распространения ненаправленного излучения в сильно рассеивающих (мутных) средах, при фиксации передаваемой от излучения в среду мощности, дают возможность неинвазивного определения параметров диффузионного приближения, обычно применяемого для нахождения светового поля в таких средах. Важным примером сильно рассеивающих мутных сред являются непрозрачные биоткани в оптическом и ближнем ИК - диапазоне, измерение параметров распространения и поглощения излучения в которых является актуальной задачей.
По результатам описанного в статье исследования можно сделать следующие выводы:
1) При моделировании методом Монте-Карло распространения ненаправленного излучения, освещающего поверхность биоткани, необходимо использовать функцию статистического представления для угла падения фотонов.
2) При фиксации передаваемой от излучения в биоткань мощности, решение в областях вне пучка (особенно вдали от пучка) практически не зависит от коэффициента преломления биоткани и поверхностного распределения входящей мощности в пучке, а определяется только полной входящей мощностью и параметрами распространения и поглощения излучения внутри биоткани. Эту особенность можно использовать для неинвазивного определения значений параметров.
3) Для однозначного восстановления двух неизвестных параметров диффузионного приближения р и ра необходимо измерить две величины, в качестве которых удобно использовать подповерхностную освещенность и диффузно отраженную мощность.
4) Использование ненаправленного (рассеянного) излучения позволяет не только фиксировать передаваемую от излучения в биоткань мощность, но и измерить значение подповерхностной освещенности и полного коэффициента отражения.
5) При известных (или уже измеренных) значениях р^ и ра знание полного коэффициента отражения позволяет восстановить значения ещё двух параметров, характеризующих рассеивающую среду - коэффициента рассеяния р и
параметра анизотропии £. Более удобно находить коэффициент рассеяния, используя значение параметра анизотропии из эмпирической зависимости от длины волны для данной биоткани.
Список литературы
1. Исимару А. Распространение и рассеяние волн в случайно неоднородных средах. В 2 т. Т. 1. Однократное рассеяние и теория переноса: пер. с англ. М.: Мир, 1981. 280 с.
2. Макаров С.Ю. Статистическое моделирование переноса излучения и световые переходные характеристики многослойной биоткани // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2014. № 12. С. 350-364. DOI: 10.7463/1214.0738862
3. Макаров С.Ю. Статистическое представление угла падения фотонов в ненаправленном излучении, проходящая и передаваемая мощность // Научный обозреватель. 2015. № 10. С. 40-41.
4. Wilson B.C., Patterson M.S., Flock S.T. Indirect Versus Direct Techniques for the Measurement of the Optical-Properties of Tissues // Photochemistry and Photobiology. 1987. Vol. 46. P. 601-608. DOI: 10.1111/j.1751-1097.1987.tb04820.x
5. Kienle A., Lilge L., Patterson M.S., Hibst R., Steiner R., Wilson B.C. Spatially resolved absolute diffuse reflectance measurements for noninvasive determination of the optical scattering and absorption coefficients of biological tissue // Applied Optics. 1996. Vol. 35, iss. 13. P. 2304-2314. DOI: 10.1364/AQ.35.002304
6. Bashkatov A.N., Genina E.A., Tuchin V.V. Ch. 5. Tissue optical properties // In: Handbook of Biomedical Optics / ed. by D.A. Boas, C. Pitris, N. Ramanujam. London: CRC Press, Taylor and Francis group, 2011. P. 67-102.
7. Sun J., Fu K., Wang A., Lin A.W.H., Utzinger U., Drezek R. Influence of fiber optic probe geometry on the applicability of inverse models of tissue reflectance spectroscopy: computational models and experimental measurements // Applied Optics. 2006. Vol. 45, iss. 31. P. 8152-8162. DOI: 1Q.1364/AO.45.008152
8. Сетейкин А.Ю., Красников И.В. Применение метода Монте-Карло для задач биофо-тоники. Благовещенск: Изд-во АмГУ, 2014. 68 с.
9. Palmer G.M., Ramanujam N. Use of genetic algorithms to optimize fiber optic probe design for the extraction of tissue optical properties // IEEE Transactions on Biomedical Engineering. 2007. Vol. 54, iss. 8. P. 1533-1535. DOI: 10.1109/TBME.2006.889779
10. Ravikant Samatham Venkata. Determination of Optical Scattering Properties of Tissues Using Reflectance-mode Confocal Microscopy: PhD Diss. Department of Biomedical Engineering, OHSU, 2012.
11. Макаров С.Ю. Теплофизическая модель биоткани и её численная реализация // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2013. № 10. С. 85-96. DOI: 10.7463/1013.0645537
12. Оптическая биомедицинская диагностика: пер. с англ. В 2 т. Т. 1. / под. ред. В.В. Тучина. М.: Физматлит, 2007. 560 с.
13. Тучин В.В. Оптика биологических тканей. Методы рассеяния света в медицинской диагностике: пер. с англ. М.: Физматлит, 2013. 812 с.
14. Haskell R.C., Svaasand L.O., Tsay T.T., Feng T.T., McAdams M.N., Tromberg B.J. Boundary conditions for the diffusion equation in radiative transfer // Journal of the Optical Society of America A. 1994. Vol. 11, iss. 10. P. 2727-2741. DOI: 10.1364/JQSAA.11.002727
15. Макаров С.Ю. Принцип для неинвазивного измерения параметров стационарного теплообмена в живых тканях // Наука и образование. МГТУ им. Н. Э. Баумана. Электрон. журн. 2014. № 2. С. 233-246. DOI: 10.7463/0214.0695233
16. Dunaev A.V., Zherebtsov E.A., Rogatkin D.A., Stewart N.A., Sokolovski S.G., Rafailov E.U. Substantiation of medical and technical requirements for noninvasive spectrophotometry diagnostic devices // Journal of Biomedical Optics. 2013. Vol. 18, iss. 10. Art. no. 107009. DOI: 1Q.1117/1.JBO.18.10.107009
Science ¿Education
of the Baumail MSTU
Science and Education of the Bauman MSTU, 2015, no. 12, pp. 87-109.
DOI: 10.7463/1215.0826701
Received: 25.10.2015
Revised: 10.11.2015
© Bauman Moscow State Technical Unversity
Non-Directional Radiation Spread Modeling and Non-Invasive Estimating the Radiation Scattering and Absorption Parameters in Biological Tissue
S.Yu. Makarov
1,*
msu-postigmailju Volgograd State University, Volgograd, Russia
Keywords: radiative transfer, diffusion approximation, noninvasive methods
The article dwells on a development of new non-invasive measurement methods of optical parameters of biological tissues, which are responsible for the scattering and absorption of monochromatic radiation. It is known from the theory of radiation transfer [1] that for strongly scattering media, to which many biological tissues pertain, such parameters are parameters of diffusion approximation, as well as a scattering coefficient and an anisotropy parameter.
Based on statistical modeling the paper examines a spread of non-directional radiation from a Lambert light beam with the natural polarization that illuminates a surface of the biological tissue. Statistical modeling is based on the Monte Carlo method [2]. Thus, to have the correct energy coefficient values of Fresnel reflection and transmission in simulation of such radiation by Monte Carlo method the author uses his finding that is a function of the statistical representation for the incidence of model photons [3]. The paper describes in detail a principle of fixing the power transmitted by the non-directional radiation into biological tissue [3], and the equations of a power balance in this case.
Further, the paper describes the diffusion approximation of a radiation transfer theory, often used in simulation of radiation propagation in strongly scattering media and shows its application in case of fixing the power transmitted into the tissue. Thus, to represent an uneven power distribution is used an approximating expression in conditions of fixing a total input power. The paper reveals behavior peculiarities of solution on the surface of the biological tissue inside and outside of the incident beam. It is shown that the solution in the region outside of the incident beam (especially far away from it), essentially, depends neither on the particular power distribution across the surface, being a part of the tissue, nor on the refractive index of the biological tissue. It is determined only by the total input power and by parameters of radiation propagation and absorption inside biological tissues. In addition, a special choice of the boundary conditions allowed us to obtain the model solution for which the value of the biological tissue output power in the region that is outside of the beam possessed this property as well. There is a proposal to use these features of the model solutions for non-invasive defining the diffusion approximation parameters by the method of calculation similarly to that used for non-invasive measurement of
parameters of stationary heat exchange in living tissues [15]. Thus, to provide unambiguous recovery of two unknown parameters of diffusion approximation in accordance with the found model solution are used measurement results of two values that are subsurface illumination and total output power. Thus, the method does not require the spatial resolution of the diffusely reflected power, which can be advantageous for measurements on real living tissues as compared with methods using "point" nature of the measurement, for example using fiber optic sensors.
In addition to the effect of fixing the power transferred into the tissue, the use of non-directional radiation enables us to measure the value of the sub-surface illumination, as well as the value of the total coefficient of reflection from biological tissue. Being known total reflection coefficient with known (or having been already measured) parameters of diffusion approximation allows us to determine additionally two more parameters characterizing the scattering medium, i.e. a scattering coefficient and an anisotropy parameter. The obtained values of the main parameters (absorption coefficient, scattering coefficient, and anisotropy parameter) provide the problem solution to find distribution of radiation of emerging heat sources for a given biological tissue already under any conditions using, for example, the Monte Carlo method.
Thus, the conducted study has shown that the non-invasive measurement of parameters of radiation propagation and absorption in the visible and near infrared range is possible and can find application in biomedical engineering using laser radiation.
References
1. 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.)
2. Makarov S.Yu. Statistical Modeling of Radiative Transfer and Transient Characteristics for Multilayer Biological Tissue. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2014, no. 12, pp. 350-364. DOI: 10.7463/1214.0738862 (in Russian).
3. Makarov S.Yu. Statistical representation of incidence angle of photons in non-directional radiation: transmittance power and transferred power. Nauchnyi obozrevatel', 2015, no. 10, pp. 40-41. (in Russian).
4. Wilson B.C., Patterson M.S., Flock S.T. Indirect Versus Direct Techniques for the Measurement of the Optical-Properties of Tissues. Photochemistry and Photobiology, 1987, vol. 46, pp. 601-608. DOI: 10.1111/j .1751-1097.1987.tb04820.x
5. Kienle A., Lilge L., Patterson M.S., Hibst R., Steiner R., Wilson B.C. Spatially resolved absolute diffuse reflectance measurements for noninvasive determination of the optical scattering and absorption coefficients of biological tissue. Applied Optics, 1996, vol. 35, iss. 13, pp. 2304-2314. DOI: 10.1364/A0.35.002304
6. Bashkatov A.N., Genina E.A., Tuchin V.V. Ch. 5. Tissue optical properties. In: Boas D.A., Pitris C., Ramanujam N., eds. Handbook of Biomedical Optics. London, CRC Press, Taylor and Francis group, 2011, pp. 67-102.
7. Sun J., Fu K., Wang A., Lin A.W.H., Utzinger U., Drezek R. Influence of fiber optic probe geometry on the applicability of inverse models of tissue reflectance spectroscopy: computational models and experimental measurements. Applied Optics, 2006, vol. 45, iss. 31, pp. 8152-8162. DOI: 10.1364/A0.45.008152
8. Seteikin A.Yu., Krasnikov I.V. Primenenie metodaMonte-Karlo dlya zadach biofotoniki [Use of Monte-Carlo method for the tasks of biophotonics]. Blagoveshensk, AmSU Publ., 2014. 68 p. (in Russian).
9. Palmer G.M., Ramanujam N. Use of genetic algorithms to optimize fiber optic probe design for the extraction of tissue optical properties. IEEE Transactions on Biomedical Engineering, 2007, vol. 54, iss. 8, pp. 1533-1535. DOI: 10.1109/TBME.2006.889779
10. Ravikant Samatham Venkata. Determination of Optical Scattering Properties of Tissues Using Reflectance-mode Confocal Microscopy. PhD Diss. Department of Biomedical Engineering, OHSU, 2012.
11. Makarov S.Yu. Thermo-physical model of bio-tissue and its numerical implementation. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2013, no. 10, pp. 85-96. DOI: 10.7463/1013.0645537 (in Russian).
12. Tuchin V.V., ed. Handbook of Optical Biomedical Diagnostics. Bellingham, Washington, SPIE Press, 2002. 1093 p. (Russ. ed.: Tuchin V.V., ed. Opticheskaya biomeditsinskaya diagnostika. V2 t. T. 1. Moscow, Fizmatlit Publ., 2007. 560 p.).
13. Tuchin V. Tissue optics. Light scattering methods and instuments for medical diagnostics. SPIE Press, 2007. (Russ. ed.: Tuchin V.V. Optika biologicheskih tkanei. Metody rasseyaniya sveta v medicynskoi diagnostike. Moscow, Fizmatlit Publ., 2013. 812 p.)
14. Haskell R.C., Svaasand L.O., Tsay T.T., Feng T.T., McAdams M.N., Tromberg B.J. Boundary conditions for the diffusion equation in radiative transfer. Journal of the Optical Society of America A, 1994, vol. 11, iss. 10, pp. 2727-2741. DOI: 10.1364/JQSAA.11.002727
15. Makarov S.Yu. A principle for the noninvasive measurement of steady-state heat transfer parameters in living tissues. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2014, no. 2, pp. 233-246. DOI: 10.7463/0214.0695233 (in Russian).
16. Dunaev A.V., Zherebtsov E.A., Rogatkin D.A., Stewart N.A., Sokolovski S.G., Rafailov E.U. Substantiation of medical and technical requirements for noninvasive spectrophotometric diagnostic devices. Journal of Biomedical Optics, 2013, vol. 18, iss. 10, art. no. 107009. DOI: 1Q.1117/1.JBO.18.10.107009
