SARレンジドップラーアルゴリズムを理解して実装する

人工衛星や航空機に載せたレーダで地上を撮影すると、雲や夜間に関係なく、数メートルの分解能で地表の構造を描き出せます。けれども、ここで一つ素朴な疑問が湧きます。衛星に載るアンテナはせいぜい数メートルしかありません。電波の回折を考えると、波長 5 cm(Cバンド)のレーダで 5 m のアンテナを使った場合、700 km 上空からの方位(進行方向)分解能は単純計算で 7 km 以上にもなってしまいます。これではビルどころか都市の区画すら見分けられません。それなのに、なぜ実際の衛星 SAR は方位方向でも数メートルの分解能を達成できるのでしょうか。

その答えが「合成開口(synthetic aperture)」という発想であり、その生データから鮮明な画像を再生する標準的な信号処理が「レンジドップラーアルゴリズム(Range-Doppler Algorithm, RDA)」です。衛星が飛びながら同じ目標を何度も照射し、受信した一連のエコーをコヒーレントに合成することで、あたかも数百メートルから数キロメートルの巨大なアンテナを使ったかのような、極めて高い方位分解能を作り出します。

このアルゴリズムを理解すると、次のような応用への見通しが一気に開けます。

  • 地球観測・防災: ALOS-2 や Sentinel-1 などの SAR 衛星は、洪水や地震による地表変動を全天候で監視します。生データから画像を作る心臓部がまさに RDA です
  • 干渉 SAR(InSAR): 2 回の観測の位相差から地盤沈下を mm 単位で計測する技術も、まず各観測を RDA で集束させることから始まります
  • 自動運転・ドローン: ミリ波レーダによる車載 SAR や、合成開口を使った高分解能イメージングレーダにも同じ原理が応用されています

本記事の内容

  • 合成開口がなぜ方位分解能を劇的に向上させるのか、その直感と導出
  • 送信チャープのレンジ圧縮(整合フィルタ)の原理
  • プラットフォーム移動によるアジマス方向のドップラー変調の導出
  • 合成開口長から方位分解能 $D/2$ を導く
  • レンジセルマイグレーション(RCM)の発生メカニズムと補正の必要性
  • レンジドップラー領域での RCM 補正(RCMC)
  • Python 実装: 点目標の生エコー生成 → レンジ圧縮 → RCMC → アジマス圧縮 → 集束画像の可視化

前提知識

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

SARの画像化とは — 2次元の集束問題

SAR の画像化を一言でいえば、「ぼやけた生データを、レンジ方向とアジマス方向の 2 方向それぞれで集束(focusing)させる」処理です。カメラで例えると、ピントの合っていない写真を後処理でシャープにするようなものですが、SAR の場合は最初から「ピントが合った像」が記録されているわけではありません。記録されているのは、目標からの反射波が時間的・空間的に広がった「生エコー(raw echo)」です。

ここで 2 つの方向を区別しておきましょう。レーダが電波を発射する向き、すなわちアンテナから地表へ斜めに伸びる方向を レンジ方向(range, スラントレンジ方向) と呼びます。一方、プラットフォーム(衛星や航空機)が飛行する向きを アジマス方向(azimuth, 方位方向) と呼びます。生データはこの 2 方向に広がった 2 次元の信号で、レンジ方向には送信パルスの長さの分だけ、アジマス方向にはアンテナのビーム幅で目標を照射し続ける時間の分だけ、エネルギーが拡散しています。

レンジドップラーアルゴリズムは、この 2 方向の拡散を別々の整合フィルタで畳み込み解除することで、各目標を本来あるべき 1 点に集束させます。レンジ方向の集束を「レンジ圧縮」、アジマス方向の集束を「アジマス圧縮」と呼びます。ただし、両者は完全に独立ではありません。衛星が飛行する間に目標までの距離が変化するため、1 つの目標のエコーが複数のレンジセルにまたがって記録される「レンジセルマイグレーション」という厄介な現象が起き、これを補正しないとアジマス圧縮がうまくいきません。

つまり RDA は次の 3 ステップから成ります。まずレンジ圧縮で距離方向を絞り込み、次にレンジセルマイグレーション補正(RCMC)でエコーの軌跡をまっすぐに直し、最後にアジマス圧縮で方位方向を絞り込みます。以降では、この 3 ステップを一つずつ、数式の導出を省略せずに丁寧に追っていきます。まずは出発点となる「レンジ圧縮」から見ていきましょう。

レンジ圧縮 — チャープと整合フィルタ

なぜチャープを使うのか

レーダの距離分解能は、送信パルスが短いほど良くなります。理想的には一瞬だけ光るような鋭いパルスを出せば、反射の戻り時間から距離をピンポイントで測れます。しかし、パルスを短くするとそのパルスに乗せられる総エネルギーが小さくなり、遠方の弱い目標を検出できなくなってしまいます。距離分解能(短いパルス)と検出能力(大きいエネルギー)は、素朴には両立しないトレードオフの関係にあるのです。

この矛盾を解くのが チャープ(chirp, 線形周波数変調) です。長い時間にわたって周波数を直線的に掃引する信号を送れば、パルス長を長く保ってエネルギーを稼ぎつつ、受信側の信号処理(パルス圧縮)で実効的に短いパルスと同等の分解能を得られます。広い帯域幅を持つ長いパルスを、後処理で鋭いピークに「圧縮」するわけです。

送信チャープ信号は、時間 $\tau$(高速時間, fast time と呼びます)に対して次のように書けます。

$$ s_{\text{tx}}(\tau) = \text{rect}\left(\frac{\tau}{T_p}\right) \exp\left(j\pi K_r \tau^2\right) $$

ここで $T_p$ はパルス幅、$K_r$ はチャープレート(周波数の掃引速度 [Hz/s])、$\text{rect}(\cdot)$ はパルス長で 1、それ以外で 0 をとる矩形窓です。瞬時周波数は位相を時間微分した $f(\tau) = K_r \tau$ であり、$\tau$ に比例して直線的に変化することが、まさに「線形周波数変調」の名の由来です。掃引する帯域幅は $B = |K_r| T_p$ となります。

整合フィルタによる圧縮

目標から戻ってきた受信信号は、送信チャープが時間遅延 $\tau_0 = 2R/c$($R$ は目標までの距離、$c$ は光速)だけ遅れたものです。

$$ s_{\text{rx}}(\tau) = A \, \text{rect}\left(\frac{\tau – \tau_0}{T_p}\right) \exp\left(j\pi K_r (\tau – \tau_0)^2\right) $$

この信号をできるだけ鋭いピークに変えるには、雑音中で信号対雑音比(SNR)を最大化するフィルタ、すなわち 整合フィルタ(matched filter) を使います。整合フィルタのインパルス応答は、送信信号の複素共役を時間反転したものです。チャープの場合、

$$ h(\tau) = s_{\text{tx}}^*(-\tau) = \text{rect}\left(\frac{\tau}{T_p}\right)\exp\left(-j\pi K_r \tau^2\right) $$

レンジ圧縮の出力は、受信信号と整合フィルタの畳み込みで与えられます。

$$ s_{\text{rc}}(\tau) = s_{\text{rx}}(\tau) * h(\tau) $$

この畳み込みを実際に計算すると何が起こるかを見てみましょう。畳み込みは積分で定義されますが、チャープの圧縮特性を理解するには周波数領域で考えるのが見通しが良いです。チャープのスペクトルは(停留位相法による近似で)帯域 $B$ にわたってほぼ平坦な振幅と、$K_r$ に対応する 2 次位相を持ちます。整合フィルタはこの 2 次位相を打ち消す(共役なので符号が逆)ため、出力スペクトルは平坦な振幅だけが残り、その逆フーリエ変換は理論上 sinc 関数になります。

$$ s_{\text{rc}}(\tau) \approx A \, T_p \, \text{sinc}\left[B(\tau – \tau_0)\right] \exp(-j 2\pi f_c \tau_0) $$

ここで $f_c$ は搬送波周波数です。重要なのは、出力が遅延 $\tau_0$ を中心とする鋭い sinc ピークになっている点です。長かった入力パルスが、幅 $\approx 1/B$ の鋭いピークに圧縮されたのです。sinc 関数の主ローブ幅から、距離分解能(スラントレンジ分解能)は次のように決まります。

$$ \rho_r = \frac{c}{2B} $$

たとえば帯域幅 $B = 30$ MHz なら $\rho_r = 5$ m となります。パルス長 $T_p$ がいくら長くても、分解能を決めるのは帯域幅 $B$ だけだという点が、チャープとパルス圧縮の核心です。エネルギーは $T_p$ で稼ぎ、分解能は $B$ で稼ぐ。両者を独立に設計できることが、この方式の最大の利点です。

圧縮利得とSNR

圧縮の御利益はもう一つあります。長さ $T_p$、帯域 $B$ のチャープを圧縮すると、SNR が時間帯域幅積 $T_p B$ 倍だけ向上します。これを パルス圧縮利得(pulse compression gain) と呼びます。たとえば $T_p = 40\,\mu$s、$B = 30$ MHz なら $T_p B = 1200$ で、約 31 dB もの SNR 改善が得られます。短い物理パルスでは到底届かない遠方の目標も、長いチャープと圧縮処理の組み合わせで検出できるのです。

ここまでで、レンジ方向に拡散したエコーを 1 本の鋭いピークに集束させる方法がわかりました。しかし SAR ではプラットフォームが飛行しているため、もう一方のアジマス方向にもエコーが拡散しています。次は、この拡散がどのように生じるのか — すなわちアジマス方向のドップラー変調を導出しましょう。

アジマス方向のドップラー変調

なぜアジマス方向に信号が広がるのか

レーダのアンテナは、ある幅を持ったビームで地表を照らしています。プラットフォームが飛行すると、ある地上の点目標は、ビームに入ってから出ていくまでの間、繰り返し照射され続けます。この間に送受信した一連のパルスを「アジマス方向の信号」と見なします。アジマス方向の時間(低速時間, slow time と呼び $\eta$ で表します)は、パルスを 1 つ送るごとに 1 ステップ進む、いわば「飛行に沿った時間軸」です。

ここで鍵になるのが、プラットフォームが目標に近づき、最接近し、遠ざかっていく過程で、目標までの距離が刻一刻と変化することです。距離が変われば往復の電波の位相が変わり、これがアジマス方向の信号に変調をかけます。近づくときは波が圧縮されて周波数が上がり、遠ざかるときは伸びて周波数が下がる — これはまさにドップラー効果です。

距離履歴の導出

プラットフォームが速度 $v$ で直線飛行し、低速時間 $\eta$ における位置を考えます。$\eta = 0$ で目標に最接近し、そのときのスラントレンジ(最短距離)を $R_0$ とします。プラットフォームが $\eta$ 秒間に進む距離は $v\eta$ なので、目標までの瞬時距離 $R(\eta)$ はピタゴラスの定理から、

$$ R(\eta) = \sqrt{R_0^2 + (v\eta)^2} $$

となります。SAR では通常 $v\eta \ll R_0$(飛行距離が最短距離よりずっと小さい)が成り立つため、平方根をテイラー展開して近似します。$\sqrt{R_0^2 + x} \approx R_0 + \frac{x}{2R_0}$($x = (v\eta)^2$)を使うと、

$$ R(\eta) \approx R_0 + \frac{(v\eta)^2}{2R_0} $$

距離が低速時間 $\eta$ の 2 次関数(放物線)で変化することがわかります。最接近点 $\eta = 0$ で最小値 $R_0$ をとり、前後に離れるにつれて 2 次的に増える、という直感に合った形です。

アジマス信号の位相

受信信号の位相は、往復距離 $2R(\eta)$ に対応します。搬送波の波長を $\lambda$ とすると、アジマス方向の信号位相 $\phi(\eta)$ は、

$$ \phi(\eta) = -\frac{4\pi}{\lambda} R(\eta) $$

です(往復で $2R$、位相は $2\pi/\lambda \times 2R = 4\pi R / \lambda$、近づくほど位相が進むので符号は負)。先ほどの距離近似を代入すると、

$$ \phi(\eta) \approx -\frac{4\pi}{\lambda}\left(R_0 + \frac{v^2 \eta^2}{2R_0}\right) = -\frac{4\pi R_0}{\lambda} – \frac{2\pi v^2}{\lambda R_0}\eta^2 $$

第 1 項は $\eta$ によらない定数位相なので、信号の形には影響しません。本質は第 2 項の $\eta^2$ に比例する 2 次位相 です。瞬時ドップラー周波数は位相を $2\pi$ で割って時間微分したものなので、

$$ f_\eta(\eta) = \frac{1}{2\pi}\frac{d\phi}{d\eta} = -\frac{2v^2}{\lambda R_0}\eta $$

ドップラー周波数が低速時間 $\eta$ に比例して直線的に変化しています。これは、レンジ方向のチャープと全く同じ「線形周波数変調(チャープ)」の構造です。つまり、アジマス方向の信号もチャープであり、その アジマスチャープレート(ドップラーレート) $K_a$ は、

$$ K_a = \frac{2v^2}{\lambda R_0} $$

と定義できます(符号は近づく→遠ざかるで負、ここでは絶対値の大きさに注目します)。ドップラー周波数が最接近時 $\eta = 0$ でゼロを横切り、近づくときは正、遠ざかるときは負になる様子は、列車が近づくと汽笛が高く、通り過ぎると低くなる日常のドップラー効果そのものです。

アジマス圧縮の発想

アジマス方向の信号がチャープであるということは、レンジ圧縮と全く同じ整合フィルタの手法でアジマス方向も圧縮できるということです。アジマス方向の整合フィルタは、ドップラーレート $K_a$ を打ち消す 2 次位相を持つフィルタ、

$$ h_a(\eta) = \exp\left(j\pi K_a \eta^2\right) $$

であり、これをアジマス方向に畳み込めば、拡散していた目標のエネルギーが 1 点に集束します。レンジ圧縮とアジマス圧縮が同じ「チャープの整合フィルタリング」という枠組みで統一的に扱えること — これが RDA の美しさです。次に、このアジマス圧縮が達成できる分解能を計算し、合成開口の威力を定量的に確かめましょう。

合成開口長と方位分解能

合成開口長

冒頭の疑問に戻ります。なぜ小さなアンテナで高い方位分解能が得られるのでしょうか。鍵は、目標が照射され続ける間にプラットフォームが移動した距離 — すなわち 合成開口長(synthetic aperture length) $L_s$ にあります。物理アンテナの長さを $D$ とすると、その方位方向のビーム幅(半値全幅の近似)は回折により、

$$ \theta_a \approx \frac{\lambda}{D} $$

です。距離 $R_0$ にある目標は、このビーム幅で照射される間ずっとレーダの視野に入っています。地表でビームが覆う方位方向の長さ、すなわち目標が照射され続ける区間の長さが合成開口長で、

$$ L_s = R_0 \theta_a = \frac{\lambda R_0}{D} $$

となります。プラットフォームはこの $L_s$ の区間を飛びながら同じ目標を観測し続け、その全エコーをコヒーレントに合成します。あたかも長さ $L_s$ の巨大なアンテナを作ったかのように振る舞うため「合成開口」と呼ぶわけです。たとえば $\lambda = 0.056$ m(Cバンド)、$R_0 = 700$ km、$D = 5$ m なら $L_s \approx 7.8$ km にもなります。

方位分解能の導出

合成開口(長さ $L_s$ の実効アンテナ)の方位方向ビーム幅は、開口が長いほど狭くなります。一般にアンテナのビーム幅は開口長に反比例しますが、合成開口の場合は 片道ではなく往復で位相が効く ため、通常の実開口アンテナの半分のビーム幅になります。実効的なビーム幅は、

$$ \theta_s = \frac{\lambda}{2 L_s} $$

です(係数 $1/2$ が往復経路の効果)。このビーム幅が距離 $R_0$ で張る地表の幅が、達成可能な方位分解能 $\rho_a$ です。

$$ \rho_a = R_0 \theta_s = R_0 \cdot \frac{\lambda}{2 L_s} $$

ここに合成開口長 $L_s = \lambda R_0 / D$ を代入します。すると、

$$ \rho_a = R_0 \cdot \frac{\lambda}{2} \cdot \frac{D}{\lambda R_0} = \frac{D}{2} $$

という、驚くほど簡潔な結果が得られます。方位分解能は物理アンテナ長 $D$ の半分 であり、波長にも距離にも飛行高度にも依存しません。これは SAR の最も衝撃的な結論です。アンテナを小さくする($D$ を小さくする)ほど分解能が良くなるという、常識に反した関係になっています。

なぜアンテナが小さいほど良いのか

直感的に説明しましょう。アンテナを小さくすると、ビーム幅 $\lambda/D$ が広がり、1 つの目標を照射し続ける時間が長くなります。すると合成開口長 $L_s$ が伸びるため、実効的により大きなアンテナを合成でき、結果として分解能が向上します。「物理アンテナが小さい → ビームが広い → 長く照射できる → 合成開口が長い → 分解能が良い」という連鎖です。

もちろん現実には、小さすぎるアンテナは利得が低くて SNR が確保できない、サンプリングレート(パルス繰り返し周波数, PRF)の制約がある、などの理由で限界があります。それでも、$\rho_a = D/2$ という関係式は、SAR がなぜ革命的だったかを端的に物語っています。たとえば $D = 10$ m のアンテナで、わずか 5 m の方位分解能が原理的に達成できるのです。

ここまでで、レンジ方向とアジマス方向それぞれの集束原理と、達成できる分解能がわかりました。しかし、レンジとアジマスを単純に独立処理するだけでは画像はボケてしまいます。その元凶が、次に説明する「レンジセルマイグレーション」です。

レンジセルマイグレーションと補正

マイグレーションの発生

アジマス方向の距離履歴 $R(\eta) = R_0 + v^2\eta^2/(2R_0)$ をもう一度見てください。これは「目標までの距離が低速時間とともに変化する」という意味です。レンジ圧縮された生データを 2 次元の表(縦軸=アジマス、横軸=レンジ)として見ると、1 つの点目標のエコーは、$R(\eta)$ に従って 放物線状に複数のレンジセルを横断 して記録されます。

これが レンジセルマイグレーション(Range Cell Migration, RCM) です。目標は地上では 1 点なのに、生データ上ではアジマス方向に湾曲した軌跡として広がってしまうのです。マイグレーション量(最接近点から開口の端までの距離変化)は、

$$ \Delta R = R(\eta_{\max}) – R_0 = \frac{v^2 \eta_{\max}^2}{2R_0} \approx \frac{L_s^2}{8 R_0} $$

程度になります($v\eta_{\max} \approx L_s/2$ を使用)。このマイグレーション量がレンジ分解能 $\rho_r$ の数分の 1 を超えると、アジマス圧縮で問題になります。なぜなら、アジマス圧縮はある固定したレンジセルに沿って整合フィルタをかける処理であり、目標のエネルギーが複数のレンジセルに散らばっていると、1 つのフィルタでは全エネルギーを集められず、ピークがボケて分解能と SNR が劣化するからです。

レンジドップラー領域での補正

ここで RDA の名前の由来である巧妙なアイデアが登場します。マイグレーションの軌跡 $R(\eta)$ は、同じレンジ $R_0$ にある目標であればアジマス位置によらず同じ形をしています。けれども異なる目標ごとに $\eta$ の中心がずれているため、時間領域(アジマス時間)のままでは目標ごとに異なる補正をしなければならず厄介です。

そこで、データをアジマス方向にフーリエ変換して レンジドップラー領域(range-Doppler domain) に移します。アジマス時間 $\eta$ をアジマス周波数(ドップラー周波数)$f_\eta$ に変換するのです。先に求めた瞬時ドップラー周波数の関係 $f_\eta = -K_a \eta$($K_a = 2v^2/(\lambda R_0)$)を使うと、$\eta = -f_\eta / K_a$ と書けます。これを距離履歴に代入すると、

$$ R(f_\eta) = R_0 + \frac{v^2 \eta^2}{2R_0} = R_0 + \frac{v^2}{2R_0}\cdot\frac{f_\eta^2}{K_a^2} $$

$K_a = 2v^2/(\lambda R_0)$ を代入して整理します。分母の $K_a^2 = 4v^4/(\lambda^2 R_0^2)$ を使うと、

$$ \Delta R(f_\eta) = \frac{v^2}{2R_0}\cdot\frac{\lambda^2 R_0^2}{4 v^4} f_\eta^2 = \frac{\lambda^2 R_0}{8 v^2} f_\eta^2 $$

という、アジマス周波数 $f_\eta$ だけで決まるマイグレーション量が得られます。ここがポイントです。レンジドップラー領域では、同じレンジ $R_0$ にあるすべての目標が、アジマス位置によらず全く同じマイグレーション軌跡を共有 します。なぜなら、目標のアジマス位置の違いはフーリエ変換で線形位相(時間シフト)に変わるだけで、振幅としての軌跡の形は共通になるからです。したがって、各レンジ・各ドップラー周波数に対して 1 つの補正量を計算すれば、その行のすべての目標を一括で補正できます。

RCMC の実装方法

補正(Range Cell Migration Correction, RCMC)は、レンジドップラー領域で各アジマス周波数ビンごとに、レンジ方向のデータを $\Delta R(f_\eta)$ に対応するサンプル数だけシフト(補間)することで行います。シフト量は一般に整数サンプルにならないため、sinc 補間や線形補間でサブサンプル単位の移動を実現します。この処理により、放物線状だったエコーの軌跡が、レンジドップラー領域でまっすぐな一直線(一定レンジ)に整列します。

軌跡がまっすぐになれば、各レンジセルに沿ってアジマス整合フィルタをかけるだけで、目標のエネルギーを 1 点に集められます。RCMC はレンジ圧縮とアジマス圧縮を橋渡しする不可欠なステップなのです。なお、補正量はレンジ $R_0$(=レンジビン)にも依存するため、厳密にはレンジごとに異なる補正をかけます。本記事の実装では、点目標の近傍に絞り、各レンジビンごとに $\Delta R(f_\eta)$ を計算して補正します。

これで RDA の 3 ステップ — レンジ圧縮、RCMC、アジマス圧縮 — の理論がすべて揃いました。次節では、これらをアルゴリズムとして整理し、Python で点目標の生エコーから集束画像を再生してみましょう。

レンジドップラーアルゴリズムの全体像

これまでの導出を、処理の流れとして整理します。入力は 2 次元の生エコー(縦:アジマスサンプル、横:レンジサンプル)、出力は集束された複素画像です。

  1. レンジ圧縮 — 各アジマス行に対し、レンジ方向のチャープ整合フィルタを畳み込む。周波数領域で高速に行うため、レンジ FFT → 参照スペクトル乗算 → レンジ IFFT の順で実行する。これで各目標はレンジ方向に集束し、放物線状の軌跡として現れる
  2. アジマス FFT — データをアジマス方向にフーリエ変換し、レンジドップラー領域へ移す。この領域では同一レンジの目標が共通のマイグレーション軌跡を持つ
  3. RCMC — 各アジマス周波数ビンごとに $\Delta R(f_\eta) = \lambda^2 R_0 f_\eta^2 / (8v^2)$ を計算し、レンジ方向に補間シフトして軌跡をまっすぐにする
  4. アジマス圧縮 — レンジドップラー領域で、各レンジビンのドップラーレート $K_a$ に整合した参照関数 $\exp(j\pi f_\eta^2 / K_a)$ を乗算する(周波数領域での整合フィルタリング)
  5. アジマス IFFT — アジマス方向に逆フーリエ変換し、時間領域(画像領域)へ戻す。これで各目標が 1 点に集束した複素画像が得られる

ステップ 1 と 4 がそれぞれレンジ・アジマスの整合フィルタ、ステップ 3 が両者をつなぐ補正です。すべて FFT を使って周波数領域で行うため、計算量は畳み込みを直接行うより圧倒的に少なくて済みます。それでは、この流れをコードに落とし込んでいきましょう。

Python実装:生エコーの生成

まず、シミュレーションのパラメータを設定し、点目標からの生エコーを合成します。実データの代わりに、距離履歴とチャープから理論的に生エコーを組み立てることで、アルゴリズムが正しく点目標を 1 点に集束させるかを検証できます。

import numpy as np
import matplotlib.pyplot as plt

# ---- レーダ・幾何パラメータ ----
c = 3e8                 # 光速 [m/s]
fc = 5.3e9              # 搬送波周波数 [Hz] (Cバンド)
lam = c / fc           # 波長 [m]
Tp = 2.5e-6            # パルス幅 [s]
B = 30e6              # チャープ帯域幅 [Hz]
Kr = B / Tp          # レンジチャープレート [Hz/s]
fs = 45e6            # レンジ方向サンプリング周波数 [Hz]

v = 150.0            # プラットフォーム速度 [m/s] (航空機SAR想定)
R0 = 15000.0         # 点目標の最短スラントレンジ [m]
PRF = 200.0          # パルス繰返し周波数 [Hz] (アジマスサンプリング)

print(f"波長 lambda      = {lam*100:.2f} cm")
print(f"レンジ分解能 rho_r = {c/(2*B):.2f} m")
print(f"アジマスチャープレート Ka = {2*v**2/(lam*R0):.2f} Hz/s")

このコードはシミュレーションの基本パラメータを定義します。Cバンド(5.3 GHz)、帯域 30 MHz でレンジ分解能は約 5 m になるはずです。出力されるアジマスチャープレート $K_a = 2v^2/(\lambda R_0)$ は、後でアジマス圧縮の参照関数に使う重要な値です。航空機 SAR を想定した穏やかなパラメータにすることで、計算量を抑えつつ本質を確認できます。

次に、観測区間(合成開口)の長さと、データ行列のサイズを決めます。アンテナのビーム幅から合成開口長を計算し、それに対応するアジマスサンプル数を確保します。

# ---- 観測ジオメトリと配列サイズ ----
D = 2.0                          # アンテナ長 [m]
theta_bw = lam / D               # アジマスビーム幅 [rad]
Ls = R0 * theta_bw               # 合成開口長 [m]
Ta = Ls / v                      # 目標照射時間 [s]
print(f"合成開口長 Ls = {Ls:.1f} m, 照射時間 Ta = {Ta:.3f} s")
print(f"方位分解能(理論) rho_a = D/2 = {D/2:.2f} m")

# アジマス(低速時間)軸: 照射時間の少し広めを確保
Na = int(np.ceil(1.4 * Ta * PRF))      # アジマスサンプル数
Na += Na % 2                            # 偶数化
eta = (np.arange(Na) - Na//2) / PRF     # 低速時間 [s], 中心が最接近

# レンジ(高速時間)軸: 目標距離の前後をカバー
Rmin, Rmax = R0 - 60, R0 + 60           # スラントレンジ窓 [m]
tau0 = 2*Rmin/c                          # 窓の開始遅延 [s]
Nr = int(np.ceil((2*(Rmax-Rmin)/c + Tp) * fs))
Nr += Nr % 2
tau = tau0 + np.arange(Nr) / fs          # 高速時間(往復遅延) [s]
print(f"データ行列サイズ: アジマス {Na} x レンジ {Nr}")

このコードは合成開口長 $L_s = \lambda R_0 / D$ と照射時間 $T_a = L_s/v$ を計算し、それに合わせて 2 次元データ行列の大きさを決めています。アジマス軸は最接近点を中心に対称にとり、レンジ軸は目標距離の前後 60 m をカバーします。出力される方位分解能 $D/2 = 1.0$ m は、後で再生画像の主ローブ幅と比較する目標値です。

いよいよ生エコーを合成します。各低速時間 $\eta$ における距離 $R(\eta)$ を計算し、その遅延位置に送信チャープを配置し、往復位相 $\exp(-j4\pi R/\lambda)$ を掛けます。

# ---- 点目標の生エコー生成 ----
raw = np.zeros((Na, Nr), dtype=complex)

for i, et in enumerate(eta):
    R_eta = np.sqrt(R0**2 + (v*et)**2)   # 瞬時スラントレンジ
    td = 2*R_eta/c                        # 往復遅延 [s]
    # レンジチャープ(エンベロープ): |tau-td| <= Tp/2 の範囲で有効
    dt = tau - td
    env = (np.abs(dt) <= Tp/2).astype(float)
    range_chirp = env * np.exp(1j*np.pi*Kr*dt**2)
    # アジマス方向の往復位相(ドップラー変調)
    az_phase = np.exp(-1j*4*np.pi*R_eta/lam)
    # ビーム重み(矩形ビーム近似): 照射区間内のみ
    beam = 1.0 if np.abs(et) <= Ta/2 else 0.0
    raw[i, :] = beam * az_phase * range_chirp

# 生エコーの振幅を可視化
plt.figure(figsize=(8, 5))
plt.imshow(np.abs(raw), aspect='auto', cmap='viridis',
           extent=[Rmin, Rmax, eta[-1], eta[0]])
plt.xlabel('Slant range [m]')
plt.ylabel('Azimuth slow-time [s]')
plt.title('Raw SAR echo (amplitude) of a point target')
plt.colorbar(label='|raw|')
plt.tight_layout()
plt.show()

生成された生エコーの振幅画像を見ると、点目標であるにもかかわらず、エネルギーがレンジ方向にもアジマス方向にも大きく広がっていることがわかります。レンジ方向の広がりはチャープのパルス長 $T_p$ による拡散、アジマス方向の広がりは照射時間 $T_a$ にわたって目標が照らされ続けることによる拡散です。さらに目を凝らすと、エコーの帯がアジマス方向に沿ってわずかに湾曲している — これがレンジセルマイグレーションの正体です。この「ぼやけた」生データを、これから 3 ステップで 1 点に集束させていきます。

Python実装:レンジ圧縮

最初のステップはレンジ圧縮です。各アジマス行に対し、レンジ方向のチャープ整合フィルタを周波数領域で適用します。レンジ FFT → 参照スペクトルの共役を乗算 → レンジ IFFT、という流れです。

# ---- レンジ圧縮(周波数領域での整合フィルタ) ----
# 参照チャープ(送信波形と同じ)を時間原点に作る
t_ref = (np.arange(Nr) - Nr//2) / fs
ref_chirp = (np.abs(t_ref) <= Tp/2) * np.exp(1j*np.pi*Kr*t_ref**2)

# 整合フィルタ = 参照チャープの複素共役のスペクトル
Href = np.conj(np.fft.fft(np.fft.ifftshift(ref_chirp)))

# 各アジマス行をレンジFFT → フィルタ乗算 → IFFT
Raw_f = np.fft.fft(raw, axis=1)
rc = np.fft.ifft(Raw_f * Href[None, :], axis=1)

# レンジ圧縮後の振幅
plt.figure(figsize=(8, 5))
plt.imshow(np.abs(rc), aspect='auto', cmap='viridis',
           extent=[Rmin, Rmax, eta[-1], eta[0]])
plt.xlabel('Slant range [m]')
plt.ylabel('Azimuth slow-time [s]')
plt.title('After range compression')
plt.colorbar(label='|rc|')
plt.tight_layout()
plt.show()

レンジ圧縮後の画像では、レンジ方向に広がっていたエネルギーが鋭い縦線(細い帯)に圧縮されています。これがパルス圧縮の効果です。一方、アジマス方向にはまだエネルギーが広がったままで、しかもその帯がアジマス方向に沿ってきれいな放物線を描いて湾曲していることがはっきり見えます。これがレンジセルマイグレーションで、最接近点(中央)で最もレンジが小さく、両端でレンジが大きくなる $R(\eta) \approx R_0 + v^2\eta^2/(2R_0)$ の形そのものです。この湾曲を放置するとアジマス圧縮が破綻するため、次に RCMC で補正します。

Python実装:アジマスFFTとRCMC

レンジドップラー領域に移るため、データをアジマス方向に FFT します。その上で、各ドップラー周波数 $f_\eta$ におけるマイグレーション量 $\Delta R(f_\eta) = \lambda^2 R_0 f_\eta^2 / (8v^2)$ を計算し、レンジ方向に補間シフトして軌跡をまっすぐにします。

# ---- アジマスFFTでレンジドップラー領域へ ----
rc_rd = np.fft.fft(rc, axis=0)            # アジマス方向FFT
f_eta = np.fft.fftfreq(Na, d=1/PRF)        # ドップラー周波数軸 [Hz]

# 各レンジビンの参照レンジ(窓内の実距離)
range_axis = Rmin + (tau - tau0) * c / 2   # 各レンジサンプルの距離 [m]

# ---- RCMC: ドップラー周波数ごとにレンジ方向シフト ----
rcmc = np.zeros_like(rc_rd)
for k in range(Na):
    # このドップラー周波数でのマイグレーション量 [m]
    dR = lam**2 * R0 * f_eta[k]**2 / (8 * v**2)
    # レンジサンプル単位のシフト量(往復なので距離→遅延→サンプル)
    shift_samp = dR * 2 / c * fs
    # sinc補間で非整数シフト(周波数領域での線形位相シフト)
    line = rc_rd[k, :]
    Nf = len(line)
    freq = np.fft.fftfreq(Nf)
    phase = np.exp(-1j*2*np.pi*freq*(-shift_samp))  # -shiftで近距離側へ詰める
    rcmc[k, :] = np.fft.ifft(np.fft.fft(line) * phase)

# RCMC後をアジマスIFFTして時間領域で確認
rc_check = np.fft.ifft(rcmc, axis=0)
plt.figure(figsize=(8, 5))
plt.imshow(np.abs(rc_check), aspect='auto', cmap='viridis',
           extent=[Rmin, Rmax, eta[-1], eta[0]])
plt.xlabel('Slant range [m]')
plt.ylabel('Azimuth slow-time [s]')
plt.title('After RCMC (range-time domain check)')
plt.colorbar(label='|rc|')
plt.tight_layout()
plt.show()

RCMC 後の画像を、レンジ圧縮直後の画像と見比べてください。先ほどまで放物線状に湾曲していたエコーの帯が、ほぼ垂直な一直線(一定レンジ)に整列しているはずです。これは、各ドップラー周波数に応じてレンジ方向のサンプルを補間シフトし、最接近点のレンジ $R_0$ に揃えた結果です。シフトには周波数領域での線形位相乗算(フーリエシフト定理)を使うことで、非整数サンプルの移動を sinc 補間と等価に実現しています。軌跡がまっすぐになったことで、ようやく各レンジセルに沿ったアジマス整合フィルタリングが正しく機能する準備が整いました。

Python実装:アジマス圧縮と集束画像

最後のステップはアジマス圧縮です。レンジドップラー領域のまま、各レンジビンのドップラーレート $K_a$ に整合した参照関数を乗算し、アジマス IFFT で画像領域に戻します。

# ---- アジマス圧縮(レンジドップラー領域での整合フィルタ) ----
# 各レンジビンごとのドップラーレート Ka = 2 v^2 / (lam * R)
Ka = 2 * v**2 / (lam * range_axis)         # 形状 (Nr,)

# アジマス整合フィルタ H_az(f_eta) = exp(+j*pi*f_eta^2 / Ka)
# (アジマスチャープの2次位相を打ち消す)
Haz = np.exp(1j*np.pi * f_eta[:, None]**2 / Ka[None, :])  # (Na, Nr)

ac_rd = rcmc * Haz                          # 整合フィルタ乗算
img = np.fft.ifft(ac_rd, axis=0)            # アジマスIFFTで画像へ
img = np.fft.fftshift(img, axes=0)          # アジマス中心を画面中央へ

# 集束画像(dBスケール)
img_db = 20*np.log10(np.abs(img)/np.abs(img).max() + 1e-6)
plt.figure(figsize=(8, 5))
plt.imshow(img_db, aspect='auto', cmap='jet', vmin=-40, vmax=0,
           extent=[Rmin, Rmax, -Na/2/PRF*v, Na/2/PRF*v])
plt.xlabel('Slant range [m]')
plt.ylabel('Azimuth [m]')
plt.title('Focused SAR image of point target [dB]')
plt.colorbar(label='Normalized power [dB]')
plt.tight_layout()
plt.show()

集束画像(dB スケール)を見ると、これまで 2 次元に大きく広がっていた点目標のエネルギーが、ついに 1 つの鋭い点に集束しています。周囲に十字状に伸びるのは sinc 関数のサイドローブで、レンジ・アジマス両方向の整合フィルタが理論通りに機能していることの証です。生エコーの「ぼやけた帯」から、はっきりとした 1 点へ — これがレンジドップラーアルゴリズムによる画像再生の全貌です。

最後に、集束した点目標のアジマス方向の断面(インパルス応答)を取り出し、達成された方位分解能が理論値 $D/2$ に一致するかを定量的に確認しましょう。

# ---- 方位分解能の検証(アジマス断面の半値全幅) ----
# ピーク位置を探す
pk = np.unravel_index(np.argmax(np.abs(img)), img.shape)
az_cut = np.abs(img[:, pk[1]])
az_cut /= az_cut.max()
az_axis = (np.arange(Na) - Na//2) / PRF * v   # アジマス座標 [m]

# -3dB(振幅0.707)幅を線形補間で推定
half = 1/np.sqrt(2)
above = np.where(az_cut >= half)[0]
res_meas = az_axis[above[-1]] - az_axis[above[0]]

plt.figure(figsize=(8, 5))
plt.plot(az_axis, 20*np.log10(az_cut + 1e-6), 'b-')
plt.axhline(-3, color='r', ls='--', label='-3 dB')
plt.xlim(-15, 15)
plt.ylim(-40, 1)
plt.xlabel('Azimuth [m]')
plt.ylabel('Normalized power [dB]')
plt.title(f'Azimuth impulse response (measured FWHM = {res_meas:.2f} m)')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

print(f"理論方位分解能 D/2 = {D/2:.2f} m")
print(f"実測 -3dB 幅       = {res_meas:.2f} m")

このグラフは集束された点目標のアジマス方向インパルス応答です。中央に鋭い主ローブがあり、その両側に対称なサイドローブが並ぶ、典型的な sinc 状の応答が得られています。-3 dB 幅(半値全幅)が理論値 $D/2 = 1.0$ m とおおむね一致していれば、アジマス圧縮が正しく機能し、合成開口によって物理アンテナ長の半分の分解能が達成されたことになります。実測値が理論値と多少ずれるのは、矩形ビーム近似や有限のドップラー帯域、補間誤差によるもので、窓関数を使えばサイドローブを下げられる一方で主ローブはやや広がります。これは実際の SAR 処理でも常に直面する分解能とサイドローブのトレードオフです。

まとめ

本記事では、合成開口レーダの画像化原理と、その標準的な処理であるレンジドップラーアルゴリズムを、数式の導出から Python 実装まで一気通貫で解説しました。要点を整理します。

  • レンジ圧縮: 長いチャープを送り、整合フィルタで鋭いピークに圧縮する。距離分解能は帯域幅で決まり $\rho_r = c/(2B)$、エネルギーはパルス長で稼ぐため両者を独立に設計できる
  • アジマスのドップラー変調: プラットフォーム移動で距離が $R(\eta) \approx R_0 + v^2\eta^2/(2R_0)$ と変化し、アジマス信号がドップラーレート $K_a = 2v^2/(\lambda R_0)$ のチャープになる。これも整合フィルタで圧縮できる
  • 方位分解能 $D/2$: 合成開口長 $L_s = \lambda R_0/D$ を合成することで、方位分解能は物理アンテナ長の半分になり、波長にも距離にも依存しない。アンテナが小さいほど分解能が良いという逆説が成り立つ
  • RCMC: 距離変化で生じるレンジセルマイグレーションを、レンジドップラー領域で一括補正する。これがレンジ圧縮とアジマス圧縮を橋渡しする
  • RDA 全体: レンジ圧縮 → アジマス FFT → RCMC → アジマス圧縮 → アジマス IFFT の 5 段階で、ぼやけた生エコーが鋭い点像に集束する

レンジドップラーアルゴリズムは、SAR 処理の中で最も古典的かつ直感的な手法です。これを足がかりに、より高精度な処理(チャープスケーリングアルゴリズム、$\omega$-K アルゴリズム)や、複数観測の位相差を使う干渉 SAR(InSAR)、偏波情報を使うポラリメトリ SAR へと発展していきます。

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