[發(fā)明專利]基于輸出相關(guān)自適應(yīng)卡爾曼濾波的彈道軌跡形成方法

有效

• 申請(qǐng)?zhí)枺?/td>: 201310722729.5
• 申請(qǐng)日: 2013-12-24
• 公開（公告）號(hào): CN103744057A
• 公開（公告）日: 2014-04-23
• 發(fā)明（設(shè)計(jì)）人: 歐陽廣帥;胡居榮;周晶;肖豐;曹寧
• 申請(qǐng)（專利權(quán)）人: 河海大學(xué)
• 主分類號(hào): G01S7/02
• 分類號(hào): G01S7/02
• 代理公司: 南京經(jīng)緯專利商標(biāo)代理有限公司 32200
• 代理人: 朱小兵

1.一種基于輸出相關(guān)自適應(yīng)卡爾曼濾波的彈道軌跡形成方法，其特征在于，具體步驟如下：

步驟1，根據(jù)彈道模型建立基于狀態(tài)變量的彈道方程及觀測方程，將狀態(tài)變量與標(biāo)稱狀態(tài)作差得到狀態(tài)偏差，從而建立基于狀態(tài)偏差的近似線性方程；

步驟2，由觀測數(shù)據(jù)計(jì)算出輸出相關(guān)函數(shù)序列；

步驟3，由輸出相關(guān)函數(shù)序列推算出卡爾曼濾波增益矩陣的最優(yōu)穩(wěn)態(tài)解；

步驟4，將卡爾曼濾波增益矩陣的最優(yōu)穩(wěn)態(tài)解代入卡爾曼濾波算法遞推方程中的狀態(tài)估計(jì)方程，得到輸出相關(guān)自適應(yīng)濾波的卡爾曼濾波算法遞推方程，對(duì)雷達(dá)觀測的一段飛行彈道參數(shù)進(jìn)行濾波，實(shí)現(xiàn)對(duì)彈道數(shù)據(jù)的自適應(yīng)卡爾曼濾波估計(jì)和修正；

步驟5，根據(jù)步驟4濾波后的飛行彈道參數(shù)，利用彈道模型外推彈道落點(diǎn)，形成彈道軌跡。

2.根據(jù)權(quán)利要求1所述的一種基于輸出相關(guān)自適應(yīng)卡爾曼濾波的彈道軌跡形成方法，其特征在于，步驟1中所述建立基于狀態(tài)偏差的近似線性方程具體如下：

在地面坐標(biāo)系中建立彈道方程：

dXdt=dxdtdydtdzdtdvxdtdvydtdvzdt=vxvyvz-KDSCx(Ma)2mρvr(vx-wx)+KDSCx(Ma)2mρvrwx′-KDSCx(Ma)2mρvrvy-g-KDSCx(Ma)2mρvr(vz-wz)+KLBZ+KDSCx(Ma)2mρvrwz′]]>

式中：X為t時(shí)刻的狀態(tài)變量；表示對(duì)時(shí)間t求導(dǎo)數(shù)；x，y，z分別為彈丸的距離、高度、側(cè)偏；v_x,v_y,v_z分別為彈丸速度的三軸分量；v_r為含風(fēng)速在內(nèi)的彈丸實(shí)際速度；w_x,w_z分別為縱風(fēng)和橫風(fēng)的風(fēng)速，在標(biāo)準(zhǔn)條件下均取為0；w'_x,w'_z分別為縱風(fēng)和橫風(fēng)的隨機(jī)量；S為彈丸最大截?cái)嗝娣e，d為彈徑；m為彈體質(zhì)量；Ma為馬赫數(shù)，C_s為聲速，其中k為比熱比，對(duì)于空氣k＝1.404；R_d為干空氣氣體常數(shù)；τ為虛溫，是把濕空氣折合成干空氣時(shí)對(duì)氣溫的修正；C_x(Ma)為空氣阻力系數(shù)，它是馬赫數(shù)Ma的函數(shù)；ρ為空氣密度，其中，p為氣體壓強(qiáng)；g為海拔y處的重力加速度，其中g(shù)₀為海平面重力加速度，r₀為地球半徑；K_D為阻力符合系數(shù)；B_Z為側(cè)向升力加速度；K_L為由動(dòng)力平衡角引起的側(cè)向加速度時(shí)彌補(bǔ)誤差的修正系數(shù)，又稱靜力矩符合系數(shù)；

令狀態(tài)變量X＝[x,y,z,v_x,v_y,v_z,c,ac]^T＝[x₁,x₂,x₃,x₄,x₅,x₆,x₇,x₈]^T，式中c為彈道系數(shù)，ac為修正的彈道系數(shù)，x₁～x₈為狀態(tài)變量X的8個(gè)狀態(tài)，分別與x,y,z,v_x,v_y,v_z,c,ac對(duì)應(yīng)；

令x_4r＝x₄-w_x表示彈丸在x方向受風(fēng)速影響后的實(shí)際速度；x_6r＝x₆-w_z，表示彈丸在，z方向的實(shí)際飛行速度，則

取作為中間變量；同時(shí)將聲速帶入，馬赫數(shù)表示為Ma=vrCs=vr20.0468τ,]]>將彈道方程改寫為：

dXdt=f(X)+ΓW=x4x5x6-x7Nx4r-x7Nx5-g-x7Nx6r+x8BZ00+03×3I3×302×3x7Nwx′0wz′]]>

式中，f(X)為狀態(tài)變量X的函數(shù)；Γ為噪聲驅(qū)動(dòng)陣，Γ＝[0_3×3?I_3×3?0_2×3]^Tx₇N；W為系統(tǒng)噪聲，W＝[w'_x?0?w'_z]^T；0_3×3表示3×3的零矩陣；I_3×3表示3×3的單位陣；0_2×3表示2×3的零矩陣；

建立觀測方程為：

Z＝h(X)+V

式中，h(X)=[x2+y2+z2arctanzxarctanyx2+z2]T,]]>V為雷達(dá)量測噪聲；

在采樣區(qū)間[t_k,t_k+1]上將f(X(t))近似展開并將彈道方程離散化，得到：

X_k+1＝G(X_k)+Γ_kq_k

Z_k+1＝h(X_k+1)+V_k+1

式中，X_k為k時(shí)刻的狀態(tài)變量；X_k+1為k+1時(shí)刻的狀態(tài)變量；G(Xk)=Xk+f(Xk)T+A(Xk)f(Xk)T22,]]>A(X)是f(X)關(guān)于X的偏導(dǎo)，即A(X)=∂f(X)∂X;]]>T為時(shí)間間隔，T＝t_k+1-t_k；模型噪聲Z_k+1為k+1時(shí)刻的觀測變量；V_k+1為k+1時(shí)刻的雷達(dá)量測噪聲；

若不考慮系統(tǒng)噪聲，則系統(tǒng)k時(shí)刻的狀態(tài)變量、觀測變量的標(biāo)稱狀態(tài)為

XK+1*=G(Xk)]]>

Zk+1*=h(Xk+1*)+Vk+1]]>

將狀態(tài)變量與標(biāo)稱狀態(tài)作差得到狀態(tài)偏差，則：

ΔX=Xk-Xk*]]>

ΔZ=Zk-Zk*]]>

圍繞標(biāo)稱狀態(tài)將函數(shù)h(X_k+1)和G(X_k)泰勒展開，并取一階近似，得到基于狀態(tài)偏差的近似線性方程為：

ΔXk+1=∂G(Xk*)∂Xk*ΔXk+Γ(Xk*)qk]]>

ΔZk+1=∂h(Xk*)∂Xk*ΔXk+1+Vk+1.]]>