干渉SAR(InSAR)の原理 — 位相差から地表変位を測る

2016年の熊本地震のあと、国土地理院は「地面が最大2メートル動いた」という変位図を数日で公開しました。地上に測量機材を置いていない山の中まで、面的に、しかもセンチメートル精度で。どうやって測ったのでしょうか。答えは、高度700〜800 kmを飛ぶ人工衛星が地震の前と後に撮った2枚のレーダー画像を「引き算」しただけ、です。使ったものさしは、レーダーの波長そのもの — わずか5.5 cmの電波の波です。

この「電波の波を定規にして、地面が何ミリ動いたかを宇宙から測る」技術が干渉SAR（InSAR: Interferometric Synthetic Aperture Radar）です。光学カメラで撮った写真をどれだけ拡大しても、地面が1 cm沈んだかどうかはわかりません。ところがSARが記録している位相という情報を使うと、電波の波長の100分の1、つまりサブミリメートルの世界にまで手が届きます。

干渉SARを理解すると、次のような分野への見通しが一気に開けます。

• 地震・火山モニタリング: 断層のすべり分布や火山のマグマだまりの膨張を、地上観測点のない場所でも面的に捉えられます。1993年に Nature に載ったランダース地震の干渉画像は、地球物理学の観測手段を根本から変えました
• 地盤沈下・インフラ監視: 地下水汲み上げによる都市の沈下、トンネル掘削の影響、堤防やダム・橋梁のミリ単位の変形を、定期的に、広域に監視できます
• 地形計測（DEM生成）: 位相差は高さにも敏感なので、標高モデルを作れます。2000年のスペースシャトルSRTMミッションは、この原理で地球の陸地の8割の地形図を11日間で作りました
• 氷河・氷床の流動速度: 極域の氷がどれだけ速く海に流れ込んでいるかを、氷の上に立たずに測れます

まず、干渉SARが何をしているのかを1枚の図にまとめておきます。

左の (a) が測定の舞台です。衛星は同じ場所を2回（たとえば12日あけて）撮影しますが、その間に地面が動くと往復距離が数mm〜数cmだけ変わります。中央の (b) が仕掛けで、2枚の複素画像 $s_1, s_2$ の複素共役積を取ると位相が引き算になり、$\Delta\varphi = \frac{4\pi}{\lambda}\Delta R$ という単純な関係だけが残ります。右の (c) が出力で、位相を色で表すと同心円状の縞（フリンジ）が現れ、縞1本が視線方向に半波長（Cバンドで2.77 cm）の変位に対応します。この記事は、この3枚の絵を数式で厳密に埋めていく作業だと思ってください。

一方で、干渉SARには独特の「クセ」があります。測れるのは視線方向の変位だけ、位相は $2\pi$ ごとに巻き戻ってしまう、地面の状態が変わると位相が意味を失う、大気の揺らぎが変位に化ける — こうした落とし穴をひとつずつ潰していくのが、この技術の学びどころです。

本記事の内容

• SAR画像の複素数表現と、干渉位相 $\Delta\varphi = \frac{4\pi}{\lambda}\Delta R$ の導出
• 基線ベクトルの幾何と、平行基線 $B_\parallel$・垂直基線 $B_\perp$ の定義
• 地形位相 $\varphi_{\text{topo}} = \frac{4\pi B_\perp h}{\lambda R \sin\theta}$ の導出と高度あいまい度 $h_a$
• 変位位相 $\varphi_{\text{disp}} = \frac{4\pi}{\lambda} d_{\text{LOS}}$ と視線方向への投影
• 干渉位相の分解（平地・地形・変位・大気・軌道・ノイズ）
• コヒーレンス $\gamma$ の定義と、時間・幾何・ボリューム非相関の3要因
• コヒーレンスと位相標準偏差を結ぶクラメール・ラオ型の関係式とその適用限界
• 位相アンラップ — 伊藤の条件、残差（residue）、枝切り法・最小コスト流れ
• Goldstein フィルタの原理
• Python実装: 合成干渉図の生成 → ノイズ注入 → フィルタ → アンラップ → 変位マップ復元と誤差評価

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 合成開口レーダー（SAR）の原理
• SARのレンジドップラーアルゴリズム
• 電磁波の伝搬 — マクスウェル方程式から波動方程式を導く

干渉SARとは — 波長を「ものさし」にする

巻尺で長さを測るとき、精度は目盛りの細かさで決まります。1 mm刻みの巻尺なら1 mmまで、1 cm刻みなら1 cmまでです。レーダーが持っている「目盛り」は何でしょうか。

普通に考えると、レーダーの距離分解能です。パルスの帯域幅 $W$ に対して、スラントレンジ分解能は $\rho_r = c/(2W)$ で、帯域50 MHzなら約3 mです。3 mの目盛りでは、地面が1 cm動いても検出できません。

ところがレーダーは、反射波の振幅だけでなく位相も記録しています。位相は、電波が往復した距離を波長で測った「小数点以下の端数」に相当します。波長5.5 cmのCバンドなら、位相が $2\pi$ 進むのは往復距離が5.5 cm変わったとき、つまり片道で2.775 cm変わったときです。位相を1度の精度で読めれば、$2.775 \text{ cm} \times (1/360) \approx 0.08$ mm の分解能になります。3 mの目盛りが、いきなり0.1 mmの目盛りに化けるわけです。

問題は、1枚のSAR画像の位相をそのまま読んでも何もわからないことです。1つの分解能セル（数m四方）の中には、小石も草も建物の壁も、無数の散乱体が入っています。受信信号はそれらからの反射の複素和なので、位相は事実上ランダムな値になります。これがスペックルの正体です。

干渉SARの発想は単純です。同じ場所を2回撮って、位相を引き算する。ランダムな散乱位相は2回とも同じなので、引き算すれば消えます。残るのは「2回の観測で往復距離がどれだけ変わったか」という情報だけです。

この「引き算で消える」という主張を、複素平面の図で確かめておきます。

(a) の灰色の折れ線が、セル内の個々の散乱体からの反射を順に足し合わせた様子で、その総和（青い太矢印）の向きは事実上ランダムです。(b) は2回目の観測ですが、散乱体の配置が変わっていないので、和のかたち自体は同じで、全体が $\Delta\varphi$ だけ回転しているだけです。(c) で2本を重ねると、ランダムだった散乱位相は共通なので消え、2本のなす角 $\Delta\varphi$ — つまり往復距離の変化だけが残ることがわかります。逆に言えば、2回の間に散乱体の配置が入れ替わると (b) の前提が崩れ、この角度は意味を失います。

日常のアナロジーで言えば、こうです。二人が同じ長さの物差しを持って、それぞれ違う場所から同じ壁までの距離を測ります。物差しの目盛りの読み方に共通の癖があっても、二人の測定値の差を取れば、癖は打ち消され、二人の位置の違いだけが残ります。干渉SARの「引き算」も、まったく同じ理屈です。

この「引き算で散乱体のランダム位相が消える」という前提こそが干渉SARの生命線であり、後で説明するコヒーレンスは、まさにこの前提がどれだけ成り立っているかを測る量になります。

では、この引き算を数式で厳密に書き下すところから始めましょう。

SAR画像の複素数表現と干渉位相

複素画像（SLC）

SARの生データを画像化した結果は、画素ごとに複素数を持つSLC（Single Look Complex）画像です。ある画素の値を

$$ s = A e^{j\varphi} $$

と書きます。$A$ は反射の強さ（後方散乱の振幅）、$\varphi$ は位相です。位相は次の2つの寄与からできています。

$$ \varphi = -\frac{4\pi R}{\lambda} + \varphi_{\text{scat}} $$

第1項が伝搬位相です。電波は衛星から地表まで距離 $R$ を進み、反射して同じ距離を戻ってくるので、往復の光路長は $2R$ です。波長 $\lambda$ の波が距離 $2R$ を進むと位相は $2\pi \cdot 2R/\lambda = 4\pi R/\lambda$ だけ遅れます。符号が負なのは、$e^{j\omega t}$ の時間規約で「遅れ」を表すためです。

第2項 $\varphi_{\text{scat}}$ が散乱位相で、分解能セル内の散乱体の配置で決まる、実質的にランダムな値です。ここが厄介です。$\varphi$ を1枚だけ見ても、$-4\pi R/\lambda$ と $\varphi_{\text{scat}}$ が混ざっていて、$R$ を取り出せません。

干渉図の定義

同じ地点を撮った2枚のSLC画像 $s_1, s_2$（幾何位置を精密に合わせる「レジストレーション」は済んでいるとします）に対し、干渉図（インターフェログラム）を次で定義します。

$$ \begin{equation} I = s_1 s_2^{*} = A_1 A_2 \, e^{j(\varphi_1 – \varphi_2)} \end{equation} $$

複素共役を掛けると、位相は引き算になります。ここで、2回の観測の間に分解能セル内の散乱体の配置が変わっていなければ $\varphi_{\text{scat},1} = \varphi_{\text{scat},2}$ なので、この項は完全に消えます。残るのは、

$$ \Delta\varphi = \varphi_1 – \varphi_2 = -\frac{4\pi R_1}{\lambda} + \frac{4\pi R_2}{\lambda} = \frac{4\pi}{\lambda}(R_2 – R_1) $$

です。$\Delta R \equiv R_2 – R_1$ とおけば、干渉SARの出発点となる関係式が得られます。

$$ \begin{equation} \boxed{\; \Delta\varphi = \frac{4\pi}{\lambda}\Delta R \;} \end{equation} $$

係数 $4\pi/\lambda$ の $4$ は往復（2倍）から来ています。片道の伝搬なら $2\pi/\lambda$ です。この因子2が、干渉SARの感度を2倍にしてくれています。

ラッピング — $2\pi$ の壁

ここで最初の落とし穴に出会います。位相は角度なので、実際に観測できるのは $[-\pi, \pi)$ に折り畳まれた値だけです。ラッピング演算子を

$$ \mathcal{W}\{\varphi\} = \varphi – 2\pi \left\lfloor \frac{\varphi + \pi}{2\pi} \right\rfloor $$

と定義すると、干渉図から読み取れるのは $\psi = \mathcal{W}\{\Delta\varphi\}$ であって、$\Delta\varphi$ そのものではありません。真の値との関係は

$$ \Delta\varphi = \psi + 2\pi k, \qquad k \in \mathbb{Z} $$

で、整数 $k$（巻き数）は未知です。干渉図をカラー表示すると縞模様（フリンジ）が現れますが、この1本の縞が $2\pi$ の巻き戻し1回に対応します。$\Delta\varphi = 2\pi$ に対応する視線方向距離の変化は

$$ \Delta R = \frac{\lambda}{4\pi} \cdot 2\pi = \frac{\lambda}{2} $$

つまり 1フリンジ = 半波長 です。Cバンド（$\lambda = 5.55$ cm）なら1フリンジ2.77 cm、Xバンド（$\lambda = 3.1$ cm）なら1.55 cm、Lバンド（$\lambda \approx 23.6$ cm）なら11.8 cmです。干渉画像を見て縞を数えれば、その半波長単位で変位がわかる — これが干渉画像の読み方の基本です。

(a) の青い直線が真の位相 $\Delta\varphi$、赤い鋸歯が実際に観測できるラップ位相 $\psi$ です。両者の差はどこでも $2\pi$ の整数倍で、鋸歯の1本1本が巻き数 $k$ の1つの増分に対応します。(b) は同じ鋸歯を色で表したもので、これが干渉図の縞そのものです — 縞を数えるという行為は、赤い鋸歯が何回リセットしたかを数えることに等しいわけです。(c) は1フリンジが表す変位量を波長ごとに並べたもので、Xバンドの1.55 cmに対しLバンドは11.8 cmと8倍近く粗くなります。細かく測りたいなら短波長ですが、そのぶん同じ変位で縞が密になり、次に述べるアンラップが難しくなるというトレードオフがあります。

そして、この未知の整数 $k$ を画像全体で矛盾なく決める作業が位相アンラップであり、干渉SAR処理で最も難しい工程です。後半でじっくり扱います。

ここまでで「位相差 $\Delta\varphi$ は視線距離の変化 $\Delta R$ を半波長単位で表したもの」ということがわかりました。では、$\Delta R$ には何が含まれているのでしょうか。次は幾何を詰めて、$\Delta R$ を成分に分解します。

干渉SARの幾何 — 基線ベクトル

基線の定義

2回の観測は、まったく同じ軌道位置からは行われません。繰り返し観測（リピートパス）では、12日後、24日後に同じ軌道を通っても、実際の位置は数十〜数百メートルずれます。この2つのアンテナ位置の差ベクトルを基線（baseline） $\bm{B}$ と呼びます。

基線には2種類あります。SRTMのように1機のプラットフォームに2つのアンテナを積むシングルパス方式では、$\bm{B}$ は機械的に固定された長さ（SRTMでは60 mのマスト）で、両者は同時刻に観測します。一方、Sentinel-1やALOS-2のようなリピートパス方式では、$\bm{B}$ は軌道のばらつきで決まり、時間差も生まれます。地表変位を測りたいなら、時間差があるリピートパスが必須です。

座標系を決めましょう。ターゲット $P$ を見込む面（レンジ方向の鉛直断面）を考え、水平方向でターゲットに向かう向きを $\hat{\bm{g}}$、鉛直上向きを $\hat{\bm{u}}$ とします。アンテナ1からターゲットを見る視線（LOS: Line Of Sight）の単位ベクトルは、天底からの視線角（look angle） $\theta$ を使って

$$ \hat{\bm{p}} = \sin\theta \, \hat{\bm{g}} – \cos\theta \, \hat{\bm{u}} $$

と書けます（衛星から地表へ向かう向き）。基線ベクトルは、水平からの傾き $\alpha$ を使って

$$ \bm{B} = B\left(\cos\alpha \, \hat{\bm{g}} + \sin\alpha \, \hat{\bm{u}}\right) $$

とします。

$\Delta R$ の導出

アンテナ1の位置を $S_1$、アンテナ2を $S_2 = S_1 + \bm{B}$、ターゲットを $P = S_1 + R_1 \hat{\bm{p}}$ とします。$S_2$ からターゲットまでの距離は

$$ R_2 = \left| P – S_2 \right| = \left| R_1 \hat{\bm{p}} – \bm{B} \right| $$

です。中身を展開します。

$$ R_2 = \sqrt{R_1^2 – 2 R_1 \, \hat{\bm{p}} \cdot \bm{B} + B^2} $$

ここで $B \ll R_1$（数百メートル対数百キロメートル、比にして $10^{-3}$ 程度）なので、$R_1$ でくくって平方根を展開します。

$$ R_2 = R_1 \sqrt{1 – \frac{2\,\hat{\bm{p}}\cdot\bm{B}}{R_1} + \frac{B^2}{R_1^2}} $$

$\sqrt{1+x} \approx 1 + x/2$ を使い、$B^2/R_1^2$ の項は $10^{-6}$ オーダーなので落とすと、

$$ R_2 \approx R_1 – \hat{\bm{p}} \cdot \bm{B} $$

内積を計算しましょう。

$$ \hat{\bm{p}} \cdot \bm{B} = B\left(\sin\theta\cos\alpha – \cos\theta\sin\alpha\right) = B\sin(\theta – \alpha) $$

加法定理でひとまとめになりました。この量を平行基線（parallel baseline）と呼びます。

$$ \begin{equation} B_\parallel = B\sin(\theta – \alpha) \end{equation} $$

したがって、

$$ \begin{equation} \Delta R = R_2 – R_1 = -B_\parallel = -B\sin(\theta – \alpha) \end{equation} $$

平行基線は、基線ベクトルを視線方向に射影した成分そのものです。視線方向に2つのアンテナがずれていれば、その分だけ距離が違う — 素直な結果です。

図では基線 $B$ を実スケールより大きく誇張して描いています（実際は $B/R \sim 10^{-3}$）。緑の矢印が基線 $\bm{B}$ で、水平から $\alpha$ だけ傾いています。これをオレンジの成分（視線方向）と紫の成分（視線に直交）に分けたものが $B_\parallel$ と $B_\perp$ です。オレンジの長さがそのまま距離差 $|\Delta R|$ になり、紫の長さは距離差が視線角 $\theta$ に対してどれだけ敏感かを決めます。$\theta – \alpha$ が小さいほどオレンジが短く紫が長い、つまり距離差は小さいのに地形には敏感、という関係になっている点に注目してください。

垂直基線 $B_\perp$ が主役になる理由

$B_\parallel$ が距離差を作るのに、なぜ干渉SARの教科書では $B_\perp$ ばかり出てくるのでしょうか。理由は、私たちが知りたいのが「$\Delta R$ の絶対値」ではなく「$\Delta R$ がターゲットの位置によってどう変わるか」だからです。

$\Delta R = -B\sin(\theta-\alpha)$ を視線角 $\theta$ で微分してみます。

$$ \frac{\partial (\Delta R)}{\partial \theta} = -B\cos(\theta – \alpha) \equiv -B_\perp $$

ここで登場した

$$ \begin{equation} B_\perp = B\cos(\theta – \alpha) \end{equation} $$

が垂直基線（perpendicular baseline）、すなわち基線ベクトルの視線に直交する成分です。つまり $B_\perp$ は、視線角のわずかな違いが距離差にどれだけ効くかを決める感度係数です。地形の高低は視線角の違いとして現れるので、$B_\perp$ が地形への感度を支配します。ここが干渉SARの核心です。

平坦な地面を仮定しても、視線角 $\theta$ はレンジ方向（近距離側から遠距離側へ）にゆっくり変化します。したがって $\Delta R$ も変化し、平坦地でも縞が現れます。これを平地位相（flat-earth phase） $\varphi_{\text{flat}}$ と呼び、軌道情報から解析的に計算して差し引きます。この操作をフラットニングと言い、これをやらないと、地形も変位も細かい縞に埋もれて見えません。

これで幾何の道具が揃いました。次は、フラットニング後に残る「地形位相」を定量化します。

地形位相と高度あいまい度

高さから位相へ

平地からの高さ $h$ を持つターゲットは、同じスラントレンジ $R$ で見たとき、平地上の点とは違う視線角で見えます。この視線角のずれが $\Delta R$ を変え、位相を変えます。順を追いましょう。

衛星の（基準面からの）高度を $H$ とすると、スラントレンジ $R$、視線角 $\theta$ のターゲットの高さは

$$ z = H – R\cos\theta $$

です。ここでスラントレンジ $R$ を固定したまま視線角を $\delta\theta$ 動かすと、高さは

$$ \delta z = R \sin\theta \, \delta\theta $$

だけ変わります。逆に言えば、高さ $h$ の違いは視線角の違い

$$ \delta\theta = \frac{h}{R\sin\theta} $$

に対応します。この関係が「高さ → 視線角」の変換です。次に「視線角 → 距離差」の変換は、前節で求めた $\partial(\Delta R)/\partial\theta = -B_\perp$ です。2つを繋ぎます。

$$ \delta(\Delta R) = \frac{\partial (\Delta R)}{\partial \theta}\,\delta\theta = -\frac{B_\perp h}{R\sin\theta} $$

最後に、$\Delta\varphi = (4\pi/\lambda)\Delta R$ で位相に直します。

$$ \varphi_{\text{topo}} = -\frac{4\pi}{\lambda}\cdot\frac{B_\perp h}{R\sin\theta} $$

符号は、どちらの画像を基準（マスタ）に取るか、$B_\perp$ をどちら向きに測るかで反転します。以後は符号を $B_\perp$ の向きに吸収させ、次の形を使います。

$$ \begin{equation} \boxed{\; \varphi_{\text{topo}} = \frac{4\pi B_\perp h}{\lambda R \sin\theta} \;} \end{equation} $$

(a) が「高さ → 視線角 → 距離差」という2段の変換を1枚にしたものです。緑の破線は衛星から等距離（同じスラントレンジ $R$）の円弧で、その上に載っている平地の点（青）と丘の上の点（赤）は、レーダーから見て同じレンジに写るのに視線角が $\delta\theta$ だけ違います。この $\delta\theta$ が $-B_\perp$ 倍されて距離差になり、さらに $4\pi/\lambda$ 倍されて位相になる、というのが $\varphi_{\text{topo}}$ の正体です。(b)(c) は実際の標高断面と、それを $B_\perp = 100$ mで位相に変換してラップしたものです。標高が $h_a = 139.7$ m変わるたびに位相が1周（縞が1本）していることが、鋸歯の数と点線の間隔の対応から読み取れます。

高さ $h$ に比例し、$B_\perp$ に比例する — この2点が本質です。$B_\perp = 0$ なら地形位相はゼロになり、干渉図は地形に対して完全に鈍感になります。変位だけを測りたいなら $B_\perp$ が小さいペアが理想的だ、という実務上の指針がここから出ます。

高度あいまい度 $h_a$

地形位相がちょうど $2\pi$（1フリンジ）変化する高さを高度あいまい度（height of ambiguity） $h_a$ と呼びます。定義から

$$ \frac{4\pi B_\perp h_a}{\lambda R\sin\theta} = 2\pi $$

これを $h_a$ について解くと、

$$ \begin{equation} \boxed{\; h_a = \frac{\lambda R \sin\theta}{2 B_\perp} \;} \end{equation} $$

$h_a$ は「干渉図の縞1本が標高何メートルに相当するか」を表す、干渉ペアの性格を一言で表すパラメータです。関連して、位相の高さに対する感度を垂直波数（vertical wavenumber）

$$ k_z = \frac{\partial \varphi_{\text{topo}}}{\partial h} = \frac{4\pi B_\perp}{\lambda R\sin\theta} = \frac{2\pi}{h_a} $$

と書くこともあります。

Sentinel-1 相当の条件（$\lambda = 5.55$ cm、$R = 800$ km、$\theta = 39^\circ$）で $h_a$ を並べてみましょう。

読み方はこうです。$B_\perp = 100$ mのペアなら、標高が140 m変わるごとに縞が1本現れます。$B_\perp = 600$ mなら23 mごとに1本 — 山岳地では縞が密集しすぎて、隣の画素との位相差が $\pi$ を超え、アンラップが破綻します。逆に、標高を精度よく測りたい（DEMを作りたい）なら $h_a$ は小さいほうがよく、長い基線が望ましい。変位計測では短基線、地形計測では長基線という設計指針が、この1つの式から出てきます。

表の最右列「平地の縞間隔」は、平坦地でフラットニングをしなかった場合に地上距離何メートルごとに縞が出るかで、$h_a/\tan\theta$ に等しくなります（$h_a = 139.7$ mなら $139.7/0.8098 = 172.5$ m）。

(a) は表を両対数プロットにしたもので、$h_a \propto 1/B_\perp$ がきれいな直線になっています。赤く塗った $h_a < 30$ mの帯に入ると、標高差30 mごとに縞が1本という密度になり、山岳地では隣接画素の位相差が $\pi$ を超えてアンラップが原理的に破綻します。(b) は同じ話を変位側から見たもので、$B_\perp = 600$ m（$h_a = 23$ m）のペアではDEM誤差わずか5 mが約6 mmもの偽の変位を生むのに対し、$B_\perp = 50$ m（$h_a = 279$ m）なら30 m誤っても3 mmに収まります。破線は本記事のシミュレーションで得られる雑音由来のRMSE（0.91 mm）で、長基線ではDEM誤差のほうが支配的な誤差源になることが一目でわかります。

DEM誤差が変位に化ける

変位を測るとき、私たちは外部のDEM（数値標高モデル）を使って $\varphi_{\text{topo}}$ を計算し、差し引きます。この操作を差分干渉SAR（DInSAR: Differential InSAR）と呼びます。しかしDEMには誤差 $\delta h$ があります。残ってしまう位相は $k_z \delta h$ で、これを変位に換算すると

$$ \delta d = \frac{\lambda}{4\pi} k_z \, \delta h = \frac{\delta h}{h_a}\cdot\frac{\lambda}{2} $$

です。$h_a = 139.7$ mのペアでDEM誤差10 mなら、$\delta d = (10/139.7)\times 2.775 \text{ cm} = 2.0$ mm の偽の変位が出ます。DEM誤差20 mなら4.0 mmです。ミリ単位の変位を議論するとき、これは無視できません。$B_\perp$ が小さいペアを選ぶと $h_a$ が大きくなり、この誤差が自動的に小さくなる — これも短基線が好まれる理由です。

地形の寄与を式で押さえられたので、次はいよいよ本命の変位位相に進みます。

変位位相 — ミリメートルを測る

視線方向変位と位相

2回の観測の間に地面が動いたとしましょう。1回目の観測時のスラントレンジを $R_1$、2回目を $R_2$ とし、その差を

$$ d_{\text{LOS}} = R_2 – R_1 $$

と定義します。$d_{\text{LOS}} > 0$ は「衛星から遠ざかった」を意味します。$\Delta\varphi = (4\pi/\lambda)\Delta R$ にそのまま入れれば、

$$ \begin{equation} \boxed{\; \varphi_{\text{disp}} = \frac{4\pi}{\lambda} d_{\text{LOS}} \;} \end{equation} $$

驚くほど単純です。地形位相と違って、基線にも視線角にも、衛星高度にも依存しません。変位位相は幾何パラメータから完全に独立しており、これが干渉SARによる変位計測が頑健である理由です。地形位相を測るには軌道を精密に知る必要がありますが、変位を測るのに必要なのは波長だけです。

符号の規約は教科書によって逆に取ることがあります。本質は $|d_{\text{LOS}}| = \lambda|\varphi_{\text{disp}}|/(4\pi)$ という絶対値の関係であり、実データを扱うときは必ず処理ソフトの規約を確認してください。

3次元変位を視線に射影する

地面の変位は3次元ベクトル $\bm{u} = (u_E, u_N, u_U)$（東・北・上）ですが、干渉SARが感じるのはその視線成分だけです。ターゲットから衛星へ向かう単位ベクトルを $\hat{\bm{q}} = -\hat{\bm{p}}$ とすると、視線距離が縮む量は $\bm{u}\cdot\hat{\bm{q}}$ なので、

$$ d_{\text{LOS}} = -\bm{u}\cdot\hat{\bm{q}} $$

$\hat{\bm{q}} = \cos\theta\,\hat{\bm{u}} – \sin\theta\,\hat{\bm{g}}$（$\hat{\bm{g}}$ は衛星直下から見てターゲット方向の水平単位ベクトル）を代入すると、

$$ \begin{equation} d_{\text{LOS}} = -u_U\cos\theta + u_g\sin\theta \end{equation} $$

ここで $u_g = \bm{u}\cdot\hat{\bm{g}}$ はレンジ方向の水平成分です。$\theta = 39^\circ$ で係数を評価すると $\cos\theta = 0.777$、$\sin\theta = 0.629$ です。つまり、

• 鉛直成分には感度0.777 — 1 cm沈下すると $d_{\text{LOS}} = 7.77$ mm、位相にして1.76 rad（0.28フリンジ）
• レンジ方向の水平成分には感度0.629 — 鉛直よりやや鈍い
• アジマス（飛行）方向の水平成分には感度ゼロ

(a) を見ると、赤い矢印（実際に地面が動いた3次元ベクトル $\bm{u}$）のうち、干渉SARが感じるのは緑の矢印、すなわち視線方向に落とした影だけだとわかります。赤と緑の長さの差が「測れない成分」で、とくに紙面の奥に向かうアジマス方向の動きは影の長さを1ミリも変えません。(b) はその係数を視線角の関数として描いたもので、$\theta = 39^\circ$ では鉛直に0.777、レンジ水平に0.629、アジマスにちょうど0です。$\theta$ を20度側に浅くすると鉛直感度が0.94まで上がるので地盤沈下向き、50度側に深くすると水平感度が0.77まで上がるので断層の水平すべり向き、という使い分けが読み取れます。

3つ目が干渉SARの構造的な弱点です。極軌道衛星は南北方向に飛ぶので、南北方向の水平変位には原理的にほぼ見えません。上昇軌道（アセンディング）と下降軌道（ディセンディング）の2つの干渉ペアを組み合わせれば、鉛直成分と東西成分の2つは分離できますが、南北成分は別の手法（ピクセルオフセット追跡、GNSS、マルチアパーチャ干渉）に頼る必要があります。

フリンジの読み方

実際の干渉画像を見て変位を読むとき、次の暗算ができると便利です。Cバンドで沈降域が同心円状の縞になっているとして、縞を中心から外側へ5本数えたとしましょう。$5 \times 2.775 \text{ cm} = 13.9$ cm の視線方向変位です。これが純粋な鉛直沈下なら、$13.9/0.777 = 17.9$ cm の沈下、ということになります。

$\varphi_{\text{topo}}$ と $\varphi_{\text{disp}}$ という2大成分が揃いました。しかし現実の干渉位相には、まだ厄介な同居人がいます。次はそれを整理します。

干渉位相の完全な分解

観測される（ラップされた）干渉位相は、次の成分の和です。

$$ \begin{equation} \psi = \mathcal{W}\left\{\varphi_{\text{flat}} + \varphi_{\text{topo}} + \varphi_{\text{disp}} + \varphi_{\text{atm}} + \varphi_{\text{orb}} + \varphi_{\text{noise}}\right\} \end{equation} $$

まず、これらの成分がどれくらいの大きさで、どんな空間的な形をしているのかを見ておきましょう。

(a) の平地位相だけが桁違いに大きく、30 kmの断面で198 rad（32フリンジ）に達します。これを取り除かないと他の成分は完全に埋もれてしまうので、フラットニングが処理の最初に来るわけです。(b) が残りの成分で、欲しい変位（赤、ピーク12.5 rad = 1.98フリンジ = 5.5 cm）と大気（青、最大4.2 rad = 0.67フリンジ）が同じオーダーで並んでいることに注目してください。大気だけで見かけ1.9 cmの変位が作れてしまうので、1枚の干渉図から数cmの変位を主張するのは危険だ、という警告がここに可視化されています。(c) は全部を足してラップした「実際に観測されるもの」で、赤い変位信号がどこにあるのか目では判別できません。干渉SARの処理とは、この黒い縞から赤い曲線を掘り出す作業だということです。

それぞれの正体と対処法を見ていきます。

$\varphi_{\text{flat}}$（平地位相）: 基準楕円体面を仮定したときの $\Delta R$ の変化による位相。軌道情報から解析的に計算でき、ほぼ完全に除去できます。除去後の干渉図を「フラットニング済み」と呼びます。

$\varphi_{\text{topo}}$（地形位相）: 前節で導いた $k_z h$。外部DEMがあれば計算して除去できますが、DEM誤差 $\delta h$ が $k_z\delta h$ の残差として残ります。DEMがない場合、あるいは意図的にDEMを作りたい場合は、この項こそが求めたい信号になります。

$\varphi_{\text{disp}}$（変位位相）: $\frac{4\pi}{\lambda}d_{\text{LOS}}$。DInSARで求めたい信号です。

$\varphi_{\text{atm}}$（大気位相）: 電波が対流圏・電離圏を通るときの伝搬遅延の、2回の観測間の差です。これが干渉SARの実用上、最も厄介な誤差源です。対流圏の水蒸気分布は時間・空間ともに大きく変動し、視線方向の遅延差で1〜2 cm、山岳部や湿潤地域では数cmに達します。1 cmの遅延差はCバンドで $4\pi \times 0.01/0.0555 = 2.26$ rad、すなわち0.36フリンジに相当します。つまり、大気は数cmの「偽の変位」を作れてしまうわけです。1枚の干渉図から「3 cm沈下した」と主張するのは危険で、大気の可能性を否定できません。

対処法は主に3つあります。(a) 気象モデルや分光観測で遅延を推定して補正する、(b) 大気は時間的に相関がない（各観測日で独立）が変位は時間的に滑らかである、という性質を使って多数の干渉図から統計的に分離する、(c) 空間的に低周波な成分としてフィルタリングする。(b) の思想を体系化したのがPS-InSAR（Ferretti ら, 2001）とSBAS（Berardino ら, 2002）で、数十枚の画像の時系列から年間ミリメートル精度の変位速度を推定する現代の標準手法です。

$\varphi_{\text{orb}}$（軌道誤差位相）: 軌道決定の誤差が $B_\perp$ の誤差となり、画像全体にゆるやかな傾き（ランプ）として現れます。GNSSベースの精密軌道が使える現在は数cm精度で軌道が決まるため、影響は小さくなりましたが、広域では低次多項式でフィッティングして除去します。

$\varphi_{\text{noise}}$（雑音位相）: 熱雑音と、散乱体の配置変化による位相の乱れ。この大きさを支配するのが次節のコヒーレンスです。

こうして並べると、干渉SARの処理は「欲しい成分だけを残して他を剥ぎ取る作業」だとわかります。そして、剥ぎ取りが成功するかどうかの前提条件が、位相にそもそも情報が残っているかどうか — コヒーレンスです。

コヒーレンス — 干渉位相の信頼度

直感

干渉SARの前提は「分解能セル内の散乱体の配置が2回の観測で同じ」でした。もし配置が変わっていたら、散乱位相 $\varphi_{\text{scat}}$ が引き算で消えず、干渉位相はランダムノイズになってしまいます。

草原を想像してください。1回目の観測から12日後、草は伸び、風で倒れ、雨で濡れています。セル内の散乱体の相対配置は完全に入れ替わっており、位相は無関係な値になります。一方、コンクリートの建物の壁や露岩は、12日程度ではまったく変わりません。位相はきれいに保存されます。

この「位相がどれだけ保存されているか」を数値化したのがコヒーレンス（複素相関係数の絶対値）です。

定義

2つのSLC画像を確率過程とみなし、複素相関係数を

$$ \begin{equation} \gamma = \frac{\left| \mathbb{E}\!\left[ s_1 s_2^{*} \right] \right|}{\sqrt{\mathbb{E}\!\left[|s_1|^2\right]\,\mathbb{E}\!\left[|s_2|^2\right]}} \end{equation} $$

と定義します。$0 \le \gamma \le 1$ で、$\gamma = 1$ は完全に位相が保存されている状態（ノイズなし）、$\gamma = 0$ は完全にランダムです。分母の正規化により、反射強度の絶対値には依らない量になっています。

実データでは期待値を取れないので、周囲 $L$ 画素の空間平均（マルチルック）で推定します。

$$ \hat{\gamma} = \frac{\left| \sum_{i=1}^{L} s_{1,i} s_{2,i}^{*} \right|}{\sqrt{\sum_{i=1}^{L} |s_{1,i}|^2 \sum_{i=1}^{L}|s_{2,i}|^2}} $$

この推定量は上に偏るという重要な性質があります。真のコヒーレンスがゼロでも、有限個のランダムな複素数のベクトル和はゼロにならないからです。後でモンテカルロで確かめますが、真値 $\gamma = 0$ に対して $L = 4$ では $\hat\gamma$ の平均が0.46にもなります。「コヒーレンスが0.4ある」と言っても、ルック数が少なければ完全な無相関と区別がつきません。ルック数を明記せずにコヒーレンスの値を語ってはいけない、ということです。

非相関の3+2要因

コヒーレンスは、独立な要因の積に分解できます。

$$ \gamma = \gamma_{\text{temporal}} \cdot \gamma_{\text{geometric}} \cdot \gamma_{\text{volume}} \cdot \gamma_{\text{thermal}} \cdot \gamma_{\text{proc}} $$

(1) 時間非相関 $\gamma_{\text{temporal}}$: 2回の観測の間に散乱体の配置が変わることによる低下です。植生、農地、雪、水面で激しく、都市部や露岩では小さい。経験的に指数モデル

$$ \gamma_{\text{temporal}}(t) = \exp\!\left(-\frac{t}{\tau}\right) $$

がよく使われます。時定数 $\tau$ は森林でCバンド数十日以下、乾燥した裸地で数百日以上と、地表被覆で桁が変わります。$\tau = 30$ 日なら12日間隔でも $\gamma = 0.67$、36日で0.30まで落ちます。波長が長いほど時間非相関に強い（Lバンドは森林でも数年オーダーで干渉が組める）というのが、L バンドSARの大きな利点です。

(2) 幾何（基線）非相関 $\gamma_{\text{geometric}}$: 視線角がわずかに違うと、地表から見た「観測している空間周波数」がずれます。フラットな地表の場合、地上レンジ方向の波数は $k_g = (4\pi f/c)\sin\theta$ です。2つの観測が同じ地上波数を見るには、視線角のずれ $\delta\theta$ を周波数のずれ $\delta f$ で補償する必要があります。$\delta k_g = 0$ の条件から、

$$ \frac{4\pi}{c}\left(\sin\theta\,\delta f + f\cos\theta\,\delta\theta\right) = 0 \;\;\Longrightarrow\;\; \delta f = -\frac{f\,\delta\theta}{\tan\theta} $$

2つの軌道のなす角は $\delta\theta = B_\perp / R$ なので、スペクトルシフトは

$$ |\delta f| = \frac{f\,B_\perp}{R\tan\theta} $$

です。このシフトが送信帯域幅 $W$ を超えると、2つの画像は共通の周波数成分をまったく持たなくなり、コヒーレンスはゼロになります。$|\delta f| = W$ から臨界基線が決まります。

$$ \begin{equation} B_{\perp,\text{crit}} = \frac{W R \tan\theta}{f} = \frac{\lambda R \tan\theta}{2\rho_r} \end{equation} $$

最後の変形では $f = c/\lambda$ と $\rho_r = c/(2W)$ を使いました。Sentinel-1 相当（$W = 56.5$ MHz、$R = 800$ km、$\theta = 39^\circ$、$\lambda = 5.55$ cm）では $B_{\perp,\text{crit}} \approx 6.8$ km となります。実運用では地形の傾斜などで実効的な臨界基線はこれより小さく、5 km程度が目安です。線形モデル $\gamma_{\text{geom}} = 1 – |B_\perp|/B_{\perp,\text{crit}}$ がよく使われます。なお、共通帯域だけを取り出すコモンバンドフィルタをかければ、この非相関はかなり回復できます。

(3) ボリューム非相関 $\gamma_{\text{volume}}$: 森林のように散乱体が高さ方向に分布していると、異なる高さからの反射が異なる位相 $k_z z$ を持ち、和を取ると打ち消しあいます。高さ方向の散乱分布を $F(z)$ とすると、

$$ \gamma_{\text{volume}} = \frac{\left|\int F(z) e^{jk_z z}\,dz\right|}{\int F(z)\,dz} $$

つまり散乱プロファイルのフーリエ変換です。厚さ $h_v$ の一様な層なら、

$$ \gamma_{\text{volume}} = \left|\mathrm{sinc}\!\left(\frac{k_z h_v}{2}\right)\right|, \quad \mathrm{sinc}(x) = \frac{\sin x}{x} $$

$h_v = 20$ mの森林で $B_\perp = 100$ m（$k_z = 0.045$）なら $\gamma_{\text{vol}} = 0.97$ とほぼ無害ですが、$B_\perp = 600$ m（$k_z = 0.270$）では $0.16$ まで落ちます。この式を逆に使って森林の高さを推定するのがPolInSAR（偏波干渉SAR）です。

3つの要因は効き方の「軸」がまったく違います。(1) は時間だけの関数で、同じ12日間隔でも森林（$\tau = 15$ 日）なら $\gamma = 0.45$ まで落ちるのに、都市・露岩（$\tau = 600$ 日）なら0.98を保ちます。地表被覆で桁が変わるという意味がこの縦の広がりです。(2) と (3) はどちらも $B_\perp$ の関数ですが、横軸のスケールが違うことに注意してください — 幾何非相関がゼロになるのは6.78 kmですが、ボリューム非相関は樹高30 mの森林なら $B_\perp = 500$ m程度でもうゼロになります。実運用で使う数百mの範囲（(2) の緑帯）では幾何非相関はほとんど無害で、森林では先にボリューム非相関が効いてくるというのが実務上の結論です。

(4) 熱雑音 $\gamma_{\text{thermal}}$: 受信機雑音による低下で、信号対雑音比 $\mathrm{SNR}$ を使って $\gamma_{\text{thermal}} = (1 + \mathrm{SNR}^{-1})^{-1}$ と書けます。

(5) 処理起因 $\gamma_{\text{proc}}$: レジストレーション誤差（画素合わせのずれ）や補間の誤差による低下。一般に、干渉SARでは0.1画素以下の精度でレジストレーションする必要があります。

コヒーレンスと位相精度を結ぶ式

コヒーレンスがわかると、位相の標準偏差が見積もれます。$L$ ルックの多視処理を行った場合、干渉位相の分散に対するクラメール・ラオ下界は次で与えられます。

$$ \begin{equation} \boxed{\; \sigma_\varphi = \sqrt{\frac{1 – \gamma^2}{2L\gamma^2}} \;} \end{equation} $$

そして変位に換算した標準偏差は $\sigma_d = \frac{\lambda}{4\pi}\sigma_\varphi$ です。数値を並べます（$\lambda = 5.55$ cm）。

「干渉SARはミリメートル精度」という宣伝文句の中身がわかります。$\gamma = 0.8$ で20ルック平均すれば、確かに0.5 mmの位相ノイズです。ただしこれはノイズだけの話であって、大気やDEM誤差といった系統誤差はこの表に入っていません。実際の1枚の干渉図の精度は、大気に支配されてcmオーダーになるのが普通です。

$\sigma_\varphi$ の式には重要な適用限界があります。これはクラメール・ラオ下界であり、位相誤差が小さい（$\sigma_\varphi \ll \pi$）ことを前提とした近似です。位相は $[-\pi,\pi)$ に閉じ込められているので、標準偏差は一様分布の値 $\pi/\sqrt{3} = 1.814$ rad（103.9度）を超えられません。ところが式は $\gamma = 0.2, L = 1$ で $3.46$ radという不可能な値を返します。低コヒーレンス・少ルックの領域では、この式を鵜呑みにしてはいけません。後でモンテカルロシミュレーションでこの破綻を実際に確認します。

コヒーレンスによって位相の信頼度がわかりました。しかし、信頼できる位相であっても、$2\pi$ の巻き数が未知という問題は残ったままです。いよいよ最難関のアンラップに進みます。

位相アンラップ

問題設定と伊藤の条件

観測されるラップ位相 $\psi(x,y) \in [-\pi,\pi)$ から、連続な真の位相 $\phi(x,y)$ を復元したい。しかし $\phi$ と $\phi + 2\pi k$ は区別できないので、このままでは解が定まりません。何らかの先験情報が必要です。

そこで使うのが、真の位相は隣接画素間で急には変わらないという仮定です。数式で書けば、

$$ \begin{equation} \left|\phi(x+1,y) – \phi(x,y)\right| < \pi \quad \text{（すべての隣接ペアで）} \end{equation} $$

これを伊藤の条件（Itoh, 1982）と呼びます。この条件が成り立てば、隣接画素の真の位相差はラップ位相差から一意に復元できます。

$$ \phi(x+1,y) – \phi(x,y) = \mathcal{W}\left\{\psi(x+1,y) – \psi(x,y)\right\} $$

なぜなら、真の位相差が $(-\pi,\pi)$ に入っており、$\mathcal{W}$ は $2\pi$ の倍数を除いた同値な値を $(-\pi,\pi]$ に返すからです。したがって、ある基準点から出発して隣接画素を辿りながら差を積分していけば、位相場全体が復元できます。

$$ \phi(\bm{r}) = \phi(\bm{r}_0) + \int_{C} \mathcal{W}\{\nabla\psi\}\cdot d\bm{l} $$

これが位相アンラップの原理です。1次元なら NumPy の np.unwrap がやっているのと同じことです。

なぜ2次元は難しいのか — 残差（residue）

1次元と2次元の決定的な違いは、経路が一意でないことです。ある点まで辿る道は無数にあり、伊藤の条件が破れている画素が1つでもあると、通る道によって答えが変わってしまいます。

これを定量化するのが残差（residue）です。$2\times 2$ の最小の閉ループを反時計回りに一周し、ラップされた位相差を足し合わせます。

$$ \begin{equation} q = \frac{1}{2\pi}\sum_{i=1}^{4} \mathcal{W}\{\Delta_i \psi\} \end{equation} $$

真の位相 $\phi$ は連続なスカラー場なので、その勾配を閉曲線に沿って一周積分すれば必ずゼロです。したがって伊藤の条件がループ上の4辺すべてで成り立っていれば $q = 0$ になります。逆に $q \ne 0$（実際には $q = \pm 1$）となるループがあれば、そこでは条件が破れています。この点を残差と呼び、$q = +1$ を正残差、$q = -1$ を負残差と言います。

残差の重要な性質を2つ挙げます。

1. 残差は必ず正負ペアで生まれる。閉領域全体で残差を足し合わせると、内部の辺の寄与はすべて隣接ループ間で打ち消し合い、外周に沿った循環だけが残ります。外周の循環がゼロなら、領域内の残差の総和もゼロです。後のシミュレーションでも、正残差101個・負残差101個・総和0という結果を確認します。
2. 残差を囲まない経路を選べば、積分は経路に依らない。正残差と負残差を線分（枝切り線, branch cut）で結び、その線を横切らないように積分すれば、矛盾は生じません。

(a) と (b) は同じ位相場の別の場所を切り出したものです。(a) では赤い矢印でたどる $2\times2$ ループの4辺のラップ位相差を足すと $2\pi$ になり、残差 $q = +1$ が立ちます。一方 (b) の滑らかな場所では和がちょうど0で、どちら回りに積分しても同じ答えになります。数値だけを見比べても、(a) のループの左辺では $-2.4 \to 2.4$ と一気に飛んでおり（真の差は $\pi$ を超えている）、(b) の $-2.4 \to -2.6$ のような緩やかな変化とは対照的です。(c) が枝切り法の発想で、正残差（赤）と負残差（青）を黒い線分で結び、緑の破線のようにその線を横切らない経路だけを使えば、どんな道順で積分しても矛盾が生じません。残差が疎ならこの戦略は非常にうまく働きます。

この観察から古典的な枝切り法（Goldstein, Zebker, Werner, 1988）が生まれました。アルゴリズムは、(a) 全残差を検出し、(b) 正負残差を最短距離で結んで枝切り線を張り、(c) 枝切り線を横切らないように領域を塗りつぶしながら積分する、という流れです。実装が軽く、残差が少ない良質な干渉図では非常にうまく動きます。弱点は、残差が密集すると枝切り線が絡み合って大きな島ができ、その領域がアンラップ不能になることです。

大域最適化としてのアンラップ

現代の主流は、アンラップを最適化問題として定式化する方法です。

最小コスト流れ（MCF）法（Costantini, 1998）は、各辺に付ける $2\pi$ の補正整数 $k_{ij}$ を未知数とし、

$$ \min \sum_{ij} c_{ij}\left|k_{ij}\right| \quad \text{s.t.} \quad (\text{各ループでの残差が解消される)} $$

という整数最小費用流問題に落とします。コスト $c_{ij}$ にコヒーレンスの逆数などを使えば、「信頼できない画素を通って補正を入れる」ことを優先させられます。ネットワークフローの多項式時間アルゴリズムで大域最適解が求まるのが強みです。

最小二乗法（Ghiglia & Romero, 1994）は、真の勾配とラップ勾配の食い違いを最小化します。

$$ \min_\phi \sum \left(\phi_{i+1,j} – \phi_{i,j} – \mathcal{W}\{\Delta_x\psi\}\right)^2 + \left(\phi_{i,j+1} – \phi_{i,j} – \mathcal{W}\{\Delta_y\psi\}\right)^2 $$

停留条件は離散ポアソン方程式 $\nabla^2\phi = \rho$ になり、DCT（離散コサイン変換）で高速に解けます。必ず解が得られるのが利点ですが、誤差を画像全体に薄く広げてしまう（局所的な破綻が広域を汚す）という欠点があります。

品質誘導法（quality-guided） は、実装が簡単で実務でもよく使われるヒューリスティックです。各画素に信頼度（コヒーレンスや位相の局所分散から作る品質マップ）を与え、品質の高い画素から順にアンラップしていきます。優先度付きキューを使えば数十行で書けて、しかも「信頼できる領域を先に確定させ、ノイズ領域は最後に回す」ため誤りが伝播しにくい。本記事の実装ではこれを採用します。

Goldstein フィルタ

アンラップの前に、干渉図のノイズを落として残差を減らしておくと成功率が劇的に上がります。ここで定番なのがGoldstein フィルタ（Goldstein & Werner, 1998）です。

アイデアは「干渉図の局所的なスペクトルにおいて、フリンジは狭帯域のピークとして現れ、ノイズは全周波数に広がっている」という観察です。そこで、小さなパッチ（32×32など）ごとに2次元FFTをかけ、スペクトルの振幅そのものを重みに使って強調します。

$$ \begin{equation} Z_{\text{filt}}(u,v) = \left|\hat{S}(u,v)\right|^{\alpha} Z(u,v) \end{equation} $$

ここで $Z$ はパッチの複素干渉図のFFT、$\hat{S}$ は $|Z|$ を平滑化して最大値1に正規化したものです。$\alpha \in [0,1]$ がフィルタの強さで、$\alpha = 0$ で無処理、$\alpha = 1$ で最強です。強いスペクトル成分（フリンジ）は $|\hat{S}|^\alpha \approx 1$ でそのまま残り、弱い成分（ノイズ）は $|\hat{S}|^\alpha \ll 1$ で押さえ込まれます。

重要なのは、これが適応的であることです。フリンジの方向や間隔がどこでどう変わっても、その場のスペクトルが自動的に重みを決めます。実装では、パッチを重ねてハニング窓で重み付き加算（オーバーラップ加算）し、ブロック境界の不連続を避けます。

注意点として、フィルタは位相を平滑化するので、断層の食い違いのような本物の位相不連続も鈍らせます。$\alpha$ を上げすぎると、フリンジは滑らかになるがピーク変位が過小評価される、というトレードオフがあります。

理論が一通り揃いました。ここからはPythonで、合成データを使って干渉SARの処理チェーン全体を組み立て、実際にミリ単位の変位が復元できるかを確かめます。

Pythonでの実装

幾何パラメータの計算と可視化

まず、これまでに導いた $h_a$、$B_{\perp,\text{crit}}$、視線感度を数値で確かめます。Sentinel-1 相当の条件を使います。

import numpy as np

c0 = 299792458.0          # 光速 [m/s]
lam = 0.0555              # Cバンド波長 [m]
R0 = 800e3                # スラントレンジ [m]
theta = np.deg2rad(39.0)  # 視線角
Wr = 56.5e6               # レンジ帯域幅 [Hz]

def h_amb(Bperp, lam=lam, R=R0, th=theta):
    """高度あいまい度 h_a = λ R sinθ / (2 B_perp)"""
    return lam * R * np.sin(th) / (2.0 * Bperp)

def B_crit(lam=lam, R=R0, th=theta, W=Wr):
    """臨界基線 B_crit = λ R tanθ W / c"""
    return lam * R * np.tan(th) * W / c0

print("1フリンジ = λ/2 = %.2f cm" % (lam / 2 * 100))
print("臨界基線 = %.0f m" % B_crit())
for Bp in [50, 100, 150, 300, 600]:
    ha = h_amb(Bp)
    print("  B_perp=%4.0f m -> h_a=%6.1f m, k_z=%.4f rad/m, 平地縞間隔=%5.1f m"
          % (Bp, ha, 2 * np.pi / ha, ha / np.tan(theta)))

print("鉛直感度 cosθ=%.3f, レンジ水平感度 sinθ=%.3f" % (np.cos(theta), np.sin(theta)))
print("1 cm の沈下 -> 位相 %.2f rad (%.2f フリンジ)"
      % (4 * np.pi * 0.01 * np.cos(theta) / lam,
         4 * np.pi * 0.01 * np.cos(theta) / lam / (2 * np.pi)))

実行すると、1フリンジ2.77 cm、臨界基線6776 m、$B_\perp = 100$ mで $h_a = 139.7$ m、といった値が出ます。表で示した数字と一致しました。注目すべきは最後の行で、1 cmの沈下がわずか0.28フリンジにしかならないことです。干渉画像で縞が1本も閉じていない微小変位を議論するには、後述する時系列解析が必要になる、という現実がここに表れています。

続いて、$B_\perp$ に対する $h_a$ と幾何非相関、そして視線感度の方位依存性を図にします。

import numpy as np
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

Bp = np.linspace(10, 2000, 400)
fig, ax = plt.subplots(1, 3, figsize=(16, 4.5))

ax[0].loglog(Bp, h_amb(Bp), lw=2, color="#1f77b4")
ax[0].axhline(30, ls="--", c="gray"); ax[0].text(12, 33, "山岳でアンラップが困難", fontsize=9)
ax[0].set_xlabel("垂直基線 $B_\\perp$ [m]"); ax[0].set_ylabel("高度あいまい度 $h_a$ [m]")
ax[0].set_title("基線が長いほど地形に敏感（$h_a$ が小さい）"); ax[0].grid(alpha=.3, which="both")

ax[1].plot(Bp, np.clip(1 - Bp / B_crit(), 0, 1), lw=2, color="#d62728")
for hv, cl in [(10, "#2ca02c"), (20, "#9467bd"), (30, "#8c564b")]:
    x = (2 * np.pi / h_amb(Bp)) * hv / 2
    ax[1].plot(Bp, np.abs(np.sin(x) / x), lw=2, color=cl, label="ボリューム（樹高 %d m）" % hv)
ax[1].plot([], [], lw=2, color="#d62728", label="幾何（臨界基線モデル）")
ax[1].set_xlabel("垂直基線 $B_\\perp$ [m]"); ax[1].set_ylabel("コヒーレンス $\\gamma$")
ax[1].set_title("基線が長いほどコヒーレンスは落ちる"); ax[1].legend(fontsize=9); ax[1].grid(alpha=.3)

th_deg = np.linspace(20, 50, 200)
ax[2].plot(th_deg, np.cos(np.deg2rad(th_deg)), lw=2, label="鉛直変位への感度 $\\cos\\theta$")
ax[2].plot(th_deg, np.sin(np.deg2rad(th_deg)), lw=2, label="レンジ水平への感度 $\\sin\\theta$")
ax[2].axvline(39, ls=":", c="k"); ax[2].text(39.5, .3, "$\\theta=39°$", fontsize=9)
ax[2].set_xlabel("視線角 $\\theta$ [deg]"); ax[2].set_ylabel("視線方向への投影係数")
ax[2].set_title("視線角と変位感度（アジマス方向は常に0）"); ax[2].legend(fontsize=9); ax[2].grid(alpha=.3)
plt.tight_layout(); plt.show()

3枚の図から、干渉ペアの設計思想が読み取れます。左は $h_a \propto 1/B_\perp$ の両対数直線で、$B_\perp$ が数百mを超えると $h_a$ が数十mまで落ち、山岳地では隣接画素の位相差が $\pi$ を超えてアンラップが破綻する領域に入ります。中央は、その長基線がコヒーレンスも同時に下げることを示しています。とくに樹高30 mの森林では $B_\perp = 300$ m付近でボリューム非相関が0.44まで落ち、幾何非相関より先に効いてきます。右のパネルは、$\theta$ が浅い（20度側）ほど鉛直変位に敏感で、深い（50度側）ほど水平変位に敏感になることを示します。地盤沈下を狙うなら浅い入射角、断層の水平すべりを狙うなら深い入射角が有利、という指針です。

これで幾何パラメータの感覚がつかめました。次に、コヒーレンスと位相ノイズの関係式が本当に正しいのか、シミュレーションで検証します。

位相標準偏差のモンテカルロ検証

$\sigma_\varphi = \sqrt{(1-\gamma^2)/(2L\gamma^2)}$ が成り立つ範囲を調べます。コヒーレンス $\gamma$ を持つ複素ガウス散乱を人工的に作るには、独立な複素ガウス $a, b$ から

$$ s_1 = a, \qquad s_2 = \gamma a + \sqrt{1-\gamma^2}\, b $$

とすればよく、$\mathbb{E}[s_1 s_2^*] = \gamma$、$\mathbb{E}[|s_i|^2] = 1$ になります。これを $L$ サンプル平均して位相のばらつきを測ります。

import numpy as np
rng = np.random.default_rng(1)

def sim_phase_std(gamma, L, M=200000, rng=rng):
    """コヒーレンスγの複素ガウス対をL多視し、干渉位相の標準偏差を実測"""
    a = (rng.standard_normal((M, L)) + 1j * rng.standard_normal((M, L))) / np.sqrt(2)
    b = (rng.standard_normal((M, L)) + 1j * rng.standard_normal((M, L))) / np.sqrt(2)
    s1 = a
    s2 = gamma * a + np.sqrt(1 - gamma**2) * b     # 相関を持つ2枚目
    ph = np.angle((s1 * np.conj(s2)).sum(axis=1))  # L多視後の干渉位相
    return np.sqrt(np.mean(ph**2))                 # 真値0まわりのばらつき

print(" γ    L   実測σ[rad]  CRB[rad]   比")
for g in [0.9, 0.8, 0.6, 0.4, 0.2]:
    for L in [1, 5, 20]:
        emp = sim_phase_std(g, L)
        crb = np.sqrt((1 - g**2) / (2 * L * g**2))
        print("%4.1f %4d %9.3f %9.3f %6.2f" % (g, L, emp, crb, emp / crb))
print("一様分布の限界 π/√3 = %.3f rad" % (np.pi / np.sqrt(3)))

出力から3つのことが読み取れます。第一に、$\gamma = 0.9$, $L = 20$ では実測0.079 rad に対しCRBが0.077 radで、比は1.03 — 高コヒーレンス・多ルックの領域ではCRBが実用的な精度予測になります。第二に、$\gamma \ge 0.4$ の範囲で $L = 5$ という少ルックにすると比は1.14〜1.33となり、CRBは真の誤差を1〜3割過小評価します。第三に、$\gamma = 0.2$, $L = 1$ では実測1.64 radに対しCRBが3.46 radで、比が0.47 — CRBが物理的にありえない値を返しています。実測値が一様分布の限界 $\pi/\sqrt3 = 1.814$ rad を超えていないことも確認できます。低コヒーレンス域でCRBを使ってはいけないという結論です。

(a) の実線がCRB、同色の白丸が実測です。$\gamma$ が高い右側では両者がぴたりと重なりますが、左に行くほど実線だけが $1/\gamma$ で発散し、白丸は $\pi/\sqrt{3} = 1.814$ radの水平線（一点鎖線）を決して超えません。位相が $[-\pi,\pi)$ に閉じ込められている以上、これは当たり前の頭打ちですが、CRBの式はその事情を知らないのです。(b) は両者の比で、$L = 20$（青）はどの $\gamma$ でも1.0〜1.25に収まる一方、$L = 1$（赤）は $\gamma = 0.9$ で2.0、$\gamma = 0.2$ で0.47と、過小評価と過大評価の両方向に大きく外れます。理論式が使えるのは高コヒーレンス・多ルックの領域だけという適用範囲が、この1枚に凝縮されています。

コヒーレンス推定量のバイアスも確認しておきましょう。

import numpy as np
rng = np.random.default_rng(7)

print("真のコヒーレンスが 0 のときの推定値の平均:")
for L in [4, 16, 64]:
    a = (rng.standard_normal((20000, L)) + 1j * rng.standard_normal((20000, L))) / np.sqrt(2)
    b = (rng.standard_normal((20000, L)) + 1j * rng.standard_normal((20000, L))) / np.sqrt(2)
    num = np.abs((a * np.conj(b)).sum(axis=1))
    den = np.sqrt((np.abs(a)**2).sum(axis=1) * (np.abs(b)**2).sum(axis=1))
    print("  L = %2d -> γ_hat の平均 = %.3f" % (L, np.mean(num / den)))

真値がゼロにもかかわらず、$L = 4$ では平均0.460、$L = 16$ で0.225、$L = 64$ で0.110 という推定値が出ます。おおむね $1/\sqrt{L}$ で減っていく上方バイアスです。コヒーレンスマップを見て「0.4以上あるから使える」と判断するときは、必ずルック数を確認しなければならない、という実務上の教訓がここにあります。

(a) のヒストグラムは、真のコヒーレンスがゼロ（左端の黒い縦線）にもかかわらず、推定値がどこに散らばるかを示しています。$L = 4$（赤）では分布が0.2〜0.8に広く伸び、0.6を超える値すら珍しくありません。$L = 64$（青）まで多視すると0.1付近に集中します。(b) はその平均を両対数でプロットしたもので、実測点が $\sqrt{\pi}/(2\sqrt{L})$ の破線にほぼ乗っており、バイアスが $1/\sqrt{L}$ で減ることが定量的に確認できます。ここから、たとえば「$L = 4$ でコヒーレンス0.46」という報告は完全な無相関と区別がつかない、と即座に判断できるわけです。

理論式の適用範囲がわかったので、いよいよ合成干渉図を作ります。

合成干渉図の生成

地形（滑らかな起伏）とガウス型の沈降域を持つ地表を作り、干渉位相を組み立てます。$B_\perp = 100$ m、画素間隔40 m、256×256画素とします。

import numpy as np

rng = np.random.default_rng(0)
N, dpix = 256, 40.0
lam, R0, theta, Bperp = 0.0555, 800e3, np.deg2rad(39.0), 100.0
ha = lam * R0 * np.sin(theta) / (2 * Bperp)          # 高度あいまい度
x = (np.arange(N) - N / 2) * dpix
X, Y = np.meshgrid(x, x)

def smooth_field(N, corr_px, rng):
    """相関長 corr_px 画素のガウス平滑ランダム場（標準化）"""
    z = rng.standard_normal((N, N))
    fx = np.fft.fftfreq(N)[:, None]
    fy = np.fft.fftfreq(N)[None, :]
    H = np.exp(-2 * (np.pi * corr_px)**2 * (fx**2 + fy**2))
    out = np.real(np.fft.ifft2(np.fft.fft2(z) * H))
    return out / out.std()

# 地形（起伏 + 1つの丘）と、視線方向に6cm遠ざかる沈降域
h = 150.0 * smooth_field(N, 18.0, rng) \
    + 200.0 * np.exp(-((X - 2000)**2 + (Y + 1500)**2) / (2 * 2500.0**2))
h -= h.mean()
d_los = 0.060 * np.exp(-((X + 1500)**2 + (Y - 1000)**2) / (2 * 1200.0**2))

phi_topo = 2 * np.pi * h / ha             # 地形位相 = k_z * h
phi_disp = 4 * np.pi * d_los / lam        # 変位位相
phi_true = phi_topo + phi_disp
wrap = lambda p: np.angle(np.exp(1j * p)) # ラッピング演算子

print("h_a = %.1f m, 標高レンジ %.0f m -> 地形位相 %.1f フリンジ"
      % (ha, h.ptp(), phi_topo.ptp() / (2 * np.pi)))
print("沈降ピーク %.1f cm -> 変位位相 %.2f フリンジ" % (d_los.max()*100, phi_disp.max()/(2*np.pi)))
gx, gy = np.diff(phi_true, axis=1), np.diff(phi_true, axis=0)
print("隣接画素の最大位相差 %.3f rad (π=%.3f) -> 伊藤の条件を満たす"
      % (max(np.abs(gx).max(), np.abs(gy).max()), np.pi))

標高レンジ884 mが地形位相6.3フリンジに、沈降ピーク6.0 cmが変位位相2.16フリンジに対応します。$h_a = 139.7$ mなので、884/139.7 = 6.3 と手計算とも一致しています。最後の行が重要で、隣接画素間の真の位相差の最大値が0.650 radと $\pi$ より十分小さく、ノイズがなければアンラップは原理的に必ず成功する設定になっていることを確認しています。逆にこの値が $\pi$ を超えていたら、どんな高級なアンラッパーを使っても正解には到達できません。

この位相場を可視化します。

import numpy as np
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

ext = [x[0]/1000, x[-1]/1000, x[0]/1000, x[-1]/1000]
fig, ax = plt.subplots(2, 2, figsize=(11, 9.5))
panels = [
    (h, "terrain", "標高 $h$ [m]", "(a) 地形（真値）"),
    (d_los * 100, "inferno", "視線方向変位 [cm]", "(b) 沈降域（真値, 遠ざかる向きが正）"),
    (phi_true, "viridis", "位相 [rad]", "(c) アンラップ済み真位相 $\\varphi_{topo}+\\varphi_{disp}$"),
    (wrap(phi_true), "hsv", "位相 [rad]", "(d) 干渉図（ラップ済み・雑音なし）"),
]
for a, (dat, cm, cb, ti) in zip(ax.ravel(), panels):
    im = a.imshow(dat, cmap=cm, extent=ext, origin="lower")
    plt.colorbar(im, ax=a, fraction=.046, label=cb)
    a.set_title(ti, fontsize=11); a.set_xlabel("東西 [km]"); a.set_ylabel("南北 [km]")
plt.tight_layout(); plt.show()

(a)の地形と(c)の真位相を見比べると、位相が標高の等高線をそのままなぞっていることがわかります。(b)の沈降域は左上にありますが、(c)ではその位置に地形とは無関係な同心円状の盛り上がりが加わっています。そして(d)が実際に観測される干渉図です。地形起源の縞と変位起源の同心円が重なり、このままでは変位だけを目で読み取れないことが一目でわかります。DInSARでDEMを使って地形位相を引く必要性が、この図から実感できるはずです。

次に、現実的なコヒーレンス低下と位相ノイズを加えます。

コヒーレンスマップとノイズ注入・残差の検出

右上に「植生に覆われて時間非相関が強い領域」を作り、コヒーレンスに応じた位相ノイズを注入します。ノイズの標準偏差にはCRBの式を使います（前節で $\gamma \gtrsim 0.15$、$L = 20$ なら妥当と確認済みです）。

import numpy as np

# コヒーレンスマップ（右上に低コヒーレンス域を作る）
coh = 0.85 - 0.25 * (0.5 * smooth_field(N, 20.0, rng) + 0.5)
coh = np.clip(coh - 0.55 * np.exp(-((X - 2500)**2 + (Y - 2500)**2) / (2 * 1500.0**2)), 0.15, 0.95)

L_looks = 20
sigma_phi = np.sqrt((1 - coh**2) / (2 * L_looks * coh**2))   # 画素ごとの位相σ
psi = wrap(phi_true + sigma_phi * rng.standard_normal((N, N)))

def residues(p):
    """2x2 ループを一周したラップ位相差の和 / 2π（残差）"""
    d1 = wrap(p[:-1, 1:] - p[:-1, :-1])
    d2 = wrap(p[1:,  1:] - p[:-1, 1:])
    d3 = wrap(p[1:, :-1] - p[1:,  1:])
    d4 = wrap(p[:-1, :-1] - p[1:, :-1])
    return np.round((d1 + d2 + d3 + d4) / (2 * np.pi)).astype(int)

r_raw = residues(psi)
print("コヒーレンス 平均 %.2f (最小 %.2f / 最大 %.2f)" % (coh.mean(), coh.min(), coh.max()))
print("残差: 正 %d 個 / 負 %d 個 / 総和 %d" % ((r_raw > 0).sum(), (r_raw < 0).sum(), r_raw.sum()))
for lo, hi in [(0.8, 0.95), (0.5, 0.8), (0.3, 0.5), (0.15, 0.3)]:
    m = (coh[:-1, :-1] >= lo) & (coh[:-1, :-1] < hi)
    if m.sum() > 100:
        print("  γ %.2f-%.2f: 残差密度 %.3f %%" % (lo, hi, 100 * np.mean(r_raw[m] != 0)))

結果は劇的です。正残差101個、負残差101個、総和はぴったり0 — 理論で述べた「残差は正負ペアで生まれる」が数値実験で確認できました。さらに残差密度をコヒーレンス別に見ると、$\gamma > 0.3$ の領域では残差密度が0.000 %、$\gamma = 0.15$〜$0.3$ の領域だけで5.04 %と、残差が完全に低コヒーレンス域に集中しています。位相ノイズはどこにでも乗っているのに、アンラップを壊すレベルの誤りは低コヒーレンス域でしか起きないわけです。これがコヒーレンスを品質マップとして使う正当性の根拠になります。

この状況を図で確認します。

import numpy as np
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

fig, ax = plt.subplots(1, 3, figsize=(16, 4.8))
im = ax[0].imshow(coh, cmap="gray", vmin=0, vmax=1, extent=ext, origin="lower")
plt.colorbar(im, ax=ax[0], fraction=.046, label="コヒーレンス $\\gamma$")
ax[0].set_title("(a) コヒーレンスマップ（右上＝植生で非相関）")

im = ax[1].imshow(psi, cmap="hsv", extent=ext, origin="lower")
plt.colorbar(im, ax=ax[1], fraction=.046, label="位相 [rad]")
ax[1].set_title("(b) 雑音を含む干渉図（%dルック相当）" % L_looks)

ry, rx = np.nonzero(r_raw)
ax[2].imshow(coh, cmap="gray", vmin=0, vmax=1, extent=ext, origin="lower", alpha=.7)
pos = r_raw[ry, rx] > 0
ax[2].scatter(x[rx[pos]]/1000, x[ry[pos]]/1000, s=6, c="red",  label="正残差 %d" % pos.sum())
ax[2].scatter(x[rx[~pos]]/1000, x[ry[~pos]]/1000, s=6, c="blue", label="負残差 %d" % (~pos).sum())
ax[2].legend(fontsize=9, loc="lower left")
ax[2].set_title("(c) 残差の分布（低コヒーレンス域に集中）")
for a in ax: a.set_xlabel("東西 [km]"); a.set_ylabel("南北 [km]")
plt.tight_layout(); plt.show()

(b)を(a)と見比べると、右上の低コヒーレンス域では縞が完全に潰れて「砂嵐」になっているのに対し、それ以外の領域では縞がきれいに読めます。(c)の残差プロットは決定的で、赤（正残差）と青（負残差）が右上の暗い（低コヒーレンス）領域にほぼ限定され、しかも赤と青がペアで近接して分布しています。枝切り法が「近い正負残差を結ぶ」という単純な戦略で機能する理由が、この図から直感的に理解できます。

残差が可視化できたので、これを減らすためにGoldsteinフィルタをかけます。

Goldstein フィルタの実装

32×32のパッチを16画素ずつずらしながら、ハニング窓でオーバーラップ加算します。

import numpy as np

def goldstein_filter(z, alpha=0.7, W=32, step=16):
    """複素干渉図 z にGoldsteinフィルタを適用（オーバーラップ加算）"""
    Nn = z.shape[0]
    out = np.zeros_like(z, dtype=complex)
    wsum = np.zeros((Nn, Nn))
    win = np.hanning(W + 2)[1:-1]
    w2 = np.outer(win, win)              # 2次元ハニング窓
    k = np.ones((3, 3)) / 9.0            # スペクトル平滑化カーネル
    for i in range(0, Nn - W + 1, step):
        for j in range(0, Nn - W + 1, step):
            Z = np.fft.fft2(z[i:i+W, j:j+W])
            # |Z| を平滑化して最大値1に正規化 → これを重みにする
            S = np.abs(np.real(np.fft.ifft2(
                np.fft.fft2(np.abs(Z)) * np.fft.fft2(k, s=Z.shape))))
            S /= (S.max() + 1e-12)
            out[i:i+W, j:j+W] += np.fft.ifft2(Z * S**alpha) * w2
            wsum[i:i+W, j:j+W] += w2
    wsum[wsum == 0] = 1.0
    return out / wsum

psi_f = np.angle(goldstein_filter(np.exp(1j * psi), alpha=0.7))
r_flt = residues(psi_f)
print("残差数: フィルタ前 %d -> 後 %d (%.1f %% 減)"
      % (np.count_nonzero(r_raw), np.count_nonzero(r_flt),
         100 * (1 - np.count_nonzero(r_flt) / np.count_nonzero(r_raw))))
print("ラップ位相RMSE: %.3f rad -> %.3f rad"
      % (np.sqrt(np.mean(wrap(psi - phi_true)**2)),
         np.sqrt(np.mean(wrap(psi_f - phi_true)**2))))

$\alpha = 0.7$ のフィルタで、残差は202個から7個へ96.5 %減、ラップ位相のRMSEは0.276 radから0.154 radへほぼ半減しました。残差がここまで減れば、あとはどんなアンラッパーでもほぼ正解に辿り着けます。位相RMSEの減少（1.8倍）に比べて残差の減少（29倍）が圧倒的に大きいことに注目してください。残差は「位相差が $\pi$ を超えた」という閾値的な事象なので、ノイズをわずかに下げるだけで発生確率が急激に落ちるのです。フィルタの効果が非線形に効く、というのが干渉SAR処理の面白いところです。

上段がフィルタの中身です。(a) の32×32パッチをフーリエ変換すると、(b) のようにフリンジ由来の1点だけが突出し、雑音は全面に薄く広がります。(c) はそれを平滑化・正規化して0.7乗した重みで、ピーク付近だけが1に近く、それ以外は0.3以下に落ちています。(d) の断面を見ると、ピーク（中央）はほぼ素通しなのに、その両側の雑音成分は1桁近く押し下げられていることがわかります。下段が結果で、(e) の低コヒーレンス域の「砂嵐」が (f) では滑らかな縞に変わり、残差が202個から7個に減りました。ノイズの振幅は2倍弱しか減っていないのに残差が29倍減るのは、残差が「位相差が $\pi$ を超えたか」という閾値事象だからです。

ノイズを落とせたので、アンラップに進みます。

品質誘導アンラップの実装

コヒーレンスを品質マップとして、優先度付きキューで高品質画素から順にアンラップします。

import numpy as np, heapq

def unwrap_quality_guided(psi, quality):
    """品質の高い画素から順に伊藤の条件で積分していくアンラップ"""
    Nn = psi.shape[0]
    unw = np.zeros_like(psi)
    done = np.zeros((Nn, Nn), bool)
    s = np.unravel_index(np.argmax(quality), quality.shape)   # 最良画素から出発
    unw[s], done[s] = psi[s], True
    heap = []
    def push(i, j, fi, fj):
        if 0 <= i < Nn and 0 <= j < Nn and not done[i, j]:
            heapq.heappush(heap, (-quality[i, j], i, j, fi, fj))  # 品質が高い順
    for di, dj in ((1,0), (-1,0), (0,1), (0,-1)):
        push(s[0]+di, s[1]+dj, s[0], s[1])
    while heap:
        _, i, j, fi, fj = heapq.heappop(heap)
        if done[i, j]:
            continue
        unw[i, j] = unw[fi, fj] + wrap(psi[i, j] - psi[fi, fj])  # 伊藤の条件
        done[i, j] = True
        for di, dj in ((1,0), (-1,0), (0,1), (0,-1)):
            push(i+di, j+dj, i, j)
    return unw

真値と比較する評価関数を用意し、フィルタの有無・走査順の違いを比較します。アンラップ結果には全体を貫く $2\pi$ の定数不定性があるので、中央値でオフセットを合わせてから比べます。

import numpy as np

def score(unw, truth):
    off = 2 * np.pi * np.round(np.median(unw - truth) / (2 * np.pi))  # 定数不定性を除去
    e = unw - off - truth
    return np.sqrt(np.mean(e**2)), 100 * np.mean(np.abs(e) > np.pi), unw - off

for name, p in [("フィルタなし", psi), ("Goldstein後", psi_f)]:
    rm, bad, u = score(unwrap_quality_guided(p, coh), phi_true)
    print("%s + 品質誘導 : RMSE %.3f rad, |誤差|>π の画素 %.2f %%" % (name, rm, bad))
    if name == "Goldstein後":
        unw_final = u

u_raster = np.unwrap(np.unwrap(psi_f, axis=1), axis=0)   # 素朴なラスタ走査
rm, bad, _ = score(u_raster, phi_true)
print("Goldstein後 + ラスタ走査 : RMSE %.3f rad, |誤差|>π の画素 %.2f %%" % (rm, bad))

3つの結果が出ます。フィルタなし＋品質誘導ではRMSE 1.072 rad、$\pi$ を超える誤り画素が1.30 %。Goldstein後＋品質誘導ではRMSE 0.292 rad、誤り画素0.15 %と9倍近く改善します。そしてGoldstein後でも素朴なラスタ走査（行方向・列方向に順に np.unwrap）ではRMSE 0.585 rad、誤り画素0.41 %と、品質誘導の2倍悪くなります。ラスタ走査は低コヒーレンス域を必ず横断してしまうので、そこで生じた $2\pi$ の誤りが下流の行すべてに伝播するからです。フィルタと走査順のどちらも効くことが定量的に確認できました。

上段のアンラップ結果は3枚とも一見似ていますが、下段の誤差図が決定的な違いを見せます。(a) のフィルタなしでは右上の低コヒーレンス域が青と赤の縦縞でびっしり埋まっており、$2\pi$ 単位の誤りが密集していることがわかります。(b) のGoldstein後は誤差がほぼ真っ白で、残っているのは細い線が1本だけです。(c) のラスタ走査は、その細い線が上方向に長く伸びているのが特徴で、これは低コヒーレンス域で1回だけ生じた $2\pi$ の誤りが、走査方向に沿って下流の画素すべてに引き継がれている様子です。品質誘導が誤りを局所に閉じ込められるのは、悪い画素を後回しにして「まわり道」できるからだ、と誤差図が語っています。

最後に、地形位相を除去して変位マップを復元します。

変位マップの復元と誤差評価

DEMが完全に既知であるとして地形位相を差し引き、$d = \frac{\lambda}{4\pi}\varphi$ で変位に換算します。

import numpy as np

d_est = lam / (4 * np.pi) * (unw_final - phi_topo)   # 地形位相を除去して変位へ
err_mm = (d_est - d_los) * 1000
mask = coh > 0.5

print("変位マップ RMSE = %.2f mm (全画素) / %.2f mm (γ>0.5 の %.0f %%)"
      % (np.sqrt(np.mean(err_mm**2)), np.sqrt(np.mean(err_mm[mask]**2)), 100 * mask.mean()))
print("最大誤差 %.1f mm（低コヒーレンス域）" % np.abs(err_mm).max())
print("推定沈降ピーク %.2f cm (真値 %.2f cm)" % (d_est.max() * 100, d_los.max() * 100))
for dh in [5, 10, 20]:
    print("DEM誤差 %2d m -> 見かけの変位 %.2f mm" % (dh, lam / (4*np.pi) * (2*np.pi*dh/ha) * 1000))

全画素でRMSE 1.29 mm、コヒーレンスが0.5を超える画素（全体の84 %）に限ればRMSE 0.91 mmです。表で予測した「$\gamma = 0.6$, $L = 20$ で $\sigma_d = 0.93$ mm」とみごとに一致しました。推定沈降ピークは6.01 cmで、真値6.00 cmとの差はわずか0.1 mm — 6 cmの沈降を0.2 %の誤差で復元できています。一方、最大誤差は34.2 mmで、これはすべて低コヒーレンス域に集中しています。干渉SARの誤差が「全体的に少し悪い」のではなく「特定の場所だけ壊滅的に悪い」という空間的に極めて不均一な性質を持つことが、この2つの数字の対比によく表れています。

DEM誤差の影響も並べてあります。$h_a = 139.7$ mのペアでは、DEM誤差10 mが2.0 mm、20 mが4.0 mmの偽の変位を生みます。今回のノイズ由来のRMSEが0.9 mmだったことを思えば、この設定ではDEM誤差のほうがノイズより大きな誤差源になり得ることがわかります。

結果を可視化して締めましょう。

import numpy as np
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

fig, ax = plt.subplots(2, 2, figsize=(12, 10))
im = ax[0,0].imshow(d_los*100, cmap="inferno", extent=ext, origin="lower", vmin=0, vmax=6)
plt.colorbar(im, ax=ax[0,0], fraction=.046, label="変位 [cm]"); ax[0,0].set_title("(a) 真の視線方向変位")

im = ax[0,1].imshow(d_est*100, cmap="inferno", extent=ext, origin="lower", vmin=0, vmax=6)
plt.colorbar(im, ax=ax[0,1], fraction=.046, label="変位 [cm]"); ax[0,1].set_title("(b) 復元された変位")

im = ax[1,0].imshow(err_mm, cmap="coolwarm", extent=ext, origin="lower", vmin=-5, vmax=5)
plt.colorbar(im, ax=ax[1,0], fraction=.046, label="誤差 [mm]")
ax[1,0].set_title("(c) 誤差（低コヒーレンス域に集中）")

i0 = N // 2 + 25   # 沈降中心を通る東西プロファイル
ax[1,1].plot(x/1000, d_los[i0]*100, lw=2.5, label="真値")
ax[1,1].plot(x/1000, d_est[i0]*100, lw=1.2, ls="--", color="crimson", label="復元値")
ax[1,1].set_xlabel("東西 [km]"); ax[1,1].set_ylabel("視線方向変位 [cm]")
ax[1,1].set_title("(d) 沈降中心を通る断面"); ax[1,1].legend(); ax[1,1].grid(alpha=.3)
for a in ax.ravel()[:3]: a.set_xlabel("東西 [km]"); a.set_ylabel("南北 [km]")
plt.tight_layout(); plt.show()

(a)と(b)は目視ではほぼ区別がつきません。6 cmのガウス型沈降が、形も振幅も忠実に復元されています。(c)の誤差図は右上の低コヒーレンス域だけが赤青にざわついており、それ以外は白（誤差ほぼゼロ）です。(d)の断面プロファイルでは、復元値（赤破線）が真値（青実線）にぴったり重なっています。ピーク付近でわずかな平滑化も見られませんから、$\alpha = 0.7$ のGoldsteinフィルタはこの空間スケールの信号を損なっていないと言えます。ただし、断層の食い違いのように空間的に急峻な信号を扱う場合は、フィルタ強度を下げるか、フィルタ後の値でピークを議論しない配慮が必要です。

まとめ

本記事では、干渉SAR（InSAR）の原理を、電波の位相から地表変位を測るという発想の出発点から、実際の処理チェーンの実装まで通して解説しました。

• 干渉位相の基本式は $\Delta\varphi = \frac{4\pi}{\lambda}\Delta R$ — 2枚のSLC画像の複素共役積を取ることで、分解能セル内のランダムな散乱位相が打ち消され、視線距離の変化だけが残ります。1フリンジが半波長（Cバンドで2.77 cm）に対応します
• 基線ベクトルの視線直交成分 $B_\perp$ が地形感度を支配する — $\Delta R = -B\sin(\theta-\alpha)$ を視線角で微分すると $-B_\perp$ が現れ、地形位相 $\varphi_{\text{topo}} = \frac{4\pi B_\perp h}{\lambda R\sin\theta}$ が導かれます
• 高度あいまい度 $h_a = \frac{\lambda R\sin\theta}{2B_\perp}$ が干渉ペアの性格を一言で表す — 変位計測なら短基線（大きい $h_a$、DEM誤差に鈍感）、地形計測なら長基線（小さい $h_a$）です
• 変位位相 $\varphi_{\text{disp}} = \frac{4\pi}{\lambda}d_{\text{LOS}}$ は幾何パラメータから独立 — ただし測れるのは視線成分だけで、鉛直に $\cos\theta = 0.78$、レンジ水平に $\sin\theta = 0.63$、アジマス（南北）方向にはゼロの感度しか持ちません
• 干渉位相は平地・地形・変位・大気・軌道・雑音の和 — なかでも大気は1 cmの遅延差が0.36フリンジ（Cバンド）に化けるため、1枚の干渉図だけで数cmの変位を主張するのは危険です。PS-InSARやSBASのような時系列解析が現代の標準です
• コヒーレンス $\gamma$ が位相の信頼度を決める — 時間・幾何・ボリューム・熱雑音・処理の5要因の積で表され、臨界基線 $B_{\perp,\text{crit}} = \lambda R\tan\theta/(2\rho_r)$ やボリューム非相関のsinc則が定量的に導けます
• $\sigma_\varphi = \sqrt{(1-\gamma^2)/(2L\gamma^2)}$ には適用限界がある — モンテカルロ検証により、高コヒーレンス・多ルック域では誤差1〜3 %の良い近似だが、$\gamma = 0.2$, $L = 1$ では実測の2倍以上の値を返して破綻することを確認しました。位相の標準偏差は $\pi/\sqrt{3} = 1.814$ radを超えられません
• 位相アンラップは残差との戦い — 伊藤の条件が破れた場所に生じる残差は必ず正負ペアで現れ（実験でも正101・負101・総和0）、しかも低コヒーレンス域に集中します。Goldsteinフィルタで残差を96.5 %減らし、品質誘導アンラップで信頼度の高い画素から積分することで、素朴なラスタ走査の2倍の精度が得られました
• 合成データでの最終性能 — 6 cmのガウス型沈降を、全画素RMSE 1.29 mm、高コヒーレンス域で0.91 mm、ピーク値の誤差0.1 mmで復元できました。これは理論表の $\sigma_d$ とよく一致します

干渉SARは「位相という捨てられがちな情報を活かすと、分解能の数千分の一の変化が見える」という、計測工学の中でも際立って美しい実例です。同時に、その精度は $2\pi$ のあいまいさ、コヒーレンスの喪失、大気の揺らぎという3つの制約と常に隣り合わせでもあります。この記事で扱ったのは1つの干渉ペアの処理でしたが、実務では数十〜数百枚の画像を積み重ね、時間方向の冗長性でこれらの制約を突破していきます。その入口に立つための道具立ては、ここまでで一通り揃ったはずです。

次のステップとして、以下の記事も参考にしてください。

• 合成開口レーダー（SAR）の原理
• SARのレンジドップラーアルゴリズム
• ISARイメージング — 回転する目標の像を作る
• レーダー方程式 — 受信電力と探知距離の関係