基于最大相关熵的高阶扩展卡尔曼滤波器设计方法

技术领域

本发明属于滤波器设计领域，具体涉及针对一类具有多项式形式的非线性非高斯系统，利用最大相关熵提出一种新型的高阶扩展卡尔曼滤波器设计方法。

背景技术

滤波器的应用在国内外各个领域都占据重要地位，其进步和发展在国民经济建设，尤其是国防建设，如实时估计和目标跟踪，发挥重要作用。1960年，针对线性系统，Kalman等人提出了最小均方误差准则下的卡尔曼滤波，并迅速得到广泛应用。为了求解非线性问题，以卡尔曼滤波为基础，相继出现了扩展卡尔曼滤波(EKF)、无迹卡尔曼滤波(UKF)和容积卡尔曼滤波(CKF)等。但是上述滤波均要求其建模误差为高斯白噪声。然而，在实际应用中，大多数的动态系统噪声大都是非高斯白噪声。本发明针对一类具有多项式形式的非线性非高斯系统，利用最大相关熵提出一种新型的高阶扩展卡尔曼滤波器设计方法(H-MCEKF)。

发明内容

本发明针对一类具有多项式形式的非线性非高斯动态系统，提出一种基于最大相关熵的高阶扩展卡尔曼滤波器设计方法。

本发明首先，给定两个一维随机变量X,Y∈R

本发明的有益效果：本发明将高阶扩展卡尔曼滤波和最大相关熵结合，可用于求解非线性非高斯系统的估计问题，可以解决在非线性非高斯系统情况下，随着非线性提高，滤波性能下降和发散的问题，可以将其应用到实时估计和目标跟踪领域中去。

附图说明

图1是本发明基于最大相关熵的高阶扩展卡尔曼滤波器设计方法的步骤流程图；

图2是本发明步骤4的流程图；

图3是实施例1情况1下位移的估计值；

图4是实施例1情况1下速度的估计值；

图5是实施例1情况1的估计误差；

图6是实施例1情况2位移的估计值；

图7是实施例1情况2速度的估计值；

图8是实施例1情况2的估计误差；

图9是实施例2情况1位移的估计值；

图10是实施例2情况1速度的估计值；

图11是实施例2情况1的估计误差；

图12是实施例2情况2位移的估计值；

图13是实施例2情况2速度的估计值；

图14是实施例2情况2的估计误差。

具体实施方式

以下结合附图对本发明作进一步说明。

本发明提出针对一类具有多项式形式的非线性非高斯系统，基于目标跟踪，利用最大相关熵提出一种新型的高阶扩展卡尔曼滤波器设计方法(H-MCEKF)。

如图1所示，本发明包括以下各步骤：

步骤(1)相关熵描述:

相关熵是描述两个随机变量之间的广义相似性度量。给定两个一维随机变量X,Y∈R

V(X,Y)＝E[κ(X,Y)]＝∫κ(X,Y)dF

其中，E为期望算子，κ(·,·)为平移不变的Mercer内核。在本文中，没有特别强调，这个核函数就是高斯核，定义如下：

上(2)式中，e＝x-y，σ>0表示核带宽。

对(2)式进行泰勒级数展开，可以得到

则(1)式的相关熵有如下表达式：

其中，E{(X-Y)

然而，在大多数实际情况下，联合分布F

则在有限数据驱动下随机变量对(X,Y)的相关熵表达式

当X,Y∈R

步骤(2)给定无人机运动的状态模型和测量模型，其均为具有强非线性特性的复杂动态系统，将状态模型中的高阶多项式定义为系统的隐性变量，将其伪线性化以及线性化表示；系统的状态模型和测量模型：

x(k+1)＝f(x(k))+w(k) (8)

y(k+1)＝h(x(k+1))+v(k+1) (9)

其中，x(k)∈R

步骤(2-1)对给定无人机运动系统的状态模型伪线性化表示：

为便于描述和理解，令m＝n＝2为例，进行示范性描述。两个状态变量x

假设状态转移函数f

其中，

定义1：

定义2：

根据定义1和定义2，式(13)的伪线性扩维形式如下

令：

则式(13)的矩阵形式为

同样，假设测量函数h

根据定义1和定义2、式(12)和(13)，可得式(15)的矩阵形式

步骤(2-2)将伪线性化后的状态模型转化成真正的线性形式，并将非线性测量模型等价改写为以状态和参数为变量的线性形式：

为了将步骤2-1节中建立的伪线性模型转化成真正的线性形式，需要建立l阶隐变量和u阶隐变量之间的动态关系。

其中，

再结合定义1和定义2、式(13)和式(14)，状态模型(8)有如下线性矩阵形式

令

X(k)＝[(x

则式(19)有如下线性化形式

X (k+1)＝A(k+1,k)X(k)+W(k) (20)

其中，W(k)为建模误差。

同理，可得到测量模型的线性矩阵形式

进一步，测量模型的线性化形式如下所示

Y(k+1)＝H(k+1)X(k+1)+V(k+1) (22)

其中，V(k+1)为建模误差。

步骤(3)对线性形式的状态模型和测量模型，利用递归滤波器设计思想得到高阶扩展卡尔曼滤波器：

对于线性模型(20)和(22)，给出基于KF的滤波器。给定初始值X(0)，当W(k)和V(k+1)都为零均值的高斯白噪声，且方差分别记为Q

递归滤波器设计思想:

P(k+1|k)＝A(k+1,k)P(k|k)A

K(k+1)＝(P(k+1|k)H

P(k+1|k+1)＝(I-K(k+1)H(k+1))P(k+1|k) (27)

其中Q(k)＝diag{Q

步骤(4)基于最大相关熵的高阶扩展卡尔曼滤波器设计，见图2：

步骤(4-1)基于k时刻得到系统状态变量X(k)估计值

记：系统状态变量X(k+1)的一步预测估计误差为

并可等价改写成一个关于系统状态变量X(k+1)的测量模型

将

其中，I为相应维数的单位阵，

这里，依据式(28)，可得系统状态变量X(k+1)的一步预测误差协方差矩阵为

上式中，

其中的Q

式(31)中

其中的v

步骤(4-2)测量模型组合中非高斯建模误差向量u(k+1)中各分量的统计独立化过程：

测量模型组合式(30)中的向量u(k+1)，是维度为

其中，B

将B

上式可被进一步简化为

D(k+1)＝S(k+1)X(k+1)+e(k+1) (36)

其中

由于

E{e(k)e

＝B

＝B

＝I (37)

因此，非高斯建模误差随机变量u(k+1)经过矩阵B

步骤(4-3)基于测量模型组合建立求解系统状态变量X(k+1)估计值的相关熵目标函数：

基于式(30)对系统状态变量X(k+1)的L维观测信息的N次采样实现的，建立求解状态变量估计值

其中，经过独立变换后，

因此，通过求取相关熵目标函数(38)极大值解，可得到系统状态变量X(k+1)在最大相关熵准则下的最优估计值

步骤(4-4)基于最大相关熵目标函数的系统状态变量估计值求解不动点，并对其迭代数值求解：

对于核函数

对目标函数J

则有

由于

X(k+1)＝f (X(k+1)) (42)

其中

上式中

这样，利用式(42)就可以基于相关熵目标函数，实现对系统状态变量估计值的不动点迭代数值求解。

步骤(4-5)利用不动点求解方程向Kalman滤波等价转换：

有式(36)中各方程，并利用著名的矩阵求逆引理，不动点求解方程(44)和(45)可等价变换如下

其中，

进一步推导过程中，利用到增益矩阵

这样，不动点数值求解方程(41)可等价改写为

从而完成了从不动点求解方程向Kalman滤波器求解的等价转换过程。

步骤(4-6)在线迭代Kalman滤波求解状态估计值

给定迭代初始值

其中

上式中

得系统状态变量X(k+1)的第t次迭代解决数值

当迭代截止于第r(k+1)步结束后，就得到X(k+1)的估计值和估计误差协方差矩阵分别为

注释：式(56)的C

式(57)的C

为了验证本发明所提方法的有效性，用两个实施例进行对比仿真实验：其中两个状态变量x

实施例1考虑状态方程为高阶多项式，测量方程为线性模型的无人机运动系统

其中，w(k+1)～0.8N(0,0.01)+0.2N(0,0.1)，且有v(k+1)～N(0,0.01)，初始值x(0)＝[1 1]

用MCEKF和所提出的H-MCEKF两种滤波方法对状态变量在两种情况下进行估计，比较位移和速度的估计值及位移和速度的估计误差。分别计算MCEKF和H-MCEKF的准确率进行比较。

实施例2状态方程和测量方程都是高阶多项式形式的无人机运动系统

其中，w

用MCEKF和所提出的H-MCEKF两种滤波方法对目标状态变量在两种情况下进行估计，比较位移和速度的估计值及位移和速度的估计误差。分别计算MCEKF和H-MCEKF的准确率进行比较。

实施例1分别以情况一σ＝0.5,ε＝10

表1 Estimation errors with different filters

实施例2分别以情况一σ＝0.5,ε＝10

表2 Estimation errors with different filters

由上述两个实施例可以看出，当测量噪声呈非高斯分布，给定过小的核带宽σ时，H-MCEKF的滤波性能稍差于MCEKF。但是，当给定一个合适的核带宽σ，H-MCEKF的滤波性能可接近MCEKF，甚至优于MCEKF。在实施例2中，尽管MCEKF对于速度的估计可达到较小的估计误差，但是对位移及整个状态变量的误差较大。此外，从表1和2可以看出，与核带宽σ相比，阈值ε对滤波性能影响不大。以上数据分析验证了所提方法的有效性。