Unscented Kalman Filter（UKF）：非線形状態推定の理論とPython実装

はじめに

カルマンフィルタは線形ガウスシステムで最適ですが、非線形システムでは直接適用できません。**拡張カルマンフィルタ（EKF）**はヤコビアンによる線形化で対処しますが、高次の非線形性では精度が低下します。

Unscented Kalman Filter（UKF）は、確率分布を代表する少数のシグマポイントを非線形関数に通すことで、ヤコビアンの計算なしに高精度な推定を実現します。

Unscented 変換

シグマポイントの生成

\(n\) 次元の状態 \(\mathbf{x}\) （平均 \(\bar{\mathbf{x}}\) , 共分散 \(\mathbf{P}\) ）に対し、\(2n+1\) 個のシグマポイントを生成します。

ここで \(\lambda = \alpha^2(n + \kappa) - n\) はスケーリングパラメータです。

重み

典型的なパラメータ: \(\alpha = 10^{-3}\) , \(\beta = 2\) （ガウス分布で最適）, \(\kappa = 0\)

重みの導出：モーメントマッチング条件

式(4)(5)の重みがなぜこの形になるのか、 Unscented変換の基礎 で扱った枠組みに沿って、モーメントマッチングの条件から導出します。

平均マッチング条件. シグマ点は中心 \(\boldsymbol{\chi}_0\) の周りに \(\pm\) 対称に配置されているため、奇数次のずれ成分は加重平均をとると打ち消し合います。したがって平均が正しく \(\bar{\mathbf{x}}\) に一致するための条件は、重みの総和が1になることだけです。式(4)(5)の重みを代入すると：

任意の \(\lambda\) （つまり任意の \(\alpha, \kappa\) の組）に対してこの恒等式が成り立ちます。つまりこの制約だけからは \(\lambda\) は一意に決まらず、\(\lambda\) はシグマ点の広がり \(\gamma = \sqrt{n+\lambda}\) をスケーリングする自由パラメータとして導入されていることが分かります。実際に \(\alpha = 10^{-3}, \kappa = 0\) で \(n = 2, 3, 4\) について数値的に確認すると、\(\sum_i W_i^{(m)} = 1.0000000000\) （浮動小数点誤差の範囲内）が常に成立します。

共分散マッチング条件と負の中心重みの意味. 次に、外側の \(2n\) 個のシグマ点だけで共分散がどこまで再現されるかを見ます。\(\mathbf{S}_i\) を \(\sqrt{\mathbf{P}}\) の第 \(i\) 列とすると、\(\boldsymbol{\chi}_i - \bar{\mathbf{x}} = \gamma \mathbf{S}_i\) （\(i \le n\) ）、\(\boldsymbol{\chi}_{i+n} - \bar{\mathbf{x}} = -\gamma \mathbf{S}_i\) なので、外側の点だけの加重共分散和は：

（最後の等号は \(\gamma^2 = n + \lambda\) と \(\sum_i \mathbf{S}_i \mathbf{S}_i^T = \mathbf{P}\) より）。重要なのは、この結果が \(W_0^{(c)}\) の値に一切依存しないことです。中心点 \(\boldsymbol{\chi}_0\) のずれは常にゼロなので、線形（あるいは恒等）変換であれば \(W_0^{(c)}\) がどんな値（負であっても）でも共分散は厳密に再現されます。これは後述のエッジケース①で数値的にも確認します。

4次モーメント（尖度）補正としての \(\beta\) . では \(W_0^{(c)}\) に加わる \((1 - \alpha^2 + \beta)\) の項は何をしているのか。これは非線形変換のもとで、事前分布の4次中心モーメント（尖度）を正しく反映させるための補正です。最も単純にこれが見える例として、スカラー（\(n=1\) ）の純粋な2次関数 \(g(x) = g_0 + \frac{g_2}{2}(x - \bar{x})^2\) を考えます。\(X \sim \mathcal{N}(\bar{x}, \sigma^2)\) なら、ガウス分布の4次モーメント \(E[(X-\bar{x})^4] = 3\sigma^4\) を使って真の分散は：

一方、UTでこの \(g\) にシグマ点を通して共分散を計算すると（\(n=1\) : \(\boldsymbol{\chi}_0 = \bar{x}\) , \(\boldsymbol{\chi}_{1,2} = \bar{x} \pm \gamma\sigma\) ）、SymPyで記号計算すると：

両者を等置すると \(\alpha^2 \kappa + \beta = 2\) 、すなわち：

典型的な設定 \(\kappa = 0\) では、\(\alpha\) の値によらず常に \(\beta = 2\) が厳密解になります。これが「ガウス分布では \(\beta = 2\) が最適」という経験則の根拠です（SymPyで \(\alpha \in \{0.001, 1\}\) , \(\kappa=0\) について数値的に検算し、いずれも \(\beta=2\) を確認済み）。\(\kappa \ne 0\) を使う流儀（例えば \(\kappa = 3-n\) ）では、この関係式が示すとおり \(\beta\) 側の補正も調整が必要になります。

変換後の統計量

シグマポイントを非線形関数 \(f\) に通した後:

UKF アルゴリズム

予測ステップ

1. 現在の状態推定 \((\hat{\mathbf{x}}_{k-1}, \mathbf{P}_{k-1})\) からシグマポイントを生成
2. 状態遷移関数 \(f\) で変換: \(\boldsymbol{\chi}_{k|k-1} = f(\boldsymbol{\chi}_{k-1})\)
3. 予測平均と共分散を計算:

更新ステップ

1. 予測シグマポイントを観測関数 \(h\) で変換
2. カルマンゲインを計算:

3. 状態と共分散を更新:

EKFとの比較

Python実装

UKFの強みは非線形の観測または状態遷移に対して発揮されます。ここでは、レーダーが原点から目標までの距離（range）と方位角（bearing）を測定する、古典的なレンジ・ベアリング観測モデルを例にとり、EKFとUKFを同一データ・同一条件で直接比較します。この観測モデルは EKFの記事 でも扱った非線形モデルです。

問題設定

状態ベクトル \(\mathbf{x} = [p_x, p_y, v_x, v_y]^T\) （位置と速度）に対して、状態遷移は等速直線運動モデルで線形です：

一方、観測モデルは、原点にいるセンサーが距離 \(r\) と方位角 \(\theta\) を計測する非線形関数です：

EKFではこの \(h\) をヤコビ行列で線形化します。\(r = \sqrt{p_x^2 + p_y^2}\) として：

UKFはこのヤコビ行列を一切使わず、シグマポイントを \(h\) に直接通すことで非線形性を近似します。状態遷移 \(f\) は線形なのでヤコビ行列 \(\mathbf{F}\) は定数行列（近似誤差なし）であり、EKFとUKFの差はもっぱら非線形観測 \(h\) の扱いから生まれます。

シミュレーション設定

目標がセンサー（原点）からわずか距離1まで接近する軌道を等速直線運動で追跡します。あえて初期推定を真値から大きくずらし、初期共分散も大きく設定しています（追尾再獲得のシナリオを想定）。これにより、観測関数 \(h\) の曲率が大きい原点近傍で線形化誤差の影響が現れやすい状況を作り、EKFとUKFの違いを際立たせます。

import numpy as np
import matplotlib.pyplot as plt

class UKF:
    def __init__(self, n, m, f, h, Q, R, alpha=1.0, beta=2, kappa=0):
        self.n = n
        self.m = m
        self.f = f
        self.h = h
        self.Q = Q
        self.R = R

        self.lam = alpha**2 * (n + kappa) - n
        self.gamma = np.sqrt(n + self.lam)

        # 重み
        self.Wm = np.full(2 * n + 1, 1 / (2 * (n + self.lam)))
        self.Wc = np.full(2 * n + 1, 1 / (2 * (n + self.lam)))
        self.Wm[0] = self.lam / (n + self.lam)
        self.Wc[0] = self.lam / (n + self.lam) + (1 - alpha**2 + beta)

    def sigma_points(self, x, P):
        L = np.linalg.cholesky(P)
        sigmas = np.zeros((2 * self.n + 1, self.n))
        sigmas[0] = x
        for i in range(self.n):
            sigmas[i + 1] = x + self.gamma * L[:, i]
            sigmas[i + 1 + self.n] = x - self.gamma * L[:, i]
        return sigmas

    def predict(self, x, P, dt):
        sigmas = self.sigma_points(x, P)
        sigmas_pred = np.array([self.f(s, dt) for s in sigmas])

        x_pred = np.dot(self.Wm, sigmas_pred)
        P_pred = self.Q.copy()
        for i in range(2 * self.n + 1):
            d = sigmas_pred[i] - x_pred
            P_pred += self.Wc[i] * np.outer(d, d)

        return x_pred, P_pred, sigmas_pred

    def update(self, x_pred, P_pred, sigmas_pred, z):
        sigmas_z = np.array([self.h(s) for s in sigmas_pred])

        z_pred = np.dot(self.Wm, sigmas_z)
        Pzz = self.R.copy()
        Pxz = np.zeros((self.n, self.m))

        for i in range(2 * self.n + 1):
            dz = sigmas_z[i] - z_pred
            dz[1] = (dz[1] + np.pi) % (2 * np.pi) - np.pi  # 方位角の正規化
            dx = sigmas_pred[i] - x_pred
            Pzz += self.Wc[i] * np.outer(dz, dz)
            Pxz += self.Wc[i] * np.outer(dx, dz)

        y = z - z_pred
        y[1] = (y[1] + np.pi) % (2 * np.pi) - np.pi

        K = Pxz @ np.linalg.inv(Pzz)
        x_new = x_pred + K @ y
        P_new = P_pred - K @ Pzz @ K.T

        return x_new, P_new


class EKF:
    """比較用のEKF（https://yuhi-sa.github.io/posts/20260224_ekf/1/ と同じ構造）"""

    def __init__(self, x0, P0, Q, R, dt):
        self.x = x0.copy()
        self.P = P0.copy()
        self.Q = Q
        self.R = R
        self.dt = dt

    def predict(self):
        F = jacobian_F(self.dt)
        self.x = state_transition(self.x, self.dt)
        self.P = F @ self.P @ F.T + self.Q

    def update(self, z):
        H = jacobian_H(self.x)
        y = z - observation(self.x)
        y[1] = (y[1] + np.pi) % (2 * np.pi) - np.pi

        S = H @ self.P @ H.T + self.R
        K = self.P @ H.T @ np.linalg.inv(S)
        self.x = self.x + K @ y
        self.P = (np.eye(4) - K @ H) @ self.P

# --- レンジ・ベアリング観測モデル（式14-16） ---
def state_transition(x, dt):
    px, py, vx, vy = x
    return np.array([px + vx * dt, py + vy * dt, vx, vy])

def jacobian_F(dt):
    return np.array([
        [1, 0, dt, 0],
        [0, 1, 0, dt],
        [0, 0, 1, 0],
        [0, 0, 0, 1],
    ])

def observation(x):
    px, py = x[0], x[1]
    r = np.sqrt(px**2 + py**2)
    theta = np.arctan2(py, px)
    return np.array([r, theta])

def jacobian_H(x):
    px, py = x[0], x[1]
    r = np.sqrt(px**2 + py**2)
    return np.array([
        [px / r,      py / r,    0, 0],
        [-py / r**2,  px / r**2, 0, 0],
    ])

# --- シミュレーション：センサー近傍を通過する目標の追跡 ---
np.random.seed(42)
dt = 1.0
T = 80

# 真の軌跡（センサーから最短距離1まで接近する等速直線運動）
t = np.arange(0, T + 1) * dt
px_true = np.full(T + 1, 1.0)
py_true = -10.0 + 0.25 * t
true_states = np.column_stack([
    px_true, py_true, np.zeros(T + 1), np.full(T + 1, 0.25)
])

R = np.diag([0.3, 0.2])   # 距離ノイズ分散, 方位角ノイズ分散(rad^2)
Q = np.diag([0.05, 0.05, 0.01, 0.01])

measurements = np.array([
    observation(true_states[k]) + np.random.multivariate_normal([0, 0], R)
    for k in range(1, T + 1)
])

# 真値から大きくずれた初期推定・大きな初期共分散から開始（追尾再獲得）
x0 = np.array([6.0, -3.0, 0.0, 0.1])
P0 = np.diag([16.0, 16.0, 1.0, 1.0])

# --- EKFの実行 ---
ekf = EKF(x0.copy(), P0.copy(), Q, R, dt)
ekf_estimates = [x0.copy()]
for k in range(T):
    ekf.predict()
    ekf.update(measurements[k].copy())
    ekf_estimates.append(ekf.x.copy())
ekf_estimates = np.array(ekf_estimates)

# --- UKFの実行（同一データ・同一初期値・同一ノイズ共分散） ---
ukf = UKF(n=4, m=2, f=state_transition, h=observation, Q=Q, R=R)
x_est, P_est = x0.copy(), P0.copy()
ukf_estimates = [x_est.copy()]
for k in range(T):
    x_pred, P_pred, sigmas_pred = ukf.predict(x_est, P_est, dt)
    x_est, P_est = ukf.update(x_pred, P_pred, sigmas_pred, measurements[k].copy())
    ukf_estimates.append(x_est.copy())
ukf_estimates = np.array(ukf_estimates)

# --- RMSE比較 ---
ekf_rmse = np.sqrt(np.mean(
    np.sum((true_states[1:, :2] - ekf_estimates[1:, :2])**2, axis=1)))
ukf_rmse = np.sqrt(np.mean(
    np.sum((true_states[1:, :2] - ukf_estimates[1:, :2])**2, axis=1)))

print(f"EKF Position RMSE: {ekf_rmse:.4f}")
print(f"UKF Position RMSE: {ukf_rmse:.4f}")

結果

同一の観測データ・同一の初期推定・同一のノイズ共分散でEKFとUKFを実行して位置RMSE（真値とのユークリッド距離の二乗平均平方根）を比較すると、EKF推定のRMSEは1.7155、UKF推定のRMSEは1.5489となり、UKFはEKFに対して約9.7%の誤差削減を達成しました（乱数シード42で再現可能）。

左図は目標の真の軌跡・観測値（デカルト座標換算）・EKF推定・UKF推定を重ねたものです。右図は各時刻での位置誤差（真値からのユークリッド距離）の推移で、破線はそれぞれのRMSEを示します。軌道序盤（追尾再獲得直後、目標がセンサー近傍を通過する高曲率領域）でUKFの誤差がEKFを下回る場面が多く見られ、シグマポイントによる非線形変換の直接近似が、ヤコビアンによる一次線形化よりも観測関数の曲率を精度良く捉えていることが分かります。一方、目標がセンサーから十分離れて観測関数がほぼ線形とみなせる終盤では、両者の誤差はほぼ一致します。前掲の「EKFとの比較」表が示す「非線形精度：中程度（EKF） vs 高い（UKF）」という差は、非線形性が強い局面でこそ顕在化することが、この数値実験からも裏付けられます。

実践上の注意点：エッジケースとチューニング

前節で導出した重みは、線形変換に対する平均・共分散の厳密な再現と、ガウス事前分布に対する4次モーメント（尖度）の補正を保証します。しかし \(\alpha\) を極端な値に設定したり、フィルタを長時間運用したりすると、理論上は無害なはずの負の重みや、共分散の数値的破綻が実際に問題化することがあります。以下、3つの典型的な落とし穴をPythonで実行し、実測した数値とともに示します。

エッジケース①：負の重みと破局的な打ち消し

\(\kappa=0\) の場合 \(\lambda = (\alpha^2-1)n\) なので、典型的な推奨範囲 \(\alpha \in [10^{-4}, 1)\) ではほぼ常に \(\lambda < 0\) となり、\(W_0^{(c)}\) も負になります。\(\alpha=0.5, \beta=2, \kappa=0\) で実際に計算すると：

式(5b)で示した通り、線形変換に対しては \(W_0^{(c)}\) が負でも共分散は厳密に再現されます。実際に \(n=2\) , \(\alpha=0.5\) （\(W_0^{(c)}=-0.25\) ）で適当な線形写像 \(\mathbf{A}\) を通した共分散をUTで計算し、解析解 \(\mathbf{A P A}^T\) と比較すると、最大絶対誤差は \(4.44 \times 10^{-16}\) （倍精度の丸め誤差の範囲内）でした。負の重み自体は問題ではありません。

本当に問題になるのは、\(\alpha\) を極端に小さくしたときの**破局的な打ち消し（catastrophic cancellation）**です。\(\alpha \to 0\) では \(n+\lambda \to 0\) となり、シグマ点の広がり \(\gamma=\sqrt{n+\lambda}\) は0に近づく一方、重みの大きさ \(\sim 1/(n+\lambda)\) は発散します。つまり「ほぼ等しい2つの数の差」（ずれ \(\boldsymbol{\mathcal{Y}}_i - \bar{\mathbf{y}}\) ）に「非常に大きい重み」を掛け合わせる、丸め誤差が増幅されやすい典型的な構成になります。単精度（float32、組み込み・GPU実装で現実的）で \(\bar{\mathbf{x}}=[1000,-500]\) , \(\mathbf{P}=\mathrm{diag}(2.0,1.5)\) に恒等変換（\(f(\mathbf{x})=\mathbf{x}\) 、理論値は \(\mathbf{P}\) そのもの）を適用する、次のコードで実際に確認できます。

import numpy as np

def ut_identity_f32(alpha, n=2, kappa=0.0, beta=2.0):
    x = np.array([1000.0, -500.0], dtype=np.float32)
    P = np.diag(np.array([2.0, 1.5], dtype=np.float32))
    lam = np.float32(alpha**2 * (n + kappa) - n)
    c = np.float32(n + lam)
    L = np.linalg.cholesky((c * P).astype(np.float64)).astype(np.float32)

    sigmas = np.zeros((2 * n + 1, n), dtype=np.float32)
    sigmas[0] = x
    for i in range(n):
        sigmas[i + 1] = x + L[:, i]
        sigmas[n + i + 1] = x - L[:, i]

    Y = sigmas  # f(x) = x (恒等変換)
    Wm0 = np.float32(lam / c)
    Wc0 = np.float32(lam / c + (1 - alpha**2 + beta))
    Wi = np.float32(1.0 / (2 * c))

    ymean = Wm0 * Y[0]
    for i in range(1, 2 * n + 1):
        ymean = ymean + Wi * Y[i]

    Py = Wc0 * np.outer(Y[0] - ymean, Y[0] - ymean)
    for i in range(1, 2 * n + 1):
        Py = Py + Wi * np.outer(Y[i] - ymean, Y[i] - ymean)

    eig = np.linalg.eigvalsh(Py.astype(np.float64))
    return Py, eig

for alpha in [1.0, 0.1, 0.01, 0.001, 0.0001]:
    with np.errstate(all="ignore"):
        try:
            Py, eig = ut_identity_f32(alpha)
        except np.linalg.LinAlgError:
            print(f"alpha={alpha}: c*P のコレスキー分解が失敗(float32丸めで半正定値性を喪失)")
            continue
        psd = "Yes" if eig.min() >= 0 else "No"
        print(f"alpha={alpha}: Py[0,0]={Py[0,0]:.6f} Py[1,1]={Py[1,1]:.6f} min_eig={eig.min():.3e} PSD={psd}")

実行結果：

alpha=1.0: Py[0,0]=2.000000 Py[1,1]=1.500008 min_eig=1.500e+00 PSD=Yes
alpha=0.1: Py[0,0]=2.000307 Py[1,1]=1.500237 min_eig=1.500e+00 PSD=Yes
alpha=0.01: Py[0,0]=3.003601 Py[1,1]=1.752060 min_eig=1.577e+00 PSD=Yes
alpha=0.001: Py[0,0]=12544.000000 Py[1,1]=65.250000 min_eig=1.811e+00 PSD=Yes
alpha=0.0001: c*P のコレスキー分解が失敗(float32丸めで半正定値性を喪失)

（真値は \(\mathrm{diag}(2.0, 1.5)\) ）\(\alpha=10^{-3}\) の時点で再構成された共分散はすでに真値から3〜4桁ずれています。\(\alpha=10^{-4}\) では \(c\mathbf{P}=(n+\lambda)\mathbf{P}\) 自体がfloat32の丸めで半正定値性を失い、コレスキー分解の時点で失敗します（実装や丸めの順序によっては、この手前で \(P_y\) が非半正定値な行列として得られたり、重み計算がオーバーフロー/アンダーフローしNaNになったりすることもあります——正確にどの \(\alpha\) でどちらの破綻モードが先に起こるかは、演算の順序や使用するBLAS実装に依存します）。倍精度でも同じ破綻モードは起こり、現れる \(\alpha\) が小さくなるだけです（同一設定で倍精度なら \(\alpha \sim 10^{-7}\text{-}10^{-8}\) 付近から数値のずれが目立ち始め、\(\alpha=10^{-9}\) では \(n+\lambda\) が完全に0にアンダーフローしゼロ除算エラーになります）。

\(\alpha\) が小さくなるほどシグマ点は中心点 \(\boldsymbol{\chi}_0\) に密集し、外側の点と中心点がほぼ区別できなくなります。このほぼ縮退した配置で、大きさが急増する重み（\(W_0^{(c)} \to -\infty\) ）を使って差を取るため、丸め誤差が増幅されるのです。実務上の指針：単精度なら \(\alpha \gtrsim 10^{-3}\) 、倍精度でも \(\alpha \gtrsim 10^{-6}\) 程度を下限の目安とし、それより小さい \(\alpha\) が必要な場面では、共分散 \(\mathbf{P}\) を毎回再構成するのではなくコレスキー因子を直接伝播する Square-Root UKF の採用を検討してください。

エッジケース②：コレスキー分解の数値不安定性

\(\alpha\) を極端にしなくても、通常運用中の丸め誤差の蓄積だけで \(\mathbf{P}\) が半正定値性を失うことがあります。とくに観測が強く相関している（あるいはほぼ重複した）場合、\(\mathbf{P}_{yy}\) がほぼ特異になり、式(13)の素朴な更新 \(\mathbf{P}_k = \mathbf{P}_k^{-} - \mathbf{K}_k \mathbf{P}_{yy} \mathbf{K}_k^T\) （Joseph形式でない更新）がわずかな丸め誤差を増幅します。次のコードは、ほぼ重複した2つの観測（片方をもう片方に \(10^{-9}\) だけずらしたもの）と小さい観測ノイズ \(\mathbf{R}=10^{-8}\mathbf{I}\) を使って、この破綻を再現したものです。

import numpy as np
from scipy.linalg import sqrtm

P = np.diag([1.0, 1.0, 0.5, 0.5])
H = np.array([[1.0, 0.0, 0.0, 0.0], [1.0 + 1e-9, 0.0, 0.0, 0.0]])
R = np.diag([1e-8, 1e-8])

failed_at = None
for k in range(200):
    Pyy = H @ P @ H.T + R
    Pxy = P @ H.T
    K = Pxy @ np.linalg.inv(Pyy)
    P = P - K @ Pyy @ K.T
    try:
        np.linalg.cholesky(P)
    except np.linalg.LinAlgError:
        failed_at = k
        break

print(f"Cholesky first failed at iteration: {failed_at}")
eig = np.linalg.eigvalsh((P + P.T) / 2)
print(f"Eigenvalues of symmetrized P: {eig}")

実行結果：コレスキー分解は最初の更新（イテレーション0）で早くも失敗し、対称化した \(\mathbf{P}\) の最小固有値は \(-6.90 \times 10^{-9}\) でした。\(\mathbf{P}\) 自体のフロベニウスノルムでの非対称性はちょうど0.0であり、非対称性が原因ではなく、ほぼ特異な \(\mathbf{P}_{yy}\) の逆行列計算が生む丸め誤差によって固有値がゼロをわずかに下回ったことが原因です。標準的な対処法は3つあります。

1. 固有値クリッピング：\(\mathbf{P}\) を対称化し固有値分解、負の固有値を小さな正の下限（例: \(10^{-10}\) ）にクリップして再構成 → クリップ後の固有値は \([0, 0.5, 0.5, 1.0]\) となり、以降のコレスキー分解は成功。
2. ジッタ加算 \(\mathbf{P} + \epsilon \mathbf{I}\) ：この例では \(\epsilon=10^{-12}, 10^{-10}\) ではまだ失敗し、\(\epsilon=10^{-8}\) で初めて成功しました。ジッタの大きさは問題依存であり固定値を過信しないこと。
3. 対称行列平方根（scipy.linalg.sqrtm）：不定行列に対しても直接計算でき、虚部は \(8.3 \times 10^{-5}\) 程度の無視できる大きさに収まるため .real を取れば実用上問題ありません。

本番運用では、素朴な更新式(13)の代わりに数値的に安定な Joseph形式

を使うか、上記のクリッピング／ジッタを共分散更新の直後に毎回挟むことを推奨します。

エッジケース③：\(\alpha, \beta, \kappa\) の選び方（数値スイープで検証）

実践的な指針は次の通りです。

• \(\alpha\) ：シグマ点の広がりを制御します。小さいほど平均近傍の情報に絞られ高次の近似誤差を抑えられますが、上記の通り小さすぎると数値的に破綻します。\(10^{-4} \sim 1\) の範囲で選び、既定値は \(10^{-3}\) です。
• \(\beta\) ：ガウス事前分布では式(5c)より \(\kappa=0\) のとき厳密に \(\beta=2\) です。観測・プロセスノイズがガウスから大きく外れると分かっている場合のみ変更を検討します。
• \(\kappa\) ：\(0\) が安全な既定値です。歴史的には \(\kappa=3-n\) （\(n\) は状態次元）という発見的な選び方もありますが、\(n \ge 3\) では \(\lambda + n = \alpha^2 \cdot 3\) となり \(n\) に対して相対的に小さくなるため、高次元では \(\kappa=0\) の方が扱いやすいのが実情です。

本記事のレンジ・ベアリング追跡の題材で、\(\beta=2, \kappa=0\) を固定し \(\alpha\) を \(0.001\) から \(1\) までスイープしてUKF位置RMSEを比較しました。

最小RMSEは \(\alpha=0.01\) で 1.510643、最大は \(\alpha=0.7\) で 1.560870 となり、その差はわずか 3.32% でした。本問題では状態遷移 \(f\) が線形で非線形性は観測 \(h\) のみに由来するため、\(\alpha\) の影響は比較的穏やかです。これは実務上重要な示唆で、多くの問題では既定値 \(\alpha=10^{-3}\) （またはその近傍の \([10^{-3}, 10^{-1}]\) ）を使えば十分であり、\(\alpha\) のチューニングに多くの時間をかける必要はありません。ただし、状態遷移そのものが強く非線形な問題や、上記エッジケース①のような純2次関数に近い局面では、\(\alpha\) の影響がより大きく現れうるため、本記事のようなRMSEスイープで実データに対して検証することを推奨します。

おすすめ書籍

カルマンフィルタの導出から実装例までを日本語で丁寧に解説した定番書です。状態空間モデルの理解が深まります。

※ 上記は Amazon アソシエイトのリンクです。

関連記事

• カルマンフィルタの理論とPython実装 - UKFの基盤となる線形カルマンフィルタの理論と実装を解説しています。
• 拡張カルマンフィルタ（EKF）の理論とPython実装 - UKFの前身であるEKFの線形化アプローチを解説しています。
• Cubature Kalman Filter（CKF）の理論とPython実装 - UKFの重み問題を解決するキュバチャカルマンフィルタを解説しています。
• パーティクルフィルタのPython実装 - UKFより高い非線形性にはパーティクルフィルタが有効です。
• カルマンスムーザ（RTS Smoother）の理論とPython実装 - UKFベースのスムーザへの拡張の基盤となるRTSスムーザを解説しています。
• カルマンスムーザの比較 - UKFベースの非線形スムーザを含むスムージング手法を横断比較しています。
• 非線形カルマンスムーザ（EKS/URTS） - 本記事のシグマポイントを後向きパスに用いるUnscented RTS Smootherを実装し、EKFベースのEKSと精度比較しています。
• Unscented変換の基礎 - UKFの核となるUnscented変換の詳細な解説です。
• ウィーナーフィルタの理論とPython実装 - 定常過程における最適線形フィルタ。UKFとの位置づけの違いを理解できます。
• 指数移動平均（EMA）フィルタの周波数特性 - 最も基本的なフィルタとの比較で、UKFの高度さが分かります。
• 機械学習による時系列予測・分類・異常検知ハブ - UKF を線形カルマン・EKF・粒子フィルタ・LSTM の系列モデル群に並べ、非線形状態空間における予測・異常検知の使い分けを整理したハブ記事。

参考文献

• Julier, S. J., & Uhlmann, J. K. (2004). “Unscented Filtering and Nonlinear Estimation”. Proceedings of the IEEE, 92(3), 401-422.
• Wan, E. A., & van der Merwe, R. (2000). “The Unscented Kalman Filter for Nonlinear Estimation”. Proceedings of the IEEE Adaptive Systems for Signal Processing, Communications, and Control Symposium.
• Simon, D. (2006). Optimal State Estimation. Wiley. Chapter 14.

近年の研究動向（2024-2025）

• Yang, L., Lin, X., Hou, Y., Ren, J., & Wang, M. (2025). “Application of an Improved Adaptive Unscented Kalman Filter in Vehicle Driving State Parameter Estimation”. International Journal of Adaptive Control and Signal Processing. カイ二乗検定とZスコア正規化を組み合わせた適応UKF（CAUKF）で異常値を補正し、車両運動状態推定の頑健性を高めています。
• Zhao, J., Zhang, Y., Li, S., Wang, J., Fang, L., Ning, L., Feng, J., & Zhang, J. (2025). “An Improved Unscented Kalman Filter Applied to Positioning and Navigation of Autonomous Underwater Vehicles”. Sensors, 25(2), 551. ローリングホライズン推定に基づく適応UKF（RHAUKF）により、AUVの測位精度を向上させています。
• Nguyen, D. V., Zhao, H., Hu, J., & Giang, L. N. (2025). “Adaptive Robust Unscented Kalman Filter for Dynamic State Estimation of Power System”. arXiv:2504.07731. 一般化最小混合誤差エントロピーに基づく適応UKFで、非ガウス性雑音下の電力系統状態推定を頑健化しています。