光の輸送現象とその逆問題
町田学

Department of Mathematics, University of Michigan

ランダム媒質：雲ランダム媒質：生体組織

ランダム媒質：星間物質

ランダム媒質中の光ℓ∗ aBB B

Q

QsQ
吸収µa = 単位長さ当りに吸収される割合

散乱µs = 単位長さ当りに散乱される割合
f (ŝ, ŝ′ ) = 散乱されたときに、ŝ′ 方向の光が ŝ 方向を向く確率

光の伝播：水光の伝播：微量

醤油

牛乳

光の伝播：少量

醤油

牛乳

光の伝播：もっと

醤油

牛乳

光の伝播：さらに

醤油

牛乳

輸送平均自由行程
頻繁に散乱される場合（µa ≪ µs ）を考える。
散乱長:ℓs =平均自由行程:∫g=S21µsℓ∗ =

1µa + µs (1 − g)(ŝ · ŝ′ )f (ŝ, ŝ′ ) dŝ′

等方散乱 =
前方散乱 =
生体組織では ℓ∗ ≫ ℓs

g=0 ⇒0≪g<1ℓ∗ < ℓs
⇒ ℓ∗ ≫ ℓs

光強度I(r, ŝ, t) = 位置 r (∈ R3 ) 時刻 t で ŝ (∈ S2 ) 方向の、単位面積当り単位時間当りの輻射ŝdŝ
dσdσ
輻射エネルギー = I(r, ŝ, t)(dσ · ŝ)dŝdt

エネルギー輸送
断面積 ∆A で奥行き ∆ℓ = c∆t (c は物質中の光速度) の微小直方体

I(r + ∆ℓŝ, ŝ, t + ∆t)∆A∆t = I(r, ŝ, t)∆A∆t − 散乱と吸収
LHS = [I(r, ŝ, t + ∆t) + ∆ℓŝ · ∇I(r, ŝ, t + ∆t)] ∆A∆t][∂I= I(r, ŝ, t) + ∆t (r, ŝ, t) + ∆ℓŝ · ∇I(r, ŝ, t) ∆A∆t∂tŝ∆t ∂I
散乱と吸収+ ŝ · ∇I = −∆ℓ ∂t∆A∆ℓ∆t∆ℓ = c∆t∆A光の吸収と散乱

吸収µa I(r, ŝ)

散乱−ŝ から別方向へ∫

別方向からŝ へ∫= µsf (ŝ′ , ŝ)I(r, ŝ) dŝ′ − µsf (ŝ, ŝ′ )I(r, ŝ′ ) dŝ′S2S2∫
= µs I(r, ŝ) − µsf (ŝ, ŝ′ )I(r, ŝ′ ) dŝ′S2

輻射輸送方程式
微小直方体でのエネルギーバランス

1 ∂I(r, ŝ, t)+ ŝ · ∇I(r, ŝ, t) + (µa + µs )I(r, ŝ, t)c∂t∫= µsf (ŝ, ŝ′ )I(r, ŝ′ , t) dŝ′ + S(r, ŝ, t)

S2S(r, ŝ, t) = 光源⇐輻射輸送方程式＝線形ボルツマン方程式

１階の微分方程式だが積分項があり変数の数は６つ (x, y, z, θ, φ, t)様々な場面で登場

ランダム媒質中の光
海洋、流氷
原子炉中の中性子やガンマ線地震波

マクスウェル方程式から輻射輸送方程式を導出する試み
摂動論
均質化法

拡散近似
µa ≪ µs でしかも r が境界や光源から遠く離れているときに∫1
I(r, ŝ, t) dŝ
u(r, t) =c S2は拡散方程式に従う。

∂u(r, t)− ∇ · (D∇)u(r, t) + cµa u(r, t) = SDA (r, t)
∂t
D=cℓ∗
3

階層
ミクロマクスウェルメソ=⇒λ ≪ ℓ∗

輻射輸送λ: 光の波長
L: 系のサイズℓ∗ マクロ=⇒ℓ ≪L∗aQBQsQ

B B
-L

拡散今日主に使う式

時間に依存しない輻射輸送方程式を考える。また、不純物は球対称で散乱は相対的な角度だけに依るとする。

∫ŝ · ∇I(r, ŝ) + µt I(r, ŝ) = µs

S2f (ŝ · ŝ′ )I(r, ŝ′ ) dŝ′ + S(r, ŝ)

全減衰： µt = µa + µs
１階の微分方程式だが積分項があり変数の数は５つ (x, y, z, θ, φ)

逆問題??

観測逆輸送問題１⇒

境界での I(r, ŝ)µa , µs , f

つまり、I|Γ+⇒µa

Γ± = {(r, ŝ) ∈ ∂Ω × S2 ; ±ŝ · ν̂(r) > 0}、ν̂(r) は r での外向き単位
法線ベクトルµt = µa + µs で∫ŝ · ∇I(r, ŝ) + µt I(r, ŝ) = µs

S2f (ŝ · ŝ′ )I(r, ŝ′ ) dŝ′ + S(r, ŝ)逆輸送問題２方程式は線形。つまり、与えられた µa に対して
（順問題）

µa7→

境界の出射光 I(r, ŝ)は線形。しかし、µa については非線形。つまり（逆問題）

境界の出射光 I(r, ŝ)

7→µaは非線形になる。

しかも、境界の I の微小変化に対して、求まる µa は大きく変化する (ill-posed)

光トモグラフィー

レーザー−→−→−→µa−→−→−→−→−→−→CCD カメラ（測定データ） =⇒CCD カメラµa

X 線 CT vs 光トモグラフィー

吸収吸収散乱逆問題の理論
順問題：?
存在（実験）
逆問題：

一意性
安定性?

復元光トモグラフィーは不安定拡散方程式 vs 輻射輸送方程式
拡散光トモグラフィー： 実験あり
輸送光トモグラフィー： （基本的に）未知
拡散光トモグラフィーによる乳がん検査

Penn System
http://www.physics.upenn.edu/yodhlab/

Choe, et al. (2005)

拡散光トモグラフィーによる脳機能計測
生後２ヶ月の幼児に母親が子供向けの本を読み聞かせる。

normal speech

reversenone
Hitachi, Ltd. (2003)

医用イメージング

機能解剖

Optical Tomography

OT:
CT:
MRI:

CT (Computed Tomography)
MRI

tissue oxygenation (Hb/HbO2 )骨軟組織OT: 安全、低価格

医用イメージングの歴史

X 線 Röntgen (1895)[1901 年ノーベル物理学賞]
超音波 水中ソナー検出 (1912)、医用超音波 (1940 年代後半)
CT Cormack (1963) and Hounsfield (1971)[1979 年ノーベル生理・医学賞]
MRI Lauterbur and Mansfield (1973)[2003 年ノーベル生理・医学賞]
光トモグラフィー 1990 年代初頭から

Part I順問題

拡散近似の破綻

媒質が薄い、境界付近、光源付近、吸収 µa が大きい、場合には拡散近似が正当化できない
マクロ領域

⇒

メソ領域

拡散方程式

⇒輻射輸送方程式

輻射輸送方程式を解く方法
解析的に解ける場合
１次元（z だけに依存）
1
等方散乱 (p = 4π) in Rn (n ≥ 2)

数値的手法
モンテカルロdiscrete ordinates
PN
有限要素法
逆問題のために、できるだけ解析的に解きたい！

１次元µ̄t と µ̄s は定数とする。
∂I(z, ŝ)cos θ+ µ̄t I(z, ŝ) = µ̄s∂z∫
S2f (ŝ · ŝ′ )I(z, ŝ′ ) dŝ′ + S(z, ŝ)

解 I(z, ŝ) は（初等関数で書ける）固有モードの重ねあわせで与えられる∑I(z, ŝ) =cν ϕν (ŝ)e−z/νν

係数 cν は境界条件と S から決まる

Case (1960)
Mika (1961), McCormick and Kuščer (1966)

３次元への拡張の試み“Although there have been many attempts to extend thesemethods (the integral transform and singular eigenfunction
methods) to two- and three-dimensional problems, these
extensions have usually encountered extreme mathematicalcomplexity and have met with only marginal success.”[Duderstadt-Martin, Transport Theory (1979), p. 122]
M. M. R. Williams, J. Math. Phys. (1968)
A. G. Gibbs, J. Math. Phys. (1969)
G. Garrettson and A. Leonard, J. Math. Phys. (1970)
K. M. Case and R. D. Hazeltine, J. Math. Phys. (1970)
C. J. Cannon, J. Quant. Spec. Rad. Trans. (1973)
G. C. Pomraning, Ann. Nucl. Energy (1996)

回転座標系の方法
Iν (z, ŝ) = ϕν (ŝ)e−z/ν⇒Iν (r, ŝ) = ϕν (ŝ(k̂))e−r·k̂/νsin θ(k̂) cos φ(k̂)
ŝ(k̂) =  sin θ(k̂) sin φ(k̂)  ,cos θ(k̂)

k̂ · k̂ = 1θ(k̂) と φ(k̂) は、z 軸が k̂ 方向を向く座標系で測った仰角と方位角。
cf., V. A. Markel (2004): Ylm (ŝ) ⇒ Ylm (ŝ; k̂)
M. M., G. Y. Panasyuk, J. C. Schotland, and V. A. Markel (2010)３次元→１次元ŝ · ∇ϕν (ŝ(k̂))e−r·k̂/νŝ · k̂ϕν (ŝ(k̂))e−r·k̂/ννcos θ(k̂)= −ϕν (ŝ(k̂))e−r·k̂/νν= −座標を回転： k̂ → ẑcos θϕν (ŝ)e−z/νν∂= cos θ ϕν (ŝ)e−z/ν∂z⇒ −３次元での解析解Theorem
Let I(r, ŝ) be the solution to the RTE in R3 . We assumelN ∑∑′fl Ylm (ŝ)Ylm (ŝ′ ). Then there exists a set off (ŝ · ŝ ) =l=0 m=−leigenmodes Iνm (r, ŝ; q) such that
I(r, ŝ) =N∑∑∫m=−NνR2mcm
ν (q)Iν (r, ŝ; q) dq,where cm
ν (q) are determined by the source term S(r, ŝ).
M.M. (2014)

エネルギー密度エネルギー密度の計算

1u(r) =c∫I(r, ŝ) dŝS2U(z)U(µt z) v / µt2 Sa

102µa=0.03, µs=100, f1=0Fourier transformµa=0.03, µs=100, f1=0.3µa=0.3, µs=100, f1=0.3

101Monte Carlo
10010-10246µt z810Part II逆問題出射光→吸収
一般的に時間に依存する場合を考える。
ランダム媒質はなめらかな境界 ∂Ω で囲まれた有界な領域 Ω にあるとする。
（逆問題）

外向きの出射光 I|Γ+ ×(0,T ) = AIΓ− ×(0,T )ここで{}
Γ± = (r, ŝ) ∈ ∂Ω × S2 ; ±ŝ · ν̂(r) > 0ν̂(r) は r での外向き単位法線ベクトル
A は反射能演算子と呼ばれる[反射能＝反射光／入射光 = albedo]
A は線形だが、I|Γ+ ×(0,T ) 7→ µa は非線形→µa一意性と安定性

一意性
Choulli-Stefanov (1996, 1999), Stefanov (2003),
Stefanov-Tamasan (2009), Bal-Jollivet (2009),
Klibanov-Pamyatnykh (2006),...
安定性
Romanov (1997, 1998), Wang (1999), Tamasan (2002),
Bal-Jollivet (2008, 2009, 2010), Bal (2008),
Bal-Langmore-Monard (2008), Stefanov-Uhlmann (2003),
Langmore (2008), McDowall-Stefanov-Tamasan (2010),...

大域リプシッツ安定性
これまで、∥A∥ に対する Hölder stability⇒ 一つの境界値に対する Global Lipschitz stability

TheoremI k = I(µkt )(r, ŝ, t) (k = 1, 2) は輻射輸送方程式の解。ある γ ∈ Rd ,
̸= 0 について V ⊂ {ŝ; γ · ŝ > 0} となるように V をとる。I k は十分
なめらかで ∥µkt ∥L∞ (Ω×V ) 、∥µs ∥L∞ (Ω×V ) ≤ M とする。T(t ∈ (0, T )) は十分に大きく、初期条件は (r, ŝ) ∈ Ω × V で正、a(r, ŝ) > 0 とする。以下の定数 C = C(M ) > 0 が存在する。
(∫ T ∫∥µ1t − µ2t ∥L2 (Ω×V ) ≤ C(ν̂(r) · ŝ)|∂t (I 1 − I 2 )|2 dSdŝdt0) 21.

Γ+M.M. and M. Yamamoto (2014)

輸送光トモグラフィーの復元アルゴリズム

Dorn (1998)
Klose and Hielscher (1999)
Wright, Schweiger, and Arridge (2007)
Tarvainen, Vauhkonen, and Arridge (2008)
Gonzalez-Rodriguez and Kim (2009)···

非接触光トモグラフィー
ここからはまた時間に依存しない場合

Schotland and Markel (2007)

輸送光トモグラフィー
η > −1, supp η ⊆ Ba ⊆ Ω とする。

∫ŝ · ∇I(r, ŝ) + µ̄t [1 + η(r)]I(r, ŝ) = µ̄sŝ · ∇I0 (r, ŝ) + µ̄t I0 (r, ŝ) = µ̄s∫I(r, ŝ) = I0 (r, ŝ) + µ̄t= I0 (r, ŝ) + µ̄t∫RS2∫S2f (ŝ · ŝ′ )I(r, ŝ′ ) dŝ′f (ŝ · ŝ′ )I0 (r, ŝ′ ) dŝ′G(r, ŝ; r′ , ŝ′ )η(r′ )I(r′ , ŝ′ ) dr′ dŝ′
3R3

G(r, ŝ; r′ , ŝ′ )η(r′ )I0 (r′ , ŝ′ ) dr′ dŝ′ + . . .

境界で (I0 − I)|Γ+ を測定⇒

未知の η を復元

ボルン級数と逆ボルン級数
データ：Φ = (I0 − I)|Γ+ とする。

ボルン級数：

Φ = K1 η + K2 η ⊗ η + K3 η ⊗ η ⊗ η + . . .
逆ボルン級数：η = K1 Φ + K2 Φ ⊗ Φ + K3 Φ ⊗ Φ ⊗ Φ + . . .
Theorem
Kj : L1 (Γ+ × · · · × Γ+ ) → L∞ (Ba ) (j = 1, 2, . . . ) であり、級数は∥K1 ∥ と ∥K1 Φ∥L∞ (Ba ) が十分に小さいときに収束する。

ボルン近似η = K1 Φ + K2 Φ ⊗ Φ + . . .
= K1+ Φ + [−K1 K2 K1 ⊗ K1 ] Φ ⊗ Φ + . . .
≈ K1+ Φill-posed⇒⇒

K1 の inverse は存在しない
K1 = K1+ は pseudoinverse
低周波フィルターのようなもの

復元画像
高速アルゴリズム：ν を求める対角化と係数を求める連立方程式、２次元逆フーリエ変換
大きなデータ：108 程度の光源と検出器の組は容易（cf., 104 ）レモン片の復元画像

レモンとレンコン

レモンとレンコン

さらに発展させるには境界がある場合のグリーン関数の解析解→連立方程式を解かずに係数が求まる
非線形項を使った η の復元→トポロジカル縮約＝各項で１つのダイヤグラムだけが寄与µt と µs の同時復元
偏光を考慮→ベクトル輻射輸送方程式