Математическое моделирование процессов переноса излучения в многослойных средах с подвижными рассеивателями

Размещено на http://www.allbest.ru/

Автореферат

диссертации на соискание ученой степени

кандидата физико-математических наук

Математическое моделирование процессов переноса излучения в многослойных средах с подвижными рассеивателями

Специальность 05.13.18 - Математическое моделирование,

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

Старухин Павел Юрьевич

Саратов - 2011



Работа выполнена в ГОУ ВПО «Саратовский государственный технический университет»

Научный руководитель: доктор физико-математических наук, профессор Клинаев Юрий Васильевич

Официальные оппоненты: доктор физико-математических наук, Заслуженный деятель науки РФ,

Лауреат государственной премии РФ, профессор Синицын Николай Иванович

доктор физико-математических наук, профессор Губатенко Валерий Петрович

Ведущая организация: Саратовский государственный университет им. Н.Г. Чернышевского

Защита состоится 15 июня 2011 года в 13-00 часов на заседании диссертационного совета Д 212.242.08 при ГОУ ВПО «Саратовский государственный технический университет» по адресу: 410054, г. Саратов, ул. Политехническая, д.77, корп. 1, ауд. 319.

С диссертацией можно ознакомиться в научно-технической библиотеке ГОУ ВПО «Саратовский государственный технический университет».

Автореферат размещён на сайте Саратовского государственного технического университета www.sstu.ru 13 мая 2011 г.

Автореферат разослан «13» мая 2011 г.

Ученый секретарь

диссертационного совета Терентьев А.А.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность исследований. В настоящее время одно из актуальных научных направлений связано с анализом влияния физиологических и структурных свойств биологических объектов на их физические и, в частности, оптические свойства [Приезжев А.В., Тучин В.В., Ульянов С.С.]. В этой связи важное значение приобретает математическое моделирование физических процессов, протекающих в объекте исследования в условиях лазерного облучения, при этом возникает задача решения интегродифференциального уравнения теории переноса излучения (ТПИ) (Чандрасекар, Исимару).

Одним из часто используемых методов численного решения задач ТПИ является метод статистического моделирования Монте-Карло (М. Кейзер, С. Жак, И.В. Ярославский, И.В. Меглинский, Ф. де Мул, М. Кулинк и др.).

К настоящему моменту разработаны программные средства для расчета распределения рассеянного излучения в неоднородных средах, в том числе с учетом локальных неоднородностей, например опухолей. Существует коммерческий программный комплекс TracePro фирмы Lambda Research Corporation для расчета оптических систем, позволяющий проводить трассировку лучей методом Монте-Карло. Изучены особенности динамики импульсов рассеянного излучения (A. Popov, A. Priezzhev, and R. Myllylд) и возможности диагностики скорости потоков жидкостей доплеровским методом (A.V. Bykov, M.Yu. Kirillin, A.V. Priezzhev, R.А. Myllylд).

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

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

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

Для достижения поставленной цели сформулированы следующие задачи: фотон излучение рассеиватель

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

· Разработка численных алгоритмов решения стационарного и нестационарного уравнений переноса излучения в непрозрачных средах на основе метода моделирования Монте-Карло, позволяющего учитывать чрезвычайно большие объёмы (~10¹⁰) рассеиваемых пакетов фотонов.

· Реализация разработанных алгоритмов моделирования в виде комплекса программ для многопроцессорных вычислительных систем.

· Анализ возможности применения современных многопоточных центральных процессоров и многопоточных графических процессоров для производства параллельных вычислений по методу Монте-Карло.

Научная новизна работы:

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

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

· Предложены эффективные подходы к реализации моделирования переноса излучения с использованием многопоточных графических процессоров, позволившие сократить время моделирования на 2 порядка.

· Разработан проблемно-ориентированный комплекс программ для моделирования распространения излучения в рассеивающих средах.

· Проведен ряд вычислительных экспериментов, в результате которых установлены следующие функциональные зависимости:

- при непрерывном характере облучения соотношение интенсивности доплеровской фракции излучения к полной интенсивности рассеянного излучения квадратично растет с увеличением объемной концентрации подвижных рассеивателей в поверхностных слоях среды;

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

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

Научная и практическая значимость.

Научная значимость работы заключается в применимости разработанного на основе методов Монте-Карло проблемно-ориентированного комплекса программ моделирования распространения излучения в мутных средах для исследования распределения интенсивности рассеянного излучения и уширения спектра рассеянного излучения вследствие эффекта Доплера.

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

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

Положения и результаты, выносимые на защиту:

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

2. Программный комплекс «Mcml-doppler» позволяет решать задачи стационарной и нестационарной ТПИ, рассчитывать пространственно-временное распределение интенсивности рассеянного излучения и характеристики доплеровского уширения спектра излучения рассеянного подвижными частицами.

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

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

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

6. Разработанный алгоритм моделирования переноса излучения позволяет распараллелить вычислительный процесс в требуемом масштабе для эффективного использования многоядерных/многопоточных CPU, и многократно сократить время вычислительного эксперимента при использовании современных графических процессоров, либо увеличить точность вычислений при заданном времени, затрачиваемом на процесс моделирования.

Апробация работы. Основные результаты диссертационной работы докладывались на российских и международных конференциях: Молодежная научная школа по оптике, лазерной физике и оптоэлектронике «Проблемы оптической физики» (Саратов, Россия, 1997), Международная научно-техническая конференция «Актуальные проблемы электронного приборостроения» (Саратов, Россия, 2010), Всероссийская научно-практическая конференция молодых ученых «Инновации и актуальные проблемы техники и технологий 2010» (Саратов, Россия, 2010), Ежегодная Всероссийская научная школа-семинар «Методы компьютерной диагностики в биологии и медицине 2010» (Саратов, Россия, 2010), XX Ежегодная международная конференция-выставка «Информационные технологии в образовании «ИТО - 2010» (Москва, Россия, 2010).

Публикации. По теме диссертации опубликовано 20 печатных работ, в том числе пять статей - в изданиях, рекомендованных ВАК РФ. Программный комплекс моделирования зарегистрирован ИНИМ РАО РФ / Объединённый фонд электронных ресурсов «Наука и образование» 02.12.2010г. Список 13 публикаций приведен в конце автореферата.

Структура и объем работы. Диссертационная работа состоит из введения, 4 глав, заключения, списка цитируемой литературы из 128 наименований и 4 приложений. Диссертация содержит 7 таблиц, иллюстрирована 43 рисунками. Общий объем диссертационной работы составляет 162 страницы.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

Во введении дана общая характеристика работы: обоснована актуальность, отмечена научная новизна и практическая значимость работы. Формулируются цели и задачи исследования. Проанализированы и обоснованы применяемые в работе методы исследования.

Рис. 1. Схема многослойной модели кожи и моделируемого эксперимента

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

мa - коэффициент поглощения (величина, обратная расстоянию, на котором интенсивность уменьшается за счет поглощения в е раз) [1/см];

мs - коэффициент рассеяния [1/см];

g - фактор анизотропии или средний косинус угла рассеяния;

n - показатель преломления;

d - толщина слоя [см]). (1)

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

где - лучевая интенсивность в точке в направлении , в момент времени t[Вт·м-²·ср-¹] (ср - стерадиан), t₂=(µ_t c)-¹ - среднее время между взаимодействиями, µ_t=µ_s+µ_a - полный коэффициент ослабления или экстинкции (в биологических средах µ_s>>µ_a), описывает временную деформацию д-образного импульса после акта рассеяния, - фазовая функция рассеяния, - единичный телесный угол в направлении , описывает функцию источника излучения в среде, c - скорость света в среде. Уравнение (2) дополняется начальными условиями:



Исключение из (2) переменной t, приводит к уравнению стационарной ТПИ.

Если процесс переноса излучения исследуется в области , а ?G - граница области, то граничные условия на ?G могут быть записаны в виде

где , n - внешняя нормаль к ?G, - лучевая интенсивность падающего излучения, - оператор отражения. При наличии в области G отражающих и преломляющих поверхностей условия, аналогичные (3') должны быть записаны для каждой из них. Данное граничное условие отражает тот факт, что входящее в G излучение складывается из излучения внешних источников и той части вышедшего из G излучения, которое снова попало в G за счет отражения, описываемого оператором (например, Френелевское отражение).

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

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

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

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

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

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

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

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

Тогда на каждом шаге алгоритма моделирования выполняются действия:

1. Формирование очередного пакета фотонов на границе среды.

Плоский возбуждающий импульс длительности T, соответствующий члену уравнения (2), моделируется путем облучения поверхности среды N пакетами фотонов, для каждого из которых время влета t₀ в среду равномерно распределено в диапазоне от 0 до T и определяется

Здесь и далее - случайное число, равномерно распределенное на отрезке [0,1].

Пространственное и угловое распределение пакетов фотонов в импульсе соответствует профилю падающего излучения. Например, для точечного источника, излучающего вдоль оси Z по нормали к поверхности среды, начальные координаты нового пакета фотонов:

.	(5)

Исходные направляющие косинусы волнового вектора пакета фотонов:

; .	(6)

Безразмерный параметр, отражающий энергию поступающего в неоднородную среду пакета фотонов, называется весом пакета и имеет начальное значение:

.	(7)

2. Расчет траектории.

Поглощение излучения в однородном слое среды описывается законом Ламберта-Бугера-Бера. Ослабление интенсивности за счет рассеяния также носит экспоненциальный характер. Поэтому ослабление интенсивности излучения в заданном направлении с учетом поглощения и рассеяния будет определяться коэффициентом полного ослабления (экстинкции) и при распространении излучения в направлении оси z можно записать

.	(8)



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

,	(9)

что отражает ослабление интенсивности излучения в соответствии с первым членом из правой части (2).

Дальнейшее направление движения будет определяться фазовой функцией рассеяния из второго члена правой части (2).

При переходе к полярным координатам направление определяется двумя углами рассеяния - полярным и и азимутальным ц. Тогда можно записать

=p(и)p(ц).	(10)

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

.	(11)

Расчет полярного угла и в соответствии с (11), осуществляется по формуле:



.	(12)

В случае, если рассеивающая среда характеризуется изотропным рассеянием, т.е. g=0, угол и можно определить на основе выражения

.	(13)

Азимутальный угол рассеяния ц в однородном слое среды считается произвольным и равновероятным в диапазоне от 0 до 2р:

.	(14)

На основе (9)-(14) формируется модель учета многократного рассеянного излучения в соответствии со вторым членом из правой части (2).

Направляющие косинусы , , единичного вектора скорости пакета фотонов k_s после рассеяния могут быть определены из соотношений:

,	(15)

где , , - направляющие косинусы единичного вектора скорости пакета фотонов k_i до рассеяния.

Координаты пакета после рассеяния определяются по формулам:



.	(16)

Время, затрачиваемое пакетом фотонов на перемещение между последовательными рассеяниями, учитывается следующим образом:

,	(17)

где l_св - длина пробега между рассеяниями, n - показатель преломления среды.

Учет доплеровского уширения спектра рассеянного излучения

При наличии в среде подвижных частиц и в случае однократного рассеяния фотона (пакета фотонов) на частице, движущейся со скоростью v, возникает эффект Доплера, учитывать который можно в форме сдвига частоты излучения следующим образом:

,	(18)

где - сдвиг частоты излучения, k_i, k_s - единичные векторы скорости фотона «до» и «после» рассеяния, v - вектор скорости движущейся частицы.

Если в процессе распространения фотон испытывает n рассеяний движущимися частицами, его суммарный сдвиг частоты вычисляется как: . В результате расчета распространения в среде заданного числа пакетов фотонов (10⁶ - 10¹⁰), детектора заданной геометрии достигнут m пакетов фотонов с общим весом: . Из них испытавших доплеровский сдвиг частоты. Тогда средний доплеровский сдвиг частоты компоненты излучения, рассеянного движущимися частицами:

.

Рис. 2. Эффект Доплера при рассеянии излучения подвижным рассеивателем

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

,	(19)

где V_подв - объем, занимаемый подвижными частицами в слое среды, V_среды - общий объем частиц.

Пространственное распределение подвижных рассеивателей в каждом слое среды считалось случайным и равномерным. Поэтому, если в процессе моделирования случайная величина о, равномерно распределенная в интервале [0,1], о < C_слоя, пакет фотонов считался рассеянным движущейся частицей и приобретал доплеровский сдвиг частоты в соответствии с (18).

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

Интенсивности отраженной и преломленной компонент излучения на границе показателя преломления рассчитываются на основе закона преломления и формул Френеля:

, ,	(20)

где - показатели преломления среды, из которой поступает излучение, и среды в которую проникает излучение, соответственно, - углы падения и преломления, - коэффициент внутреннего отражения.

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

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

3. Поглощение.

На входе в непрозрачную среду вес каждого пакета фотонов задается в соответствии с (7). За счет поглощения и рассеяния после каждого i-го акта взаимодействия со средой вес пакета уменьшается в соответствии с поглощающими и рассеивающими свойствами среды по формуле

.	(21)

Поглощенная доля энергии пакета, равная , учитывается при расчете пространственного распределения поглощенного излучения. Таким образом, рассчитывается ослабление потока излучения за счет поглощения и рассеяния, что соответствует первому члену из правой части уравнения (2).

4. Прекращение отслеживания траектории пакета.

Если после i-го рассеяния вес пакета фотонов становится меньше некоторого минимума, считается, что он полностью поглощен средой и осуществляется возврат к п. 1. алгоритма, иначе производится переход к п. 2.

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

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

После повторения шагов 1-5 для достаточного количества N траекторий пакетов фотонов карта распределения энергии отраженного, пропущенного и поглощенного излучения рассчитывается и накапливается в компьютере.

С учетом нормировки на суммарную энергию исходного возбуждающего импульса NW₀ можно записать



,	(22)

где A_m -суммарная доля поглощенной средой энергии из всех N построенных траекторий пакетов, T_d и R_d - суммарные доли диффузно-отраженной и диффузно-пропущенной энергии соответственно.

Выражение (22) отражает суть закона сохранения энергии. В результате проведения ряда статистических испытаний, формируется пространственно-временное распределение интенсивности и среднего значения доплеровского сдвига частоты рассеянного излучения.

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

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

С прикладной точки зрения особый интерес представляют методы восстановления физических свойств среды по результатам измерения характеристик рассеянного излучения. Для нахождения решения обратных задач ТПИ использовались результаты решения прямой задачи. В общем случае постановка обратной задачи ТПИ является некорректной, поскольку приводит к целому классу решений. Однако для узкопоставленной задачи со всеми известными, кроме одной переменной, можно свести задачу к однопараметрической. Тогда метод решения обратной задачи, например определения объемной концентрации подвижных рассеивателей C в среде по результатам измерения соотношения интенсивности доплеровской фракции рассеянного излучения I_D к полной интенсивности I излучения достигшего детектора, можно представить в виде последовательности шагов:

1. Определение основных характеристик объекта согласно (1);

2. Численный поиск зависимости соотношения I_D/I от концентрации подвижных рассеивателей C методом решения последовательности прямых задач ТПИ в виде I_D/I =f(C) в заданном интервале C;

3. Получение обратной зависимости C =f ^-1(I_D/I) в аналитической или численной форме на основе аппроксимации функции f;

4. Использование f ^-1 для определения концентрации подвижных рассеивателей C по экспериментально измеряемому значению соотношения I_D/I.

Во второй главе на основе разработанного алгоритма для решения прямых задач стационарной и нестационарной ТПИ реализован проблемно-ориентированный комплекс программ Mcml-doppler.

В отличие от известных программ mcml и gcmcml, предназначенных для расчета распределения интенсивности рассеянного излучения, предложенный комплекс позволяет: учитывать подвижность рассеивателей, рассчитывать доплеровское уширение спектра рассеянного излучения, использовать генераторы случайных чисел с большим периодом (~10⁹ чисел), распараллеливать процесс численного моделирования с применением многопоточных центральных и графических процессоров, управлять временем решения на основе выбора технических средств вычислительного комплекса.

Обобщённая архитектура программного комплекса представлена на рис. 3.

В главе 3 для иллюстрации возможностей разработанного программного комплекса рассмотрен ряд моделей рассеивающих сред.

Рис. 3. Схема программного комплекса Mcml-doppler

3.1. В целях верификации программного комплекса решена прямая задача стационарной ТПИ по расчету коэффициентов отражения R_d и пропускания T_d однослойной рассеивающей среды со следующими параметрами: µ_a=10см-¹, µ_s=90 см-¹, g=0.75, n=1, толщина среды d=0.02 см, показатели преломления для окружающей среды n_вн=1. Сравнение полученных результатов с данными более ранних исследований представлено в табл. 1. В столбцах R_d, среднее и Ошибка R_d табл. 1 указаны среднее значение и стандартная ошибка для коэффициента отражения соответственно, в столбцах T_d, среднее и Ошибка T_d для коэффициента пропускания соответственно.

Из табл. 1 видно, что результаты, полученные с использованием комплекса Mcml-doppler, незначительно отличаются от результатов, полученных другими исследователями, и их можно считать достоверными.

Таблица 1. Результаты верификации комплекса Mcml-doppler

Источник	Метод расчета	Rd, средн.	Ошибка Rd	Td, средн.	Ошибка Td
де Хюльст, 1980	Теория Ми	0.09739		0.66096
Пралл, 1989	Монте-Карло	0.09711	0.00033	0.66159	0.00049
Mcml-Doppler	Монте-Карло	0.09750	0.00163	0.66091	0.00041

3.2. Решена задача по моделированию распространения непрерывного излучения в двухслойной среде с расположенным во втором слое сосудом цилиндрической формы с подвижными рассеивателями. Цилиндр диаметром 0.06-0.10 см располагался на глубине 0.05-0.25 см параллельно внешней поверхности среды, ось пучка падающего излучения проходила через центр цилиндра. Значения параметров слоев среды выбраны близкими к значениям характерным для нормальной кожи человека на длине волны 633 нм.: µ_a1=0.8 см-¹, µ_s1=300 см-¹, g₁=0.8, n₁=1.4, d₁=0.02 см, С₁=0.0, µ_a2=1.0 см-¹, µ_s2=200 см-¹, g₂=0.9, n₂=1.3, d₂=1.0 см. Параметры для подвижных частиц выбраны близкими к значениям обычным для крови человека: µ_a=2.0см-¹, µ_s=200 см-¹, g=0.95, n=1.3. В рамках решения прямой задачи стационарной ТПИ получено, что контур равной интенсивности доплеровской компоненты диффузно отраженного излучения имеет форму эллипса, параметры которого описываются значениями полуосей (a - большая полуось, b - малая полуось). Коэффициент сжатия эллипса (эллиптичность) определяется как:. Для всех рассмотренных значений диаметра с увеличением глубины залегания цилиндра эллиптичность контура равной интенсивности снижается в зависимости, близкой к линейной. С точки зрения постановки обратной задачи ТПИ, измерение может быть положено в основу методики томографии кровеносных сосудов.

3.3. Численно исследованы возможности диагностики состояния кровеносной системы человека на примере семислойной модели кожного покрова, все слои которой (кроме двух верхних, моделирующих роговой слой и эпидермис) содержали сеть микроциркуляции крови. Значения параметров слоев среды выбраны близкими к значениям, характерным для нормальной кожи человека на длине волны 633 нм.

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

Нормальные значения объемной концентрации крови в каждом отдельном слое кожного покрова составляют от 4% для верхней дермы с капиллярными петлями до 14% в нижнем слое дермы с глубинными сплетениями сосудов кровеносной системы. При таких значениях концентрации крови можно говорить о нормальном состоянии объемного кровенаполнения микроциркуляторного русла V_кр. При изменениях состояния кровенаполнения, соответственно будут изменяться и значения объемной концентрации крови в каждом слое.

Задача состоит в том, чтобы связать измеряемые характеристики рассеянного излучения со значением V_кр для решения задач диагностики кровеносной системы.

Использовав нормальные значения объемной концентрации крови в каждом слое как базовые (V_кр=1.0), численно исследовано влияние объемного кровенаполнения микроциркуляторного русла кожи человека V_кр на соотношение интенсивности доплеровской компоненты к полной интенсивности излучения достигшего детектора в диапазоне V_кр от 0.5 до 1.5.

На рис. 4 приведено построение изолиний соотношения интенсивности доплеровской компоненты к полной интенсивности излучения, попавшего на детектор, при различных значениях V_кр и расстояния между источником и детектором излучения.

С использованием методики решения обратной задачи, описанной в гл. 1, и полученного решения прямой задачи ТПИ найдена зависимость объемного кровенаполнения микроциркуляторного русла V_кр от измеряемого соотношения интенсивностей доплеровской и полной интенсивности рассеянного излучения:

,	(23)

где I_D и I - интенсивность доплеровской компоненты и полная интенсивность рассеянного излучения, соответственно. Результат решения обратной задачи для одного значения расстояния источник-детектор излучения приведены на рис. 5. Полученный результат может быть полезен при разработке методов диагностики биологических тканей.

3.4. Возможности разработанного комплекса по решению прямых задач нестационарной ТПИ изучены на примере задачи о прохождении короткого импульса излучения в многослойной среде, содержащей подвижные рассеиватели. В качестве модели выбрана двухслойная среда, со статическим верхним слоем и вторым слоем с сетью микроциркуляции. Моделировалось распространение в среде прямоугольных импульсов монохроматического излучения, длительностью от 0.1 до 10 пс.

Рис. 4. Изолинии соотношения интенсивности доплеровской компоненты I_D к суммарной интенсивности I рассеянного излучения при различных V_кр и расстояния источник-детектор

Рис.5. Зависимость V_кр от измеряемого отношения интенсивности доплеровской компоненты к суммарной интенсивности рассеянного излучения

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

Рис. 6. Импульс длительностью 5пс и его доплеровская компонента

Рис. 7. Влияние толщины статического слоя на задержку доплеровской компоненты импульса

3.4а) Задержка между появлением несмещенной и доплеровской компонентами сигнала возникает вследствие прохождения излучением толщи первого статического слоя среды без подвижных рассеивателей. В результате исследования для модели двухслойной среды выявлена квадратичная зависимость временной задержки доплеровской компоненты рассеянного излучения от толщины статического слоя среды (рис. 7).

3.4.б) Верхние слои кожного покрова содержат сеть микроциркуляции биожидкости, в то время как нижележащие слои содержат крупные сосуды с более скоростными частицами.

Рис. 8. Влияние строения среды на поведение доплеровского уширения спектра рассеянного сигнала

Тогда более реалистичная математическая модель биоткани будет представлена тремя слоями: верхний слой статический толщиной 0,02 см, средний с 10% содержанием частиц со скоростью 0,1 мм/c, и нижний толщиной 0,5 см с 10% содержанием частиц со скоростью 1,0 мм/c. На рис. 8 представлены результаты моделирования для трех случаев: в первом рассматривалась 2-слойная модель со вторым слоем толщиной 5мм и абсолютным значением скорости рассеивателей 0,1 мм/с, во 2 и 3 случаях трехслойная среда, для которой толщина 2-го слоя (со скоростью рассеивателей 0,1 мм/c) составляла 0,09 и 0,11 см соответственно, толщина 3-го слоя со скоростью рассеивателей 1,0 мм/c 5 мм. При исследовании зависимости среднего доплеровского сдвига частоты рассеянного излучения от времени наблюдения обнаружен излом данной зависимости (рис. 8), связанный с проникновением излучения в слой с высокоскоростными рассеивателями. Положение излома (обозначено окружностями на рис. 8) определяется толщиной второго слоя с низкоскоростными частицами. Область I на рис. 8 соответствует толщине второго слоя 0,09 см, область II - 0,11см. В диапазоне толщин второго слоя от 0,08 до 0,15 см, (соответствующем поверхностным слоям кожного покрова) зависимость положения излома от толщины второго слоя носит характер, близкий к линейному.

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

В главе 4 проведено исследование производительности многопроцессорных вычислительных систем в задачах моделирования переноса излучения методом Монте-Карло. Рассмотрены как многоядерные/многопоточные центральные процессоры общего назначения - CPU, так и многоядерные графические процессоры (GPU - Graphics Processing Units).

Алгоритмы реализации метода Монте-Карло позволяют распараллелить процесс поиска решения прямой задачи ТПИ на этапе моделирования случайного блуждания пакетов фотонов в среде. При проведении статистического моделирования с использованием многоядерных CPU общего назначения сокращение времени вычислений составило до 40-45% на каждое дополнительное процессорное ядро и до 20% на дополнительный поток по технологии Intel Hyper-Threading. Современные графические процессоры Nvidia позволяют проводить массивно-параллельные математические вычисления для решения неграфических задач с использованием интерфейса прикладного программирования CUDA (Compute Unified Device Architecture). Для оценки эффективности использования технологии CUDA проведена серия численных экспериментов по моделированию распространения излучения в трехслойной модельной среде, основные оптические параметры которой близки к характеристикам биологических сред. В задаче использовались программные комплексы Mcml-Doppler, MCML (http://omlc.ogi.edu/software/mc/) а также GCMCML (http://bmp.hust.edu.cn/ GPU_Cluster/GPU_Cluster_MCML.HTM). Количество моделируемых траекторий фотонов составило 3·10⁷. Время моделирования заданного количества траекторий с использованием CPU Intel Core i3 330M (2133 МГц) составила 1700с. Характеристики использованных GPU и результаты измерения времени моделирования представлены в табл. 2. Переход к использованию GPU сократил время моделирования на 2 порядка (в ~140 раз) по сравнению с CPU общего назначения, для поставленной задачи. Учитывая сильную взаимосвязь между производительностью и двумя параметрами - шириной шины данных и числом потоковых процессоров GPU, предложен оценочный параметр производительности:

,	(24)

где N_p - число потоковых процессоров, Q_d - ширина шины данных в битах.

Полученная зависимость времени моделирования от характеристического параметра p аппроксимируется экспоненциальной функцией:

.	(25)

Присутствие в параметре p двух характеристик, влияющих на производительность GPU, позволяет управлять временем моделирования, подбирая устройства с различным числом потоковых процессоров и шириной шины данных (24), обеспечивающих либо значение параметра p по требуемому времени моделирования (25), либо сохраняя значение р как инвариант.



Таблица 2. Характеристики использованных GPU и измеренное время моделирования тестовой задачи

	Тип GPU	Потоковых процессоров, шт.	Шина данных, бит.	Параметр p	Время расчета, сек.
1	Nvidia GT210	16	64	1024	63,28
2	Nvidia GT9500	32	128	4096	28,05
3	Nvidia GT240M	48	128	6144	22,11
4	Nvidia GT240	96	128	12288	20,86
5	Nvidia GT9800	112	256	28672	12,12

При реализации методов Монте-Карло для моделирования рассеяния излучения с применением CPU общего назначения, как правило, используется представление чисел с плавающей точкой с двойной точностью, в соответствии со стандартом IEEE 754. Архитектура современных GPU позволяет эффективно выполнять обработку только чисел с одинарной точностью. Исследование сходимости результатов моделирования с использованием GPU показало, что для достижения дисперсии вычисляемых величин (например, коэффициентов отражения R_d и пропускания T_d) на уровне результатов получаемых при моделировании на CPU с двойной точностью, необходимо увеличение числа экспериментов в методе Монте-Карло приблизительно в 4 раза.

В целом при использовании современных GPU время моделирования методами Монте-Карло сокращается многократно.

Приложения 1-4 содержат: техническое описание комплекса Mcml-doppler, программный код для моделирования переноса излучения (при стационарном и нестационарном режимах облучения), свидетельство о регистрации программного комплекса, перечень из 20 публикаций автора по теме работы.



ЗАКЛЮЧЕНИЕ

1. На основе метода Монте-Карло разработаны модели рассеивающих сред и алгоритмы решения прямых задач стационарной и нестационарной ТПИ для сложных, оптически неоднородных объектов.

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

3. Построенная математическая модель позволяет распараллелить процесс вычислений и при использовании графических процессоров, с поддержкой технологии Nvidia CUDA, многократно сократить время моделирования методами Монте-Карло.

4. Разработан комплекс программ Монте-Карло моделирования для использования на компьютерах с ОС Microsoft Windows и Linux.

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

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



Список основных публикаций

1. публикации в изданиях, рекомендованных ВАК РФ:

1. Старухин П.Ю. Моделирование доплеровского уширения спектра рассеянного лазерного излучения и диагностика кровотока в биологических тканях / П.Ю. Старухин, Ю.В. Клинаев // Вестник Саратовского государственного технического университета. 2010. № 4 (49). C. 28-35.

2. Старухин П.Ю. Математическое моделирование уширения спектров рассеянного лазерного излучения при распространении в биоткани методом Монте-Карло / П.Ю. Старухин, Ю.В. Клинаев, Д.В. Терин // Вестник Саратовского государственного технического университета. 2010. №3(48). С. 24-32.

3. Старухин П.Ю. Эффективность использования параллельных вычислений по технологии CUDA в решении задач биомедицинской физики / П.Ю. Старухин, Ю.В. Клинаев // Вестник Саратовского государственного технического университета. 2010. №3(48). С. 126-130.

4. Старухин П.Ю. Применение метода Монте-Карло для моделирования прохождения сверхкоротких лазерных импульсов через неоднородную среду с подвижными рассеивателями / П.Ю. Старухин, Ю.В. Клинаев // Журнал прикладной спектроскопии. 2011. Т. 78. № 2.С. 277-281.

5. Starukhin P. Blood-flow measurements with a small number of scattering events / P. Starukhin, S. Ulyanov, E. Galanzha, V. Tuchin // Applied Optics. 2000. Vol.39. № 16. P.2823-2830.

2. Другие публикации:

6. Старухин П.Ю. Расчет доплеровского сдвига методом Монте-Карло / П.Ю. Старухин // Проблемы оптической физики - 1997 : материалы молодежной научной школы по оптике, лазерной физике и оптоэлектронике, Саратов. ноябрь 1997 г. Саратов: Изд-во Сарат. ун-та, 1997. С. 170-171. ISBN 5-292-01579-2

7. Старухин П.Ю. Исследование влияния геометрии детекторов излучения в задачах лазерного доплеровского измерения скорости потока крови в биоткани методом математического моделирования Монте-Карло / П.Ю. Старухин, Ю.В. Клинаев // Актуальные проблемы электронного приборостроения АПЭП-2010 : материалы Междунар. науч.-техн. конф., Саратов, 22-23 сентября 2010 г. Саратов, 2010. C.163-169.

8. Старухин П.Ю. Моделирование применения эффекта Доплера для диагностики параметров микроциркуляции крови в биотканях / П.Ю. Старухин, Ю.В. Клинаев // Методы компьютерной диагностики в биологии и медицине - 2010 : материалы ежегодной Всерос. школы-семинара, Саратов, 9-11 ноября 2010 г. Саратов: Изд-во Сарат. ун-та, 2010. C. 27-35.

9. Старухин П.Ю. О сходимости результатов моделирования методами Монте-Карло при вычислениях по технологии CUDA / П.Ю. Старухин, Ю.В. Клинаев, Д.В. Терин // Инновации и актуальные проблемы техники и технологий.: материалы Всероссийской научно-практической конференции молодых ученых. Саратов, 26-29 октября 2010 г. Саратов: Сарат. гос. техн. ун-т, 2010. Т. 1. C.19-22.

10. Старухин П.Ю. Эффективность использования параллельных вычислений по технологии CUDA в задачах моделирования методами Монте-Карло / П.Ю. Старухин, Ю.В. Клинаев, Д.В. Терин // Информационные технологии в образовании. ХХ Международная конференция-выставка : сб. трудов. Ч.III. М. : МИФИ, 2010. С.68-70. ISBN 978-5-7262-1363-7

11. Старухин П.Ю., Клинаев Ю.В. // Программа моделирования процесса распространения лазерного излучения в мутных средах Mcml-doppler // ХРОНИКИ ОБЪЕДИНЕННОГО ФОНДА ЭЛЕКТРОННЫХ РЕСУРСОВ «НАУКА И ОБРАЗОВАНИЕ». №12 (19), декабрь 2010. С. 1.

12. Старухин П.Ю. Моделирование методом Монте-Карло задач теории переноса излучения через биологические многослойные среды с подвижными рассеивателями / П.Ю Старухин., Ю.В. Клинаев, О.А. Монахова // Исследования в области естественных наук и методики их преподавания: сб. науч. статей. Саратов : ООО Изд. Центр «Наука», 2011. С. 102-111. ISBN 978-5-9999-0773-8

3. Патентные документы:

13. Свидетельство о регистрации электронного ресурса № 16458 ИНИМ РАО РФ. Объединённый фонд электронных ресурсов «Наука и образование»: Программа моделирования процесса распространения лазерного излучения в мутных средах Mcml-doppler / П.Ю. Старухин, Ю.В. Клинаев. Инв. № 50201050204, дата регистрации 02.12.2010.

Размещено на Allbest.ru

...