Весовые оценки метода Монте-Карло для неинвазивного измерения оптических параметров биологических тканей Текст научной статьи по специальности «Математика»

АПВПМ-2019

ВЕСОВЫЕ ОЦЕНКИ МЕТОДА МОНТЕ-КАРЛО ДЛЯ НЕИНВАЗИВНОГО ИЗМЕРЕНИЯ ОПТИЧЕСКИХ ПАРАМЕТРОВ БИОЛОГИЧЕСКИХ ТКАНЕЙ

А. В, Лаппа1, А. Е, Анчугова1, Д. Ю, Шакаева1

1 Челябинский государственный университет, 454001, Челябинск

УДК 519.245

Б01: 10.24411/9999-016А-2019-10042

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

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

Введение

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

1 Методы

1.1 Общие оценки в общей модели

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

ISBN 978-5-901548-42-4

1.1.1 Уравнение и функционал

В кинетической модели световое поле в мутной среде от стационарных монохромных источников в пренебрежении неупругими взаимодействиями и поляризацией описываются скалярной лучевой интенсивностью I = I(х) с фазовыми координатами х = (г, П), описывающими положение точки в пространстве R3, направление лучей на единичной сфере О С R3. Параметрами модели являются показатель преломления п, коэффициенты поглощения и рассеяния и индикатриса (фазовая функция) рассеяния т(П ^ П') (зависит от угла рассеяния arccos(nn') ).

Объект, вообще говоря, гетерогенный, т.е. занимаемая им область G С R3 состоит из одной или нескольких попарно непересекающихся открытых зон G\, G2,Gm и кусочно-гладкой граннцы Г, состоящей из внутренней границы Г^п (между зонами) и внешней границы Гои4. Не ограничивая общности, будем полагать область G ограниченной, одпосвязпой и обобщенно выпуклой. Последнее означает, что G невогнута и что любая прямая, проходящая через точку г G G = G\ + ... + Gm, имеет конечное число общих точек (пересечений и касаний) с границей Г. Дадее, без оговорок, термин граница и символы Г Fin, Гои4 применяются только в отношении гладких (открытых) частей соответствующих границ. Таким образом, ребра исключаются, и в любой точке границы г G Г определена нормаль nr; условимся считать, что па внутренней границе она направлена в сторону зоны с меньшим номером, а на внешней — из среды.

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

Q = (I, v) + [I,b]> = i I(x)v(x)dx + i dF i dtt |Ппг| I(x)b(x) (1)

JX JT JOr>

Здесь и далее: X = {G x О}, dx = drdn, Or<(Or>) — множество направлений П таких, что Ппг = 0 при г G Fin и Ппг > 0(Ппг < 0) щи г G rout. Принято, что интенсивность на границе есть интенсивность падающего на нее излучения. Функционал (1) задает показания детекторов двух типов: объемного и граничного. Показание первого формируются светом внутри среды, второго — светом, падающим на границу и покидающим ее; v,b — соответствующие функции отклика. Предполагается, что граничные детекторы, покидающего границу излучения, сосредоточены в фазовой области, свободной от поверхностных источников (детектирующая граница не излучает).

В сопряженной форме этот функционал имеет вид:

Q = (s,J) + [a,J]< = s(x)J(x)dx + dF dn |Ппг| a(x)J(x), (2)

Jx J Г Jor<

Где j — функция ценности, имеющая смысл показания детектора, создаваемого лучом света единичной мощности, испускаемым из фазовой точки х, s — функция (объемных) источников внутри объекта (фазовая плотность мощности), а — интенсивность покидающего границу света, создаваемого поверхностными источниками (поверхностная плотность источников есть а (г, П) |Ппг|).

Прямое представление (1) и сопряженное представление (2) определяют одну величину, но условия их применимости различны. Для применимости прямого требуется регулярность источников, сопряженного — детекторов. Первое условие исключает, в частности, точечные, мононаправленные, и другие сингулярные источники, являющиеся естественными, широко используемыми моделями для множества реальных задач. Второе условие — гораздо менее ограничительно, оно отражает реальность и является естественным для алгоритмов метода Монте-Карло. Поэтому, далее мы развиваем сопряженный подход, допуская сингулярность источников, т.е. полагая функции «sb представлении (2), вообще говоря, обобщенными. Функция ценности и интегро-дифференциальное сопряженное уравнение, при этом, остаются классическими.

Интенсивность же при сингулярных источниках является обобщенной функцией и уравнение переноса необходимо рассматривать в обобщенном смысле. Разрывы коэффициентов приводят к затруднениям такого рассмотрения, для преодоления которых нужно обобщать обобщенный подход в духе теории Коломбо [6] и других, что громоздко и не сделано для уравнения переноса. Более адекватный подход требует введения мер для описания поля излучения вместо интенсивности. Не все физики готовы к такому переходу, в частности из-за потери уравнения переноса в интегро-дифференциальной форме. Желание избежать всех этих усложнений, явилось важной причиной выбора сопряженного подхода.

В работе [1] мы получили интегральное уравнение для функции ценности, корректно учитывающие внутренние и внешние граничные условия:

J = (LS + &PX)J +(Lv+@b), (3)

где сопряженные операторы пробега, рассеяния, достижения границы и отражения-преломления определяются, соответственно, выражениями:

Ь/(х) = [ * е-т(х'1}/(г + П1, П)М, Б/(х)=^е(г) [ /(г, П')т(г, П, ^ П')(!П', 0/(х) = е-т(х'1 /(гп, П),

■10 ■! О

Рх = Хгп(Р< + Р>) + ХоыР<, Р<> д(х) = |пгП'| д(г, П')рь(г, П ^ П')йП'.

«/ ( Пг

Здесь т(х, /) = /0 ¡(г + П 1')Ш' — оптическое расстояние между точками г и г + П1; ¡л = ¡а + — коэффициент взаимодействия, Iх — расстояние вдоль луча, выходящего из точки г в направлении П до первого достижения границы, г^ = г + ПIх — точка этого достижения, рь(г, П ^ П') — дифференциальный коэффициент отражения-преломления (при отрицательных (пгП)(пгП') — отражение, при положительных — преломление), Хгп(г),Хош(г) — индикаторы внутренней и внешней границ.

Это уравнение определено и справедливо внутри объекта, т.е. в X, и на границах для «пересекающих» направлений: входящих та внешней границе и всех на внутренних, т.е. на Х< = {х : г € Г, О € Ог<}. Объединение X = X У Х< совпадает с областью интегрирования функции ценности в (2), так что X является естественным фазовым пространством для уравнения (3) и задачи в целом. Естественным пространством для решений уравнения (3) является банахово пространство ограниченных (строго) в X функций. Уравнение имеет в нем единственное итерационное решение, если функции отклика V, Ь и введенные выше операторы

ограничены, а среда всюду имеет ненулевое поглощение: М |ча(г)| > 0. Для существования функционала

геа

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

1.1.2 Простые оценки сигналов

Оператору уравнения (3) и заданным источникам отвечает аналоговая («физическая») марковская цепь в X. Цепь эта однородная и обрывающаяся: п < то. Реализацию этой цепи х0, ...,хп естественно интерпретировать как последовательность фазовых координат событий на траектории частицы: старт, столкновения двух типов (рассеяние и поглощение) и достижение границы двух типов (с вылетом из среды и без тако-п

Дополним реализацию этими событиями, приписав им помер п + 1 и направления до поглощения и после вылета.

Физическая цепь полностью задается геометрией объекта, функциями источников и оптическими параметрами, т.е. набором исходных характеристик: а, а, ¡а, ¡е,т,рь. Введем другую цепь, модифицированную, отвечающую той же геометрии объекта, но другим источникам и оптическим параметрам, т.е. определенную в том же фазовом пространстве X, но обусловленную набором модифицированных характеристик: Ч,Ч,ча,ч8,п),рь. Потребуем, чтобы выполнялись условия регулярности: каждая модифицированная характеристика переходит в исходную посредством умножения на некоторую классическую функцию: ^з(г) = ^3(г)р,3(г) и т.д. Эти функции будем обозначать символами характеристик, подчеркнутых волной. Определим на реализациях модифицированной цепи 6 функционалов:

п

Фаф = 53 ® (V? + Ч^ ), а = Р = е, а, (4)

г=0

где вес Цг — е TгYlm=0 Um, ит — Хг(г)р(гт; Пт—1 ^ Пт Пт — 1 ^ П^m), т >

Дтг = ^;=0 /03 Дч(гу + П 1)й1 — разность оптических путей от старта до г-го события, вычисленных с коэф-

фициентами ч и Дч = ¡V — ¡), 1г = |гг+1 — гг|, \г и ха — индикаторы Г и С, Дт (хг, I) = /0 Дч(П + Ог 1')М'

VI = 1^(хг), VI = [1" еАт(х*>1 + ПI, П)М, Г1=ха(гг+1)eiv(г i+1, П)/ч(П+1), 0

Ч = 0Ъ(хг), Ч = хг(гг+1)<цЬ(п+1, Пг), а = еАт(х-Ч

Каждый из этих функционалов является несмещенной оценкой функции ценности, объемной и граничной частей искомого функционала (2) : J(х) = Мхфа,р, (в, ^ = Ш"ф>аф, [а, 1\< = Мьфа,/з на траекториях,

начинающейся, соответственно, из неслучайной точки х и го случайных точек в X и Х<, распределенных с модифицированными нормированными плотностями И(х) и а(х) |Ппг| ((я, 1) = [а, 1]< = 1). Начальные веса мо оценок равны 1, в(хо), а(хо), соответственно. Символы М с индексами обозначают математические ожидания на модифицированных траекториях с описанными стартами.

Оценки фаф отличаются степенью аналитического осреднения «физических» вкладов, она возрастает в порядке а = с,1, е; 0 = а, е. При этом дисперсия оценок уменьшается, но увеличивается их сложность. Поэтому, вообще говоря, нельзя заранее выбрать лучшую из них. Но, в важном частном случае, когда фазовый объем детектора (области в XX, где отличаются от нуля функции отклика) мал, и вероятность попадания в него траекторией мала, следует ожидать увеличения эффективности с увеличением степени аналитического осреднения: фе<е тогда, вероятно, будет лучшей в семействе (4). При V = 0 (имеются только граничные детекторы) семейство состоит из двух оценок: «по достижениям» границ фа (0 = а) и «по событиям» фс (0 = с). Количество вкладов в последней гораздо больше (из-за аналитического осреднения), что и обуславливает ее большую эффективность.

Еще больше можно увеличить эффективность оценок, проведя дополнительное аналитическое осреднение вкладов в духе локальных оценок [7, 8, 9].

1.1.3 "Локальные" оценки сигналов

Идея перехода к таким оценкам заключается в применении весовой оценки (4) к вычислению итерированной ценности .1 = (ЬЯ + 0РХ)удовлетворяющей уравнению (3) с другим свободным членом: (ЬЯ + 0Рх)(Ьу(х)+0Ь(х)). Его носитель шире носителя исходного свободного члена: Ьу(х)+0Ь(х), и потому оценка .1 с вкладами (ЬЯ + 0Рх)(^+0Ь) будет иметь большую вероятность ненулевых реализаций, чем оценка . с вкладами ^+06. В результате новый алгоритм может быть существенно эффективнее старого. Вместо (4) получаем семейство функционалов на траекториях модифицированной цепи:

п

£а,13 = ^3 Ф^Н + а = е,г,с, 0 = е,а, (6)

¿=о

. 1 (х) =

Мх£,аф, (в,.1) = £„,/3, [а,.1]< = Мь£,аф при тех же стартах и весах, что и в (4). Вклады Ь^ рассчитываются по формулам (5) с заменой V, Ъ на v1 = БЬу + Я0Ъ,Ъ 1 = РхЬу+Рх0Ь. Переход к искомым (не итерированным) ценностям и функционалам осуществляется с помощью связи: . = .1 + Ьу+0Ь. Заметим, что при граничных детекторах (V = 0) даже простейшая оценка из (6) имеет ненулевые вклады от событий, причем, гораздо чаще, чем любая из (4).

1.1.4 Оценки производных

Пусть г — некоторый параметр и существуют и ограничены производные по нему от всех заданных функций, определяющих систему: источники в, а, среду р, рв, и, ръ и детекторы V, Ь. Тогда, при выполнении принятых ранее допущений, обеспечивающих существование ценности и искомой характеристики, также существуют и их производные по этому параметру: .', Q', причем операцию дифференцирования можно ввести под интегралы в (2).

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

п п

ф'а,р ^[ф(Vf + Ц)' + ^ + Ц)], = + ьЫ + ч'Ж + Цг)], (7)

¿=о ¿=о

являются несмещенными оценками производных от функций ценности . и .1, соответственно: .'(х) =

Мхф'аф,л[(х) = Мх^'аф.

Оценки для ^^ ^^^^^^^^^^^ выражения (2) и применением для оценок ф'а ^ или

+ (^(хо)+0Ъ(х0)) . Если функции источников в, а те зависят от параметра г/, то получаются оценки вида (4), (6), в которых вместо величин ф, V, Ь стоят их производныеф', V', Ъ'.

1.2 Оптические параметры 1.2.1 Модель и оптические параметры

Оптические параметры биологической ткани это параметры модели, принятой для описания распространения света в ткани. В кинетической модели это ¡а, п. Ясно, что объем информации, содержащейся в измеряемых характеристиках, недостаточен для определения всех этих параметров, и потому необходимо перейти к редуцированной модели с меньшим количеством параметров. Редуцированная кинетическая модель (сформулированная в [3], фактически используемая и в других работах) получается из общей простым фиксированием индикатрисы, т.е. априорным назначением ее (индикатрисы) формы и параметров. Тогда в модели остаются 3 варьируемых кинетических параметра: п,^а,ч3, однозначно связанных с известным 3-х параметрическим набором, используемым для описания полей излучения: п, ¡а, где = ¡3(1 — д) —

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

1.2.2 Схема

Вместо пары ¡а, удобнее определять пару ¡в,р = Ча/Ч однозначно связанную с искомой: ¡а = ¡8р/(1 — р), = Ч' (1 — д). Они определяются в диффузно-отражательной схеме (рис. 1): Облучение и регистрация от-

Рис. 1: Схема измерений оптических параметров

раженного излучения осуществляется с помощью специального многоволоконного зонда. Он состоит из 1 п

чающего. По излучающему световоду подводится свет от монохромного источника к исследуемому образцу среды, по считывающим свет отводится от образца к детекторам мощности, в качестве которых мы используем фотодиоды. Сначала измеряется сигнал ¿0 от излучающего световода, затем зонд прикладывается к образцу, и измеряются сигналы ¿1,..., ¿п от считывающих световодов. Оптические параметры ¡3,р извлекаются из п те зависящих от мощности источника отношений ((г = ^/¿0, г = 1, ...,п, с использованием таблиц теоретических оценок этих относительных сигналов и их производных по этим параметрам для заданного достаточно широкого набора пар ¡3,р в описанной выше модели.

1.2.3 Оценки сигналов и производных

Оценка сигналов производилась в следующих, достаточно реалистических допущениях:

Теоретическая оценка относительных сигналов производилась в следующих, достаточно реалистических

допущениях:

А1. Среда однородна и полу бесконечна.

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

A3. Излучающий световод испускает свет равномерно по торцу его сердцевины и азимутальному углу. Угловое распределение мощности Sa(r&), измеренное в воздухе (проинтегрированное по торцу), нормированное на 1, определяет источник для оценки Qi.

А4. Мощность света, выходящего из детектирующего световода, не зависит от азимутального угла и точки падения луча па торец его сердцевины. Зависимость мощности Bai ($) от полярного угла падения луча единичной мощности, измеренная в воздухе, определяет соответствующий детектор.

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

А6. Все световоды цилиндричны с одинаковым радиусом сердцевины го-

В этих предположениях относительный сигнал Qi для заданной пары ps,p может быть представлен в виде второго (граничного) члена функционала (1), [I, bi]>, в котором граница Г = rout — плоская поверхность раздела среды и зонда, 1(г, П) = I(p, z, П) — интенсивность, рассчитанная с единичным коэффициентом

границе), имеющего угловое распределение, отвечающее «воздушному» распределению Sa(r&). Зависимость от ps и радиуса физического источника г0 перенесена па функцию отклика:

bi(p, П) = Bi(0)<pu/ro(p/(psго)), Р = |p|, cose = nrn, Bi(d) = 2n-2[Tmf (0)/Tat(6ma)]Bai(9ma),

( f min( t+l,A+l) /12 _ l+t/2 \ (\2- 1+t '2 \

<pA(t) = < arceos - arceos --- t'dt', |i - A| < 2; 0, |i - A| > 2

[Jmax(t-1,A-1) V 2t>t J V 2t'A J

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

Оценку такого функционала будем строить на модифицированной марковской цепи, отвечающей среде

с прежней (физической) индикатрисой рассеяния, единичным коэффициентом рассеяния и назначенной

р

Сначала выбирается полярный угол го «воздушного» распределения Ба, затем он пересчитывается как = 1 т. е. выполняется преломление из воздуха в среду. Азимутальный угол старта полагается пулевым. Веса в оценках (4), (6) при такой модификации достаточно просты: Як = (Т{т(-&™*)/Т{а(-&™1)) ехр {[и - V]Ьк}, где V = (1 -р)~\ V = (1 -р)-1, Ьк — путь, пройденный частицей до к-то события. Сами оценки (4), (6) в простейших вариантах принимают вид:

фа = ^2 1аЪг(Ра, 0а), Сс,а ^V-1/ ¿Пи(П'к ^ П) eXp{-lкv}Bi(в)<Plí/ro (|Гк + Пк/к| / (.Р^о)) (8)

а к -1ПпГ>°

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

Первая оценка строится по достижениям границы с углом падения внутри апертуры детектора А = тах{в : В^(в) > 0}. Вторая оценка формируется из гораздо более частых событий — столкновений, дающих больше ненулевых вкладов, так как для их расчёта не требуется попадание в угловую апертуру. Дисперсия этой оценки намного ниже, но ее сложность и время вычисления больше. Можно упростить оценку, ран-домизируя интеграл, например, вычисляя его методом Монте-Карло по одному случайному направлению внутри апертуры (0, А) х (0, 2-к).

Дифференцируя почленно оценки по параметрам, получаем оценки для производных Qi^s и Qlip. Из £с,аполучаются оценки производных по столкновениям:

^ =12 дк ^Рэ)-1 ¿Пи(П'к ^ П)ехр{- 1к v}Bi(в)фli/r0 (|гк + Пк1к |/(ра г о)) (9)

к ./Ппг>о

= Чк V2 ¡>-1 ¿Пт(П'к ^ П)(Ьк + 1к) ехр{—1к 1у}В0)рк/Го (|гк + Пк1к |/(¡3 Г0), в) (10)

к ¿Ппг>0

( т!п( 4+1,Л+1) .,(.2 ,1_ ,,2) /2 + \2 _ 1 \ 1

фл&) = { -—-)-Ш атссо8( - Л', — Л| < 2; 0, \Ъ — А| > 2} .

I «/ тах( £— 1,л-1) ((2Ы')2 — (12 — 1+1 '2)2) ' \ 2 Л / )

Здесь также для упрощения можно применить рандомизацию интегралов.

Поскольку модифицированная цепь не связана с параметрами среды и геометрией зонда, оценку всех сигналов и их производных для множества заданных наборов параметров можно выполнять в одном расчете на одном ансамбле траекторий. Функции рл(£), Фл&) имеет смысл предварительно затабулировать.

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

2 Результаты

На основе «локальных» оценок по столкновениям (8) был разработан метод получения оптических параметров биологических тканей. Определение пары происходит в два этапа: предварительного и рабочего. Первый этап включает

1. измерение угловых распределений интенсивности источника ва(^) и чувствительности детектора в воздухе Ваг (вта), характеристик зонда а, I;

2. расчет и табулирование показаний детекторов (г = 1,4 — номер детектора) и их производных Я^3, ((' Для сред с различными заданными оптическими параметрами;

На рабочем этапе выполняется:

1. освещение исследуемой среды низкоинтенсивным монохроматическим излучением и измерение сигналов ¿0, ¿г",

2. решение обратной задачи по поиску оптических параметров исследуемой среды ¡3,р и пересчет па искомые ¡а, ч'3.

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

( г = 1 , . . . , 4

их производных по оптическим параметрам , для 17 х 21 различных тар параметров ¡3,р, заданных в диапазоне [1, 500] х [0.001, 0.1]. Все 4284 расчетов были проведены на одном ансамбле из 107 траекторий модифицированной марковской цепи менее чем за несколько минут, выбрана оценка по столкновениям. Аппаратная часть комплекса включает в себя используемые для проведения эксперимента приборы и инструменты, где основные из них это оптоволоконный зонд (для измерения диффузно отраженного света), источник монохроматического света, детектор излучения, моторизированная поворотная платформа (для снятия угловых распределений).

Оптические параметры не являются строгими константами среды (зависят от принятой модели), и проверка метода состоит в оценке точности прогнозирования полей излучения по определенным параметрам. Мы провели такие испытания в эксперименте по прохождению широкого пучка света.

Излучение с помощью световода направляли на смесь чернил, интралипида и воды и регистрировали с помощью оптического зонда, погруженного в смесь и расположенного на заданном расстоянии от торца световода. Ось пучка совпадала с осью излучающего световода зонда. Мы подготовили две смеси с разной толщиной слоя жидкости (3.85 и 4 мм) и измерили их параметры методом диффузного отражения с использованием 2, 3 и 4 детекторов. На рис. 3 измеренные относительные сигналы в этих смесях сравниваются с расчетными. Расчеты проводились с использованием той же программы Монте-Карло, которая применялась для подготовительных расчетов. Мы оцениваем результаты тестирования как хорошие. Ошибка предсказания полей внутри ткани на глубинах до 4 транспортных длин не превосходит 10%, а на глубинах до 2 транспортных длин — не превосходит 5%.

Рис. 2: Сигналы четырех детекторов и их производных по оптнческнм параметрам полученные одновременно в одном расчете. Разные цвета линий соответствуют разным расстояниям между центрами излучающего и считывающих волокон. Значения расстояний указаны у кривых в мм.

Рис. 3: Измеренные сигналы от четырех детекторов, расположенных на разных расстояниях от оси пучка, в эксперименте на прохождение и их предсказания, рассчитанные из измерений по диффузному отражению (по данным 2. 3 и 4 детекторов) для двух смесей

В качестве примера реального применения метода мы изморили оптические параметры цельной неразбавленной венозной крови человека. Полученные результаты согласуются с данными других авторов (Таблица 1). В ней использованы следующие обозначения: Ш.с гоматокрит отношение объема эритроцитов к объему жидкой части крови. БаЮ2 сатурация уровень насыщения кислородом гемоглобина крови, выраженный в процентах, 7 — скорость кровотока.

Таблица 1: Оптические параметры человеческой крови

Статья Образец А, нм 9 p's, см 1 Ра, СМ 1

Наша работа Цельная кровь (Ни- 0.52). 8аЮ2 = 98%, без течения 657 0.98 11.3 12.7

Уагов^-вку и др. [9] Цельная кровь (Шс=0.45-0.46), разбавленная до Шс=0.001 с помощью фосфатного буферного раствора, 8аЮ2 > 98%, без течения 633 0.982 11.6 15.5

Метке и др. [10] (дан- 4.8x106 ГШСв/пЛ в плазме (Шс=0.4), 650 0.988 12 1

ные из графиков) 8аЮ2 >98%, 7=600 с-1

Ш^ап и др. [11] (дан- Эритроцитарный концентрат, разбав- 633 0.978 20 8

ные из графиков) ленный до 11|с 0.5 с помощью фосфатного буферного раствора, 8аЮ2 >98%, 7=500 с-1

Заключение

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

Список литературы

fl] Lappa A.V. and Anchugova А.Е., Monte Carlo weighted algorithms for calculation of radiation characteristics and their derivatives in the biomedical optics problems// Proc. SPIE 10876, 1087626, 2019.

[2] Lappa A.V., Anchugova A.E. and Shakaeva D.Yu. Extraction of tissue optical parameters from diffuse reflectance measurements with a new able to count derivatives inverse Monte Carlo method// Proc. SPIE 10876, 108768, 2019.

[3] Lappa A.V., Kulikovskiy Ar.N., Kulikovskiy An.N., Makarova T.A. A non-destructive diffuse reflectance calibration-free method for determine optical parameters of biological tissues," Proc. SPIE 8579, 857912, 2013.

[4] Лаппа А.В., Анчугова A.E. Весовые оценки метода Монте-Карло для решения прямых и обратных задач оптики мутных сред// Труды Международной конференции "Марчуковские научные чтения-2017 с. 521-527.

[5] Henyey L. G., and Greenstein J. L. Diffuse radiation in the galaxy// Astrophys. J.„V. 93, p. 70-83, 1941.

[6] Colombeau J. F. New Generalized Functions and Multiplication of the Distributions// North Holland, Amsterdam, 1984.

[7] Михайлов Г. А., Войтишек А. В. Численное статистическое моделирование. Методы Монте-Карло// Изд.: Академия, 2006.

[8] Kalos М.Н. On the estimation of the flux at a point by Monte Carlo// IEEE Transactions on Nuclear Science,V. 16, iss. 3, 1963.

[9] Yaroslavsky A.N., Yaroslavsky I.V., Goldbach T. and Schwarzmaier, H.J. The optical properties of blood in the near infrared spectral range// Proc SPIE 2678, 314-324, 1996.

[10] Meinke М. С., Miller G. J., Helfmann J. and Friebel M. Optical properties of platelets and blood plasma and their influence on the optical behavior of whole blood in the visible to near infrared wavelength range// J. Biomed. Opt., V. 12, iss. 1, 014024, 2007.

[11] Roggan A., Friebel M., Doerschel K., Hahn A. and Mueller G. J. Optical properties of circulating human blood in the wavelength range 400-2500 urn// J. Biomed. Opt., V. 4, iss. 1, 1999.

Лаппа Александр Владимирович — д.ф.-м.н., проф. кафедры теоретической физики

Челябинского государственного университета;

e-mail: lappa@csu.ru;

Анчугова Анастасия Евгеньевна — аспирант кафедры теоретической физики

Челябинского государственного университета;

e-mail: anchugova.ae@gmail.com; Шакаева Дарина Юсуповна — аспирант кафедры теоретической физики

Челябинского государственного университета;

e-mail: sh.darina@mail.ru.

Дата поступления — 30 июня 2019 г.