ЖЭТФ, 2022, том 161, вып. 6, стр. 779-787

МОДЕЛИРОВАНИЕ ПЕРЕНОСА ИЗЛУЧЕНИЯ В ТЕРМИНАХ

УРАВНЕНИЯ БЕТЕ - СОЛПИТЕРА ДЛЯ ДВУХСЛОЙНЫХ

СИСТЕМ БИОТКАНЕЙ

В. Л. Кузьминa*, Ю. А. Жаворонковb**, С. В. Ульяновb***, А. Ю. Вальковa,b****

^a Санкт-Петербургский политехнический университет Петра Великого

195251, Санкт-Петербург, Россия

^b Санкт-Петербургский государственный университет

198504, Санкт-Петербург, Россия

Поступила в редакцию 23 декабря 2021 г.,

после переработки 29 декабря 2021 г.

Принята к публикации 30 декабря 2021 г.

Интенсивность обратного рассеяния излучения ближнего инфракрасного диапазона рассчитана для двух-

слойной модели сильно неоднородной среды, которую можно рассматривать как систему биотканей

«череп-мозг». На основе уравнения Бете - Солпитера описана процедура моделирования многократного

рассеяния методом Монте-Карло для двухслойной случайно-неоднородной среды. В качестве индика-

трисы однократного рассеяния использована фазовая функция Хеньи - Гринштейна. Рассчитаны зависи-

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

излучения и приемником. Вид этих зависимостей оказался чувствительным к изменению таких парамет-

ров системы, как анизотропия индикатрисы рассеяния, толщина слоев, длина волны лазерного излуче-

ния. Эта особенность может быть использована в медицинской диагностике. Предложен альтернативный

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

что, начиная с расстояния между источником и приемником порядка нескольких транспортных длин,

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

DOI: 10.31857/S0044451022060013

чатки, диагностики венозных повреждений, травм

EDN: DTRSKT

мозга и мягких тканей [2, 9, 11-14].

1. ВВЕДЕНИЕ

Для корректного описания состояния мозга по

Методы рассеяния света в настоящее время ак-

данным обратно рассеянного лазерного излучения

тивно используются для исследования биотканей,

обязательно требуется дополнительный учет нали-

прежде всего в связи с медицинской диагности-

чия черепа, т. е. для теории — как минимум двух-

кой [1-4]. Важную область здесь составляет быст-

слойная модель головы «череп-мозг» [15, 16]. За-

рая «полевая» диагностика различных травм голов-

дача определения оптических параметров модели

ного мозга [5-11]. Наличие так называемого «ок-

«череп-мозг» по данным эксперимента в биофанто-

на прозрачности» биологических тканей в ближнем

мах и теоретического описания интенсивности как

инфракрасном диапазоне привело к развитию ме-

аналитически — на основе уравнения переноса, так и

тода диффузной спектроскопии ближнего инфра-

численно — моделированием методом Монте-Карло

красного диапазона DNIRS (Diffuse Near-infrared

(МК), рассматривалась в работе [17]. Аналогич-

Spectroscopy, см. [1-4,7-9]). Указанный метод широ-

ная модель была описана в [18] и использовалась

ко используется для определения глубины и степени

для неинвазивного определения оптических свойств

повреждения кожных покровов и подкожной клет-

мозга в рамках DNIRS. В работе [19] показано, что

оптические параметры мозга можно получить, если

^* E-mail: kuzmin_vl@mail.ru

расстояние между источником и детектором вдоль

^** E-mail: zhavoronkov95@gmail.com

^*** E-mail: ulyanov_sv@mail.ru

поверхности черепа больше некоторого характерно-

^**** E-mail: alexvalk@mail.ru

го значения.

779

В. Л. Кузьмин, Ю. А. Жаворонков, С. В. Ульянов, А. Ю. Вальков

ЖЭТФ, том 161, вып. 6, 2022

Данная работа посвящена МК-моделированию

нечным k_f волновыми векторами, k_ij — волновой

обратного рассеяния в двухслойной среде. Выбор

вектор, направленный из точки r_j в точку r_i, k_ij =

для исследования двухслойной модели продиктован

= k₀r_ij/r_ij, rij = ri - rj. Здесь k0 = 2π/λ — вол-

желанием избежать громоздкости в теоретическом

новое число в вакууме, λ — длина волны. Пропа-

описании, эта модель является удобной базой для

гатор однократного рассеяния Λ(r) = r-2 exp(-μr)

верификации результатов для реальных многослой-

возникает из произведения двух комплексно-сопря-

ных моделей, учитывающих дополнительно слой

женных средних функций Грина скалярного поля,

скальпа, прозрачную цереброспинальную жидкость

μ = μ_s + μ_a — коэффициент экстинкции, μ_s и μ_a —

и сложную структуру мозга (см., например, [11]).

коэффициенты рассеяния и абсорбции соответствен-

Наш расчет интенсивности обратно рассеянного

но, p(k_f - k_i) — фазовая функция,

лазерного излучения основан на итерационном ре-

G(k_f - k_i)

шении уравнения Бете - Солпитера. Многократные

p (k_f - k_i) =

∫

(2)

интегралы, являющиеся членами разложения реше-

dΩ_f G(k_f - k_i)

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

дом Монте-Карло. Широко применяемая процедура

G(k)

— преобразование Фурье корреляционной

MCML (Monte Carlo Multi Layer) [20] основана на

функции диэлектрической проницаемости,

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

∫

зультате случайных блужданий. MCML-подход мо-

G(k) = d(r - r₀)e-ik·(r-r0)〈δε(r) δε^∗(r₀)〉.

(3)

дифицирован нами так, чтобы вклады в интенсив-

ность обратного рассеяния определялись на каждом

акте рассеяния [21,22], это приводит к существенно-

Уравнение (1) написано в лестничном приближе-

му снижению времени расчета. Также мы модифи-

нии, подразумевающем условие слабого рассеяния:

λ ≪ l_s = μ-1s. В формуле (1) мы использовали оп-

цируем процедуру обратного преобразования в ме-

тоде МК, явно учитывая наличие абсорбции в каж-

тическую теорему для скалярного поля,

дом порядке рассеяния.

∫

k⁴⁰

В разд. 2 данной работы приведено решение

μ_s =

dΩ_f G(k_f - k_i).

(4)

(4π)²

уравнения Бете - Солпитера в лестничном прибли-

жении для интенсивности обратного рассеяния в ви-

Пусть z — декартова координата, нормальная к

де ряда по кратностям рассеяния. В разд. 3 описа-

границе полубесконечной среды, r = (r_⊥, z), z > 0.

на процедура МК-вычисления интегралов методом

Для исходящего поля мы используем приближение

обратного преобразования. В разд. 4 в явном ви-

Фраунгофера, когда поле в дальней зоне является

де приведен метод MCML обратного преобразова-

произведением сферической волны и плоской вол-

ния кумулятивной функции распределения вероят-

ны, направленной в точку наблюдения [22]. Тогда

ности длины пробега фотона в двухслойной случай-

основная, некогерентная часть интенсивности об-

но неоднородной среде. В разд. 5 описана предло-

ратного рассеяния [23] представляется в виде

женная нами модификация метода обратного пре-

образования. В разд. 6 представлены результаты

∫

∞

∫

расчетов интенсивности обратного рассеяния, а в

J (s_i, s_f ) = 4π dz₁

dr₂Γ(r₂, r₁ | k_f , k_i) ×

разд. 7 проведен анализ полученных результатов и

z2>0

сделаны выводы.

× exp(-μ(s_f z₂ + s_iz₁)),

(5)

где s_i = 1/ cos θ_i, s_f = 1/ cos θ_f , θ_i — угол падения, а

2. УРАВНЕНИЕ БЕТЕ - СОЛПИТЕРА

θ_f — угол обратного рассеяния, отсчитываемый от

Мы описываем перенос излучения в случайной

обратного направления.

среде с помощью уравнения Бете - Солпитера [21,22]

Проводя итерации в уравнении Бете - Солпите-

ра, представляем интенсивность рассеяния как ряд

Γ(r₂, r₁ | k_f , k_i) = μ_sp (k_f - k_i) δ(r₂ - r₁) +

по кратностям рассеяния [22, 24]:

∫

+ μ_s dr₃p(k_f-k₂₃)Λ(r₂-r₃)Γ(r₃,r1 |k₂₃,k_i),

(1)

∑

J (s_i, s_f ) =

J(n)(s_i, s_f ),

(6)

n=1

где пропагатор Γ(r₂, r₁ | k_f , k_i) соответствует излу-

чению из точки r₁ в точку r₂ с начальным k_i и ко-

где вклад в рассеяние n-го порядка

780

ЖЭТФ, том 161, вып. 6, 2022 Моделирование переноса излучения в терминах уравнения Бете - Солпитера. . .

∫

∞

∫

J(n)(s_i, s_f ) = 4πμⁿ

dz₁

dr₂ . . . dr_nΛ(r₂₁) ×

drΛ(r)p(t) =

dξ dχ dφ.

(10)

2πμ

∏

× p(k₂₁ - k_i)

Λ(rj+1j )p(kj+1j - kjj-1) ×

После этого интеграл вычисляется как среднее по

j=2

выборке трех равномерно распределенных перемен-

ных ξ, χ, φ, где первые две принадлежат интервалу

× H(z_j)H(z_n)p(k_f - kn n-1)e-μ(siz¹+sfzⁿ).

(7)

[0, 1], а азимутальный угол φ — интервалу [0, 2π].

Функции Хевисайда H(z) обеспечивают исчезнове-

Приближая член n-го порядка J(n)(1, s_f ) сред-

ние этого вклада при вылете фотона из среды.

ним по выборке из N_ph падающих фотонов, имеем

N_ph

)

∑

J(n)(1, s_f) ≃

ⁿ p k_f -k(i)

e-μsf zni) ,

(11)

3. МОДЕЛИРОВАНИЕ МЕТОДОМ

n n-1

N_ph

МОНТЕ-КАРЛО

i=1

Опишем алгоритм моделирования для однород-

где веса

n для n > 1 задаются формулой

ной полубесконечный среды. Метод основан на из-

(

) (

)

вестной процедуре обратного преобразования [25],

⁽μ_s)ⁿ ∏

W(i)n =

H z(i)

H T -z(i)

(12)

который представляет пространственные интегралы

j≤n

итерационного ряда (6) по полубесконечному интер-

валу в виде интегралов по единичным интервалам.

Функция H(z)H(T -z) учитывает, что реальная сре-

Трехмерный пространственный интеграл в декарто-

да заполняет конечный слой толщиной T , 0 ≤ z ≤ T .

вых координатах r_j = (x_j , y_j , z_j) преобразуется в ин-

Отметим, что интенсивность рассеяния, рассчитан-

теграл по сферическим координатам (r, θ, φ), с нача-

ная с помощью (11), может быть интерпретирована

лом в rj-1. Ось z находится под прямым углом к по-

как среднее значение экспоненты exp(-μs_f zni)), ко-

верхности образца и уходит вглубь него. Форма про-

торая описывает затухание фотона, возвращающе-

пагатора однократного рассеяния Λ(r) показывает,

гося из среды к границе после n актов рассеяния.

что длина свободного пробега фотона имеет экспо-

Вес

n представляет собой случайное значение

ненциальное распределение вероятности с плотно-

многократного пространственного интеграла, полу-

стью f(r) = μ exp(-μr) на интервале r ∈ [0, ∞), где

ченного в результате итерации n-го порядка урав-

r = |rj - rj-1| — расстояние между точками j-го

нения Бете - Солпитера. Вычисляя его, моделируем

и (j - 1)-го порядков актов рассеяния. Для экспо-

стохастическую последовательность, или траекто-

ненциального распределения кумулятивная функ-

рию, точек рассеяния r₁, . . . , r_n. Переменная z(i)j —

ция распределения ξ = F (r) находится элементарно,

это расстояние до границы от j-го события рассея-

ния. Функция φ_BLB(zni)) = exp(-μs_f zni)) возникает

∫r

вследствие затухания рассеянного излучения, рас-

ξ = F(r) = f(r^′)dr^′ = 1 - exp(-μr).

(8)

пространяющегося от случайной точки n-го события

рассеяния zni) к границе. Оно зависит от локаль-

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

Обратное преобразование r = F-1(ξ) дает

жущегося к границе, и имеет вид в приближении

r = -μ-1 ln(1 - ξ) = -μ-1 lnξ^′,

(9)

Фраунгофера закона Бугера - Ламберта - Бера. Для

z). В фор-

однородной среды φ_BLB(z) = exp(-μs_f

где ξ и ξ^′ = 1 - ξ — случайные величины, равно-

муле (12) мы пренебрегли потерями энергии на от-

мерно распределенные в единичном интервале [0, 1].

ражение на границах между слоями и с вакуумом.

Аналогичное обратное преобразование выполняется

Эти потери можно учесть, если весовые коэффици-

с косинусом угла рассеяния: от t = cos θ переходим к

енты

n умножить на коэффициенты отражения

Френеля (см. [26]).

∫t

Метод МК широко используется для моделиро-

χ = 2π p(t^′)dt^′.

вания миграции фотонов в тканях и тканевых фан-

-1

томах, в основном в рамках известного алгоритма

Таким образом, трехмерный пространственный ин-

MCML [20]. В рамках MCML в сигнал вносят вклад

теграл по относительной координате r = r_j - rj-1

фотоны, выходящие из рассеивающей среды, что

преобразуется как

требует довольно большой выборки из-за того, что

781

В. Л. Кузьмин, Ю. А. Жаворонков, С. В. Ульянов, А. Ю. Вальков

ЖЭТФ, том 161, вып. 6, 2022

число случайно вышедших фотонов с заданной гео-

метрией задачи могут составлять очень малую до-

лю падающего света. В развитой в настоящей рабо-

те модификации каждый фотон вносит свой вклад

в сигнал при каждом акте рассеяния, пока не по-

кинет среду. Таким образом, объем выборки и, со-

ответственно, время вычислений, необходимое для

получения результатов, существенно уменьшается.

4. МЕТОД ОБРАТНОГО ПРЕОБРАЗОВАНИЯ

ДЛЯ ДВУХСЛОЙНОЙ СРЕДЫ

Рассмотрим неоднородную среду, в которой оп-

тические параметры зависят от положения фотона,

а именно — от декартовой координаты z по оси, нор-

мальной к границам.

В общем случае положим μ = μ(z). Мы предпо-

Рис. 1. Схематическое представление случайного пути фо-

лагаем, что на границе раздела слоев отражение от-

тона из начальной точки rj (черный кружок) в точку сле-

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

дующего рассеяния rj+1 (белый кружок); варианты дви-

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

жения фотона вглубь образца, cos θ > 0: (1a) — из среды

A в A, (1b) — из A в B, (2) — из B в B; варианты движения

ной среде, получаем, что экспоненциальная функ-

фотона в направлении поверхности, cos θ < 0: (3a) — из

ция распределения, определяющая затухание фото-

среды B в B, (3b) — из B в A, (4) — из A в A

на, движущегося из точки r₀ в точку r, может быть

представлена в следующем виде:

⎛

⎞

∫

а также с учетом того, движется ли фотон в слое A

exp(-μ|r - r₀|) → exp^⎝-

μ(z^′) dz^′⎠ .

(13)

cosθ

или B, либо пересекает границу этих слоев. Таким

образом, нам необходимо построить алгоритм опре-

Здесь θ — угол между направлением движения фо-

деления последовательных шагов, которые должен

тона и осью z, который должен быть определен за-

пройти фотон, начиная с точки z₀.

ранее. Заметим, что распределение зависит от на-

Пусть начальная точка z₀ выбрана в слое A для

чального положения фотона r₀. Определим плот-

фотона, движущегося вглубь среды, cosθ > 0, что

ность вероятности для пространственной координа-

соответствует путям (1a) или (1b) среди шести пу-

ты z, задающей новое положение фотона, формулой

тей, показанных на рис. 1. Для описания случай-

⎛

⎞

∫z

ного расстояния, пройденного фотоном между дву-

f (z, z₀) = C-10 exp ⎝-

μ(z^′) dz^′⎠ ,

(14)

мя последовательными актами рассеяния, опреде-

cosθ

лим функцию плотности вероятности [20]:

где C₀ — нормировочная постоянная.

Рассматривается двухслойная модель, состоя-

f (z, z₀) =

⎧

(

)

щая из слоев A и B. Слой A занимает область 0 <

⎪μ(A)

μ(A)

⎨

exp

(z - z₀)

z≤T_A,

< z < T_A, слой B — область T_A < z < T, где

cosθ

T = T_A + T_B — толщина данной двухслойной сис-

(

)

(15)

⁼⎪

⎪ μ(B)

μ(B)

темы. При построении обратного преобразования

⎩ξ_A

exp

(z - T_A)

z>T_A,

cosθ

мы рассматриваем среду B как полубесконечный

слой T_A < z, а в численных расчетах будем по-

лагать T_B конечным, но T_B ≫ T_A. Нам необходи-

где параметр ξ_A = exp (-μ(A)(T_A - z₀)/ cos θ) со-

мо получить кумулятивную функцию, зависящую

ответствует движению фотона от точки z₀ до гра-

от знака cosθ, и выполнить обратное преобразова-

ницы z = T_A. Для краткости мы опускаем зави-

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

симость введенных параметров от z₀. Интегрируя

видно из рис. 1, для каждого возможного направле-

плотность вероятности (15), получаем кумулятив-

ния движения фотона, вверх или вниз на рисунке,

ную функцию распределения

782

ЖЭТФ, том 161, вып. 6, 2022 Моделирование переноса излучения в терминах уравнения Бете - Солпитера. . .

F (z) =

МК-моделировании более точно учесть наличие аб-

⎧

(

)

сорбции и, в частности, отличие альбедо μ_s/μ от 1.

μ(A)

⎪

⎨1 - exp

(z - z₀)

z≤T_A,

Движение вглубь среды: z > z₀, cos θ > 0.

cosθ

(

)

(16)

Пусть фотон начинает двигаться в слое A, z₀ < T_A,

⎪

μ(B)

⎩1- ξ_A exp

(z - T_A)

, z>T_A.

в положительном направлении, cos θ > 0. Плотность

cosθ

вероятности для z₀ < T_A имеет вид

Значение кумулятивной функции рассматривается

далее как равномерно распределенная случайная ве-

f⁽⁺⁾(z, z₀) =

⎧

(

)

личина ξ = F (z, z₀). Выполняя обратное преобра-

μ_s(A)

μ(A)

⎪

exp

(z-z₀)

z≤T_A,

зование, т. е. определяя пространственную перемен-

⎨

C⁽⁺⁾ cosθ

cosθ

(

)

(21)

ную z как обратную функцию z = F-1(ξ), получаем

⎪

ξ_Aμ_s(B)

μ(B)

⎧

⎩

exp

(z-T_A)

, z>T_A.

cosθ

C⁽⁺⁾ cosθ

cosθ

⎪

z₀ -

ln(1 - ξ),

ξ≤1-ξ_A,

⎨

μ(A)

(+)

(17)

В нормировочную постоянную C⁽⁺⁾ = C₁

+C⁽⁺⁾²

⎪

cosθ

⁽1-ξ⁾

⎩T_A -

ξ>1-ξ_A.

вносят вклад два члена,

μ(B)

ξ_A

μ_s(B)

В случае z₀ > T_A и положительного направления

C⁽⁺⁾¹ = (1 - ξ_A)μs(A),

C⁽⁺⁾

=ξ_A

(22)

μ(A)

μ(B)

движения фотона, cos θ > 0, коэффициент рассея-

ния не меняется, μ(z) = μ(B), и применение обрат-

Теперь найдем обратную кумулятивную функ-

ного преобразования дает

цию z = (F⁽⁺⁾)-1(ξ). Выполняя схему обратного

преобразования, находим переменную z как функ-

cosθ

z=z₀ -

ln(1 - ξ).

(18)

цию случайной величины ξ:

μ(B)

Аналогично для фотона, движущегося из слоя B в

отрицательном направлении, cosθ < 0, получаем

⎧

(

)

cosθ

μ(A)

⎪

⎧

⎨z0-

1-C⁽⁺⁾ξ

, ξ≤ξ⁽⁺⁾,

cosθ

μ(A)

μ_s(A)

⎪

z₀ -

ln(1 - ξ),

ξ≤1-ξ_B,

(

)

⎨

(23)

μ(B)

⁼⎪

cosθ

C⁽⁺⁾

(19)

⎪

(1-ξ)

ξ>ξ⁽⁺⁾.

⎪

cosθ

⁽1-ξ⁾

⎩TA-

μ(B)

⎩T_A -

ξ>1-ξ_B,

C(+)2

μ(A)

ξ_B

Здесь

где ξ_B = exp (-μ(B)(T_A - z₀)/ cos θ). Для фотона,

движущегося в сторону границы z = 0 из слоя A,

μ_s(A)

ξ⁽⁺⁾ = F⁽⁺⁾(T_A) = (1 - ξ_A)

(24)

аналогично получаем

μ(A)C⁽⁺⁾

cosθ

— значение кумулятивной функции на границе сло-

z=z₀ -

ln(1 - ξ).

(20)

μ(A)

ев A и B. При z₀ > T_A, при движении вглубь, cosθ >

> 0, обратное преобразование дает

Таким образом, формулы (17)-(20) определяют из-

менение координаты z фотона в результате одного

cosθ

z=z₀ -

ln(1 - ξ).

(25)

акта рассеяния как функцию равномерно распреде-

μ(B)

ленной величины ξ для двухслойной среды. Число

Движение к поверхности: z₀ > z, cos θ < 0.

возможных случаев увеличивается как k(k+1) с рос-

Двигаясь в отрицательном направлении из слоя B,

том числа слоев k; так, в трехслойной среде имеется

cosθ < 0, мы получаем аналогичный результат:

12 различных вариантов движения фотона.

⎧

(

)

5. МОДИФИКАЦИЯ АЛГОРИТМА MCML

cosθ

μ(B)

⎪

⎨z0-

1+C(-)ξ

, ξ≤ξ(-),

μ(B)

μ_s(B)

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

(

)

(26)

⁼⎪

cosθ

C(-)

в выражении (7) для интенсивности рассеяния по-

⎪

⎩TA-

(1-ξ)

ξ>ξ(-),

μ(A)

рождает коэффициент μ_s(A) или μ_s(B). Мы разра-

C(-)

ботали процедуру обратного преобразования, вклю-

где

чающую эти множители в явном виде в функ-

μ_s(B)

ξ(-) = (1 - ξ_B)

(27)

цию плотности вероятности. Это позволит нам при

μ(B)C(-)

783

В. Л. Кузьмин, Ю. А. Жаворонков, С. В. Ульянов, А. Ю. Вальков

ЖЭТФ, том 161, вып. 6, 2022

Таблица. Типичные коэффициенты рассеяния и абсорбции рассмотренных биотканей μa и μ^s (в мм-1) для раз-

личных длин волн λ (в нм)

λ = 750

λ = 850

λ = 950

λ = 1050

μ_a

μ^′s

μ_a

μ^′s

μ_a

μ^′s

μ_a

μ^′s

Мозг

0.036

0.859

0.106

0.762

0.114

0.622

0.118

0.525

Череп

0.006

1.974

0.013

1.876

0.019

1.757

0.019

1.665

а нормировочная постоянная C(-) = C(-)1 + C(-)2

1-g²

p_HG(cos θ) =

(30)

представляет собой сумму двух слагаемых,

4π (1 + g² - 2g cosθ)3/2

μ_s(A)

где g = 〈cos θ〉 — параметр анизотропии однократ-

C(-)1 = (1 - ξ_B)μs(B),

C(-)

=ξ_B

(28)

μ(B)

μ(A)

ного рассеяния. Удобство фазовой функции ХГ свя-

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

В случае z₀ < T_A, при движении к поверхности,

выражается обратная кумулятивная функция рас-

cosθ < 0, обратное преобразование дает

пределения вероятности по углу рассеяния θ.

cosθ

Во многих оптических исследованиях биотканей

z=z₀ -

ln(1 - ξ).

(29)

μ(A)

и биофантомов используется приведенный коэффи-

циент рассеяния μ^′s (обратная транспортная длина),

связанный с коэффициентом рассеяния μ_s (обрат-

6. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

ной длиной рассеяния) соотношением μ^′s = (1-g)μ_s.

В данном разделе приводятся результаты рас-

На рис. 2 показаны зависимости интенсивности об-

четов интенсивности обратного рассеяния лазерно-

ратно рассеянного излучения от расстояния меж-

го излучения ближнего инфракрасного диапазона

ду источником фотонов и приемником на плоскости

двухслойной средой.

z = 0 при изменении параметра анизотропии g. На

Существуют обширные данные по оптическим

рис. 2а приведены результаты расчетов для различ-

параметрам биотканей (см., например, [27-30]). Па-

ных значений параметра g при одном и том же зна-

раметры для системы «череп-мозг», использован-

чении коэффициента рассеяния μ_s. В соответствии с

ные нами при моделировании, приведены в таблице.

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

Моделирование проводилось методом МК с ис-

менной остается интегральная интенсивность одно-

пользованием как традиционного метода обратно-

кратного рассеяния, в то время как приведенный ко-

го преобразования кумулятивной функции распре-

эффициент рассеяния μ^′s с изменением параметра g

деления вероятности длины свободного пробега фо-

меняется.

тона [20], так и предложенной в разд. 5 данной ра-

Как видно из рис. 2а, с ростом параметра g, т.е.

боты его модификации. В расчетах полагалось, что

с увеличением вытянутости индикатрисы однократ-

в среду, занимающую полупространство z ≥ 0, ла-

ного рассеяния вперед, растет и интенсивность об-

зерное излучение попадает вдоль направления оси

ратного рассеяния на любом расстоянии ρ между

z. Далее, в результате многократного рассеяния в

источником и приемником, исключая самую близ-

двухслойной среде возникает обратное рассеяние,

кую к источнику область. Такое поведение находит-

интенсивность которого рассчитывалась в точке на

ся в согласии с аналитическими расчетами [32, 33],

поверхности z = 0 на расстоянии ρ от точки входа

выполненными для изотропной и сильно анизотроп-

лазерного луча в среду. В расчетах кратность рас-

ной фазовых функций. Отметим, что на рис. 2а для

сеяния ограничивалась числом n = 3 · 10⁵, а объем

каждого значения g показаны результаты вычисле-

выборки был выбран N = 10⁶.

ний, выполненных по двум алгоритмам нахожде-

Результаты многократного рассеяния существен-

ния длины свободного пробега фотона, приведен-

но зависят от анизотропии рассеяния. Для изучения

ным в разд. 4 и 5 данной работы: а именно — за-

этой зависимости нормированную фазовую функ-

крашенные символы соответствуют классическому

цию p(cos θ), где θ — угол однократного рассеяния,

алгоритму MCML из работы [20], а белые симво-

мы описывали широко применяемой модельной фа-

лы — нашему алгоритму. На рис. 2б приведены ре-

зовой функцией Хеньи- Гринштейна (ХГ) [3,31]:

зультаты МК-моделирования для другого варианта

784

ЖЭТФ, том 161, вып. 6, 2022 Моделирование переноса излучения в терминах уравнения Бете - Солпитера. . .

lg(J)

–2

–4

–6

–8

–10

–12

–10

, мм

lg(J)

Рис. 3. Зависимости интенсивности обратного рассеяния

–2

в двухслойной модели череп (A) и мозг (B) от расстояния

источник-приемник для различных толщин слоя A: 0 (■),

3 мм (+), 5 мм (▴), 7 мм (×), 10 мм (♦), 150 мм (∗).

–4

Длина волны λ = 750 нм. Параметр анизотропии g = 0.9

–6

ком для системы «череп-мозг» с различной толщи-

ной верхнего слоя, т. е. черепа. Кривые для T_A = 0 и

–8

T_A = 150 мм показывают результаты обратного рас-

сеяния на полупространстве, содержащем лишь ве-

–10

щество мозга и черепа соответственно. Другие кри-

вые на рис. 3 показывают, как с увеличением тол-

, мм

щины верхнего слоя, т. е. черепа, изменяется на-

клон кривых. По результатам МК-моделирования

Рис. 2. Зависимости интенсивности обратного рассеяния

излучения, обратно рассеянного полупространством

в двухслойной модели череп (A) и мозг (B) от расстояния

только из одной биоткани, определялся угловой ко-

источник-приемник. Закрашенные символы соответству-

эффициент μ_eff в уравнении [34]

ют классическому алгоритму, белые — модифицированно-

му. Значения параметра анизотропии g = 0.8 (□ и ■), 0.9

ln[ρ²I(ρ)] = -μ_eff ρ + I₀,

(31)

(△ и ▴), 0.95 (♦ и ♦). Длина волны λ = 750 нм. Толщина

слоя A — TA = 5 мм. Графикам а и б отвечают варьируе-

справедливом при ρ ≫ ltr = 1/μ′s, т. е. в случае

мые параметры μ^′s и μs соответственно

применимости диффузионного приближения пере-

носа излучения. По данным кривой на рис. 3 с

T_A = 0, т. е. для полупространства с параметра-

роста параметра анизотропии g, а именно — неиз-

ми ткани мозга, было получено значение μ_eff (B) =

менной оставалась величина приведенного коэффи-

= 0.31 мм-1, которое очень близко к результату

циента рассеяния μ^′s, а коэффициент μ_s изменялся

μ_eff (B) = 0.305 мм-1, найденному по данным таб-

согласованно с изменением параметра анизотропии

лицы:

g. Из совпадения графиков на этом рисунке следует,

√

μ_eff =

3μ^′sμ_a.

(32)

что зависимость интенсивности обратного рассея-

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

По данным кривой на рис. 3 с T_A = 150 мм бы-

в рассмотренной двухслойной системе определяет-

ло найдено значение μ_eff (A) = 0.197 мм-1, кото-

ся приведенными коэффициентами рассеяния μ^′s(A)

рое также близко к значению μ_eff (A) = 0.188 мм-1,

и μ^′s(B), что свидетельствует о сформировавшемся

найденному по формуле (32) с данными из таблицы.

диффузионном режиме переноса излучения.

Таким образом, для данных тканей при моделирова-

На рис. 3 представлены результаты моделирова-

нии обратного рассеяния на расстояниях от источ-

ния зависимости интенсивности обратного рассея-

ника излучения больших 10 мм можно пользоваться

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

диффузионным приближением.

785

В. Л. Кузьмин, Ю. А. Жаворонков, С. В. Ульянов, А. Ю. Вальков

ЖЭТФ, том 161, вып. 6, 2022

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

lg(J)

ком излучения. Изучено влияние на интенсивность

–2

рассеяния изменения параметра анизотропии фазо-

–4

вой функции, толщины черепа и длины волны. Чув-

ствительность интенсивности обратного рассеяния к

–6

изменению параметров биоткани позволяет исполь-

зовать данные обратного рассеяния в медицинской

–8

диагностике. В расчетах мы применяли нашу мо-

–10

дификацию известной процедуры MCML [20], от-

личающуюся способом регистрации фотонов [22], а

–12

для моделирования длины свободного пробега фо-

тона мы пользовались как традиционным методом

–14

MCML, так и предложенной в данной работе его мо-

дификацией. Модифицированный метод позволил

, мм

существенно сократить время вычислений, а также

Рис. 4. Зависимости интенсивности обратного рассеяния

явно учесть в моделировании альбедо μ_s/μ.

в двухслойной модели череп (A) и мозг (B) от рассто-

яния источник-приемник для различных длин волн λ:

750 нм (+), 850 нм (♦), 950 нм (•), 1050 нм (×). Пара-

ЛИТЕРАТУРА

метр анизотропии g = 0.9. Толщина слоя A — TA = 3 мм

S. L. Jacques, Phys. Med. Biol. 58, R37 (2013).

На рис. 4 показано, как меняется зависимость

D. J. Davies, Z. Su, M. T. Clancy et al., J. Neuro-

интенсивности обратного рассеяния на двухслойной

trauma 32, 933 (2015).

среде от расстояния источник-приемник для раз-

В. В. Тучин, Оптика биологических тканей. Ме-

ных длин волн. Значительное отличие кривой, по-

тоды рассеяния света в медицинской диагности-

строенной для λ = 750 нм, от трех других кривых

ке, IPR Media, Москва (2021).

связано с большим отличием коэффициента абсорб-

ции и согласуется с формулой (32) для μ_eff .

А. Н. Башкатов, А. В. Приезжев, В. В. Тучин, КЭ

41, 283 (2011).

D. K. Joseph, T. J. Huppert, M. A. Franceschini,

7. ЗАКЛЮЧЕНИЕ

and D. A. Boas, Appl. Opt. 45, 8142 (2006).

В работе проведен расчет интенсивности обрат-

M. Dehaes, P. E. Grant, D. D. Sliva et al., Biomed.

ного рассеяния лазерного инфракрасного излучения

Opt. Express 2, 552 (2011).

на двухслойной случайно неоднородной биоткани. В

J. Selb, D. A. Boas, S.-T. Chan et al., Neurophoton.

качестве основной модели рассматривалась система

1, 015005 (2014).

«череп-мозг». Для сравнения с реальными биомеди-

цинскими данными по обратному рассеянию инфра-

A. Sabeeh and V. V. Tuchin, J. Biomed. Photon. Eng.

красного излучения головой человека двухслойная

6, 040201 (2020).

модель является очень упрощенной. Однако она мо-

жет служить предельным случаем для верификации

R. Francis, B. Khan, G. Alexandrakis et al., Biomed.

результатов, учитывающих слой скальпа, прозрач-

Opt. Express 6, 3256 (2015).

ную цереброспинальную жидкость и более слож-

10.

S. Mamani, L. Shi, T. Ahmed et al., J. Biophotonics

ную структуру мозга (см., например, [11]). Моде-

11, e201800096 (2018).

лирование переноса излучения в биоткани основы-

валось на уравнении Бете - Солпитера в лестнич-

11.

A. P. Tran, S. Yan, and Q. Fang, Neurophoton. 7,

ном приближении, решение которого представлено

015008 (2020).

в виде разложения по кратностям рассеяния. Каж-

12.

E. Zinchenko, N. Navolokin, A. Shirokov et al.,

дый член разложения является многократным ин-

Biomed. Opt. Express 10, 4003 (2019).

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

Монте-Карло. Были рассчитаны зависимости интен-

13.

E. S. Papazoglou, M. D. Weingarten, S. Michael et

сивности обратного рассеяния от расстояния вдоль

al., J. Biomed. Opt. 13, 044005 (2008).

786

ЖЭТФ, том 161, вып. 6, 2022 Моделирование переноса излучения в терминах уравнения Бете - Солпитера. . .

14. E. S. Papazoglou, M. T. Neidrauer, L. Zubkov et al.,

24. V. L. Kuzmin, M. T. Neidrauer, D. Diaz et al.,

J. Biomed. Opt. 14, 064032 (2009).

J. Biomed. Opt. 20, 105006 (2015).

25. L. Devroye, Non-Uniform Random Variate Genera-

15. S. Mahmoodkalayeh, M. A. Ansari, and V. V. Tuchin,

tion, Springer-Verlag, New York (1986).

Biomed. Opt. Express 10, 2795 (2019).

26. T. H. Pham, O. Coquoz, J. B. Fishkin et al., Rev.

16. M. S. Cano-Velazquez, N. Davoodzadeh, D. Halaney

Sci. Instrum. 71, 2500 (2000).

et al., Biomed. Opt. Express 10, 3369 (2019).

27. A. N. Bashkatov, E. A. Genina, V. I. Kochubey et

17. A. Kienle, M. S. Patterson, N. Dögnitz et al., Appl.

al., J. Phys. D Appl. 38, 2543 (2005).

Opt. 37, 779 (1998).

28. A. N. Bashkatov, E. A. Genina, V. I. Kochubey et

18. J. H. Choi, W. Martin, V. Yu. Toronov et al.,

al., Proc. SPIE 6163, 616310 (2006).

J. Biomed. Opt. 9, 221 (2004).

29. J. D. Johansson, J. Biomed. Opt. 15, 0570059 (2010).

19. M. A. Franceschini, S. Fantini, L. A. Paunescu et al.,

30. E. A. Genina, A. N. Bashkatov, D. K. Tuchina et al.,

Appl. Opt. 37, 7447 (1998).

Biomed. Opt. Express 10, 5182 (2019).

20. L. Wang, S. L. Jacques, and L. Q. Zheng, Comput.

31. T. Durduran, R. Choe, W. B. Baker et al., Rep.

Meth. Prog. Bio. 47, 131 (1995).

Progr. Phys. 73, 076701 (2010).

21. В. Л. Кузьмин, А. Ю. Вальков, Письма в ЖЭТФ

32. T. M. Nieuwenhuizen and J. M. Luck, Phys. Rev.

105, 261 (2017).

E 48, 569 (1993).

22. В. Л. Кузьмин, А. Ю. Вальков, Л. А. Зубков,

33. V. L. Kuzmin and A. Yu. Valkov, J. Quant. Spectrosc.

ЖЭТФ 155, 460 (2019).

Radiat. Transf. 272, 107760 (2021).

23. V. L. Kuzmin, V. P. Romanov, and E. V. Aksenova,

34. D. Tamborini, P. Farzam, B. B. Zimmermann et al.,

Phys. Rev. E 65, 016601 (2001).

Neurophoton. 5, 011015 (2017).

787