左瑞芹

(中國空空導(dǎo)彈研究院，河南 洛陽 471009)

Kalman 濾波是一種建立在最小均方誤差估計(jì)基礎(chǔ)上的時(shí)域?yàn)V波方法。Kalman 濾波采用狀態(tài)描述系統(tǒng)，其算法采用遞推形式，數(shù)據(jù)存儲量小，便于實(shí)現(xiàn)在線實(shí)時(shí)濾波。自Kalman 濾波器問世以來，其在科學(xué)和工程上已得到了廣泛應(yīng)用。

以雷達(dá)尋的制導(dǎo)系統(tǒng)而論，由于雷達(dá)目標(biāo)的運(yùn)動(dòng)具有隨機(jī)性，故雷達(dá)導(dǎo)引頭輸出的代表導(dǎo)彈-目標(biāo)視線角速度的電信號中，一般都混雜有測量噪聲，通常需要設(shè)法抑制這些噪聲，求取有用信號的最佳估計(jì)值，以削弱其對導(dǎo)引精度的影響。

本文主要介紹了帶有噪聲的導(dǎo)引系統(tǒng)方程，建立了具有Kalman 濾波器的導(dǎo)彈比例導(dǎo)引的系統(tǒng)方程，給出了濾波增益的計(jì)算方法，并對不同條件進(jìn)行了仿真分析，對工程實(shí)現(xiàn)具有一定的參考意義。

1 系統(tǒng)方程

對線性系統(tǒng)而言，Kalman 濾波器能夠提供最優(yōu)估計(jì)。為了說明Kalman 濾波理論，這里用矩陣微分方程的形式來描述系統(tǒng)

其中:x 是描述系統(tǒng)的列向量;F 為系統(tǒng)的動(dòng)態(tài)矩陣;u 為已知的控制向量;w 為白色噪聲。

定義1 個(gè)與噪聲向量相關(guān)的矩陣Q

從另一方面講，Q 為白色噪聲在時(shí)域變換的期望。

如果觀測量與系統(tǒng)狀態(tài)向量線性相關(guān)，則有

其中:z 為觀測向量;H 是觀測矩陣;v 是白色噪聲觀測值。

定義1 個(gè)與噪聲觀測向量v 相關(guān)的噪聲觀測矩陣R

式(1)～式(4)描述的都是連續(xù)系統(tǒng)，對于離散系統(tǒng)需要經(jīng)過1 個(gè)采樣時(shí)間Ts后系統(tǒng)才會更新信息，因此需要把系統(tǒng)模型進(jìn)行離散化。

定義與系統(tǒng)動(dòng)態(tài)矩陣相關(guān)的矩陣為

其中:I 是辨識矩陣;L-1為拉普拉斯逆變換。

式(5)的離散變換矩陣為

式(3)、式(4)的離散形式為:

2 Kalman 濾波器方程

離散形式的Kalman 濾波系統(tǒng)方程以遞推的形式給出:

其中:Kk代表Kalman 增益矩陣，從如下的Ricatti 矩陣方程可計(jì)算得到:

其中:Pk是狀態(tài)估計(jì)更新前的誤差協(xié)方差矩陣;Mk是狀態(tài)估計(jì)更新后的誤差協(xié)方差矩陣;Qk為離散噪聲矩陣。

3 比例導(dǎo)引系統(tǒng)

導(dǎo)彈的末制導(dǎo)，一般采用比例導(dǎo)引規(guī)律引導(dǎo)導(dǎo)彈飛向目標(biāo)，其一般形式為

式中:nc為導(dǎo)彈加速度，m/s2;N'為比例導(dǎo)引常數(shù)，或?qū)Ш奖?Vc為導(dǎo)彈與目標(biāo)的相對速度，m/s;λ 為導(dǎo)彈對目標(biāo)視線與參考線的夾角或視線角，rad;λ˙ 為視線角的導(dǎo)數(shù)，也即視線角速度。

圖1 為在同一平面內(nèi)導(dǎo)彈和目標(biāo)的相對運(yùn)動(dòng)圖。圖1中:VT為目標(biāo)速度，m/s;VM為導(dǎo)彈速度，m/s;nT為目標(biāo)加速度，m/s2;RTM為彈目相對距離;y 為彈目相對高度。

圖1 在同一平面內(nèi)導(dǎo)彈和目標(biāo)的相對運(yùn)動(dòng)

從圖1 中可以看出，彈目沿y 方向的相對加速度為

如果β、λ 為很小，則有

同時(shí)有

假設(shè)彈目相對速度Vc為常數(shù)，則有

定義

其中:tF為導(dǎo)彈總共飛行時(shí)間;t 為當(dāng)前時(shí)間;tgo為截獲時(shí)間。

由式(18)～式(20)可得

對式(20)兩邊求導(dǎo)，可得

通過上面的分析，可以得到比例導(dǎo)引的制導(dǎo)模型框如圖2 所示。

圖2 比例導(dǎo)引律制導(dǎo)模型

假設(shè)已知導(dǎo)彈加速度nc，且目標(biāo)加速度nT是定常的，目標(biāo)隨機(jī)機(jī)動(dòng)，且起動(dòng)時(shí)刻在導(dǎo)彈的攔截工作時(shí)間內(nèi)均勻分布，目標(biāo)機(jī)動(dòng)的隨機(jī)性可以看成通過1 個(gè)積分器的白色噪聲過程，其噪聲的譜密度為［6］

圖2 中的模型可以表示為

把式(24)中F 的值帶入式(5)，可得

寫成離散形式為

由式(9)、式(24)和式(26)可得離散形式的G 為

由式(7)可取系統(tǒng)離散形式的觀測方程

由式(9)、式(26)～式(28)可得圖2 的帶Kalman 濾波器的形式，表示成矩陣為:

由式(15)和式(22)可得與狀態(tài)估值相關(guān)的制導(dǎo)指令為

由式(18)、式(29)和式(30)可知，系統(tǒng)可用圖3 所示的框圖形式表示。圖3 中為引入噪聲之后的視線角;z-1為純延遲;z-1yk表示yk-1。圖3 中Kalman 濾波器給出了相對位置、相對速度和目標(biāo)加速度的最優(yōu)估計(jì)。

圖3 加入了Kalman 濾波器制導(dǎo)系統(tǒng)

4 Kalman 增益的計(jì)算

圖3 所示的系統(tǒng)需要先計(jì)算濾波器的增益Kk，這些增益可通過一系列Ricatti 方程計(jì)算得到。對于圖2 所示的三狀態(tài)系統(tǒng)，第1 個(gè)Ricatti 方程Mk可通過代入式(2)、式(14)和式(26)后整理得到

第2 個(gè)Ricatti 方程用來計(jì)算Kalman 增益。對于圖2 來說，Rk是1 ×1，因此有:

那么，第3 個(gè)Ricatti 方程的解為

為了求解Ricatti 方程，式(31)需要初始化系數(shù)矩陣P0。1 種經(jīng)驗(yàn)的初始化方法為

5 仿真實(shí)例

采用圖2 的制導(dǎo)系統(tǒng)模型和圖3 的濾波框圖，取nT=3 g，Vc=1 200 m/s，VM=850 m/s，VT=350 m/s，σnoise=0.001，HE=20°，tF=10 s，ts=0.1 s 來進(jìn)行制導(dǎo)系統(tǒng)仿真。圖4所示為取不同σnoise時(shí)的加速度估計(jì)值曲線，圖5 為加速度估計(jì)誤差的曲線。

從圖4 可以看出，σnoise取值較小時(shí)，系統(tǒng)估計(jì)值的上升時(shí)間短。從圖5 中可看出，加速度估計(jì)誤差很快進(jìn)入了σ 的理論范圍之內(nèi)。

圖4 不同σnoise的估計(jì)值比較

圖5 加速度誤差

6 結(jié)束語

本文主要介紹了Kalman 濾波器及其在導(dǎo)彈自動(dòng)導(dǎo)引中的應(yīng)用，給出了Kalman 濾波器方程及其增益計(jì)算方法，推導(dǎo)了Kalman 濾波器在比例導(dǎo)引中的工程應(yīng)用方法，對給定的飛行條件進(jìn)行了仿真分析，仿真結(jié)果證明了方法的有效性。

［1］Yuan P J，Hsu S C.Rendezvous Guidance with Proportional Navigation［J］. Guidance，Control and Dynamics，1994，17(2):409-411.

［2］Ludeman L C. Filtering，Estimation and Detection，Wiley［M］.Random processes，2003.

［3］Becker K.Closed-Form Solution of Pure Proportional Navigation［J］.Transnactions on Aerospace and Electronic Systems，1990，20(3):526-533.

［4］Ghose D.On the generalization of TPN［J］.Transaction on Aerospace and Electronic Systems，1994，30(2):545-555.

［5］Lee G T，Lee J G.Improved Command to Line-of-Sight for Homing Guidance［J］.Transactions on Aerospace and Electronic Systems，1995，31(1):506-510.

［6］錢輝，丁永忠.大航程AUVSINS/DVL 組合導(dǎo)航定位精度研究［J］.兵工自動(dòng)化，2010(2):46-48.