秦月梅,周倩倩,楊衍波,梁 彥,潘 泉

(1.西安郵電大學(xué)自動化學(xué)院,陜西 西安 710121;2.西北工業(yè)大學(xué) 自動化學(xué)院,陜西 西安 710129)

1 引言

離散時間動態(tài)系統(tǒng)狀態(tài)估計問題廣泛存在于自動控制[1]、導(dǎo)航制導(dǎo)[2]、目標(biāo)跟蹤[3-5]等軍、民用領(lǐng)域[6].針對僅受高斯白噪聲擾動的線性系統(tǒng),卡爾曼濾波給出了最小線性均方誤差下的最優(yōu)估計.然而在諸如非合作目標(biāo)跟蹤的典型應(yīng)用場景[7-8],高斯白噪聲的假設(shè)很難成立.例如,由于持續(xù)干擾或快速采樣等,導(dǎo)致量測值受有色噪聲[9-10]的影響.又如,博弈對抗環(huán)境下,量測的先驗標(biāo)稱模型與實際模型往往失配,使獲得的量測數(shù)據(jù)表征為先驗標(biāo)稱模型上施加了缺少先驗信息的未知擾動.

一般地,針對有色噪聲處理方法有自適應(yīng)濾波法[11]、狀態(tài)增廣法[12-13]和量測差分法[14].常見的自適應(yīng)濾波算法[11],利用歷史殘差序列來估計有色噪聲的協(xié)方差矩陣實現(xiàn)對有色噪聲的擬合和修正,但是增大了濾波算法內(nèi)存和計算量.狀態(tài)增廣法[13]將有色噪聲擴維生成新的狀態(tài)方程,一方面增加了計算負擔(dān);另一方面,重構(gòu)的系統(tǒng)系數(shù)矩陣在濾波過程中可能帶來數(shù)值不穩(wěn)定.量測差分法通過相鄰歷元的觀測值線性組合構(gòu)造出只含有高斯白噪聲和一步延遲狀態(tài)的偽量測模型,有效避免上述兩種算法帶來的數(shù)值問題和病態(tài)計算.針對含有色量測噪聲的線性系統(tǒng),文獻[14]對量測差分后的偽量測模型進行3 種等效性研究分析,提出了一種更易實現(xiàn)、更易推廣到非線性系統(tǒng)的差分法.

實際工業(yè)過程中,受目標(biāo)偽裝、欺騙,以及復(fù)雜地球磁場干擾影響,量測值獲取過程中受有色噪聲和未知擾動共存影響.上述算法不具備抵制觀測異常和動力學(xué)擾動異常影響[15-16].實際非合作目標(biāo)跟蹤過程中,復(fù)雜外部環(huán)境和系統(tǒng)偏差不僅需要考慮有色噪聲,還需要進一步考慮廣義未知擾動(generalized unknown disturbance,GUD)的存在.這類擾動無任何先驗信息故無法利用噪聲進行有效建模.針對特性未知GUD的自適應(yīng)濾波器[17-19]將系統(tǒng)辨識與濾波有機地結(jié)合為一體,對未知參數(shù)進行辨識和修正.近年來,針對含GUD的線性系統(tǒng),文獻[4-5,19]通過構(gòu)造上限濾波器,自適應(yīng)地最小化狀態(tài)預(yù)測誤差、新息協(xié)方差和估計誤差協(xié)方差矩陣上界,有效抑制濾波誤差的峰值.

現(xiàn)有針對有色噪聲和廣義未知擾動下的濾波算法研究,不是將復(fù)雜干擾建模為協(xié)方差未知的零均值隨機噪聲就是假設(shè)GUD 有確定的系統(tǒng)模型.例如,文獻[13]將未知擾動建模為協(xié)方差未知的不確定噪聲,通過確定噪聲協(xié)方差的擾動域計算未知參數(shù),提出了一種混合不確定系統(tǒng)魯棒Kalman濾波器.實際上,廣義未知擾動具有間斷性且無任何先驗信息,實際干擾狀態(tài)模型難以建立.目前,文獻[5]所提的最小上限濾波器針對具有統(tǒng)計約束擾動的離散時間線性隨機系統(tǒng),具有較好的濾波效果,可有效解決廣義未知擾動,但無法適用于有色量測噪聲下帶廣義未知擾動的隨機動態(tài)系統(tǒng).針對該受干擾系統(tǒng),一方面有色噪聲的存在導(dǎo)致相鄰時刻間量測相關(guān),不滿足擴展上限濾波要求噪聲獨立的基礎(chǔ),不再適用于已有的上限濾波框架;另一方面,廣義未知擾動的存在導(dǎo)致理論的估計誤差協(xié)方差很難在線計算,無法利用有色量測噪聲下類卡爾曼濾波框架獲得高精度在線估計.

基于上述討論,本文基于上限濾波框架提出一種有色噪聲下最小上限遞推濾波框架.首先,在有色噪聲白化框架下,借助多拍優(yōu)化參數(shù)構(gòu)造新息協(xié)方差上限以表征重構(gòu)的量測方程中多拍廣義未知擾動的影響,以在線計算估計誤差協(xié)方差的上限代替估計誤差協(xié)方差,提出具有遞推結(jié)構(gòu)的有色噪聲下上限濾波(upper-bound filter with colored noises,CUBF);其 次,將估計誤差協(xié)方差上限最小化轉(zhuǎn)化為在線標(biāo)量參數(shù)優(yōu)化,通過求解線性矩陣不等式獲得最優(yōu)濾波參數(shù)以得到有色噪聲下最小上限濾波(minimum CUBF,CMUBF).

符號=:表示定義,(·)表示與相鄰公式相同的內(nèi)容.I和O分別為合適維數(shù)的單位陣和零矩陣.對于任意兩個(半)正定陣A和B,若A-B為正定陣或半正定陣,則可記為A>B或A≥B.Rn代表n維實向量所構(gòu)成的空間.

2 問題描述

考慮如下離散時間隨機線性系統(tǒng)動態(tài)方程:

其中:xk∈Rn和zk∈Rm分別為系統(tǒng)狀態(tài)向量和量測向量;Fk為系統(tǒng)傳遞矩陣;Hk為觀測矩陣;Ak+1為擾動系數(shù)矩陣,且各矩陣具有合適的維數(shù);wk∈Rn為高斯白噪聲,其協(xié)方差矩陣Qk≥0;vk∈Rm為服從下式(2)的有色噪聲,即

其中:Mk為噪聲系數(shù)矩陣;ek為高斯白噪聲,且其協(xié)方差矩陣Rk>0.wk和ek互不相關(guān).初始狀態(tài)x0為均值ˉx0,協(xié)方差P0的高斯隨機向量與wk,vk,ek均不相關(guān).

參考文獻[3],GUD統(tǒng)計模型如下所示:

本文擬在有色噪聲白化框架下,以在線計算估計誤差協(xié)方差的上限代替理論估計誤差協(xié)方差,利用線性矩陣不等式求解最優(yōu)濾波參數(shù),實現(xiàn)有色噪聲下最小上限濾波的遞推.擬設(shè)計的CMUBF具有如下優(yōu)點:

1) 以估計誤差協(xié)方差上限代替理論估計誤差協(xié)方差來刻畫未知擾動對濾波的影響,相比于未知擾動在線辨識更容易實現(xiàn),放寬了濾波器設(shè)計條件;

2) 估計誤差協(xié)方差上限的構(gòu)造不需要依賴過多參數(shù),相應(yīng)濾波器的遞推實現(xiàn)具有更強的可設(shè)計性和靈活性;

3) 用于求解最優(yōu)濾波參數(shù)的線性矩陣不等式依賴當(dāng)前量測,既保證濾波充分挖掘有效量測信息,又使得濾波器設(shè)計對建模誤差有一定魯棒性.

3 有色量測噪聲下遞推上限濾波

在有色噪聲白化框架下,以在線計算狀態(tài)估計誤差協(xié)方差的上限代替理論估計誤差協(xié)方差.通過構(gòu)造多拍濾波參數(shù)并在線優(yōu)化,實現(xiàn)有色量測噪聲下遞推上限濾波.

3.1 有色噪聲下遞推上限濾波框架設(shè)計

利用量測差分將有色噪聲白化以構(gòu)建新量測方程,即

其中殘差rk+1中含無先驗信息的δk+1和δk,因此無法通過在線計算其協(xié)方差來對新息協(xié)方差進行在線補償.文獻[4-5,19]僅針對量測模型中含有單一時刻擾動的遞推上限濾波問題,并未涉及量測模型中含有多拍廣義未知擾動的情況.針對此情況,本文引入雙濾波參數(shù)以構(gòu)造狀態(tài)估計誤差協(xié)方差上限,實現(xiàn)有色量測噪聲下上限濾波設(shè)計.考慮到雙濾波參數(shù)求解的復(fù)雜性,擬將雙濾波參數(shù)求解近似為單一參數(shù)凸優(yōu)化,以簡化濾波實現(xiàn)過程.

定義1基于式(5)-(8),定義如下滿足式(10)-(13)的有色噪聲下上限濾波器,即CUBF:

根據(jù)上述分析,可得CUBF濾波框架如下所示:

3.2 有色噪聲下遞推最小上限濾波器實現(xiàn)

理論上,通過引入雙濾波參數(shù)構(gòu)造誤差協(xié)方差上限,進行參數(shù)求解時應(yīng)同時優(yōu)化濾波參數(shù)?k和?k+1.然而,這個求解過程較為復(fù)雜且需要更多的先驗信息.故可將k時刻求解獲得的“最優(yōu)”濾波參數(shù)作為已知參數(shù)直接代入k+1時刻的參數(shù)優(yōu)化求解中,以簡化濾波實現(xiàn)過程.如下定理1討論了有色噪聲下遞推最小上限濾波的存在性條件.

定理1若滿足

定理1的證明過程見附錄.

定理1需同時滿足式(19)-(20)所示的初值問題及參數(shù)優(yōu)化問題.這里,初始估計誤差協(xié)方差往往可以設(shè)置較大數(shù)值作為誤差協(xié)方差上限.對參數(shù)優(yōu)化而言,通過求解線性矩陣不等式(24)求得滿足的濾波參數(shù)?k+1,進一步地,結(jié)合式(23)-(24)解出使得線性矩陣不等式(21)-(22)成立的最優(yōu)參數(shù)及增益,得到CMUBF的遞推實現(xiàn).

將式(18)代入式(24),結(jié)合式(23),濾波參數(shù)?k+1的求解可簡化為下述線性矩陣不等式(26),即

注1為簡化上述雙參數(shù)優(yōu)化下的線性矩陣不等式(26)的求解,在k+1 時刻,將k時刻已求解獲得的濾波參數(shù)作為已知值代入式(26),將上述雙參數(shù)優(yōu)化的線性矩陣不等式求解簡化為單參數(shù)優(yōu)化下的線性矩陣不等式求解,以保證其求解的穩(wěn)定性以及所設(shè)計的CMUBF的遞推計算實時性.

注2考慮到諸如非合作目標(biāo)跟蹤等對實時性要求較高的實際應(yīng)用場景,為進一步保證所設(shè)計的CMUBF遞推求解的實時性,可通過如下特征值求解近似上述線性矩陣不等式(26)的參數(shù)優(yōu)化,即

注3當(dāng)Mk-1=O2時,所提CMUBF可退化為廣義未知擾動下遞推上限濾波方法[5];當(dāng)?k+1=0時,所提CMUBF可退化為有色量測噪聲下離散時間動態(tài)系統(tǒng)濾波方法[14].當(dāng)Mk-1=O2,?k+1=0時,所提CMUBF可退化為傳統(tǒng)卡爾曼濾波方法.事實上,所提CMUBF可同時處理有色量測噪聲和廣義未知擾動,其遞推框架具有更加一般的意義.

為方便理解,表1和圖1分別給出了所提CMUBF的算法步驟和遞推實現(xiàn)框圖.

圖1 CMUBF遞推結(jié)構(gòu)框圖Fig.1 The recursive block diagram of CMUBF

表1 CMUBF算法步驟Table 1 The computational procedure of the CMUBF

4 仿真分析

考慮博弈對抗環(huán)境下伴隨有色量測噪聲的非合作目標(biāo)快速跟蹤,目標(biāo)在o-ξη平面進行勻速直線運動,其初始位置為(100000 150000)Tm、初始速度為(200 150)Tm/s.目標(biāo)動力學(xué)模型如式(1)所示,“?”為矩陣Kronecker積運算.Qk=I4.采樣步長T=1 s,采樣時刻k≤200 s.

該目標(biāo)被位于坐標(biāo)原點的高頻采樣雷達所探測,并獲得受有色噪聲和未知擾動污染的量測信息.相應(yīng)量測模型如式(1)所示,且有色量測噪聲建模為式(2).Hk=[1 0]?I2,Ak+1=I2,Rk=1002×I2,Mk=m×I2,其中“m”為可變的噪聲系數(shù)矩陣參數(shù),令m=0.7.

當(dāng)k<50 s 時,GUD 表現(xiàn)為周期為25 s 的梯形波,相應(yīng)起始值為350 m、峰值為700 m;當(dāng)50 s ≤k<150 s時,表現(xiàn)為形如600 sin(2πk/50)m的正弦波;當(dāng)150 s ≤k≤200 s時,GUD 表現(xiàn)為周期為25 s的梯形波,相應(yīng)起始值為-350 m、峰值為-700 m.

將所提CMUBF與卡爾曼濾波(簡記為KF)、有色量測噪聲下卡爾曼濾波(簡記為CbKF)[14]、廣義未知擾動下最小上限濾波(簡記為MUBF)[15]進行對比.所有對比算法初始狀態(tài)估計、初始誤差協(xié)方差以及未知擾動干擾保持一致,即=(100200 215,150250 170)T,P0=104×diag{5,0.1,5,0.1},且未知擾動干擾如下圖2所示.

圖2 廣義未知擾動示意圖Fig.2 The schematic diagram of GUD

圖3給出了1000次蒙特卡洛(Monte Carlo)仿真下各方法在ξ方向和η方向的位置和速度估計的均方根誤差(root mean square errors,RMSEs)對比.明顯地,從圖3可以看出,由于既沒有考慮有色噪聲,又沒有考慮GUD,KF算法的RMSEs最大;CbKF和MUBF要么僅考慮有色噪聲,要么僅考慮GUD,因此兩者RMSEs居中;所提CMUBF同時考慮了有色噪聲和GUD,相應(yīng)的RMSEs在絕大多數(shù)時刻明顯低于KF,CbKF和MUBF等對比算法的RMSEs.這種現(xiàn)象說明所提CMUBF在這種情況下取得最好的濾波精度,驗證了所提CMUBF的有效性.

圖3 各對比算法均方根誤差Fig.3 RMSEs of compared algorithms

同時,基于1000次Monte Carlo實現(xiàn),圖4給出了對比方法在不同m值(有色噪聲的相關(guān)性強度不同)下ξ方向和η方向的位置和速度RMSEs對比.在所考慮的有色噪聲系數(shù)矩陣參數(shù)m取值范圍內(nèi),由于分別考慮了有色噪聲或廣義未知擾動,CbKF和MUBF的RMSEs均值皆小于KF的RMSEs均值.進一步地,當(dāng)m較小時,MUBF的RMSEs均值小于Cb-KF的RMSEs均值,表明此時廣義未知擾動對濾波性能的影響占主導(dǎo)地位;當(dāng)m較大時,MUBF的RMSEs均值大于CbKF的RMSEs均值,表明此時有色噪聲對濾波性能的影響占主導(dǎo)地位.特別地,由于同時考慮了廣義未知擾動和有色噪聲,所提CMUBF的RMSEs均值在所考慮的m取值范圍內(nèi)皆小于所有對比算法的RMSEs均值,且其RMSEs均值隨m值的增大而基本保持不動.這表明所提CMUBF相比于所對比算法可以有效解決有色量測噪聲下帶廣義未知擾動的隨機動態(tài)系統(tǒng)濾波問題,且對m的取值具有一定的魯棒性.

圖4 不同m值下各對比算法均方根誤差均值Fig.4 Average RMSEs of compared algorithms under different values of m

圖5給出了所提CMUBF在不同正定矩陣D(即:D1,D2和D3)下基于1000次Monte Carlo 實現(xiàn)的目標(biāo)位置估計和速度估計的RMSEs曲線.這里,D1=diag{1,1},D2=diag{100,10},D3=diag{10,100}.不同正定矩陣下所提CMUBF的RMSEs曲線基本一致,表明該算法對于正定矩陣D的取值不敏感.

圖5 不同D值下CMUBF算法均方根誤差Fig.5 RMSEs of CMUBF with different values of D

最后,圖6給出了所提CMUBF 在不同初始估計誤差協(xié)方差矩陣P0(即:P01,P02,P03,P04和P05)下基于1000次Monte Carlo 仿真實現(xiàn)的目標(biāo)位置估計和速度估計的RMSEs 曲線.P01=103×diag{1,0.02,1,0.02},P02=103× diag{5,0.1,5,0.1},P03=104×diag{5,0.1,5,0.1},P04=105×diag{5,0.1,5,0.1},P05=105×diag{25,0.5,25,0.5}.當(dāng)系統(tǒng)運行時間較短時,不同初始估計誤差協(xié)方差P0下所提CMUBF 的RMSEs 波動較大(即P0較大時,速度估計的RMSEs較大;P0較小時,位置估計的RMSEs 較大),這是因為此時初始估計誤差協(xié)方差矩陣對當(dāng)前濾波影響較大;當(dāng)系統(tǒng)運行時間足夠長時,不同初始估計誤差協(xié)方差P0下所提CMUBF 的RMSEs 曲線基本一致,此時,初始估計誤差協(xié)方差陣對當(dāng)前濾波的影響幾乎可以忽略.該現(xiàn)象表明,在系統(tǒng)短時運行時,所提濾波方法的有效執(zhí)行依賴初始估計誤差協(xié)方差矩陣的設(shè)定,此時,為獲得高精度估計結(jié)果,應(yīng)取較為準(zhǔn)確的初始估計誤差協(xié)方差矩陣(比理論值適當(dāng)偏大以滿足式(19),即≥P0|0);當(dāng)系統(tǒng)長時間運行時,所提濾波方法的估計精度不依賴于初始估計誤差協(xié)方差矩陣.因此,從濾波精度角度而言,為獲得較高精度的估計結(jié)果,在條件允許的前提下,建議初始估計誤差協(xié)方差矩陣的取值略大于理論初始估計誤差協(xié)方差矩陣.

圖6 不同P0值下CMUBF算法均方根誤差Fig.6 RMSEs of CMUBF with different values of P0

5 結(jié)論

針對有色量測噪聲下帶廣義未知擾動的隨機動態(tài)系統(tǒng),本文提出了有色噪聲下遞推最小上限濾波.利用量測差分對有色噪聲白化處理,構(gòu)造了包含相鄰兩個時刻系統(tǒng)狀態(tài)和未知擾動的量測方程.在此基礎(chǔ)上,定義了遞推上限濾波框架,以在線遞推估計誤差協(xié)方差的上限取代理論估計誤差協(xié)方差的遞推,以彌補未知擾動的影響.進一步地,將估計誤差協(xié)方差上限的計算轉(zhuǎn)化為標(biāo)量參數(shù)的凸優(yōu)化處理,以獲得估計誤差協(xié)方差最小上限的有效近似,從而實現(xiàn)有色噪聲下最小上限濾波遞推.后續(xù)研究應(yīng)進一步考慮多拍量測噪聲相關(guān)性,并利用多參數(shù)優(yōu)化以獲取更加準(zhǔn)確的估計誤差協(xié)方差的最小上限.同時,考慮實際系統(tǒng)在諸多場合建模為連續(xù)時間隨機系統(tǒng),可類比離散時間卡爾曼濾波和連續(xù)時間卡爾曼濾波之間的聯(lián)系,在重點關(guān)注連續(xù)時間未知擾動建模的基礎(chǔ)上,進一步研究將所提有色噪聲下最小上限濾波推廣到相應(yīng)連續(xù)系統(tǒng).

附錄 定理1的證明

1)CUBF存在性證明.

因此,式(13)成立.

2)CMUBF存在性證明.

式(21)成立.