1993・粒子濾波：Gordon–Salmond–Smith 的自助濾波器#

「Kalman 用兩個數字追蹤世界，粒子用一群螞蟻淹沒世界。」本案時間：1993 年，地點：英國皇家信號與雷達研究院，報案人：雷達追蹤官。

案件檔案	內容
案號	PF-1993-GSS
年份	1993 年
主角	Neil Gordon、David Salmond、Adrian Smith
關鍵文獻	Novel Approach to Nonlinear/Non-Gaussian Bayesian State Estimation（1993）
涉案對象	狀態 \(x\) 、觀測 \(y\) 、粒子 \(N\) 、權重 \(w\)
核心證物	bootstrap 濾波器、SIR 重採樣、雷達追蹤實測
關聯舊案	1960 年 Kalman 濾波、1992 年數值聖經、1998 年粗糙路徑
遺產	SLAM、基因定序、資料同化的蒙地卡羅正統

案發現場#

1993 年的雷達螢幕前，追蹤官盯著一個亂跳的光點。

目標飛機忽快忽慢，還會機動轉彎，量測夾雜雜波與漏檢。1960 年的 Kalman 濾波優雅無比，前提卻是線性高斯：一旦運動非線性、噪聲非高斯，兩個數字（均值與協方差）就裝不下真實的後驗分佈。多假設追蹤（MHT）能硬撐，計算量卻隨假設指數爆炸。

當時的貝葉斯遞迴公式人人都會寫：預測用 Chapman–Kolmogorov 積分，更新用 Bayes 定理相乘。但寫得出不代表算得出，高維積分卡死了所有人。擴展 Kalman（EKF）用線性化矇混，無味 Kalman（UKF）用 sigma 點取巧，遇到強非線性照樣翻車。

Gordon 等人是做雷達出身的實戰派。他們的想法近乎無賴：既然解析分佈裝不下，那就養一群粒子 \(\{x^i\}\) ，每顆粒子帶一個權重 \(w^i\) ，用經驗分佈硬逼近後驗。名字取得很挑釁——bootstrap 濾波器，自己提著自己的鞋帶升天。

偵查過程#

偵探先把狀態空間模型釘在牆上。設隱狀態為 \(x_t\) ，觀測為 \(y_t\) ：

\[ x_t \sim p(x_t \mid x_{t-1}) \]

\[ y_t \sim p(y_t \mid x_t) \]

目標是遞迴估計濾波分佈 \(p(x_t \mid y_{1:t})\) 。預測步把上一代粒子向前傳播，更新步按似然加權：

\[ w^i_t \propto w^i_{t-1}\,p(y_t \mid x^i_t) \]

權重歸一後，加權粒子集即為後驗的蒙地卡羅近似。期望估計化為加權平均：

\[ \mathbb{E}[g(x_t) \mid y_{1:t}] \approx \sum_{i=1}^N w^i_t\,g(x^i_t) \]

麻煩很快現形：跑幾十步之後，權重退化到只剩一顆粒子獨活，其餘皆死。這叫權重簡併，有效樣本量 \(N_{\text{eff}}\) 崩盤：

\[ N_{\text{eff}} = \frac{1}{\sum_i (w^i_t)^2} \]

解藥是重採樣（resample）：按權重有放回地複製粒子，殺死輕權重，大量繁殖重權重，之後把權重置平為 \(1/N\) 。SIR（序貫重要性重採樣）三部曲就此定型：sample、importance-weight、resample。

下表是審訊室裡的流程筆錄，每一步都有明確口供：

步驟	動作	數學	目的
傳播	由 \(x^i_{t-1}\) 抽 \(x^i_t\)	\(p(x_t \mid x_{t-1})\)	實現預測積分
加權	算似然乘舊權	\(w \propto w\,p(y_t \mid x_t)\)	實現 Bayes 更新
診斷	算 \(N_{\text{eff}}\)	上式倒數平方和	判斷是否退化
重採樣	複製與淘汰	多項式／系統／分層重抽	重置權重，延命粒子

1993 年論文的殺手鐧，是拿雷達機動目標開刀。狀態含位置與速度，轉彎率突變，量測為帶噪的距離與方位。bootstrap 濾波器的追蹤誤差顯著低於 EKF，且不需線性化雅可比。Kalman 的線性王座，從此多了一位非線性接班人。

理論側寫隨後跟上：在適當條件下，粒子估計以 \(1/\sqrt{N}\) 速率收斂，偏差隨 \(N\) 消失。代價是維度詛咒：狀態維度一高，所需粒子數指數膨脹。這一點，至今仍是資料同化與 SLAM 的心頭大患。

結案報告#

1993 年一案，把貝葉斯濾波從解析牢籠裡放了出來。

只要能模擬狀態轉移、能算觀測似然，粒子濾波就能開工。不問線性非線性，不問高斯非高斯。這種「即插即用」的野蠻，正是它橫掃機器人定位、語音增強、流行病追蹤的原因。同一家族後來長出輔助粒子濾波、Rao–Blackwellized 粒子濾波、粒子 MCMC，族譜極盛。

在歷史年表裡，本案是 1960 年 Kalman 的叛逆之子：父親用解析封閉傳家，兒子用蒙地卡羅闖蕩。兩者在線性高斯案裡握手言和——粒子數趨於無窮時，粒子均值回到 Kalman 解；實務上更有混合打法，線性部分 Kalman 化，非線性部分粒子化。

它也反哺了數值學：每顆粒子的傳播都是一次 SDE 離散，腳下踩的是 1992 年聖經的格式。沒有穩定的傳播格式，再多的重採樣也救不了發散的粒子群。

案件狀態：已破案，雷達光點歸位，粒子群就地解散又隨時集結。

證據與工具#

證物一，bootstrap 濾波器極簡虛擬碼（粒子數 \(N\) ）：

\[ \text{init: } x^i_0 \sim p(x_0), \quad w^i_0 = 1/N \]

\[ \text{loop: propagate} \to \text{weight by } p(y_t \mid x^i_t) \to \text{normalize} \to \text{resample if } N_{\text{eff}} < N/2 \]

證物二，重採樣方案對照表：

方案	做法	特點
多項式	按 \(w\) 有放回抽 \(N\) 次	最直觀，變異數最大
系統	一把等距梳子一次定位	變異數小，實務最愛
分層	每層獨立均勻抽樣	兼顧簡單與穩定
殘差	先確定複製整數份再抽餘數	大權重粒子不浪費

證物三，驗案工具。線性高斯案以 Kalman 為真值，粒子均值應逼近 Kalman 均值；對應程式 _code/1993-particle_filter.py 示範 RMSE 達到 Kalman 九成水準。非線性案則檢查 \(N_{\text{eff}}\) 曲線：長期貼地，表示建議分佈太差，該換輔助或橋接格式。

探員守則：權重退化是病，重採樣是藥，樣本枯竭是藥的副作用。連續重採樣而不加入抖動（jitter），粒子會近親繁殖到全員同臉。

補充：程式實作#

對應程式#

1993-particle_filter.py：以一維非線性二次觀測模型實作 bootstrap 粒子濾波，並以線性子問題的 Kalman 解為基準對照。

理論呼應#

本文遞迴目標為濾波分佈 \(p(x_t \mid y_{1:t})\) ，更新步按 \(w^i_t \propto w^i_{t-1}\,p(y_t \mid x^i_t)\) 加權，期望以 \(\sum_i w^i_t\,g(x^i_t)\) 逼近。退化診斷用有效樣本量 \(N_{\mathrm{eff}} = 1/\sum_i (w^i_t)^2\) ，程式在每步後以系統重採樣把權重置平，正是 SIR 三部曲的實作。理論保證粒子估計以 \(1/\sqrt{N}\) 速率收斂，程式以 RMSE 比值驗收：在二次觀測的非線性案中追上線性 Kalman 基準，即為收斂的具體兌現。

執行方式#

在 隨機微積分 目錄下執行 python3 _code/1993-particle_filter.py 。

實測輸出#

RMSE_pf = 0.484537, RMSE_kal = 0.562554, ratio = 0.861316
VERIFY rmse_pf=0.484537 rmse_kal=0.562554 ratio=0.861316

讀者實驗#

將粒子數 Np 由 2000 改為 500 再重跑，觀察 RMSE 比值是否變差並體會 \(1/\sqrt{N}\) 收斂速率的實際含義。

完整程式如下：

"""1993 粒子濾波（一維非線性、二次觀測）vs Kalman 基準。
對應 wiki：Gordon–Salmond–Smith (1993) bootstrap filter；Kalman (1960) 線性最優濾波。
狀態：x_{k+1} = a x_k + b + √Q ε（a=0.9, b=0.5, Q=0.5），均值 5 恆正，
  使二次觀測 y_k = x_k^2/10 + √R η（R=0.5）免於符號不可辨識。
粒子濾波：bootstrap（提議=轉移先驗），2000 粒子，系統重採樣。
同場 Kalman：相同動態、相同真值軌跡，但線性觀測 y_lin = x + √R η'，
  Kalman 為該線性子問題的最優解，作為 RMSE 基準。
驗證：RMSE_pf / RMSE_kal < 1.5（九成水準內、容忍 50% 退化），
  或 RMSE_pf < 2.0。閾值 2.0 的合理性：狀態平穩標準差
  √(Q/(1-a²)) ≈ 1.62，2.0 約 1.2 倍先驗 spread，代表濾波有效收斂。
只用 numpy，固定種子，不畫圖。
"""
import numpy as np

rng = np.random.default_rng(2)
a, b = 0.9, 0.5
Q, R = 0.5, 0.5
sQ, sR = np.sqrt(Q), np.sqrt(R)
T = 100
Np = 2000


def systematic_resample(w, rng):
    N = w.size
    u0 = rng.random() / N
    u = u0 + np.arange(N) / N
    cum = np.cumsum(w)
    idx = np.searchsorted(cum, u)
    return idx


# 真值與兩種觀測
x_true = np.zeros(T + 1)
x_true[0] = 5.0
for k in range(T):
    x_true[k + 1] = a * x_true[k] + b + sQ * rng.standard_normal()
y_nl = x_true[1:] ** 2 / 10.0 + sR * rng.standard_normal(T)
y_lin = x_true[1:] + sR * rng.standard_normal(T)

# --- Kalman（線性子問題基準） ---
xk, Pk = 5.0, 1.0
xhat_kal = np.zeros(T)
for k in range(T):
    xp = a * xk + b
    Pp = a * a * Pk + Q
    Kg = Pp / (Pp + R)
    xk = xp + Kg * (y_lin[k] - xp)
    Pk = (1.0 - Kg) * Pp
    xhat_kal[k] = xk
rmse_kal = float(np.sqrt(np.mean((xhat_kal - x_true[1:]) ** 2)))

# --- Bootstrap 粒子濾波（二次觀測） ---
parts = 5.0 + 1.0 * rng.standard_normal(Np)
xhat_pf = np.zeros(T)
for k in range(T):
    parts = a * parts + b + sQ * rng.standard_normal(Np)
    pred = parts ** 2 / 10.0
    # 對數權重避免下溢
    logw = -0.5 * ((y_nl[k] - pred) ** 2) / R
    logw -= np.max(logw)
    w = np.exp(logw)
    w /= np.sum(w)
    xhat_pf[k] = float(np.sum(parts * w))
    parts = parts[systematic_resample(w, rng)]

rmse_pf = float(np.sqrt(np.mean((xhat_pf - x_true[1:]) ** 2)))
ratio = rmse_pf / rmse_kal if rmse_kal > 0 else np.inf

print(f"RMSE_pf = {rmse_pf:.6f}, RMSE_kal = {rmse_kal:.6f}, ratio = {ratio:.6f}")
print(f"VERIFY rmse_pf={rmse_pf:.6f} rmse_kal={rmse_kal:.6f} ratio={ratio:.6f}")
assert np.isfinite(rmse_pf) and np.isfinite(rmse_kal)
assert (ratio < 1.5) or (rmse_pf < 2.0), f"粒子濾波未達標: ratio={ratio}, rmse={rmse_pf}"