光を勉強するブログ

光と流体に関して。

光の伝搬現象における輸送方程式の概要と導出

概要と方針

 輸送方程式(Radiative Transfer Equation : RTE)とはボルツマン方程式の一種である偏微分方程式です。光計測における基礎方程式として知られ、光CTや拡散光トモグラフィに利用されていものです。

   さて、この記事の方針ですが、題目の通り輸送方程式を光輸送の観点から解析学に知見の浅い方でも分かるよう導出していきます。というのも、以前、ネットで検索したら中性子関連のものだらけで光に関するものが(あまり)見つからなく、あったとしても難解なものだらけでした。初学者であり、解析学に知見の浅い自分でも分かるような、直観でイメージできる記事があったら…と考え、ないなら自分で書いてしまおうと思ったことが動機です。

 

導出

 輸送方程式はエネルギー保存則を考えることで導くことが叶います。これは光の伝搬現象を導く際にも同様です。具体的には光強度の空間的・時間的な変化は光の吸収と散乱によっておこる。ということを記述します。式で表すと以下のようになります。

\begin{align}
\displaystyle I(\boldsymbol{r}+\Delta l\boldsymbol{s}, \boldsymbol{s}, t+\Delta t) = I(\boldsymbol{r}, \boldsymbol{s}, t) + \Delta I_{吸収} + \Delta I_{散乱}  \tag{1}
\end{align}

  {I(\boldsymbol{r}, \boldsymbol{s}, t)}はある位置\boldsymbol{r} \in \boldsymbol{R}^3、方向 \boldsymbol{s}\in \boldsymbol{S}^2、時刻tにおける光強度を表しており、 \Delta I_{吸収} + \Delta I_{散乱}はそれぞれ吸収と散乱による光強度の変化を表します。この式の左辺では\boldsymbol{s}方向に微小距離\Delta l移動し、時刻\Delta t経過したときの光強度の微小変化を示し、右辺で、その変化量は吸収と散乱によると表しています。まず、この左辺についてテイラー展開をもちいて一次近似を施します。すると

\begin{align}
\displaystyle L.H.S &=  I(\boldsymbol{r}+\Delta l\boldsymbol{s}, \boldsymbol{s}, t+\Delta t) \\ \displaystyle &= \sum_{k=0}^1 \frac{1}{k!} (\Delta l \boldsymbol{s} \cdot \nabla + \Delta t \frac{\partial}{\partial t} )^k I(\boldsymbol{r}, \boldsymbol{s}, t) \\ \displaystyle  &= I(\boldsymbol{r}, \boldsymbol{s}, t)+\Delta l \boldsymbol{s} \cdot \nabla  I(\boldsymbol{r}, \boldsymbol{s}, t) + \Delta t \frac{\partial}{\partial t}I(\boldsymbol{r}, \boldsymbol{s}, t) \tag{2}
\end{align}

そして(1)式右辺の I(\boldsymbol{r}, \boldsymbol{s}, t)を導出された(2)式左辺に移項することによって(1)式は次のようになります。

\begin{align}
\displaystyle \Delta l \boldsymbol{s} \cdot \nabla  I(\boldsymbol{r}, \boldsymbol{s}, t) + \Delta t \frac{\partial}{\partial t}I(\boldsymbol{r}, \boldsymbol{s}, t) = \Delta I_{吸収} + \Delta I_{散乱}  \tag{3}
\end{align}

 さて、次に光の吸収による光強度の変化I_{吸収}をみます。

               f:id:Iwalabo:20180830112043p:plain

上図のように、ある媒体に光強度 I(\boldsymbol{r}, \boldsymbol{s}, t)で照射された光は、微小距離進む中で媒質による吸収の影響を受け、dIだけ変化します。
 また、I_{吸収}\Delta lの関係は以下のような式で記述されることが知られています。

\begin{align}
\displaystyle \Delta I_{吸収}=-\mu_a I(\boldsymbol{r}, \boldsymbol{s}, t) \Delta l \tag{4}
\end{align}

 これは、吸収による光強度の変化量が変化前の光強度に比例し、その比例定数が吸収係数\mu_a[mm^{-1}]であることを示しています。そして符号が負なのは吸収によって光強度が減少していることを示しています。

…と書きましたがこの(4)は正確ではありません。

 正確にはI_{吸収}dIとおき、\Delta ldlとおきかえることが叶うときに成り立ちます。また、そうすることでこの(4)は単純な微分方程式となり、変数分離型として解くことで光トポグラフィにおける重要な式、Beer Lambert則を得ることができます。

 

iwalabo.hatenablog.com

 

 

 次に散乱による光強度の変化I_{散乱}をみます。散乱の影響は吸収によるものと比べて少し複雑になります。また、輸送方程式において散乱による項は2つ現れます。散乱による光強度の増加項と散乱による光強度の減少光の2つです。

              f:id:Iwalabo:20180830112114p:plain

  Fig.2は散乱によって生ずる光強度の減少を表しています。\boldsymbol{s}方向からの光の一部が媒質内で粒子に衝突し、\boldsymbol{s}’方向へと経路を変えることで本来の\boldsymbol{s}方向の光が減少します。この現象は以下のように記述されます。

\begin{align}
\displaystyle \Delta I_{散乱-}=-\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s},\boldsymbol{s}')I(\boldsymbol{r},\boldsymbol{s},t)d\boldsymbol{s}' \Delta l \tag{5}
\end{align}

 ここで、符号の負は先ほどと同様に光強度の減少を示し、\mu_sは散乱係数[mm^{-1}]を示します。そして、p(\boldsymbol{s},\boldsymbol{s}')は方向\boldsymbol{s}から方向\boldsymbol{s}'への散乱確率を示す位相関数になります。この位相関数と光強度の積をとり、\boldsymbol{s}'に関して全立体角で積分することで\boldsymbol{s}方向の光強度が別方向にどれだけ散乱されるかわかります。また、位相関数p(\boldsymbol{s},\boldsymbol{s}')は全立体角で積分すると1になるよう定義されており、これを用いてこの積分式を整理すると、

\begin{align}
\displaystyle \Delta I_{散乱-} &= -\mu_sI(\boldsymbol{r},\boldsymbol{s},t) \int_{\boldsymbol{s}^2} p(\boldsymbol{s},\boldsymbol{s}')d\boldsymbol{s}' \Delta l \\ \displaystyle &= -\mu_s I(\boldsymbol{r},\boldsymbol{s},t)\Delta l \tag{6}
\end{align}

 次に散乱による光強度の増加をみます。 

                  f:id:Iwalabo:20180830114246p:plain

  Fig.3は散乱によって生ずる光強度の増加を表しています。減少によるものと同様に、\boldsymbol{s}'方向からの光の一部が媒質内で粒子に衝突し、\boldsymbol{s}方向へと経路を変えることで本来の\boldsymbol{s}方向の光が増加します。この現象は以下のようになります。

\begin{align}
\displaystyle \Delta I_{散乱-}=\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'} \Delta l \tag{7}
\end{align}
  先ほどと異なるのは位相関数の変数の順序と、光強度Iの方向変数が\boldsymbol{s}'となっており、増加現象を示しているため符号が正となっている点です。

 最後に、これまで導出した吸収・散乱項(4)、(6)、(7)式を(3)式に代入すると、
\begin{align}
\displaystyle \Delta l \boldsymbol{s} \cdot \nabla I(\boldsymbol{r},\boldsymbol{s}, t) + \Delta t \frac{\partial}{\partial t} I(\boldsymbol{r}, \boldsymbol{s}, t) =  -\mu_a I(\boldsymbol{r}, \boldsymbol{s}, t) \Delta l -\mu_s I(\boldsymbol{r},\boldsymbol{s},t)\Delta l +\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'} \Delta l \tag{8}
\end{align}
となります。そして全項を\Delta lで割り、媒質中の光速度c=\frac{\Delta l}{\Delta t}の関係式を用いると、
\begin{align}
\displaystyle \boldsymbol{s} \cdot \nabla I(\boldsymbol{r},\boldsymbol{s}, t) + \frac{1}{c} \frac{\partial}{\partial t} I(\boldsymbol{r}, \boldsymbol{s}, t) =  -\mu_a I(\boldsymbol{r}, \boldsymbol{s}, t)-\mu_s I(\boldsymbol{r},\boldsymbol{s},t)+\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'}\tag{9}
\end{align}
となり、整理することで輸送方程式を得ることができます。
\begin{align}
\displaystyle\frac{1}{c} \frac{\partial}{\partial t} I(\boldsymbol{r}, \boldsymbol{s}, t)+ \boldsymbol{s} \cdot \nabla I(\boldsymbol{r},\boldsymbol{s}, t)+(\mu_a +\mu_s)I(\boldsymbol{r}, \boldsymbol{s}, t) =\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'}\tag{10}
\end{align}

 

  さて、輸送方程式を導くことが出来ました。なるべく直感的にわかりやすく、かつ曖昧な議論を避けたつもりですが、良く分からない点や間違いがあったら申し訳ありません。

  今後は輸送方程式を用いて光拡散近似を行っていこうと考えております。興味があったらご一読お願いします。

 

iwalabo.hatenablog.com