【2005年_H Hammer】Reconstruction of spatially inhomogeneous dielectric tensors through optical tomography

https://opg.optica.org/viewmedia.cfm?seq=0&uri=josaa-22-2-250 透明な媒体内部に存在する弱い異方性をもつ不均質誘電テンソルを再構成する方法が提案されている。積分幾何学の数学理論を実用的な枠組みに落とし込み、6つの平面断層撮影サイクルから得られる偏光変換データに対してスカラーRadon逆変換を行うことで、誘電テンソル場を完全に決定できるようにした。さらに、材料の不均質性や異方性に関する経験的な長さスケールを用いて、通常の積分型光弾性法の方程式を丁寧に導出しており、本論文は任意の三次元的で弱い異方性をもつ誘電テンソル場の再構成について自己完結的に記述している。  

I. 序論

固定周波数における外部測定から媒質の異方的かつ空間的に変化する電磁特性を復元する境界値逆問題は、逆問題の中でも最も数学的に困難なものの一つである。静的近似が有効な低周波の電磁測定では、異方的な誘電テンソルまたは導電率テンソルは、あるゲージ条件を除いて完全な表面電気測定によって一意に決定されることが知られている[1]。等方性[2]およびカイラル等方性媒質[3]の場合、境界上の接線方向の電場と磁場の全ての組み合わせを単一の例外的でない周波数で知ることが、材料特性を復元するのに十分であることが知られている。異方性材料は、流れ、変形、結晶または液晶構造、あるいは繊維質や層状複合材料の均質化から生じる有効異方性特性など、多くの実際的な問題において重要である。一般の異方性媒質に対して、復元に必要なデータが十分であるかどうかの問題は依然として未解決である。 本論文では、実用的に重要性の高い特定の高周波ケースに焦点を当てる。ここでは、材料が非磁性、すなわち透磁率が真空値 µ0 と等しい一様かつ等方的な値を持ち、導電率は無視でき、誘電率すなわち誘電テンソルが弱い異方性を持つ場合を仮定する（その意味は後で定義する）。本論文では、光線が光学的異方性材料を様々な角度から通過する際に得られるトモグラフィ測定データから、（空間的に変化する）誘電テンソル εij の6つの独立成分を逆算するための数学的枠組みを提示する。弱い異方性の条件下では、光線が媒質を走査する平面を6つ適切に選んで測定すれば、本手法によって全てのテンソル成分を完全に決定できることを示す。 不均質かつ弱異方性媒質中を光が通過する際の方程式は文献[4–6]で定式化されており、第II章では、幾何光学的アプローチから出発し、媒質の不均質性と異方性を記述する適切な長さスケールで電場を展開することによりこれらの方程式を導出する方法を示す。また、「弱い」異方性の条件を定義し、その条件の下で光線に沿った偏光伝達行列を決定する方程式を線形化する（第IV章）。続いて第V章では、6つの異なる平面で試料を走査して得られる偏光変換データに対してスカラーRadon変換を行うことにより、線形化極限で誘電テンソルを再構成する手法を提示する。第VI章では、軸方向荷重を受ける円柱棒内部の誘電テンソルを再構成した視覚的な例を示す。この場合、透明媒質内の応力テンソルが光学的異方性を生じさせ、本手法が弱異方性の極限で有効に機能することが示される。本手法の数学的基盤はSharafutdinov[7]の先駆的な研究に見られ、彼の著書ではn次元ユークリッド空間におけるテンソル場の「光線変換」の一般理論と、それを逆変換して基礎となるテンソル場を再構成する可能性について議論されている。本研究は、この高度に数学的な枠組みを、不均質な誘電テンソル場を光学トモグラフィによって再構成するという具体的目的に適応・部分的に再定式化したものである。 さらに、本研究で提示する異方性誘電テンソルの決定法は、荷重を受けた透明材料内部の応力テンソル場を再構成する問題とも密接に関連している。もともと等方的であった媒質が歪みによって光学的異方性を示す現象は光弾性[8–12]と呼ばれ、トモグラフィ法で得られる偏光変換データから、歪んだ媒質内部の応力情報を得るために利用できる。光弾性効果を応力情報の取得に利用する研究は広く行われており、「積分型光弾性法」[13–17]として知られる方法がある。また、主応力成分の差については、偏光変換データに対する適切なRadon変換から得られることが指摘されている[18–21]。しかし、これらの方法では応力成分を個別に再構成することができず、そのため全応力テンソルおよび線形に関連する誘電テンソルは、軸対称性のような特定の対称性を持つ系でのみ完全に決定できる。他の再構成法として、三光束測定スキーム[22]（軸対称系に対してはオニオンピーリング法[23, 24]が提案されている）や、一般的な三次元誘電テンソルを決定できる可能性のある「負荷増分法」[25]などもある。後者では、対象にかかる応力を小刻みに増加させ、各段階で測定を行う。 本論文における新しい成果は二つある。第一に、物体を走査する特定の光線に対して偏光伝達行列を決定する方程式群を導出した点である。これらの方程式は積分型光弾性法[14]の標準方程式に関連するが、本研究では経験的な長さスケール[6]とそれに基づく近似をより厳密に展開し、物質媒質におけるマクスウェル方程式から通常これらの方程式が導出される過程を整理した（第II–IV章）。第二に、複雑に不均質な誘電テンソルであっても、「弱い」複屈折という条件（第V・VI章で定義）に従う限り、試料内部で再構成する新しい方法を提示した点である。

II. ヒューリスティックな長さスケール

材料は光学波長に対して非吸収性であり、さらに非磁性であると仮定する。すなわち、磁気透磁率は $\mu_{ij} = \mu_0 \delta_{ij}$、ここで $\mu_0$ は真空の磁気透磁率である。誘電率については、誘電テンソルが大域的な空間平均から「弱く」逸脱していると仮定する： $$\bar{\epsilon} = \frac{1}{\text{vol}B} \int_B d^3x \, \frac{1}{3}\text{tr}(\epsilon) , \tag{1}$$ ここで $B$ は体積 $\text{vol}B$ を持つ物体、$\text{tr}\epsilon = \sum_{i=1}^3 \epsilon_{ii}$ はトレース演算、$\epsilon = \epsilon(x)$ は誘電テンソルである。零次近似では、この材料は均質等方性とみなされ、式 (1) で定義されるスカラー誘電定数 $\bar{\epsilon}$ を持つ。一般に、このような均質等方的背景からの弱い逸脱は、ガラスや特定の樹脂が中程度の内部応力や外部荷重を受けた場合に成立する。スカラー定数 $\bar{\epsilon}$ は、後に式 (8) において無次元の異方性度を規定する際の基準量となる。 本研究に関連する実際的な問題においては、材料内の不均質を特徴づける長さスケールは、物体を通過する（単色）光の波長よりもはるかに大きい。したがって幾何光学近似の使用は正当化される。経験的に、次の2つの長さスケールが考えられる[6]。

1. 不均質性を特徴づけるスケール $l_0$ は次式で導入できる：

$$l_0 \, |\kappa \cdot \nabla \text{tr}\epsilon(x)| \sim |\text{tr}\epsilon(x)| , \tag{2}$$ ここで $\kappa$ は波の伝搬方向の単位ベクトル。

2. 任意の異方性媒質では、各伝搬方向ごとに優先する2つの偏光方向 $e_p(x), p=1,2$ が存在する[26–30]。これらの偏光方向の変化率を測るスケール $l_p$ は次で与えられる：

$$l_p \, |\kappa \cdot \nabla e_p(x)| \sim |e_p(x)| . \tag{3}$$ これらのスケールは、媒質を通過する単色光の平均波長 $\lambda$ と比較されるべきであり、 $$l_p, l_0 \gg \lambda$$ が満たされるとき、電場は以下の形で表せる： $$\mathbf{E}(x,t) = E(x) e^{i\phi(x) - i\omega t} , \tag{4–5}$$ ここで、エイコナル $\phi(x + \lambda s) = \phi(x) + \lambda s \cdot \nabla \phi(x) + O(\lambda/l_0)$ は波数ベクトル $\nabla \phi(x)$ を持つ局所平面波を記述し、振幅は $\lambda \kappa \cdot \nabla E \sim O(\lambda/l_p)E$ のオーダーで変化する。 Fuki, Kravtsov, Naida (FKN)[6] は、次のような無次元スケールを導入した： $$\alpha = \max \left( \frac{\lambda}{l_0}, \, \frac{\lambda}{l_p} \right). \tag{6}$$ 幾何光学の極限は $\alpha \ll 1$ で与えられる。 次に、異方性の明示的な尺度が必要である： $$A_{ij} = \frac{\epsilon_{ij} - \delta_{ij}\bar{\epsilon}}{\bar{\epsilon}}, \quad \beta \sim \max A_{ij}, \tag{7–8}$$ ここで $\beta$ は無次元異方性テンソル $A_{ij}$ の成分の大きさを特徴づける指標（例えば最大値）である。異方性が強い場合、任意の伝搬方向において2つの異なる位相速度と2つの異なる光線速度が存在するため、媒質内の各点で光線の連続的分裂が生じる。したがって、弱異方性の条件とは、光の伝搬が単一の光線によって記述でき、その光線は局所的な光学テンソルの変動によって偏光方向が回転する形でのみ影響を受ける、という状況を意味する。この状況は光弾性で典型的に現れ、我々の研究にとって最も重要である。 FKN[6] は、次の条件が満たされると光線分裂を無視できることを示した： $$\frac{\beta}{\alpha} \lesssim 1 . \tag{9}$$ この場合、「分裂した光線を実験的に区別することは不可能」であり、誘電テンソルの等方成分から得られる単一の等方光線で置き換えることができる。これは光弾性実験において光線分裂が観測されない領域であり、我々の関心の中心である。実際、対象とする応用では光は通常、試料を直線的に通過するので、式 (5) の試行解はさらに強い仮定で置き換えられる： $$\mathbf{E}(x,t) = E(x) e^{i k \cdot x - i \omega t} , \tag{10}$$ ここで $k$ は定数の波数ベクトルであり、平面波と同様である。この波数に関連する位相速度は、式 (1) で定義された平均誘電率 $\bar{\epsilon}$ によって決まる： $$k = \frac{\omega}{u}, \quad u = \frac{1}{\sqrt{\mu_0 \bar{\epsilon}}}, \quad \lambda = \frac{2\pi}{k}. \tag{11}$$ ただし、振幅 $E(x)$ は光線に沿った偏光方向の変化を反映して空間依存性を持つ。 条件 (7) および (9) の下で、電場 $E$ と電束密度 $D$ は次のように振る舞う： $$\kappa \cdot D \sim O\!\left(\frac{\lambda}{l_p}\right)D , \quad \kappa \cdot E \sim O\!\left(\frac{\lambda}{l_p}\right)E + \beta E . \tag{12}$$ つまり、$D$ はほぼ横波であり、$E$ についても異方性が十分小さい場合に限り同様のことが成り立つ： $$\beta \ll 1 .$$ これが「弱い異方性」の条件であり、本研究の手法はこの「準等方的」領域を前提として定式化されている。  

III. 伝達行列が満たす方程式

材料を通過する光線の偏光状態の変化に関する情報は、2次元のユニタリ伝達行列に符号化される。光線に沿った伝達行列が満たす方程式は、光弾性トモグラフィに関する多くの研究[18–24]で示されているが、$\lambda / l_0$ や $\lambda / l_p$ の比の高次項を無視する過程での近似は必ずしも明示されていない。そこで、ここでは必要なステップを簡単にまとめる。 式 (10) をマクスウェル方程式に代入すると、 $$k \times (k \times E) + \omega^2 \mu_0 \epsilon E - i \nabla \times (k \times E) + k \times (\nabla \times E) - \nabla \times \nabla \times E = 0 \tag{14}$$ を得る。 次を示すのは容易である： $$\nabla \times \nabla \times E \sim O\!\left( \left(\frac{\lambda}{l_p}\right)^2 \right) E . \tag{15}$$ したがって、幾何光学的極限 (7) ではこの項を無視できる。すると式 (14) は次の形になる： $$\kappa \times (\kappa \times E) - \frac{i}{k} \Big[ \nabla \times (\kappa \times E) + \kappa \times (\nabla \times E) \Big] + \mu_0 u^2 \epsilon E = 0 , \tag{16}$$ ここで $u$ は式 (11) で与えられた量である。式 (16) の縦成分（光線の伝搬方向 $\kappa = k/k$ に投影した成分）は、 $$O\!\left(\frac{\lambda}{l_p}\right) E \tag{17}$$ のオーダーであり、幾何光学的極限では無視される。したがって、保持するのは波の伝搬方向 $\kappa$ に垂直な $E$ と $D$ の横成分のみである。 次に、$\kappa$ が $z$ 軸に沿っている座標系で式 (16) を考えると、 $$\frac{d}{dz} \begin{pmatrix} E_1 \\ E_2 \end{pmatrix} = i \frac{\pi}{\lambda} \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} \begin{pmatrix} E_1 \\ E_2 \end{pmatrix} , \tag{18}$$ が得られる。ここで $A_{ij}$ は式 (8a) で定義されたものである。式 (18) の解は伝達行列を用いて表すことができる： $$\begin{pmatrix} E_1(z) \\ E_2(z) \end{pmatrix} = U(z, z_0) \begin{pmatrix} E_1(z_0) \\ E_2(z_0) \end{pmatrix} , \tag{19}$$ ここで $U$ は次の常微分方程式を満たす： $$\frac{d}{dz} U(z, z_0) = i \frac{\pi}{\lambda} A_\perp(z) U(z, z_0), \quad U(z_0, z_0) = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} . \tag{20}$$ 式 (20) は積分方程式の形でも書ける： $$U(z, z_0) = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} + i \frac{\pi}{\lambda} \int_{z_0}^z dz_1 \, A_\perp(z_1) U(z_1, z_0) , \tag{21}$$ ここで $A_\perp$ は式 (18) に現れる $A_{ij}$ の横成分の行列である。この形式的解はボルン＝ノイマン展開によって与えられる： $$\begin{aligned} U(z, z_0) &= \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} + i \frac{\pi}{\lambda} \int_{z_0}^z dz_1 A_\perp(z_1) \\ &\quad + \left( i \frac{\pi}{\lambda} \right)^2 \int_{z_0}^z dz_1 A_\perp(z_1) \int_{z_0}^{z_1} dz_2 A_\perp(z_2) + \cdots . \end{aligned} \tag{22}$$ 伝達行列 $U$ はユニタリであり、したがって複素電場ベクトルのノルムを保存する。物理的にはこれは強度が保存されることを意味し、ユニタリ性は光線のエネルギー保存を表している。実際、ここでは媒質を非吸収性と仮定しているため、この条件は当然満たされる。

IV. 線形化された逆問題

媒質を走査する十分多くの光線について伝達行列 $U$ が測定できたと仮定する。このとき、関連する逆問題は、これらの伝達行列の集合から異方性テンソル $A_{ij}$ を再構成することにある。この逆問題は式 (21) および (22) から明らかなように $A_{ij}$ に対して非線形である。完全な非線形問題の解はまだ知られていない。しかし、準等方的領域では線形化された逆問題を扱うことができる。これは、式 (22) のボルン＝ノイマン級数を一次項で打ち切ったものとして定義される： $$U(z,z_0) = I_2 + i \frac{\pi}{\lambda} \int_{z_0}^z dz_1 \, A_\perp(z_1) . \tag{23}$$ 例えば、相対異方性 $\nu \sim 10^{-9}$、波長 $\lambda \sim 0.5 \times 10^{-6}\,\text{m}$、試料長さ $L \sim 1\,\text{m}$ を仮定すると、式 (23) の一次項は $10^{-3}$ 程度であり、この場合、線形化は良い近似となる。 伝達行列 $U$ は、媒質を様々な角度で通過する光線に沿った偏光変化の適切な測定によって決定されなければならない。その際、干渉計測なども利用できる。いわゆる3つの特性パラメータ[13]を測定することで、ポアンカレ等価定理[31]を用いて伝達行列 $U$ の SU(2) 成分を計算できる。この定理は行列分解定理であり、特性パラメータを、遅相角 $\Delta$、速軸の角度 $\theta$、回転角 $\delta$ として解釈でき、これらを持つ直線リターダーとローテータからなる等価光学モデルに対応させることができる。ポアンカレ等価定理は、ジョーンズ行列やポアンカレ球上のストークスパラメータによって定式化でき、最近の解説は文献[32]に与えられている。 したがって、$\Delta,\theta,\delta$ の測定から、次を満たすユニモジュラ行列 $S(\Delta,\theta,\delta)$ が決定される： $$U = S(\Delta,\theta,\delta) \, e^{i\Phi}, \quad S(\Delta,\theta,\delta) \in SU(2) , \tag{24}$$ ここで $\Phi$ は伝達行列 $U$ のグローバル位相である。一般の場合、異方性に制限がないとき、この位相 $\Phi$ は任意に大きくなり得て、特性パラメータだけからは決定できないため、各光線ごとに干渉計測などで測定する必要がある。しかし弱異方性の極限では、グローバル位相は特性パラメータだけで実質的に決定される。すなわち、$\Delta,\theta,\delta$ を測定して $S$ を求めると、式 (23) から、未知の位相 $\Phi$ は次の条件を満たすように選ばれるべきである： $$\Re \left(e^{i\Phi} S_{11}\right) = \Re \left(e^{i\Phi} S_{22}\right) = 1, \quad \Re \left(e^{i\Phi} S_{12}\right) = \Re \left(e^{i\Phi} S_{21}\right) = 0 . \tag{25}$$ したがって、$\Phi$ はこれらの方程式のいずれか、または数値的安定性のためにはすべてから平均値をとって決定できる。したがって弱異方性の極限では、伝達行列に含まれる実数の自由度は一般の場合の4ではなく、実質的に3となる。

---

V. 六つのスカラー Radon 逆変換による線形化逆問題の解法

線形化された逆問題は、偏光変換データに対して6つのスカラー Radon 逆変換を行うことで解けることを示す。まず、試料を横切る点 $y$ を含む平面 $P(y,\eta)$ を指定する。この平面の向きは、平面の法線である単位ベクトル $\eta$ によって決まる。単位ベクトル $\kappa$ が $P(y,\eta)$ 内にあり、光線 $t \mapsto y + t\kappa$ が試料を通過するとする。さらに、$\kappa$ と $\eta$ に直交する3つ目の単位ベクトル $\xi$ を導入し、$(\xi,\eta,\kappa)$ が右手系をなすとする。このとき光線に沿った偏光変換は、式 (23) に類似する形で次のように表される： $$U(t,t_0) = I_2 + i\frac{\pi}{\lambda} \int_{t_0}^t dt_1 \begin{pmatrix} A_{\xi\xi}(t_1) & A_{\xi\eta}(t_1) \\ A_{\eta\xi}(t_1) & A_{\eta\eta}(t_1) \end{pmatrix}. \tag{26}$$ 左辺の伝達行列は特性パラメータの測定から決定される。これを平面 $P(y,\eta)$ 内の全ての光線について繰り返すと、法線成分 $A_{\eta\eta}$ に関する線積分の集合が得られる： $$\int dt_1 \, A_{\eta\eta}(t_1) = \frac{\lambda}{i\pi} \left[ U_{\eta\eta}(+\infty,-\infty) - 1 \right]. \tag{27}$$ この積分集合は、$\eta$ に関する $A_{ij}$ の横光線変換[7]と呼ばれる。特に $A_{\eta\eta}$ は平面 $P$ に垂直であり、SO(2) 回転に対してスカラー関数のように振る舞う。したがって (27) の集合はスカラー関数 $A_{\eta\eta}(x)$ の標準的な2次元 Radon 変換であり、数値的な逆 Radon 変換手法（例えばフィルタ付き逆投影[34]）によって復元可能である。これにより、平面 $P(y,\eta)$ 上の全点で $A_{\eta\eta}(x)$ が再構成される。 この手続きを平行な全ての平面に繰り返すことで、試料全体にわたり $A_{\eta\eta}$ が復元される。 さらに、次の6つの $\eta$ の選択について同様の手続きを行う： $$\eta_1 = e_1, \quad \eta_2 = e_2, \quad \eta_3 = e_3 , \tag{28a}$$ $$\eta_{12} = \tfrac{1}{\sqrt{2}} (e_1 + e_2), \quad \eta_{23} = \tfrac{1}{\sqrt{2}} (e_2 + e_3), \quad \eta_{31} = \tfrac{1}{\sqrt{2}} (e_3 + e_1). \tag{28b}$$ 式 (28a) に対応する走査では $A_{11}, A_{22}, A_{33}$ が得られる。一方、例えば $\eta_{12}$ の場合、再構成されるのは $A(\eta_{12},\eta_{12})$ であり、テンソルの線形性と対称性から： $$A(\eta_{12},\eta_{12}) = \tfrac{1}{2}(A_{11} + A_{22}) + A_{12} . \tag{30}$$ したがって、既に $A_{11}, A_{22}$ を得ていれば、この式から $A_{12}$ を計算できる。残り2つの $\eta$ についても同様にすれば、全6成分 $A_{ij}$ が復元される。 もし平均誘電定数 $\bar{\epsilon}$ が既知であれば、定義式 (8a) を用いて完全な誘電テンソル $\epsilon_{ij}$ を直ちに求められる。未知の場合でも、 $$\epsilon_{ij} = \bar{\epsilon}(A_{ij} + \delta_{ij}) \tag{31}$$ と書けるため、$\bar{\epsilon}$ のスケール因子までで $\epsilon_{ij}$ を復元できる。これは依然として媒質内部構造に関する有用な情報を提供する。

---

VI. 数値例

ここでは、物体を斜めに横切る平面におけるテンソル成分 $A_{\eta\eta}$ の再構成例を示す。偏光変換データは、円柱棒が軸方向荷重を受けて中央部が膨らむ人工的な応力モデルから得られる（図1(a)）。このデータを用い、本手法で円柱の中心を通り対称軸と $22^\circ$ の角度をなす平面における $A_{\eta\eta}$ を再構成した（図1(b)）。元の $A_{\eta\eta}$ の分布（図1(c)）と比較すると、再構成画像（図1(d), 254×254ピクセル）はラドン逆変換に典型的なアーティファクトを含むが、走査回数を増やせば改善される。例えば、1°ごとに計180回走査すれば、図1(c) とほとんど区別がつかない結果となる。 逆 Radon 変換には、フィルタ付き逆投影アルゴリズムを用いた Matlab 関数 iradon を使用した[34]。

---

VII. まとめ

本論文では、透明で非吸収性の媒質内部において、均質等方的平均からの弱い逸脱として現れる複屈折条件の下で、任意に不均質な異方性誘電テンソルを再構成する新しい方法を提示した。光学トモグラフィによって得られる偏光変換データから誘電テンソルを復元する線形化逆問題は、6回のスカラー Radon 逆変換に帰着でき、これにより誘電テンソルを完全に決定できることを示した。また、誘電テンソルの逆問題を定義する従来の積分型光弾性の方程式について、マクスウェル方程式からの導出に含まれる各種近似を詳細に記述し、より厳密な定式化を与えた。