论动态自适应滤波

第30卷　第4期测　绘　学　报

V ol.30,N o.4　2001年11月

ACT A G EOD AET ICA et CAR T O GR AP HICA SI NI CA

N ov.,2001

文章编号:1001-1595(2001)04-0293-06中图分类号:P207 文献标识码:A

论动态自适应滤波

杨元喜1

,何海波2

,徐天河

(1.西安测绘研究所,陕西西安710054; 2.解放军信息工程大学测绘学院,河南郑州450052)

Adaptive Robust Filtering for Kinematic GPS Positioning

YANG Yuan -x i 1,HE Hai -bo 2,XU T ian -he

(1.X i 'an Research I nstitute of Surv ey ing and M app ing ,X i 'an 710054,China ; 2.I nstitute of Sur vey ing and M app ing ,I nf or mation Engineer ing U niv er sity of PL A ,Zheng z hou 450052,China )

Abstract :A fter a br ief rev iew o f Sag e adaptiv e filter ing ,the relatio ns of analy tical ex pressio ns and cova ria nce

matr ices between the basic r ando m vecto r ,such as r esidual vecto r ,innov atio n vecto r and cor rectio n vecto r of predicted st ate ,ar e der iv ed and discussed .T he sho rt co ming s o f co var iance mat rices by windo wing the r esidual vecto rs,inno vat ion v ect or s a nd cor rection vecto rs o f kinem atic state,are analyzed.T he r ecent ly dev elo ped ro -bust filter ing ,Sag e adaptive filter ing and the adaptively ro bust filter ing a re compar ed.It is sho w n,by der iv a-tio ns and calculat ions ,t hat the new adaptiv e filt ering is not o nly simple in calculation but also r obust in con-tro lling the measurement o ut lier s and kinemat ic st ate disturbing .

Key words :Sage filt er ing ;adaptiv e est imation ;r obust estimation ;kinemat ic system

摘　要:动态导航与定位的质量取决于对动态载体扰动和观测异常扰动的认知和控制。本文首先介绍了目前广泛使用的Sage 自适应滤波,讨论了自适应滤波的残差向量、新息向量及状态参数预报值残差向量的解析关系,以及它们之间的协方差矩阵之间的关系;分析了基于新息向量、残差向量和状态参数预报值残差向量的自适应协方差估计存在的问题。对新近发展起来的抗差滤波、Sag e 自适应滤波及抗差自适应滤波进行了综合比较与分析,结果表明抗差自适应滤波解算理论与方法除自适应地估计载体状态预报向量的协方差矩阵外,还能自适应地估计任意历元观测量的权。计算结果证实,抗差自适应滤波不仅计算简单,而且能有效地控制观测异常和载体状态扰动异常对动态系统参数估值的影响。收稿日期:2000-12-18;修回日期:2001-05-21

基金项目:国家杰出青年科学基金资助项目(49825107)

作者简介:杨元喜(1956-),男,江苏泰县人,教授,主要研究方向为大地测量数据处理理论与方法。

关键词:Sag e 滤波;自适应估计;抗差估计;动态系统

1　引　言

运动定位与导航一般应用Kalm an 滤波算法。可靠的Kalm an 滤波算法要求有可靠的函数模型、随机模型以及合理的估计方法。然而,运动

物体一般难以确保规则运动,因而精确的函数模型的构造十分困难;随机模型先验信息的获取一般都是基于验前统计信息,而任何统计信息都难免不失真,尤其是难以精确表征当前物理现实和观测现实。于是围绕如何利用当前观测信息和状

态估值更新先验信息和补偿运动方程误差,成为广大工程应用领域尤其是动态定位领域研究的热点问题之一。在Kalm an滤波中,这种由估计过程自适应地调整、更新先验信息的算法称为自适应Kalman滤波。

在各类工程应用中,针对不同的工程应用背景已构造出多种自适应Kalm an滤波算法。Jazw inski提出了模型方差自适应补偿法[1],即在Kalman滤波过程中,利用观测信息自适应地生成模型误差的协方差矩阵,使预测残差与其相应的统计量保持良好的一致性。M ehra提出了利用新息序列的自适应估计(innovatio n-based adap-tiv e estimation,简称IAE)开窗逼近法[2],即要求观测方程协方差阵及状态误差协方差阵随时自适应于观测信息[3,4]。这种自适应滤波一般称为Sage滤波。它采用前m步残差序列估计观测向量的协方差阵和状态误差协方差阵。近年许多学者提出了利用抗差估计自适应地求解观测权阵和状态权阵[5～9]。为了适应于动态载体的导航,中国学者又提出了利用抗差估计和状态协方差矩阵膨胀模型相结合的自适应滤波[10]。此外还有基于多模型的自适应滤波方法(m ultiple-m odel-based adaptive estimation,简称M M AE)[3]。若将上述各自适应滤波进行分类,则M M AE滤波为函数模型自适应滤波,而其余的算法均属于随机模型自适应滤波。关于状态噪声的各类函数模型补偿法已另文作过分析和比较[11]。本文试图将目前动态GPS定位与导航中经常使用的随机模型Sage 自适应滤波与抗差自适应滤波进行理论和实用方面的比较,期望能对实际应用提供实用参考信息。2　观测值残差与新息向量之间的解析与统计关系

自适应Kalm an滤波的构造大多数都基于观测残差序列或新息(innov ation)序列。为分析各类自适应滤波的有效性,下面将详细分析残差与新息序列的解析关系及统计关系。

设有动力模型方程

X k= k,k-1X k-1+W k(1)式中,X k为t k时刻的u维状态向量; k,k-1为u×u维状态转换矩阵;W k为状态噪声向量。

t k时刻的观测模型为

L k=A k X k+ k(2)式中,L k为n k维观测向量;A k为n k×u维设计矩阵; k为观测噪声向量。设W k和 k的协方差矩阵分别为∑W k和∑k,W k,W j, k, j互不相关。又设状态预测向量为

X k= k,k-1X

k-1+W k(3)解向量为

k=X k+K k(L k-A k X k)(4)其中

K k=∑X k A T k(A k∑X k A T k+∑k)-1(5)

∑X

= k,k-1∑X k-1 T k,k-1+∑W k(6)残差向量和新息向量(也称预测残差向量)分别为

V k=A k X

k-L k(7)

V k=A k X k-L k(8) V k与V k的信息量不同。新息向量V k由t k时刻的预测状态X k推算,而残差向量V k是由t k时刻的状态估值向量求得,它含有t k时刻观测向量L k的信息;此外,由于V k相对于已经L k修正过的状态, V k相对于未经L k修正过的状态,所以,新息向量V k更能反映动态系统的扰动。

此外,由式(4)、式(7)及式(8)可得V k与V k 的解析关系

V k=(I-A k K k)V k(9) V k与V k的协方差分别为

∑V

=∑k+A k∑X k A T k(10)∑V

=∑k-A k∑X k A T k(11)并可推得

∑V

=(I-A k K k)∑V

(I-K T k A T k)(12)

注意到A k∑X

A T k与A k∑X

A T k都是半正定矩阵,显然有

tr(∑V k)

3　观测噪声协方差矩阵的自适应开窗估计法存在的问题

观测噪声协方差阵的自适应估计一般采用开窗法,即由m个历元新息向量或残差向量估计当前观测残差的协方差矩阵。前者简称为IAE滤波。如此,后者可称为RAE(residual-based adaptive estim ation)滤波。这2种开窗自适应估计法已在多篇论文中作过介绍[3,4,12,13]。

294测　绘　学　报第30卷

3.1　IAE 开窗法

设观测误差近于正态分布,若取计算窗口为

m ,则∑V k

的估值∑^

V k

可取为

∑

=1m ∑m

j =0

V k -j V T k -j

(14)

将式(14)代入式(10),则可求得t k 时刻观测信息的协方差矩阵∑^

∑

-A k ∑X k

T k

(15)

3.2　RAE 开窗法

类似于式(11),观测残差向量V k 的协方差矩阵∑V k 的估值可表示成

∑^

V k =1

m ∑m

j =0

V k -j V T k -j (16)再由式(11)可求得t k 时刻观测向量的协方差矩阵∑k 的估值,即

∑

+A k ∑X k A T

(17)

注意到由式(17)自适应地估计∑^

k ,需要用到t k 历元的∑X k 和残差向量V k ,而求解∑X k 和V k 应先有∑^

k 。为此,求解t k 历元的观测向量的协方差阵可用t k -1之前的m 个历元的信息求解,即

∑

k -1

m ∑m +1

j =1

V k -j V T k -j (18)

式(17)变为

∑

k -1

+A k -1∑X k -1

A T k -1

(19)

求得∑^

k 即可求得t k 时刻观测向量的权矩阵P k 。比较IAE 估式(14)、式(15)和RAE 估式(18)、式(19)容易看出:

1.IAE 法估计的∑^

k 中含有状态预报值X k

的误差,若X k 的误差较大,则V k 的误差必然大,由V k 计算的∑^V k 以及∑^

k 的可靠性就差。

2.RAE 法估计的∑^

k 实为t k -1时刻的

∑

^k -1

,因为欲用残差向量V k 计算∑^V k 及∑^

k 必

须要先知道X k 和V k ,而欲求X k 、V k 又必须先解

∑^

,所以只能由t

k -1

时刻以前的观测残差预测

当前历元的观测噪声协方差矩阵。这种预测的可靠性取决于当前历元观测精度与历史历元观测精度的一致性。否则这种预测协方差阵的代表性和可靠性就很难保证。

3.IAE 估计∑^

k 还有可能出现负定现象,即

∑

可能小于A k ∑^

X k

A T k 。

4.一般情况下,由式(18)求得的∑^

V k

远小于由式(14)求得的∑^

V k

。

5.由IAE 和RAE 法估算∑^

k 都要利用m 个历元的残差和新息向量,一方面增加了历史信息的存贮,而且窗口宽度m 也是一个难以确定的变数。

6.无论由式(15)还是式(19)计算∑^

k ,都是对历史精度信息的平均,这种计算结果很难表征当前历元的观测信息的离散程度。于是,这种自适应估计,很难真正实现“自适应”。

7.由IAE 还是RAE 法估算∑^

k ,要求各历元残差向量同类、同分布、且同维,否则根本不可

能由式(18)或式(19)计算∑^

k 。这一问题正是IAE 和RAE 法估算∑^

k 存在的最本质的缺陷。这一问题导致IAE 和RAE 开窗法几乎不可能应用于动态导航定位,因为载体在运动过程中,不可能保证每一历元载体上观测设备获取相同分布的信息。

4　动态模型噪声矩阵∑W 的开窗估

计存在的问题

令状态预报向量的改正数向量为 X k

X k =X k -X k

(20)容易推得

∑

- k ,k -1∑X

k -1

T k ,k -1

(21)

注意到E ( X k )=0,则 X k 的协方差估值可取为

∑^ X =1

m ∑m

j =0

X k -j X T k -j (22)将式(21)改成估值形式,即可求得∑W k 的估值∑^

W k 。

∑

- k ,k -1∑X k -1 T

k ,k -1(23)

295

第4期杨元喜等:论动态自适应滤波

然而,由式(28)估计∑^

W k

存在如下问题:

1.表达式中含有t k 历元的状态参数估值的协方差矩阵∑X k ,而∑X k 的求解往往需要先求得∑

。

∑

k 中含有∑^ X k ,而∑^

X k

是由m 个历元的 X j =X

j -X j 求得,其中包含t k 历元的X k

-X k ,这也要求先解∑^

W k

。

3.即使可以用t k -1之前m 个历元的 X j 估算

∑^ X k -1

和∑^W k -1

,并由此作为∑^

W k

的近似,但由

于t k -1历元状态有时不能可靠地反映t k 历元的状态扰动,且由m 个历元 X j X T

j 的平均值来估计∑^

X k -1

也不能反映t k -1历元的状态噪声水平,尤

其当运动状态产生大的扰动时,这种∑W k 的“自适应估计”,很难保证∑W k “适应”于实际运动载体的动态噪声水平。

为了避免上述前两个问题,可以直接估计

∑

。考虑式(4)和式(25), X k 可表示成

X k =-K k V k

(24)则有

∑

=K k ∑V k

K T k

(25)

若由式(14)求得∑^

V k

的估计,则可求得∑ X k

的

估计

∑

=K k ∑^V k K T

(26)

在稳态情况下,可以直接由∑^ X k 近似代替∑^

W k ,

即有

∑

=K k ∑^V k K T

(27)

最后应该指出,上述自适应滤波过程还存在着另一对矛盾,即若利用自适应估计法增大状态噪声协方差矩阵,则相当于减小观测噪声协方差,反之亦相反。但若同时增大或减少∑W k 和∑k ,则会出现矛盾,有时甚至会造成死循环,导致滤波结果发散。

5　抗差自适应滤波

利用抗差估计原理也可实现动态系统的自适

应滤波。基于观测误差方程(7)和状态误差方程及抗差估计极值条件[9],可求得抗差滤波解为

X k =X k +K k (L k -A k X k )(28)

K k =∑

k A T k (A k ∑X k A T

k +∑k

)

-1

(29)

式中,

∑

为L k 的等价协方差矩阵,它是观测协

方差阵的自适应估计,

∑

为状态预报值X k 的

协方差阵的自适应估计。等价协方差阵∑k 和

∑

可以由等价权阵P k 和P X k 求逆获得,而P k ,

P X k 可参考文献[8,9,14,15]。

由抗差估计方法自适应确定观测噪声协方差和状态噪声协方差,可以控制观测异常和动态模型噪声异常对状态参数估值的影响。但是由抗差估计方法进行自适应滤波也存在一定问题,因为

抗差滤波通常需要迭代,而且在求观测噪声的等价协方差阵时,需要可靠的状态估值,若状态估值受异常影响,则很难求得可靠的、表征观测噪声水平的∑k 。

为此,我们提出了在观测历元单步抗差解的基础上,整体控制状态噪声协方差阵的方法[10]

即若状态预报值X k 与历元抗差解X ~

k 相差较大

时,可将∑X k

变换成 ∑X k

(其中 <1),即

∑

/ (30)

式中, 可取为

[5,12]

　1 X ~k ≤c c

~k X

~k >c (31)

这里c 可取为2.5, X

~k 为 X ~k =‖X ~

k -X k ‖/

tr(∑X k )

(32)

6　计算与比较

本文算例所用数据与文献[12]中相同,即拟对一组机载动态GPS 观测数据进行处理。有两台T rim ble 4000SSE 型接收机,一台固定于参考站,另一台安置于飞机上。经过约10min 的初始化后,飞机起飞,整个飞行时间约90min 。

为了分析和比较多种Kalman 滤波结果中飞机非平稳扰动对状态参数的影响,我们用高精度的载波相位双差解作为参考值,而用双差单P 码进行各种滤波解算,其结果与参考值作差进行比

296

测　绘　学　报第30卷

较。滤波模型采用常速度模型。位置、速度和P 码的初始方差分别取为0.2m 2,9×10-6m 2?s -2和1m 2。而速度的谱密度取为0.2m 2s -2。状态模型方差-协方差阵直接采用文献[4]中的模型。由于不可能利用IAE 和RAE 法求解观测向量的协方差矩阵,故本算例只考虑对载体运动状态噪声协方差进行自适应估计比较。如图1,图2,图3所示,共采用如下3种方案进行计算比较: 经典Kalman 滤波; Sage 自适应滤波(开窗法);

抗差自适应滤波。

图1　经典Kalman 滤波F ig.1　Classical Kalman filter

ing

图2　Sag e 自适应滤波(m =10)Fig .2　Sag e ada ptive estim atio n (m =

10)

图3　抗差自适应滤波

Fig.3　T he new adaptive K alman filtering

而方案2中又分别取窗口宽度m =5,10,15三种情况进行比较。为节省篇幅,文中仅给出经多次试算找到的最优结果,即m =10的计算结果。

由计算结果,我们可以看出:

1.飞机飞行状态有两处显著非平稳状态,一处在大约1000历元,即飞机起飞时;另一处在

3000历元到4000历元之间,即飞机拐弯时。这两处的“抖动”对经典Kalman 滤波有明显的影响,见图1。

2.Sage 自适应滤波在一定程度上能抵制动态模型误差的影响,其结果优于经典的Kalm an 滤波,但结果并不十分明显。因为Sag e 自适应滤波并不能真正实现“自适应”,而且窗口的宽度m 也是一个难以确定的变数。若m 取得过小,所用的历史信息不足,不能求得表征飞机运动状态扰动的动态模型的方差-协方差矩阵。若m 过大,则由于使用了大量历史信息的平滑,很难反映瞬时历元运动状态出现的扰动。计算结果表明,m =10优于m =5和m =15。

3.在稳态情况下,Sage 自适应滤波有较好的滤波效果,其结果优于一般的Kalman 滤波。但实际的运动中,很难保证载体状态完全平稳。当出现剧烈扰动时,Sag e 自适应滤波结果会很差,甚至出现发散。此外,当观测有异常时,Sage 自适应滤波未采取任何措施抵制异常观测的影响。

4.抗差自适应滤波能实现真正意义上的“自适应”,特别是在载体“抖动”处,其结果明显优于经典Kalman 滤波和Sag e 自适应滤波,它不但可灵活地控制动态模型噪声异常的影响,而且可控制观测异常对状态参数估值的影响。因文献[12]已证明了抗差自适应滤波对观测异常的抵制能力,故本算例中未在观测样本中施加任何观测异常。

7　结　论

1.GPS 动态导航定位系统的Kalm an 滤波不具有抵制观测异常和载体动态非平稳扰动影响的能力。

2.Sage 自适应滤波利用IAE 和RAE 法对观测噪声协方差矩阵进行估计,但这两种估计都隐含一个前提条件,即载体各历元观测信息不仅同类、同分布,而且同维(即历元观测向量应具有相同的观测个数),否则无论基于新息或残差向量均无法估计观测向量的协方差矩阵。此外,IAE 和RAE 开窗法估计观测向量协方差矩阵还要求当前历元观测噪声与窗口内各历元观测平均噪声相近,否则IAE 和RAE 开窗法不可能可靠地估计当前历元观测噪声的协方差矩阵。

3.Sage 自适应滤波对状态扰动的“自适应”取决于窗口内各历元状态参数与当前载体状态的

297

第4期杨元喜等:论动态自适应滤波

适应程度,若所开的窗口内的状态参数改正数向量正好反映载体当前的状态扰动水平,则由开窗法求得的状态方程协方差矩阵能调节状态方程的噪声水平;反之,若窗口内各历元的状态参数改正向量与当前历元扰动水平不一致,则由于对历史状态信息的平滑,反而会歪曲当前历元的载体的扰动水平,从而使状态参数估值偏离。

4.Sag e自适应滤波需对历史信息进行存贮,一方面窗口大小具有明显的不确定性,而且窗口过大,会大大增加计算存贮。

5.抗差自适应滤波不仅能抵制观测异常的影响,而且对状态扰动具有较强的控制能力。抗差自适应滤波在计算方面具有很强的灵活性。因为它只基于状态预报值与当前历元的观测信息的不符值确定状态扰动水平。此外,抗差自适应滤波可在LS估计、Kalm an滤波、自适应滤波、抗差Kalman滤波和抗差自适应滤波中自适应地选择最适应的算法,因而做到了真正的“自适应”。

参考文献:

[1]　JA ZW IN SKI A H.Sto chast ic Pr ocesses and F ilter-

ing T heo ry[C].New Y or k:A cademic Pr ess,

1970.

[2]　M EHRA R K.On the Ident ification of V ariances

and A dapt ive K alman F ilter ing[J].I EEE T ransac-

tio ns o n A uto matic Contr ol,1970,15(2):175-184.

[3]　M OHA M ED A H,SCHWA RZ K P.A dapt ive K al-

man Filt ering for IN S/G PS[J].Jo urnal of

G eodesy.1999,73(4):193-203.

[4]　W AN G J,ST EWA RT M P,T SA K IRI M.Adap-

tiv e Kalman Filter ing fo r Integ rat ion of G PS w ith

G LO N A SS and IN S[R].Birming ham:IU GG/

IA G,1999.18-31.

[5]　K OCH K R,YA N G Y uan-xi.Robust K alman F il-

ter for Rank Deficient O bser vat ion M odel[J].Jour-

nal o f Geodesy,1998,72(8):436-441.

[6]　W AN G Y J,KU BIK K K.Ro bust Ka lman Filter

and Its G eodetic A pplicatio ns[J].M anuscr ipta

G eodaetica,1993,18(6):349-354.[7]　YA N G Y ua n-xi.Ro bust Bay esian Estimatio n[J].

Bulletin G eodesique,1991,65(3):145-150.

[8]　YA N G Yuan-x i.T he T heo r y and Applicatio n of

Robust Estimatio n[M].Beijing:Bayi Publishing

Ho use,1993.(in Chinese)

[9]　Y A N G Y uan-xi.Ro bust K alm an F ilter fo r Dy nam-

ical Sy st em[J].Jour na l o f the Z heng zho u Inst itute

o f Surv eying and M apping,1997,14(2):79-84.(in

Chinese)

[10]　Y A N G Yuan-x i,HE Hai-bo,X U G uo-chang.A

N ew A daptiv ely R obust F ilter ing for Kinematic

G eodetic Po sitioning[J].Jour nal of Geo desy,

2001,75(2):109-116.

[11]　Y A NG Yuan-x i.T he Basis in T heo ry and Co m-

pariso ns o f V ar io us Ro bust Filter M odels[C].

Beijing:Publishing Ho use of Sur vey ing and M ap-

ping,1999.51-57.(in Chinese)

[12]　HU Cong-w ei.T he Basic M odels for Sho rt G PS

Baseline a nd D ata Pr ocess M etho ds for G PS K ine-

mat ic P osit ioning[D].Wuhan:Instit ut e of

G eodesy and Geo phy sics o f Chinese Academy of

Sciences,1999.(in Chinese)

[13]　HU Guo-r ong,O U Ji-kun.T he Impro ved M ethod

of A daptiv e Kalman F iltering fo r G PS Hig h K ine-

matic Po sitio ning[J].A cta G eodaetica et Car to-

gr aphica Sinica,1999,28(4):290-294.(in Chi-

nese)

[14]　ZHOU Jiang-w en.Classical T heo r y of Err or s and

Ro bust Estimat ion[J].A cta Geo daetica et Car to-

gr aphica Sinica,1989,18(2):115-120.(in

Chinese)

[15]　ZHOU Jiang-w en,HU A N G Y ou-cai,YA N G Y uan-x i,O U Ji-kun.Robustified L east Squares M et hods[M].W uhan:P ress of Huazhong U niver-

sity of Science and T echno lo gy,1997.(in Chinese)

[16]　SCHW A RZ K P,CA N N ON M E,W ON G R V C.

A Co mpariso n o f GPS K inematic M o dels for De-

term inatio n of Po sitio n and V elo city along a T r a-

jecto ry[J].M anuscr ipta Geo daetica,1989,14(6):

345-353.

298测　绘　学　报第30卷