位相アンラップとは？2π折り返しを復元するアルゴリズムの理論と実装

77 GHzのミリ波レーダーで、目標までの距離をミリメートル単位で測りたいとします。反射波の位相 $\phi$ は距離 $R$ に対して $\phi = 4\pi R/\lambda$ で変化するので、位相を読めば $\lambda/2 = 1.95$ mm の分解能どころか、その1/100の精度で距離が求まります。ところが実際に受信信号から位相を取り出すと、値は必ず $-\pi$ から $\pi$ の間に収まってしまいます。目標が 2 mm 動いただけで位相計は一周し、また元の値に戻ってしまう。位相そのものは物理量を高精度に運んでいるのに、観測できるのはその「時計の針の位置」だけなのです。

この、針の位置から「何周したか」を復元する処理が 位相アンラップ（phase unwrapping） です。地味な名前ですが、これができないと成立しない技術は驚くほど多くあります。

上のグラフの青い線が本当に知りたい位相で、距離や標高に比例してどこまでも増えていきます。しかし赤い線、つまり実際に測定器から得られる値は、$-\pi$ に達するたびに $+\pi$ へジャンプするノコギリ波にしかなりません。緑の矢印が示す「失われている $2\pi$ の整数倍」を、下段のノコギリ波だけを手がかりに埋め戻すのがこの記事のテーマです。

• 干渉SAR（InSAR）による標高・地殻変動計測 — 2枚のSAR画像の位相差から地形を測ります。Sentinel-1（C帯、$\lambda = 5.55$ cm）で入射角35度・基線長150 mなら、位相が $2\pi$ 進むごとに標高が約 85 m 変わります。数千メートルの山を測るには縞を何十本も数え直す必要があります。差分干渉（DInSAR）で地盤沈下を測る場合はさらにシビアで、$2\pi$ が 2.77 cm に相当します。
• FMCWレーダーの高精度測距・振動計測 — ビート周波数から得られる粗い距離に、位相から得られる細かい距離を足し込みます。位相の折り返しを解かない限り、$\lambda/2$ 以上の変位は追跡できません。
• 光干渉計測・ホログラフィ・OCT — 干渉縞の位相が表面形状や屈折率分布に対応します。
• MRIの磁場不均一補正・水脂肪分離、ToFカメラの距離画像、GNSS搬送波位相測位のアンビギュイティ決定 — いずれも本質的に同じ「未知の整数を決める」問題です。

本記事では、この問題を数学的にきちんと定式化し、1次元での完全解から、2次元で必ず現れる「レジデュー」という障害、そして実務で最も広く使われる最小二乗アンラップまでを、導出を省略せずに追いかけます。

本記事の内容

• 位相が $[-\pi, \pi)$ に折り返される仕組みと、ラップ作用素 $\mathcal{W}$ の性質
• 観測位相 $\psi = \phi + 2\pi k$ における未知整数 $k$ の決定問題としての定式化
• 1次元アンラップの原理と、それが成立する条件（Itohの定理）— 標本化定理との等価性
• 雑音による位相跳躍の発生確率と、ガウス近似が破綻する理由
• 2次元で生じる経路依存性と、その正体であるレジデュー（離散版のStokesの定理）
• Goldsteinの枝切り法の考え方
• 最小二乗アンラップの導出（離散ポアソン方程式）と、DCTによる高速解法
• 重み付き最小二乗とPCG（前処理付き共役勾配法）
• Pythonでの実装：チャープの群遅延推定、レジデュー可視化、SNR掃引によるRMS誤差比較

前提知識

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

• 解析信号と瞬時周波数
• 干渉SAR（InSAR）の原理
• FMCWレーダーの原理
• サンプリング定理

位相はなぜ折り返されるのか

アナログ時計を思い浮かべてください。針が「3時の方向」を指しているとき、それが今日の3時なのか、15時なのか、あるいは3日と3時間経った後なのかは、針を見ただけでは分かりません。針の位置が持っている情報は、経過時間を12時間で割った余りだけです。位相もまったく同じで、複素数の偏角として測る限り、$2\pi$ の整数倍の情報は原理的に失われます。

もう少し正確に言うと、こうです。私たちが実際に手にするのは複素振幅 $z$ であって、位相 $\phi$ そのものではありません。$z$ から位相を取り出すには

$$ \psi = \arctan_2\!\left(\operatorname{Im} z,\ \operatorname{Re} z\right) $$

を計算します。しかし $e^{j\phi} = e^{j(\phi + 2\pi k)}$ が任意の整数 $k$ について成り立つので、複素数 $z$ を見ても $\phi$ と $\phi + 2\pi$ は完全に同一です。$\arctan_2$ はこの多価性を「主値」に押し込めて、必ず $[-\pi, \pi)$ の値を返します。つまり 折り返しは計算機の実装上の都合ではなく、複素指数関数の周期性そのものに由来する本質的な情報欠落 なのです。

ラップ作用素の定義と性質

この「主値に押し込める」操作に名前を付けておくと、以後の議論がすっきりします。実数 $x$ に対して ラップ作用素 $\mathcal{W}$ を

$$ \begin{equation} \mathcal{W}(x) = \bmod(x + \pi,\ 2\pi) – \pi \end{equation} $$

と定義します。$\bmod(\cdot, 2\pi)$ は $[0, 2\pi)$ を返す剰余演算です。$x + \pi$ を $2\pi$ で割った余りは $[0, 2\pi)$ に入り、そこから $\pi$ を引くので、結果は必ず $[-\pi, \pi)$ に収まります。等価な表現として、四捨五入 $\operatorname{round}$ を使った

$$ \mathcal{W}(x) = x – 2\pi \operatorname{round}\!\left(\frac{x}{2\pi}\right) $$

も便利です。こちらの形からは、$\mathcal{W}$ が「$x$ に最も近い $2\pi$ の整数倍を引き算する」操作であることが一目で分かります。

このグラフは $\mathcal{W}$ の全体像を1枚に凝縮しています。出力は必ず $[-\pi, \pi)$ の帯に収まり、入力を $2\pi$ ずらしても波形がそのまま繰り返される（＝出力が変わらない）ことが読み取れます。そして濃い緑で塗った $|x| < \pi$ の帯では、赤い実線が灰色の破線 $y = x$ と完全に重なっており、この範囲では $\mathcal{W}$ が何の情報も壊していないことが視覚的に確認できます。

以降で繰り返し使う性質を2つ、証明とともに挙げます。

性質1（$2\pi$ 周期性）: 任意の整数 $m$ に対して $\mathcal{W}(x + 2\pi m) = \mathcal{W}(x)$。

これは定義から直ちに従います。$\bmod((x + 2\pi m) + \pi, 2\pi) = \bmod(x + \pi, 2\pi)$ なので、$\pi$ を引いた結果も同じです。この性質が意味するのは、ラップ作用素は「$2\pi$ の整数倍だけ違う情報」を完全に潰す ということです。

性質2（小さい引数への恒等作用）: $|x| < \pi$ ならば $\mathcal{W}(x) = x$。

$|x| < \pi$ なら $x + \pi \in (0, 2\pi)$ なので剰余演算は何もせず、$\mathcal{W}(x) = (x + \pi) - \pi = x$ です。つまり 絶対値が $\pi$ 未満の量に対しては $\mathcal{W}$ は透明で、何の情報も失いません。

この2つの性質を組み合わせると、後で出てくるアンラップ原理の核心が見えてきます。「$2\pi$ の整数倍のズレを含む量に $\mathcal{W}$ をかけると、ズレだけが消えて、真の量が $\pi$ 未満なら真の量がそのまま残る」。これがアンラップの唯一にして最大の武器です。次節ではこれを問題の定式化に接続します。

問題の定式化：未知整数を決める

観測された位相を $\psi_n$（$n$ は標本番号）、私たちが本当に欲しい連続的な位相を $\phi_n$ とすると、両者の関係は

$$ \begin{equation} \psi_n = \phi_n + 2\pi k_n, \qquad k_n \in \mathbb{Z} \end{equation} $$

と書けます。$\psi_n$ は既知、$\phi_n$ と $k_n$ は未知です。位相アンラップとは、この式から整数列 $\{k_n\}$ を決定して $\phi_n = \psi_n – 2\pi k_n$ を復元する問題です。

ここで即座に気づくべき点があります。この問題は、そのままでは解が無数にあります。$N$ 個の標本に対して未知数は $N$ 個の整数、観測は $N$ 個の実数ですが、観測値はすでに $\psi_n$ という形で与えられており、$k_n$ をどう選んでも $\psi_n – 2\pi k_n$ は式(2)を満たしてしまうからです。言い換えると、$\phi_n$ の候補は各点で独立に $\{\dots, \psi_n – 2\pi, \psi_n, \psi_n + 2\pi, \dots\}$ という格子上を自由に動けます。データだけからは何も決まりません。

灰色の丸が、各標本で式(2)を満たす候補（$2\pi$ 間隔の格子点）です。青い直線と赤い破線はどちらもこの格子点だけを通っており、観測データとの整合性という意味では完全に同格です。それでも私たちが青を選びたくなるのは、赤がジグザグに暴れているから — つまり選択の根拠はデータではなく「なめらかさ」という外から持ち込んだ仮定なのだ、ということがこの図から見て取れます。

したがって、アンラップは 必ず何らかの事前知識（正則化）を必要とする不良設定問題 です。実務で使われる仮定はほぼ例外なく次の一つです。

なめらかさ仮定: 真の位相 $\phi_n$ は隣接標本間で急激には変わらない。具体的には $|\phi_{n+1} – \phi_n| < \pi$。

この仮定は物理的には「十分に細かく標本化している」という主張であり、後で見るように標本化定理そのものです。逆に言えば、この仮定が破れる場面（急峻な地形、位相の不連続、レイオーバー）では、どんなアルゴリズムを持ってきても原理的に正解には到達できません。アンラップアルゴリズムの良し悪しとは、仮定が破れた場所で被害をどこまで局所化できるか の勝負なのです。

問題の形が見えたところで、まずは仮定が素直に効く1次元の場合から始めましょう。

1次元アンラップ：差分をラップして積み直す

直感から入ります。時計の針を1分おきに見ているとします。3時→3時1分と動いたとき、針は少しだけ右に回りました。この「少しだけ」という差分の情報は折り返しの影響を受けていません。だから、絶対的な時刻は分からなくても、差分を順に足し上げていけば経過時間は完全に復元できるわけです。位相アンラップの1次元アルゴリズムは、文字通りこれだけです。

アルゴリズム

1. 隣接標本の観測位相差 $D_n = \psi_{n+1} – \psi_n$ を計算する。これは $[-2\pi, 2\pi]$ の範囲に散らばる。
2. これをラップして $d_n = \mathcal{W}(D_n) \in [-\pi, \pi)$ とする。
3. 積分（累積和）する：$\hat{\phi}_0 = \psi_0$、$\hat{\phi}_{n} = \hat{\phi}_0 + \sum_{m=0}^{n-1} d_m$。

たったこれだけです。ステップ2が「差分は小さいはずだから、$\pi$ を超える見かけの差分は折り返しのせいだ」という仮定を注入している唯一の場所です。

中段が手続きの核心です。生の差分（灰色）は折り返しのたびに $-2\pi$ 近くまで落ち込む深い谷を作りますが、$\mathcal{W}$ をかけた差分（緑）はどこでも $+0.16$ rad 前後の一定値に戻っています。この「本来の小さな位相進み」を積み上げた下段の結果は、真値との差の標準偏差が $10^{-14}$ rad 台、つまり浮動小数点の丸め誤差レベルで真の位相と一致します。

Itohの定理：いつ厳密に復元できるか

上の手続きが本当に真の位相を再現するのか、きちんと確かめておきましょう。

定理（Itoh, 1982）: すべての $n$ について $|\phi_{n+1} – \phi_n| < \pi$ が成り立つならば、

$$ \mathcal{W}(\psi_{n+1} – \psi_n) = \phi_{n+1} – \phi_n $$

が成り立つ。したがって上のアルゴリズムは $\phi_n$ を（全体の定数 $2\pi k_0$ を除いて）厳密に復元する。

証明: 式(2)を差分に代入します。

$$ \psi_{n+1} – \psi_n = (\phi_{n+1} + 2\pi k_{n+1}) – (\phi_n + 2\pi k_n) = (\phi_{n+1} – \phi_n) + 2\pi (k_{n+1} – k_n) $$

$k_{n+1} – k_n$ は整数なので、ラップ作用素の性質1（$2\pi$ 周期性）により、この整数倍の項は $\mathcal{W}$ で消えます。

$$ \mathcal{W}(\psi_{n+1} – \psi_n) = \mathcal{W}(\phi_{n+1} – \phi_n) $$

ここで仮定 $|\phi_{n+1} – \phi_n| < \pi$ を使うと、性質2（小さい引数への恒等作用）により $\mathcal{W}$ は何もしません。

$$ \mathcal{W}(\phi_{n+1} – \phi_n) = \phi_{n+1} – \phi_n $$

以上の2式をつなげば定理の主張が得られます。あとは累積和が望遠鏡和になることを確認すれば十分です。

$$ \hat{\phi}_n – \hat{\phi}_0 = \sum_{m=0}^{n-1} d_m = \sum_{m=0}^{n-1} (\phi_{m+1} – \phi_m) = \phi_n – \phi_0 $$

初期値を $\hat{\phi}_0 = \psi_0 = \phi_0 + 2\pi k_0$ と置いているので、$\hat{\phi}_n = \phi_n + 2\pi k_0$。すなわち全体に共通の $2\pi$ の整数倍だけの不定性を残して一致します。$\square$

この「全体定数の不定性」は多くの応用では無害です。InSARなら別途の高度基準点で、FMCWレーダーならビート周波数から得た粗距離で、絶対値を後から拘束できるからです。

なぜこれは標本化定理と同じことなのか

Itohの条件 $|\phi_{n+1} – \phi_n| < \pi$ を、周波数の言葉に翻訳してみましょう。標本化周波数 $f_s$ で瞬時周波数 $f$ の信号を観測しているなら、1標本あたりの位相進みは

$$ \Delta\phi = 2\pi \frac{f}{f_s} $$

です。これを条件に代入すると

$$ \left| 2\pi \frac{f}{f_s} \right| < \pi \quad \Longleftrightarrow \quad |f| < \frac{f_s}{2} $$

となり、ナイキスト条件そのものが出てきます。つまり位相アンラップとは、周波数領域のエイリアシングを位相領域で言い直したものにほかなりません。

左のグラフでは、推定周波数（オレンジ）が 500 Hz までは真値（青破線）にぴったり重なり、500 Hz を越えた瞬間に $-500$ Hz 側へ跳んで $f – f_s$ の直線に乗り換えます。右のグラフはまったく同じ境界を「1標本あたりの位相進み $\Delta\phi$ が $\pi$ を超えるか」という位相の言葉で描いたもので、緑と赤の境目が左図の折り返し点と同じ 500 Hz に来ています。2つの条件が別物ではなく同一の不等式の2つの表現であることが、この対比から読み取れます。条件が破れると、ナイキスト周波数とエイリアシングで見た折り返しとまったく同じ形で、真の周波数 $f$ が $f – f_s$ に化けます。この対応は後でPythonで実測します。

誤差伝搬：一度の失敗が永久に残る

累積和という構造には、見逃せない弱点があります。ある1点 $m$ でラップの判定を誤って $d_m$ が真値より $2\pi$ ずれたとすると、その誤差は $\hat{\phi}_{m+1}$ 以降のすべての標本に $2\pi$ のオフセットとして残り続けます。局所的な誤りが大域的な誤りに増幅されるのです。

$$ \hat{\phi}_n = \phi_n + 2\pi k_0 + 2\pi \varepsilon \cdot \mathbb{1}[n > m] $$

これは1次元アンラップの構造的な性質で、どんな工夫をしても「積分する」という枠組みの中では避けられません。この性質は2次元に拡張したときに、さらに厄介な形で顔を出します。

まずは、この失敗がどのくらいの確率で起きるのかを、雑音の観点から評価しておきましょう。

雑音があると位相はどれだけ暴れるか

実際の観測では、複素振幅は $z = A e^{j\phi} + \nu$ の形をしています。$\nu$ は循環対称な複素ガウス雑音で、分散を $\sigma_\nu^2$（実部・虚部それぞれ $\sigma_\nu^2/2$）とします。SNRを $\gamma = A^2/\sigma_\nu^2$ と定義しましょう。

高SNRのとき、位相誤差は次のように見積もれます。信号ベクトル $A e^{j\phi}$ に対して、雑音を「信号と同じ向きの成分」と「直交する成分」に分解します。位相を回すのは直交成分だけで、その分散は $\sigma_\nu^2/2$ です。誤差角が小さければ

$$ \delta\phi \approx \frac{\nu_\perp}{A}, \qquad \sigma_\phi \approx \frac{\sigma_\nu/\sqrt{2}}{A} = \frac{1}{\sqrt{2\gamma}} $$

$\gamma = 10$ dB（$\gamma = 10$）なら $\sigma_\phi \approx 0.224$ rad $= 12.8^\circ$、$\gamma = 20$ dB なら $\sigma_\phi \approx 0.0707$ rad です。実際にモンテカルロで測ると 20 dB で 0.0709 rad、10 dB で 0.230 rad と、高SNR側ではこの近似がよく当たります。一方 0 dB では近似値 0.707 rad に対して実測 0.869 rad と外れます。SNRが下がると位相は $[-\pi, \pi)$ 全体に広がっていき、そもそもガウス分布では近似できなくなるからです。

隣接2標本の位相差の揺らぎは、雑音が独立なら $\sigma_\Delta = \sqrt{2}\,\sigma_\phi$ です。アンラップが失敗するのは、この揺らぎが真の位相差を $\pm\pi$ の外まで押し出したときです。素朴にガウス近似すると失敗確率は $2Q(\pi/\sigma_\Delta)$ 程度になりそうですが、この見積もりは実際より何桁も小さく出ます。

6 dB で6桁もの乖離です。原因は、位相誤差の真の分布が「中心のガウス的な山」と「雑音が信号を上回った標本による一様に近い裾」の混合になっていることにあります。裾の重さが $\pi$ を超える事象の確率を支配するので、ガウス近似は使い物になりません。アンラップの失敗確率を見積もるときは、必ず実測かちゃんとした位相分布（ライス位相分布）で評価する べきです。

左の対数目盛のヒストグラムを見ると、SNR 20 dB では実測（青）が破線のガウス曲線とぴったり重なるのに、0 dB では中心付近こそ山になっているものの $\pm\pi$ の端まで確率密度が $10^{-2}$ 程度で寝そべっており、ガウスの破線が急降下するのと対照的です。この寝そべった裾が $\pi$ を超える事象を支配するため、右のグラフのように 8 dB では実測がガウス近似の $3\times10^{9}$ 倍にもなります。ガウス近似は SNR が高いほど当てにならなくなるという、直感に反する挙動がはっきり見えます。

この裾の重さこそが、実務で「SNRが十分に見えるのに位相がときどき飛ぶ」現象の正体です。では実際に、1次元アンラップを動かして確かめてみましょう。

Pythonでの実装1：チャープ信号の位相と群遅延

最初の実験では、線形チャープ（周波数が時間とともに直線的に上がる信号）を題材にします。チャープはFMCWレーダーの送信波そのものであり、位相が二次関数になるので理論値との比較が容易です。まずラップ作用素と1次元アンラップを自前で実装します。

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

def wrap(x):
    """ラップ作用素 W(x) = mod(x+pi, 2pi) - pi。結果は [-pi, pi)"""
    return np.mod(x + np.pi, 2 * np.pi) - np.pi

def unwrap_1d(psi):
    """差分をラップして積分する1次元アンラップ"""
    d = wrap(np.diff(psi))            # 隣接差分をラップ
    return np.concatenate([[psi[0]], psi[0] + np.cumsum(d)])

# 線形チャープ: f0 -> f1 に 1 秒で掃引
fs, T, f0, f1 = 2000.0, 1.0, 50.0, 450.0
t = np.arange(0, T, 1 / fs)
k = (f1 - f0) / T                     # チャープ率 [Hz/s]
phi_true = 2 * np.pi * (f0 * t + 0.5 * k * t**2)   # 真の位相
psi = wrap(phi_true)                  # 観測できるのはこれだけ

phi_hat = unwrap_1d(psi)
print("自作アンラップと真値の差の標準偏差:", np.std(phi_hat - phi_true))
print("numpy.unwrap との最大差:", np.max(np.abs(phi_hat - np.unwrap(psi))))

真の位相 $\phi$ は $0$ から $2\pi(50 + 200) \approx 1571$ rad まで単調に増加します。ラップされた $\psi$ は $[-\pi,\pi)$ の間を250回近く往復するノコギリ波になりますが、差分をラップして積分し直すと真の位相と完全に一致します（差の標準偏差は数値誤差レベル）。numpy.unwrap とも一致するので、実装の正しさが確認できます。ここで大事なのは、折り返しで一見めちゃくちゃになったノコギリ波から、1571 rad という大きな値が誤差ゼロで復元されている という点です。差分だけで全体が決まるという構造の強さがよく分かります。

瞬時周波数の推定とItoh条件の破れ

次に、アンラップした位相を微分して瞬時周波数を求め、さらに標本化周波数を変えてItoh条件が破れる様子を見ます。

import numpy as np

# --- 瞬時周波数の推定 ---
f_inst = np.diff(phi_hat) / (2 * np.pi) * fs      # アンラップ後の位相を微分
f_naive = np.diff(psi) / (2 * np.pi) * fs         # アンラップせずに微分（誤り）
m = (t[:-1] > 0.05) & (t[:-1] < 0.95)
A = np.vstack([t[:-1][m], np.ones(m.sum())]).T
slope, intercept = np.linalg.lstsq(A, f_inst[m], rcond=None)[0]
print(f"チャープ率  真値 {k:.1f} Hz/s -> 推定 {slope:.2f} Hz/s")
print(f"開始周波数  真値 {f0:.1f} Hz  -> 推定 {intercept:.2f} Hz")
print(f"アンラップなしの平均瞬時周波数: {f_naive[m].mean():.2f} Hz（真値 {f_inst[m].mean():.2f} Hz）")

# --- Itoh条件の破れ: 単一正弦波の周波数を変えて推定 ---
fs2 = 1000.0
t2 = np.arange(0, 1, 1 / fs2)
print("\n真の周波数 | 1標本あたりΔφ | アンラップ推定")
for f in [100, 300, 480, 520, 700]:
    p = wrap(2 * np.pi * f * t2)
    est = np.polyfit(t2, unwrap_1d(p), 1)[0] / (2 * np.pi)
    print(f"  {f:4d} Hz  |  {2*np.pi*f/fs2:.3f} rad  |  {est:8.2f} Hz")

出力を見ると、チャープ率は真値 400.0 Hz/s に対して推定 400.00 Hz/s、開始周波数は 50.0 Hz に対して 50.10 Hz と、ほぼ完璧に一致します。一方アンラップせずに位相を微分すると、平均瞬時周波数は $-0.04$ Hz になってしまいました。ノコギリ波の急落部分が巨大な負のスパイクを作り、上りの寄与とほぼ打ち消し合うためです。位相微分は必ずアンラップ後に行う という鉄則が数値で確認できます。

後半の表はさらに示唆的です。$f_s = 1000$ Hz のとき、480 Hz（$\Delta\phi = 3.016 < \pi$）までは完璧に推定できるのに、520 Hz（$\Delta\phi = 3.267 > \pi$）になった途端に推定値が $-480$ Hz に化けます。700 Hz は $-300$ Hz になります。これは $f \to f – f_s$ という、まさに標本化定理の折り返しです。Itoh条件は工学的なヒューリスティックではなく、ナイキスト条件の言い換えなのだと納得できます。

周波数軸のアンラップ：群遅延の推定

アンラップが必要になるのは時間軸だけではありません。信号のフーリエ変換 $X(f) = |X(f)| e^{j\Phi(f)}$ の位相 $\Phi(f)$ をアンラップして微分すると、群遅延

$$ \tau_g(f) = -\frac{1}{2\pi} \frac{d\Phi(f)}{df} $$

が得られます。線形チャープの場合、定常位相法から $\Phi(f) \approx -\pi(f – f_0)^2/k – \pi/4$ なので、理論群遅延は $\tau_g(f) = (f – f_0)/k$、つまり「その周波数が送信された時刻」です。ここでItoh条件は周波数軸上の条件になり、$|d\Phi/df| \cdot \Delta f < \pi$ すなわち $\tau_g \Delta f < 1/2$ を要求します。周波数分解能 $\Delta f$ が粗いとアンラップが破綻するわけです。

import numpy as np

N = len(t)
tau_theory = lambda f: (f - f0) / k
print("ゼロ詰め率 | 周波数分解能 | 群遅延のRMS誤差")
for pad in [1, 2, 4, 8]:
    Nf = N * pad
    X = np.fft.rfft(np.cos(2 * np.pi * (f0 * t + 0.5 * k * t**2)) * np.hanning(N), Nf)
    f = np.fft.rfftfreq(Nf, 1 / fs)
    Phi = unwrap_1d(np.angle(X))                  # 周波数軸に沿ってアンラップ
    tau = -np.gradient(Phi, f) / (2 * np.pi)      # 群遅延
    sel = (f > 100) & (f < 400)
    rms = np.sqrt(np.mean((tau[sel] - tau_theory(f[sel]))**2))
    dphi_max = 2 * np.pi * tau_theory(400.0) * (f[1] - f[0])
    print(f"   x{pad:<5d} | {f[1]-f[0]:6.3f} Hz  | {rms:.5f} s  (bin間の最大Δφ = {dphi_max:.2f} rad)")

結果は劇的です。ゼロ詰めなし（$\Delta f = 1$ Hz）では群遅延のRMS誤差が 0.707 秒、つまり測定区間の長さ 0.75 秒とほぼ同じで、完全に無意味な出力です。ところが2倍にゼロ詰めして $\Delta f = 0.5$ Hz にすると、誤差は $4 \times 10^{-5}$ 秒に激減し、4倍・8倍にしても改善しません。

境目の理由もはっきりしています。400 Hz における理論群遅延は 0.875 秒なので、bin間の位相変化は $\Delta f = 1$ Hz で $2\pi \times 0.875 \times 1 = 5.50$ rad となり $\pi$ を大きく超えます。$\Delta f = 0.5$ Hz なら 2.75 rad で、ぎりぎり $\pi = 3.14$ を下回ります。ゼロ詰めは新しい情報を足さないのに結果を劇的に改善する という一見不思議な現象は、それが周波数軸の標本化を細かくしてItoh条件を満たさせているからだ、と理解できます。

左のグラフで、ゼロ詰めなし（赤）の群遅延は 250 Hz 付近で $+0.5$ 秒から $-0.5$ 秒へ落下し、そこから理論値（黒破線）とは無関係な直線を描いています。これは周波数軸上でItoh条件が破れた結果で、$\times 2$（オレンジ）は黒破線に完全に乗っています。右の棒グラフは $\Delta\phi$ が $\pi$ を下回った $\times 2$ で誤差が4桁落ち、それ以上細かくしても頭打ちになることを示しており、必要なのは「たくさんゼロ詰めすること」ではなく「$\pi$ を切ること」だと分かります。

ここまでは1次元、しかも「積分経路が一本しかない」世界の話でした。2次元に移った瞬間、この牧歌的な状況は崩壊します。

2次元の壁：経路依存性とレジデュー

干渉SARの位相画像やホログラフィの干渉縞は2次元の配列です。1次元アルゴリズムを素直に拡張しようとすると、「まず1列目を縦にアンラップして、各行を横にアンラップする」といったラスタ走査を考えることになります。しかしここで根本的な問いが生まれます — どの経路を通ってその画素に到達するかで、答えは変わらないのか？

1次元では経路は一本しかないので、この問いは存在しませんでした。2次元では左回りにも右回りにも行けます。もし答えが経路に依存するなら、そもそも「正しいアンラップ結果」という概念自体が怪しくなります。

巡回積分とレジデューの定義

答えを出すために、最小の閉ループ、すなわち隣り合う4画素からなる $2\times2$ のループを考えます。画素 $(i,j) \to (i,j+1) \to (i+1,j+1) \to (i+1,j) \to (i,j)$ と一周する間に、ラップした差分を足し上げます。

$$ \begin{aligned} \Delta_1 &= \mathcal{W}(\psi_{i,j+1} – \psi_{i,j}) \\ \Delta_2 &= \mathcal{W}(\psi_{i+1,j+1} – \psi_{i,j+1}) \\ \Delta_3 &= \mathcal{W}(\psi_{i+1,j} – \psi_{i+1,j+1}) \\ \Delta_4 &= \mathcal{W}(\psi_{i,j} – \psi_{i+1,j}) \end{aligned} $$

この4つの和を $2\pi$ で割った量

$$ \begin{equation} q_{i,j} = \frac{1}{2\pi}\left(\Delta_1 + \Delta_2 + \Delta_3 + \Delta_4\right) \end{equation} $$

を、そのループの レジデュー（residue、残差点） と呼びます。

$q_{i,j}$ が整数になることを確認しましょう。ラップ前の差分 $\psi_{i,j+1} – \psi_{i,j}$ などをそのまま4つ足すと、一周して戻るので恒等的に $0$ です。各 $\Delta$ はラップ前の差分から $2\pi$ の整数倍を引いたものなので、和は $2\pi \times (\text{整数})$ になります。したがって $q_{i,j} \in \mathbb{Z}$。さらに各 $\Delta \in [-\pi, \pi)$ なので和は $[-4\pi, 4\pi)$ に入り、$q_{i,j} \in \{-2, -1, 0, +1\}$ に限られます。$|q| = 2$ は4辺すべてが同時に $\pm\pi$ 近くになる必要があり、実際にはまず起きません。実験でも $-5$ dB という劣悪なSNRまで下げても $|q| \le 1$ しか観測されませんでした。そこでレジデューは通常 電荷 $+1$（正のレジデュー）、$-1$（負のレジデュー）、$0$（レジデューなし） の3値として扱います。

右の図は、左の分布から実際のレジデューを1つ取り出して4辺を一周したものです。各辺のラップ差分はどれも $[-\pi, \pi)$ に収まった「もっともらしい」値なのに、$0.422 + 1.778 + 2.216 + 1.867 = 6.283$ となって $2\pi$ ちょうどになっています。閉ループを一周して元の点に戻ったのに $2\pi$ 増えてしまう — この局所的な矛盾がレジデューの正体で、左図ではそれが位相の乱れた中央部に集中して現れていることが分かります。

レジデューが意味すること

$q_{i,j} = 0$ なら、その $2\times2$ ループを一周してもラップ差分の和はゼロ、つまり どちら回りに行っても同じ答え になります。$q_{i,j} \neq 0$ なら、回り方を変えると $2\pi q_{i,j}$ だけ違う答えが出ます。

この局所的な性質は、大きなループにそのまま持ち上がります。任意の閉曲線 $C$ に沿った巡回積分は、$C$ が囲む領域内のレジデューの総和に等しくなります。

$$ \begin{equation} \frac{1}{2\pi}\oint_C \mathcal{W}(\nabla\psi) \cdot d\bm{l} = \sum_{(i,j) \in \operatorname{int}(C)} q_{i,j} \end{equation} $$

証明の骨子: 領域内部を $2\times2$ ループで敷き詰めます。隣り合う2つのループが共有する辺は、一方のループでは順方向、他方では逆方向に横断されます。$\mathcal{W}(-x) = -\mathcal{W}(x)$（$x = \pi$ ちょうどという零集合を除く）なので、内部の辺の寄与はすべて相殺し、外周の辺の寄与だけが残ります。これは連続系におけるStokesの定理（$\oint \bm{F}\cdot d\bm{l} = \iint \nabla\times\bm{F}\, dS$）の離散版で、レジデューはラップ勾配場の「回転（curl）」に対応します。$\square$

ここから重要な帰結が2つ出ます。

帰結1: 画像全体にレジデューが1つも無ければ、アンラップ結果は経路に依らず一意に定まる。これは「回転がゼロのベクトル場はポテンシャル場である」という事実の離散版です。

帰結2: レジデューの総電荷は、画像外周を一周した巡回積分に等しい。外周上でItoh条件が満たされていればこれはゼロなので、正のレジデューと負のレジデューは必ず同数現れます。実験でも、SNR 8 dB で $+31$ 個・$-31$ 個、SNR 0 dB で $+718$ 個・$-717$ 個（境界の影響で微差）と、見事に対で発生します。物理的には正負の電荷対、あるいは渦・反渦の対と見なせます。

レジデューはなぜ現れるのか

原因は2つあります。

• 雑音: 位相の揺らぎが $\pi$ を超えると、局所的にItoh条件が破れてレジデューになります。前節の表のとおり、SNRが下がると急激に増えます。実験では 20 dB で 1 個、15 dB で 7 個、10 dB で 38 個、0 dB で 1371 個（$128\times128$ 画像、$127\times127 = 16129$ ループ中）でした。
• 位相の急峻さ（エイリアシング）: 雑音がゼロでも、真の位相勾配が $\pi$ を超えるとレジデューが出ます。ガウス丘の高さを変えて実験すると、最大勾配 2.96 rad（$<\pi$）ではレジデュー 0 個で完全復元できるのに、勾配 4.27 rad にした途端に 62 個のレジデューが発生し、ラスタ走査アンラップのRMS誤差は 19.8 rad まで悪化します。SARにおけるレイオーバーや急斜面がこれに当たります。

この3枚は雑音を一切加えずに、丘の高さだけを変えたものです。最大勾配 2.96 rad の左端はレジデューが 0 個でラスタ走査の誤差もゼロですが、$\pi$ を超えた中央（3.48 rad）で 26 個、右端（4.27 rad）で 62 個のレジデューが縞の詰まった部分に湧き出し、誤差はそれぞれ 9.42 rad、19.81 rad まで悪化します。レジデューは雑音だけの産物ではなく、位相そのものが急峻すぎるときにも必ず生じることが、この対比からはっきり分かります。

つまりレジデューは「アンラップが失敗しかけている場所」を教えてくれる 診断量 であり、同時に「単純な経路積分では手に負えない」という 警告 でもあります。次に見るべきは、この警告にどう対処するかです。

Goldsteinの枝切り法

レジデューがある場所では経路積分の答えが経路に依存する — ならば、答えが変わってしまうような経路を最初から禁止すればよい。これがGoldstein、Zebker、Wernerが1988年に提案した 枝切り法（branch cut algorithm） の発想です。複素解析で多価関数を扱うときに分岐切断（branch cut）を入れて一価にするのと、まったく同じアイデアです。

手順は次のとおりです。

1. すべての $2\times2$ ループのレジデュー $q_{i,j}$ を計算する。
2. 正のレジデューと負のレジデューを線分（枝切り）で結び、各枝切りグループの総電荷がゼロになるようにする。電荷が余る場合は画像境界に接続して中和する。
3. 積分経路が枝切りを横切ることを禁止したうえで、通常のフラッドフィル（塗りつぶし）で全画素をアンラップする。

ステップ2でなぜ「総電荷ゼロ」が要るのかは、式(4)から分かります。枝切りを横切らない任意の閉ループは、正負が打ち消し合った電荷しか囲めないため、巡回積分がゼロになります。つまり 枝切りを入れた領域の中では、ラップ勾配場は回転ゼロの場に戻るわけです。これでステップ3の経路積分は一意に定まります。

Goldsteinのオリジナルアルゴリズムでは、ステップ2を貪欲法で実現します。未処理のレジデューを1つ選び、その周囲の $3\times3$ の探索窓の中で最も近い反対電荷のレジデューを探して結びます。見つからなければ窓を $5\times5$、$7\times7$ と広げていきます。窓が境界に達したら境界に接続して中和します。枝切りの総延長をできるだけ短くする という方針は、「レジデューは対で近くに生じやすいので、近いもの同士を結ぶのが誤りの局在化に最も効く」という経験則に基づいています。

枝切り法の長所と短所を整理しておきます。

• 長所: レジデューが無い（＝データが健全な）領域では、真の位相を厳密に復元します。誤差が枝切りの近傍に閉じ込められ、健全な領域に漏れ出しません。位相の不連続を保存できます。
• 短所: レジデューが密集すると枝切りが絡み合い、フラッドフィルが到達できない孤立領域（アイランド）が生じます。この領域はアンラップ不能として欠測扱いになります。また枝切りの張り方は貪欲法であって最適ではなく、結果がレジデューを処理する順序に依存します。

「どうしても答えを出さないといけない、欠測は許されない」という要求に対しては、まったく別の設計思想が必要です。それが次節の最小二乗アプローチです。

最小二乗アンラップ：ポアソン方程式に帰着させる

枝切り法は「レジデューを避けて通る」戦略でした。最小二乗アンラップは正反対に、レジデューの存在を認めたうえで、矛盾を全体に薄く均して丸ごと飲み込む 戦略を取ります。

考え方はシンプルです。ラップされた勾配 $\mathcal{W}(\nabla\psi)$ は、レジデューがある場所では「どんな滑らかな関数の勾配でもない」矛盾したベクトル場です。ならば、その矛盾場に 最も近い「本物の勾配場」を持つ関数 $\phi$ を探せばよい。

目的関数と正規方程式

観測位相 $\psi_{i,j}$（$M \times N$ 格子）から、ラップした前進差分を作ります。

$$ \Delta^x_{i,j} = \mathcal{W}(\psi_{i,j+1} – \psi_{i,j}), \qquad \Delta^y_{i,j} = \mathcal{W}(\psi_{i+1,j} – \psi_{i,j}) $$

境界の外では $\Delta^x = \Delta^y = 0$ と定義します。求める $\phi_{i,j}$ の勾配がこれらにできるだけ近くなるよう、次の $L^2$ 誤差を最小化します。

$$ \begin{equation} E(\phi) = \sum_{i,j}\left(\phi_{i,j+1} – \phi_{i,j} – \Delta^x_{i,j}\right)^2 + \sum_{i,j}\left(\phi_{i+1,j} – \phi_{i,j} – \Delta^y_{i,j}\right)^2 \end{equation} $$

これは $\phi$ について二次形式なので、偏微分をゼロと置けば線形方程式（正規方程式）が得られます。$\phi_{i,j}$ が現れる項は4つあることに注意してください — 水平方向では $(i,j)$ 番目の項に $-\phi_{i,j}$ として、$(i,j-1)$ 番目の項に $+\phi_{i,j}$ として現れ、垂直方向も同様です。それぞれを微分すると

$$ \frac{1}{2}\frac{\partial E}{\partial \phi_{i,j}} = -\left(\phi_{i,j+1} – \phi_{i,j} – \Delta^x_{i,j}\right) +\left(\phi_{i,j} – \phi_{i,j-1} – \Delta^x_{i,j-1}\right) -\left(\phi_{i+1,j} – \phi_{i,j} – \Delta^y_{i,j}\right) +\left(\phi_{i,j} – \phi_{i-1,j} – \Delta^y_{i-1,j}\right) $$

これをゼロと置き、$\phi$ の項を左辺、$\Delta$ の項を右辺に集めます。$\phi$ の項をまとめると $-\phi_{i,j+1} – \phi_{i,j-1} – \phi_{i+1,j} – \phi_{i-1,j} + 4\phi_{i,j}$ となるので、全体の符号を反転させて

$$ \begin{equation} \underbrace{\phi_{i,j+1} + \phi_{i,j-1} + \phi_{i+1,j} + \phi_{i-1,j} – 4\phi_{i,j}}_{\text{5点離散ラプラシアン } \nabla^2\phi} = \underbrace{\left(\Delta^x_{i,j} – \Delta^x_{i,j-1}\right) + \left(\Delta^y_{i,j} – \Delta^y_{i-1,j}\right)}_{\rho_{i,j}\ =\ \operatorname{div}\mathcal{W}(\nabla\psi)} \end{equation} $$

最小二乗アンラップは離散ポアソン方程式 $\nabla^2 \phi = \rho$ を解くことと等価 だと分かりました。右辺 $\rho$ は「ラップ勾配場の発散」です。連続系の言葉なら、$\phi$ は電位、$\rho$ は電荷分布に対応します。

ここで前節の議論が効いてきます。ヘルムホルツ分解によれば、任意のベクトル場は「回転のない成分（勾配場）」と「発散のない成分（回転場）」に分けられます。$\rho$ は発散だけを拾っているので、最小二乗解はラップ勾配場の勾配成分だけを取り出し、レジデュー（＝回転成分）を完全に捨てているのです。だからこの方法は絶対に失敗しません — 矛盾は最初から解の外に置かれるからです。同時に、捨てられた矛盾のぶん、誤差が画像全体に薄く広がるという代償を負います。この性質は後で数値実験ではっきり見えます。

境界条件

境界の外側で $\Delta = 0$ と置いたことは、境界で「差分がゼロ」すなわち $\partial\phi/\partial n = 0$ を課したことに相当します。これは ノイマン境界条件 です。ノイマン問題は解が定数の加算に対して不定なので、$\phi$ の平均をゼロと決めて一意化します。もともとアンラップの答えには全体定数の不定性があるので、これは実害になりません。

DCTによる高速解法

$M \times N$ 画像なら未知数は $MN$ 個です。$128\times128$ でも 16384 元連立方程式で、素朴に解くのは非現実的です。しかし係数行列は5点ラプラシアンという極めて構造的な行列なので、適切な直交変換で対角化できます。それが 離散コサイン変換（DCT） です。

なぜDCTなのかを1次元で確認します。ノイマン条件下のラプラシアンの固有ベクトル候補として

$$ c_i^{(m)} = \cos\!\left(\frac{\pi(2i+1)m}{2M}\right), \qquad i = 0,\dots,M-1 $$

を取ります。これはDCT-IIの基底関数そのものです。差分を計算しましょう。三角関数の和積公式 $\cos(A+B) + \cos(A-B) = 2\cos A \cos B$ を、$A = \pi(2i+1)m/(2M)$、$B = \pi m/M$ として使うと

$$ c_{i+1}^{(m)} + c_{i-1}^{(m)} = \cos\!\left(\frac{\pi(2i+3)m}{2M}\right) + \cos\!\left(\frac{\pi(2i-1)m}{2M}\right) = 2\cos\!\left(\frac{\pi m}{M}\right) c_i^{(m)} $$

したがって

$$ c_{i+1}^{(m)} – 2c_i^{(m)} + c_{i-1}^{(m)} = \left(2\cos\frac{\pi m}{M} – 2\right) c_i^{(m)} $$

となり、$c^{(m)}$ は固有値 $2\cos(\pi m/M) – 2$ の固有ベクトルです。境界（$i=0$ と $i=M-1$）でも、DCT-IIが暗黙に行う半標本ずらしの鏡像拡張のおかげでノイマン条件が自動的に満たされ、同じ関係が成り立ちます。

2次元では変数分離により、$x$ 方向と $y$ 方向の固有値が足し算になります。

$$ \begin{equation} \lambda_{m,n} = 2\cos\!\left(\frac{\pi m}{M}\right) + 2\cos\!\left(\frac{\pi n}{N}\right) – 4 \end{equation} $$

これで解法が完成します。$\rho$ の2次元DCTを $\hat{\rho}_{m,n}$ とすると、ポアソン方程式は変換領域で単なる割り算になります。

$$ \hat{\phi}_{m,n} = \frac{\hat{\rho}_{m,n}}{\lambda_{m,n}}, \qquad \hat{\phi}_{0,0} = 0 $$

$\lambda_{0,0} = 2 + 2 – 4 = 0$ はまさに定数モード（ノイマン問題の零空間）に対応するので、$\hat{\phi}_{0,0} = 0$ と置くことで平均ゼロ解を選びます。最後に逆DCTで $\phi$ に戻します。計算量は $O(MN\log MN)$ で、$1000\times1000$ の画像でも一瞬です。

重み付き最小二乗

最小二乗法の弱点は、信頼できないデータと信頼できるデータを同列に扱ってしまう ことです。InSARなら水面や植生でコヒーレンスが低く、そこの位相はほぼ雑音です。それを健全な領域と同じ重みで拘束すると、雑音領域の矛盾が画像全体に染み出します。

そこで、各差分に重み $w$ を付けます。

$$ \begin{equation} E_w(\phi) = \sum_{i,j} w^x_{i,j}\left(\phi_{i,j+1} – \phi_{i,j} – \Delta^x_{i,j}\right)^2 + \sum_{i,j} w^y_{i,j}\left(\phi_{i+1,j} – \phi_{i,j} – \Delta^y_{i,j}\right)^2 \end{equation} $$

重みにはコヒーレンス $\gamma$ の関数（$w = \gamma^2$ など）や、レジデュー近傍をゼロにするマスクを使います。同じように偏微分をゼロと置くと、正規方程式は

$$ -\operatorname{div}\left(w\,\nabla\phi\right) = -\operatorname{div}\left(w\,\mathcal{W}(\nabla\psi)\right) $$

という 係数が空間変化するポアソン方程式 になります。係数が一定でないので、もはやDCTでは対角化できません。

ここでGhigliaとRomero（1994）の実用的な解法が効いてきます。この方程式の係数行列は対称正定値（定数モードを除けば半正定値）なので 共役勾配法（CG） が使えます。しかも、重みを外した定数係数版（＝DCTで一発で解ける問題）を前処理行列として使う と、収束が劇的に速くなります。重みが $0$ と $1$ の間で緩やかに変化する限り、重み付き作用素は非重み付き作用素と「近い」ので、良い前処理になるわけです。実験では $128\times128$ の問題が11回のPCG反復で残差 $10^{-9}$ まで収束しました。1回の反復はDCT2回ぶんなので、非重み付き解法の20倍程度のコストで済みます。

なお、$L^2$ ノルムの代わりに $L^1$ ノルム（$\sum w|\nabla\phi – \mathcal{W}(\nabla\psi)|$）を最小化する定式化もあり、これは最小費用流問題に帰着します（Costantiniのネットワークフロー法）。$L^1$ 解は誤差を薄く広げず不連続を保存するので品質は高い一方、計算コストは跳ね上がります。InSARの標準ツールSNAPHUは、この系列に統計的コスト関数を組み合わせたものです。

理論が揃ったので、2次元アンラップを実際に動かして比較しましょう。

Pythonでの実装2：2次元のラップ位相とレジデュー

まず、真の位相として「一様な傾斜＋中央のガウス丘」を作ります。傾斜は干渉SARにおける地球楕円体成分や軌道縞、ガウス丘は局所的な地形隆起や地殻変動に対応する典型的なモデルです。

import numpy as np

def true_phase(M=128, N=128, amp=45.0):
    """真の位相：一様傾斜 + 中央のガウス丘"""
    y, x = np.mgrid[0:M, 0:N]
    X, Y = (x - N / 2) / (N / 2), (y - M / 2) / (M / 2)
    return 12 * np.pi * X + 8 * np.pi * Y + amp * np.exp(-(X**2 + Y**2) / (2 * 0.18**2))

phi2 = true_phase()
grad_max = max(np.abs(np.diff(phi2, axis=1)).max(), np.abs(np.diff(phi2, axis=0)).max())
print(f"位相のレンジ: {phi2.ptp():.1f} rad = {phi2.ptp()/(2*np.pi):.1f} 縞")
print(f"最大勾配: {grad_max:.3f} rad  (Itoh条件の上限 pi = {np.pi:.3f})")

def residues(psi):
    """2x2ループのレジデュー（電荷）を計算。戻り値は (M-1, N-1) の整数配列"""
    a = wrap(psi[:-1, 1:] - psi[:-1, :-1])    # 右へ
    b = wrap(psi[1:,  1:] - psi[:-1,  1:])    # 下へ
    c = wrap(psi[1:, :-1] - psi[1:,   1:])    # 左へ
    d = wrap(psi[:-1, :-1] - psi[1:,  :-1])   # 上へ
    return np.round((a + b + c + d) / (2 * np.pi)).astype(int)

psi2 = wrap(phi2)
r0 = residues(psi2)
print(f"無雑音でのレジデュー数: {np.abs(r0).sum()}")

rng = np.random.default_rng(11)
snr = 10 ** (6 / 20)          # SNR 6 dB
z = np.exp(1j * phi2) + (rng.standard_normal(phi2.shape)
                         + 1j * rng.standard_normal(phi2.shape)) / (np.sqrt(2) * snr)
psi_n = np.angle(z)
r = residues(psi_n)
print(f"SNR 6 dB でのレジデュー: 正 {(r>0).sum()} 個, 負 {(r<0).sum()} 個, 総電荷 {r.sum()}")

位相のレンジは 124.7 rad、約 19.8 本の干渉縞に相当します。最大勾配は 2.957 rad で $\pi$ をわずかに下回っているため、無雑音ではレジデューが 1 個も出ません。この設定は「雑音さえなければ完璧に解ける、ぎりぎり健全な」データです。ここに SNR 6 dB の雑音を加えると、正のレジデューと負のレジデューがきっちり同数（総電荷ゼロ）現れます。これは帰結2で証明した保存則の実測確認になっています。

左端のなめらかな真の位相（レンジ 124.7 rad）が、中央ではラップによって約 19.8 本の干渉縞に化けています。中央と右端は同じ位相なのに、SNR 6 dB の雑音が乗った右端では縞のエッジがざらつき、それだけで 122 個のレジデューが生まれます。見た目の劣化は軽微に見えても、$2\times2$ ループ単位ではこれだけの矛盾が仕込まれている、というのがこの比較の要点です。

経路依存性を実測する

レジデューが本当に経路依存性を生むのか、2つの経路で同じ画素に到達して確かめます。

import numpy as np

def path_row_first(psi, i, j):
    """(0,0) -> (0,j) -> (i,j) の経路で積分（先に横、次に縦）"""
    return psi[0, 0] + wrap(np.diff(psi[0, :j+1])).sum() + wrap(np.diff(psi[:i+1, j])).sum()

def path_col_first(psi, i, j):
    """(0,0) -> (i,0) -> (i,j) の経路で積分（先に縦、次に横）"""
    return psi[0, 0] + wrap(np.diff(psi[:i+1, 0])).sum() + wrap(np.diff(psi[i, :j+1])).sum()

# 矩形 [0:i, 0:j] が囲む総電荷を累積和で求める
C = np.cumsum(np.cumsum(r, axis=0), axis=1)
print(f"経路で答えが変わる画素の割合: {100*np.mean(C != 0):.1f} %")

for (i, j) in [(8, 121), (63, 61), (68, 46), (127, 92)]:
    a, b = path_row_first(psi_n, i, j), path_col_first(psi_n, i, j)
    print(f"画素({i:3d},{j:3d})  経路A={a:8.3f}  経路B={b:8.3f}  "
          f"差/2pi={(a-b)/(2*np.pi):+.3f}  囲む総電荷={C[i-1, j-1]:+d}")

出力は式(4)の離散Stokesの定理を寸分違わず再現します。たとえば画素 $(63,61)$ では経路A が 88.811 rad、経路B が 82.528 rad で、差はぴったり $2\pi$。そしてその矩形が囲むレジデューの総電荷はちょうど $+1$ です。他の画素でも「差 $/2\pi$ ＝ 囲む総電荷」が完全に一致します。

さらに衝撃的なのは、全画素の約15%で経路によって答えが変わる という事実です。SNR 6 dB という決して劣悪ではない条件で、$16129$ ループ中レジデューはわずか 122 個しかないにもかかわらず、その影響は画像の広い範囲に及びます。レジデューが1つあれば、それを分断する境界の「向こう側」全体に $2\pi$ の段差が入るからです。単純なラスタ走査アンラップが実データで惨敗する理由が、これで腑に落ちます。

左は「先に横、次に縦」（水色）と「先に縦、次に横」（黄）という2本の経路で同じ星印の画素に到達したもので、答えは 88.811 rad と 82.528 rad、差はぴったり $2\pi$ です。右は全画素について「左上からその画素までの矩形が囲む総電荷」を色で示したもので、赤や青に染まった領域が経路によって答えが変わる場所を表します。レジデューはたった 122 個なのに、そこから水平・垂直に帯状の汚染が広がって全画素の 14.8% を覆っている点に注目してください。

Pythonでの実装3：DCT最小二乗アンラップとSNR掃引比較

いよいよ最小二乗アンラップを実装します。式(6)のポアソン方程式を、式(7)の固有値を使ってDCT領域で解きます。

import numpy as np
from scipy.fft import dct, idct

def dct2(a):
    return dct(dct(a, axis=0, norm='ortho'), axis=1, norm='ortho')

def idct2(a):
    return idct(idct(a, axis=0, norm='ortho'), axis=1, norm='ortho')

def ls_unwrap(psi):
    """DCTによる非重み付き最小二乗アンラップ（Ghiglia & Romero 1994）"""
    M, N = psi.shape
    # ラップした前進差分（境界外は0＝ノイマン条件）
    dx, dy = np.zeros((M, N)), np.zeros((M, N))
    dx[:, :-1] = wrap(psi[:, 1:] - psi[:, :-1])
    dy[:-1, :] = wrap(psi[1:, :] - psi[:-1, :])
    # 右辺 rho = div(W(grad psi))
    rho = np.zeros((M, N))
    rho[:, 0] = dx[:, 0]
    rho[:, 1:] = dx[:, 1:] - dx[:, :-1]
    rho[0, :] += dy[0, :]
    rho[1:, :] += dy[1:, :] - dy[:-1, :]
    # DCT領域で割り算
    R = dct2(rho)
    m, n = np.arange(M)[:, None], np.arange(N)[None, :]
    lam = 2 * np.cos(np.pi * m / M) + 2 * np.cos(np.pi * n / N) - 4
    lam[0, 0] = 1.0                      # ゼロ除算回避
    P = R / lam
    P[0, 0] = 0.0                        # 平均ゼロ解を選ぶ
    return idct2(P)

def raster_unwrap(psi):
    """比較用：1列目を縦にアンラップ→各行を横にアンラップ（ラスタ経路積分）"""
    u = np.unwrap(psi, axis=1)
    col = np.unwrap(psi[:, 0])
    return u + (col - psi[:, 0])[:, None]

def rms_err(est, truth):
    """全体定数の不定性を除いたRMS誤差"""
    d = est - truth
    return np.sqrt(np.mean((d - d.mean())**2))

print(f"無雑音: ラスタ {rms_err(raster_unwrap(psi2), phi2):.2e} rad, "
      f"最小二乗 {rms_err(ls_unwrap(psi2), phi2):.2e} rad")

無雑音ではどちらも $10^{-11}$ rad 以下、つまり数値誤差レベルで完全に一致します。レジデューがゼロなら経路積分は一意で厳密、そして最小二乗解も同じ答えに収束するからです。手法の差は「データが健全なとき」には現れません。差が出るのはレジデューが生まれてからです。

雑音を入れると、2つの手法の誤差の「形」がまったく違うことが見えてきます。左のラスタ経路積分は、レジデューから右へ伸びる濃紺の帯（$-2\pi$ や $-4\pi$ の段差）を作り、それ以外の場所はほぼ無誤差という局所的だが致命的な壊れ方をします。中央の最小二乗は段差が一切なく、代わりに薄い誤差が画像全体になだらかに広がる大域的だが穏やかな壊れ方です。右のヒストグラムでも、ラスタだけが $-2\pi \approx -6.3$ や $-4\pi \approx -12.6$ に離散的な山を作っているのが確認できます。

SNRを振ってRMS誤差を比較する

import numpy as np

rng = np.random.default_rng(0)
print(f"{'SNR[dB]':>8} {'レジデュー数':>12} {'DCT最小二乗':>14} {'ラスタ経路積分':>16}")
for snr_db in [30, 25, 20, 15, 12, 10, 8, 6, 4, 2, 0]:
    a = 10 ** (snr_db / 20)
    e_ls, e_ra, n_res = [], [], []
    for _ in range(12):                       # 12試行の平均
        z = np.exp(1j * phi2) + (rng.standard_normal(phi2.shape)
                                 + 1j * rng.standard_normal(phi2.shape)) / (np.sqrt(2) * a)
        p = np.angle(z)
        n_res.append(np.abs(residues(p)).sum())
        e_ls.append(rms_err(ls_unwrap(p), phi2))
        e_ra.append(rms_err(raster_unwrap(p), phi2))
    print(f"{snr_db:>8d} {np.mean(n_res):>12.0f} "
          f"{np.mean(e_ls):>14.3f} {np.mean(e_ra):>16.3f}")

得られる結果は次のとおりです（RMS誤差の単位は rad）。

この表からは3つのことが読み取れます。

第一に、レジデューがゼロの領域（25 dB以上）では両者の性能は完全に同一 です。誤差 0.022〜0.040 rad は雑音そのものによる位相揺らぎで、アルゴリズムの責任ではありません。レジデューが 0 個である限り、どんなアンラップ法も同じ答えを出します。

第二に、レジデューが1個でも現れた瞬間から差が開き始めます。20 dB でレジデュー1個のとき、既に最小二乗が 0.074 rad に対しラスタ経路積分は 0.190 rad と 2.6 倍です。たった1個のレジデューが、ラスタ走査では画像の一部に $2\pi$ の段差を作ってしまうためです。

第三に、最小二乗の誤差はレジデュー数に対してなだらかに増える のに対し、経路積分の誤差はより急激に増えます。10 dB では 0.575 rad 対 2.288 rad で4倍の開きです。最小二乗が「矛盾を全体に薄く均す」設計であることの効果が、そのまま数字に出ています。ただし 0 dB まで下がると 8.33 rad 対 13.05 rad と差が縮まります。誤差が $2\pi$ を超えるレベルになると、もはやどちらも意味のある答えを出していません。アンラップは魔法ではなく、データがItoh条件を大きく破っていれば手法選択では救えない ことを、この行が正直に示しています。

左のグラフでは、緑（最小二乗）と紫（ラスタ）が緑色に塗ったレジデュー0の領域では完全に重なり、25 dB を割った途端に離れ始めます。右のグラフは同じデータを「レジデュー数 対 誤差比」に描き替えたもので、レジデューが 7 個前後のところで比が 5.0 倍と最大になり、そこから増えるほど 1.0 に向かって落ちていきます。つまり 最小二乗の優位が最も大きいのはレジデューが少しだけあるときで、大量に出る劣悪な条件ではどちらも同じくらい役に立たなくなる、という非自明な構図が読み取れます。

Pythonでの実装4：重み付き最小二乗とPCG

最後に、局所的にコヒーレンスが崩れた場合を扱います。画像の一部だけSNRを 0 dB に落とし、それ以外は 25 dB という状況を作ります。InSARで湖や植生域が画像中にある状況の単純化です。

import numpy as np

def div_w(rx, ry):
    """重み付き残差場の発散"""
    M, N = rx.shape
    o = np.zeros((M, N))
    o[:, 0] += rx[:, 0]; o[:, 1:] += rx[:, 1:] - rx[:, :-1]
    o[0, :] += ry[0, :]; o[1:, :] += ry[1:, :] - ry[:-1, :]
    return o

def apply_A(f, wx, wy):
    """重み付きラプラシアン作用素 A f = -div(w grad f)"""
    M, N = f.shape
    gx, gy = np.zeros((M, N)), np.zeros((M, N))
    gx[:, :-1] = f[:, 1:] - f[:, :-1]
    gy[:-1, :] = f[1:, :] - f[:-1, :]
    return -div_w(wx * gx, wy * gy)

def wls_unwrap(psi, wx, wy, iters=500, tol=1e-9):
    """重み付き最小二乗アンラップ（DCT前処理付き共役勾配法）"""
    M, N = psi.shape
    dx, dy = np.zeros((M, N)), np.zeros((M, N))
    dx[:, :-1] = wrap(psi[:, 1:] - psi[:, :-1])
    dy[:-1, :] = wrap(psi[1:, :] - psi[:-1, :])
    b = -div_w(wx * dx, wy * dy)
    m, n = np.arange(M)[:, None], np.arange(N)[None, :]
    lam = 2 * np.cos(np.pi * m / M) + 2 * np.cos(np.pi * n / N) - 4
    lam[0, 0] = 1.0
    def precond(v):                      # 非重み付きポアソン解を前処理に使う
        P = -dct2(v) / lam; P[0, 0] = 0.0
        return idct2(P)
    f = np.zeros((M, N)); res = b - apply_A(f, wx, wy)
    z = precond(res); p = z.copy(); rz = np.sum(res * z)
    for it in range(iters):
        Ap = apply_A(p, wx, wy)
        alpha = rz / np.sum(p * Ap)
        f += alpha * p; res -= alpha * Ap
        if np.sqrt(np.mean(res**2)) < tol:
            break
        z = precond(res); rz_new = np.sum(res * z)
        p = z + (rz_new / rz) * p; rz = rz_new
    return f, it + 1

前処理関数 precond の中で符号を反転しているのは、作用素を $A = -\nabla^2$（正定値になる向き）と定義したためです。DCT領域の固有値 $\lambda$ はラプラシアンのものなので、$-1$ を掛けて $A$ の固有値に合わせます。あとは教科書どおりのPCG反復です。

import numpy as np

rng = np.random.default_rng(3)
M, N = phi2.shape
yy, xx = np.mgrid[0:M, 0:N]
blob = ((yy - 40)**2 + (xx - 90)**2) < 22**2      # 低コヒーレンス領域

snr_map = np.full((M, N), 10 ** (25 / 20.))
snr_map[blob] = 10 ** (0 / 20.)                   # ここだけ SNR 0 dB
z = np.exp(1j * phi2) + (rng.standard_normal((M, N))
                         + 1j * rng.standard_normal((M, N))) / (np.sqrt(2) * snr_map)
psi_b = np.angle(z)

coh = snr_map**2 / (1 + snr_map**2)               # 疑似コヒーレンス
w = coh ** 2                                      # 重み = コヒーレンスの2乗
wx, wy = np.zeros((M, N)), np.zeros((M, N))
wx[:, :-1] = np.minimum(w[:, 1:], w[:, :-1])      # 辺の重みは両端の小さい方
wy[:-1, :] = np.minimum(w[1:, :], w[:-1, :])

u_plain = ls_unwrap(psi_b)
u_weight, n_it = wls_unwrap(psi_b, wx, wy)
good = ~(((yy - 40)**2 + (xx - 90)**2) < 30**2)   # 健全域だけで評価

def rms_masked(est, truth, mask):
    d = est - truth
    return np.sqrt(np.mean((d - d[mask].mean())[mask]**2))

print(f"レジデュー総数: {np.abs(residues(psi_b)).sum()},  PCG反復回数: {n_it}")
print(f"健全域RMS誤差   非重み付き {rms_masked(u_plain, phi2, good):.3f} rad, "
      f"重み付き {rms_masked(u_weight, phi2, good):.3f} rad")

PCGはわずか11反復で収束します。DCT前処理が効いている証拠で、前処理なしのCGなら数百反復かかる規模です。誤差の改善も明確で、健全域のRMS誤差は非重み付きの 0.844 rad から重み付きの 0.365 rad へと 2.3 倍改善しました。健全域における最大誤差も 2.34 rad から 1.05 rad に下がります。

3枚目と4枚目を見比べると、重み付けが何をしているのかが一目で分かります。非重み付き（3枚目）では、緑の破線で囲んだ低コヒーレンス円の外側にまで赤と青の広い濃淡が染み出し、画像の左下から右上にかけて誤差の勾配ができています。重み付き（4枚目）では円の中の誤差はそのまま残る一方、外側の色がぐっと薄くなり、汚染が円の内側に閉じ込められています。復元できない領域を復元するのではなく、その害を他所に広げないのが重み付けの役割だということです。

読み取るべき本質は、重み付けは低コヒーレンス領域の位相を良くするのではなく、その汚染が健全域へ漏れ出すのを止める という点です。低コヒーレンス領域そのものは、どんな手法でも復元不能です（情報が無いのだから当然です）。しかし非重み付き最小二乗は、その領域の矛盾を「全体に薄く均す」設計ゆえに、まったく関係ない遠くの画素まで巻き込んで歪ませてしまいます。重みを入れることで、$L^2$ 最小化の柔軟さを保ったまま被害を局所化できるわけです。これが、実務のInSAR処理でコヒーレンス重み付けが標準となっている理由です。

対数目盛の縦軸に対して残差がほぼ直線的に落ちており、1反復あたり約1桁というきれいな線形収束を示しています。11反復で $10^{-9}$ に到達するということは、DCT2回ぶんのコストを22回払うだけで済むということで、非重み付き解法（DCT2回）に対しておよそ20倍のコストです。重み付けという実用上不可欠な拡張が、この程度の追加コストで手に入る点が、GhigliaとRomeroの解法が広く使われている理由です。

実務で気をつけること

理論と実験を踏まえて、位相アンラップを実務で使うときの指針を整理しておきます。

1. アンラップの前に、まずアンラップしやすくする。 アルゴリズムを凝るより、入力の質を上げるほうが圧倒的に効果的です。複素領域での多視化（マルチルック平均）やコヒーレンス適応フィルタ（Goldsteinフィルタなど）でSNRを上げれば、レジデュー数は激減します。SNRが 6 dB から 15 dB に上がるだけで、実験ではレジデューが 122 個から 7 個へ、RMS誤差が 1.41 rad から 0.18 rad へ改善しました。位相をラップした実数値のまま平均してはいけません。必ず複素数のまま平均して偏角を取ります。$-3.1$ rad と $+3.1$ rad の平均は $0$ ではなく $\pi$ 付近であるべきだからです。

2. 既知の縞成分は先に引く。 InSARなら地球楕円体成分と軌道縞、外部DEMによる地形縞を先に差し引く（フラットニングする）ことで、残差位相の勾配が小さくなり、Itoh条件が満たされやすくなります。FMCWレーダーなら粗距離から予測される位相を引きます。アンラップは「残差が小さい」ほど簡単になる ので、前段で減らせるものはすべて減らしておくのが鉄則です。

3. レジデューを必ず数える。 レジデュー数はアンラップの難易度を表す最良の単一指標です。処理の前後で数えておけば、結果を信じてよいかの判断材料になります。ゼロならどの手法でも同じ答えが出るので、最速の手法を選べばよいことになります。

4. 手法は「欠測を許すか」で選ぶ。 枝切り法や品質誘導法（quality-guided path following）は、健全な領域で厳密な答えを出す代わりに、アンラップ不能な領域を残します。最小二乗系は必ず全画素に値を出す代わりに、誤差を薄く広げます。信頼できない値を出すくらいなら欠測にしたい 用途（変動量の計測）と、穴のない面が必要な 用途（DEM生成）とでは、正しい選択が逆になります。

5. 結果の検証を忘れない。 アンラップ後の位相を再びラップして、元の観測位相と一致するかを確認します（$\mathcal{W}(\hat{\phi}) = \psi$ が成り立つべきです）。最小二乗系はこの一致性を保証しないので、不一致が大きい画素は誤差の集中を示す有力な手がかりになります。また複数の干渉ペアがあるなら、閉ループ位相（三角関係の一貫性）で検証する方法も強力です。

まとめ

本記事では、位相アンラップの理論とアルゴリズムを、1次元から2次元まで導出を追いながら解説しました。

• 折り返しは本質的: $e^{j\phi} = e^{j(\phi+2\pi k)}$ という複素指数の周期性により、観測位相は $\psi = \phi + 2\pi k$ の形でしか得られない。ラップ作用素 $\mathcal{W}(x) = \bmod(x+\pi, 2\pi) – \pi$ は $2\pi$ 周期性と、$|x| < \pi$ での恒等作用という2つの性質を持ち、これがアンラップの原理をすべて支えている。
• アンラップは不良設定問題: データだけからは解が定まらず、「隣接位相差が $\pi$ 未満」というなめらかさ仮定が必ず必要になる。
• 1次元は完全解: Itohの定理により、$|\phi_{n+1}-\phi_n| < \pi$ なら差分をラップして積分するだけで真の位相が厳密に復元される。この条件は $|f| < f_s/2$ というナイキスト条件そのもので、破れると周波数が $f - f_s$ に折り返される。実験では $f_s = 1000$ Hz で 480 Hz は正しく、520 Hz は $-480$ Hz と推定された。
• 周波数軸にも同じ条件: スペクトル位相から群遅延を求めるときも $\tau_g \Delta f < 1/2$ が必要で、ゼロ詰めによる周波数分解能の向上がRMS誤差を 0.707 秒から $4\times10^{-5}$ 秒へ激変させた。
• 2次元ではレジデューが現れる: $2\times2$ ループの巡回積分が $\pm 2\pi$ になる点がレジデューで、これはラップ勾配場の回転成分に対応する。任意の閉曲線の巡回積分は囲むレジデューの総電荷に等しい（離散Stokesの定理）。数値実験でも「経路差 $/2\pi$ ＝ 囲む総電荷」が厳密に成立し、SNR 6 dB では全画素の約15%で答えが経路に依存した。レジデューは正負が必ず対で生じる。
• 枝切り法: レジデューを反対電荷同士で結んで積分経路を遮断し、健全な領域では厳密な解を得る。誤差を局在化できる代わりに、アンラップ不能な領域が残る。
• 最小二乗アンラップ: $\min\sum|\nabla\phi – \mathcal{W}(\nabla\psi)|^2$ の正規方程式は離散ポアソン方程式 $\nabla^2\phi = \operatorname{div}\mathcal{W}(\nabla\psi)$ になり、ノイマン境界条件のもとDCTで対角化できる。固有値は $\lambda_{m,n} = 2\cos(\pi m/M) + 2\cos(\pi n/N) – 4$ で、$O(MN\log MN)$ で解ける。レジデュー（回転成分）を捨てるので必ず解が得られるが、誤差は全体に広がる。
• 重み付き最小二乗: コヒーレンスに応じた重みを入れると係数可変のポアソン方程式になり、DCT解を前処理としたPCGで11反復程度で解ける。低コヒーレンス領域の汚染が健全域に漏れるのを防ぎ、健全域のRMS誤差を 0.844 rad から 0.365 rad へ改善した。
• 前処理が最も効く: 多視化やコヒーレンスフィルタでSNRを上げ、既知の縞成分を先に引くことで、そもそもレジデューを作らないのが最良の戦略。

位相アンラップは、複素解析の多価性、標本化定理、ポアソン方程式、離散微分幾何（Stokesの定理）、最適化と、驚くほど多くの分野が交差する場所です。「たかが $2\pi$ の整数倍」と侮れない奥行きがあり、だからこそSARやレーダー、光計測の現場で今なお改良が続いています。

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

• 解析信号と瞬時周波数
• ヒルベルト変換
• 干渉SAR（InSAR）の原理
• 合成開口レーダー（SAR）の原理
• FMCWレーダーの原理
• 群遅延とは
• コヒーレンス関数とスペクトル解析
• ラプラシアン