基于Gauss-Hermite逼近的非线性加权观测融合无迹Kalman滤波器

李云 孙书利 郝钢

滤波算法在定位、目标跟踪、导航和故障诊断等方面发挥着重要作用[1−3].然而,单个传感器难以满足高精度、高容错性等要求,因此,多传感器融合估计技术应运而生.在过去的几十年里,线性系统的融合估计理论已经有了一系列完整的理论基础[3].目前常用的信息融合估计方法主要包括两个基本的结构:集中式融合估计和分布式融合估计.集中式融合估计将所有传感器信息进行增广,并基于增广的观测设计融合状态估计[4−5].该算法没有信息丢失,当所有传感器没有故障时,估计精度具有全局最优性,可作为其他融合算法在精度上的衡量标准,也是现在多传感器系统经常采用的融合方式之一[6−7].然而,由于集中式融合算法计算量大,在传感器数量较多的情况下,集中式融合算法会导致整个系统实时性差.特别是当存在故障传感器时可能导致滤波器发散.分布式融合算法是把各个局部状态估计送入融合中心,根据一定的融合准则进行加权得到融合估计[3,8−9].分布式融合方式具有良好的鲁棒性,计算量小且容错性强,估计精度是局部最优、全局次优的.

加权观测融合算法根据加权最小二乘准则,将集中式融合系统增广的高维观测进行压缩处理,得到降维的观测,基于降维观测设计的滤波器可以明显地减小计算负担.对于线性系统,加权观测融合算法在最小方差意义下和集中式融合算法具有数值等价性,因而具有重要的应用价值[10].然而,绝大多数系统具有非线性特性,例如,大多数定位系统观测方程是在球面坐标系下建立的,而估计和分析状态时往往又是在笛卡尔坐标系下进行的,这使得观测方程具有某种非线性特性[6−7].

近些年,基于贝叶斯估计框架和采样逼近的非线性滤波算法得到了广泛研究,例如无迹Kalman滤波器 (Unscented Kalman filer,UKF)[11−12]、容积滤波器 (Cubature Kalman filer,CKF)[13−14]、粒子滤波器(Particle filter,PF)[15],以及其他一些非线性滤波器[16].这些非线性滤波器都可以统一处理非线性滤波问题,但各具优缺点.UKF与CKF具有相同的滤波精度,区别在于粒子权值的计算上存在差异.PF在有充足粒子条件下具有较高的滤波精度,精度普遍要高于UKF与CKF,但是较大的计算负担成为了PF的一大缺点.事实上,以上提到的滤波器都可以与本文提出的加权观测融合算法相结合,形成加权观测融合滤波算法,本文将以UKF滤波器为例,给出一种非线性加权观测融合滤波算法.

非线性滤波算法的大量涌现表明了学者们对非线性问题的关注.涉及到非线性系统的融合方法也层出不穷[17−20].近年来,有学者通过随机集、人工神经网络、模糊逻辑、粗糙集、D-S证据理论等非概率方法提出了非线性融合方法[21−23].这些方法可实现非线性系统的信息融合以及决策级融合,但这些方法普遍存在信息丢失等情况,所以这些算法不具有最优性或渐近最优性.文献[24]提出了一种在线性最小方差意义下最优非线性加权观测融合UKF滤波器.该算法要求传感器观测方程是相同的,因此具有较大的局限性.文献[25]中,基于Taylor级数和UKF,提出了加权观测融合无迹Kalman滤波器.该算法可以统一处理非线性融合估计问题,但该算法需要实时计算Taylor级数展开项系数,这将带来一定的在线计算负担,而且在展开点(状态预报)偏离过大,或者Taylor级数展开项较少的时候,滤波精度难以保证.

Gauss-Hermite逼近方法[26−28]可以通过固定点采样、Gauss函数和Hermite多项式逼近任意初等函数,且具有较好的拟合效果.为了降低该逼近方法的计算负担,本文采用了分段处理方法,即将状态区间进行分段逼近,并离线计算每段的加权系数矩阵.本文主要创新点及工作如下:首先,利用分段的Gauss-Hermite逼近方法将系统观测方程统一处理,得到近似的中介函数以及系数矩阵.进而基于此中介函数、系数矩阵以及加权最小二乘法,提出了非线性加权观测融合算法.该融合算法可对增广的高维观测进行压缩降维,为后续滤波等工作降低计算负担.最后,结合UKF滤波算法,提出了非线性加权观测融合UKF滤波算法(Weighted measurement fusion UKF,WMF-UKF).该算法可以处理非线性多传感器系统的融合估计问题.与集中式融合UKF(Centralized measurement fusion UKF,CMF-UKF)算法相比,WMF-UKF具有与之逼近的估计精度,但计算量明显降低,并且随着传感器数量的增加,该算法在计算量上的优势将更加明显.本文为非线性多传感器系统信息融合估计提供了一个有效途径.在定位、导航、目标跟踪、通信和大数据处理等领域具有潜在应用价值[29−31].

1 问题阐述

考虑一个非线性多传感器系统

其中,E为均值号,上标T为转置号,δtt=1,δtk=0().

在传感器网络中,传感器的能量是有限的,为了节省能量,假设分布在空间上的传感器之间没有通信,传感器的观测数据通过网络传输给融合中心,在融合中心对数据进行压缩和滤波处理.而在工程中经常遇到的未知参数问题[32−33]、相关性问题[34−35]、传感器分布及管理[36]等问题,本文没有涉及.

本文将从集中式融合结构入手,引出本文所提出的基于Gauss-Hermite逼近的加权观测融合方法.该融合方法将观测函数分解成Gauss函数和Hermite多项式的组合形式,利用其系数矩阵对集中式融合系统观测方程进行降维,得到一个维数较低的加权融合观测方程.对加权融合观测方程与状态方程形成的加权观测融合系统进行滤波器设计,可获得与集中式融合逼近的估计精度,并降低了集中式融合估计算法的计算量.

引理1[4−5].对系统式(1)和式(2),全局最优集中式融合系统的观测方程为:

其中

并且v(0)(k)的协方差矩阵由下式给出:

对系统式(1)和式(4),应用非线性滤波算法(例如扩展Kalman滤波器(Extended Kalman filter,EKF),UKF,CKF,PF等),可得到相应的全局最优集中式融合非线性滤波器.但由于集中式融合的观测方程式(4)是观测增广扩维形成的,使得基于该高维观测的估计算法的计算负担随着传感器数量的增加而迅速增加.因此,找到等效的或者近似的融合方法来降低计算量是十分必要的.下面本文将解决非线性系统增广观测的降维问题.

定理1.对系统式(1)和式(2),若存在一个中介函数ψ(x(k),k)∈Rψ,使得局部观测函数h(j)(x(k),k)(j=1,2,···,L)满足h(j)(x(k),k)=H(j)ψ(x(k),k),其中矩阵H(j)∈Rmj×ψ,则加权观测融合系统的观测方程可由下式给出:

其中

其中,R(0)−1=(R(0))−1,并且v(I)(k)的协方差矩阵为:

其中,M(列满秩)和(行满秩)是的满秩分解矩阵:

其中,M,H(I)可以用Hermite规范形得到[25].

证明.由于M和H(I)为H(0)的满秩分解,则有:

由于M为列满秩,因而MTR(0)−1M为非奇异矩阵.令H(I)ψ(x(k),k)为观测对象,应用加权最小二乘法,则H(I)ψ(x(k),k)的最优Gauss-Markov估计为式(9)所示.□

对加权观测融合系统式(1)和式(9),应用非线性滤波算法,可得到全局最优加权观测融合非线性滤波算法.

2 Gauss-Hermite逼近

本节将引入一种函数逼近方法,该方法借由Gauss函数和Hermit多项式的组合形式逼近任意初等函数.通过此逼近方法,可得到的近似函数,进而可将统一转化为的形式,其中,ψ(x(k),k)由Gauss函数和Hermit多项式构成,H(j)为系数矩阵.非线性多传感器系统观测函数经过转换,将满足定理1中要求.

引理 2[26].设在区间[a,b]中存在一个点集,对于任意点存在yi,满足,其中y(x)是一个确定的函数.进而y(x)的近似函数可由Gauss-Hermite折叠函数得出:

其中,γ是一个与∆xi(i=1,···,S)有关的常系数,为一系列Hermite多项式的组合:

由式(17)和式(18)有:

其中,‘!’表示阶乘,双阶乘‘m!!’表示不超过自然数m且与m有相同奇偶性的所有正整数的乘积.

注1.对于多维情况,假设是一个采样集合,对于集合中每一个点存在点满足,其中Y(·)是确定的多维函数.那么Gauss-Hermite折叠函数如下:

其中,n维函数为函数Y(·) 的近似函数. 引理2给出了一种利用Gauss函数和Hermite多项式组合的逼近方法,该方法可以利用较少的函数项获得很好的逼近效果.如果将引理1中的视为定理 1 中的中介函数ψ(x(k),k),将视为H(j),则定理1可以得以实施.

3 基于Gauss-Hermite逼近的加权观测融合算法

由文献[26]和大量仿真试验表明,在p=0,2,4等情况下,合理的选择 ∆xiµ和γµ(i=1,···,S;µ=1,···,n) 即可很好地逼近任意初等连续函数. 本文选取p=2,∆xiµ=1,γµ=γ(i=1,···,S;µ=1,···,n), 则由式 (18) 和式(19)有C2=−1/4,H2(u)=4u2−2,进而有f2(u)=1.5−u2.令

定理2.对系统式(1)和式(2),基于Gauss-Hermite逼近的近似加权观测融合方程为:

证明.利用式(22)将集中式融合系统观测方程式(6)进行近似,得到近似的集中式融合观测方程:

其中

注2.定理2通过Gauss-Hermite逼近构建了一个近似的中介函数.它使得形如式(1)和式(2)的任意非线性多传感器系统的局部观测函数具有了定理1中所阐述的关系,可使定理1得以实施.

注3.如果状态范围过大,拟合采样点数量会急剧增加,导致计算量增加,因此本文采取分段的处理方法.例如,对一维状态系统,可以将状态的范围划分成多个区间,对二维状态系统,可以将状态的范围分成若干小的区域.在每个区间或区域分别进行Gauss-Hermite逼近.逼近过程中形成的中介函数及其满秩分解矩阵可离线计算,在线调用,减少了在线计算负担.

4 WMF-UKF算法

对加权观测融合系统式(1)和式(23),应用非线性滤波算法(EKF、UKF、PF、CKF等),可得加权观测融合非线性滤波算法.本文将以UKF为例,给出一种基于Gauss-Hermite逼近和UKF滤波算法的非线性加权观测融合估计算法.

4.1 基于Gauss-Hermite逼近的WMF-UKF算法

本文UKF采样策略选用比例对称抽样,即Sigma采样点可由式(31)计算.

且有粒子权值如式(32)和式(33)所示.

其中,α>0是比例因子,λ=α2(n+κ)−n,κ是比例参数,通常设置κ=0或者κ=3−n,β=2.下面给出WMF-UKF算法.

WMF-UKF算法.对非线性系统式(1)和式(2),基于定理2的WMF-UKF算法如下:

步骤1.设置初始值

基于多传感器的观测数据z(j)(0)∼z(j)(k)(j=1,2,···,L),加权观测融合系统Sigma采样点可以计算为:

其中初值条件为:

步骤2.预测方程

预测Sigma采样点:

状态预报:

状态预测误差方差阵:

观测预报Sigma采样点:

观测预报:

观测预报误差方差阵:

协方差矩阵由下式计算:

步骤3.更新方程

滤波增益由下式计算:

滤波误差协方差矩阵为:

4.2 WMF-UKF算法的时间复杂度分析

算法1中的式(45)出现了矩阵求逆运算,因此该算法的时间复杂度由决定[37],即WMF-UKF的时间复‡杂度为O(r·3),而CMFUKF的时间复杂度为.由定理2知,所以WMF-UKF的时间复杂度小于CMF-UKF.

另外,随着传感器数量L的增加,将不断增加.而在拟合采样点数S不改变的情况下,由于,故r将保持在Sn(或者更小)不改变.因此随着传感器数量的增加,WMFUKF较CMF-UKF在计算量上的优势将更加明显.

5 仿真例子

例1.考虑一个带有4传感器的非线性系统[38]

其中

w(k)和v(j)(k)(j=1,···,4)是相互独立的白噪声,方差分别为:,.状态初值为x(0)=0.由于状态x(k)介于−1～4,因此选取拟合采样点集为:{−2,−1,···,5}(8 个等间隔点),相应的系数选取为:γ=1.选择p=2,则中介函数为:

系数矩阵H(0),M和H(I)分别为:

最后得到基于Gauss-Hermite逼近的WMFUKF估计曲线和真实曲线如图1所示.

本例采用k时刻累积均方误差(Accumulated mean square error,AMSE)[24,39]作为衡量估计准确性的指标函数如式(55)所示.

其中,xi(t)是t时刻第i次Monte Carlo实验的真实值,是t时刻第i次Monte Carlo实验的估计值.独立进行20次Monte Carlo实验,得到的AMSE曲线如图2所示,其中本例选取局部UKF估计AMSE曲线(Local filter 1～4,LF 1～4)、集中式融合UKF估计AMSE曲线(CMF-UKF)以及本文提出的加权观测融合UKF估计AMSE曲线(WMF-UKF)进行对比.由图2可以看出CMFUKF与WMF-UKF具有接近的估计精度,而高于局部UKF.在计算量方面,由于本文压缩后的观测为3维,因此WMF-UKF滤波过程中的时间复杂度为O(33).而集中式融合系统观测方程为4维,因此时间复杂度为O(43).因此,WMF-UKF计算量要低于CMF-UKF.

图1 真实状态及WMF-UKF估计曲线Fig.1 Curves of the true state and the WMF-UKF estimate

图2 局部UKF,WMF-UKF以及CMF-UKF的AMSE曲线Fig.2 AMSE curves of local UKF,WMF-UKF and CMF-UKF

例2.考虑一个带有8传感器的平面跟踪系统,在笛卡尔坐标下的状态方程和观测方程如下:

经测试,本例选取Gauss-Hermite系数γ=1.04.为了减少计算量,本例将目标移动范围划分成了16个1平方公里的子区域,如图3(a)所示.每个子区域采用以该区域为中心,外扩2点的方法避免边缘拟合效果不良.以子区域7为例,以点(0,0),(0,1),(1,1)和(1,0)所围区域为中心,外扩2点得到该子区域的拟合采样点如图3(b)所示.计算该区域的系数矩阵,如图3(c)所示.不难看出,由于8个传感器位于4个点,这里至少可以将16维的集中式融合观测方程压缩成8维的加权观测融合方程.将16个区域对应的与中介函数离线计算存储并形成数据库.根据每时刻状态预报,在数据库中选取相应的以及可减少在线计算负担.

为了对比分析WMF-UKF的精度和计算量,本文选取了8传感器集中式融合UKF(8-CMF-UKF),5传感器集中式融合UKF(5-CMF-UKF)以及3传感器集中式融合UKF(3-CMF-UKF).传感器的选择原则是尽量的分散,例如,3-CMF-UKF选择的是1,3和5传感器,5-CMF-UKF选择的是1,3,5,7和8传感器.各种融合系统的滤波跟踪轨迹曲线如图4所示.

本例采用k时刻位置(x(k),y(k))的累积均方误差(AMSE)作为指标函数,如式(58)所示.

其中,(xi(t),yi(t))是t时刻第i次Monte Carlo实验的真实值,是t时刻第i次Monte Carlo实验的估计值.独立进行20次Monte Carlo实验,得到的AMSE曲线如图5所示.

图4 真实轨迹和WMF-UKF,8-CMF-UKF和5-CMF-UKF的估计曲线Fig.4 True and estimated tracks using WMF-UKF,8-CMF-UKF and 5-CMF-UKF

在精度方面,由图5可以看到AMSE由低到高依次是8-CMF-UKF,WMF-UKF,5-CMF-UKF和3-CMF-UKF.实验说明,随着传感器数量的增加,集中式融合算法的精度不断提高,而本文提出的WMF-UKF算法的精度接近全观测集中式融合8-CMF-UKF.

在计算量方面,加权观测融合系统观测方程为8维,因此时间复杂度为O(83).3传感器集中式融合系统观测方程为6维,因此时间复杂度为O(63).5传感器集中式融合系统观测方程为10维,因此时间复杂度为O(103).8传感器集中式融合系统观测方程为16维,因此时间复杂度为O(163).因此,时间复杂度由高到低依次为:8-CMF-UKF,5-CMFUKF,WMF-UKF和3-CMF-UKF.

图5 位置融合估计的AMSE曲线Fig.5 AMSE curves of position fusion estimates

此外,为了比较分析,本例应用文献[25]中的Taylor级数逼近方法得到的WMF-UKF的AMSE曲线也绘于图5中,这里我们采用2阶Taylor级数逼近.由于Taylor级数展开阶数以及展开点等原因,使得其精度低于其他融合算法.而且与本文的不需要在线计算融合矩阵的WMF-UKF算法相比,文献[25]的WMF-UKF(2-order Taylor)算法需要根据在线预报值实时计算融合参数矩阵,因而具有更大的在线计算负担.

本例根据不同Hermite多项式(p=0,2,4)情形进行了仿真分析.经离线测试,选取Gauss-Hermite系数分别为:γ=0.83(p=0),γ=1.04(p=2),γ=1(p=4),其他参数不变.得到Monte Carlo实验的AMSE曲线如图6所示.图6中可以看到,Hermite多项式的数量与函数逼近效果并无直接关系,得到融合估值精度间也不存在渐近最优性.因此,根据被逼近函数形式,离线测试逼近函数效果,对本文所提出WMF-UKF算法的精度起到非常关键的作用.

图6 带不同Hermite多项式的WMF-UKF位置AMSE曲线Fig.6 AMSE curves of WMF-UKFs with different Hermite polynomials for position

综上,合理的选择Gauss-Hermite逼近函数以及相应的系数γ,可使本文提出的WMF-UKF在精度方面接近集中式融合算法,而减少计算量.

6 结论

本文首先基于Gauss-Hermite逼近方法和加权最小二乘法,提出了一种具有普适性的非线性加权观测融合算法.进而结合UKF算法,提出了非线性加权观测融合UKF(WMF-UKF)算法.与CMFUKF算法相比,WMF-UKF具有与之逼近的估计精度,但计算量明显降低,并且随着传感器数量的增加,该算法在计算量上的优势将更加明显.本文通过仿真实例对比已有的相关算法,说明了本算法的有效性.