盧春光， 周中良， 劉宏強(qiáng)， 寇添， 楊遠(yuǎn)志

(空軍工程大學(xué)航空工程學(xué)院, 西安 710038)

蛇形機(jī)動(dòng)是戰(zhàn)斗機(jī)飛行員在隱蔽接敵、機(jī)動(dòng)規(guī)避及協(xié)同探測(cè)等戰(zhàn)術(shù)動(dòng)作中經(jīng)常采用的戰(zhàn)術(shù)機(jī)動(dòng)類型。在一定的空戰(zhàn)態(tài)勢(shì)下，飛行員進(jìn)行蛇形機(jī)動(dòng)反映了其戰(zhàn)術(shù)意圖，因此在空戰(zhàn)中如何快速而又準(zhǔn)確地識(shí)別出目標(biāo)的蛇形機(jī)動(dòng)模式對(duì)于明確其戰(zhàn)術(shù)意圖以及評(píng)估當(dāng)前的戰(zhàn)場(chǎng)態(tài)勢(shì)具有十分重要的意義。目前，機(jī)動(dòng)模式的識(shí)別主要從機(jī)動(dòng)模式的幾何特征(宏觀)以及運(yùn)動(dòng)參量特征(微觀)2個(gè)方面著手：文獻(xiàn)[1]通過(guò)將戰(zhàn)術(shù)機(jī)動(dòng)進(jìn)行分割并按時(shí)間序列進(jìn)行編碼，采用隱馬爾可夫模型進(jìn)行訓(xùn)練、識(shí)別；文獻(xiàn)[2]通過(guò)提取目標(biāo)航向角變化率、高度變化率等運(yùn)動(dòng)參量特征，設(shè)計(jì)了一個(gè)2級(jí)識(shí)別方案，采用模糊推理和時(shí)間自動(dòng)機(jī)方法實(shí)現(xiàn)了機(jī)動(dòng)模式的識(shí)別。對(duì)于蛇形機(jī)動(dòng)目標(biāo)而言，轉(zhuǎn)彎角速度是其最為關(guān)鍵的運(yùn)動(dòng)參量，因此通過(guò)載機(jī)雷達(dá)的量測(cè)數(shù)據(jù)精確辨識(shí)出轉(zhuǎn)彎角速度對(duì)于提高戰(zhàn)斗機(jī)蛇形機(jī)動(dòng)模式的識(shí)別率具有重要意義。

目前，機(jī)動(dòng)目標(biāo)轉(zhuǎn)彎角速度辨識(shí)的問(wèn)題已經(jīng)得到了許多學(xué)者的關(guān)注。文獻(xiàn)[3-5]將轉(zhuǎn)彎角速度作為目標(biāo)狀態(tài)變量，通過(guò)重新構(gòu)造狀態(tài)方程，分別使用擴(kuò)展卡爾曼濾波(EKF)/無(wú)跡卡爾曼濾波(UKF)、基于5階球徑容積規(guī)則的高階容積卡爾曼濾波(HCKF)等非線性濾波算法估計(jì)角速度和其他狀態(tài)變量，其缺點(diǎn)是角速度估計(jì)的精度依賴于目標(biāo)初始狀態(tài)及協(xié)方差、過(guò)程噪聲等因素。文獻(xiàn)[6]根據(jù)角速度與目標(biāo)速度方向角之間的物理關(guān)系，通過(guò)采用濾波算法對(duì)目標(biāo)速度方向角進(jìn)行估計(jì)，間接求出了角速度的估計(jì)值，該方法的缺點(diǎn)是采用方向角等中間變量間接地估計(jì)角速度，將引入新的估計(jì)誤差。文獻(xiàn)[7]將角速度當(dāng)作系統(tǒng)的未知參數(shù)，采用期望最大化(EM)算法將角速度與目標(biāo)狀態(tài)進(jìn)行聯(lián)合估計(jì)與辨識(shí)，該方法的缺點(diǎn)是使用的是全量測(cè)數(shù)據(jù)，算法運(yùn)行時(shí)間長(zhǎng)，不能滿足空戰(zhàn)中實(shí)時(shí)識(shí)別的要求。

在實(shí)際空戰(zhàn)中，目標(biāo)的角速度是未知且時(shí)變的，并且與目標(biāo)狀態(tài)相互耦合，因此在進(jìn)行聯(lián)合狀態(tài)估計(jì)與角速度辨識(shí)時(shí)，為了獲得角速度辨識(shí)的解析解，需要解除這種耦合關(guān)系。對(duì)在傳統(tǒng)意義上的狀態(tài)估計(jì)與參數(shù)辨識(shí)中，通常認(rèn)為量測(cè)噪聲之間是相互獨(dú)立的，但是在實(shí)際跟蹤過(guò)程中，雷達(dá)的采樣頻率過(guò)高[8]，使得相鄰量測(cè)噪聲之間的相關(guān)性往往不能夠被忽略，因此，在對(duì)角速度辨識(shí)和目標(biāo)狀態(tài)估計(jì)的過(guò)程中必須考慮這種相關(guān)性。針對(duì)上述問(wèn)題，本文基于EM算法框架提出一種帶有色量測(cè)噪聲的聯(lián)合估計(jì)與辨識(shí)算法,并引入滑窗思想，以提高辨識(shí)的實(shí)時(shí)性：E-step利用當(dāng)前估計(jì)的角速度以及經(jīng)過(guò)帶有色量測(cè)噪聲的高階容積卡爾曼平滑(HCKS)獲得的后驗(yàn)平滑概率密度和聯(lián)合概率密度，計(jì)算完備數(shù)據(jù)似然函數(shù)的條件期望；M-step通過(guò)使期望最大化進(jìn)而更新獲得下一次迭代的角速度估計(jì)量。

1 問(wèn)題描述

蛇形機(jī)動(dòng)可以分解成若干個(gè)具有不同轉(zhuǎn)彎角速度的圓弧轉(zhuǎn)彎?rùn)C(jī)動(dòng)，因此蛇形機(jī)動(dòng)目標(biāo)轉(zhuǎn)彎角速度的辨識(shí)問(wèn)題轉(zhuǎn)化成圓弧轉(zhuǎn)彎?rùn)C(jī)動(dòng)的角速度辨識(shí)問(wèn)題。在有色量測(cè)噪聲背景下，假設(shè)目標(biāo)在二維平面中運(yùn)動(dòng)，轉(zhuǎn)彎?rùn)C(jī)動(dòng)模型狀態(tài)方程及量測(cè)方程為

(1)

f(xk,Ωk)=

vk=ψk,k-1vk-1+ζk-1

(2)

由式(1)可知，蛇形機(jī)動(dòng)目標(biāo)的轉(zhuǎn)彎角速度Ωk這一未知參數(shù)耦合在狀態(tài)轉(zhuǎn)移矩陣之中，為了獲取轉(zhuǎn)彎角速度的解析解，需要解除轉(zhuǎn)彎角速度與狀態(tài)轉(zhuǎn)移矩陣之間的非線性耦合關(guān)系，從而將狀態(tài)方程等價(jià)轉(zhuǎn)換成如式(3)形式：

xk+1=F(xk)θ+xk+wk

(3)

式中：

θ1=sin(ΩkT)/Ωk

θ2=(1-cos(ΩkT))/Ωk

θ3=1-cos(ΩkT)

θ4=sin(ΩkT)

由于在有色噪聲條件下，傳統(tǒng)的高斯近似濾波器和平滑器估計(jì)性能不佳，進(jìn)而影響EM算法的辨識(shí)效果，所以為了實(shí)現(xiàn)目標(biāo)狀態(tài)的精確估計(jì)和轉(zhuǎn)彎角速度的準(zhǔn)確辨識(shí)，需要解除相鄰量測(cè)之間相互耦合的關(guān)系，實(shí)現(xiàn)有色量測(cè)噪聲的白化[10]。本文采取量測(cè)差分的方法實(shí)現(xiàn)有色量測(cè)噪聲的白化，重新構(gòu)造量測(cè)方程如下：

(4)

依據(jù)式(1)和式(2)可以得到具有一步狀態(tài)延遲的量測(cè)方程：

(5)

2 基于HCKS-EM的聯(lián)合估計(jì)與辨識(shí)算法

基于極大似然估計(jì)準(zhǔn)則，本文提出了一種帶有色量測(cè)噪聲的聯(lián)合估計(jì)與辨識(shí)算法對(duì)未知角速度進(jìn)行辨識(shí)，具體的算法框架如圖1所示。在第t次迭代時(shí)，通過(guò)將帶有色噪聲的量測(cè)序列經(jīng)過(guò)HCKS，獲得目標(biāo)狀態(tài)的平滑估計(jì)量，并采用EM算法進(jìn)行參數(shù)θ的辨識(shí)，辨識(shí)后的參數(shù)θ用于第t+1次的狀態(tài)估計(jì)，不斷迭代直到滿足設(shè)定的要求為止[12-13]。

圖1 聯(lián)合估計(jì)與辨識(shí)算法Fig.1 Joint estimation and identification algorithm

(6)

(7)

式中：

(8)

dxk-idxk-i-1

(9)

(10)

其中：l為滑動(dòng)窗口的長(zhǎng)度，其長(zhǎng)度大小可依據(jù)實(shí)際需求自行設(shè)定。當(dāng)l=k-1時(shí)，I1是一個(gè)常數(shù)，與待辨識(shí)量θ無(wú)關(guān)；當(dāng)0≤l

pθ(xk-i|xk-i-1)～N(xk-i;F(xk-i-1)θ+xk-i-1,Q)

(11)

ψh(xk-i-1),R)

(12)

2.1 E-step

xk-i-1)TQ-1F(xk-i-1)·

(13)

式中：后驗(yàn)平滑概率密度和基于滑窗內(nèi)數(shù)據(jù)的xk-i和xk-i-1的聯(lián)合概率密度滿足如下分布：

為了求解I2，首先定義如下2種積分函數(shù)：

(14)

(15)

然后分別將F(xi)和xi按xi的分量進(jìn)行分解，可得

Γ(xk-i-1)和Δ(xk-i,xk-i-1)可分別表示為

(16)

(17)

(18)

(19)

則

(20)

(21)

2.2 M-step

(22)

當(dāng)I2取得極大值時(shí)滿足

(23)

則可求得參數(shù)θ的迭代表達(dá)式為

(24)

2.3 帶有色量測(cè)噪聲的HCKS算法的設(shè)計(jì)

眾所周知，平滑是在濾波的基礎(chǔ)之上進(jìn)行的，因此首先設(shè)計(jì)一種帶有色量測(cè)噪聲的HCKF算法，算法如下：

算法1帶有色量測(cè)噪聲的高斯濾波算法框架。

步驟1狀態(tài)預(yù)測(cè)。

(25)

步驟2狀態(tài)修正。

(26)

(27)

帶有色量測(cè)噪聲的HCKS算法是一類前向后向平滑器，其主要包括2個(gè)部分：前向一步平滑和后向固定區(qū)間平滑，具體算法如下：

算法2帶有色量測(cè)噪聲的高斯平滑算法框架。

步驟1前向一步平滑。

(28)

步驟2后向固定區(qū)間平滑。

(29)

2.4 帶有色量測(cè)噪聲的HCKS算法的實(shí)現(xiàn)

本文中采用了Jia等[5]提出的5階球徑容積規(guī)則計(jì)算式(25)、式(27)和式(28)中的非線性高斯積分：

(30)

ξi=

(31)

式中：ei為n維單位向量，且其第i個(gè)元素為1。

(32)

容積點(diǎn)所對(duì)應(yīng)的權(quán)重ωi為

(33)

可將本文算法總結(jié)如下：

算法3基于HCKS-EM的聯(lián)合估計(jì)與辨識(shí)算法

步驟1E-step。已知第t次迭代后參數(shù)θ的估計(jì)值為θt，由式(28)和式(29)可以求得后驗(yàn)平滑概率密度和基于滑窗內(nèi)數(shù)據(jù)的xk-i和xk-i-1的聯(lián)合概率密度，進(jìn)而可以求解得到I2的值。

遞歸：t=0，k←k+1，并將參數(shù)θ的初值設(shè)為上次迭代結(jié)束時(shí)的辨識(shí)值，進(jìn)入下一輪迭代。

3 仿真分析

本文利用水平方向上的轉(zhuǎn)彎?rùn)C(jī)動(dòng)非線性動(dòng)態(tài)模型，仿真出一條蛇形機(jī)動(dòng)軌跡。假設(shè)機(jī)動(dòng)目標(biāo)在1～80 s以Ω1=-2(°)/s作轉(zhuǎn)彎運(yùn)動(dòng)，在k=81 s時(shí)轉(zhuǎn)彎角速度突變?yōu)棣?=2(°)/s并持續(xù)到300 s。設(shè)置轉(zhuǎn)彎角速度初始值為Ω0=-1(°)/s，采樣周期T=1 s，q1=0.1 m2/s3，過(guò)程噪聲wk的協(xié)方差為Q。

Q=diag(q1M,q1M)

通過(guò)載機(jī)雷達(dá)可以獲得目標(biāo)與載機(jī)之間的相對(duì)距離r、方向角φ的信息，則可獲得系統(tǒng)的非線性量測(cè)方程為

基于HCKS-EM聯(lián)合估計(jì)與辨識(shí)算法和擴(kuò)維法的初始狀態(tài)以及協(xié)方差分別為

擴(kuò)維法初始狀態(tài)和協(xié)方差的分別設(shè)置為

100 mrad/s]

為了評(píng)估分析本文聯(lián)合估計(jì)和辨識(shí)算法的性能，設(shè)置了如下3組仿真：

仿真1窗口長(zhǎng)度l分別設(shè)置為5、10、15、20，最大迭代次數(shù)均為5次，各執(zhí)行100次蒙特卡羅仿真。圖2為不同窗口長(zhǎng)度下角速度辨識(shí)結(jié)果，圖3為不同窗口長(zhǎng)度下角速度辨識(shí)均方根誤差(RMSE)。從圖2和圖3可以看出，隨著窗口長(zhǎng)度的增大，該算法在第一階段收斂于真實(shí)值的時(shí)刻就越早，精度越高，并且越穩(wěn)定，但是當(dāng)角速度發(fā)生突變時(shí)，對(duì)于突變的角速度反應(yīng)的就越慢。并且從圖4和表1可以看出，窗口長(zhǎng)度越大，該算法估計(jì)的目標(biāo)狀態(tài)整體精度就越高，但是消耗的時(shí)間越長(zhǎng)，這顯然是時(shí)間與精度之間的“博弈”問(wèn)題，當(dāng)窗口大于10時(shí)，由窗口長(zhǎng)度即量測(cè)數(shù)據(jù)帶來(lái)的精度收益優(yōu)勢(shì)不太明顯，相反運(yùn)算時(shí)間的增加確實(shí)比較可觀的，也就是說(shuō)算法的執(zhí)行效率降低了。并且在角速度發(fā)生突變時(shí)刻及其附近，窗口長(zhǎng)度為10時(shí)的估計(jì)效果明顯好于其他3個(gè)窗口。

圖2 不同窗口長(zhǎng)度下轉(zhuǎn)彎角速度辨識(shí)結(jié)果Fig.2 Turn rate identification results with different window lengths

圖3 不同窗口長(zhǎng)度下轉(zhuǎn)彎角速度辨識(shí)的均方根誤差Fig.3 RMSE of turn rate identification with different window lengths

圖4 不同窗口長(zhǎng)度下位置和速度估計(jì)均方根誤差Fig.4 RMSE of position and velocity estimation with different window lengths

滑動(dòng)窗口長(zhǎng)度時(shí)間/s50.0465100.1048150.1222200.1649

仿真2窗口長(zhǎng)度設(shè)置為10，最大迭代次數(shù)r分別設(shè)置為3、5、10、15，各執(zhí)行100次蒙特卡羅仿真。從圖5和圖6可以看出，隨著迭代次數(shù)的增加，該算法在初始時(shí)刻收斂于真實(shí)值的時(shí)刻就越早，并且對(duì)于角速度突變反應(yīng)的也比較靈敏。從圖7目標(biāo)狀態(tài)4個(gè)分量的RMSE可以看出，迭代次數(shù)越大，該算法估計(jì)的目標(biāo)狀態(tài)整體精度就越高，但是根據(jù)表2中的運(yùn)算時(shí)間可知花費(fèi)的時(shí)間也會(huì)相應(yīng)地增大。尤其是當(dāng)最大迭代次數(shù)大于5時(shí)，最大迭代次數(shù)增加所帶來(lái)的精度收益越發(fā)的不明顯。

圖5 不同最大迭代次數(shù)下轉(zhuǎn)彎角速度辨識(shí)結(jié)果Fig.5 Turn rate identification result with different numbers of maximum iterations

圖6 不同迭代次數(shù)下轉(zhuǎn)彎角速度辨識(shí)的均方根誤差Fig.6 RMSE of turn rate identification with different numbers of iterations

仿真3將本文聯(lián)合估計(jì)與辨識(shí)算法窗口長(zhǎng)度設(shè)置為10，最大迭代次數(shù)設(shè)置為5，將該算法與帶有色的UKF算法和帶有色的HCKF算法各執(zhí)行100次蒙特卡羅仿真。從圖8和圖9可以看出，本文EM算法收斂于角速度真實(shí)值的時(shí)刻比帶有色的UKF算法和帶有色的HCKF算法都要早，并且穩(wěn)定、精度高，但是當(dāng)角速度發(fā)生突變時(shí)，對(duì)于突變的角速度反應(yīng)的速度比較慢。從圖10目標(biāo)狀態(tài)4個(gè)分量的RMSE可以看出，除了在角速度突變時(shí)刻及其鄰近時(shí)刻之外，本文EM算法估計(jì)的精度都要明顯高于帶有色的UKF算法和帶有色的HCKF算法估計(jì)的精度，帶有色的HCKF算法比帶有色的UKF算法估計(jì)的精度要高。

圖7 不同迭代次數(shù)下位置和速度估計(jì)均方根誤差Fig.7 RMSE of position and velocity estimation with different numbers of iterations

最大迭代次數(shù)時(shí)間/s30.049550.0827100.1599150.2577

圖8 EM算法與擴(kuò)維法之間的對(duì)比Fig.8 Comparison between EM algorithm and augmentation method

圖9 轉(zhuǎn)彎角速度辨識(shí)均方根誤差Fig.9 RMSE of turn rate identification

圖10 位置和速度估計(jì)均方根誤差Fig.10 RMSE of position and velocity estimation

4 結(jié) 論

針對(duì)蛇形機(jī)動(dòng)目標(biāo)角速度辨識(shí)與目標(biāo)狀態(tài)聯(lián)合估計(jì)的問(wèn)題，基于EM算法框架提出一種帶有色量測(cè)噪聲的聯(lián)合估計(jì)與辨識(shí)算法。主要結(jié)論有：

1) EM算法設(shè)計(jì)：在E-Step，通過(guò)將完備似然函數(shù)進(jìn)行分解，從而將其轉(zhuǎn)換成帶參數(shù)θ的解析表達(dá)式，提高了辨識(shí)與估計(jì)的精度，降低了由狀態(tài)擴(kuò)維法擴(kuò)維所帶來(lái)的復(fù)雜性；在M-Step，通過(guò)極大化完備似然函數(shù)求得高精度的解析解。

2) 帶有色量測(cè)噪聲的HCKS算法的設(shè)計(jì)：考慮到在實(shí)際空戰(zhàn)中，相鄰噪聲序列之間的相關(guān)性對(duì)于目標(biāo)狀態(tài)估計(jì)與轉(zhuǎn)彎角速度辨識(shí)的影響，設(shè)計(jì)了帶有色量測(cè)噪聲的HCKS算法用于消除這種影響，提高了聯(lián)合辨識(shí)與估計(jì)的精度。