胡浩然,陳樹新,吳昊,2,何仁珂,吳強(qiáng),張喜慶

(1.空軍工程大學(xué)信息與導(dǎo)航學(xué)院,710077,西安;2.地理信息工程國家重點(diǎn)實(shí)驗室,710054,西安;3.復(fù)雜航空系統(tǒng)仿真重點(diǎn)實(shí)驗室,100076,北京;4.解放軍95894部隊,100076,北京;5.華恒智禾電子信息技術(shù)有限公司,710077,西安)

在目標(biāo)跟蹤[1-2]領(lǐng)域,離散-離散時間系統(tǒng)將目標(biāo)運(yùn)動模型及其測量模型均視為離散模型。但是,考慮目標(biāo)運(yùn)動的實(shí)際情況,運(yùn)動模型在本質(zhì)上應(yīng)是連續(xù)的,而其有間隔的測量是離散的,這樣的系統(tǒng)稱為連續(xù)-離散時間系統(tǒng)[3]。相比較于傳統(tǒng)跟蹤算法常采用的離散-離散系統(tǒng)模型,連續(xù)-離散系統(tǒng)的數(shù)學(xué)解算方式精度更高,在目標(biāo)跟蹤領(lǐng)域更具優(yōu)勢[4]。

文獻(xiàn)[4]將連續(xù)-離散時間系統(tǒng)與容積準(zhǔn)則相結(jié)合,形成了連續(xù)-離散容積卡爾曼濾波(continuous-discrete cubature Kalman filter,CD-CKF)算法。該算法以隨機(jī)微分方程(stochastic differential equation,SDE)描述連續(xù)時間目標(biāo)運(yùn)動模型,而SDE的求解精度又很大程度上影響著濾波的性能。文獻(xiàn)[5]提出了基于1.5階的It-Taylor算法,有效解決了連續(xù)時間的預(yù)測問題,但同時存在精度上限不高的缺點(diǎn)。為進(jìn)一步提高狀態(tài)估計的精度,文獻(xiàn)[6]將基于SDE的連續(xù)模型描述為狀態(tài)期望和協(xié)方差的形式,進(jìn)而用高階數(shù)值近似算法對連續(xù)模型進(jìn)行了求解。文獻(xiàn)[7]采用了一種嵌套式的隱式高階龍格-庫塔算法,解決了連續(xù)-離散擴(kuò)展卡爾曼濾波(continuous-discrete extended Kalman filter,CD-EKF)的狀態(tài)估計問題,該算法也成功應(yīng)用于CD-CKF算法,有效提高了濾波的精度[8]。文獻(xiàn)[9]采用自適應(yīng)步長對上述數(shù)值近似算法完成了優(yōu)化,進(jìn)一步提高了連續(xù)-離散濾波算法的性能和效率。

以上提高濾波精度的算法都是在假設(shè)測量噪聲為高斯型的理想條件下建立的。但是,在實(shí)際測量環(huán)境中,信號干擾等惡劣條件會導(dǎo)致非高斯噪聲或異常誤差的出現(xiàn)。如果不加以處理,濾波就會不準(zhǔn)確甚至發(fā)散。對于此類測量不確定問題,有關(guān)學(xué)者在離散-離散時間系統(tǒng)做了充分的研究,提出了多種魯棒濾波算法。文獻(xiàn)[10]采用Huber代價函數(shù)來降低異常數(shù)據(jù)的權(quán)重;Wu等[11]提出了基于廣義M估計的魯棒CKF算法,無需考慮降權(quán)函數(shù)的經(jīng)驗閾值和異常誤差的統(tǒng)計特征信息。近年來,相關(guān)熵準(zhǔn)則[12]因能夠有效抑制非高斯噪聲而受到廣泛關(guān)注。作為一種局部相似性度量,其與僅使用測量量二階矩信息的經(jīng)典最小均方誤差準(zhǔn)則相比,能夠充分利用測量量的更高階信息。因此,相關(guān)熵理論適用于受非高斯噪聲污染的隨機(jī)系統(tǒng)。最大相關(guān)熵準(zhǔn)則(maximum correntropy criterion,MCC)已成功應(yīng)用于離散-離散時間系統(tǒng)濾波中[13-14],提高了算法的魯棒性。

為了進(jìn)一步提升目標(biāo)跟蹤精度和魯棒性,本文將改進(jìn)的MCC代價函數(shù)(傳統(tǒng)的MCC濾波算法存在數(shù)值問題)擴(kuò)展到連續(xù)-離散時間系統(tǒng),同時引入自適應(yīng)因子和平方根濾波技術(shù),提出了一種具有魯棒機(jī)制的連續(xù)-離散目標(biāo)跟蹤算法。該算法使用連續(xù)-離散模型描述目標(biāo)運(yùn)動,利用所提MCC算法處理測量噪聲,由不同測量環(huán)境確定自適應(yīng)因子的取值,并通過自適應(yīng)因子對觀測噪聲的協(xié)方差矩陣在線調(diào)整,進(jìn)而完成對狀態(tài)的估計以及協(xié)方差的實(shí)時更新。同時,為保證誤差協(xié)方差矩陣的正定性和對稱性,將其與平方根濾波技術(shù)相結(jié)合,從而得到了魯棒平方根連續(xù)-離散自適應(yīng)最大相關(guān)熵容積卡爾曼濾波(robust square-root continuous-discrete adaptive maximum correntropy cubature Kalman filter,RSRCD-AMCCKF)算法。仿真驗證結(jié)果表明,本文算法能夠有效處理測量中可能出現(xiàn)的高斯、非高斯噪聲以及測量突發(fā)異常的情況。

1 連續(xù)-離散系統(tǒng)目標(biāo)跟蹤模型和平方根連續(xù)-離散容積卡爾曼濾波算法

1.1 連續(xù)離散系統(tǒng)目標(biāo)跟蹤模型

在連續(xù)-離散目標(biāo)跟蹤系統(tǒng)中,目標(biāo)運(yùn)動模型為連續(xù)時間模型,測量模型是離散時間模型?？紤]到隨機(jī)擾動,連續(xù)時間目標(biāo)運(yùn)動模型可以描述為隨機(jī)微分方程的形式

(1)

式中:x(t)是n維的狀態(tài)向量;f為狀態(tài)轉(zhuǎn)移函數(shù);w(t)是n維標(biāo)準(zhǔn)布朗運(yùn)動;Q為布朗運(yùn)動的協(xié)方差矩陣。

依據(jù)狀態(tài)期望和協(xié)方差進(jìn)行狀態(tài)估計有利于提高連續(xù)-離散濾波算法的精度,期望和協(xié)方差可表示為

(2)

(3)

離散時間測量模型為

zk=h(xk,k)+vk

(4)

式中:zk是實(shí)際觀測值;h為觀測函數(shù);vk是觀測噪聲,指代高斯或非高斯噪聲;k是指離散的時間測量點(diǎn)。

1.2 平方根連續(xù)離散容積卡爾曼濾波算法

1.2.1 時間更新

在平方根連續(xù)-離散容積卡爾曼濾波(square-root continuous-discrete cubature Kalman filter,SRCD-CKF)算法中,狀態(tài)容積點(diǎn)可定義為

(5)

式中:ξi是容積點(diǎn)集;S(t)為協(xié)方差矩陣P(t)的下三角矩陣,滿足P(t)=S(t)ST(t)。

基于容積準(zhǔn)則,期望和協(xié)方差可重新表示為

(6)

(7)

式中ε和W是容積參數(shù),定義為

(8)

其中,I2n表示維數(shù)為2n的單位矩陣,1是單位列向量,?代表Kronecker積。

對協(xié)方差傳播的計算,本文采用文獻(xiàn)[15]提出的平方根協(xié)方差計算算法,將其與容積準(zhǔn)則結(jié)合,并用高階數(shù)值近似算法進(jìn)行求解

(9)

(10)

式中:Bi,j是Φ(B)的第i行、第j列的元素;Φ(B)作為判斷矩陣,完成平方根形式的誤差協(xié)方差矩陣更新。B(t)的定義為

B(t)=S-1(t)[X(t)WFT(X(t))+

F(X(t))XT(t)+Q(t)]S-T(t)

(11)

(2)協(xié)方差分解,計算式為

P(tk)=S(tk)ST(tk)

(12)

(3)計算狀態(tài)容積點(diǎn),公式為

(13)

(4)狀態(tài)容積點(diǎn)傳播,計算式為

(14)

(5)計算狀態(tài)預(yù)測值,公式為

(15)

(6)求解預(yù)測平方根協(xié)方差,計算式為

(16)

1.2.2 測量更新

測量更新流程如下。

(1)計算狀態(tài)容積點(diǎn),公式為

Sk|k-1=S(tk)

(17)

(18)

(2)容積點(diǎn)的測量傳播,計算式為

Zi,k|k-1=h(Xi,k|k-1,k)

(19)

(3)計算測量預(yù)測值,公式為

(20)

(4)構(gòu)建測量加權(quán)中心矩陣,計算式為

(21)

(5)計算新息協(xié)方差矩陣,公式為

(22)

(6)構(gòu)建狀態(tài)加權(quán)中心矩陣,計算式為

(23)

(7)計算交叉協(xié)方差矩陣,公式為

(24)

(8)計算連續(xù)-離散的容積增益,公式為

(25)

(9)計算狀態(tài)估計值,公式為

(26)

(10)更新協(xié)方差矩陣,計算式為

(27)

2 基于連續(xù)-離散系統(tǒng)的自適應(yīng)最大相關(guān)熵平方根容積卡爾曼濾波算法

2.1 最大相關(guān)熵準(zhǔn)則

相關(guān)熵是用于衡量兩個隨機(jī)變量之間相似性的一種新方法。對于給定的兩個隨機(jī)變量X和Y,相關(guān)熵定義為[12]

V(X,Y)=E[κ(X,Y)]=?κ(x,y)dFX,Y(x,y)

(28)

式中:E[·]為數(shù)學(xué)期望;FX,Y(x,y)為聯(lián)合概率密度函數(shù);κ(·)為梅瑟型正定核函數(shù)。核函數(shù)的選擇決定了相關(guān)熵的求解方式,本文選取高斯核作為相關(guān)熵的核函數(shù),其定義為

(29)

式中e=x-y,σ>0是核函數(shù)的核帶寬。

實(shí)際中,聯(lián)合概率密度函數(shù)通常是難以獲取的,且所獲得的樣本數(shù)據(jù)采樣點(diǎn)也是有限的。因此,相關(guān)熵可以通過T個樣本數(shù)據(jù)采樣點(diǎn)進(jìn)行估計

(30)

(31)

由式(31)可以看出,隨機(jī)相關(guān)熵信息包含了變量誤差X-Y的所有偶數(shù)階矩。因此,只要選擇合適的核帶寬,就可以獲取更高階的信息。對受高斯噪聲污染的系統(tǒng),與使用測量值二階矩信息的最小均方誤差準(zhǔn)則相比,相關(guān)熵是完全可以用來描述其統(tǒng)計信息的。對于測量中出現(xiàn)的非高斯噪聲,最大相關(guān)熵準(zhǔn)則能夠利用測量量的更高階信息,對處理非高斯隨機(jī)系統(tǒng)也是十分有效的。

2.2 魯棒平方根連續(xù)離散自適應(yīng)最大相關(guān)熵容積卡爾曼濾波算法

在連續(xù)時間狀態(tài)模型中,其過程噪聲是高斯形式的,但離散的測量噪聲是非合作的。因此,本文將加權(quán)最小二乘(weighted least squares,WLS)代價函數(shù)與最大相關(guān)熵準(zhǔn)則相結(jié)合,定義新的相關(guān)熵代價函數(shù)

(32)

(33)

(34)

考慮式(32)所示的代價函數(shù),濾波算法的最終目的是獲得xk的最優(yōu)估計,即

(35)

式(35)的最優(yōu)解可通過求解極值獲得

(36)

式中自適應(yīng)因子Lk的表達(dá)式為

(37)

這里由最大相關(guān)熵準(zhǔn)則得到的自適應(yīng)因子Lk是一個標(biāo)量,可通過式(29)完成計算,它很容易分離并且避免了零矩陣求逆的數(shù)值計算問題。

為了保證當(dāng)σ→∞時,此算法能夠收斂于傳統(tǒng)的CD-CKF算法,令γ=1、β=-2σ2,則式(36)可重新表述為

(38)

兩邊移項,提取公共項可以得到

(39)

由此獲得的方程式(39)是一個不動點(diǎn)方程,因為Lk與xk相關(guān),可以通過固定點(diǎn)迭代算法求解。經(jīng)一次迭代后即可得到所需的狀態(tài)更新值與濾波增益

(40)

(41)

相關(guān)的誤差協(xié)方差矩陣的更新計算式為

(42)

(43)

(1)當(dāng)測量值發(fā)生較大突變或受噪聲影響較大,即‖zk‖∞→∞時,傳統(tǒng)的相關(guān)熵濾波算法[16]會出現(xiàn)零矩陣的求逆計算,這會導(dǎo)致原有的濾波算法失效,而本文提出的RSRCD-AMCCKF算法則可以避免上述數(shù)值計算問題。

(2)高斯核帶寬σ在很大程度上影響濾波算法的精度。已有研究表明,σ較小時對非高斯測量噪聲具有更強(qiáng)的魯棒性[16]。但是,當(dāng)σ過小時,算法的濾波性能會退化較為嚴(yán)重,甚至?xí)霈F(xiàn)濾波發(fā)散的情況。針對連續(xù)-離散時間系統(tǒng)的實(shí)際情況,根據(jù)σ的取值對自適應(yīng)因子進(jìn)行調(diào)整,以便更加合理地處理不同測量情況也是十分重要的。

3 仿真與分析

仿真中考慮典型的空中交通管制場景,跟蹤目標(biāo)在水平面內(nèi)以一定的轉(zhuǎn)彎率做協(xié)調(diào)轉(zhuǎn)彎運(yùn)動[17-18],分析本文提出的目標(biāo)跟蹤算法在測量噪聲為高斯噪聲、測量噪聲為非高斯噪聲、測量值發(fā)生突變這3種場景下的性能變化。目標(biāo)的協(xié)調(diào)轉(zhuǎn)彎運(yùn)動是一種典型的非線性運(yùn)動,其連續(xù)時間運(yùn)動模型可以用隨機(jī)微分方程表示。

采用CKF算法[18]、MCCKF算法、SRCD-CKF算法、基于Huber兩段權(quán)函數(shù)的連續(xù)-離散魯棒(Huber robust square-root continuous-discrete cubature Kalman filter,HRSRCD-CKF)算法、基于P-Huber三段權(quán)函數(shù)的連續(xù)-離散魯棒(M-estimation based robust square-root continuous-discrete cubature Kalman filter,MRSRCD-CKF)算法作為對比算法。HRSRCD-CKF算法和MRSRCD-CKF算法可以理解為Huber魯棒函數(shù)在本文提出的連續(xù)-離散系統(tǒng)框架中的應(yīng)用,其較離散-離散系統(tǒng)的濾波精度更高[19]。MCCKF算法采用本文提出的最大相關(guān)熵算法,σ取值參考文獻(xiàn)[20],HRSRCD-CKF算法中異常誤差的判別門限值取為1.345,而MRSRCD-CKF算法中的α分位點(diǎn)取為α1=0.1,α2=0.1%。

離散時間的非線性測量模型為

zk=arctan(y(tk)/x(tk))+vk

(44)

式中vk為測量噪聲,具體表達(dá)形式視仿真場景進(jìn)行定義。

克拉美羅下界是目標(biāo)參數(shù)跟蹤誤差的下界,它提供了濾波算法理論達(dá)到的最好均方誤差?？死懒_下界在k時刻的定義為

(45)

文獻(xiàn)[21]已證明,Fisher信息矩陣Mk的遞歸計算式為

(46)

考慮到本文采用純方位跟蹤系統(tǒng),Mk在該系統(tǒng)中的表達(dá)式[21]為

(47)

式中

(48)

ψk是觀測函數(shù)的雅可比矩陣,其具體定義為

(49)

為了評估各算法的性能,分別定義位置量、速度量的均方根誤差(root mean square error,RMSE),計算式為

(50)

(51)

(1)仿真場景1,測量噪聲為高斯噪聲。高斯條件下的濾波問題是目標(biāo)跟蹤算法的基礎(chǔ),目標(biāo)跟蹤過程常常討論該條件下濾波算法的精度問題[22-23]。測量噪聲vk滿足0均值的高斯分布,即vk～N(0,Rk),本次仿真中Rk的初始值為一維常值0.001。當(dāng)高斯核帶寬σ取不同值時,所提RSRCD-AMCCKF算法的eRMS,x和eRMS,v如表1所示。

由表1可知:當(dāng)高斯核帶寬σ不同時,RSRCD-AMCCKF算法表現(xiàn)出了不同的性能,σ過小或過大都會降低濾波的精度;當(dāng)σ=0.1時,濾波發(fā)散;當(dāng)σ越來越大時,RSRCD-AMCCKF濾波性能逐漸趨近于SRCD-CKF算法,符合之前的理論推導(dǎo);當(dāng)σ=50時,RSRCD-AMCCKF濾波精度仍高于SRCD-CKF算法。

表1 高斯條件下不同核帶寬時的算法精度

高斯核帶寬σ=5時,各算法的eRMS,x和eRMS,v分別如圖1、圖2所示。

圖1 高斯噪聲下各算法的位置均方根誤差比較Fig.1 Comparison of position root mean square error of each algorithm under Gaussian noise

圖2 高斯噪聲下各算法的速度均方根誤差比較Fig.2 Comparison of velocity root mean square error of each algorithm under Gaussian noise

由圖1和圖2可知:當(dāng)測量噪聲為高斯噪聲時,各算法前期走勢大致相同,都能夠隨著測量次數(shù)的增加趨于收斂,其中,連續(xù)-離散濾波算法較離散-離散濾波算法表現(xiàn)出了更高的精度;MCCKF算法較CKF算法略有優(yōu)勢,其精度比較接近連續(xù)-離散算法,HRSRCD-CKF和MRSRCD-CKF較原始的SRCD-CKF算法有小幅波動,在測量正常的情況下精度基本持平;RSRCD-AMCCKF相比較于以上算法最接近于CRLB,濾波性能更優(yōu),這是因為最大相關(guān)熵準(zhǔn)則能夠利用測量量的更高階信息,將高斯核應(yīng)用于估計誤差矩陣的每個元素,與SRCD-CKF的結(jié)合效果更好。

(2)仿真場景2,測量噪聲為非高斯噪聲。在目標(biāo)跟蹤過程中,由于復(fù)雜目標(biāo)不同部位的散射強(qiáng)度和相對相位的隨機(jī)變化,造成回波相位波前面的畸變,波前在接收天線口徑面上的傾斜和隨機(jī)擺動必然引起測量誤差,這種現(xiàn)象引起的測量噪聲稱之為閃爍噪聲。閃爍噪聲分布與高斯分布的主要差別在于尾部較長,而在中心區(qū)域則類似于高斯形狀,常表現(xiàn)出非高斯特征。目前,針對被污染目標(biāo)的這一特性,常采用具有不同方差的高斯噪聲加權(quán)對閃爍噪聲建模,以此來表示測量中出現(xiàn)的非高斯誤差[24],即測量噪聲vk為高斯混合噪聲,vk～0.99N(0,Rk)+0.01N(0,500Rk)。當(dāng)σ不同時,RSRCD-AMCCKF算法在非高斯條件下的eRMS,x和eRMS,v如表2所示。

表2 非高斯條件下不同核帶寬時的算法精度

由表2可知:能夠?qū)Ψ歉咚乖肼暺鸬揭种谱饔玫摩野l(fā)生了變化,相比于高斯條件下RSRCD-AMCCKF算法的濾波精度有所降低,但合適的σ仍能對測量中出現(xiàn)的非高斯噪聲起到一定的抑制作用。因此,根據(jù)不同測量環(huán)境選擇是比較重要的。

在非高斯條件下,σ=5時各算法的eRMS,x和eRMS,v分別如圖3和圖4所示。

圖3 非高斯噪聲下各算法的位置均方根誤差比較Fig.3 Comparison of position root mean square error of each algorithm under non-Gaussian noise

圖4 非高斯噪聲下各算法的速度均方根誤差比較Fig.4 Comparison of velocity root mean square error of each algorithm under non-Gaussian noise

由圖3和圖4可知:當(dāng)測量噪聲為非高斯噪聲時,各算法的誤差均出現(xiàn)了不同程度的增大,但在濾波精度上連續(xù)-離散算法仍表現(xiàn)出了相對較高的精度;HRSRCD-CKF算法和MRSRCD-CKF算法不能夠?qū)Ψ歉咚箿y量噪聲起到抑制作用,算法基本失效;MCCKF算法具備一定程度的抑制效果,但穩(wěn)定性不高,這可能與核帶寬的取值有關(guān);RSRCD-AMCCKF算法能夠?qū)Ψ歉咚乖肼暺鸬捷^好的抑制作用,最終表現(xiàn)出相對較高的收斂精度,體現(xiàn)了MCC算法在處理非高斯噪聲問題上所具有的優(yōu)勢。

圖5 測量異常時各算法的位置均方根誤差比較Fig.5 Comparison of position root mean square error of each algorithm under abnormal measurement

圖6 測量異常時各算法的速度均方根誤差比較Fig.6 Comparison of velocity root mean square error of each algorithm under abnormal measurement

由圖5和圖6可知:當(dāng)測量發(fā)生異常時,不具備魯棒性的CKF算法和SRCD-CKF算法誤差陡然增加,CKF算法較SRCD-CKF算法精度稍低,而其他算法均有效抵制了異常測量。RSRCD-AMCCKF算法在受到干擾時表現(xiàn)出了更好的魯棒性,對異常測量抑制的效果最為明顯。HRSRCD-CKF算法、MRSRCD-CKF算法和MCCKF算法抑制效果稍差。同時,以連續(xù)-離散算法跟蹤目標(biāo)運(yùn)動的RSRCD-AMCCKF算法較MCCKF算法擁有更高的精度,其收斂精度更趨向于CRLB曲線,說明本文提出的RSRCD-AMCCKF算法能夠較好地實(shí)現(xiàn)濾波精度和魯棒性的統(tǒng)一。

4 結(jié) 論

為了能夠在提升目標(biāo)跟蹤精度的同時更好地應(yīng)對測量中可能出現(xiàn)的異常情況,本文提出了一種魯棒平方根連續(xù)-離散自適應(yīng)最大相關(guān)熵CKF算法。通過仿真實(shí)驗可得到以下結(jié)論。

(1)RSRCD-AMCCKF算法基于連續(xù)-離散算法跟蹤目標(biāo)實(shí)際運(yùn)動,較傳統(tǒng)目標(biāo)跟蹤算法中采用的離散-離散模型更準(zhǔn)確,跟蹤精度更高。

(2)RSRCD-AMCCKF算法將自適應(yīng)因子與改進(jìn)的最大相關(guān)熵準(zhǔn)則相結(jié)合,增強(qiáng)了算法的環(huán)境適應(yīng)性。該算法能夠有效限制測量中的可能發(fā)生的不利情況,并且比傳統(tǒng)魯棒算法精度更高,在目標(biāo)跟蹤系統(tǒng)中優(yōu)勢更為明顯。