基于Masreliez-Martin的鲁棒分数阶容积卡尔曼滤波算法及应用

穆 静, 严东升, 蔡远利, 王长元

(1. 西安工业大学计算机科学与工程学院, 陕西 西安 710021; 2. 北京航天长征飞行器研究所,北京 100076; 3. 西安交通大学电子与信息学部, 陕西 西安 710049)

0 引 言

国内外学者研究了许多非线性系统的免微分滤波算法,如无迹卡尔曼滤波(unscented Kalman filter, UKF)采用了一组确定Sigma点的Unscented变换[1]、离差差分卡尔曼滤波算法(divided difference Kalman filter, DDKF)利用斯特林插值公式计算均值和协方差[2]、容积卡尔曼滤波(cubature Kalman filter, CKF)利用球面径向容积准则逼近状态后验分布[3]。CKF克服了UKF和DDKF在高维非线性状态估计中的缺点,同时能够具有更高的滤波精度,获得了广泛的应用[4-5]。

在实际应用中,被估计系统的函数模型和噪声统计特性常常不能够准确获得, 直接采用上述的免微分滤波算法会导致估计精度大幅度降低。为了提高估计精度,许多学者提出了鲁棒滤波算法。文献[6]提出了基于非高斯量测噪声下的Huber自适应UKF算法,并讨论了关键参数对算法的影响。文献[7]提出了基于Huber方法的高阶容积卡尔曼跟踪算法,使用统计线性回归模型近似非线性量测模型的基础上, 利用Huber M估计算法实现状态的量测更新,结合高阶容积准则构成了所提的算法,同时分析了Huber代价函数的调节因子对算法跟踪性能的影响。文献[8]研究了基于鲁棒CKF的自适应目标跟踪算法,借鉴Huber等价权函数的思想, 构造了基于平方根平滑逼近函数的修正因子以抑制观测野值的影响, 并结合CKF器求解框架推导出该算法。文献[9]讨论了基于Huber的鲁棒广义高阶CKF算法,采用 Huber方法处理观测量,进行标准的广义CKF量测更新,从而实现算法的鲁棒化,充分利用了容积变换的优势,无需通过统计线性回归模型对系统的非线性量测模型进行近似。文献[10]提出了非高斯噪声下的非线性随机系统的基于Masreliez-Martin(M-M)方法的鲁棒扩展卡尔曼滤波算法。进一步的,文献[11]提出了一种基于M-M方法的自适应和鲁棒UKF算法,可以在存在未知过程噪声和测量噪声协方差矩阵的情况下提供可靠的状态估计。

另一方面,实际应用中对有些系统的建模本身需要用分数阶来描述,而有些系统使用分数阶模型可以更好建模所研究的系统,因此分数阶系统成为了越来越多学者的研究热点。针对高斯环境下的离散分数阶随机系统,文献[12]提出了分数阶(扩展)卡尔曼滤波并应用到目标的状态估计。文献[13]以概率论为基础,从分数阶的非线性动态系统模型出发,推导了分数阶UKF(fractional UKF, FUKF),并应用于再入飞行器的状态估计上。文献[14]使用统计线性方法和容积变换,提出了分数阶CKF(fractional CKF, FCKF)算法,相比于分数阶扩展卡尔曼滤波算法,仿真实验证明了该算法的有效性和优越性。文献[15]提出了分数阶插值CKF(fractional interpolatory CKF, FICKF)算法,并讨论了FUKF和FCKF是FICKF算法的两种特殊形式。文献[16]提出了未知先验信息状态下的分数阶中心差分滤波算法,并且,从贝叶斯滤波的视角分析了分数阶非线性系统的状态估计问题[17]。而文献[18-19]分别提出了基于新息的自适应分数阶卡尔曼滤波算法和基于模糊逻辑方法的FUKF算法。在高斯分布噪声环境下,上述方法都获得了有益的结果。进一步的,一些学者研究了某些特殊噪声下分数阶系统的状态估计问题[20-21]。

针对连续分数阶系统的状态估计问题,文献[22-23]讨论了基于Tustin 生成函数的分数阶卡尔曼滤波器设计方法, 解决了相互关联的分数阶有色过程噪声和分数阶有色测量噪声下的连续时间线性分数阶系统的状态估计问题。文献[24]讨论了连续分数阶系统的FUKF算法,相比于扩展卡尔曼滤波算法,提出的算法具有较准确的估计效果。文献[25-26]讨论了使用Tustin生成函数和分数阶导数平均值方法,提出了针对具有相关和非相关过程噪声和量测噪声的连续分数阶系统下的分数阶扩展卡尔曼滤波算法和FUKF算法。进一步,文献[27-28]研究了针对非线性连续时间系统的CKF算法的初始参数问题和未知参数和分数阶问题。

本文针对具有非高斯量测噪声的分数阶非线性系统的状态估计问题进行研究,以分数阶动态随机系统为模型,结合三阶容积原则和M-M方法,改变其量测更新过程,建立基于M-M方法的鲁棒分数阶CKF算法(M-M based robust factional CKF, MMFCKF)。由于分数阶微积分的固有特性,FCKF器具有“记忆”特性,其下一时刻的估计值不仅与系统当前时刻的状态有关,还依赖于当前时刻之前的系统的运行情况。在观测噪声污染率较高的情况下能够有效地抑制污染噪声的影响, 具有较强抗差能力。最后,将MMFCKF算法应用到再入目标的状态估计中,同时分析不同程度的污染噪声对算法性能的影响。仿真实验表明,MMFCKF优于FUKF和FCKF,同时验证了新算法的有效性和鲁棒性。

1 分数阶非线性系统和容积原则

1.1 分数阶非线性随机系统

分数阶微积分是非整数阶微积分或者是任意阶微积分,是整数阶微积分的推广[29]。自从莱布尼兹等人提出后,数学家们从不同角度给出了分数阶微积分不同的定义, 每个定义的合理性和科学性在实践中得到了广泛检验。为了易于处理状态估计问题,本文使用Grünwald-Letnikov定义[12]:

(1)

(2)

使用式(1),非高斯量测噪声的分数阶非线性离散随机系统[12,15]可表示为

(3)

(4)

zk=h(xk)+ek

(5)

(6)

式中:xk∈Rnx为状态向量;f(·)和h(·)为非线性状态方程和量测方程。式(4)计算机实现时需要使用有限个过去状态(L),即xk-L+1,xk-L+2,…,xk。zk∈Rm是量测值,α1,α2,…,αn,i∈{1,2,…,n}是系统的分数阶;wk是均值为零和方差为Qk的高斯白噪声;ek是非高斯的量测噪声,其概率分布函数如下:

p(ek)=(1-)pN(ek)+q(ek)

(7)

式中:pN(ek)服从高斯分布函数N(ek,0,R1,k),称为标称量测噪声,q(ek)服从N(ek,0,R2,k)是污染噪声,参数∈(0,1)代表了污染噪声对系统的污染度。记是k时刻的量测集合。

1.2 三阶容积原则

在贝叶斯框架下,非线性滤波问题可归结为以高斯概率分布为权值的非线性函数的积分问题[3]。首先考虑如下标准正态分布的积分问题:

(8)

式中:Integral(g)为所求积分;x为滤波估计状态向量,x∈Rnx,Rnx为积分域;g(x)为非线性函数; N(x;0,I)为标准正态分布,I为单位矩阵。使用基于三阶容积原则的数值积分方法可得式(8)的近似解[3]:

(9)

式中:ξi为基本容积点;ωi为容积点对应的权值。ξi和ωi的取值分别为

(10)

实际应用中,需要求解权值为非标准正态分布的非线性函数的积分。对于一个非标准正态分布函数N(x;μ,Σ),μ为均值,Σ为协方差矩阵,那么:

(11)

(12)

利用式(9),可得到式(12)的近似解:

(13)

2 MMFCKF方法

本文提出的算法使用文献[12]中的两个假设,其合理性是显然的[15-16]。

2.1 时间更新

(14)

(15)

(16)

(17)

已知噪声wk是均值为零的白噪声,在假设1下,可得式(16)的第2项和第3项:

(18)

(19)

将式(17)～式(19)代入式(16),可得状态预测为

(20)

(21)

(22)

(23)

(24)

(25)

(26)

2.2 量测更新

(1) 计算容积点

(27)

(2) 计算通过非线性量测方程传播的容积点

γj,k=h(Xj,k)

(28)

(3) 计算量测预测、新息方差和协方差估计

(29)

(30)

(31)

(4) 使用M-M方法,计算状态估计和协方差

(32)

(33)

其中,增益Kk和vk定义如下:

(34)

(35)

(36)

(37)

(38)

式中:L为常数。

3 仿真实验

为了验证本文所提的MMFCKF的有效性,将算法应用到再入弹道目标(re-entry ballistic target, RBT)跟踪问题[30-31]中的状态估计,并与FUKF和FCKF算法对比,进一步分析量测污染噪声因子对算法的性能影响。

3.1 再入目标跟踪模型

图1描绘了再入弹道目标在一维坐标系下雷达对再入弹道目标的跟踪过程。弹体垂直下落,雷达跟踪弹体并测量了弹体和雷达之间的距离。

图1 RBT和雷达的几何关系Fig.1 Geometry of RBT and radar

根据内大气阶段的阻力,RBT的状态方程[30-31]可以建模为

(39)

式中:状态值x=[x1(t),x2(t),x3(t)]T分别表示垂直落体的高度、向下速度和弹道系数,这里取弹道系数为β=5×10-5。式(39)是连续时间模型,对应采样时间T的离散化模型f(xk-1)为

x1(k)=x1(k-1)-tx2(k-1)+w1,k-1

(40)

x2(k)=x2(k-1)-te-βx1(k-1)x2(k-1)2x3(k-1)+w2,k-1

(41)

x3(k)=x3(k-1)+w3,k-1

(42)

这里假设wk=[w1,k,w2,k,w3,k]T是均值为零和协方差为Qk的高斯噪声。使用式(40)～式(42)直接表示弹道目标的运动模型时,会损失掉有用的状态信息,使用分数阶系统描述系统能获得更为丰富的信息[32]。定义f(xk)=[x1(k),x2(k),x3(k)]T,再入目标的分数阶状态方程可表示为式(4)的形式。

根据RBT与雷达的几何关系,雷达以采样时间间隔t测量与再入目标的距离r,再入目标在k时刻相对雷达的距离rk由测量方程给出:

(43)

在实际应用中,ek为非高斯噪声,往往伴随高尾分布而出现,独立于状态噪声wk,表达式如式(7)所示。M为雷达到落体路径的最短水平距离,H为雷达距离地面的高度。本文取M=100 000 ft,H=100 000 ft[30-31]。

蒙特卡罗仿真实验中,使用均方根误差(root mean square error, RMSE)验证算法的有效性。k时刻的位置RMSE定义为

(44)

为了进一步分析算法的有效性,本文使用累加RMSE(accumulated RMSE, ARMSE),位置ARMSE定义为

(45)

式中:N为测量数据点个数。用类似的方法定义速度和弹道系数公式的ARMSE。仿真实验中蒙特卡罗次数为100次,使用式(44)和式(45)计算RMSE和ARMSE。为了更清晰地显示各种算法的结果,后文仿真结果图中y轴上使用对数坐标。

3.2 仿真实验分析

3.2.1 MMFCKF与FUKF, FCKF算法性能比较

标称测量噪声协方差选择R1,k=1002ft2,污染测量噪声协方差选取R2,k=5002ft2。图2～图4分别给出了污染噪声因子=0.5的MMFEKF,FUKF和FCKF算法的位置、速度和弹道系数的均方根误差。

图2 目标位置估计的RMSEFig.2 RMSE of target position estimation

从图2中可以看出,相比于FUKF和FCKF,MMFCKF的位置状态估计RMSE大幅度下降,滤波器收敛速度较快。图3和图4显示了MMFCKF的速度和弹道系数的RMSE从32 s时开始大幅下降,说明MMFCKF对于状态突变具有较强的跟踪能力,这是由于当存在目标状态突变时,虽然系统模型失准,M-M方法能根据测量残差的变化修正一步预测协方差阵,进而在线实时调整增益矩阵Kk,这样就能充分提取残差序列中的有效信息,使跟踪精度有了明显提高。同时,滤波收敛速度较快,说明MMFCKF对于状态突变具有较强的跟踪能力。

图3 目标速度估计的RMSEFig.3 RMSEs of target velocity estimation

图4 目标弹道系数估计的RMSEFig.4 RMSEs of target ballistic coefficient

使用ARMSE进一步对提出的算法进行性能分析,3种算法的ARMSE的柱状图如图5所示。图5显示,相比于FUKF和FCKF算法,MMFCKF算法的位置、速度和弹道系数的ARMSE大幅度减小。相比于FUKF 和FCKF 算法,MMFCKF算法的位置ARMSE分别降低了35.47%和35.22%,速度ARMSE分别降低了20.98%和20.93%,弹道系数ARMSE分别降低了6.23%和6.06%,再次表明了MMFCKF 具有很强的跟踪能力。从图2～图5也表明了FUKF 和FCKF 算法对目标的跟踪性能相似。

图5 不同算法的目标ARMSEFig.5 Target ARMSEs of different algorithms

3.2.2 噪声污染因子对算法性能的影响

图6 不同的MMFCKF的目标位置RMSEFig.6 Target position RMSEs for MMFCKF with various

图7 不同的MMFCKF的目标速度RMSEFig.7 Target velocity RMSEs for MMFCKF with various

图8 不同的MMFCKF的目标弹道系数RMSEFig.8 Target ballistic coefficient RMSEs for MMFCKF with various

从图6和图7可以看出,随着污染噪声因子增加,MMFCKF算法的位置和速度RMSE稍微在增加,但增加时相对稳定,对于弹道系数的RMSE并没有随着污染噪声因子的增加出现明显的增加,意味着MMFCKF算法可以有效抑制污染的测量噪声,M-M方法可以有效抑制污染噪声,大大的提高了估计精度,体现了MMFCKF算法在污染噪声环境下对目标跟踪的优越性。

4 结 论

本文研究了分数阶离散非线性系统状态估计问题,利用M-M方法对量测更新过程进行改进,构造了基于M-M方法的鲁棒FCKF器。该算法充分使用了量测数据的有用信息,抑制了量测噪声中的污染噪声对算法性能的影响,具有很好的稳定性。大量仿真结果表明,与FUKF算法、FCKF算法相比,MMFCKF算法能够很好地抑制污染噪声,具有更好的性能。最后,分析了不同程度的量测污染噪声对算法的性能影响,进一步验证了算法具有较高的鲁棒性。