Основным инструментом описания процесса распространения оптического излучения через сильно рассеивающие среды (СРС), такие как биологические ткани, является уравнение переноса излучения (УПИ). Оно представляет собой интегро-дифференциальное уравнение баланса энергии в среде, содержащей частицы [7]. Главная проблема заключается в том, что УПИ в общем виде не имеет аналитического решения, так как рассеяние фотонов является случайным. По этой причине для описания распространения оптического излучения в СРС с коэффициентами поглощения и рассеяния в ближнем инфракрасном диапазоне как у биологических тканей используются различные приближённые модели [1?4]. Чаще всего используются многопотоковые модели (метод Кубелки ? Мунка), моделирование методом статистических испытаний (Монте-Карло) и диффузионное приближение к УПИ.

Главным допущением метода Кубелки ?Мунка является то, что лучевая интенсивность считается диффузной. Внутри ткани (при одномерной геометрии) поток фотонов разделяется на два: в направлении падающего излучения и в обратном направлении (рассеяние назад). Лучевая интенсивность в каждом направлении испытывает два акта снижения (из-за поглощения и рассеяния) и один акт усиления (из-за рассеяния фотонов с противоположного направления). Основным недостатком метода является адекватность только в случаях, когда рассеяние многократно превышает поглощение [2, 10].

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

Главным допущением диффузионного приближения является то, что ключевым фактором ослабления света признаётся процесс рассеяния (диффузии фотонов). УПИ упрощается посредством разложения в ряд сферических гармоник. Результатом упрощения является система связанных дифференциальных уравнений в частных производных, которая может быть сведена к одному дифференциальному уравнению. Основным достоинством диффузионного приближения является то, что оно хорошо описывает распространение излучения в толще биологической ткани, а недостатком – то, что оно справедливо не во всех случаях, а лишь при больших альбедо и небольших значениях фактора анизотропии рассеяния [8].

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

Материалы и методы исследования

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

Для моделирования миграции фотонов в исследуемом объекте используется диффузионное приближение к УПИ следующего вида [2, 6]:

где  – скорость света в среде; c0  – ско рость света в вакууме; ?object – относительный коэффициент преломления моделируемого объекта (?) и его границы (??); x, y, z – координаты всех точек конечной моделируемой области; и ?a(x, y, z) – коэффициент диффузии и коэффициент поглощения в точках с координатами x, y, z; ?s(x, y, z) – коэффициент рассеяния в точках с координатами x, y, z; g – параметр анизотропии; ?(x, y, z, t) – фотонная плотность в точке с координатами x, y, z в момент времени t;  – функция источника фотонов; ? – дельта-функция,  – средняя длина рассеяния, т.е. глубина на которой возникает точечный виртуальный изотропный источник.

Для описания распространения фотонов на границе [8] моделируемого объекта ? используется граничное условие третьего типа (Робина):

(2)

где n(x, y, z) – направление внешней нормали к границе ?? в точке с координатами x, y, z. F – коэффициент френелевского отражения [2], вычисляемый как

где R0 и Qc – коэффициенты, соответственно равные

и  

где ?medium – относительный коэффициент преломления для окружающей объект среды.

Численное решение уравнения (1) с граничным условием (2) было выполнено по семиточечному шаблону. Начальное приближение функции ?(x, y, z, t) во всех узлах сетки генерируется с учетом геометрии объекта, позиции источника фотонов и количества инжектируемых в исследуемый объект в течение одиночного импульса фотонов. Количество фотонов рассчитывается на основе средней мощности используемого фемтосекундного импульсного лазера, а также его длины волны и длительности одиночного импульса. Критерием окончания итерационного процесса служит истечение заданного времени.

После завершения итерационного процесса для получения НМФП функция ?(x, y, z, t) нормируется [9] относительно своего максимума ?max(x, y, z, t):

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

где P – экспериментально найденный минимальный уровень фотонной плотности НМПФ, 0 < P ? 1.

Результаты исследования и их обсуждение

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

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

В результате компьютерного моделирования (для P = 0,997), проведённого с помощью диффузионной модели миграции фотонов в биологических тканях, было установлено, что все выявленные ранее закономерности миграции НМФП в СРС с цилиндрической геометрией [7, 9] справедливы и для СРС со сложной структурой, с той лишь разницей, что в этом случае их следует использовать в комплексе.

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

Таким образом, ключевое влияние на движение НМФП в конкретном маммографическом объекте оказывают железистая ткань, как самая крупная по площади и в тоже время наиболее поглощающая структура, и карцинома протока как крупная по площади и в то же время наиболее рассеивающая структура.

Заключение и выводы

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