CN106487358A - 一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法 - Google Patents

Abstract

本发明属于信号处理技术领域，提供了一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法，经过初始化滤波方法、时间预测、定点迭代和更新误差协方差矩阵等步骤后，能够对非线性系统状态进行估计。该方法在最大相关熵卡尔曼滤波器基础上，利用统计学中的统计线性回归方法，将最大相关熵代价函数融入容积卡尔曼滤波器的框架中，从而应用到非线性系统状态估计上。该方法在做非线性系统状态估计的时候较传统的容积卡尔曼方法，有明显的鲁邦性和和抗噪能力的提升。可广泛用于非线性系统，尤其系统的噪声项是在非高斯情况下。

Description

一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法

技术领域

本发明涉及信号处理技术领域，具体涉及一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法。

背景技术

系统的状态估计是信号处理中的一个重要问题。卡尔曼滤波方法是解决线性系统的状态估计的主要方法之一，其充分利用了系统的状态空间模型和观测数据，通过求解优化问题，使得状态估计的误差最小，从而来获得系统状态的最优估计。非线性卡尔曼滤波器包括无迹卡尔曼滤波器，容积卡尔曼滤波器，以及稀疏网格求积分滤波器，这三种滤波器都是解决非线性系统的状态估计的有效工具。

状态空间模型为：

其中f(·)、h(·)分别表示状态向量、观测向量、状态转移函数和观测函数，m表示状态向量的维数，n表示观测向量的维数，v_k和w_k为随机噪声项，i、j表示变量。

估计效果的好坏可以用RMSE指标来衡量，其定义为：

其中，x_i,k表示真实的状态，表示对状态的估计，下标k表示第k时刻，下标i表示第i次蒙特卡罗，各个时刻RMSE的数值越小表示估计效果越好。

容积卡尔曼滤波器，作为状态估计中的一种有效方法，广泛应用于跟踪、导航。然而从贝叶斯滤波框架中得出的容积卡尔曼滤波，为了应用数值积分方法，需要假设噪声项为高斯分布，这严重影响了容积卡尔曼滤波在实际应用中的抗噪能力。同时也由于系统内部的变化不定，目前考虑了鲁邦性的滤波方法都是对特定应用定制特定的方法。所以，提出一种适用于非线性系统的鲁棒性和抗噪能力都提升的容积卡尔曼滤波方法显得尤为重要。

发明内容

针对现有技术中的缺陷，本发明提供一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法，适用于非线性系统，提高了系统的鲁棒性和抗噪能力。

一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法，包括以下步骤：

S1：初始化滤波：获取从高斯分布中随机产生的初始状态估计；

S2：时间预测：根据上一个时刻的状态估计和状态估计误差协方差矩阵计算第一采样点，并根据第一采样点计算经过状态转移函数变换后得到的第一变换采样点、状态预测以及预测误差协方差矩阵；

S3：定点迭代：基于时间预测得到的状态预测以及预测误差协方差矩阵计算第二采样点，并根据第二采样点计算经过观测函数变换后得到的第二变换采样点、观测预测值和观测估计误差协方差矩阵；利用统计线性回归和定点迭代的方法计算状态估计；

S4：更新状态估计误差协方差矩阵，并在获得带噪观测值后返回步骤S2循环。

优选地，所述步骤S1具体为：从均值为协方差矩阵为P₀的高斯分布中随机产生初始状态估计令k←1,熵计算所用核函数取高斯核函数其中σ为高斯核函数的核宽，k为离散时刻。

优选地，所述步骤S2包括：

S2a：获取上一个时刻的状态估计和状态估计误差协方差矩阵P_k-1，并计算第一采样点X_i,k-1：其中S_k-1是P_k-1的平方根矩阵，e_j是第j个元素为1的单位列向量；n为k时刻下的状态向量x_k的维数，i取值1～2n；

S2b：根据下式计算经过状态转移函数x_k＝f(x_k-1)变换后的第一变换采样点和状态预测

S2c：计算状态预测对应的预测误差协方差矩阵P_k|k-1：

其中，Q_k-1为状态转移噪声的协方差矩阵。

优选地，所述步骤S3包括：

S3a：基于时间预测得到的状态预测和对应的预测误差协方差矩阵P_k|k-1重新计算第二采样点X_i,k|k-1，其中S_k|k-1是P_k|k-1的平方根矩阵；

S3b：根据下式计算经过观测函数z_k＝h(x_k)变换后的第二变换采样点和观测预测值

S3c：计算观测估计误差协方差矩阵P_xz；

S3d：计算统计线性回归系数矩阵并定义以下关系：

d_k＝W_kx_k+e_k；

其中r_k为观测噪声，z_k为k时刻的带噪观测值；I为单位矩阵；另外φ_k＝diag([φ_p,k|k-1,φ_r]),φ_p,k|k-1和φ_r分别表示P_k|k-1和观测噪声的协方差矩阵R_k的平方根矩阵；

S3e：令j←1,采用定点迭代的方法计算状态估计：

其中，

上标j表示第j次迭代，对角系数矩阵和中对角元素由误差元素的第l个元素通过核函数κ_σ计算到，其中l＝1～(n+m)，

S3f：计算终止条件判断是否小于常数ε，如果是，记录迭代次数为N_k，并终止迭代，执行步骤S4；否则继续S3e。

优选地，所述步骤S4包括：

S4a：更新状态估计误差协方差矩阵

S4b：令k＝k+1,返回步骤S2。

由上述技术方案可知，本发明提供的基于统计线性回归的最大相关熵容积卡尔曼滤波方法，在最大相关熵卡尔曼滤波器基础上，利用统计学中的统计线性回归方法，将最大相关熵代价函数融入容积卡尔曼滤波器的框架中，从而应用到非线性系统状态估计上。该方法在做非线性系统状态估计的时候较传统的容积卡尔曼方法，有明显的鲁邦性和和抗噪能力的提升。可广泛用于非线性系统，尤其系统的噪声项是在非高斯情况下。

附图说明

图1为本发明所提出的最大相关熵容积卡尔曼滤波方法的流程图。

图2为实施例1中CKF和MCCKF-SLR对转弯机动目标跟踪效果图。

图3为实施例1中CKF和MCCKF-SLR对位置估计误差的对比图。

图4为实施例2中高斯噪声项增加了新噪声项后仿真得到的时间序列的估计误差对比图。

图5为实施例2中高斯噪声项增加了噪声项后仿真得到的时间序列的估计误差对比图。

具体实施方式

下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例仅用于更加清楚地说明本发明的技术方案，因此只作为示例，而不能以此来限制本发明的保护范围。

需要注意的是，除非另有说明，本申请使用的技术术语或者科学术语应当为本发明所属领域技术人员所理解的通常意义。

由于现有卡尔曼滤波器存在抗噪能力低、鲁邦性差等缺陷，因此考虑替换卡尔曼滤波器求解优化过程的代价函数，是一个很好的解决这两个问题的办法。最近基于最大相关熵代价函数提出的最大相关熵卡尔曼滤波器，作为新的线性滤波器较传统的卡尔曼滤波器在鲁邦性和抗噪能力方面提高明显，而迄今还没有人在适用于线性系统的最大相关熵卡卡尔曼方法基础上，提出适用于非线性系统的最大相关熵卡尔曼方法。基于现有的容积卡尔曼方法(CKF)不能很好处理被非高斯噪声干扰的非线性系统状态估计问题和现有的最大相关熵卡尔曼滤波器(MCKF)只能对线性系统做状态估计的缺陷，本发明提供一种新的非线性滤波方法。

一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法，包括以下步骤：

S1：初始化滤波：从均值为协方差矩阵为P₀的高斯分布中随机产生初始状态估计令k←1,熵计算所用核函数取高斯核函数其中σ为高斯核函数的核宽，k为离散时刻。

S2：时间预测：

S2a：获取上一个时刻的状态估计和状态估计误差协方差矩阵P_k-1，并计算第一采样点X_i,k-1：其中S_k-1是P_k-1的平方根矩阵，e_j是第j个元素为1的单位列向量；n为k时刻下的状态向量x_k的维数，i取值1～2n；

S2b：根据下式计算经过状态转移函数x_k＝f(x_k-1)变换后的第一变换采样点和状态预测

S2c：计算状态预测对应的预测误差协方差矩阵P_k|k-1：

其中，Q_k-1为状态转移噪声的协方差矩阵；

S3：定点迭代：

S3a：基于时间预测得到的状态预测和对应的预测误差协方差矩阵P_k|k-1重新计算第二采样点X_i,k|k-1，其中S_k|k-1是P_k|k-1的平方根矩阵；

S3b：根据下式计算经过观测函数z_k＝h(x_k)变换后的第二变换采样点和观测预测值

S3c：计算观测估计误差协方差矩阵P_xz；

S3d：计算统计线性回归系数矩阵并定义以下关系：

d_k＝W_kx_k+e_k；

其中r_k为观测噪声，z_k为k时刻的带噪观测值；I为单位矩阵；另外φ_k＝diag([φ_p,k|k-1,φ_r]),φ_p,k|k-1和φ_r分别表示P_k|k-1和观测噪声的协方差矩阵R_k的平方根矩阵；

S3e：令j←1,采用定点迭代的方法计算状态估计：

其中，

上标j表示第j次迭代，对角系数矩阵和中对角元素由误差元素的第l个元素通过核函数κ_σ计算到，其中l＝1～(n+m)，

S3f：计算终止条件判断是否小于常数ε，如果是，记录迭代次数为N_k，并终止迭代，执行步骤S4；否则继续S3e。

S4：更新状态估计误差协方差矩阵：

S4a：更新状态估计误差协方差矩阵

S4b：令k＝k+1,返回步骤S2。

该方法在最大相关熵卡尔曼滤波器基础上，利用统计学中的统计线性回归方法，将最大相关熵代价函数融入容积卡尔曼滤波器的框架中，从而应用到非线性系统状态估计上。能够在状态估计过程中提供的很好的鲁邦性，而且在非高斯噪声环境中也表现明显的优势，性能较传统的容积卡尔曼算法有很大的提升，可广泛用于非线性系统，尤其系统的噪声项是在非高斯情况下。

该方法利用了统计学中的统计线性回归(SLR)做线性化，将最大相关熵代价函数融入了卡尔曼滤波器(CKF)的框架中，因此该算法继承CKF的优点，可以直接应用于非线性系统的状态估计问题中。

算法分为时间预测和观测更新两部分，观测更新部分选取了定点迭代的方法来多步计算状态更新，直到满足迭代收敛的条件。定点迭代过程，选取预测更新的结果作为观测更新部分的定点迭代算法的初值，当满足终止条件终止迭代；经验表明，一般3到4步迭代即可。

在实际情况中，很多突变是突发的，而且噪声的统计特性很难被精确建模，因而我们的发明相比传统的容积卡尔曼滤波器有更好实际应用价值。

实施例1：

如图1-3所示，将上述卡尔曼滤波方法用于对机动目标转弯的跟踪，离散状态空间模型如下：

其中x_k和y_k是位置，和代表x和y方向上的速度，r_k和θ_k分别表示雷达测得的距离和角度，转弯速度ω_k仿真中设置为：

仿真中的其他参数设置为T＝1，(s_x,s_y)＝(0,0)。v_k和w_k是均值为0，协方差矩阵为Q和R的高斯噪声，其中，

而

根据上述基于统计线性回归的最大相关熵容积卡尔曼滤波方法，首先经过初始化滤波算法。状态实际初值和协方差矩阵初值P₀分别设置为P₀＝diag([100,10,100,10,1×10^-4])，状态估计的初值从高斯分布中随机选取。高斯核函数中的核宽σ设置为5。然后依次经过时间预测、定点迭代和观测更新处理后，得到图2、3的仿真图，本发明提供的基于统计线性回归的最大相关熵容积卡尔曼滤波方法的仿真结果见图中曲线MCCKF-SLR，可见较传统CKF，性能优势明显。其中常数ε＝1×10^-6。

实施例2：

更换仿真模型为单变量非平稳增长模型。状态空间模型如下：

其中q_k＝1。根据上述基于统计线性回归的最大相关熵容积卡尔曼滤波方法处理后进行仿真，即可得到图4和5。状态的初始值和方差分别设置为P₀＝1。和实施例1相比，实施例2改变了模型中的噪声r_k的统计特性。在图4的仿真中，以0.5的概率在原来的模型的噪声基础上增加了新的噪声项σ_r＝1。用相同方法在图5的仿真中加的噪声项为从图4和5中可以看出，随着混合高斯噪声的非高斯程度的增加，该算法的性能较CKF有越来越大的提升。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。

Claims (5)

1.一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法，其特征在于，包括以下步骤：

S1：初始化滤波：获取从高斯分布中随机产生的初始状态估计；

S2：时间预测：根据上一个时刻的状态估计和状态估计误差协方差矩阵计算第一采样点，并根据第一采样点计算经过状态转移函数变换后得到的第一变换采样点、状态预测以及预测误差协方差矩阵；

S3：定点迭代：基于时间预测得到的状态预测以及预测误差协方差矩阵计算第二采样点，并根据第二采样点计算经过观测函数变换后得到的第二变换采样点、观测预测值和观测估计误差协方差矩阵；利用统计线性回归和定点迭代的方法计算状态估计；

S4：更新状态估计误差协方差矩阵，并在获得带噪观测值后返回步骤S2循环。

2.根据权利要求1所述的基于统计线性回归的最大相关熵容积卡尔曼滤波方法，其特征在于，所述步骤S1具体为：从均值为协方差矩阵为P₀的高斯分布中随机产生初始状态估计令k←1,熵计算所用核函数取高斯核函数其中σ为高斯核函数的核宽，k为离散时刻。

3.根据权利要求1所述的基于统计线性回归的最大相关熵容积卡尔曼滤波方法，其特征在于，所述步骤S2包括：

S2a：获取上一个时刻的状态估计和状态估计误差协方差矩阵P_k-1，并计算第一采样点X_i,k-1：其中S_k-1是P_k-1的平方根矩阵，e_j是第j个元素为1的单位列向量；n为k时刻下的状态向量x_k的维数，i取值1～2n；

S2b：根据下式计算经过状态转移函数x_k＝f(x_k-1)变换后的第一变换采样点和状态预测

$X_{i,k-1}^{*}=f\left(X_{i,k-1}\right);$

${\hat{x}}_{k|k-1}=\frac{1}{2n}\underset{i=1}{\overset{2n}{\Σ}}{X}_{i,k-1}^{*};$

S2c：计算状态预测对应的预测误差协方差矩阵P_k|k-1：

$P_{k|k-1}=\frac{1}{2n}\underset{i=1}{\overset{2n}{\Σ}}{X}_{i,k-1}^{*}{X}_{i,k-1}^{*T}-{\hat{x}}_{k|k-1}{\hat{x}}_{k|k-1}^T+Q_{k-1};$

其中，Q_k-1为状态转移噪声的协方差矩阵。

4.根据权利要求1所述的基于统计线性回归的最大相关熵容积卡尔曼滤波方法，其特征在于，所述步骤S3包括：

S3a：基于时间预测得到的状态预测和对应的预测误差协方差矩阵P_k|k-1重新计算第二采样点X_i,k|k-1，其中S_k|k-1是P_k|k-1的平方根矩阵；

S3b：根据下式计算经过观测函数z_k＝h(x_k)变换后的第二变换采样点和观测预测值

$Z_{i,k|k-1}^{*}=h\left(X_{i,k|k-1}\right);$

${\hat{z}}_{k|k-1}=\frac{1}{2n}\underset{i=1}{\overset{2n}{\Σ}}{Z}_{i,k|k-1}^{*};$

S3c：计算观测估计误差协方差矩阵P_xz；

$P_{xz}=\frac{1}{2n}\underset{i=1}{\overset{2n}{\Σ}}{X}_{i,k|k-1}{Z}_{i,k|k-1}^{*T}-{\hat{x}}_{k|k-1}{\hat{z}}_{k|k-1}^T;$

S3d：计算统计线性回归系数矩阵 并定义以下关系：

d_k＝W_kx_k+e_k；

其中r_k为观测噪声，z_k为k时刻的带噪观测值；I为单位矩阵；另外φ_k＝diag([φ_p,k|k-1,φ_r]),φ_p,k|k-1和φ_r分别表示P_k|k-1和观测噪声的协方差矩阵R_k的平方根矩阵；

S3e：令j←1,采用定点迭代的方法计算状态估计：

${\hat{x}}_k^{j+1}={\hat{x}}_{k|k-1}+K_k^j(z_k-{\hat{z}}_k);$

其中，

$\left\{\begin{array}{c}{K}_k^j={\stackrel{\‾}{P}}_{k|k-1}^j{H}_k^T{(H_k{\stackrel{\‾}{P}}_{k|k-1}^j{H}_k^T+{\stackrel{\‾}{R}}_k^j)}^{-1}\\ {\stackrel{\‾}{P}}_{k|k-1}^j={\mathrm{\φ}}_{p,k|k-1}{\left(C_{n,k}^j\right)}^{-1}{\mathrm{\φ}}_{p,k|k-1}^T\\ {\stackrel{\‾}{R}}_k^j={\mathrm{\φ}}_r{\left(C_{m,k}^j\right)}^{-1}{\mathrm{\φ}}_r^T\end{array}\right.;$

上标j表示第j次迭代，对角系数矩阵和中对角元素由误差元素

的第l个元素通过核函数κ_σ计算到，其中l＝1～(n+m)，

$\left\{\begin{array}{c}{C}_{n,k}^j=diag\left({\mathrm{\κ}}_{\mathrm{\σ}}\right(e_{1,k}^j),{\mathrm{\κ}}_{\mathrm{\σ}}(e_{2,k}^j),...,{\mathrm{\κ}}_{\mathrm{\σ}}(e_{n,k}^j\left)\right)\\ C_{m,k}^j=diag\left({\mathrm{\κ}}_{\mathrm{\σ}}\right(e_{n+1,k}^j),{\mathrm{\κ}}_{\mathrm{\σ}}(e_{n+2,k}^j),...,{\mathrm{\κ}}_{\mathrm{\σ}}(e_{n+m,k}^j\left)\right)\end{array}\right.;$

S3f：计算终止条件判断是否小于常数ε，如果是，记录迭代次数为N_k，并终止迭代，执行步骤S4；否则继续S3e。

5.根据权利要求1所述的基于统计线性回归的最大相关熵容积卡尔曼滤波方法，其特征在于，所述步骤S4包括：

S4a：更新状态估计误差协方差矩阵

S4b：令k＝k+1,返回步骤S2。