空港のレーダー画面には、1回転ごとに無数の輝点が現れます。そのうちのいくつかは旅客機ですが、残りは鳥だったり、雨雲だったり、単なる受信機の熱雑音がしきい値をたまたま超えただけの「誤警報」だったりします。しかもレーダーは「この輝点は日本航空123便です」という名札を付けてくれません。返ってくるのは、ただのラベルのない点の集合です。
ここで、追尾する側は毎スキャン次の問題を突きつけられます。今スキャンで得られた $m$ 個の観測のうち、どれが追尾中のどの航跡の続きなのか？ これがデータアソシエーション（data association、観測・航跡の対応付け）の問題です。カルマンフィルタは「観測がどの目標のものかわかっている」という前提で組まれた道具なので、この前提が崩れた瞬間に、どんなに精緻なフィルタ設計も無意味になります。実際、多目標追尾システムの性能を最終的に決めるのはフィルタのチューニングではなく、このアソシエーションの巧拙です。
データアソシエーションを理解すると、次のような場面で何が起きているのかが見通せるようになります。
- 航空管制・空域監視: 二次監視レーダー（SSR）が使えない非協力目標を、一次レーダーの点群だけから追尾し続ける処理の中核です
- 自動運転のミリ波レーダー / LiDAR: 前方の車列を個別に追い続けるとき、隣り合う車両の反射点が入れ替わると衝突判断そのものが破綻します
- 宇宙状況認識（SSA）: 地上レーダーがデブリを観測すると、1パスで数十個の未識別点が得られ、これをカタログ上の物体と結び付ける必要があります
- 海上監視・ソナー: 波浪クラッタが密な環境で小型船舶を追尾する典型的な応用先です
本記事では、この問題に対する古典的かつ現役の3つの解法 — 最近傍法（NN）、確率的データ関連（PDA）、同時確率的データ関連（JPDA） — を、数式の導出を省略せずに追いかけます。最後にPythonで、交差する2目標と誤警報が混在するシナリオを作り、NNが航跡の入れ替わりを起こす一方でPDA・JPDAが持ちこたえる様子を定量的に確認します。
本記事の内容
- Track-While-Scan（TWS）における観測の曖昧さと、3つの困難
- イノベーション $\bm{v}$ とその共分散 $\bm{S}$ の導出
- マハラノビス距離 $\bm{v}^\top \bm{S}^{-1} \bm{v}$ がカイ二乗分布に従う理由と、ゲーティングの設計
- ゲート確率 $P_G$、ゲート体積、ゲート内の期待誤警報数のトレードオフ
- 最近傍法（NN）と大域最近傍法（GNN）— 尤度からコスト関数を導く
- PDA の関連確率 $\beta_i$ をベイズ則から導出する
- PDA の共分散更新に現れる「イノベーションの散らばり」項の導出
- PDA の弱点（航跡の併合）と、JPDA の同時事象確率による解決
- Python 実装: 交差2目標シナリオでの NN / GNN / PDA / JPDA の比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
Track-While-Scan とデータアソシエーション問題とは
まず、この記事がこれから扱う問題の全体像を1枚の絵にしておきます。
図の ✚ は追尾中の2本の航跡が「次のスキャンではここに現れるはず」と予言している位置、破線の楕円はその予言の不確かさから決まる探索範囲（ゲート）です。× で描かれた検出点には名札がなく、$z_2$ のように両方のゲートに同時に入ってしまう点が存在します。この $z_2$ を航跡1に渡すか航跡2に渡すか、あるいはどちらでもない誤警報とみなすか — この判断こそが本記事の主題であり、$z_5$ のようにどちらのゲートにも入らない点を機械的に捨てる仕組み（ゲーティング）から順に組み立てていきます。
昔のレーダーは、アンテナを1つの目標に向けたまま追い続ける「単一目標追尾（STT）」が主流でした。これなら「今アンテナが向いている方向にいるのが、追いかけている目標だ」と決め打ちできるので、アソシエーションの問題は生じません。しかしアンテナを占有してしまうため、1度に1機しか追えません。
そこで登場したのが Track-While-Scan（TWS、捜索しながら追尾） です。アンテナは通常どおり空域を回転走査し、1回転（1スキャン、周期 $T$）ごとに、その回転中に検出されたすべての目標候補をまとめて信号処理部に渡します。これで同時に何十もの目標を追えるようになりました。その代償として払うことになったのが、まさにデータアソシエーションのコストです。
1スキャンで何が渡されるのか
1スキャンで得られるのは、位置（と場合によってはドップラー速度）を表すベクトルの集合です。
$$ Z_k = \{\bm{z}_{k,1}, \bm{z}_{k,2}, \dots, \bm{z}_{k,m_k}\} $$
ここが決定的に重要なのですが、この集合には順序に意味がありません。$\bm{z}_{k,1}$ が目標Aで $\bm{z}_{k,2}$ が目標Bだという保証はどこにもなく、番号は単に検出処理が出力した順番にすぎません。つまり $Z_k$ は「順序なしの点集合」であって、通常のカルマンフィルタが想定する「時刻 $k$ の観測ベクトル $\bm{z}_k$」ではないのです。この構造のギャップを埋めるのがアソシエーションの役目です。
3つの困難
観測の曖昧さは、次の3つの要因が絡み合って生じます。
3つのパネルは、どれも「ゲートの中に何が入っているか」という同じ視点で描いてあります。左は目標が確かにいるのに検出が1つも出ない状態、中央はゲート内外に誤警報が散らばって本物が埋もれる状態、右は2目標のゲートが重なって観測の帰属が決められない状態です。この3つは独立ではなく同時に起こりうるので、実際のゲート内は「本物が入っているかどうかすらわからない数個の点」という状況になります。
(1) 検出漏れ（$P_D < 1$） — CFAR検出器は、しきい値を超えた信号だけを検出として出力します。目標のRCSがフェージングで落ち込んだスキャンでは検出が出ません。実用的な設定では $P_D = 0.7 \sim 0.95$ 程度で、「目標は確かにそこにいるのに観測がない」スキャンが必ず混じります。
(2) 誤警報（クラッタ） — 逆に、目標がいないのに検出が出ることもあります。熱雑音がしきい値を超えるケース、地面・海面・降雨からの反射、鳥や昆虫。これらは監視領域内にほぼ一様に散らばり、その個数はポアソン分布に従うとモデル化するのが標準です。単位面積あたりの期待個数をクラッタ密度 $\lambda$ と呼びます。
(3) 目標同士の混同 — もっとも厄介なのがこれです。2機の航空機が接近・交差すると、両者の予測位置が観測誤差の大きさと同程度まで近づきます。このとき「左の点が機体A、右の点が機体B」と言い切る根拠はほとんどありません。ここで誤ると、以後フィルタは別の目標を追い続けることになります。これを航跡の入れ替わり（track swap）と呼びます。
なぜ「入れ替わり」は致命的なのか
一度アソシエーションを誤ると、フィルタは間違った観測に引っ張られて状態推定を更新します。カルマンフィルタは自分の推定を信頼しているので、誤った更新の後は誤った予測を出し、その予測の周りにゲートを張り、また誤った観測を取り込みます。誤りが自己強化されるわけです。しかも共分散 $\bm{P}$ は観測を取り込むたびに小さくなるので、フィルタは「自信満々に間違える」状態に陥ります。これが、アソシエーションが単なる前処理ではなく多目標追尾の本丸である理由です。
ここまでで、なぜアソシエーションが必要かがわかりました。次は、この問題を扱うための土台となる追尾サイクルの構造と、そこで使う記法を整理しておきましょう。
追尾サイクルの全体像と記法
多目標追尾の1スキャンは、次の4段階で進みます。
- 予測（prediction） — 前スキャンの推定を運動モデルで1周期進め、「次はここに現れるはず」という予測 $\hat{\bm{x}}_{k|k-1}$ と、その不確かさ $\bm{P}_{k|k-1}$ を得る
- ゲーティング（gating） — 予測位置の周りに「ありえる範囲」を張り、その外の観測を候補から機械的に除外する（計算量の削減と誤対応の防止）
- アソシエーション（association） — ゲート内に残った候補の中から、どれをどの航跡に割り当てるか（あるいはどう重み付けるか）を決める
- 更新（update） — 決まった対応に基づいて状態推定を更新する
NN・PDA・JPDA が違うのは、このうち3番目だけです。1・2・4は共通の枠組みで動きます。
この図で押さえてほしいのは、赤く塗った第3段階だけが手法によって差し替わり、その前後は共通の部品だという構造です。観測集合 $Z_k$ が入ってくるのは第2段階のゲーティングで、そこから第4段階の更新までの間に「どの観測を、どれだけの重みで使うか」が決まります。だからこそ本記事のコードでも、予測・ゲーティング・更新の関数は1つずつしか書かず、アソシエーション部分だけを4通り用意する構成になります。
状態空間モデル
以降、目標は等速直線運動（CV モデル）に従うとします。状態を位置と速度からなる4次元ベクトルとし、
$$ \bm{x}_k = \begin{bmatrix} x_k & \dot{x}_k & y_k & \dot{y}_k \end{bmatrix}^\top $$
状態遷移とプロセス雑音を次のように書きます。
$$ \bm{x}_k = \bm{F}\bm{x}_{k-1} + \bm{w}_k, \quad \bm{w}_k \sim \mathcal{N}(\bm{0}, \bm{Q}) $$
観測は位置のみが得られるとし、
$$ \bm{z}_k = \bm{H}\bm{x}_k + \bm{n}_k, \quad \bm{n}_k \sim \mathcal{N}(\bm{0}, \bm{R}) $$
行列は次のとおりです。$T$ はスキャン周期です。
$$ \bm{F} = \begin{bmatrix} 1 & T & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & T \\ 0 & 0 & 0 & 1 \end{bmatrix}, \quad \bm{H} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} $$
そして、$Z^{k} = \{Z_1, Z_2, \dots, Z_k\}$ で時刻 $k$ までの全スキャンの累積を表します。追尾の目的は、事後密度 $p(\bm{x}_k \mid Z^k)$ を（近似的に）求めることです。
この枠組みの中で、まず「予測がどれだけ外れうるか」を定量化する量 — イノベーションとその共分散 — を導きましょう。ゲーティングもアソシエーションも、すべてこの2つの量の上に建てられます。
イノベーションとその共分散の導出
追尾フィルタが持っている「予言」は、次スキャンで観測がどこに現れるかという予測観測値です。
$$ \hat{\bm{z}}_{k|k-1} = \bm{H}\hat{\bm{x}}_{k|k-1} $$
実際に得られた観測 $\bm{z}_k$ と、この予言のズレをイノベーション（innovation、新機軸）と呼びます。
$$ \bm{v}_k = \bm{z}_k – \hat{\bm{z}}_{k|k-1} $$
「イノベーション」という名前は、この差分こそが観測にしか含まれていない新情報だ、という意味です。予測と一致した部分はフィルタがすでに知っていたことなので、情報量はゼロです。
共分散 $\bm{S}$ を導く
イノベーションがどれくらいばらつくかを知りたいので、共分散を計算します。まず、予測誤差を $\tilde{\bm{x}}_{k|k-1} = \bm{x}_k – \hat{\bm{x}}_{k|k-1}$ と置き、これを使ってイノベーションを書き換えます。$\bm{z}_k = \bm{H}\bm{x}_k + \bm{n}_k$ を代入すると、
$$ \bm{v}_k = \bm{H}\bm{x}_k + \bm{n}_k – \bm{H}\hat{\bm{x}}_{k|k-1} $$
右辺の $\bm{H}$ でくくれる部分をまとめると、予測誤差が現れます。
$$ \bm{v}_k = \bm{H}(\bm{x}_k – \hat{\bm{x}}_{k|k-1}) + \bm{n}_k = \bm{H}\tilde{\bm{x}}_{k|k-1} + \bm{n}_k $$
つまりイノベーションは「予測誤差が観測空間に写ったもの」と「観測雑音」の和です。ここで共分散を取ります。$\mathbb{E}[\bm{v}_k] = \bm{0}$ なので、
$$ \bm{S}_k = \mathbb{E}\left[\bm{v}_k \bm{v}_k^\top\right] = \mathbb{E}\left[(\bm{H}\tilde{\bm{x}}_{k|k-1} + \bm{n}_k)(\bm{H}\tilde{\bm{x}}_{k|k-1} + \bm{n}_k)^\top\right] $$
括弧を展開すると4つの項が出ます。交差項 $\mathbb{E}[\bm{H}\tilde{\bm{x}}_{k|k-1}\bm{n}_k^\top]$ は、時刻 $k$ の観測雑音が過去の情報から作られた予測誤差と独立なのでゼロになります。これが効いて、残るのは対角的な2項だけです。
$$ \bm{S}_k = \bm{H}\,\mathbb{E}[\tilde{\bm{x}}_{k|k-1}\tilde{\bm{x}}_{k|k-1}^\top]\,\bm{H}^\top + \mathbb{E}[\bm{n}_k\bm{n}_k^\top] $$
期待値の中身はそれぞれ予測共分散 $\bm{P}_{k|k-1}$ と観測雑音共分散 $\bm{R}$ そのものなので、
$$ \begin{equation} \bm{S}_k = \bm{H}\bm{P}_{k|k-1}\bm{H}^\top + \bm{R} \end{equation} $$
が得られました。この式は「観測が現れうる散らばりは、自分の予測がどれだけ怪しいか（$\bm{H}\bm{P}\bm{H}^\top$）とセンサがどれだけ雑なのか（$\bm{R}$）の足し算だ」と読めます。追尾を始めたばかりで $\bm{P}$ が大きいときは $\bm{S}$ も大きく、収束後は $\bm{S} \to \bm{R}$ に近づきます。
左の図では、青い楕円（予測の怪しさ）とオレンジの楕円（センサの雑さ）を足し合わせた赤い楕円が $\bm{S}$ にあたり、実際の観測 $\bm{z}_k$ はこの赤い楕円の広がりの中に現れます。緑の矢印がイノベーション $\bm{v}_k$ です。右の図は本記事のパラメータ（$q=0.05$、$\sigma_m=2$ m）で $\sqrt{\det\bm{S}}$ の推移を追ったもので、10スキャンほどで定常値 6.42 に落ち着きます。この値がセンサ単体の下限 $\sqrt{\det\bm{R}} = 4$ より大きいのは、プロセス雑音のぶんだけ予測に不確かさが残り続けるからです。
イノベーションの分布
線形ガウスモデルの下では、$\tilde{\bm{x}}_{k|k-1}$ も $\bm{n}_k$ もガウス分布に従い、その和もガウス分布です。したがって、観測 $\bm{z}_k$ が本当にこの目標由来であるという条件の下で、
$$ \bm{v}_k \mid Z^{k-1} \sim \mathcal{N}(\bm{0}, \bm{S}_k) $$
が成り立ちます。この「条件の下で」が重要です。誤警報や他目標由来の観測に対しては、この分布はまったく成り立ちません。だからこそ、$\bm{v}$ の大きさを測れば「この観測はこの目標のものらしいか」を判定できるわけです。この発想を具体化したのが次節のゲーティングです。
ゲーティング — マハラノビス距離とカイ二乗分布
なぜユークリッド距離ではダメなのか
「予測位置から半径 $r$ 以内の観測を候補にする」という素朴な方法を考えてみます。しかしこれは2つの意味で不適切です。第一に、$r$ をどう決めればよいかの根拠がありません。第二に、そして本質的に、予測の不確かさは方向によって違うからです。等速直線運動する目標を追尾していると、進行方向の不確かさ（速度誤差が距離誤差に化ける）は横方向より大きくなります。予測共分散の等高線は円ではなく、進行方向に長い楕円になっています。
この楕円の形を考慮して距離を測るのがマハラノビス距離です。イノベーション共分散 $\bm{S}$ を使って、
$$ \begin{equation} d^2 = \bm{v}^\top \bm{S}^{-1} \bm{v} \end{equation} $$
と定義します。共分散の大きい方向では $\bm{S}^{-1}$ の固有値が小さくなるため、同じ物理距離でも $d^2$ は小さく評価されます。「ばらつきやすい方向のズレは大目に見る」という、極めて自然な重み付けです。
この比較が効いています。観測Aはユークリッド距離で 7.5 m、観測Bは 5.3 m ですから、素朴な円ゲートではBの方が「近い」ことになります。ところが $\bm{S}$ の長軸方向を考慮した $d^2$ で測ると、Aが 2.2 でゲート内、Bが 16.2 でゲート外と順位がまるごと逆転します。予測が不確かな方向に沿ったズレは目標由来としてありふれており、逆に自信のある方向へのズレは異常だ — という判断が、$\bm{S}^{-1}$ という重みだけで自動的に実現されているわけです。
$d^2$ がカイ二乗分布に従うことの導出
さて、この $d^2$ の分布がわかれば、しきい値を確率的な根拠を持って決められます。導出しましょう。
$\bm{S}$ は共分散行列なので対称正定値です。したがってコレスキー分解 $\bm{S} = \bm{L}\bm{L}^\top$（$\bm{L}$ は下三角行列で対角成分が正）が一意に存在します。ここで白色化された変数を導入します。
$$ \bm{u} = \bm{L}^{-1}\bm{v} $$
この $\bm{u}$ の分布を調べます。$\bm{v} \sim \mathcal{N}(\bm{0}, \bm{S})$ で $\bm{u}$ は $\bm{v}$ の線形変換なので、$\bm{u}$ もガウス分布に従い、平均は $\bm{L}^{-1}\bm{0} = \bm{0}$ です。共分散は線形変換の公式 $\mathrm{Cov}(\bm{A}\bm{v}) = \bm{A}\,\mathrm{Cov}(\bm{v})\,\bm{A}^\top$ を使うと、
$$ \mathrm{Cov}(\bm{u}) = \bm{L}^{-1}\bm{S}\bm{L}^{-\top} $$
ここに $\bm{S} = \bm{L}\bm{L}^\top$ を代入すると、$\bm{L}^{-1}$ と $\bm{L}$、$\bm{L}^\top$ と $\bm{L}^{-\top}$ がそれぞれ打ち消し合います。
$$ \mathrm{Cov}(\bm{u}) = \bm{L}^{-1}\bm{L}\bm{L}^\top\bm{L}^{-\top} = \bm{I} $$
つまり $\bm{u} \sim \mathcal{N}(\bm{0}, \bm{I})$ — 各成分が独立な標準正規分布です。次に $d^2$ を $\bm{u}$ で書き直します。$\bm{v} = \bm{L}\bm{u}$ と $\bm{S}^{-1} = (\bm{L}\bm{L}^\top)^{-1} = \bm{L}^{-\top}\bm{L}^{-1}$ を $d^2$ の定義に代入すると、
$$ d^2 = (\bm{L}\bm{u})^\top \bm{L}^{-\top}\bm{L}^{-1} (\bm{L}\bm{u}) = \bm{u}^\top \bm{L}^\top \bm{L}^{-\top}\bm{L}^{-1}\bm{L}\,\bm{u} $$
中央の4つの行列が単位行列に潰れるので、
$$ d^2 = \bm{u}^\top\bm{u} = \sum_{i=1}^{n_z} u_i^2 $$
これは独立な標準正規変数の二乗和そのものです。カイ二乗分布の定義そのままなので、
$$ \begin{equation} d^2 = \bm{v}^\top\bm{S}^{-1}\bm{v} \sim \chi^2_{n_z} \end{equation} $$
が示せました。$n_z$ は観測ベクトルの次元です。2次元の位置観測なら $n_z = 2$、距離・方位・仰角なら $n_z = 3$ です。
ゲート確率 $P_G$ とゲートしきい値 $\gamma$
しきい値 $\gamma$ を決めて、$d^2 \le \gamma$ を満たす観測だけを候補に残すことにします。目標由来の観測がこのゲートに入る確率をゲート確率 $P_G$ と呼び、カイ二乗分布の累積分布関数で与えられます。
$$ P_G = P(d^2 \le \gamma) = F_{\chi^2_{n_z}}(\gamma) $$
$n_z = 2$ の場合、カイ二乗分布の累積分布関数は初等関数で書けます。自由度2のカイ二乗分布は平均2の指数分布と一致するので、
$$ P_G = 1 – e^{-\gamma/2} \quad \Longleftrightarrow \quad \gamma = -2\ln(1 – P_G) $$
具体的な値を見ておきましょう。
| $P_G$ | $\gamma$（$n_z=2$） | $\gamma$（$n_z=3$） |
|---|---|---|
| 0.90 | 4.605 | 6.251 |
| 0.95 | 5.991 | 7.815 |
| 0.99 | 9.210 | 11.345 |
| 0.997 | 11.618 | 13.932 |
実務では $P_G = 0.99$ 前後がよく使われます。ここを欲張って $P_G = 0.999$ にすると、後述するようにゲートに入る誤警報が増えてしまいます。
ゲート体積と、その中に紛れ込む誤警報の数
$d^2 \le \gamma$ という条件は、観測空間で超楕円体の内部を表します。この体積を計算します。$\bm{u}$ 座標では条件は $\|\bm{u}\|^2 \le \gamma$、つまり半径 $\sqrt{\gamma}$ の $n_z$ 次元球です。$\bm{v} = \bm{L}\bm{u}$ という線形変換のヤコビアンは $|\det \bm{L}| = \sqrt{\det \bm{S}}$ なので、元の空間での体積は球の体積を $\sqrt{\det\bm{S}}$ 倍したものになります。
$$ \begin{equation} V_G = c_{n_z}\,\gamma^{n_z/2}\sqrt{\det \bm{S}} \end{equation} $$
ここで $c_{n_z}$ は単位球の体積で、$c_1 = 2, c_2 = \pi, c_3 = 4\pi/3$ です。2次元では $V_G = \pi\gamma\sqrt{\det\bm{S}}$ となります。
クラッタ密度 $\lambda$ の環境では、ゲート内に入る誤警報の期待個数は
$$ \mathbb{E}[\text{ゲート内誤警報数}] = \lambda V_G $$
です。ここに、ゲーティング設計の核心的なトレードオフが現れます。
- $\gamma$ を大きくする → $P_G$ が上がり、真の観測を取りこぼしにくくなる
- しかし $V_G \propto \gamma^{n_z/2}$ で膨らむ → 紛れ込む誤警報が増え、誤対応の確率が上がる
しかも $\bm{S}$ が大きい（＝予測に自信がない）ときほど $V_G$ も大きくなるので、フィルタが迷子になりかけているときほど誤警報を拾いやすいという悪循環の構造があります。これが、追尾の失敗が急激に進行する理由の1つです。
左の図で、$P_G$ の曲線（青）は $\gamma = 6$ を超えたあたりから寝てしまうのに、期待誤警報数（橙の破線）は $\gamma$ に線形に増え続けます。得られるものが飽和する一方でコストだけが伸び続ける、という典型的な収穫逓減の形です。右の図はこれを $P_G$ 側から見たもので、$P_G$ を 0.99 から 0.999 に上げるとゲート体積は 1.50 倍に膨らむのに、取りこぼしの改善は 1 % → 0.1 % にすぎません。$P_G = 0.99$ 前後が実務の定番になっているのは、まさにこの膝の位置だからです。
ゲートで候補を絞り込むところまでは、これから紹介する3手法すべてに共通です。ここから先、絞り込んだ候補をどう使うかで手法が分かれます。まずはもっとも素直な発想 — 一番近いものを選ぶ — から見ていきましょう。
最近傍法（NN）— 一番近いものを選ぶ
直感とアルゴリズム
いちばん自然な方針は「ゲート内でもっとも予測に近い観測が、その目標のものだろう」です。これが最近傍法（Nearest Neighbor, NN）です。ただし「近い」はユークリッド距離ではなくマハラノビス距離で測ります。
航跡 $t$ について、
$$ j^\star = \arg\min_{j:\, d^2_{t,j} \le \gamma} d^2_{t,j} $$
を選び、通常のカルマンフィルタ更新をそのまま適用します。
$$ \hat{\bm{x}}_{k|k} = \hat{\bm{x}}_{k|k-1} + \bm{W}_k \bm{v}_{k,j^\star}, \quad \bm{P}_{k|k} = \bm{P}_{k|k-1} – \bm{W}_k\bm{S}_k\bm{W}_k^\top $$
ここで $\bm{W}_k = \bm{P}_{k|k-1}\bm{H}^\top\bm{S}_k^{-1}$ はカルマンゲインです。ゲート内に観測が1つもなければ更新をスキップし、予測値をそのまま推定値とします（$\hat{\bm{x}}_{k|k} = \hat{\bm{x}}_{k|k-1}$）。
実装が単純で計算も軽く、$P_D$ が高くクラッタが薄く目標が十分離れている状況では、これで十分機能します。
尤度から見た NN — なぜ $d^2$ が正しいコストなのか
「一番近い」という直感的な基準は、実は最尤推定として正当化できます。観測 $\bm{z}_j$ が航跡 $t$ に由来するとしたときの尤度は、ガウス密度そのものです。
$$ \Lambda_{t,j} = \mathcal{N}(\bm{z}_j;\, \hat{\bm{z}}_t, \bm{S}_t) = \frac{1}{(2\pi)^{n_z/2}\sqrt{\det\bm{S}_t}} \exp\left(-\tfrac{1}{2}d^2_{t,j}\right) $$
最大化したいのは $\Lambda$ ですが、扱いやすいように対数を取って符号を反転し、コスト（最小化したい量）にします。
$$ -2\ln\Lambda_{t,j} = d^2_{t,j} + \ln\det\bm{S}_t + n_z\ln(2\pi) $$
右辺第3項は定数なので無視できます。第2項 $\ln\det\bm{S}_t$ は、航跡が1つだけなら（$\bm{S}_t$ が共通なので）これも定数となり、結局 $d^2$ を最小化することと最尤選択が一致します。これが「マハラノビス距離最小＝最尤」の正体です。
一方、複数の航跡を横断して比較するときには $\ln\det\bm{S}_t$ を落としてはいけません。共分散の大きい航跡（自信がない航跡）は、そもそもどの観測にもそこそこ高い $d^2$ を許してしまうので、$\ln\det\bm{S}_t$ のペナルティで補正する必要があります。これを忘れると、迷走している航跡が近くの観測を片っ端から奪う挙動になります。
大域最近傍法（GNN）— 排他性を入れる
素朴な NN には明白な欠陥があります。各航跡が独立に「自分に一番近い観測」を取るので、2つの航跡が同じ観測を取ってしまうことがあるのです。物理的には、1つの検出が2つの物体に同時に対応することはありえません。この矛盾を放置すると、2つの航跡が同じ点に吸い寄せられて重なってしまいます（航跡の併合、track coalescence）。
そこで、「1つの観測は高々1つの航跡にしか割り当てない」という排他制約を課し、全航跡・全観測をまとめて最適化するのが大域最近傍法（Global Nearest Neighbor, GNN）です。コスト行列を
$$ C_{t,j} = d^2_{t,j} + \ln\det\bm{S}_t \quad (\text{ゲート内のみ、外は} +\infty) $$
と定め、割当 $\pi$ の総コスト $\sum_t C_{t,\pi(t)}$ を最小化します。これは線形割当問題なので、ハンガリアン法（Munkres法）や Jonker-Volgenant 法で $O(n^3)$ で厳密に解けます。「非割当」（検出漏れ）を許すために、しきい値 $\gamma$ をコストとするダミー列を用意するのが定石です。
この配置がまさに NN の弱点を突いています。マハラノビス距離は航跡1から見て $z_1$ が 2.67、$z_2$ は 19.26 でゲート外。航跡2から見ると $z_1$ が 1.63、$z_2$ が 2.25 で、どちらもゲート内だが $z_1$ の方がわずかに近い。そのため NN では両航跡が $z_1$ を選び、1つの検出が2つの物体に対応するという物理的にありえない結論になります。GNN は総コスト $2.67 + 2.25 = 4.92$ の組み合わせを選ぶので、航跡2はわずかに遠い $z_2$ を引き受けて排他性が保たれます。「個々に最適」と「全体で最適」が食い違う典型例です。
NN 系の根本的な弱み
GNN で排他性を入れても、まだ本質的な問題が残ります。ハードな決定をしていることです。2つの候補の尤度が 0.51 対 0.49 でも、NN/GNN は 0.51 の方を確率1で採用し、0.49 の可能性を完全に捨てます。そして捨てた情報は二度と戻ってきません。
さらに悪いことに、更新後の共分散は $\bm{P}_{k|k} = \bm{P}_{k|k-1} – \bm{W}\bm{S}\bm{W}^\top$ と、「正しい観測を使った」前提の値まで縮みます。つまり間違えたときほど、間違いに自信を持つのです。目標が交差する数スキャンの間だけ判断を保留できれば、離れた後に自然に解けるはずの問題を、NN は無理に即断して固定してしまいます。
では、「決めない」という選択肢はあるのでしょうか。実は、あります。次節で見る PDA は、まさに「決めずに確率で重み付ける」という発想です。
確率的データ関連（PDA）— 決めずに重み付ける
直感
PDA（Probabilistic Data Association）の発想は明快です。どれが正しいかわからないなら、それぞれの「正しさの確率」で重み付けて全部使う。0.51 対 0.49 なら、0.51 と 0.49 の重みで両方の情報を取り込みます。判断を保留したまま前に進めるので、後のスキャンで曖昧さが解ければ自然に正解へ収束します。
ベイズ的に言えば、これはアソシエーションという離散的な隠れ変数を周辺化（marginalize）していることに他なりません。
前提と事象の定義
PDA は「追尾対象は1つ、それ以外の観測はすべてクラッタ」という前提で組まれています。時刻 $k$ でゲート内に残った観測を $\bm{z}_1, \dots, \bm{z}_m$（$m = m_k$）とし、次の $m+1$ 個の互いに排反な事象を考えます。
$$ \theta_i = \begin{cases} \text{観測 } \bm{z}_i \text{ が目標由来}, & i = 1, \dots, m \\ \text{目標由来の観測はゲート内に存在しない}, & i = 0 \end{cases} $$
$\theta_0$ は「検出漏れした」または「検出はされたがゲートの外に落ちた」の両方を含みます。求めたいのは、これらの事後確率 — 関連確率（association probability） — です。
$$ \beta_i = P(\theta_i \mid Z^k) $$
$\beta_i$ をベイズ則から導く
ベイズの定理を、時刻 $k$ の観測集合 $Z_k$ を条件として適用します。
$$ \beta_i = P(\theta_i \mid Z_k, Z^{k-1}) = \frac{p(Z_k \mid \theta_i, Z^{k-1})\, P(\theta_i \mid Z^{k-1})}{\sum_{j=0}^{m} p(Z_k \mid \theta_j, Z^{k-1})\, P(\theta_j \mid Z^{k-1})} $$
分子の2つの因子をそれぞれ評価していきます。
尤度 $p(Z_k \mid \theta_i, Z^{k-1})$ の評価。 $\theta_i$（$i \ge 1$）の下では、$\bm{z}_i$ は目標由来なので密度 $\mathcal{N}(\bm{z}_i; \hat{\bm{z}}, \bm{S})$ に従います。ただしゲート内に入ったという条件付きなので、正規化のため $P_G$ で割ります。残りの $m-1$ 個はクラッタで、ゲート内に一様分布すると仮定するので、それぞれ密度 $1/V_G$ を持ちます。
$$ p(Z_k \mid \theta_i, Z^{k-1}) = \frac{1}{P_G}\mathcal{N}(\bm{z}_i;\hat{\bm{z}},\bm{S}) \cdot V_G^{-(m-1)}, \quad i \ge 1 $$
$\theta_0$ の下では $m$ 個すべてがクラッタなので、
$$ p(Z_k \mid \theta_0, Z^{k-1}) = V_G^{-m} $$
事前確率 $P(\theta_i \mid Z^{k-1})$ の評価。 クラッタ数がパラメータ $\lambda V_G$ のポアソン分布に従うとする「パラメトリック PDA」で計算します。少し場合分けが要りますが、結論だけ先に言うと、$i \ge 1$ に対して事前確率は $i$ によらず等しくなり、
$$ P(\theta_i \mid Z^{k-1}) = \frac{P_D P_G}{m}\mu_F(m-1)\Big/ \mathcal{C}, \quad P(\theta_0 \mid Z^{k-1}) = (1 – P_D P_G)\mu_F(m)\big/ \mathcal{C} $$
の形になります。ここで $\mu_F(n) = e^{-\lambda V_G}(\lambda V_G)^n/n!$ はクラッタ数の確率、$\mathcal{C}$ は正規化定数です。「目標が検出されかつゲートに入る」確率が $P_D P_G$、そのとき残りのクラッタは $m-1$ 個、そしてどの番号に付くかが $1/m$ — という読み方をすれば自然な形です。
両者を掛けて整理する。 分子（$i \ge 1$）を書き下すと、$\mathcal{C}$ は分母と約分されて消えます。
$$ p(Z_k\mid\theta_i)P(\theta_i) \propto \frac{P_D P_G}{m}\,\frac{\mathcal{N}(\bm{z}_i;\hat{\bm{z}},\bm{S})}{P_G}\, V_G^{-(m-1)}\,\frac{e^{-\lambda V_G}(\lambda V_G)^{m-1}}{(m-1)!} $$
$V_G^{-(m-1)}$ と $(\lambda V_G)^{m-1}$ を掛け合わせると $V_G$ が完全に打ち消えて $\lambda^{m-1}$ だけが残ります。ここが導出の山場です。
$$ p(Z_k\mid\theta_i)P(\theta_i) \propto \frac{P_D}{m}\,\mathcal{N}(\bm{z}_i;\hat{\bm{z}},\bm{S})\,\frac{e^{-\lambda V_G}\lambda^{m-1}}{(m-1)!} $$
同様に $i = 0$ の分子は、
$$ p(Z_k\mid\theta_0)P(\theta_0) \propto (1-P_D P_G)\,V_G^{-m}\,\frac{e^{-\lambda V_G}(\lambda V_G)^{m}}{m!} = (1-P_D P_G)\,\frac{e^{-\lambda V_G}\lambda^{m}}{m!} $$
共通因子を落とす。 両者に共通する $e^{-\lambda V_G}\lambda^{m-1}/m!$ で割ると（$1/(m-1)! = m/m!$ を使います）、驚くほどきれいな形になります。
$$ \begin{equation} \beta_i = \frac{P_D\,\mathcal{N}(\bm{z}_i;\hat{\bm{z}},\bm{S})}{\lambda(1 – P_D P_G) + P_D\sum_{j=1}^{m}\mathcal{N}(\bm{z}_j;\hat{\bm{z}},\bm{S})}, \quad i \ge 1 \end{equation} $$
$$ \begin{equation} \beta_0 = \frac{\lambda(1 – P_D P_G)}{\lambda(1 – P_D P_G) + P_D\sum_{j=1}^{m}\mathcal{N}(\bm{z}_j;\hat{\bm{z}},\bm{S})} \end{equation} $$
もちろん $\sum_{i=0}^{m}\beta_i = 1$ です。
$\beta_i$ の式を読む
左の図で、予測位置から各候補に伸びる矢印の太さが $\beta_i$ の大きさです。予測にもっとも近い $z_1$（$d^2 = 0.2$）が $\beta_1 = 0.550$ と最大の重みを得る一方、少し離れた $z_2, z_3$ も 0.280 / 0.169 とそれなりの重みを保っています。NN ならここで $z_1$ 以外を完全に捨てるところを、PDA は捨てずに全部使うわけです。ゲート外の $z_4$（$d^2 = 12.5$）だけが $\beta_4 = 0$ になり、ゲーティングが確率の世界でもそのまま機能していることがわかります。赤い太矢印が、これらを重み付き平均した $\bar{\bm{v}}$ です。
この式には、追尾の常識がすべて詰まっています。
- 分子の $\mathcal{N}(\bm{z}_i;\hat{\bm{z}},\bm{S})$ — 予測に近い観測ほど大きな重みを得ます。これは NN と同じ直感です
- 分母の $\lambda(1-P_DP_G)$ — クラッタ密度 $\lambda$ が高いほど、$\beta_0$（＝どれも目標ではない）が大きくなります。汚い環境では、フィルタは自動的に観測を信用しなくなり、予測を維持する方向に振れます
- $P_D$ が高いほど $\beta_0$ が小さくなります。「検出されるはずの目標が検出されないのは変だ」という推論が働きます
- $m = 0$（ゲート内に観測なし）のときは $\beta_0 = 1$ となり、更新は行われません。式が自動的に正しい極限を与えます
左の図は $\beta_0$ をクラッタ密度の関数として描いたもので、候補が予測から離れている（$d^2 = 6$、点線）ほど早く $\beta_0$ が立ち上がります。本記事の設定である 5.8 個/スキャンの位置ではまだ $\beta_0$ はほぼゼロですが、100個を超えると 0.3〜0.5 に達し、フィルタは観測より予測を信じるようになります。右の図は同じ量を $P_D$ 別に描いたもので、$P_D = 0.99$ の曲線が一貫して最下段にあります。「検出されるはずの目標が見つからないのはおかしい」という推論が、$1 – P_D P_G$ という係数ひとつで表現されているわけです。
$\lambda$ を推定できない場合は、$\lambda$ を $m/V_G$（実際に観測された個数から密度を推定）で置き換えるノンパラメトリック PDA を使います。実装上はこちらの方が扱いやすいことも多いです。
状態の更新
各仮説 $\theta_i$ の下での更新推定を $\hat{\bm{x}}_i$ とします。$i \ge 1$ なら通常のカルマン更新、$i = 0$ なら予測そのままです。
$$ \hat{\bm{x}}_i = \hat{\bm{x}}_{k|k-1} + \bm{W}\bm{v}_i \ (i\ge1), \qquad \hat{\bm{x}}_0 = \hat{\bm{x}}_{k|k-1} $$
全確率の定理により、事後平均はこれらの $\beta_i$ による混合になります。
$$ \hat{\bm{x}}_{k|k} = \sum_{i=0}^{m}\beta_i \hat{\bm{x}}_i = \beta_0\hat{\bm{x}}_{k|k-1} + \sum_{i=1}^{m}\beta_i\left(\hat{\bm{x}}_{k|k-1} + \bm{W}\bm{v}_i\right) $$
$\sum_i \beta_i = 1$ を使って $\hat{\bm{x}}_{k|k-1}$ でくくると、非常にコンパクトな形にまとまります。
$$ \begin{equation} \hat{\bm{x}}_{k|k} = \hat{\bm{x}}_{k|k-1} + \bm{W}\bar{\bm{v}}, \qquad \bar{\bm{v}} = \sum_{i=1}^{m}\beta_i\bm{v}_i \end{equation} $$
つまり PDA の状態更新は、「重み付き平均イノベーション $\bar{\bm{v}}$ を1個の観測とみなした普通のカルマン更新」にすぎません。カルマンゲイン $\bm{W}$ も通常どおりです。実装上、既存のカルマンフィルタに数行足すだけで PDA になるのは、この構造のおかげです。
なお $\beta_0$ が大きい（どれも目標らしくない）とき、$\sum\beta_i \ (i\ge1)$ が小さくなるので $\bar{\bm{v}}$ は自動的に縮み、更新量が小さくなります。信用できない観測は、それ自体が弱められるという良い性質です。
ただし、平均だけを見ていると重大な見落としをします。次に共分散を正しく扱いましょう。ここに PDA のもう1つの本質があります。
PDA の共分散 — 「イノベーションの散らばり」項の導出
なぜ普通の共分散ではダメか
もし $\bar{\bm{v}}$ を「1個の確かな観測」として扱い、$\bm{P}_{k|k} = \bm{P}_{k|k-1} – \bm{W}\bm{S}\bm{W}^\top$ としたら何が起きるでしょうか。フィルタは「観測を1つ確実に取り込んだ」ことになり、共分散が正規の量だけ縮みます。しかし実際には、どれが正しいかわかっていないのです。この曖昧さの分だけ、推定はもっと不確かなはずです。この「関連の不確かさ」を共分散に反映させないと、フィルタは過信して発散します。
全分散の法則による分解
事後共分散を素直に定義から計算します。
$$ \bm{P}_{k|k} = \mathbb{E}\left[(\bm{x}_k – \hat{\bm{x}}_{k|k})(\bm{x}_k – \hat{\bm{x}}_{k|k})^\top \mid Z^k\right] $$
アソシエーション事象 $\theta_i$ で条件付けて、全分散の法則（law of total variance）を使います。混合分布の共分散は「各成分の共分散の重み付き平均」＋「各成分の平均のばらつき」で書けるという、よく知られた分解です。
$$ \bm{P}_{k|k} = \underbrace{\sum_{i=0}^{m}\beta_i \bm{P}_i}_{\text{各仮説内のばらつき}} + \underbrace{\sum_{i=0}^{m}\beta_i\left(\hat{\bm{x}}_i – \hat{\bm{x}}_{k|k}\right)\left(\hat{\bm{x}}_i – \hat{\bm{x}}_{k|k}\right)^\top}_{\text{仮説間のばらつき}} $$
第1項を計算する。 $\bm{P}_i$ は仮説 $i$ の下での共分散です。$i \ge 1$ では通常のカルマン更新の共分散 $\bm{P}^c = \bm{P}_{k|k-1} – \bm{W}\bm{S}\bm{W}^\top$、$i = 0$ では更新しないので $\bm{P}_{k|k-1}$ です。したがって、
$$ \sum_{i=0}^{m}\beta_i\bm{P}_i = \beta_0\bm{P}_{k|k-1} + (1-\beta_0)\bm{P}^c $$
「$\beta_0$ の確率で予測共分散のまま、$1-\beta_0$ の確率で縮んだ共分散」という、直感どおりの補間です。
第2項を計算する。 ここが「イノベーションの散らばり（spread of the innovations）」項です。まず $\hat{\bm{x}}_i – \hat{\bm{x}}_{k|k}$ を具体的に書きます。$i \ge 1$ では $\hat{\bm{x}}_i = \hat{\bm{x}}_{k|k-1} + \bm{W}\bm{v}_i$、そして $\hat{\bm{x}}_{k|k} = \hat{\bm{x}}_{k|k-1} + \bm{W}\bar{\bm{v}}$ なので、$\hat{\bm{x}}_{k|k-1}$ が打ち消し合って、
$$ \hat{\bm{x}}_i – \hat{\bm{x}}_{k|k} = \bm{W}(\bm{v}_i – \bar{\bm{v}}), \quad i \ge 1 $$
$i = 0$ では $\hat{\bm{x}}_0 = \hat{\bm{x}}_{k|k-1}$ なので、$\hat{\bm{x}}_0 – \hat{\bm{x}}_{k|k} = -\bm{W}\bar{\bm{v}}$ となります。これは上の式で $\bm{v}_0 = \bm{0}$ と置いたものと一致するので、$\bm{v}_0 \equiv \bm{0}$ という規約を導入すれば $i=0$ も同じ式に含められます。すると第2項は、
$$ \sum_{i=0}^{m}\beta_i \bm{W}(\bm{v}_i – \bar{\bm{v}})(\bm{v}_i – \bar{\bm{v}})^\top\bm{W}^\top = \bm{W}\left[\sum_{i=0}^{m}\beta_i(\bm{v}_i-\bar{\bm{v}})(\bm{v}_i-\bar{\bm{v}})^\top\right]\bm{W}^\top $$
括弧内を展開します。$\sum_i\beta_i = 1$ と $\sum_i \beta_i\bm{v}_i = \bar{\bm{v}}$（$\bm{v}_0=\bm{0}$ なので $i\ge1$ の和と同じ）を使うと、交差項が $-2\bar{\bm{v}}\bar{\bm{v}}^\top$ にまとまり、
$$ \sum_{i}\beta_i(\bm{v}_i-\bar{\bm{v}})(\bm{v}_i-\bar{\bm{v}})^\top = \sum_{i=1}^{m}\beta_i\bm{v}_i\bm{v}_i^\top – \bar{\bm{v}}\bar{\bm{v}}^\top $$
という、分散の定番の形（二乗の平均マイナス平均の二乗）に落ち着きます。以上をまとめると、PDA の共分散更新式が得られます。
$$ \begin{equation} \bm{P}_{k|k} = \beta_0\bm{P}_{k|k-1} + (1-\beta_0)\bm{P}^{c} + \tilde{\bm{P}} \end{equation} $$
$$ \begin{equation} \tilde{\bm{P}} = \bm{W}\left[\sum_{i=1}^{m}\beta_i\bm{v}_i\bm{v}_i^\top – \bar{\bm{v}}\bar{\bm{v}}^\top\right]\bm{W}^\top \end{equation} $$
$\tilde{\bm{P}}$ の意味
$\tilde{\bm{P}}$ は必ず半正定値です（分散の形をしているので）。つまり共分散を必ず増やす方向にはたらきます。その大きさは、ゲート内の候補イノベーションがどれだけばらついているかで決まります。
- 候補が1つだけ、あるいは候補が全部ほぼ同じ場所にある → $\bm{v}_i \approx \bar{\bm{v}}$ となり $\tilde{\bm{P}} \approx \bm{0}$。曖昧さがないので普通のカルマンフィルタと同じ
- 候補が離れた2点にあり、$\beta$ がほぼ半々 → $\tilde{\bm{P}}$ が大きくなり、共分散が大きく膨らむ
3つのパネルは同じ予測共分散から出発し、ゲート内の候補配置だけを変えたものです。左（候補1個）では $\tilde{\bm{P}}$ の位置トレースが 0.00 m$^2$ で、赤い PDA 更新後の楕円は青い通常カルマン更新の楕円とぴたり重なります。中央（候補2個が近接）でも 0.11 m$^2$ とわずかです。ところが右（候補2個が離れて $\beta$ が半々）では 8.95 m$^2$ に跳ね上がり、更新後の楕円が予測の楕円よりもむしろ大きくなっている方向すらあります。しかも膨らむ方向が2候補を結ぶ向きに揃っている点が重要で、「どちらか判断がつかない軸」に沿ってだけ不確かさが増しているのです。
「どっちかわからない」という認識そのものが、$\tilde{\bm{P}}$ という形で共分散に定量的に記録されるわけです。共分散が膨らめば次スキャンのゲートも広がり、フィルタは慎重になります。曖昧なときに慌てて確定させず、判断を先送りしながら不確かさを正直に持ち越す — これが PDA が NN より頑健である理由です。
さて、ここまでの PDA は「目標は1つ」という前提の上に建っていました。目標が複数あって近接するとき、この前提が破れるとどうなるでしょうか。次にその失敗と、その修正である JPDA を見ていきます。
PDA の限界 — 航跡の併合（track coalescence）
複数の目標がいるとき、素朴には「各目標に PDA を1つずつ走らせればよい」と考えたくなります。実際、目標同士が十分離れていればこれで問題なく動きます。しかし接近すると破綻します。
問題は、PDA が他の目標の存在を知らないことです。PDA から見れば、隣の目標が生んだ観測も「クラッタ」の一種です。ところがこの「クラッタ」は一様分布ではなく、まさに自分の予測位置の近くに集中して現れます。結果として、両方の航跡が両方の観測に高い $\beta$ を与え、どちらの航跡も2つの観測の中間点に引き寄せられます。両者が同じ場所へ引かれるので、2本の航跡は1本に重なってしまいます。これが航跡の併合（track coalescence）です。
しかも併合は入れ替わりより厄介です。入れ替わりなら少なくとも各航跡は「どこかの本物の目標」を追っており、ID を付け替えれば復旧できます。併合の場合はそうはいきません。典型的には次の2つの形で現れます。ひとつは、両方の航跡がどの目標でもない中間の幻に取り残されて、そのまま両方とも失探する形。もうひとつは、両方の航跡が同じ1つの目標に吸着し、もう一方の目標が誰にも追われないまま置き去りにされる形です。どちらの場合も「追われていない目標」が生まれるので、監視システムとしては最悪の失敗です。
根本原因ははっきりしています。「1つの観測は高々1つの目標から生じる」という排他性が、独立な PDA には組み込まれていないのです。これを入れるのが JPDA です。
JPDA — 同時事象で排他性を組み込む
同時事象（joint event）
$T$ 本の航跡と $m$ 個の観測があるとき、すべての観測の帰属を同時に指定した1つの割当を同時事象 $\theta$ と呼びます。$\theta$ は写像 $\theta: \{1,\dots,m\} \to \{0, 1, \dots, T\}$ として書けます（$\theta(j) = t$ は「観測 $j$ は航跡 $t$ 由来」、$\theta(j) = 0$ は「観測 $j$ はクラッタ」）。
実行可能（feasible）な同時事象とは、次の2条件を満たすものです。
- 1つの観測は高々1つの航跡に由来する（写像として定義されている時点で自動的に満たされます）
- 1つの航跡は高々1つの観測しか生まない — これが排他制約の本体です
条件2により、$\theta(j) = \theta(j’) = t \neq 0$（$j \neq j’$）は禁止されます。PDA を独立に走らせると、この禁止された組み合わせに確率を割り当ててしまっていたわけです。
同時事象の確率
実行可能事象 $\theta$ の事後確率は、PDA と同じ論法で計算できます。結果だけ書くと、
$$ \begin{equation} P(\theta \mid Z_k) = \frac{1}{c}\prod_{j:\,\theta(j)\neq 0} \frac{P_D^{(t_j)}\,\mathcal{N}\!\left(\bm{z}_j; \hat{\bm{z}}_{t_j}, \bm{S}_{t_j}\right)}{\lambda} \prod_{t:\,\text{未検出}}\left(1 – P_D^{(t)}P_G\right) \end{equation} $$
ここで $t_j = \theta(j)$、$c$ は正規化定数です。式の読み方は次のとおりです。
- 第1の積は、$\theta$ が「観測 $j$ は航跡 $t_j$ 由来」と主張している分の尤度比です。分母の $\lambda$ は「その観測がクラッタだったとしたときの密度」で、目標仮説とクラッタ仮説の尤度比になっています。$\lambda$ が大きい（クラッタが濃い）ほど、目標に帰属させる主張は弱くなります
- 第2の積は、$\theta$ の中で1つも観測を割り当てられなかった航跡に対するペナルティです。$1-P_D P_G$ は「検出されないかゲート外に落ちる確率」で、$P_D$ が高いほどこの値は小さく、「検出されないのはおかしい」という強いペナルティになります
周辺化して $\beta$ を得る
同時事象の確率が求まったら、個々の航跡についての関連確率は周辺化で得られます。
$$ \begin{equation} \beta_j^{t} = \sum_{\theta:\,\theta(j) = t} P(\theta \mid Z_k) \end{equation} $$
つまり「観測 $j$ を航跡 $t$ に割り当てているすべての実行可能事象の確率を足し上げる」わけです。$\beta_0^t$（航跡 $t$ が未検出である確率）も同様に、$t$ に何も割り当てていない事象の総和として得られます。
そして、ここからの更新式は PDA とまったく同じです。 $\beta$ を上の周辺化で計算するか PDA の式で計算するかだけが違い、$\bar{\bm{v}}$ による平均更新も $\tilde{\bm{P}}$ を含む共分散更新も一字一句同じ形を使います。JPDA の本質は「$\beta$ の計算方法を、排他制約を尊重したものに差し替える」ことに尽きます。
なぜこれで併合が防げるのか
2航跡・2観測の場面を考えてみましょう。実行可能な事象は主に「$\bm{z}_1 \to$ 航跡1、$\bm{z}_2 \to$ 航跡2」と「$\bm{z}_1 \to$ 航跡2、$\bm{z}_2 \to$ 航跡1」の2つです（未検出の組み合わせを除けば）。排他制約により、「両航跡ともに $\bm{z}_1$ に由来する」という事象は最初から存在しません。
したがって、航跡1が $\bm{z}_1$ に高い確率を与えるなら、その同じ事象の中で航跡2は $\bm{z}_2$ に高い確率を与えることになります。つまり 2本の航跡は互いに反発します。この「押し合い」が、両者が同じ中間点に落ちるのを防ぐのです。独立 PDA にはこの反発力がまったくありません。
具体例で確かめる
数値で見てみましょう。$P_D = 0.9$, $P_G = 0.99$, $\lambda = 10^{-3}$、そして2つの観測が両航跡のゲートに入っており、ガウス尤度が次のとおりだとします。
$$ \mathcal{N}_{1,1} = 0.020,\ \mathcal{N}_{1,2} = 0.010,\ \mathcal{N}_{2,1} = 0.010,\ \mathcal{N}_{2,2} = 0.020 $$
（$\mathcal{N}_{t,j}$ は航跡 $t$ に対する観測 $j$ の尤度。航跡1は $\bm{z}_1$ 寄り、航跡2は $\bm{z}_2$ 寄りという素直な配置です。）
独立 PDA では、航跡1について $\lambda(1-P_DP_G) = 10^{-3}\times(1-0.891) = 1.09\times10^{-4}$ が非検出項なので、
$$ \beta_1^{(1)} = \frac{0.9\times0.020}{1.09\times10^{-4} + 0.9\times(0.020+0.010)} = \frac{0.0180}{0.0271} = 0.664 $$
同様に $\beta_2^{(1)} = 0.332$、$\beta_0^{(1)} = 0.004$ です。航跡2は対称に $\beta_1^{(2)} = 0.332,\ \beta_2^{(2)} = 0.664$ となります。ここで注目すべきは、航跡1も航跡2も、$\bm{z}_1$ と $\bm{z}_2$ の両方から 0.33 以上の重みで引っ張られている点です。両者が中間へ寄る力が働いています。
JPDA では、実行可能事象のうち支配的な2つを比べます。事象 $\theta_A$（$\bm{z}_1\to$航跡1, $\bm{z}_2\to$航跡2）の重みは、
$$ w_A = \frac{P_D\mathcal{N}_{1,1}}{\lambda}\cdot\frac{P_D\mathcal{N}_{2,2}}{\lambda} = \frac{0.9\times0.020}{10^{-3}}\cdot\frac{0.9\times0.020}{10^{-3}} = 18\times18 = 324 $$
事象 $\theta_B$（$\bm{z}_1\to$航跡2, $\bm{z}_2\to$航跡1）は、
$$ w_B = \frac{0.9\times0.010}{10^{-3}}\cdot\frac{0.9\times0.010}{10^{-3}} = 9\times9 = 81 $$
この2つだけで正規化すると $P(\theta_A) = 324/405 = 0.80$、$P(\theta_B) = 0.20$ です。周辺化すると航跡1については $\beta_1^{(1)} = 0.80,\ \beta_2^{(1)} = 0.20$ となります。
左のパネルは実行可能な同時事象を列挙したもので、「両航跡がともに $z_1$ 由来」という組み合わせが最初から候補に存在しないことが排他制約の実体です。中央の棒グラフは同じ尤度から得られる $\beta$ を並べたもので、誤った対応 $z_2$ への重みが PDA の 0.332 から JPDA の 0.200 へ約4割減っています。右のパネルは尤度比を横軸に振ったもので、JPDA（青）の曲線が常に PDA（緑）の上にあり、とくに尤度比が 2〜4 の「判断に迷う領域」で差が最大になります。曖昧なときほど JPDA の方が決められるという、実務でもっとも効く性質がここに出ています。
PDA の $(0.664, 0.332)$ と比べて、JPDA の $(0.80, 0.20)$ は正しい対応にはっきり傾いています。理由は明快で、$\theta_B$ という「両方とも遠い観測を取る」不自然な事象は、2つの小さい尤度の積になるためダブルでペナルティを受けるからです。独立 PDA では航跡ごとに1回しかペナルティが効かないので、この差が生じます。この鋭さの違いが、併合するかしないかを分けます。
計算量という代償
美しい定式化ですが、実行可能事象の数は航跡数・観測数に対して組合せ的に爆発します。実際、JPDA の $\beta$ を厳密に計算する問題は行列パーマネントの計算と等価で、#P完全であることが知られています。数十目標が密集する現場では厳密計算は不可能です。
実用では次のような近似が使われます。
- クラスタリング — ゲートを共有する航跡だけをグループ化し、グループごとに厳密計算する。ほとんどの場面でグループは2〜4航跡程度に収まります
- Murty のアルゴリズム — 上位 $K$ 個の割当だけを列挙し、それらで近似する
- Cheap JPDA（Fitzgerald） — 尤度の行和・列和から $\beta$ を近似的に構成し、事象列挙を完全に回避する
- ベリーファンクション / MCMC — パーマネントを確率的に近似する
さらに、JPDA でも併合が完全になくなるわけではありません。JPDA* と呼ばれる変種は、$\beta$ の混合による平滑化を抑えるために支配的な仮説だけを残すことで、併合をさらに抑制します。
理論の準備が整いました。ここからは実際にコードを書いて、これらの手法が交差目標のシナリオでどう振る舞うかを定量的に確かめましょう。
Python 実装 — シナリオとフィルタ
モデルとパラメータの設定
まず状態空間モデル、レーダーのパラメータ、ゲートしきい値を定義します。監視領域は 240 m × 120 m、観測雑音の標準偏差は各軸 2 m とします。
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.stats import chi2
# 日本語フォント設定
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
T = 1.0 # スキャン周期 [s]
F = np.array([[1, T, 0, 0], [0, 1, 0, 0], [0, 0, 1, T], [0, 0, 0, 1]], float)
H = np.array([[1, 0, 0, 0], [0, 0, 1, 0]], float)
q = 0.05 # プロセス雑音強度
Qb = np.array([[T**3 / 3, T**2 / 2], [T**2 / 2, T]]) * q
Q = np.zeros((4, 4)); Q[:2, :2] = Qb; Q[2:, 2:] = Qb
sigma_m = 2.0 # 観測雑音の標準偏差 [m]
R = np.eye(2) * sigma_m**2
PD, PG = 0.9, 0.99 # 検出確率・ゲート確率
gamma = chi2.ppf(PG, 2) # ゲートしきい値（自由度2）
lam = 2e-4 # クラッタ密度 [個/m^2]
K = 40 # スキャン数
XLIM, YLIM = (-120, 120), (-60, 60)
V_REGION = (XLIM[1] - XLIM[0]) * (YLIM[1] - YLIM[0])
print(f"ゲートしきい値 gamma = {gamma:.3f}")
print(f"監視領域 {V_REGION} m^2, 1スキャンあたり期待誤警報数 = {lam * V_REGION:.2f} 個")
実行すると gamma = 9.210、期待誤警報数 5.76 個と表示されます。$P_G = 0.99$ に対応する自由度2のカイ二乗分位点が 9.21 で、先ほどの表と一致します。1スキャンあたり平均6個ほどの誤警報が監視領域に散らばる、そこそこ厳しい環境です。
真値の生成 — 交差する2目標
2機を、$x$ 方向には同じ速度で並走させ、$y$ 方向には互いに逆向きに動かして交差させます。ただし $x$ 方向に 6 m のオフセットを入れ、最接近距離が 6 m（＝観測雑音の3倍）になるようにします。「近すぎて区別が難しいが、原理的には区別可能」という、アソシエーションが試される絶妙な設定です。
def make_truth(dx=6.0):
"""2目標の真値軌跡（k=20 付近で最接近）"""
k = np.arange(1, K + 1)
t1 = np.stack([-80 + 4 * k, np.full(K, 4.0), -20 + 1.0 * k, np.full(K, 1.0)], 1)
t2 = np.stack([-80 + 4 * k + dx, np.full(K, 4.0), 20 - 1.0 * k, np.full(K, -1.0)], 1)
return t1, t2
def make_scan(rng, k, t1, t2):
"""スキャン k の観測集合（検出漏れ + 一様クラッタ）"""
Z = []
for t in (t1, t2):
if rng.random() < PD: # 検出確率 PD で検出
Z.append(t[k, [0, 2]] + rng.normal(0, sigma_m, 2))
for _ in range(rng.poisson(lam * V_REGION)): # 誤警報はポアソン個
Z.append(np.array([rng.uniform(*XLIM), rng.uniform(*YLIM)]))
return np.array(Z) if Z else np.zeros((0, 2))
t1, t2 = make_truth()
sep = np.linalg.norm(t1[:, [0, 2]] - t2[:, [0, 2]], axis=1)
print(f"目標間距離: 初期 {sep[0]:.1f} m / 最小 {sep.min():.1f} m (k={sep.argmin()+1}) / 終端 {sep[-1]:.1f} m")
出力は「初期 38.5 m / 最小 6.0 m (k=20) / 終端 40.4 m」となります。開始時は 38.5 m 離れているので楽に区別でき、$k=20$ で 6 m まで接近し、その後また離れていきます。$\sqrt{\det\bm{S}}$ が定常状態で約 6.4 なので、最接近時の分離はイノベーションの標準偏差とほぼ同じオーダーです。ここが正念場になります。
シナリオの可視化
真値、観測、誤警報を1枚の図にしてみます。
rng = np.random.default_rng(2000)
scans = [make_scan(rng, k, t1, t2) for k in range(K)]
fig, ax = plt.subplots(figsize=(11, 5.5))
allZ = np.vstack(scans)
ax.scatter(allZ[:, 0], allZ[:, 1], s=9, c="0.65", label="観測（目標+誤警報が混在）")
ax.plot(t1[:, 0], t1[:, 2], "-", lw=2.2, color="#0072B2", label="目標1の真値")
ax.plot(t2[:, 0], t2[:, 2], "-", lw=2.2, color="#D55E00", label="目標2の真値")
kx = sep.argmin()
ax.scatter(t1[kx, 0], t1[kx, 2], s=120, marker="o", facecolors="none",
edgecolors="k", lw=1.8)
ax.annotate("最接近（分離 6 m）", xy=(t1[kx, 0], t1[kx, 2]), xytext=(-20, 40),
arrowprops=dict(arrowstyle="->"), fontsize=11)
ax.set_xlabel("x 座標 [m]"); ax.set_ylabel("y 座標 [m]")
ax.set_title("交差する2目標と一様クラッタ（40スキャン分を重ね描き）")
ax.legend(loc="upper left"); ax.grid(alpha=0.3)
plt.tight_layout(); plt.show()
この図から、追尾側が置かれている状況の厳しさがわかります。灰色の点だけを見せられて「この中から2本の滑らかな軌跡を取り出せ」と言われているわけで、目標由来の点と誤警報を見た目で区別する手がかりはありません。しかも中央付近では2本の真値が交差しており、その周辺の点がどちらの目標のものかは点の位置だけからは決められません。時間方向の連続性（＝運動モデル）を使って初めて分離できる、というのが追尾フィルタの存在意義です。
ゲーティングの検証
理論で導いた「$d^2 \sim \chi^2_2$」を数値で確かめておきます。共分散 $\bm{S}$ を持つガウス乱数を大量に発生させ、マハラノビス距離の二乗のヒストグラムを理論分布と重ねます。
rng = np.random.default_rng(0)
S_demo = np.array([[9.0, 2.0], [2.0, 4.0]]) # わざと相関のある S
L = np.linalg.cholesky(S_demo)
v = (L @ rng.standard_normal((2, 200000))).T # v ~ N(0, S)
d2 = np.einsum("ij,jk,ik->i", v, np.linalg.inv(S_demo), v)
print(f"経験的 P(d2 <= 9.21) = {(d2 <= 9.21).mean():.4f} (理論 {chi2.cdf(9.21, 2):.4f})")
print(f"経験平均 {d2.mean():.3f} (理論 2.0), 経験分散 {d2.var():.3f} (理論 4.0)")
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].hist(d2, bins=120, range=(0, 20), density=True, color="0.75",
label="経験分布（数値実験）")
xs = np.linspace(0, 20, 400)
axes[0].plot(xs, chi2.pdf(xs, 2), lw=2.2, color="#D55E00", label=r"理論: $\chi^2_2$")
axes[0].axvline(gamma, color="#0072B2", ls="--", lw=2, label=f"ゲート $\\gamma$={gamma:.2f}")
axes[0].set_xlabel(r"マハラノビス距離の二乗 $d^2$"); axes[0].set_ylabel("確率密度")
axes[0].set_title("イノベーションのマハラノビス距離はカイ二乗分布に従う")
axes[0].legend(); axes[0].grid(alpha=0.3)
g = np.linspace(0.5, 20, 300)
axes[1].plot(g, chi2.cdf(g, 2), lw=2.2, color="#0072B2", label=r"$P_G$（真の観測を捕らえる確率）")
ax2 = axes[1].twinx()
ax2.plot(g, lam * np.pi * g * np.sqrt(np.linalg.det(S_demo)), lw=2.2,
color="#D55E00", ls="--", label="ゲート内の期待誤警報数")
axes[1].set_xlabel(r"ゲートしきい値 $\gamma$"); axes[1].set_ylabel(r"ゲート確率 $P_G$")
ax2.set_ylabel("ゲート内の期待誤警報数 $\\lambda V_G$")
axes[1].set_title("ゲートを広げるほど取りこぼしは減るが誤警報が増える")
axes[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
上の図は、導出した白色化の議論を3段階でそのまま可視化したものです。左は相関のある $\bm{v} \sim \mathcal{N}(\bm{0},\bm{S})$ の散布で、ゲートは傾いた楕円になります。中央は $\bm{u} = \bm{L}^{-1}\bm{v}$ を施した後で、点群が等方的になると同時に楕円ゲートがぴったり半径 $\sqrt{\gamma}$ の円に変わっているのが見えます。この「円になる」ことこそが $d^2 = \|\bm{u}\|^2$ を独立な標準正規変数の二乗和として読める根拠です。右のヒストグラムは経験分布が理論の $\chi^2_2$ 密度に完全に重なることを示しており、数値でも $P(d^2 \le 9.21)$ が経験値 0.9897 に対して理論値 0.9900、平均は 2.006（理論 2）、分散は 4.043（理論 4）と一致しました。
本文コードの右の図はゲーティング設計のトレードオフをそのまま表しています。$P_G$ の曲線は $\gamma$ が 6 を超えたあたりから寝てしまい、9 以降はほとんど平らです。一方、期待誤警報数は $\gamma$ に線形に（2次元では $V_G \propto \gamma$）増え続けます。つまり $\gamma$ を 9.21 より大きくしても、得られる $P_G$ の改善はわずか 1 % 未満なのに、拾う誤警報は際限なく増えます。$P_G = 0.99$ 前後が実務の定番になっている理由が、この図に集約されています。
4手法を1つの関数にまとめる
予測・ゲーティング・更新は共通なので、アソシエーション部分だけを分岐させる形で実装します。まず予測とゲーティングまでです。
def gate_and_likelihood(states, Z):
"""予測状態リストと観測集合から、ゲート判定・マハラノビス距離・尤度を返す"""
pred = []
for x, P in states:
S = H @ P @ H.T + R
Si = np.linalg.inv(S)
pred.append((H @ x, S, Si, P @ H.T @ Si, np.linalg.det(S)))
nT, m = len(states), len(Z)
d2 = np.full((nT, m), np.inf)
L = np.zeros((nT, m))
for i in range(nT):
zh, S, Si, W, detS = pred[i]
for j in range(m):
v = Z[j] - zh
d2[i, j] = v @ Si @ v
val = d2 <= gamma # ゲート内フラグ
for i in range(nT):
for j in range(m):
if val[i, j]:
L[i, j] = np.exp(-0.5 * d2[i, j]) / (2 * np.pi * np.sqrt(pred[i][4]))
return pred, d2, val, L
ゲート内かどうかの判定は $d^2 \le \gamma$ の一行で済み、尤度はガウス密度をそのまま計算しています。ゲート外の尤度を 0 にしておくことで、以降の $\beta$ 計算で自動的に除外されます。
続いて NN と GNN の割当部分です。GNN は2航跡なので総当たりで厳密解が求まります（一般には scipy.optimize.linear_sum_assignment を使います）。
def assign_nn(d2, val, method):
"""NN: 各航跡が独立に最近傍を取る / GNN: 排他制約つき総コスト最小"""
nT = d2.shape[0]
a = [-1] * nT
if method == "NN":
for i in range(nT):
idx = np.where(val[i])[0]
if len(idx):
a[i] = idx[np.argmin(d2[i, idx])]
return a
cand = [np.concatenate([[-1], np.where(val[i])[0]]).astype(int) for i in range(nT)]
best = np.inf
for a0 in cand[0]:
for a1 in cand[1]:
if a0 >= 0 and a0 == a1: # 同じ観測を2航跡に割り当てるのは禁止
continue
c = (d2[0, a0] if a0 >= 0 else gamma) + (d2[1, a1] if a1 >= 0 else gamma)
if c < best:
best, a = c, [a0, a1]
return a
非割当（検出漏れ）のコストを $\gamma$ に設定しているのがポイントです。$d^2 > \gamma$ の観測を取るより「取らない」方が安い、という自然な基準になります。
PDA と JPDA の $\beta$ 計算
導出した式をそのままコードにします。PDA は航跡ごとに独立、JPDA は実行可能事象を列挙して周辺化します。
def beta_pda(L, nT, m):
"""PDA: 航跡ごとに独立に関連確率を計算"""
B = np.zeros((nT, m + 1)) # 列0 が「非検出」に対応
for i in range(nT):
w = np.concatenate([[lam * (1 - PD * PG)], PD * L[i]])
B[i] = w / w.sum()
return B
def beta_jpda(L, val, nT, m):
"""JPDA: 実行可能な同時事象を列挙し、周辺化して関連確率を得る"""
cand = [np.concatenate([[-1], np.where(val[i])[0]]).astype(int) for i in range(nT)]
events, weights = [], []
for a0 in cand[0]:
for a1 in cand[1]:
if a0 >= 0 and a0 == a1: # 排他制約：同一観測の共有は不可
continue
p = 1.0
for i, a in enumerate((a0, a1)):
p *= (1 - PD * PG) if a < 0 else PD * L[i, a] / lam
events.append((a0, a1)); weights.append(p)
weights = np.array(weights); weights /= weights.sum()
B = np.zeros((nT, m + 1))
for e, (a0, a1) in enumerate(events):
B[0, a0 + 1] += weights[e] # 周辺化：該当事象の確率を足し込む
B[1, a1 + 1] += weights[e]
return B
beta_jpda の二重ループが同時事象の列挙で、if a0 >= 0 and a0 == a1: continue の一行が排他制約そのものです。この一行があるかないかで、航跡が併合するかしないかが決まります。列挙後は事象確率を正規化し、各航跡について「その観測を選んでいる事象」の確率を足し込むだけで周辺化が完了します。
更新の実装
PDA / JPDA の更新は、$\beta$ の出所によらず共通です。導出した $\bar{\bm{v}}$ と $\tilde{\bm{P}}$ をそのまま書きます。
def update_pda(states, pred, Z, B):
"""重み付きイノベーションによる更新 + イノベーション散らばり項"""
out = []
for i, (x, P) in enumerate(states):
zh, S, Si, W, detS = pred[i]
b = B[i]
vbar = np.zeros(2)
Spread = np.zeros((2, 2))
for j in range(len(Z)):
v = Z[j] - zh
vbar += b[j + 1] * v # 重み付き平均イノベーション
Spread += b[j + 1] * np.outer(v, v) # 二乗の重み付き平均
Spread -= np.outer(vbar, vbar) # 平均の二乗を引く
Pc = P - W @ S @ W.T # 通常のカルマン更新共分散
Pnew = b[0] * P + (1 - b[0]) * Pc + W @ Spread @ W.T
out.append([x + W @ vbar, Pnew])
return out
Spread の2行が「イノベーションの散らばり」項 $\tilde{\bm{P}}$ で、まさに「二乗の平均 − 平均の二乗」という分散の形になっています。$\beta$ が1つの観測に集中していればここはほぼゼロ行列になり、複数に分散していれば大きくなって共分散を膨らませます。この項を落とすとフィルタは過信して発散するので、実装上の要注意点です。
追尾ループの全体
上の部品を組み合わせて、1回のシミュレーションを走らせる関数にします。
def track(method, rng, t1, t2):
"""method: 'NN' / 'GNN' / 'PDA' / 'JPDA'。誤差と推定軌跡を返す"""
P0 = np.diag([sigma_m**2, 4.0, sigma_m**2, 4.0])
states = []
for t in (t1, t2): # 真値近傍から初期化（追尾開始済みを仮定）
x = t[0].astype(float).copy()
x[[0, 2]] += rng.normal(0, sigma_m, 2)
states.append([x, P0.copy()])
err = np.zeros((K, 2)); hist = np.zeros((K, 2, 2))
for k in range(K):
Z = make_scan(rng, k, t1, t2)
states = [[F @ x, F @ P @ F.T + Q] for x, P in states] # 予測
pred, d2, val, L = gate_and_likelihood(states, Z) # ゲーティング
m = len(Z)
if method in ("NN", "GNN"):
a = assign_nn(d2, val, method)
for i in range(2):
if a[i] < 0:
continue
zh, S, Si, W, detS = pred[i]
x, P = states[i]
states[i] = [x + W @ (Z[a[i]] - zh), P - W @ S @ W.T]
else:
B = beta_pda(L, 2, m) if method == "PDA" else beta_jpda(L, val, 2, m)
states = update_pda(states, pred, Z, B)
for i in range(2):
hist[k, i] = states[i][0][[0, 2]]
err[k, i] = np.linalg.norm(hist[k, i] - (t1, t2)[i][k, [0, 2]])
return err, hist
初期化は真値の近傍から行い、「追尾はすでに確立している」状況を仮定しています（航跡開始のロジックは本記事の範囲外です）。予測は全手法共通、ゲーティングも共通で、分岐しているのはアソシエーションと更新の部分だけです。この構造は理論の整理そのままなので、他の手法（MHT など）に差し替えるのも同じ場所です。
ここまでで道具が揃いました。実際に走らせて、何が起きるかを見ましょう。
実験 — 交差で何が起きるか
単一ランの軌跡比較
まず1回分の追尾を可視化します。seed=2000 は NN が失敗し JPDA が持ちこたえる典型的なランです。
fig, axes = plt.subplots(2, 2, figsize=(13, 8), sharex=True, sharey=True)
titles = {"NN": "最近傍法（NN）", "GNN": "大域最近傍法（GNN）",
"PDA": "確率的データ関連（PDA）", "JPDA": "同時確率的データ関連（JPDA）"}
for ax, meth in zip(axes.ravel(), ["NN", "GNN", "PDA", "JPDA"]):
e, h = track(meth, np.random.default_rng(2000), t1, t2)
ax.plot(t1[:, 0], t1[:, 2], "-", color="0.6", lw=3, label="真値")
ax.plot(t2[:, 0], t2[:, 2], "-", color="0.6", lw=3)
ax.plot(h[:, 0, 0], h[:, 0, 1], "o-", ms=3.2, color="#0072B2", label="航跡1の推定")
ax.plot(h[:, 1, 0], h[:, 1, 1], "s-", ms=3.2, color="#D55E00", label="航跡2の推定")
ax.set_title(f"{titles[meth]} 終端誤差 {e[-10:].mean():.1f} m")
ax.grid(alpha=0.3)
if meth == "NN":
ax.legend(loc="upper left", fontsize=9)
for ax in axes[1]:
ax.set_xlabel("x 座標 [m]")
for ax in axes[:, 0]:
ax.set_ylabel("y 座標 [m]")
plt.tight_layout(); plt.show()
このランでは NN の終端誤差が 17.1 m、GNN が 1.5 m、PDA が 1.6 m、JPDA が 1.5 m でした。図を見ると違いは一目瞭然です。NN では交差点を過ぎたあと2本の推定軌跡が同じ方向へ折れ曲がってしまい、片方の目標を2本の航跡が追いかける状態（併合）に陥っています。真値が上下に分かれていくのに対して、推定は分かれられていません。一方 JPDA では、交差点でいったん推定が真値からずれるものの、数スキャンで正しく分離して各目標に復帰しています。排他制約による「反発」が効いた結果です。
モンテカルロによる定量比較
1ランでは偶然の可能性があるので、乱数を変えて500回繰り返し、統計を取ります。
def evaluate(method, nmc=500, seed0=2000):
E = np.zeros((nmc, K, 2)); Hh = np.zeros((nmc, K, 2, 2))
for mc in range(nmc):
E[mc], Hh[mc] = track(method, np.random.default_rng(seed0 + mc), t1, t2)
# 入れ替わり: 終盤10スキャンで航跡1が目標2の方に近いか
d_self = np.linalg.norm(Hh[:, -10:, 0] - t1[None, -10:, [0, 2]], axis=-1).mean(1)
d_cross = np.linalg.norm(Hh[:, -10:, 0] - t2[None, -10:, [0, 2]], axis=-1).mean(1)
swap = (d_cross < d_self).mean()
# 併合: 終盤10スキャンで2航跡の推定間距離が真の距離の半分未満
est_sep = np.linalg.norm(Hh[:, -10:, 0] - Hh[:, -10:, 1], axis=-1).mean(1)
coal = (est_sep < 0.5 * sep[-10:].mean()).mean()
return E, swap, coal
results = {}
for meth in ["NN", "GNN", "PDA", "JPDA"]:
E, swap, coal = evaluate(meth)
results[meth] = (np.sqrt((E**2).mean(axis=(0, 2))), swap, coal)
print(f"{meth:5s} 全RMSE {np.sqrt((E**2).mean()):5.2f} m | "
f"交差前 {np.sqrt((E[:, :15]**2).mean()):5.2f} | "
f"交差後 {np.sqrt((E[:, 25:]**2).mean()):5.2f} | "
f"入替率 {100*swap:4.1f}% | 併合率 {100*coal:4.1f}%")
結果は次のとおりでした。
| 手法 | 全RMSE [m] | 交差前 ($k\le15$) | 交差後 ($k\ge26$) | 入れ替わり率 | 併合率 |
|---|---|---|---|---|---|
| NN | 13.22 | 4.79 | 19.27 | 14.0 % | 26.6 % |
| GNN | 11.36 | 4.66 | 16.08 | 2.0 % | 0.8 % |
| PDA | 9.77 | 3.53 | 14.60 | 14.2 % | 30.0 % |
| JPDA | 5.90 | 3.09 | 8.16 | 0.6 % | 0.2 % |
3つのパネルは同じ500ランを別の角度から見たものです。左は交差前後の RMSE で、灰色の棒（交差前）はどれも 3〜5 m とほぼ横並びなのに、色付きの棒（交差後）は 8.16〜19.27 m と大きく開きます。中央は入れ替わり率と併合率で、RMSE が NN より小さい PDA の併合率が 30.0 % と最悪という逆転が見えます。右は終盤の2航跡の推定間距離のヒストグラムで、真の距離 31.6 m の位置に山があるのが正常、0 m 付近の山が併合したランです。NN と PDA には 0 m 付近に明確な山があり、GNN と JPDA にはほとんどありません。
この表と図からいくつもの事実が読み取れます。第一に、交差前の性能差は小さい（4.79 対 3.09 m）のに、交差後は3倍近い差（19.27 対 8.16 m）に開きます。目標が離れている限り手法の違いはほとんど問題にならず、接近した瞬間にすべてが決まるという、多目標追尾の本質がそのまま数字に出ています。
第二に、NN と PDA は入れ替わり・併合の両方を高い率で起こしています。 NN は排他制約がないので2航跡が同じ観測を取り、PDA は隣の目標をクラッタとみなして両方に重みを与えます。原因は違いますが、症状（併合率 26.6 % と 30.0 %）は同じです。実際、終盤の2航跡の推定間距離を調べると、真の距離が平均 31.6 m あるのに対し、PDA では 30 % のランで 10 m 未満に潰れていました。JPDA では 0 % です。
第三に、PDA は NN より RMSE が小さい（9.77 対 13.22）にもかかわらず、併合率はむしろ高い。 一見矛盾しますが、これは PDA が「間違えるときも中間で間違える」ためです。ハードに決める NN は間違えると大きく外れ、確率で重み付ける PDA は誤りが平均化されて誤差の絶対値は小さくなります。RMSE だけを見ていては併合を見逃すという、追尾の評価における重要な教訓です。
第四に、GNN は NN より明確に良い（11.36 対 13.22、併合率 26.6 %→0.8 %）。 排他制約を入れるだけでこれだけ改善します。しかし RMSE では PDA に負けます。ハードな決定が持つ「間違いを丸ごと取り込む」性質は、排他制約では救えないからです。排他性（GNN の強み）と柔らかい重み付け（PDA の強み）の両方を持つのが JPDA であり、実際 JPDA だけがすべての指標で最良になっています。
誤差の時間推移
RMSE の時間変化を見ると、何が起きているかがさらにはっきりします。
fig, ax = plt.subplots(figsize=(11, 5))
colors = {"NN": "#D55E00", "GNN": "#E69F00", "PDA": "#009E73", "JPDA": "#0072B2"}
for meth in ["NN", "GNN", "PDA", "JPDA"]:
ax.plot(np.arange(1, K + 1), results[meth][0], lw=2.2,
color=colors[meth], label=titles[meth])
ax.axvline(sep.argmin() + 1, color="k", ls="--", lw=1.5)
ax.annotate("最接近スキャン", xy=(sep.argmin() + 1, ax.get_ylim()[1] * 0.9),
xytext=(6, 0), textcoords="offset points", fontsize=11)
ax2 = ax.twinx()
ax2.plot(np.arange(1, K + 1), sep, color="0.6", ls=":", lw=2)
ax2.set_ylabel("真の目標間距離 [m]", color="0.4")
ax.set_xlabel("スキャン番号 $k$"); ax.set_ylabel("位置 RMSE [m]（500回平均）")
ax.set_title("交差時に誤差が急増し、NN系は回復しない")
ax.legend(loc="upper left"); ax.grid(alpha=0.3)
plt.tight_layout(); plt.show()
このグラフには、素朴な予想を裏切る事実が写っています。$k = 10$ あたりまではすべての手法が 3〜5 m で横並びですが、その後 RMSE は最後まで単調に増え続け、最接近の $k = 20$ にピークらしいピークは現れません。$k = 20$ の時点で NN 9.88 m / JPDA 4.79 m だったものが、$k = 40$ では NN 24.38 m / JPDA 9.77 m へと、目標が離れていく区間でこそ差が開いていきます。JPDA も例外ではなく、増え方がもっとも緩やかというだけで、絶対値は増加します。
種明かしが右のパネルです。位置誤差の中央値を取ると、どの手法も全区間で 1.4〜2.0 m に収まっています。 つまり大半のランは最接近を越えても正しく追い続けており、左の右上がりの曲線は少数の失敗ランが二乗平均を支配した結果にすぎません。交差後（$k \ge 26$）に両航跡を 10 m 以内で保てなかったランの割合は NN 33.0 % / GNN 5.0 % / PDA 33.4 % / JPDA 2.0 %。そして失敗しなかったランだけで RMSE を取り直すと、4手法とも 1.70〜1.82 m とほぼ差がありません。
ここから読み取るべきは、手法の違いは「うまくいったときの精度」ではなく「失敗する確率」に現れるということです。アソシエーションの巧拙は推定を少しずつ改善するのではなく、追尾が破綻するかしないかという離散的な分岐を左右します。
では失敗ランの誤差はなぜ時間とともに増えるのか。併合が起きると、2本の航跡は中間の幻に取り残されるか、どちらか一方の目標に2本とも吸着します。いずれの場合も「誰にも追われていない目標」が1つ生まれ、その目標に対する誤差は目標間距離そのものになります。目標が離れるほど距離は開くので、誤差は機械的に増大するわけです。実際、上で見た seed=2000 の NN のランでは、終盤10スキャンで2本の航跡の間隔が平均 0.74 m（最終スキャンでは 0.35 m）まで潰れて両方とも目標2を追っており、目標1に対する誤差は 32.5 m に達していました。「一度失った航跡は自然には戻らない」という追尾の鉄則が、右上がりの曲線として可視化されています。
なお、右のパネルで NN と PDA の中央値が $k = 20$ 前後にだけ小さな盛り上がり（1.9〜2.2 m）を見せているのは、最接近の数スキャンでは失敗していないランでも推定が隣の目標に引かれるためです。この一時的な劣化はすぐ収まりますが、その数スキャンの間に判断を誤ったランだけが二度と戻れなくなる — というのが、左の曲線の裏側で起きていることです。
クラッタ密度を変えたときの頑健性
最後に、交差ではなく「クラッタの濃さ」に対する頑健性を調べます。$\lambda$ を 5 倍・20 倍に上げて RMSE を比較します（200回平均）。
rows = []
for l in [2e-4, 1e-3, 4e-3]:
globals()["lam"] = l # クラッタ密度を差し替え
row = {}
for meth in ["NN", "GNN", "PDA", "JPDA"]:
E = np.zeros((200, K, 2))
for mc in range(200):
E[mc], _ = track(meth, np.random.default_rng(3000 + mc), t1, t2)
row[meth] = np.sqrt((E**2).mean())
rows.append((l, row))
print(f"lam={l:.0e} (1スキャン {l*V_REGION:5.1f}個) " +
" ".join(f"{m} {row[m]:5.2f}" for m in row))
globals()["lam"] = 2e-4 # 元に戻す
fig, ax = plt.subplots(figsize=(9, 5))
lams = [r[0] for r in rows]
for meth in ["NN", "GNN", "PDA", "JPDA"]:
ax.plot([l * V_REGION for l in lams], [r[1][meth] for r in rows],
"o-", lw=2.2, ms=7, color=colors[meth], label=titles[meth])
ax.set_xscale("log")
ax.set_xlabel("1スキャンあたりの誤警報数 [個]"); ax.set_ylabel("位置 RMSE [m]（200回平均）")
ax.set_title("クラッタが濃くなるとNN系は破綻し、PDA系は耐える")
ax.legend(); ax.grid(alpha=0.3, which="both")
plt.tight_layout(); plt.show()
図を見ると、NN と GNN（橙・黄）の2本がぴたりと重なって急勾配で上昇し、PDA と JPDA（緑・青）の2本が緩やかな傾きで下に留まるという、手法が2つのグループにきれいに分かれる構図になっています。分かれ目は排他制約の有無ではなく「ハードに1つ選ぶか、確率で重み付けるか」です。
結果は劇的でした。1スキャンあたりの誤警報数が 5.8 個のとき NN 12.98 m / PDA 10.24 m と大差ありませんが、28.8 個になると NN 30.05 m に対し PDA 12.55 m、115.2 個では NN 49.17 m に対し PDA 16.30 m と、3倍の開きになります。GNN も NN とほぼ同じように破綻しており（28.62 m、49.19 m）、排他制約はクラッタ耐性には寄与しないことがわかります。
理由は $\beta_i$ の式にあります。PDA では分母の $\lambda(1-P_DP_G)$ がクラッタ密度に比例して増えるので、$\lambda$ が上がると $\beta_0$（どれも目標ではない）が自動的に大きくなり、フィルタは観測を信用しなくなって予測を維持します。濃いクラッタの中では「動かない」ことが最善の戦略であり、PDA はその判断を数式から自動的に導いています。対する NN は、ゲート内にクラッタしかなくても必ずどれかを選び、その1点を全面的に信じて更新してしまいます。ここが両者を分ける決定的な設計思想の差です。
なお、クラッタ密度が高い領域では PDA と JPDA の差は小さくなります（16.30 対 15.90 m）。クラッタが支配的な状況では失敗の主要因が「目標同士の混同」から「クラッタへの誤同定」に移るため、排他制約の恩恵が相対的に薄れるからです。JPDA が真価を発揮するのは、クラッタは中程度で目標が密集する状況だと言えます。
こうして、交差シナリオとクラッタ環境の両方で、理論から予想された振る舞いがそのまま数値に現れることを確認できました。最後に、ここまでの議論を要点として整理しておきましょう。
まとめ
本記事では、レーダーの Track-While-Scan における中心的な課題であるデータアソシエーションを、理論の導出から Python 実装まで通して解説しました。
- イノベーション共分散 $\bm{S} = \bm{H}\bm{P}\bm{H}^\top + \bm{R}$ は「予測の怪しさ」と「センサの雑さ」の和である。交差項が消えるのは、観測雑音が過去情報と独立だから
- ゲーティングは、コレスキー分解による白色化 $\bm{u} = \bm{L}^{-1}\bm{v}$ を経て $d^2 = \bm{v}^\top\bm{S}^{-1}\bm{v} \sim \chi^2_{n_z}$ を示すことで、確率的根拠を持ったしきい値 $\gamma$ を与える。$P_G = 0.99$ なら $n_z=2$ で $\gamma = 9.21$
- ゲート体積 $V_G = c_{n_z}\gamma^{n_z/2}\sqrt{\det\bm{S}}$ とゲート内期待誤警報数 $\lambda V_G$ の間にトレードオフがあり、$\gamma$ を広げすぎても $P_G$ の改善はわずかで誤警報だけが増える
- NN / GNN はマハラノビス距離最小＝最尤選択という正当化を持つが、ハードな決定ゆえに誤ると回復できない。GNN の排他制約は併合を防ぐが（併合率 26.6 %→0.8 %）、誤りの丸ごと取り込みは防げない
- PDA は関連確率 $\beta_i$ で重み付け、$\bar{\bm{v}} = \sum\beta_i\bm{v}_i$ を1つの観測とみなす普通のカルマン更新に帰着する。$\lambda$ が分母に現れるため、クラッタが濃いほど自動的に慎重になる
- イノベーションの散らばり項 $\tilde{\bm{P}} = \bm{W}[\sum\beta_i\bm{v}_i\bm{v}_i^\top – \bar{\bm{v}}\bar{\bm{v}}^\top]\bm{W}^\top$ が、関連の不確かさを共分散に正直に記録する。これを落とすとフィルタは過信して発散する
- PDA は単一目標前提なので、近接目標では両航跡が中間に引かれて併合する（実験で 30.0 %）。JPDA は実行可能な同時事象を列挙して排他制約を課し、周辺化で $\beta$ を得ることでこれを解消する（0.2 %）
- 実験では、交差前の RMSE 差はわずか（4.79 対 3.09 m）でも交差後は 19.27 対 8.16 m まで開いた。RMSE だけを見ると併合を見逃すので、入れ替わり率・併合率も併せて評価する必要がある
- 誤差の中央値はどの手法も全区間で 1.4〜2.0 m にとどまり、交差後に追尾を失うランの割合だけが NN 33.0 % / GNN 5.0 % / PDA 33.4 % / JPDA 2.0 % と大きく違う。アソシエーションの巧拙は「うまくいったときの精度」ではなく「破綻する確率」に現れる
データアソシエーションは、追尾の「入口」であると同時に性能の上限を決める部分です。ここを押さえたうえで、非線形観測モデル（距離・方位で観測する現実のレーダー）、機動する目標への対応、そして仮説を複数スキャンにわたって保持する多仮説追尾（MHT）へと進むと、実用的な追尾システムの全体像が見えてきます。
次のステップとして、以下の記事も参考にしてください。
- 拡張カルマンフィルタ — 距離・方位で観測する現実のレーダーモデルへの拡張
- アンセンテッドカルマンフィルタ — 強い非線形性に対する代替アプローチ
- パーティクルフィルタ — 非ガウス・多峰性の事後分布を扱う
- カルマンフィルタと逐次ベイズ推定 — アソシエーションの周辺化を統一的に見る視点
- ポアソン過程 — クラッタモデルの数学的基礎
- 多変量正規分布 — 尤度計算とゲーティングの土台
- モノパルス追尾レーダー — 単一目標追尾における角度測定の高精度化