梁 浩 張子劍 賈 睿 崔利軍

1.北京宇航系統(tǒng)工程研究所，北京 100076 2.太原衛(wèi)星發(fā)射中心，山西太原 030027

0 引言

星光導(dǎo)航系統(tǒng)(CNS)具有精度高、自主隱蔽的特點(diǎn)，能夠?yàn)檩d體提供高精度姿態(tài)信息[1]；合成孔徑雷達(dá)(Synthetic Aperture Radar，SAR)是一種高精度微波成像雷達(dá)，它具有全天時(shí)、全天候、遠(yuǎn)距離工作的優(yōu)點(diǎn)[2]；CNS、SAR和捷聯(lián)慣導(dǎo)(SINS)組成的CNS/SAR/SINS組合導(dǎo)航系統(tǒng)是一種性能優(yōu)異的無(wú)依托自主導(dǎo)航方式[3-4]。

Kalman濾波限制條件比較苛刻，要求系統(tǒng)模型精確已知，然而在實(shí)際應(yīng)用中，受各種條件的限制，很難獲得噪聲準(zhǔn)確的先驗(yàn)統(tǒng)計(jì)特性；另一方面，系統(tǒng)的噪聲特性也是不穩(wěn)定的，噪聲特性因內(nèi)部器件特性變化或外部力學(xué)沖擊等環(huán)境因素，可能產(chǎn)生不可預(yù)知的變化。

為了處理這種存在偏差的系統(tǒng)，魯棒濾波被研究并發(fā)展應(yīng)用起來(lái)，M估計(jì)濾波和H∞濾波是2種最典型的魯棒濾波，被廣泛用來(lái)處理高斯分布受到污染或干擾的系統(tǒng)。M估計(jì)濾波對(duì)應(yīng)用系統(tǒng)的噪聲特性變化較遲鈍，能夠解決系統(tǒng)噪聲和測(cè)量噪聲統(tǒng)計(jì)特性不確定的問(wèn)題[5-6]，H∞估計(jì)濾波將模型誤差看作未知但有界的噪聲，使噪聲污染情況下的估計(jì)誤差最小[7]，但當(dāng)系統(tǒng)出現(xiàn)尖銳的野值時(shí)，M估計(jì)濾波和H∞濾波均會(huì)出現(xiàn)較大誤差以至發(fā)散。

本文針對(duì)組合導(dǎo)航系統(tǒng)中出現(xiàn)多種尖銳野值的情況，引入GM估計(jì)野值檢測(cè)抑制方法，構(gòu)建了魯棒PM估計(jì)濾波，應(yīng)用于高空飛行器CNS/SAR/SINS組合導(dǎo)航中，并進(jìn)行了仿真分析。

1 基于M估計(jì)的線性濾波

1964年，Huber經(jīng)過(guò)嚴(yán)格的推導(dǎo)[5]，提出了廣義極大似然估計(jì)，即M估計(jì)魯棒濾波。針對(duì)高斯噪聲受到污染的系統(tǒng)，M估計(jì)結(jié)合了l1和l2范數(shù)構(gòu)造代價(jià)函數(shù)，其魯棒性優(yōu)于l2范數(shù)估計(jì)，其優(yōu)點(diǎn)是保證最大漸近估計(jì)方差最小以及純高斯分布時(shí)l2范數(shù)估計(jì)的效率。

線性離散系統(tǒng)如下

xk=Fk-1xk-1+wk-1

(1)

yk=Hkxk+vk

(2)

式中：xk為狀態(tài)量，yk為觀測(cè)量，F(xiàn)k-1為狀態(tài)轉(zhuǎn)移矩陣；系統(tǒng)噪聲wk-1和量測(cè)噪聲vk為高斯白噪聲，方差為Qk-1和Rk。

1)預(yù)測(cè)

(3)

(4)

2)更新

(5)

針對(duì)量測(cè)更新過(guò)程，構(gòu)造線性回歸模型以提高系統(tǒng)的冗余性：

(6)

定義如下變量：

(7)

(8)

(9)

(10)

則線性回歸模型可以轉(zhuǎn)化為

zk=Mkxk+ξk

(11)

M估計(jì)求解如下代價(jià)函數(shù)的極小值，通過(guò)求解代價(jià)函數(shù)獲得狀態(tài)估計(jì)值：

(12)

式中：ei=(zk-Mkxk)i；ρ為任意函數(shù)，通過(guò)選擇ρ函數(shù)，可以使估計(jì)器獲得某些特定的性質(zhì)，Huber建議的ρ函數(shù)形式如下

(13)

式中：μ為可調(diào)節(jié)參數(shù)。而當(dāng)殘差較大時(shí)，ρ函數(shù)具有l(wèi)1范數(shù)的性質(zhì)，相反當(dāng)殘差較小時(shí)，ρ函數(shù)具有l(wèi)2范數(shù)的性質(zhì)。通過(guò)配置μ，可以使濾波器在高斯分布下具有較高的估計(jì)效率，并且能夠抑制濾波噪聲污染對(duì)系統(tǒng)的影響。

Huber在公開文獻(xiàn)中已證明，當(dāng)選擇式(13)作為ρ函數(shù)的形式時(shí)，M估計(jì)濾波對(duì)受污染的高斯系統(tǒng)具有漸近最優(yōu)魯棒性。

2 魯棒PM估計(jì)

2.1 系統(tǒng)野值的類別及影響

在實(shí)際系統(tǒng)中除了存在受污染的高斯噪聲，還會(huì)存在各種干擾野值，一般認(rèn)為在組合導(dǎo)航系統(tǒng)中出現(xiàn)的野值主要是系統(tǒng)狀態(tài)野值和量測(cè)野值2種，然而狀態(tài)轉(zhuǎn)移矩陣和量測(cè)矩陣中也會(huì)出現(xiàn)野值，稱為結(jié)構(gòu)野值，結(jié)構(gòu)野值是由時(shí)變系統(tǒng)模型匹配錯(cuò)誤和計(jì)算浮點(diǎn)數(shù)錯(cuò)誤等引起的，3種類型野值列表如下。

在回歸模型中野值點(diǎn)遠(yuǎn)離其它正常點(diǎn)，會(huì)對(duì)M估計(jì)產(chǎn)生較大的影響，M估計(jì)也僅能夠抑制系統(tǒng)的部分量測(cè)野值，當(dāng)系統(tǒng)同時(shí)存在表1中的3種野值時(shí)，M估計(jì)濾波會(huì)出現(xiàn)偏差。

表1 組合導(dǎo)航系統(tǒng)中3種野值類型

近年來(lái)，在電力系統(tǒng)狀態(tài)估計(jì)中，出現(xiàn)了一種廣義M估計(jì)(GM估計(jì))濾波算法[8]，通過(guò)構(gòu)造回歸模型并采用迭代方法求解最優(yōu)估計(jì)，GM估計(jì)對(duì)量測(cè)野值、狀態(tài)野值和結(jié)構(gòu)野值均具有魯棒性。GM估計(jì)包括Mallows型GM估計(jì)和Schweppe型GM估計(jì)2種，Mallows型GM估計(jì)在目標(biāo)函數(shù)中引入權(quán)函數(shù)以降低所有量測(cè)點(diǎn)的權(quán)重，Schweppe型GM估計(jì)只針對(duì)野值點(diǎn)進(jìn)行衰減，后者具有更高的效率。

2.2 GM估計(jì)中的PS野值檢測(cè)方法

圖1 導(dǎo)航系統(tǒng)中被污染的樣本點(diǎn)

圖2 預(yù)白化處理后的樣本點(diǎn)分布

MD(Mahalanobis Distances)是一種經(jīng)典的基于樣本點(diǎn)均值和方差的野值檢測(cè)方法，定義樣本點(diǎn)集的中心為c=[p1,p2]T，同一組狀態(tài)信息樣本點(diǎn)形成一擬合圓，令di=[di1,di2]T表示樣本點(diǎn)xi(i=1,…,m)到擬合圓的矢量，xi到擬合圓的距離為p3，那么di又可表示為di=(xi-c)(1-p3/‖xi-c‖)。MD野值檢測(cè)方法表示如下

(14)

(15)

(16)

為此，人們提出了一種魯棒的野值檢測(cè)方法，魯棒PS(Projection Statistics)檢測(cè)法[9]，PS利用樣本點(diǎn)中值和中值絕對(duì)偏差(Median Absolute Deviation, MAD)來(lái)構(gòu)造檢測(cè)方程：

(17)

式中：uk表示PSi檢測(cè)的基準(zhǔn)方向，每個(gè)點(diǎn)都需對(duì)m個(gè)方向的投影進(jìn)行檢測(cè)；PS方法受野值點(diǎn)干擾較小，當(dāng)同時(shí)存在多處野值時(shí)也能夠準(zhǔn)定位，其計(jì)算速度也較快。

2.3 魯棒PM估計(jì)濾波構(gòu)建

將GM估計(jì)中的PS野值檢測(cè)方法應(yīng)用到M估計(jì)濾波中，構(gòu)建一種新的魯棒PM(Projection M-estimation)估計(jì)濾波。構(gòu)建方法為：

(18)

3 CNS/SAR/SINS組合導(dǎo)航模型

3.1 捷聯(lián)慣導(dǎo)誤差傳播模型

捷聯(lián)慣性導(dǎo)航系統(tǒng)是飛行器CNS/SAR/SINS組合導(dǎo)航系統(tǒng)的核心，SINS的誤差傳播方程也是組合濾波器模型的基礎(chǔ)。根據(jù)SINS誤差傳播方程構(gòu)建組合系統(tǒng)的狀態(tài)方程，利用CNS和SAR的外部導(dǎo)航量構(gòu)建量測(cè)方程，設(shè)計(jì)濾波器對(duì)慣導(dǎo)誤差進(jìn)行估計(jì)。

以東北天地理坐標(biāo)系為導(dǎo)航坐標(biāo)系，CNS/SAR/SINS組合導(dǎo)航系統(tǒng)的慣導(dǎo)誤差傳播模型為[10-11]：

(19)

(20)

(21)

(22)

3.2 觀測(cè)方程

星光導(dǎo)航系統(tǒng)中星敏感器實(shí)際測(cè)得的星光矢量可表示為

(23)

星光矢量的理論測(cè)量值為：

(24)

星敏感器對(duì)星光矢量的實(shí)際測(cè)量值rb與理論測(cè)量值pb有一定的偏差：

(25)

星光導(dǎo)航系統(tǒng)中一般將2個(gè)星敏儀正交安裝，把星敏感器實(shí)際測(cè)得的2個(gè)星光矢量rb1、rb2與理論測(cè)量值pb1，pb2相減，得到CNS/SINS的觀測(cè)量，INS的水平位置與SAR輸出的水平位置之差構(gòu)成SAR/SINS量測(cè)信息[12]，則CNS/SAR/SINS系統(tǒng)的量測(cè)方程為

(26)

其中：vk為量測(cè)噪聲。

3 仿真校驗(yàn)

將常規(guī)Kalman濾波、M估計(jì)線性濾波以及魯棒PM濾波分別應(yīng)用于CNS/SAR/SINS組合導(dǎo)航系統(tǒng)進(jìn)行數(shù)學(xué)仿真實(shí)驗(yàn)。

仿真參數(shù)設(shè)置為：飛行器初始位置為東經(jīng)116°，北緯38°，速度180m/s，方位角40°，俯仰和滾轉(zhuǎn)角為0°，飛行器有0.3m/s2的加速度和0.2(°)/s的轉(zhuǎn)彎角速度。捷聯(lián)慣導(dǎo)陀螺零漂不穩(wěn)定值0.05(°)/h，加表零偏不穩(wěn)定值100μg；星敏感器導(dǎo)航精度為20″(1σ)，SAR導(dǎo)航位置精度為20m(1σ)。仿真時(shí)間100s。

由于高空飛行器導(dǎo)航系統(tǒng)長(zhǎng)時(shí)間工作，導(dǎo)航系統(tǒng)不可避免會(huì)出現(xiàn)各種野值、噪聲特性產(chǎn)生變化。下面通過(guò)設(shè)置不同的仿真條件，對(duì)3種線性濾波算法進(jìn)行對(duì)比分析。

仿真條件1：

各導(dǎo)航設(shè)備的噪聲統(tǒng)計(jì)特性為受到污染的混合高斯分布，方差是原高斯分布的3倍，污染率為ε=0.2，混合高斯分布描述如下

(33)

圖4 水平位置估計(jì)總誤差

仿真條件2：

圖6 水平位置估計(jì)總誤差

仿真結(jié)果圖3～4表明，當(dāng)導(dǎo)航設(shè)備的噪聲特性為混合高斯分布時(shí)，常規(guī)Kalman濾波呈發(fā)散趨勢(shì)、估計(jì)誤差較大，M估計(jì)濾波和魯棒PM濾波估計(jì)精度較高，驗(yàn)證了M估計(jì)濾波和魯棒PM濾波對(duì)受污染的高斯分布噪聲具有魯棒性。

圖3 三向姿態(tài)失準(zhǔn)角估計(jì)總誤差

仿真結(jié)果圖5～6表明，當(dāng)導(dǎo)航系統(tǒng)中存在尖銳野值和結(jié)構(gòu)野值時(shí)，常規(guī)Kalman濾波和M估計(jì)濾波對(duì)這些野值無(wú)法抑制，而本文提出的魯棒PM濾波使用了PS野值檢測(cè)方法和魯棒構(gòu)架，可以對(duì)污染噪聲和多種野值進(jìn)行有效抑制。

圖5 三向姿態(tài)失準(zhǔn)角估計(jì)總誤差

4 結(jié)論

分析了傳統(tǒng)的M估計(jì)濾波在處理系統(tǒng)野值上的缺陷，借鑒了電力系統(tǒng)中的魯棒PS野值檢測(cè)方法，將魯棒PM估計(jì)濾波應(yīng)用到組合導(dǎo)航系統(tǒng)中，以CNS/SAR/SINS非線性組合導(dǎo)航為應(yīng)用背景，對(duì)常規(guī)Kalman濾波、M估計(jì)線性濾波以及魯棒PM濾波進(jìn)行了對(duì)比分析。仿真結(jié)果表明，魯棒PM濾波能夠同時(shí)對(duì)受污染的高斯噪聲和多種系統(tǒng)野值進(jìn)行抑制，具有較強(qiáng)的結(jié)構(gòu)抗差能力。