趙 巖,王 寧,韋道知
(空軍工程大學 防空反導學院, 陜西 西安 710051)
狀態估計是從已有信息中提取有用信息[1],并對系統模型的某些參數進行校正的方法。狀態估值的精度和可靠性與系統模型是否準確描述直接相關[2]。實際系統的非線性模型可以有效還原系統的真實性,但過于復雜的模型,計算大,甚至難以計算[3]。因此,系統模型應對實際系統進行適當簡化,降低系統的非線性程度。與此同時,系統不可避免地引入了一定干擾。
目前,常用的非線性系統狀態估計方法有無跡卡爾曼濾波(Unscented Kalman Filtering, UKF)方法、粒子濾波(Particle Filtering, PF)方法以及相關改進方法。UKF通過Unscented變換,獲得狀態一步預測的一階、二階統計信息,能夠對較高非線性度系統進行有效估計,估值精度達到二階以上[4]。但UKF方法對初始噪聲相對敏感,不準確的初始誤差模型會嚴重影響UKF方法的濾波性能[5-7]。PF方法能夠處理強非線性系統[8],PF方法利用蒙特卡洛方法,隨機產生一組含有權值的粒子,近似狀態的后驗概率分布。但PF方法也會面臨粒子匱乏等問題,導致狀態估值陷入局部最優。這些非線性方法都存在計算量大、實時性相對較差的缺陷。對于線性系統,卡爾曼濾波(Kalman Filtering, KF)算法無疑是最實用和成熟的方法,廣泛應用于航空航天、工業控制、通信、經濟金融等領域。為了在非線性系統中應用KF算法,通過改進得到擴展卡爾曼濾波(Extended Kalman Filtering, EKF)算法[9],該方法利用Taylor級數,將非線性系統線性化,然后進行KF算法迭代,能夠獲得非線性系統的狀態估值。
除了系統模型,系統受到的干擾也是影響狀態估值精度和可靠性的一個重要因素[10]。一般的方法是將干擾(包括系統內部和外部干擾)轉換成隨機高斯噪聲(Stochastic Gaussian Noise, SGN)的形式。但實際系統的干擾往往相當復雜[11],單一的SGN并不能較好地模擬系統干擾。受技術條件限制,早期濾波難以處理非高斯噪聲(Non-Gaussian Noise, NGN),但隨著濾波技術的發展,產生了多種處理不同形式噪聲的方法。
文獻[12]研究了含有互相關噪聲的多傳感器系統,利用全局最優分布式算法對狀態進行估計,并將設計方法應用到目標跟蹤模型中進行仿真,驗證了方法的估計效果。文獻[13]利用濾波技術研究了含有自回歸噪聲的雙線性系統的參數估計問題,根據辨識后的模型采用改進最小二乘迭代方法獲得了有效估計值。文獻[14]通過系統辨識處理含有未知有色噪聲的雙線性系統,再利用KF算法對系統參數進行估計,通過數字仿真證實了這種方法的有效性。文獻[15]為了克服KF算法缺陷,采用最大相關熵準則對KF算法進行改進,并對非高斯噪聲系統進行狀態估計。文獻[16]采用一種回歸方法對未知統計特性的有界噪聲系統進行參數估計。上述文獻利用不同方法對線性和非線性系統的參數進行估計,都將干擾噪聲按照同一形式進行處理,沒有考慮系統噪聲的異質性對濾波帶來的影響。
綜上,本文針對基于KF框架算法在處理NGN系統時的局限性,提出一種能夠同時處理含有隨機高斯噪聲和未知有界噪聲兩種異類噪聲的改進算法。該方法利用最小均方誤差估計原則,在吸收集員估計處理未知有界噪聲能力的基礎上,調節卡爾曼濾波增益,使改進后的算法具有處理非高斯噪聲的能力。
KF算法是針對線性系統狀態的一種有效估計方法。通過賦予狀態量和協方差矩陣初值,利用不斷更新的觀測量,獲得下一時刻的狀態估值。但采用標準KF算法計算系統準確的狀態估值是有嚴格要求的。首先,系統是線性的,而在現實中,幾乎不存在真正的線性系統。然后,系統模型精確已知,而實際建立的系統模型總存在不同程度的誤差和干擾,很難給出準確的模型。另外,對系統噪聲也有一定的要求。但是,KF算法是一種遞推算法,特別適合計算機應用,被廣泛應用于航空航天、工業控制、通信、經濟金融等領域。
標準KF算法的計算流程如圖1所示。KF框架是根據KF算法流程,利用狀態量的一、二階矩統計信息,代入觀測信息后,通過狀態方程和觀測方程,實現時間更新和量測更新的過程。采用KF框架的濾波算法主要有KF算法、EKF算法及其改進算法。

圖1 標準卡爾曼濾波算法流程Fig.1 Flow chart of the standard Kalman filtering algorithm
EKF算法是KF算法在非線性領域的推廣。采用KF框架,通過將非線性模型進行Taylor展開,對二階以上項(含二階項)進行截斷,獲得狀態矩陣和觀測矩陣的Jacobian矩陣,然后再利用KF框架進行計算。雖然線性化過程產生了誤差,但是解決了KF算法在非線性系統中的應用問題。對于非線性度較大的系統模型,采用EKF算法進行計算會導致狀態估值發散。
基于KF框架的濾波算法對系統的噪聲有嚴格的要求。
當噪聲是已知的高斯噪聲時,基于KF框架的算法是一種遞推最小均方誤差(Mean Square Error, MSE)估計準則。當噪聲分布為非高斯、且前二階矩已知時,該估計器為最優線性估計器。但在實際情況中,系統噪聲分布具有不確定性,且統計信息難以準確獲得,因此很難滿足對噪聲的這一要求[17]。盡管基于KF框架的算法在實際應用的條件很苛刻,但適合計算機運算,且計算量相對較小。因此,將實際模型簡化后,KF算法仍然是目前應用最為廣泛的、獲得次優狀態估值的高效方法。
對于某一系統模型,采用基于KF框架的濾波算法進行解算時,系統的過程噪聲Wk和觀測噪聲Vk應同時滿足:
(1)
式中:Qk為系統噪聲序列的方差陣,Rk為觀測噪聲序列的方差陣,且Qk與Rk相互獨立;δkj為Kronecher函數。
如果系統噪聲分布不能滿足高斯分布,基于KF框架的濾波算法性能會顯著下降,甚至發散。在實際應用中,大多數噪聲是不滿足這種假設的,有的噪聲甚至難以用數學語言描述。
針對這類隨機噪聲,為了獲得更加準確、可靠的系統狀態估值,集員濾波(Set Membership Filtering, SMF)算法進入了研究人員的視野。SMF算法是一種在未知但有界(Unknown But Bounded, UBB)噪聲假設下的估計方法,其估計目的是尋求所有與模型結構、觀測數據和噪聲有界假設相容的狀態或參數組成的可行集合(Feasible Set, FS)[16]。通過該方法得到的估計結果不是一個準確值,而是一個集合,即可行集。狀態參數可行集中每個元素都是系統有效的狀態估值。一般將可行集的中心作為狀態的點估計,可行集的測度作為衡量SMF算法有效性的標準。
如前,集員濾波器是將UBB噪聲轉換到橢球集合中,用橢球集合表示這些參數的可能分布。
定義1一個橢球集合可表示為:
E(a,P,σ)={x∈Rn:(x-a)T(P)-1(x-a)≤σ2}
(2)
式中:a為橢球E(a,P,σ)的中心;P表示橢球E(a,P,σ)的形狀矩陣,且為正定矩陣;σ為橢球E(a,P,σ)半徑。
對于系統的UBB過程噪聲wk和觀測噪聲vk,用橢球集合可以表示為:
(3)

SMF算法分為時間更新和量測更新兩步,如圖2所示。時間更新即通過前一時刻狀態橢球和由狀態方程獲得的橢球(含過程噪聲橢球集合)進行Minkowski和計算,得到的狀態一步預測橢球集合的過程。量測更新即把得到的狀態一步預測橢球和觀測方程橢球(含觀測噪聲橢球集合)進行交集計算,得到當前時刻狀態橢球集合的過程。

圖2 集員濾波算法流程Fig.2 Flow chart of the set-membership filtering
狀態xk-1屬于橢球集合E(k-1),表示為:
(4)
狀態一步預測xk,k-1的集合為:
xk,k-1∈Fk-1E(k-1)+Wk-1
(5)
式中,Fk為系統的狀態矩陣,Wk-1為過程噪聲橢球集合。
定義2兩個橢球集合E(a1,P1,σ1)與E(a2,P2,σ2)的Minkowski和包含于的外定界橢球集合為:
E(a,P,σ)?E(a1,P1,σ1)+E(a2,P2,σ2)
(6)
式中:
其中:p的值需要根據選擇最小化橢球集合的方法而定。常用的方法有最小容積橢球和最小跡橢球[18],但無論采用哪種最小化橢球的方式,σ2最終可被消除。式(6)可最終表示為:
E(a,P*,1)?E(a1,P1,σ1)+E(a2,P2,σ2)
(7)
式中:
(8)
橢球集合Fk-1E(k-1)和過程噪聲橢球集合Wk-1的Minkowski和一般不為橢球集合,為使xk,k-1包含于一個標準的橢球集合,根據定義2,狀態一步預測xk,k-1的集合為:
xk,k-1∈E(k,k-1)?Fk-1E(k-1)+Wk-1

(9)
觀測橢球集合為:
(10)
式中,Hk為系統的觀測矩陣,Vk觀測噪聲橢球集合。
量測更新橢球集合為:
xk∈E(k)=E(k,k-1)∩S(k)
(11)
定義3兩個橢球集合E(a1,P1,σ1)與E(a2,P2,σ2)的交集包含于的外包界橢球集合為:
E(a,P,σ)=E(a1,P1,σ1)∩E(a2,P2,σ2)
(12)
式中:
根據式(12),式(11)可表示為:
(13)
針對EKF算法在處理非線性系統時,受到系統模型準確性和高斯噪聲模型已知的要求限制,設計一種基于KF框架的異類噪聲處理方法。首先區分系統中的異類噪聲,然后對濾波增益進行優化,最后給出改進算法。
模型是真實系統的高度抽象,建立與實際非線性系統一致的數學模型是極其困難的。即使能夠根據物理特性構建系統的模型,但系統仍然會受到外界環境的未知干擾。其中有些干擾能夠獲得統計分布,而大多數干擾難以獲得。加之非線性系統線性化的截斷誤差,使得EKF算法對非線性系統狀態進行估計時,無法獲得狀態的最優估值。因此,將系統無法準確描述的線性化截斷誤差、外界未知干擾和系統噪聲進行處理,轉變成UBB噪聲,與系統中可以描述的高斯噪聲共同構成本文研究異類噪聲的主體。因此對系統可能出現的異類噪聲進行分類,如表1所示。

表1 系統異類噪聲分類Tab.1 Classification of the system heterogeneous noise
通過調整濾波增益,使EKF算法能夠同時處理這兩類異類噪聲。
引理1若存在矩陣A、B和X,則下列矩陣函數的跡對X求導的等式成立。
(14)
定理1若EKF算法的濾波增益Kk滿足
(15)


證明:
假設系統狀態估值為:
(16)
類似地,假設觀測量為:
(17)
則狀態的MSE為:
(18)
存在
(19)

根據標準KF狀態估值方程,得到:
(20)
計算MSE,并將式(20)代入,則:
(21)
再將式(17)代入式(21),可得:

(22)
(23)
并將式(23)代入式(22),得到:

(24)

(25)
式(25)右邊第三項利用式(7)可得:
?E(0,Pk(p),1)
(26)
根據式(8)可得:
(27)
因為
(28)
所以
(29)
根據式(18)和式(25)可知:
(30)

(31)

(32)
□
綜上,基于KF框架的優化增益異類噪聲處理方法如式(33)所示。
(33)

建立粗糙的車輛運動模型作為導航系統的狀態方程。
xk=Fk-1xk-1+Wk-1
(34)

其中,sE,sN,vE,vN,aE,aN分別表示東向位置、北向位置、東向速度、北向速度、東向加速度和北向加速度,t為采樣時間。
zk=hk(xk)+vk
(35)
其中:zk=[d,θ]T;d,θ分別為里程計輸出的距離信息和方位陀螺輸出的方位角信息;hk為非線性觀測函數,其Hessian矩陣為Hk。


采用本文設計的濾波方法和EKF算法對4.1節中的車輛導航模型進行仿真,經處理得到東向、北向的位置和速度誤差結果,如圖3~8所示。
圖3和圖4為兩種濾波算法對系統位置誤差的估計曲線。由于前30個歷元兩種算法估值接近,為了更加清楚地比較位置誤差估計結果,分別在圖3和圖4左下角繪制出前30個歷元東向和北向的位置誤差估計曲線。由圖3和圖4可知,EKF算法可將位置誤差均值控制在15 m以內,且估計精度逐漸變差;本文算法具有較穩定的估計精度,能夠使位置誤差均值控制在4 m以內。造成這樣的原因是里程計和陀螺儀傳感器的誤差積累使算法計算精度有所下降,而本文算法具有較好的克服系統誤差的能力。

圖3 東向位置誤差曲線Fig.3 Eastern position error curves

圖4 北向位置誤差曲線Fig.4 Northern position error curves
圖5和圖6為兩種濾波算法對系統速度誤差的估計曲線。同樣,分別在圖5和圖6左下角繪制出前30個歷元東向和北向的速度誤差估計曲線。由圖5和圖6可知,EKF算法可將速度誤差均值控制在3 m/s以內,且估計精度逐漸變差;本文算法具有較穩定的估計精度,能夠使速度誤差均值控制在0.3 m/s以內。
圖5和圖8為50次蒙特卡洛的后輸出的位置和速度的均方誤差曲線,從圖中可以看出,本文提出的基于異類噪聲處理機制的KF算法精度優于EKF算法的計算精度。仿真得到的位置和速度誤差的均值和標準差如表2所示。本文算法的狀態估值均值和標準差均小于EKF算法的計算結果。EKF算法的計算時間為0.165 s,本文算法的計算時間為0.132 s。算法耗時是在Windows 10系統主頻2.2 GHz處理器條件下,由MATLAB R2015a軟件計算得到的。通過比較可以看出,兩種算法的計算量相當,本文算法的計算時間略短于EKF算法。因此,本文算法的估計精度高于EKF算法,且具有更好的抗誤差擾動能力,一定程度上能夠提高車輛導航定位的精度。

圖5 東向速度誤差曲線Fig.5 Eastern speed error curves

圖6 北向速度誤差曲線Fig.6 Northern speed error curves

圖7 位置均方誤差Fig.7 Position of the RMSE

圖8 速度均方誤差Fig.8 Speed of the RMSE

表2 水平位置誤差比較Tab.2 Horizontal position error compare
本文提出一種基于卡爾曼濾波框架的優化異類噪聲處理算法,在吸收集員濾波處理UBB噪聲優點的基礎上,通過調整濾波增益使系統狀態估值均方誤差達到最小,且能夠同時處理含有高斯噪聲和未知有界噪聲的系統。并將該算法和EKF算法應用到車輛導航系統模型中,在給定相同噪聲的條件下,當UBB噪聲上限相對較小時,通過數字仿真實驗發現本文方法對系統噪聲且具有更好的抗誤差擾動能力。