簡介面向高超聲速滑翔飛行器軌跡預測研究和防御決策需求本資料為目標識別、意圖推斷與軌跡預測相關科研人員和工程師提供完整復現(xiàn)方案。內(nèi)容圍繞貝葉斯推斷框架展開先利用攻擊意圖和戰(zhàn)場態(tài)勢信息構(gòu)建意圖代價函數(shù)再遞推機動模式與運動狀態(tài)并結(jié)合蒙特卡洛序貫濾波計算目標狀態(tài)分布與攻擊概率粒子濾波、測量更新等核心模塊均有對應代碼與逐段解釋便于理解非高斯非線性條件下的預測建模和實現(xiàn)流程并附仿真測試思路。資源包共1個文件為docx文檔大小53KB內(nèi)容集中易查閱閱讀與實踐都很方便。已有78人學習適合具備一定編程和數(shù)學基礎、希望結(jié)合代碼快速掌握貝葉斯推斷與粒子濾波應用的讀者。1. 高超聲速滑翔目標軌跡預測為什么難貝葉斯推斷從哪一步開始起作用雷達屏幕上一個高超聲速滑翔目標HGRVHypersonic Glide Reentry Vehicle進入滑翔段以后常規(guī)EKF給出的落點預報往往在第一輪橫向機動后就偏出幾十公里。這不是濾波代碼寫得差而是模型假設不匹配HGRV的攻角和傾側(cè)角都能在飛行中改變狀態(tài)分布會裂成多個峰高斯假設扛不住?；谪惾~斯推斷的軌跡預測方法放開這條限制用粒子集直接逼近后驗分布落點區(qū)間的可靠性明顯提升。這篇文章把它拆成四件事運動建模、貝葉斯框架、粒子濾波代碼實現(xiàn)、仿真驗證與調(diào)參適合做雷達數(shù)據(jù)處理和飛行器制導仿真的工程師逐步復現(xiàn)。2. 先立住運動模型HGRV的氣動方程、平衡滑翔近似與三段式運動2.1 HGRV與彈道目標的本質(zhì)區(qū)別升力改變了整個軌跡預測的邊界常規(guī)再入目標的飛行軌跡由重力和稀薄大氣阻力決定幾乎是一條確定的可預測彈道濾波器只需要在線修正初始速度偏差。HGRV不同它在進入滑翔段后靠氣動升力長時間維持高度可以在約25到60公里高度之間做跳躍滑翔并用傾側(cè)角改變航向。升力讓預測問題從“單峰參數(shù)估計”變成了“多峰行為推斷”目標一旦切換傾側(cè)角方向未來的軌跡可能左偏也可能右偏。雷達能看到的只是當前狀態(tài)看不到目標內(nèi)部的制導策略所以任何單一標稱軌跡假設都會在下一次機動時失效。這是HGRV軌跡預測的第一條邊界你預測的不是一條軌跡而是一族軌跡的分布。貝葉斯推斷天然處理這種分布這也是后文所有推導的出發(fā)點。2.2 點質(zhì)量模型的三個核心方程阻力、升力、重力與傾側(cè)角的作用工程里做HGRV軌跡預測很少直接用六自由度剛體模型計算量太大且攻角、舵偏等輸入根本拿不到。常見做法是使用速度坐標系下的點質(zhì)量模型把氣動力合成為升力加速度和阻力加速度。模型寫作dv/dt -D - g·sin(γ)dγ/dt L·cos(φ)/v - (g - v2/r)·cos(γ)/vdψ/dt L·sin(φ)/(v·cos(γ))其中 v 為速度γ 為航跡角ψ 為航向角φ 為傾側(cè)角r 為地心距。升力加速度 L 和阻力加速度 D 的表達式為L 0.5·ρ·v2·S·C_L / mD 0.5·ρ·v2·S·C_D / mρ 是大氣密度隨高度按負指數(shù)近似衰減S 是參考面積m 是飛行器質(zhì)量C_L 和 C_D 是升力、阻力系數(shù)。C_L 與 C_D 不僅隨攻角變化還受馬赫數(shù)和高度耦合影響而這個攻角恰恰是外部觀測拿不到的。所以 L/D 升阻比在這個場景里是個黑匣子典型值在2到4之間游走濾波器必須對它保留足夠的余量。傾側(cè)角 φ 的作用值得單獨說明。縱向通道由 L·cos(φ) 調(diào)制橫向通道由 L·sin(φ) 驅(qū)動。φ 為0時目標做縱向跳躍滑翔φ 為正負值時軌跡開始側(cè)向轉(zhuǎn)彎。制導策略里的“橫向機動”在雷達看來就是傾側(cè)角符號連續(xù)切換這一條是卡爾曼類方法翻車的主要來源。模型里沒有推力項因為滑翔段無動力這一假設在大多數(shù)HGRV場景成立。2.3 運動分段與過程噪聲設置滑翔段、橫向機動段、末段快降的區(qū)別把HGRV的飛行過程按可觀測特征拆成三段有助于分別設置過程噪聲和預報策略。階段典型高度范圍典型速度范圍機動特征建模重點滑翔巡航段40~60 kmMa 12~20縱向跳躍橫向機動少平衡滑翔近似過程噪聲可收緊橫向機動段25~45 kmMa 8~15傾側(cè)角頻繁切換橫向大范圍轉(zhuǎn)彎多峰分布過程噪聲必須放大末段快降段10~25 kmMa 4~8高度快速下降軌跡趨于陡直幾何外推為主落點約束生效在滑翔巡航段目標近似滿足平衡滑翔條件升力垂直分量約等于重力與離心力之差縱向加速度趨近于零。這個條件可以作為先驗約束把粒子的初始散布限制在物理可行的包線內(nèi)避免濾波前期粒子亂飛。橫向機動段是預測誤差被拉大的主戰(zhàn)場我一般會把過程噪聲中的航向角方差放大3到5倍給粒子足夠的自由度去覆蓋左偏和右偏兩簇軌跡。末段快降階段的落點預測不再是單純的狀態(tài)遞推而是要結(jié)合雷達測距幾何做約束。目標高度快速下降速度攀升軌跡接近陡直彈道此時量測更新率如果還停留在1秒1次橫向誤差會被放大建議把雷達數(shù)據(jù)率提高到5赫茲以上。過程噪聲設置的總原則是機動越強噪聲越大寧可讓粒子集散一點也不要讓濾波器過早鎖死在一條錯誤軌跡上。3. 貝葉斯推斷的高斯困境后驗遞推公式與粒子濾波的五個操作3.1 為什么在HGRV軌跡預測場景下卡爾曼的高斯假設會翻車標準的卡爾曼濾波假設系統(tǒng)噪聲和量測噪聲都是高斯分布EKF在線性化點附近做一階展開UKF用sigma點傳播非線性但最后仍然把后驗分布壓縮成一個高斯。HGRV的問題在于機動切換會產(chǎn)生真正的多峰分布目標到達某個航路點時向左轉(zhuǎn)和向右轉(zhuǎn)的概率相近狀態(tài)后驗會裂成兩簇。高斯近似會把兩個峰壓成一個峰均值落在兩峰之間的空白地帶那個位置實際上沒有任何物理軌跡會經(jīng)過。我見過不少用UKF做高超聲速滑翔目標預測的方案前幾秒跟蹤精度尚可第一次大傾側(cè)切換后協(xié)方差橢圓開始變得特別大預報落點偏向兩峰中間。這不是調(diào)參數(shù)能救的是分布假設錯了。貝葉斯推斷在這里的價值不是“換個更高級的濾波公式”而是放棄對后驗分布形式的預先設定用加權(quán)樣本集去逼近真實分布。粒子濾波就是這套思路最直接的落地工具。3.2 遞歸貝葉斯后驗公式與粒子濾波五個操作采樣、預測、更新、歸一化、重采樣貝葉斯推斷在軌跡預測語境下就是遞歸貝葉斯狀態(tài)估計。給定直到當前時刻的所有量測 Z_{1:k}目標狀態(tài) x_k 的后驗概率密度可以寫成p(x_k | Z_{1:k}) c·p(z_k | x_k)·∫ p(x_k | x_{k-1})·p(x_{k-1} | Z_{1:k-1}) dx_{k-1}其中 p(x_k | x_{k-1}) 來自上一章的運動方程p(z_k | x_k) 來自雷達量測方程c 是歸一化常數(shù)。這個積分對HGRV的非線性模型沒有解析解粒子濾波用 N 個帶權(quán)重的樣本 {x_i, w_i} 近似它。粒子濾波的標準循環(huán)是五個操作。第一步采樣從建議分布中生成新粒子第二步預測用HGRV的運動方程把每個粒子向前推進一個步長第三步更新按雷達量測似然調(diào)整每個粒子的權(quán)重第四步歸一化讓所有權(quán)重之和為1第五步重采樣按照權(quán)重重新抽取粒子把資源集中到高概率區(qū)域。工程里常用有效粒子數(shù) N_eff 1 / Σ(w_i2) 來判斷是否需要重采樣閾值取 0.3N 到 0.5N 之間。建議分布我直接用先驗轉(zhuǎn)移分布 p(x_k | x_{k-1})實現(xiàn)最簡單代價是過程噪聲較大時粒子會發(fā)散需要配合重采樣閾值控制。3.3 狀態(tài)向量怎么選在位置速度上再加一個氣動修正維狀態(tài)向量設計直接決定濾波器能不能收斂?；A六維狀態(tài)是本地切平面坐標下的位置和速度加兩個姿態(tài)角x [px, py, pz, v, γ, ψ]。其中 px、py、pz 是雷達站ENU坐標系下的目標位置v 是速度大小γ 是航跡角ψ 是航向角。這個狀態(tài)向量覆蓋了HGRV軌跡預測所需的全部可觀測信息。實際問題比這更麻煩運動方程里的 C_L 和 C_D 是未知的模型失配會讓濾波器長期預報系統(tǒng)性偏置。我常用的補法是給狀態(tài)向量增加一個升力系數(shù)修正維 c_l它的量綱是乘子初始散布在0.8到1.2之間用隨機游走描述變化。這樣濾波器在量測更新時會自動修正氣動偏差相當于給模型失配留了一顆后悔藥。c_l 的過程噪聲方差給太小濾波器不敢認錯給太大氣動修正會被噪聲淹沒。工程經(jīng)驗值取 1e-4 到 1e-3 量級具體數(shù)值需要結(jié)合仿真調(diào)校。如果想同時修正升力和阻力可以把升阻比 L/D 也作為狀態(tài)維形成八維狀態(tài)。但粒子濾波的維度越高需要的粒子數(shù)越多八維往往要三倍以上的粒子才能維持同樣精度。實際工程里我更傾向于固定 L/D只修正乘性升力系數(shù)用一個簡化的模型換回計算余量。4. 把貝葉斯環(huán)套到HGRV上狀態(tài)向量、量測方程與一次完整遞推4.1 狀態(tài)方程與微分方程把氣動不確定性設計成修正維為了和后面的代碼對齊這里用ENU本地切平面坐標給出完整的連續(xù)狀態(tài)方程。狀態(tài)向量 x [px, py, pz, v, γ, ψ]位置和姿態(tài)的導數(shù)如下dpx/dt v·cos(γ)·cos(ψ)dpy/dt v·cos(γ)·sin(ψ)dpz/dt v·sin(γ)dv/dt -D - g·sin(γ)dγ/dt L·cos(φ)/v - (g - v2/r)·cos(γ)/vdψ/dt L·sin(φ)/(v·cos(γ))大氣密度按指數(shù)近似 ρ ρ?·exp(-pz/H)升力和阻力加速度引用第二張公式。這里的簡化是沒考慮地球曲率對側(cè)向運動的影響短時預報誤差在可接受范圍內(nèi)如果做長時程滑翔段預報需要把經(jīng)緯高坐標系下的曲率項補進去。離散化時我使用四階龍格庫塔積分步長取0.1到0.2秒。這個選擇的原因很直接歐拉法在疊加過程噪聲后容易漂二階方法在HGRV這種強非線性方程上又不劃算RK4在0.1秒步長下數(shù)值表現(xiàn)可靠代碼也短。每個量測幀之間做一至五次積分雷達幀率常見1赫茲也就是每幀之間走10個0.1秒步長。4.2 雷達量測方程與噪聲矩陣距離、方位角、俯仰角的典型配置雷達量測通常提供目標相對于雷達站的斜距 r、方位角 az、俯仰角 el。量測方程寫作r sqrt(px2 py2 pz2)az atan2(py, px)el asin(pz / r)量測噪聲矩陣 R 按雷達精度配置典型工程值如下表。注意角度量測要經(jīng)過 ±π 環(huán)繞處理否則粒子濾波更新時角度差在跨越零度線時會產(chǎn)生虛假大新息。量測量符號典型噪聲標準差說明斜距r50~150 m高超聲速目標回波信噪比波動大方位角az0.1°~0.5°取決于雷達波束寬度俯仰角el0.1°~0.5°低空目標需注意多路徑效應粒子濾波相比EKF有個明顯的實現(xiàn)優(yōu)勢更新時不需要計算雅可比矩陣只需要做一次量測方程的函數(shù)調(diào)用。復雜轉(zhuǎn)角公式也不用求導改模型時省大量開發(fā)時間。代價是每個粒子都要獨立算出量測預測計算量正比于粒子數(shù)。4.3 一次完整遞推的七個步驟從初始化到落點預報把前面的方程組裝成一整套遞推流程實際代碼執(zhí)行的是下面七步。第一步初始化根據(jù)雷達首次探測的位置和速度散布生成 N 個粒子N 取1000到2000太多實時性扛不住太少多峰分布表達不出來。第二步預測對每個粒子做RK4積分推進到當前量測時刻并疊加過程噪聲。第三步更新計算每個粒子的量測預測值和雷達實測值之間的新息用高斯似然更新權(quán)重。第四步歸一化計算有效粒子數(shù) N_eff。第五步條件重采樣N_eff 低于閾值時執(zhí)行系統(tǒng)重采樣。第六步輸出狀態(tài)統(tǒng)計量粒子集的加權(quán)均值和協(xié)方差就是當前目標狀態(tài)估計。第七步是軌跡預測的關鍵最后一公里把當前粒子集繼續(xù)向前推進若干秒取每個時刻位置分布的2.5%和97.5%分位數(shù)得到預測航跡的置信走廊。傾側(cè)角 φ 在遞推中是外部參數(shù)濾波器不知道目標當前到底在左傾還是右傾。我一般跑多路并行一路假設 φ 0.3 弧度一路 φ -0.3一路 φ 0各跑一個粒子濾波器最后按每個濾波器近幾幀的平均似然做加權(quán)輸出。這種做法可以理解為工程化的多模型粒子濾波用三簇粒子的組合覆蓋HGRV最常見的三類機動策略。5. 算法實現(xiàn)與參數(shù)調(diào)優(yōu)粒子濾波核心代碼和四個必踩的坑5.1 系統(tǒng)傳播與四階龍格庫塔一段可直接搬的Python函數(shù)整個實現(xiàn)我拆成三塊第一塊是運動模型和數(shù)值積分。下面這段代碼直接定義了HGRV的點質(zhì)量運動方程并封裝了RK4步進函數(shù)。import numpy as np # 工程演示參數(shù)量級參考典型高超聲速滑翔體 G0 9.80665 # 海平面重力加速度m/s^2 RE 6371000.0 # 地球半徑m RHO0 1.225 # 海平面大氣密度kg/m^3 H_SCALE 7200.0 # 密度標高m S_REF 1.2 # 參考面積m^2 MASS 1000.0 # 飛行器質(zhì)量kg CL0 0.45 # 標稱升力系數(shù) L_D 3.0 # 標稱升阻比 def atmosphere_density(alt): # 負高度時按海平面密度處理避免數(shù)值溢出 return RHO0 * np.exp(-np.clip(alt, 0.0, None) / H_SCALE) def motion_derivative(state, phi, lift_factor1.0): HGRV 點質(zhì)量模型導數(shù)。 state [px, py, pz, v, gamma, psi] phi 傾側(cè)角lift_factor 升力系數(shù)乘性修正 px, py, pz, v, gamma, psi state rho atmosphere_density(pz) q 0.5 * rho * v * v L q * S_REF * CL0 * lift_factor / MASS D L / L_D g G0 * (RE / (RE pz)) ** 2 dpx v * np.cos(gamma) * np.cos(psi) dpy v * np.cos(gamma) * np.sin(psi) dpz v * np.sin(gamma) dv -D - g * np.sin(gamma) dgamma L * np.cos(phi) / v - (g - v * v / (RE pz)) * np.cos(gamma) / v dpsi L * np.sin(phi) / (v * np.cos(gamma)) return np.array([dpx, dpy, dpz, dv, dgamma, dpsi]) def rk4_step(state, phi, lift_factor, dt): 四階龍格庫塔單步積分dt 單位秒。 def f(s): return motion_derivative(s, phi, lift_factor) k1 f(state) k2 f(state 0.5 * dt * k1) k3 f(state 0.5 * dt * k2) k4 f(state dt * k3) return state dt / 6.0 * (k1 2 * k2 2 * k3 k4)這段代碼的關鍵在于把氣動模型壓縮成了三個可調(diào)參數(shù)CL0、L_D、lift_factor。前兩個是氣動外形參數(shù)仿真時可以從公開資料量級推算lift_factor 是留給濾波器的修正旋鈕。傾側(cè)角 phi 沒有進入狀態(tài)向量而是作為外部輸入傳入這在第4.3節(jié)說過的多模型并行結(jié)構(gòu)里可以直接復用。RK4步長 dt 我建議取0.1秒配合1赫茲的量測幀率正好每幀10步既能反映機動變化又不會把計算量頂上去。5.2 粒子濾波主循環(huán)預測-更新-重采樣的最小實現(xiàn)第二塊是粒子濾波的單次遞推包含預測、更新、重采樣三個環(huán)節(jié)。測量函數(shù)假設雷達站位于ENU原點如果雷達站不在原點先做目標坐標平移。def measurement_function(state, radar_posnp.zeros(3)): 雷達量測方程返回 [斜距, 方位角, 俯仰角]。 px, py, pz state[0] - radar_pos[0], state[1] - radar_pos[1], state[2] - radar_pos[2] r np.sqrt(px * px py * py pz * pz) if r 1e-6: return np.array([0.0, 0.0, 0.0]) az np.arctan2(py, px) el np.arcsin(pz / r) return np.array([r, az, el]) def angle_diff(a, b): 方位角/俯仰角差值的 ±pi 環(huán)繞處理。 return (a - b np.pi) % (2.0 * np.pi) - np.pi def systematic_resample(particles, weights): 系統(tǒng)重采樣方差低于多項式重采樣適合粒子數(shù)少的場景。 N len(particles) positions (np.arange(N) np.random.uniform(0.0, 1.0)) / N cumulative np.cumsum(weights) new_particles np.empty_like(particles) i, j 0, 0 while i N: if positions[i] cumulative[j]: new_particles[i] particles[j] i 1 else: j min(j 1, N - 1) return new_particles, np.ones(N) / N def pf_predict_update(particles, weights, z_meas, Q, R, phi_est, dt): 粒子濾波單幀遞推。 z_meas [斜距, 方位角, 俯仰角]; phi_est 當前傾側(cè)角假設 N len(particles) R_inv np.linalg.inv(R) # 1) 預測每個粒子獨立傳播并加過程噪聲 for i in range(N): s rk4_step(particles[i], phi_est, 1.0, dt) s s np.random.multivariate_normal(np.zeros(6), Q) s[3] max(s[3], 200.0) # 速度下限避免動壓為負 s[4] np.clip(s[4], -np.pi / 3, np.pi / 3) # 航跡角限幅防止奇點 particles[i] s # 2) 更新計算每個粒子的量測似然 for i in range(N): z_pred measurement_function(particles[i]) innov np.array([ z_meas[0] - z_pred[0], angle_diff(z_meas[1], z_pred[1]), angle_diff(z_meas[2], z_pred[2]), ]) weights[i] * np.exp(-0.5 * innov R_inv innov) # 3) 歸一化與有效粒子數(shù)檢查 w_sum np.sum(weights) if w_sum 1e-12: weights[:] 1.0 / N else: weights / w_sum N_eff 1.0 / np.sum(weights * weights) if N_eff 0.5 * N: particles, weights systematic_resample(particles, weights) state_mean np.sum(particles * weights[:, None], axis0) return particles, weights, state_mean, N_effQ 和 R 矩陣的取值決定了濾波行為的走向這套代碼里我給兩組經(jīng)驗配置。過程噪聲 Q 取 np.diag([502, 502, 302, 202, (0.5度)2, (0.5度)2])位置噪聲幾十米量級速度噪聲20米每秒兩個角度各0.5度對應中等機動工況。如果目標正在橫向機動段把航向角方差加大到2度以上粒子才有能力覆蓋左偏和右偏兩簇軌跡。量測噪聲 R 取 np.diag([802, (0.2度)2, (0.2度)2])和雷達精度匹配。R 給得太小會讓少數(shù)粒子權(quán)重迅速變成1過早退化R 給得太大則量測失去約束力預報區(qū)間的寬度失去意義。5.3 輸出軌跡預報與95%置信區(qū)間分位數(shù)統(tǒng)計代碼第三塊是預報輸出。粒子濾波的價值在預報不在濾波本身。我從當前粒子集中按權(quán)重抽取一部分粒子用同一個運動方程向前推演然后對所有粒子軌跡做逐時刻的分位數(shù)統(tǒng)計。def predict_uncertainty(particles, weights, horizon_s, dt, phi_seqNone): 將當前粒子集向前推進輸出航跡均值與 95% 置信區(qū)間。 phi_seq: 未來每個積分步的傾側(cè)角序列None 時按水平飛行外推 n_use min(200, len(particles)) idx np.random.choice(len(particles), sizen_use, replaceTrue, pweights) steps int(horizon_s / dt) if phi_seq is None: phi_seq np.zeros(steps) traj np.zeros((n_use, steps, 3)) for m, pi in enumerate(idx): s particles[pi].copy() for k in range(steps): s rk4_step(s, phi_seq[k], 1.0, dt) traj[m, k] s[:3] mean_traj traj.mean(axis0) lower np.percentile(traj, 2.5, axis0) upper np.percentile(traj, 97.5, axis0) return mean_traj, lower, upper這段代碼的細節(jié)在于未來操縱假設。如果只預報10到20秒用當前傾側(cè)角外推足夠預報30秒以上就必須給 phi_seq 多樣性。我常用的方法是并行跑三條外推路徑phi 恒為0.3弧度、恒為-0.3弧度、恒為0然后把三組預報區(qū)間合并得到的置信走廊比單一路徑寬但更真實。分位數(shù)用2.5%和97.5%而非標準差是因為粒子集不是高斯分布用標準差會低估不對稱的多峰散布。5.4 四個必踩的坑粒子退化、模型失配、雷達野值和數(shù)值奇點踩坑一粒子退化后濾波器變成“獨苗游戲”?，F(xiàn)象是有效粒子數(shù) N_eff 降到幾百甚至幾十重采樣卻遲遲不觸發(fā)輸出均值開始抖動或漂移。原因通常是重采樣閾值設得太低或者量測噪聲 R 給得太小少數(shù)粒子權(quán)重在幾幀內(nèi)被拉到接近1。解決方法是把重采樣閾值提高到 0.5N并在更新前檢查權(quán)重的最大占比如果某個粒子權(quán)重超過0.8就直接強制重采樣。踩坑二模型失配導致落點預報系統(tǒng)性偏向一側(cè)?，F(xiàn)象是濾波跟隨量測很好但預報落點始終偏同一方向幾十公里。原因多半是標稱氣動參數(shù)和真實目標不一致代碼里 lift_factor 固定為1.0等于用標稱參數(shù)外推。解決方法是把 lift_factor 放進狀態(tài)向量做在線估計或者跑多路粒子濾波每路固定一個不同的 lift_factor取0.8、1.0、1.2三檔按后驗加權(quán)輸出。這個改動能把系統(tǒng)性偏差明顯壓下來。踩坑三雷達野值一幀拉垮全部粒子權(quán)重?，F(xiàn)象是某幀測距值突然跳變幾公里更新之后幾乎所有權(quán)重歸零重采樣后粒子聚集到野值附近。原因是粒子濾波用乘積權(quán)重累積似然野值產(chǎn)生的極小似然會把之前積累的權(quán)重一并抹掉。解決方法是更新前加一個門限判別計算新息歸一化距離 gate innov^T R^{-1} innov如果 gate 超過卡方分布95%分位三維量測約7.8本次跳過更新或給權(quán)重乘一個0.9的衰減因子。這個邏輯雖然簡單卻是工程實現(xiàn)里最容易被忽略的防翻車措施。踩坑四數(shù)值奇點讓積分結(jié)果變成NaN。現(xiàn)象是粒子被過程噪聲推到大航跡角或低速度區(qū)域cos(γ) 接近零航向角導數(shù)爆炸。原因在模型本身ψ? 的分母是 v·cos(γ)粒子漫游到物理不可達區(qū)域時數(shù)值沒有保護。解決方法是傳播后對狀態(tài)做限幅速度不低于200米每秒航跡角限制在±60度高度低于0時強制按0處理更嚴格的做法是傳播后做物理可行性檢查不滿足平衡滑翔包線的粒子直接給零權(quán)重。6. 仿真驗證與評估指標用RMSE、NEES和覆蓋率判斷預測質(zhì)量6.1 三個驗收指標RMSE、NEES與置信區(qū)間覆蓋率仿真驗證我建議跑蒙特卡洛至少100次每次生成一條帶隨機過程噪聲的真值軌跡再疊加雷達量測噪聲喂給粒子濾波器。三個指標是必看的。第一個是位置和速度的RMSE衡量點預測精度。第二個是歸一化估計誤差平方NEES衡量濾波器協(xié)方差是否可信長期小于2說明協(xié)方差給得太大預測區(qū)間寬得沒價值長期大于10說明濾波器過度自信真實誤差經(jīng)常跑出預報走廊。第三個是置信區(qū)間覆蓋率取95%預報區(qū)間統(tǒng)計真實軌跡落在區(qū)間內(nèi)的比例工程上在92%到98%之間都算合理。6.2 我保留的一個驗證習慣離線平滑對照粒子濾波在線只能看到當前和過去的數(shù)據(jù)預報誤差里既有過程噪聲的影響也有濾波收斂慢的影響。我最后再看一個指標——離線固定區(qū)間平滑結(jié)果和在線濾波結(jié)果對比。如果平滑器能明顯修正在線濾波的軌跡說明在線過程噪聲給得太緊或者野值門限沒有生效如果平滑器和在線結(jié)果幾乎一致說明當前參數(shù)已經(jīng)接近這個模型的能力上限。這個對照能在你面對“濾波器似乎不準但不知道哪里不準”時快速定位是模型問題還是實現(xiàn)問題。做完這套驗證這套基于貝葉斯推斷的HGRV軌跡預測方案才算真正可信希望幫到你。本文還有配套的精品資源點擊獲取