空間時間適応処理(STAP)を理解して実装する

地上の固定レーダーであれば、地面からの反射(クラッタ)はドップラー周波数がほぼゼロに集中するため、ゼロ近傍の周波数だけを除去すれば動目標を浮かび上がらせることができます。ところが、レーダーを航空機や衛星に載せて飛ばした瞬間、状況は一変します。プラットフォームが移動すると、レーダーから見える地面の各点はアンテナに対して相対速度を持ち、その相対速度は「アンテナから見た角度」に依存して変化します。正面の地面は速く近づき、真横の地面はほとんど近づきません。結果として、地面という本来静止している物体からのクラッタが、角度に応じて連続的に異なるドップラー周波数を持つようになります。これが「クラッタリッジ」と呼ばれる、角度ドップラー平面上に斜めに走る帯状のクラッタです。

このクラッタリッジは厄介です。単純なドップラーフィルタ(MTI)でゼロ周波数だけを落としても、角度によってはクラッタが目標と同じドップラー周波数に重なってしまい、消し切れません。逆に、空間フィルタ(ビームフォーミング)で特定の角度だけを見ても、その角度のクラッタは目標と同じ方向から来るため分離できません。鍵は、角度(空間)とドップラー(時間)を同時に使った2次元のフィルタリングにあります。これが空間時間適応処理(Space-Time Adaptive Processing, STAP)です。

STAPを理解すると、現代レーダーの中核技術への見通しが一気に開けます。

  • 早期警戒機(AEW&C)の地上移動目標検出(GMTI): 高空から地表を走る車両を、強大な地表クラッタの中から検出するにはSTAPが不可欠です
  • 宇宙搭載レーダー(SBR): 衛星速度が約 7.5 km/s と非常に大きく、クラッタリッジが角度ドップラー平面全体に大きく広がるため、STAPの威力が最も発揮される応用です
  • 干渉・妨害(ジャミング)対策: 特定方向から来るジャマーは角度ドップラー平面上では「縦の帯」として現れ、STAPはこれもクラッタと同時に抑圧できます

本記事の内容

  • 移動レーダーでクラッタが角度依存のドップラーを持つ仕組み(クラッタリッジ)
  • アンテナ素子 × パルスの「空間時間スナップショット」というデータ構造
  • 空間時間共分散行列 $\bm{R}$ の構成と意味
  • 出力電力を最小化しつつ目標利得を保つ最適ウェイト $\bm{w} = \bm{R}^{-1}\bm{s}/(\bm{s}^H \bm{R}^{-1}\bm{s})$ の導出
  • 性能指標としての改善係数(SINR利得)の定義
  • Python実装: クラッタリッジ、適応ヌルの形成、改善係数の角度ドップラー可視化

前提知識

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

STAPは「空間(アンテナアレイ)の適応ビームフォーミング」と「時間(パルス列)のドップラー処理」を融合した技術です。前者はMVDRビームフォーミング、後者はパルスドップラー・MTIの記事で扱った内容が土台になります。STAPはこの2つを別々に行うのではなく、空間と時間を一体の長いベクトルとして扱う点に本質があります。

STAPが必要な理由 — 動くレーダーとクラッタリッジ

静止レーダーと移動レーダーの決定的な違い

まず、なぜ移動レーダーでは話が難しくなるのかを、具体的なイメージから掴みましょう。あなたがヘリコプターに乗って前方を見ているとします。真正面の景色は猛烈な速さで近づいてきますが、真横の景色はほとんど横に流れるだけで、近づいても遠ざかってもいません。レーダーが地面を「見る」ときも全く同じことが起こります。レーダーから放射された電波が地面で反射して戻ってくるとき、地面のその点とレーダーの相対速度(視線方向速度)が、反射波のドップラー周波数を決めます。

プラットフォームの速度を $v_p$、レーダーから地面のある点を見る方向の方位角を $\theta$(アンテナの進行方向=正面を $\theta=0$ とする)とすると、その点の視線方向速度は $v_p \cos\theta$ です。正面($\theta=0$)では $v_p$ そのもの、真横($\theta = 90^\circ$)では $0$ になります。ドップラー周波数 $f_d$ は視線方向速度に比例するので、

$$ f_d = \frac{2 v_p \cos\theta}{\lambda} $$

ここで $\lambda$ は波長、係数 $2$ は往復(送受信両方でドップラーを受ける)に由来します。重要なのは、地面という静止物体からのクラッタなのに、ドップラー周波数 $f_d$ が見る角度 $\theta$ の関数になっているという点です。静止レーダー($v_p=0$)なら全角度で $f_d=0$ ですが、動くレーダーでは角度ごとにクラッタのドップラーがずれていきます。

クラッタリッジ — 角度ドップラー平面の斜めの帯

この関係を「横軸に空間周波数(角度に対応)、縦軸に正規化ドップラー周波数」を取った平面で描くと、クラッタは1本の斜めの線(厳密には少し幅のある帯)に乗ります。これがクラッタリッジです。

レーダー工学では角度を直接使う代わりに、アンテナ素子間の位相差に対応する空間周波数 $f_s$ を使うと見通しが良くなります。素子間隔を $d$ とすると、方位 $\theta$ の方向から来る波は隣り合う素子間で位相差 $2\pi (d/\lambda)\sin\theta$ を生むので、正規化空間周波数を

$$ f_s = \frac{d}{\lambda}\sin\theta $$

と定義します。一方、正規化ドップラー周波数はパルス繰り返し周波数 $f_r$(PRF)で割って

$$ f_t = \frac{f_d}{f_r} = \frac{2 v_p \cos\theta}{\lambda f_r} $$

と定義します。ここで $\beta$ というパラメータを導入します。1パルスの間(時間 $1/f_r$)にプラットフォームが進む距離を素子間隔 $d$ の何個分かで測った量で、

$$ \beta = \frac{2 v_p}{f_r d} $$

と置くと、$f_t = \beta (d/\lambda)\cos\theta$ と書けます。$\sin^2\theta + \cos^2\theta = 1$ の関係から空間周波数とドップラーの間に

$$ \left(\frac{f_s}{d/\lambda}\right)^2 + \left(\frac{f_t}{\beta\, d/\lambda}\right)^2 = 1 $$

という楕円の関係が成り立ちます。特に $\beta=1$(1パルスでちょうど素子間隔1個分進む配置)のときは $f_t = \beta f_s$、つまり $f_t = f_s$ となり、クラッタリッジは角度ドップラー平面上で傾き1の直線になります。この「クラッタが斜めの直線に乗る」という性質こそ、STAPが活躍する舞台です。

なぜ1次元処理では消せないのか

ここで「なぜ空間か時間の片方だけでは不十分か」を明確にしておきましょう。仮に目標が空間周波数 $f_{s0}$、ドップラー $f_{t0}$ にいるとします。

  • ドップラーフィルタ(時間のみ): ドップラー $f_{t0}$ を通すフィルタを使うと、同じドップラーを持つクラッタ(リッジ上で $f_t=f_{t0}$ の点、対応する角度)も一緒に通ってしまいます。
  • ビームフォーマ(空間のみ): 目標方向 $f_{s0}$ にビームを向けると、その方向のクラッタ(リッジ上で $f_s=f_{s0}$ の点)も一緒に通ります。

つまり、クラッタリッジは角度ドップラー平面に広がった「2次元的な存在」なので、片方の軸だけを切るフィルタでは必ずリッジと交差してしまいます。リッジを避けて目標だけを通すには、角度ドップラー平面上でリッジに沿った細いヌル(除去帯)を作るしかありません。これを実現するのが、空間と時間を同時に扱うSTAPです。

クラッタリッジという「敵の正体」が見えたので、次はそれを攻略するためのデータ構造、すなわち空間時間スナップショットを定義しましょう。

空間時間スナップショットとデータモデル

データの並べ方

STAPの出発点は、受信データをどう並べるかにあります。$N$ 素子のフェーズドアレイアンテナで、1回のコヒーレント処理間隔(CPI)に $M$ パルスを送受信するとします。あるレンジゲート(注目する距離セル)に対して、各素子・各パルスの複素受信信号が1つずつ得られます。つまりデータは $N \times M$ 個の複素数です。

これを縦に積み上げて長さ $NM$ の1本の列ベクトルにしたものを、空間時間スナップショット $\bm{x}$ と呼びます。

$$ \bm{x} = \begin{bmatrix} x_{1,1} \\ x_{2,1} \\ \vdots \\ x_{N,1} \\ x_{1,2} \\ \vdots \\ x_{N,M} \end{bmatrix} \in \mathbb{C}^{NM} $$

ここで $x_{n,m}$ は第 $n$ 素子・第 $m$ パルスの受信信号です。素子インデックスを内側、パルスインデックスを外側にして並べています(逆順でも理論は同じですが、後のステアリングベクトルの形が決まります)。このように空間(素子)と時間(パルス)を1本のベクトルに「縦結合」することが、STAPの根幹です。$N$ 次元の空間問題と $M$ 次元の時間問題を別々に解くのではなく、$NM$ 次元の1つの問題として解くわけです。

目標信号のステアリングベクトル

次に、ある角度・ドップラーから来る信号がこのスナップショット上でどう見えるかを表すベクトルを作ります。これを空間時間ステアリングベクトル(または操舵ベクトル)$\bm{s}$ と呼びます。

空間方向を考えます。方位 $\theta$ から平面波が来ると、素子 $n$ には基準素子に対して位相 $2\pi (n-1)(d/\lambda)\sin\theta$ の遅れが生じます。正規化空間周波数 $f_s = (d/\lambda)\sin\theta$ を使うと、空間ステアリングベクトルは

$$ \bm{a}(f_s) = \begin{bmatrix} 1 \\ e^{j2\pi f_s} \\ e^{j2\pi \cdot 2 f_s} \\ \vdots \\ e^{j2\pi (N-1) f_s} \end{bmatrix} \in \mathbb{C}^N $$

となります。これは通常のアレイのステアリングベクトルそのものです。

時間方向も同じ構造です。パルスごとに目標のドップラーによる位相回転が積み重なります。パルス間隔は $1/f_r$ なので、パルス $m$ には位相 $2\pi (m-1) f_t$($f_t = f_d/f_r$ は正規化ドップラー)が乗ります。時間(ドップラー)ステアリングベクトルは

$$ \bm{b}(f_t) = \begin{bmatrix} 1 \\ e^{j2\pi f_t} \\ \vdots \\ e^{j2\pi (M-1) f_t} \end{bmatrix} \in \mathbb{C}^M $$

です。

空間時間スナップショットは「素子内側・パルス外側」で並べたので、全体のステアリングベクトルは時間ベクトルと空間ベクトルのクロネッカー積で書けます。

$$ \bm{s}(f_s, f_t) = \bm{b}(f_t) \otimes \bm{a}(f_s) \in \mathbb{C}^{NM} $$

クロネッカー積 $\bm{b}\otimes\bm{a}$ は、$\bm{b}$ の各成分に $\bm{a}$ 全体を掛けて縦に並べたものです。第 $m$ パルスのブロックが $e^{j2\pi(m-1)f_t}\bm{a}(f_s)$ になっており、まさに「各パルスで空間ステアリングベクトルがドップラー位相だけ回転している」というデータの作りと一致します。この $\bm{s}(f_s, f_t)$ が、角度ドップラー平面上の1点 $(f_s, f_t)$ に対応する「探したい信号の形」です。

受信信号の総和モデル

実際のスナップショット $\bm{x}$ は、目標・クラッタ・ジャマー・雑音の和としてモデル化します。

$$ \bm{x} = \alpha_t\, \bm{s}(f_{s0}, f_{t0}) + \bm{c} + \bm{j} + \bm{n} $$

各項の意味は次の通りです。第1項 $\alpha_t \bm{s}(f_{s0}, f_{t0})$ が複素振幅 $\alpha_t$ の目標信号で、角度ドップラー $(f_{s0}, f_{t0})$ に位置します。$\bm{c}$ がクラッタ、$\bm{j}$ がジャマー、$\bm{n}$ が熱雑音です。クラッタ $\bm{c}$ は、クラッタリッジ上に分布した無数の地面散乱体からの寄与の和です。すなわち、リッジ上の各角度 $\theta_i$(対応する $f_{s,i}, f_{t,i}$)から来る散乱の合計として

$$ \bm{c} = \sum_i \gamma_i\, \bm{s}(f_{s,i}, f_{t,i}) $$

と書けます。$\gamma_i$ は各散乱体のランダムな複素反射係数です。これらの「目標以外の成分(クラッタ+ジャマー+雑音)」を干渉と総称し、その統計的性質を1つの行列に凝縮したものが、次に説明する空間時間共分散行列です。

データの形と信号モデルが揃ったので、いよいよ干渉の構造を表す共分散行列を構成します。

空間時間共分散行列

共分散行列の定義と意味

適応処理の心臓部は、干渉成分がどの方向・どのドップラーにどれだけのエネルギーを持ち、互いにどう相関しているかを記述する空間時間共分散行列 $\bm{R}$ です。目標を含まない干渉のみのスナップショット $\bm{x}_{\text{in}} = \bm{c}+\bm{j}+\bm{n}$ に対して、

$$ \bm{R} = \mathbb{E}\!\left[\bm{x}_{\text{in}}\, \bm{x}_{\text{in}}^H\right] \in \mathbb{C}^{NM \times NM} $$

と定義します。ここで $(\cdot)^H$ は共役転置(エルミート転置)、$\mathbb{E}[\cdot]$ は期待値です。$\bm{R}$ は $NM \times NM$ のエルミート半正定値行列で、対角成分は各空間時間チャネルの平均パワー、非対角成分はチャネル間の相関を表します。

なぜこの行列が重要なのか、直感的に言えば「干渉のエネルギーがどの方向に集中しているか」を表すからです。ベクトル $\bm{u}$ に沿った干渉のパワーは $\bm{u}^H \bm{R}\, \bm{u}$ で測れます。クラッタリッジの方向にあたる $\bm{u}$ ではこの値が非常に大きく、リッジから外れた方向では小さくなります。STAPはこの $\bm{R}$ を逆に使って、エネルギーが大きい方向を自動的に「避ける」ように動きます。

各成分の共分散構造

$\bm{R}$ を構成する各成分の寄与を見ておきましょう。互いに無相関と仮定すると、共分散は和になります。

$$ \bm{R} = \bm{R}_c + \bm{R}_j + \bm{R}_n $$

雑音 は通常、空間時間の各チャネルで独立同分布の白色雑音として $\bm{R}_n = \sigma_n^2 \bm{I}$ と表されます($\bm{I}$ は $NM$ 次の単位行列)。これは角度ドップラー平面全体に一様に広がるフロアです。

クラッタ は、リッジ上の散乱体の寄与の和なので、各散乱体パワーを $\sigma_{c,i}^2$ とすると

$$ \bm{R}_c = \sum_i \sigma_{c,i}^2\, \bm{s}(f_{s,i}, f_{t,i})\, \bm{s}(f_{s,i}, f_{t,i})^H $$

となります。リッジ上の各点が「ランク1」の寄与 $\bm{s}\bm{s}^H$ を出し、それらの和がクラッタ共分散になります。重要な性質として、$\beta$ が整数のとき、クラッタ共分散の実効ランクは約 $N + (M-1)\beta$ に収まることが知られています(Brennanの法則)。$NM$ 次元の空間の中で、クラッタは低次元の部分空間にしか広がらないのです。この「クラッタが低ランク」という事実が、STAPで少ない自由度を残してクラッタを除去できる理由です。

ジャマー は特定方向 $\theta_J$ から来る広帯域雑音で、全ドップラーにエネルギーが広がる一方、空間方向には1つの角度に集中します。角度ドップラー平面では「縦の帯」、共分散としては

$$ \bm{R}_j = \sigma_J^2\, \left(\bm{I}_M \otimes \bm{a}(f_{sJ})\bm{a}(f_{sJ})^H\right) $$

のような構造(時間方向には相関がなく、空間方向にランク1)を持ちます。

実際の推定 — サンプル共分散行列

理論上の $\bm{R}$ は期待値ですが、実際には得られません。代わりに、目標が存在しないと考えられる近隣のレンジゲート(訓練データ)$\bm{x}_k$, $k=1,\dots,K$ を集めて標本平均で推定します。

$$ \hat{\bm{R}} = \frac{1}{K}\sum_{k=1}^{K} \bm{x}_k\, \bm{x}_k^H $$

これをサンプル共分散行列と呼びます。Reed–Mallett–Brennan(RMB)の規則として、検出損失を 3 dB 以内に抑えるには訓練サンプル数 $K$ が自由度 $NM$ の約 2 倍以上必要、という経験則が知られています。$NM$ が大きいと膨大な均質訓練データが要るため、後述の部分空間STAPなどで自由度を減らす工夫が実用上重要になります。

共分散行列という「干渉の地図」が手に入りました。次は、この地図を頼りに目標だけを通すための最適なウェイトを導出します。

最適ウェイトの導出

問題設定 — MVDR的な定式化

STAPの出力は、スナップショット $\bm{x}$ にウェイトベクトル $\bm{w}\in\mathbb{C}^{NM}$ を掛けて和を取ったスカラー

$$ y = \bm{w}^H \bm{x} $$

です。私たちが望むのは2つのことです。第一に、目標方向 $(f_{s0}, f_{t0})$ から来る信号は歪みなく通したい。第二に、それ以外の干渉(クラッタ・ジャマー・雑音)の出力パワーはできるだけ小さくしたい。これはMVDR(最小分散無歪応答)ビームフォーミングと全く同じ発想を、空間時間の $NM$ 次元に拡張したものです。

第一の要請(無歪)は、目標ステアリングベクトル $\bm{s}\equiv\bm{s}(f_{s0},f_{t0})$ に対する応答を $1$ に固定する制約として書けます。

$$ \bm{w}^H \bm{s} = 1 $$

第二の要請(干渉最小化)は、干渉出力パワー

$$ \mathbb{E}[|y_{\text{in}}|^2] = \mathbb{E}[|\bm{w}^H \bm{x}_{\text{in}}|^2] = \bm{w}^H \bm{R}\, \bm{w} $$

を最小化することです。まとめると、次の制約付き最小化問題になります。

$$ \min_{\bm{w}}\ \bm{w}^H \bm{R}\, \bm{w} \quad \text{subject to}\quad \bm{w}^H \bm{s} = 1 $$

ラグランジュの未定乗数法による解

この問題をラグランジュの未定乗数法で解きます。$\bm{w}$ は複素ベクトルなので、ウィルティンガー微分($\bm{w}^*$ に関する微分)を使うのが手早い方法です。複素ラグランジアンを

$$ \mathcal{L}(\bm{w}, \lambda) = \bm{w}^H \bm{R}\, \bm{w} – \lambda\left(\bm{w}^H \bm{s} – 1\right) $$

と置きます。ここで $\bm{R}$ はエルミートなので $\bm{w}^H\bm{R}\bm{w}$ は実数です。$\bm{w}^*$ について偏微分してゼロと置くと、第1項の微分は $\bm{R}\bm{w}$、第2項の微分は $\lambda \bm{s}$ となるので、

$$ \frac{\partial \mathcal{L}}{\partial \bm{w}^*} = \bm{R}\, \bm{w} – \lambda\, \bm{s} = 0 $$

が得られます。ここから $\bm{w}$ について解くと、$\bm{R}$ が正定値(雑音項 $\sigma_n^2\bm{I}$ のおかげで可逆)であることを使って両辺に左から $\bm{R}^{-1}$ を掛け、

$$ \bm{w} = \lambda\, \bm{R}^{-1}\bm{s} $$

となります。残るは未定乗数 $\lambda$ の決定です。これを制約 $\bm{w}^H\bm{s}=1$ に代入します。$\bm{w}^H = \lambda^* \bm{s}^H \bm{R}^{-H} = \lambda^* \bm{s}^H \bm{R}^{-1}$($\bm{R}$ はエルミートなので $\bm{R}^{-H}=\bm{R}^{-1}$)に注意すると、

$$ \bm{w}^H \bm{s} = \lambda^*\, \bm{s}^H \bm{R}^{-1}\bm{s} = 1 $$

ここで $\bm{s}^H \bm{R}^{-1}\bm{s}$ は実の正数($\bm{R}^{-1}$ が正定値だから)です。したがって $\lambda^* = 1/(\bm{s}^H\bm{R}^{-1}\bm{s})$、すなわち $\lambda$ も同じ実数値です。これを $\bm{w}=\lambda\bm{R}^{-1}\bm{s}$ に戻すと、STAPの最適ウェイトが得られます。

$$ \boxed{\ \bm{w}_{\text{opt}} = \frac{\bm{R}^{-1}\bm{s}}{\bm{s}^H \bm{R}^{-1}\bm{s}}\ } $$

これがSTAPの中心公式であり、MVDR解の空間時間版です。分子の $\bm{R}^{-1}\bm{s}$ が「干渉を避けつつ目標方向を向く」方向を決め、分母 $\bm{s}^H\bm{R}^{-1}\bm{s}$ が無歪制約 $\bm{w}^H\bm{s}=1$ を満たすための正規化です。

ウェイトの直感 — なぜクラッタにヌルが立つのか

公式の中身を直感的に解釈しましょう。もし干渉が雑音だけ($\bm{R}=\sigma_n^2\bm{I}$)なら、$\bm{R}^{-1}\propto\bm{I}$ なので $\bm{w}\propto\bm{s}$、つまり通常の整合フィルタ(目標方向にビームを向けるだけ)になります。

ところがクラッタやジャマーが存在すると、$\bm{R}$ はそれらの方向に大きな固有値を持ちます。$\bm{R}^{-1}$ はその大きな固有値を小さくします(逆数になる)。したがって $\bm{R}^{-1}\bm{s}$ は、干渉の強い方向の成分を抑え込んだウェイトになります。結果として、出来上がったビームパターン(角度ドップラー応答)はクラッタリッジやジャマー方向に深いヌルを形成します。これが「適応的にヌルを置く」というSTAPの本質的な働きです。データ($\bm{R}$)が干渉の場所を教え、ウェイトがそこを自動的に避けるのです。

最適ウェイトが求まったので、次に「どれだけ性能が上がったか」を定量化する指標、改善係数を定義します。

改善係数(SINR利得)

出力SINRの定義

STAPの性能を測る基本量は、出力における信号対干渉雑音比(SINR)です。ウェイト $\bm{w}$ を通した後、目標信号のパワーと干渉のパワーの比として定義します。目標成分は $\bm{w}^H(\alpha_t\bm{s}) = \alpha_t\,\bm{w}^H\bm{s}$ なので、その平均パワーは $\sigma_t^2 |\bm{w}^H\bm{s}|^2$($\sigma_t^2=\mathbb{E}[|\alpha_t|^2]$)です。干渉成分のパワーは $\bm{w}^H\bm{R}\bm{w}$ です。したがって、

$$ \text{SINR}_{\text{out}} = \frac{\sigma_t^2\, |\bm{w}^H \bm{s}|^2}{\bm{w}^H \bm{R}\, \bm{w}} $$

最適ウェイト $\bm{w}_{\text{opt}}=\bm{R}^{-1}\bm{s}/(\bm{s}^H\bm{R}^{-1}\bm{s})$ を代入してみましょう。まず分子の $\bm{w}^H\bm{s}$ は無歪制約により $1$ です。分母は

$$ \bm{w}_{\text{opt}}^H \bm{R}\, \bm{w}_{\text{opt}} = \frac{\bm{s}^H\bm{R}^{-1}\,\bm{R}\,\bm{R}^{-1}\bm{s}}{(\bm{s}^H\bm{R}^{-1}\bm{s})^2} = \frac{\bm{s}^H\bm{R}^{-1}\bm{s}}{(\bm{s}^H\bm{R}^{-1}\bm{s})^2} = \frac{1}{\bm{s}^H\bm{R}^{-1}\bm{s}} $$

と簡単になります(途中で $\bm{R}^{-1}\bm{R}\bm{R}^{-1}=\bm{R}^{-1}$ を使いました)。よって最適SINRは

$$ \text{SINR}_{\text{opt}} = \sigma_t^2\, \bm{s}^H \bm{R}^{-1}\bm{s} $$

という非常にきれいな形になります。スカラー $\bm{s}^H\bm{R}^{-1}\bm{s}$ が大きいほど、その角度ドップラー点での検出性能が良いことを意味します。

改善係数の定義

実用上は、「STAPによって雑音だけのときと比べてどれだけSINRが改善したか」を表す改善係数(Improvement Factor, IF)を使います。これは出力SINRと入力SINRの比として定義されます。

$$ \text{IF} = \frac{\text{SINR}_{\text{out}}}{\text{SINR}_{\text{in}}} $$

入力SINRは1チャネルあたりの目標パワーと干渉パワーの比で、慣例的に $\text{SINR}_{\text{in}}=\sigma_t^2/\sigma_n^2$ を基準に取ります(雑音支配の参照)。すると、最適STAPの改善係数は

$$ \text{IF} = \sigma_n^2\, \bm{s}^H \bm{R}^{-1}\bm{s} $$

と書けます。さらに比較しやすいよう、雑音のみの理想ケース($\bm{R}=\sigma_n^2\bm{I}$)で得られる最大SINR $\text{SINR}_0 = \sigma_t^2\, NM/\sigma_n^2$ で正規化したSINR損失 $L_{\text{SINR}}$ もよく使われます。

$$ L_{\text{SINR}}(f_s,f_t) = \frac{\text{SINR}_{\text{opt}}(f_s,f_t)}{\text{SINR}_0} = \frac{\sigma_n^2\,\bm{s}^H\bm{R}^{-1}\bm{s}}{NM} $$

$L_{\text{SINR}}$ は $0$ から $1$($0$ dB)の値を取り、$1$ に近いほど「干渉がない理想状態と同じ性能」を意味します。クラッタリッジ上やジャマー方向では $L_{\text{SINR}}$ が深く落ち込み(検出できない)、リッジから離れたドップラーでは $1$ に近づきます。この $L_{\text{SINR}}$ を角度ドップラー平面(あるいは固定角度でドップラー軸)にプロットした曲線が、STAPの性能を一目で示す「SINR損失曲線」です。

最小可検出速度(MDV)との関係

SINR損失曲線がクラッタリッジ付近で落ち込む幅は、レーダーが検出できる目標の最小速度、すなわち最小可検出速度(MDV)を決めます。落ち込みが鋭く狭いほど、クラッタに近い遅い目標まで検出できることになります。空間時間の自由度 $NM$ を大きくするほど、リッジ周りのヌルを鋭くできるため、MDVを小さく(性能を良く)できます。これがSTAPで多素子・多パルスを使う動機です。

性能指標が定義できたので、いよいよPythonで全体を実装し、クラッタリッジ・適応ヌル・SINR損失を可視化して理論を確かめましょう。

Pythonによる実装と可視化

シナリオ設定とステアリングベクトル

まず、移動レーダーのパラメータと、空間時間ステアリングベクトルを作る関数を用意します。$\beta=1$ となるよう PRF・速度を選び、クラッタリッジが傾き1の直線になる典型シナリオを作ります。

import numpy as np
import matplotlib.pyplot as plt

# --- レーダーパラメータ ---
N = 8        # アンテナ素子数(空間自由度)
M = 8        # CPI内のパルス数(時間自由度)
d_over_lam = 0.5   # 素子間隔/波長(半波長間隔)
beta = 1.0   # クラッタリッジの傾きを決めるパラメータ(=1で傾き1)

def steering(fs, ft, N, M):
    """空間時間ステアリングベクトル s = b(ft) ⊗ a(fs) を返す(長さ NM)"""
    n = np.arange(N)
    m = np.arange(M)
    a = np.exp(1j * 2 * np.pi * fs * n)   # 空間ステアリング
    b = np.exp(1j * 2 * np.pi * ft * m)   # 時間(ドップラー)ステアリング
    return np.kron(b, a)                  # クロネッカー積で NM 次元に

# 動作確認: 長さが NM か
s_test = steering(0.1, 0.2, N, M)
print("steering vector length:", s_test.shape[0], " (expected", N * M, ")")

このコードは、空間ステアリング $\bm{a}(f_s)$ と時間ステアリング $\bm{b}(f_t)$ をクロネッカー積で結合し、長さ $NM=64$ の空間時間ステアリングベクトルを作っています。出力の長さが $64$ と表示されれば、データモデルがコード上で正しく構成できていることが確認できます。この steering 関数が以降のすべての処理(共分散構成・ウェイト計算・可視化)の基礎部品になります。

クラッタ共分散行列の構成

次に、クラッタリッジ上に散乱体を多数並べてクラッタ共分散 $\bm{R}_c$ を作り、ジャマーと雑音を加えて全干渉共分散 $\bm{R}$ を構成します。

import numpy as np

# --- 共分散行列の構成 ---
NM = N * M
CNR_dB = 50.0    # クラッタ対雑音比 [dB]
JNR_dB = 30.0    # ジャマー対雑音比 [dB]
sigma_n2 = 1.0   # 雑音電力(基準)
cnr = 10 ** (CNR_dB / 10)
jnr = 10 ** (JNR_dB / 10)

# クラッタリッジ: 角度 theta を一様にサンプリングし、各点の (fs, ft) を計算
Nc = 200                                  # クラッタ散乱体の数(リッジ上)
theta = np.linspace(-np.pi/2, np.pi/2, Nc)
fs_c = d_over_lam * np.sin(theta)          # 空間周波数
ft_c = beta * d_over_lam * np.cos(theta)   # 正規化ドップラー(ft = beta*(d/lam)*cosθ)

R = sigma_n2 * np.eye(NM, dtype=complex)   # まず白色雑音 σn^2 I
clutter_power = cnr * sigma_n2 / Nc        # 各散乱体に均等配分
for fsi, fti in zip(fs_c, ft_c):
    s = steering(fsi, fti, N, M)
    R += clutter_power * np.outer(s, s.conj())   # ランク1寄与 s s^H の和

# ジャマー: 角度 30度から、全ドップラーに広がる(空間ランク1)
theta_J = np.deg2rad(30.0)
fsJ = d_over_lam * np.sin(theta_J)
aJ = np.exp(1j * 2 * np.pi * fsJ * np.arange(N))
Rj = jnr * sigma_n2 * np.kron(np.eye(M), np.outer(aJ, aJ.conj()))
R += Rj

print("R shape:", R.shape, " Hermitian error:", np.max(np.abs(R - R.conj().T)))

このコードは、リッジ上の200個の散乱体それぞれについてランク1の寄与 $\bm{s}\bm{s}^H$ を積算してクラッタ共分散を作り、$30^\circ$ 方向のジャマー(全ドップラーに広がる空間ランク1の干渉)と白色雑音を加えています。「Hermitian error」がほぼ $0$ と表示されることから、$\bm{R}$ がエルミート行列として正しく構成できていることが確認できます。クラッタ対雑音比を 50 dB と非常に大きく設定しているのは、現実の地表クラッタが目標よりもはるかに強いという状況を反映するためです。

クラッタの角度ドップラー分布(適応前)

STAPをかける前に、干渉がどの角度ドップラー点に集中しているかを見ておきましょう。各点 $(f_s,f_t)$ でのクラッタパワー $\bm{s}^H\bm{R}\bm{s}$ をマップ表示します。

import numpy as np
import matplotlib.pyplot as plt

# 角度ドップラー平面のグリッド
Ns_grid, Nt_grid = 120, 120
fs_grid = np.linspace(-0.5, 0.5, Ns_grid)   # 空間周波数
ft_grid = np.linspace(-0.5, 0.5, Nt_grid)   # 正規化ドップラー

power_map = np.zeros((Nt_grid, Ns_grid))
for i, ft in enumerate(ft_grid):
    for j, fs in enumerate(fs_grid):
        s = steering(fs, ft, N, M)
        # s^H R s = この方向に存在する干渉パワー
        power_map[i, j] = np.real(s.conj() @ (R @ s))

plt.figure(figsize=(7, 6))
plt.pcolormesh(fs_grid, ft_grid, 10 * np.log10(power_map / power_map.min()),
               cmap='jet', shading='auto')
plt.colorbar(label='Interference power [dB]')
plt.xlabel('Normalized spatial frequency  $f_s$')
plt.ylabel('Normalized Doppler  $f_t$')
plt.title('Interference power before STAP (clutter ridge + jammer)')
plt.tight_layout()
plt.savefig('stap_clutter_ridge.png', dpi=150, bbox_inches='tight')
plt.show()

このマップから、STAPが立ち向かう干渉の構造が一目で読み取れます。第一に、左下から右上へ斜めに走る明るい帯がクラッタリッジで、$\beta=1$ のために傾きがほぼ1の直線になっています。これは「角度が変わるとクラッタのドップラーも比例して動く」という移動レーダーの性質そのものです。第二に、$f_s = 0.5\sin30^\circ = 0.25$ 付近に縦に走る帯が見え、これが全ドップラーに広がるジャマーです。クラッタが斜めの帯、ジャマーが縦の帯という2次元的な広がりを持つため、横軸だけ・縦軸だけのフィルタでは除去できないことが視覚的に納得できます。

最適ウェイトの計算と適応ビームパターン

いよいよ最適ウェイトを計算します。目標を $f_s=0$(正面方向)に固定し、各ドップラー $f_t$ ごとにウェイトを作って、その適応応答(ビームパターン)を見ます。

import numpy as np

# 共分散の逆行列(数値安定化のため微小対角ロードを加える)
Rinv = np.linalg.inv(R + 1e-3 * np.eye(NM))

def stap_weight(fs0, ft0):
    """目標 (fs0, ft0) に対する最適ウェイト w = Rinv s / (s^H Rinv s)"""
    s = steering(fs0, ft0, N, M)
    Rinv_s = Rinv @ s
    return Rinv_s / (s.conj() @ Rinv_s)

# 例: 目標を fs0=0, ft0=0.3 に置いたときのウェイト
w_example = stap_weight(0.0, 0.3)
# 目標応答が 1 (無歪制約) になっているか確認
s0 = steering(0.0, 0.3, N, M)
print("w^H s =", np.round(w_example.conj() @ s0, 6), " (expected ~1)")

w^H s1 と表示されることから、最適ウェイトが無歪制約 $\bm{w}^H\bm{s}=1$ を厳密に満たしていることが確認できます。微小な対角ロード($10^{-3}\bm{I}$)を加えているのは、クラッタが低ランクで $\bm{R}$ が悪条件になりがちなため、逆行列計算を数値的に安定させる実務的な工夫です。続いて、このウェイトが角度ドップラー平面にどんなヌルを作るかを描きます。

import numpy as np
import matplotlib.pyplot as plt

# 目標を fs0=0, ft0=0.3 に固定したウェイトの 2D 応答パターン
fs0, ft0 = 0.0, 0.3
w = stap_weight(fs0, ft0)

resp = np.zeros((Nt_grid, Ns_grid))
for i, ft in enumerate(ft_grid):
    for j, fs in enumerate(fs_grid):
        s = steering(fs, ft, N, M)
        resp[i, j] = np.abs(w.conj() @ s) ** 2   # |w^H s|^2 = 適応応答

resp_dB = 10 * np.log10(resp / resp.max())

plt.figure(figsize=(7, 6))
plt.pcolormesh(fs_grid, ft_grid, resp_dB, cmap='viridis',
               shading='auto', vmin=-60, vmax=0)
plt.colorbar(label='Adapted response [dB]')
plt.plot(fs0, ft0, 'r*', markersize=16, label='target')
plt.xlabel('Normalized spatial frequency  $f_s$')
plt.ylabel('Normalized Doppler  $f_t$')
plt.title('STAP adapted beampattern (target at $f_s$=0, $f_t$=0.3)')
plt.legend(loc='upper right')
plt.tight_layout()
plt.savefig('stap_adapted_pattern.png', dpi=150, bbox_inches='tight')
plt.show()

この適応ビームパターンから、STAPの働きが鮮やかに見えます。赤い星印で示した目標位置($f_s=0,\ f_t=0.3$)では応答が $0$ dB(最大)になっており、目標が歪みなく通っています。一方、先ほどクラッタリッジがあった斜めの帯に沿って、応答が $-40$ dB 以下に深く落ち込んだ適応ヌルが形成されています。ジャマーがあった $f_s\approx0.25$ の縦の帯にも同様に深いヌルが見えます。つまり、ウェイトは $\bm{R}$ という干渉の地図を読み取り、自分で立ち向かうべき場所(クラッタリッジとジャマー)にだけピンポイントでヌルを置いているのです。誰もヌルの位置を手で指定していないのに、データから自動的に最適なヌルが立つ点がSTAPの真価です。

SINR損失曲線 — 改善係数の可視化

最後に、性能指標であるSINR損失 $L_{\text{SINR}}(f_s,f_t)$ をドップラー軸に沿ってプロットし、クラッタリッジ付近で性能がどれだけ落ちるか(=最小可検出速度)を確認します。

import numpy as np
import matplotlib.pyplot as plt

# 目標角度を fs0=0 に固定し、ドップラー ft を掃引して SINR損失を計算
fs0 = 0.0
ft_sweep = np.linspace(-0.5, 0.5, 300)
Lsinr = np.zeros_like(ft_sweep)

for k, ft in enumerate(ft_sweep):
    s = steering(fs0, ft, N, M)
    sinr_opt = np.real(s.conj() @ (Rinv @ s))   # σt^2 を1とした s^H Rinv s
    Lsinr[k] = sinr_opt / NM                     # 雑音のみ理想 (NM) で正規化

Lsinr_dB = 10 * np.log10(Lsinr)

plt.figure(figsize=(8, 5))
plt.plot(ft_sweep, Lsinr_dB, lw=2, color='darkorange')
plt.axhline(0, color='gray', ls='--', lw=1, label='ideal (no interference)')
plt.xlabel('Normalized Doppler  $f_t$')
plt.ylabel('SINR loss  $L_{SINR}$ [dB]')
plt.title('STAP SINR loss vs. Doppler (target at $f_s$=0)')
plt.ylim(-50, 5)
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('stap_sinr_loss.png', dpi=150, bbox_inches='tight')
plt.show()

このSINR損失曲線から、STAPの定量的な性能が読み取れます。第一に、ドップラーの大部分(クラッタリッジから離れた領域)で $L_{\text{SINR}}$ が $0$ dB 付近に張り付いており、干渉がない理想状態とほぼ同じSINRが得られています。第二に、$f_s=0$ に対応するクラッタのドップラー($\beta=1,\ \theta=90^\circ$ 付近で $f_t\approx0$)近傍に鋭く深い落ち込みがあり、ここが「クラッタに埋もれて検出できない速度域」、すなわち最小可検出速度に対応します。落ち込みの幅が狭いほど、遅い目標まで検出できることを意味します。第三に、もし素子数 $N$ やパルス数 $M$ を増やせば、この落ち込みはさらに鋭く狭くなり、MDVが改善します。これは空間時間の自由度 $NM$ がクラッタヌルの鋭さを決めるという理論と一致しています。

これら3つの可視化(クラッタリッジ、適応ヌル、SINR損失)を合わせて見ることで、STAPが「干渉の2次元分布を測り、そこにだけヌルを置き、目標は無歪で通す」という一貫した動作をしていることが、理論と実装の両面から確認できました。

まとめ

本記事では、移動レーダーのクラッタとジャミングを角度ドップラー2次元で抑圧する空間時間適応処理(STAP)について、原理から実装まで解説しました。

  • クラッタリッジ: プラットフォームが動くと、地面クラッタのドップラーが見る角度に依存し、角度ドップラー平面上で斜めの帯($\beta=1$ で傾き1の直線)になります。これが1次元フィルタでは消せない理由です。
  • 空間時間スナップショット: $N$ 素子 × $M$ パルスのデータを長さ $NM$ の1本のベクトルに縦結合し、ステアリングベクトルはクロネッカー積 $\bm{s}=\bm{b}(f_t)\otimes\bm{a}(f_s)$ で表します。
  • 共分散行列 $\bm{R}=\mathbb{E}[\bm{x}_{\text{in}}\bm{x}_{\text{in}}^H]$ が干渉の「地図」であり、クラッタは低ランク(Brennanの法則)に収まります。実用ではサンプル共分散 $\hat{\bm{R}}$ で推定します。
  • 最適ウェイト $\bm{w}_{\text{opt}}=\bm{R}^{-1}\bm{s}/(\bm{s}^H\bm{R}^{-1}\bm{s})$ をMVDR的な制約付き最小化から導出しました。$\bm{R}^{-1}$ が干渉方向を抑え、クラッタリッジとジャマーに適応ヌルを自動形成します。
  • 改善係数とSINR損失: 最適SINRは $\sigma_t^2\,\bm{s}^H\bm{R}^{-1}\bm{s}$ という簡潔な形になり、$L_{\text{SINR}}$ の落ち込み幅が最小可検出速度を決めます。

STAPは計算量が大きく($NM\times NM$ の逆行列)、訓練データも自由度の約2倍必要なため、実用では空間や時間の自由度を部分的に絞る「部分空間STAP」(要素空間後ドップラーSTAP、ビーム空間STAPなど)が広く使われます。これらは本記事のフルSTAPを土台に、性能と計算量のトレードオフを取る発展形です。

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