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

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

Однородный случай ^[2]

Поглощающая неоднородность ^[2]

Рассеивающая неоднородность ^[2]

Определения [ править ]

Рисунок 1: Схема потока энергии через элемент дифференциальной площади в положении внутри элемента дифференциального телесного угла .${\ displaystyle dA}$${\ displaystyle {\ vec {r}}}$${\ displaystyle d \ Omega}$

RTE может математически моделировать передачу энергии при движении фотонов внутри ткани. Поток энергии излучения через элемент малой площади в поле излучения можно характеризовать яркостью . Сияние определяется как поток энергии на единицу нормальной площади на единицу телесного угла в единицу времени. Здесь обозначает положение, обозначает единичный вектор направления и обозначает время (рисунок 1). Несколько других важных физических величин основаны на определении сияния: ^[1]${\displaystyle L({\vec {r}},{\hat {s}},t)({\frac {W}{m^{2}sr}})}$${\displaystyle {\vec {r}}}$${\displaystyle {\hat {s}}}$${\displaystyle t}$

• Скорость или интенсивность флюенса ${\displaystyle \Phi ({\vec {r}},t)=\int _{4\pi }L({\vec {r}},{\hat {s}},t)d\Omega \left({\frac {W}{m^{2}}}\right)}$
• Флюенс ${\displaystyle F({\vec {r}})=\int _{-\infty }^{+\infty }\Phi ({\vec {r}},t)dt\left({\frac {J}{m^{2}}}\right)}$
• Плотность тока ( поток энергии ) . Это векторный аналог скорости потока энергии, указывающий в преобладающем направлении потока энергии.${\displaystyle {\vec {J}}({\vec {r}},t)=\int _{4\pi }{\hat {s}}L({\vec {r}},{\hat {s}},t)d\Omega \left({\frac {W}{m^{2}}}\right)}$

Уравнение переноса излучения [ править ]

RTE - это дифференциальное уравнение, описывающее яркость . Это может быть получено путем сохранения энергии . Вкратце, RTE утверждает, что луч света теряет энергию из-за расходимости и ослабления (включая как поглощение, так и рассеяние от луча) и получает энергию от источников света в среде и рассеяния, направленного к лучу. Когерентность , поляризация и нелинейность не учитываются. Оптические свойства, такие как показатель преломления , коэффициент поглощения μ _a , коэффициент рассеяния μ _s и анизотропия рассеяния${\displaystyle L({\vec {r}},{\hat {s}},t)}$ ${\displaystyle n}$${\displaystyle g}$считаются неизменными во времени, но могут варьироваться в пространстве. Рассеяние считается упругим. Таким образом, RTE ( уравнение Больцмана ) записывается как: ^[1]

${\displaystyle {\frac {\partial L({\vec {r}},{\hat {s}},t)/c}{\partial t}}=-{\hat {s}}\cdot \nabla L({\vec {r}},{\hat {s}},t)-\mu _{t}L({\vec {r}},{\hat {s}},t)+\mu _{s}\int _{4\pi }L({\vec {r}},{\hat {s}}',t)P({\hat {s}}',{\hat {s}})d\Omega '+S({\vec {r}},{\hat {s}},t)}$

где

• ${\displaystyle c}$ скорость света в ткани, определяемая относительным показателем преломления
• μ _t μ _a + μ _s - коэффициент экстинкции${\displaystyle =}$
• ${\displaystyle P({\hat {s}}',{\hat {s}})}$- фазовая функция, представляющая вероятность рассеяния света с направлением распространения в телесный угол вокруг . В большинстве случаев, фазовая функция зависит только от угла между разрозненными и инцидентами направлениями, то есть . Анизотропия рассеяния может быть выражена как${\displaystyle {\hat {s}}'}$${\displaystyle d\Omega }$${\displaystyle {\hat {s}}}$${\displaystyle {\hat {s}}'}$${\displaystyle {\hat {s}}}$${\displaystyle P({\hat {s}}',{\hat {s}})=P({\hat {s}}'\cdot {\hat {s}})}$${\displaystyle g=\int _{4\pi }({\hat {s}}'\cdot {\hat {s}})P({\hat {s}}'\cdot {\hat {s}})d\Omega }$
• ${\displaystyle S({\vec {r}},{\hat {s}},t)}$ описывает источник света.

Теория диффузии [ править ]

Предположения [ править ]

В RTE, шесть различных независимые переменных определяют блеск в любой пространственной и временной точке ( , и от , полярного угла и азимутального угла от и ). Делая соответствующие предположения о поведении фотонов в рассеивающей среде, можно уменьшить количество независимых переменных. Эти предположения приводят к теории диффузии (и уравнению диффузии) для переноса фотонов. Два предположения позволяют применить теорию диффузии к RTE:${\displaystyle x}$${\displaystyle y}$${\displaystyle z}$${\displaystyle {\vec {r}}}$ ${\displaystyle \theta }$${\displaystyle \phi }$${\displaystyle {\hat {s}}}$${\displaystyle t}$

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

Оба эти предположения требуют среды с высоким альбедо (преимущественно рассеивающей). ^[1]

УПИ в диффузионном приближении [ править ]

Сияние может быть расширено на основе набора сферических гармоник _{n, m} . В теории диффузии яркость считается в значительной степени изотропной, поэтому используются только изотропные и анизотропные члены первого порядка: где _{n, m} - коэффициенты расширения. Сияние выражается 4 членами; один для n = 0 (изотропный член) и 3 члена для n = 1 (анизотропные члены). Используя свойства сферических гармоник и определения плотности потока энергии и тока , изотропный и анизотропный члены могут быть соответственно выражены следующим образом:${\displaystyle Y}$${\displaystyle L({\vec {r}},{\hat {s}},t)\approx \ \sum _{n=0}^{1}\sum _{m=-n}^{n}L_{n,m}({\vec {r}},t)Y_{n,m}({\hat {s}})}$${\displaystyle L}$${\displaystyle \Phi ({\vec {r}},t)}$${\displaystyle {\vec {J}}({\vec {r}},t)}$

• ${\displaystyle L_{0,0}({\vec {r}},t)Y_{0,0}({\hat {s}})={\frac {\Phi ({\vec {r}},t)}{4\pi }}}$
• ${\displaystyle \sum _{m=-1}^{1}L_{1,m}({\vec {r}},t)Y_{1,m}({\hat {s}})={\frac {3}{4\pi }}{\vec {J}}({\vec {r}},t)\cdot {\hat {s}}}$

Следовательно, мы можем аппроксимировать яркость как ^[1]

${\displaystyle L({\vec {r}},{\hat {s}},t)={\frac {1}{4\pi }}\Phi ({\vec {r}},t)+{\frac {3}{4\pi }}{\vec {J}}({\vec {r}},t)\cdot {\hat {s}}}$

Подставляя указанное выше выражение для яркости, RTE можно соответственно переписать в скалярной и векторной формах следующим образом (член рассеяния RTE интегрируется по полному телесному углу. Для векторной формы RTE перед оценкой умножается на направление ). : ^[1]${\displaystyle 4\pi }$${\displaystyle {\hat {s}}}$

${\displaystyle {\frac {\partial \Phi ({\vec {r}},t)}{c\partial t}}+\mu _{a}\Phi ({\vec {r}},t)+\nabla \cdot {\vec {J}}({\vec {r}},t)=S({\vec {r}},t)}$

${\displaystyle {\frac {\partial {\vec {J}}({\vec {r}},t)}{c\partial t}}+(\mu _{a}+\mu _{s}'){\vec {J}}({\vec {r}},t)+{\frac {1}{3}}\nabla \Phi ({\vec {r}},t)=0}$

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

Уравнение диффузии [ править ]

Используя второе предположение теории диффузии, отметим, что частичное изменение плотности тока на одной транспортной длине свободного пробега незначительно. Векторное представление теории диффузии RTE сводится к закону Фика , который определяет плотность тока в терминах градиента плотности потока энергии. Подстановка закона Фика в скалярное представление RTE дает уравнение диффузии: ^[1]${\displaystyle {\vec {J}}({\vec {r}},t)}$ ${\displaystyle {\vec {J}}({\vec {r}},t)={\frac {-\nabla \Phi ({\vec {r}},t)}{3(\mu _{a}+\mu _{s}')}}}$

${\displaystyle {\frac {1}{c}}{\frac {\partial \Phi ({\vec {r}},t)}{\partial t}}+\mu _{a}\Phi ({\vec {r}},t)-\nabla \cdot [D\nabla \Phi ({\vec {r}},t)]=S({\vec {r}},t)}$

${\displaystyle D={\frac {1}{3(\mu _{a}+\mu _{s}')}}}$- коэффициент диффузии, а μ ' _s μ _s - приведенный коэффициент рассеяния. Примечательно, что в уравнении диффузии нет явной зависимости от коэффициента рассеяния. Вместо этого в выражении для появляется только приведенный коэффициент рассеяния . Это приводит к важным отношениям; диффузия не изменяется, если анизотропия рассеивающей среды изменяется, в то время как приведенный коэффициент рассеяния остается постоянным. ^[1]${\displaystyle =(1-g)}$
${\displaystyle D}$

Решения уравнения диффузии [ править ]

Для различных конфигураций границ (например, слоев ткани) и источников света уравнение диффузии может быть решено путем применения соответствующих граничных условий и определения источника в зависимости от ситуации.${\displaystyle S({\vec {r}},t)}$

Точечные источники в бесконечных однородных средах [ править ]

В этом разделе представлено решение уравнения диффузии для простого случая точечного источника с короткими импульсами в бесконечной однородной среде. Член источника в уравнении диффузии принимает вид , где - позиция, в которой измеряется плотность потока энергии, а - позиция источника. Пик импульса во время . Уравнение диффузии решается для плотности потока энергии, чтобы дать${\displaystyle S({\vec {r}},t,{\vec {r'}},t')=\delta ({\vec {r}}-{\vec {r'}})\delta (t-t')}$${\displaystyle {\vec {r}}}$${\displaystyle {\vec {r'}}}$${\displaystyle t'}$

${\displaystyle \Phi ({\vec {r}},t;{\vec {r'}},t)={\frac {c}{[4\pi Dc(t-t')]^{3/2}}}\exp \left[-{\frac {\mid {\vec {r}}-{\vec {r'}}\mid ^{2}}{4Dc(t-t')}}\right]\exp[-\mu _{a}c(t-t')]}$

Термин представляет собой экспоненциальный спад флюенса из-за поглощения в соответствии с законом Бера . Остальные члены представляют собой уширение из-за рассеяния. При таком решении произвольный источник можно охарактеризовать как суперпозицию короткоимпульсных точечных источников. Исключение изменения во времени из уравнения диффузии дает следующее для не зависящего от времени точечного источника :${\displaystyle \exp \left[-\mu _{a}c(t-t')\right]}$${\displaystyle S({\vec {r}})=\delta ({\vec {r}})}$

${\displaystyle \Phi ({\vec {r}})={\frac {1}{4\pi Dr}}\exp(-\mu _{\mathrm {eff} }r)}$

${\displaystyle \mu _{\mathrm {eff} }={\sqrt {\frac {\mu _{a}}{D}}}}$- эффективный коэффициент затухания и указывает скорость пространственного затухания флюенса. ^[1]

Граничные условия [ править ]

Скорость потока на границе [ править ]

Учет граничных условий позволяет использовать уравнение диффузии для характеристики распространения света в среде ограниченного размера (где необходимо учитывать границы раздела между средой и окружающей средой). Чтобы начать рассмотрение границы, можно рассмотреть, что происходит, когда фотоны в среде достигают границы (т. Е. Поверхности). Интегрированная по направлениям яркость на границе и направленная в среду равна интегрированной по направлениям яркости на границе и направленной из среды, умноженной на коэффициент отражения :${\displaystyle R_{F}}$

${\displaystyle \int _{{\hat {s}}\cdot {\hat {n}}<0}L({\vec {r}},{\hat {s}},t){\hat {s}}\cdot {\hat {n}}d\Omega =\int _{{\hat {s}}\cdot {\hat {n}}>0}R_{F}({\hat {s}}\cdot {\hat {n}})L({\vec {r}},{\hat {s}},t){\hat {s}}\cdot {\hat {n}}d\Omega }$

где перпендикулярно границе и направлено от нее. Диффузионное приближение дает выражение для энергетической яркости через флюенс и плотность тока . Оценка вышеуказанных интегралов после подстановки дает: ^[3]${\displaystyle {\hat {n}}}$${\displaystyle L}$${\displaystyle \Phi }$${\displaystyle {\vec {J}}}$

${\displaystyle {\frac {\Phi ({\vec {r}},t)}{4}}+{\vec {J}}({\vec {r}},t)\cdot {\frac {\hat {n}}{2}}=R_{\Phi }{\frac {\Phi ({\vec {r}},t)}{4}}-R_{J}{\vec {J}}({\vec {r}},t)\cdot {\frac {\hat {n}}{2}}}$

• ${\displaystyle R_{\Phi }=\int _{0}^{\pi /2}2\sin \theta \cos \theta R_{F}(\cos \theta )d\theta }$
• ${\displaystyle R_{J}=\int _{0}^{\pi /2}3\sin \theta (\cos \theta )^{2}R_{F}(\cos \theta )d\theta }$

Подстановка закона Фика ( ) дает на расстоянии от границы z = 0, ^[3]${\displaystyle {\vec {J}}({\vec {r}},t)=-D\nabla \Phi ({\vec {r}},t)}$

${\displaystyle \Phi ({\vec {r}},t)=A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}}$

• ${\displaystyle A_{z}=2D{\frac {1+R_{\mathrm {eff} }}{1-R_{\mathrm {eff} }}}}$
• ${\displaystyle R_{\mathrm {eff} }={\frac {R_{\Phi }+R_{J}}{2-R_{\Phi }+R_{J}}}}$

Экстраполированная граница [ править ]

Желательно определить границу нулевого флюенса. Однако коэффициент плотности потока энергии на физической границе, как правило, не равен нулю. Экстраполированная граница при _b, для которой плотность потока энергии равна нулю, может быть определена для установления источников изображения. Используя приближение ряда Тейлора первого порядка ,${\displaystyle \Phi (z=0,t)}$${\displaystyle z}$

${\displaystyle \left.\Phi (z=-A_{z},t)\approx \Phi (z=0,t)-A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}\right|_{z=0}}$

который равен нулю, поскольку . Таким образом, по определению, _b должно быть _z, как определено выше. Примечательно, что когда показатель преломления одинаков с обеих сторон границы, _F равно нулю, а экстраполированная граница находится в точке _b . ^[3]${\displaystyle \Phi ({\vec {r}},t)=A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}}$${\displaystyle z}$${\displaystyle -A}$${\displaystyle R}$${\displaystyle z}$${\displaystyle =-2D}$

Луч карандаша, обычно падающий на полубесконечную среду [ править ]

Используя граничные условия, можно приближенно охарактеризовать диффузное отражение для стержневого пучка, обычно падающего на полубесконечную среду. Луч будет представлен как два точечных источника в бесконечной среде следующим образом (Рисунок 2): ^{[1] [4]}

1. Установите анизотропию рассеяния ₂ для рассеивающей среды и установите новый коэффициент рассеяния μ _{s2 равным} исходному μ _s1, умноженному на ₁ , где ₁ - исходная анизотропия рассеяния.${\displaystyle g}$${\displaystyle =0}$${\displaystyle (1-g}$${\displaystyle )}$${\displaystyle g}$
2. Преобразуйте карандашный луч в изотропный точечный источник на глубине, равной одной транспортной длине свободного пробега «ниже поверхности и мощность = ».${\displaystyle l}$${\displaystyle a}$
3. Реализуйте экстраполированное граничное условие, добавив источник изображения противоположного знака над поверхностью в точке ' _b .${\displaystyle l}$${\displaystyle +2z}$

Два точечных источника можно охарактеризовать как точечные источники в бесконечной среде через

${\displaystyle \Phi _{\infty }(r,\theta ,z;r',\theta ',z')={\frac {1}{4\pi D\rho }}\exp(-\mu _{\mathrm {eff} }\rho )}$

${\displaystyle \rho }$- расстояние от точки наблюдения до местоположения источника в цилиндрических координатах. Линейная комбинация вкладов плотности потока энергии от двух источников изображения равна${\displaystyle (r,\theta ,z)}$${\displaystyle (r',\theta ',z')}$

${\displaystyle \Phi (r,\theta ,z;r',\theta ',z')=a'\Phi _{\infty }(r,\theta ,z;r',\theta ',z')-a'\Phi _{\infty }(r,\theta ,z;r',\theta ',-z'-2z_{b})}$

Это можно использовать для получения коэффициента диффузного отражения _{d с} помощью закона Фика:${\displaystyle R}$${\displaystyle (r)}$

${\displaystyle \left.R_{d}(r)=D{\frac {\partial \Phi }{\partial z}}\right|_{z=0}={\frac {a'z'(1+\mu _{\mathrm {eff} }\rho _{1})\exp(-\mu _{\mathrm {eff} }\rho _{1})}{4\pi \rho _{1}^{3}}}+{\frac {a'(z'+4D)(1+\mu _{\mathrm {eff} }\rho _{2})\exp(-\mu _{\mathrm {eff} }\rho _{2})}{4\pi \rho _{2}^{3}}}}$

${\displaystyle \rho _{1}}$- расстояние от точки наблюдения до источника в точке, а - расстояние от точки наблюдения до источника изображения в точке _b . ^{[1] [4]}${\displaystyle (r,0,0)}$${\displaystyle (0,0,z')}$${\displaystyle \rho _{2}}$${\displaystyle (0,0,-z'-2z}$${\displaystyle )}$

Решения теории диффузии против моделирования Монте-Карло [ править ]

Моделирование переноса фотонов методом Монте-Карло, хотя и требует много времени, позволяет точно предсказать поведение фотонов в рассеивающей среде. Допущения, связанные с характеристикой поведения фотонов с помощью уравнения диффузии, порождают неточности. Как правило, диффузионное приближение становится менее точным, поскольку коэффициент поглощения μ _a увеличивается, а коэффициент рассеяния μ _s уменьшается. ^{[5] [6]} Для пучка фотонов, падающего на среду ограниченной глубины, ошибка из-за приближения диффузии наиболее заметна в пределах одной транспортной длины свободного пробега от места падения фотона (где яркость еще не изотропна) (рис. ).
Среди шагов в описании карандашного луча, падающего на полубесконечную среду с помощью уравнения диффузии, преобразование среды из анизотропной в изотропную (шаг 1) (рисунок 4) и преобразование луча в источник (шаг 2) (рисунок 5) генерировать больше ошибок, чем преобразование из одного источника в пару источников изображения (шаг 3) (рисунок 6). Шаг 2 генерирует наиболее значительную ошибку. ^{[1] [4]}

• Рис. 3: Зависимость диффузного отражения от радиуса падающего пучка-карандаша, определенная с помощью моделирования методом Монте-Карло (красный цвет), и зависимости диффузного отражения от радиуса от двух изотропных точечных источников, как определено решением теории диффузии для RTE (синий). Транспортная длина свободного пробега 0,1 см.

• Рис. 4. Зависимость коэффициента диффузного отражения от радиуса падающего пучка для анизотропной (синий) и изотропной (красный) среды.

• Рисунок 5: Зависимость диффузного отражения от радиуса от источника фотонов для стержневого луча (синий) и точечного изотропного источника (красный).

• Рис. 6: Зависимость диффузного отражения от радиуса от источника фотонов для изотропного точечного источника, как показано на основе решения RTE (синий) и моделирования Монте-Карло (красный).

См. Также [ править ]

• Метод Монте-Карло для переноса фотонов
• Радиационный перенос

Ссылки [ править ]

Дальнейшее чтение [ править ]

• Л.В. Ван и Х.И. Ву (2007). Биомедицинская оптика . Вайли. ISBN 978-0-471-74304-0.

• С.Г. Проскурин (2011). «Квантовая электроника. 41 402». Квантовая электроника . 41 (5): 402–406. DOI : 10.1070 / QE2011v041n05ABEH014597 . (2011)