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

ПРИБЛИЖЕННЫИ МЕТОД РАСЧЕТА РАСПРЕДЕЛЕНИЯ ЭНЕРГИИ ОПТИЧЕСКОГО ИЗЛУЧЕНИЯ В МНОГОКРАТНО РАССЕИВАЮЩИХ СРЕДАХ

И.А. Братченко, В.П. Захаров Самарский государственный аэрокосмический университет им. С. П. Королева

Аннотация

В работе представлен метод приближенного решения уравнения переноса для неоднородной многократно рассеивающей среды, описываемой анизотропной фазовой функцией Хени-Гринштейна. Поток распространяющегося в среде излучения представлен в виде коллимированной и диффузной компонент, в которых транспортное уравнение распадается на две части. Уравнение для коллимированной компоненты легко разрешается и приводит к виду, описываемому законом Бугера, а интегро-дифференциальное уравнение для диффузной компоненты сведено к системе зацепляющихся уравнений в частных производных для порядков разложения диффузного потока по малому параметру. Каждое уравнение системы фактически является задачей Коши с граничными условиями равенства производных соответствующих порядков. Проведено сравнение с численным решением исходной задачи методом Монте-Карло, показана быстрая сходимость предложенного приближенного метода.

Ключевые слова: транспортное уравнение переноса излучения, метод Монте-Карло, приближенные методы, многократно рассеивающая среда.

Введение

В наши дни проблема исследования свойств сильно рассеивающих сред относится к одному из наиболее быстро развивающихся разделов физической оптики. Интерес к данным исследованиям обусловлен важностью определения параметров среды и распределения электромагнитного поля в нем для целого ряда приложений: распространение излучения в аэрозольной среде и создание лидаров [1], кристаллооптика [2], биологические ткани и физическая медицина [3]. Следует отметить, что задача распространения излучения в многократно рассеивающих средах существенно усложняется при учете структурных микро- и макронеоднородностей среды. Фактически реконструкция распределения неоднородностей по результатам физического эксперимента является основой 3Б визуализации для всех экспериментальных методов диагностики многократно рассеивающих сред.

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

Для изучения распространения оптического излучения в многократно рассеивающих средах наиболее часто используются численные методы, в частности, метод Монте-Карло [4]. Метод Монте-Карло хорошо исследован применительно к одномерным и/или двумерным представлениям много-

кратно рассеивающей среды [5], что заведомо накладывает довольно существенные ограничения на применимость полученных результатов. Трехмерная задача в данном подходе решена только для случая кусочно-неоднородной среды [6]. Усложнение топологии исследуемой среды или ее существенная неоднородность приводит к значительному увеличению количества итераций, требуемых для получения корректного результата. Это, в свою очередь, требует постоянного контроля сходимости численных результатов, что во многих случаях труднореализуемо.

Альтернативным подходом является использование транспортного уравнения переноса излучения [7], являющегося выражением закона сохранения энергии поля в многократно рассеивающей среде. Однако данное уравнение допускает точное решение только для частного случая однородной среды с малой анизотропией. В связи с этим наиболее часто используют диффузионное приближение, позволяющее получить аналитическое решение, которое является по сути оценкой сверху для рассматриваемой задачи. Ценность такого решения по сравнению с методом Монте-Карло заключается в возможности быстрого анализа экспериментальных результатов. Однако применительно к проблеме 3Б визуализации и реконструкции распределения неоднородностей такой подход является заведомо неприемлемым.

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

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

= -ц I(г, /) +

Ц Г л л, л. . л

+ I р(5 , 5 )I(г, 5 ),Ю +х(г, 5 ) .

4п 4.

Здесь изменение потока света 1(г,8) в точке объема среды г с направлением распространения потока 5 обусловлено тремя причинами: затуханием в результате рассеяния (ц5) и поглощения (ца) вдоль пути ds (первое слагаемое) с коэффициентом затухания ц = ца + ц5; приростом энергии (второе слагаемое) за счет рассеянного в объеме излучения и пропорционального коэффициенту рассеяния ц5; вкладом в распространяющееся в среде излучение от любых источников х(г, 5 ) .

Вид решения уравнения (1) во многом определяется заданием граничных условий и анизотропией исследуемой среды, за которую отвечает фазовая функция рассеяния р^’). Выберем в качестве фазовой функции р^^’) эмпирическую функцию Хе-ни-Гринштейна, адекватность применения которой для рассматриваемых сред исследована, в частности, в работе [9]:

р(0) =

(1 + я 2 - 2я ос8(0))3

где 9 - угол, а g - параметр анизотропии рассеяния. Если пренебречь эффектами люминесценции (X = 0 ), то, разлагая поток I в ряд Фурье, можно переписать уравнение (1) для спектральной интенсивности оптического поля. Вид уравнения (1) при этом не изменится ввиду отсутствия взаимодействия различных спектральных компонент поля в отсутствии источника. Аналогичный вид транспортное уравнение будет иметь и для монохроматического лазерного излучения.

В этом случае граничное условие на поверхности раздела воздух-среда можно записать в виде

10 =

где г5 - френелевское отражение, \ = cos(9), а Б(г) -поперечное распределение интенсивности электромагнитной волны, падающей нормально на поверхность раздела.

Выделим коллимированный поток 1соП, вклад в который вносят только фотоны, не испытывающие значительного рассеяния. Тогда полный поток I в объеме среды можно представить как

1

где 1а - диффузная компонента электромагнитного поля. При подстановке (4) в (1) транспортное уравнение распадается на две части:

д/

д5

• = -Ц 1со

^ = -Ц I, + 4^ | р(5,5')IddЮ + / (г, 0) .

Здесь функция

/с (г, 0) = I Р(^ -О1со1^Ю '

4п

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

Уравнение (5) легко решается и с учетом граничного условия (3) приводит к экспоненциальной зависимости интенсивности 1соП от эффективной

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

Ц 2

1сои =

выражающей закон Бугера.

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

1

ді,

ді,

дг г бш(0) д0

=-ц1, +-7- I p(0)Idd 0+

+ Ц- (1 - г, ) Р (г С08( 0)) (1 -2Я )ехр(-ц г ) . 4 (1 + Я2-2ЯСО8(0))3/2

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

а

уравнение (8) преобразуется к системе зацепляющихся дифференциальных уравнений для порядков разложения:

д1<0) + 1 д1Т

— = -ц С +

дг г sin(9) д9

+4- (1 - г5) ^ (г cos(9)) р(9) ехр(-цг sin 9)),

1

а = -ц 1Т +

дг г sin(9) д9

+— [р(9-9') 1(- -1)(9')а9' 4п;

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

д!Г

дг

= 0.

Для гауссового распределения падающего излучения:

Р

решение уравнения (9) для нулевого порядка диффузной составляющей имеет вид:

1

4ц П (1+Я2 -2я cos(9))3/2

ехр

Если подставить полученное решение нулевого порядка в (10), уравнение для первого порядка может быть найдено по методу Коши [6]:

q = q0e—(<) + е~ ¥ Aj ер (т) / (т)а т,

' (')=1 ^

t 1п t

р = рйеТ ^) + е-(t -11 еТ (т) 5(т)а т, Т^) = 2ц t,

где

р=

д1ао)

дг

q =

д1ао)

д9

(t) ехр ^ 1 +1п (t) - - ] ] р ( ),

5(t) = ц^^---^ ^2at 1п2 (t)-ц^ 1 + 1п(t)-“;

х^ () ехр ^ц/1 ^1 +1п (г1)- П ^ р ()).

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

Анализ решений для однородной среды

Для исследования сходимости предложенного приближенного метода рассмотрим распространение лазерного излучения с длиной волны X = 633нм в по-лубесконечной однородной среде, характеризуемой коэффициентом анизотропии рассеяния g = 0,8, коэффициентом поглощения = 0,03 мм-1, коэффициентом рассеяния ^ = 30 мм-1. Данные оптические параметры соответствуют дермису кожи [13].

Распределение интенсивности излучения в исследуемой среде с использованием предложенного метода и учетом первых двух порядков разложения представлено на рис. 1а. Здесь же для сравнения приведены расчетные данные, полученные на основе аналитического решения [10] для фазовой функции Эддингтона

р(^(9)) = -1 (2 / 5(1- ^(^ + (1- / )(1 + 3я ^(ОД, 4п

в которой учитывается только линейная часть анизотропии рассеяния. Видно, что использование фазовой функции Хени-Гринштейна (2) приводит к существенному расплыванию волнового пакета по поперечному сечению, что коррелирует с экспериментальными данными работы [11].

№111!

0 2 4 6 0 2 4 6

г, мм 2, мм

Рис. ¡.Распределение поглощенной энергии: а - предложенный метод, б - приближение Эддингтона

Следует отметить, что для анализируемой среды вблизи оси падающего излучения наибольший вклад (70%) в суммарную величину поглощенной энергии (4) вносит коллимированный поток излучения. Соответственно нулевой порядок диффузной компоненты составляет около 25% от суммарного потока, а первый порядок диффузной части не превышает 2%. Учет больших порядков не приводит к существенному изменению картины распределения потока поглощенного в среде излучения. По мере удаления от оси вклад от коллимационной части экспоненциально убывает, а соотношение между порядками диффузного потока медленно перераспределяется между нулевым и первым порядком. Так, для поперечного радиуса, соответствующего полуширине

п

г=+0

диффузного расплывания пакета, вклад нулевого порядка составляет 94%, а первого порядка - 5% от общего суммарного потока.

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

Сравнение приближенного и численного методов Для проверки адекватности полученных решений задача в аналогичной постановке была решена с помощью численного метода Монте-Карло в программной среде TracePro Expert [12]. Отметим, что как для приближенного решения транспортного уравнения, так и для метода Монте-Карло все вычисления производились без учета влияния внешней среды - воздуха. Распределение интенсивности в среде при использовании метода Монте-Карло представлено на рис. 2.

У>лш I' мВт _____^

0 5 г, мм

Рис. 2. Распределение поглощенной энергии

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

Сравнительный анализ полученных зависимостей показывает, что приближенное решение транспортного уравнения (12) - (14) при больших значениях анизотропии рассеяния ^ = 0,7.1) согласуется с численным решением Монте-Карло как качественно, так и количественно, в то время как при малой анизотропии (^ ~ 0,4^0,7) лучшее согласие наблюдается при использовании линейной фазовой функции Эддингтона.

При существенном уменьшении анизотропии среды и/или коэффициента рассеяния наблюдается только качественное согласие между прибли-

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

Рис. 3. Глубина распространения (а) и полуширины диффузного расплывания (б) рассеянного излучения в зависимости от коэффициента анизотропии:

1 - приближение Эддингтона, 2 - метод Монте-Карло,

3 - предлагаемый метод

Анализ решений для неоднородной среды

Сравнительный анализ численного и приближенного решений для случая неоднородной среды проведен на примере кожи человека, представленной в виде модельной многослойной среды с оптическими параметрами, заданными в соответствии с работой [13].

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

Зависимость полуширины диффузного рассеяния в модельной среде от глубины ъ представлена на рис. 4.

Видно, что, как и для численных экспериментов с однородной средой, сохраняется хорошее совпадение результатов расчетов по обоим методам. Отклонение составляет не более 10%. При этом временные затраты на достижение сравнимой точности вычислений для выделенного набора параметров при тактовой частоте процессора 2,81 ГГц снижается до 4-5 минут по сравнению с несколькими часами, требуемыми при использовании метода Монте-Карло.

Рис. 4. Полуширина рассеянного излучения в зависимости от глубины: 1 - численное решение Монте-Карло,

2 - приближенное решение

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

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

Благодарности Работа выполнена при поддержке аналитической

ведомственной целевой программы «Развитие научного потенциала высшей школы»

и гранта РФФИ 08-02-99038.

Литература

1. Уравнения множественной регрессии для фракционных концентраций атмосферного аэрозоля / М. М Ку-гейко. С. А. Лысенко //Журнал прикладной спектроскопии, 73. -2006, -№6. -С. 807-812

2. Агранович, В. М. Кристаллооптика с учетом пространственной дисперсии и теория экситонов./ В.М. Агранович, В.Л. Гинзбург // -М. Наука, -1979г. -С.432.

3. Сравнительный спектральный анализ обратного рассеяния излучения растительной и живой тканью / Воробьева Е.В. [и другие]// Физика волновых процессов и радиотехнические системы. -2005. -№3. -Т.8. -С.69-74.

4. MCML - Monte Carlo modeling of light transport in multi-layered tissues / L. Wang, S.L. Jacques, L. Zheng // Computer methods and programs in Biomedicine -1995. -v.47, -P.131-146.

5. Моделирование распространения оптического излучения в сложной случайно-неоднородной среде методом Монте-Карло / С.Д. Словецкий // Радиотехника -1994. -№7. -С.654-671.

6. 3D визуализация многократно рассеивающих сред / В. П. Захаров, А. Р. Синдяева // Компьютерная оптика, -2007. -Т.31. -№4. -С. 44-52.

7. Численное решение краевых задач для уравнения переноса излучения в оптическом диапазоне./ И. П. Яро-венко // Вычислительные методы и программирование: новые вычислительные технологии, -7 -2006. -№1 С.93-104.

8. Исимару, А. Распространение и рассеяние волн в случайно - неоднородных средах./ А. Исимару // - М: Мир. -1981. -С.384.

9. Устойчивость фазовой функции Хенье-Гринштейна и быстрое интегрирование по путям в условиях многократного рассеяния света / В.М. Петникова, Е.В. Третьяков, В.В. Шувалов // Квантовая электроника. -2006. -Т.36. -№11. -С.1043-1047.

10. Light transport in tissue./ S. A. Prahl // -1988. - P.221.

11. Использование спектроскопии обратного диффузного отражения света для мониторинга состояния тканей при фотодинамической терапии, / А.А. Стратоников [и другие] //Квантовая электроника, -36. -2006, -№12. -С.1103-1110

12. TracePro user’s manual supplement, release 3.1. Lambda Research Corporation.

13. Tuchin, V. V. Tissue optics./ Tuchin V. V.// -2000. P.353.