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

УДК 535.361

МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ ОПТИЧЕСКОГО ИЗЛУЧЕНИЯ С МНОГОКРАТНО РАССЕИВАЮЩИМИ СРЕДАМИ

© 2010 В.П.Захаров, И.А Братченко, Е. В. Тимченко

Самарский государственный аэрокосмический университет

Поступила в редакцию 11.01.2010

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

ВВЕДЕНИЕ

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

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

Для изучения распространения оптического излучения в многократно рассеивающих средах наиболее часто используются численные методы, в частности метод Монте-Карло [4]. Метод Монте-Карло хорошо исследован применительно к одномерным и/или двумерным представлениям многократно рассеивающей среды [5], что Захаров Валерий Павлович, доктор физико-математических наук, профессор кафедры "Автоматические системы энергетических установок". Е -mail: Zakharov@ssau.ru Братченко Иван Алексеевич, аспирант, инженер кафедры "Автоматические системы энергетических установок". E-mail: Ud_liche@mail.ru

Тимченко Елена Владимировна, кандидат физико-математических наук, ассистент кафедры "Автоматические системы энергетических установок". E-mail: Vorobieva.82@mail.ru

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

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

ОСНОВНЫЕ УРАВНЕНИЯ

Для описания распространения излучения в многократно рассеивающей среде воспользуемся транспортным уравнением переноса излучения, которое в случае отсутствии выделенной поляризации может быть представлено в виде [8]:

dI (r,s, Л)

ds

Ms (r ,Л) 4n

= -p(r, Л)1 (r, s, Л) +

+

{ps (s, s')i (г, гла+хл(г, s), (1)

4n

где ¡л = ¡ла + ¡л,., ¡ла и ¡л, - коэффициенты поглощения и рассеяния, 8 и 8' - направления падающего и исходящего лучей, □ ' - телесный

угол в направлении 8', х(г, 5) - мощность излучения в единице объема.

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

Р(в) =

(1 - Е2)

(2)

(1 + е2 - 2есоз(в))3/2'

где в- угол рассеяния, а g - параметр анизотро пии. В качестве граничных условий на поверх ности раздела воздух - среда Г принималось:

I(г,5,Л)1 гП) 0 =

4 ' ' /\( 5П )<0

= 1К (г, 5, Л) + ЯI (г, 5, Л)

г еГ

. (3)

(5п) >0

где Я - оператор Френеля, а п - нормаль к поверхности раздела. Считая процесс флуоресценции изотропным, функцию источника %Л(Т, 5) можно представить в виде:

4п к

=

| алпЬЛЛпШЬЛп) \^1(ЫЛп), (4)

Л

4п

нат, слабо зависящей от азимутального угла. В этом случае диффузные потоки Iмогут быть разложены в ряд по малому параметру [10]:

=14 ^ 5).

(6)

Подставляя разложения (5-6) в уравнение (1), получаем решение в виде:

!10)(г,5) = (1 -г)п 1± ехр-тех)(г • (5-

т (0) 1/к

1

(г, 5) =—|с5К, (г')ехр(-(Тд -Т)) 4п 0 1 5

Ф] (r, 5) = —|д (г', Л) ехр(-(т^ - т}()) X, (7) X | (г ' • г 'у£"-1)(г', г '')

4п

где индекс £ = (ех,/к), 50 - определяет направление падающего потока излучения с поперечным распределением интенсивности I , г5 - коэффициент френелевского отражения от границы раздела сред, а оптическая толщина определена как:

л

= |д(г ,Л?)сЬ.

где - квантовый выход флуоресценции к-го флуорофора при облучении среды излучением с длиной волны Лп .

МЕТОД МАЛОГО ПАРАМЕТРА

Выделим диффузную часть интенсивности

излучения в среде:

I = I (0) +1

ех ех ех,С/

I = I (0) +1 , (5)

1 /к 11к

где IX определяет вклад в интенсивность на длине волны Лп "баллистических" фотонов, не испы-

т(0)

тывающих значительного рассеяния; Iу - определяет вклад в интенсивность на длине волны л, за счет переизлучения поглощенной энергии (флуоресценции) в направлении г; а - диффузная часть интенсивности излучения на длине волны Л^. Учитывая, что в многократно рассеивающей среде распределение неоднородностей и флуорофоров носит стохастически-случайный характер, без существенной потери общности можно считать распределение оптических параметров среды медленно меняющейся функцией коорди-

Задавая поперечное распределение I^ падающего излучения и используя рекуррентные соотношения (7), можно найти аналитическое решение для диффузного потока первого порядка IX и всех последующих порядков разложения (6), что в конечном итоге позволяет определить ^ . Подстановка данного распределения в выражение (4) дает решение для потока 7(0), который рекуррентным образом определяет все последующие порядки разложения /^к""1.

ПРИБЛИЖЕНИЕ СФЕРИЧЕСКИХ ГАРМОНИК

Одним из альтернативных подходов к решению уравнения переноса является использование разложения интенсивностей и фазовой функции в ряд Лапласа по сферическим гармоникам УМ (5) [11]. Фактически такое разложение интенсивности эквивалентно учету дипольного, квадрупольного и последующих порядков рассеяния и флуоресценции. В этом случае, применяя теорему сложения сферических гармоник, из уравнения (1) получаем:

» I ( » г ]

(5 -У+(д +дж -хУ"(5)- д I ^Г (5)Е ЕрД^С [=0

1=0 т=-1[ 4п I=0" =-1 ]

"=1

где 1Ы и ру1 - коэффициенты разложения интенсивности и фазовой функции соответственно. Умножая (8) на УМ (у), интегрируя по 8, используя ортогональность сферических гармоник и ограничиваясь дипольным приближением (Ь < 1) уравнение (8) может быть сведено к уравнению диффузии для спектральной плотности энергии и (Г, Л):

является плоскопараллельной многослойной многократно рассеивающей средой, основным флуорофором которой является хлорофилл. Морфология и оптические параметры среды задавались в соответствии с данными работы [12], а в качестве критериев исследования использовались дифференциальные коэффициенты обЛ

V 2и (г, Л) - ки (г, Л) = -цх(Г, Л) , (9) ратного рассеяния (К = —— ) и флуоресценции

где

к=л

в

п = Л в =-1-

В 3(ЛУ + Ла )

г 1 4п

и (Г ,Л) = [-1 (Г, У Л)^ = — 100(г ,Л), а

и

скорость света в среде.

Считая инициирующее излучение гауссовс-ким и применяя к (9) прямое и обратное преобразование Ханкеля нулевого порядка, в конечном итоге получаем:

да к да I-

и(Г,Л) = \акг 1 Л(к^)\сЫк±,д,Л)вхр-1к +к1 2 -ф

о к + к 1 0 '

(10)

где ^ - функция Бесселя нулевого порядка, О (кТ ,2 - д) - функция Грина.

СРАВНИТЕЛЬНЫЙ АНАЛИЗ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ И ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ

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

( F =

Л

Л

■), где ЛЛ - полный поток излучения

на длине волны Л. Зависимости данных коэффициентов от суммарной концентрации хлорофилла СсЫ представлено на рис. 1. При использовании метода малого параметра распределение интенсивности рассенного излучения и флуоресценции в среде определялось с учетом первых десяти порядков разложения (7). Функциональная зависимость К(СсЫ), полученная с использованием метода малого параметра, имеет линейный характер, однако если численная модель Монте-Карло дает заниженное значение ( < 4% в области Сы <0.025 моль/л и 10.8% при максимальном значении концентрации хлоро-филлов) коэффициента обратного рассеяния, то метод малого параметра - завышенное (от 3.4% при СсЫ=0.03 моль/л до 11.5% при СсЫ=0.04 моль/л). Для коэффициента флуоресценции F метод малого параметра лишь качественно правильно описывает экспериментальные зависимости, а ошибке вычислений по сравнению с численным решением возрастает с 16% до 24.4% при увеличении концентраций хлорофилла.

а б

Рис. 1. Зависимость дифференциальных оптических коэффициентов К (а) и F (б)

от концентрации хлорофилла СсЫ: 1 - численное решение Монте-Карло; 2 - экспериментальные данные; 3 - решение уравнения переноса на основе метода малого параметра; 4 - решение уравнение переноса на основе приближения сферических гармоник

Для коэффициента К приближение сферических гармоник дает несколько заниженное значение по сравнению с экспериментальными данными. Количественное различие при этом варьируется от 3% (в областях с низкой концентрацией хлорофил-лов) до 11% (в областях с высокой концентрацией). Отличие же от результатов моделирования методом Монте-Карло не превышают 5%. Для коэффициента F максимальные различия данных приближения сферических гармоник с экспериментальными данными не превышают 12,6%, а с результатами метода Монте-Карло - 5,3%.

ЗАКЛЮЧЕНИЕ

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

БЛАГОДАРНОСТИ

Работа выполнена при поддержке Федеральной целевой программы "Научные и научно-пе-

дагогические кадры инновационной России" (2009-2013гг.), выполняемой в рамках мероприятия 1.3.1. "Проведение научных исследований

молодыми учеными - кандидатами наук"

СПИСОК ЛИТЕРАТУРЫ

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

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

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

4. Wang L, Jacques S.L., Zheng L. MCML - Monte Carlo modeling of light transport in multi-layered tissues // 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. Яровенко И.П. Численное решение краевых задач для уравнения переноса излучения в оптическом диапазоне // Вычислительные методы и программирование: новые вычислительные технологии. 2006. №1. С. 93-104.

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

9. Tuchin V.V. Tissue optics. 2000. 353 p.

10. PrahlSA. Light transport in tissue. 1988. 221p.

11. Kim A.D., Ishimary A. Optical Diffusion of Focused Beam Wave Pulses in Discrete Random Media // J. of App. Opt. A. 2000. №2. P.321-334.

12. Захаров В.П., Братченко И.А., Синдяева А.Р., Тимченко Е.В. Моделирование распределения энергии оптического излучения в растительной ткани // Оптика и спектроскопия. 2009. Т.9. №6. С. 957-962.

MODELING OF OPTICAL RADIATION INTERACTION WITH MULTIPLE SCATTERING MEDIUMS

© 2010 V. P. Zakharov, I. A. Bratchenko, E.V. Timchenko Samara State Aerospace University

Mathematical model of optical radiation interaction with multiple scattering medium is build with consideration of medium structural heterogeneity, its spectral properties and fluorescence effects. Developed model is implemented for Henyey-Greenstein phase function. The approximate solution based on expansion in perturbation series of diffuse and fluorescent fluxes of transport equation is proposed. The application of spherical harmonics approximation of first order for description of optical radiation interaction with multiple scattering mediums is investigated. The comparative analysis of approximate model shows its correspondence with physical experimental data.

Key words: mathematical model, multiple scattering medium, approximate methods.

Valery Zakharov, Doctor of Physics and Mathematics, Professor at the Automatic Systems of Power Plants Department. E-mail: Zakharov@ssau.ru Ivan Bratchenko, Graduate Student, Engineer at the Automatic Systems of Power Plants Department. E-mail: Ud_liche@mail.ru Elena Timchenko, Candidate of Physics and Mathematics, Assistant Lecturer at the Automatic Systems of Power Plants Department. E-mail: Vorobieva.82@mail.ru