非线性、非高斯系统滤波的最优估算

陈 蕊

（芜湖职业技术学院，安徽 芜湖 241000）

0 引言

在实际工程应用中有很多非线性、非高斯系统，如目标跟踪[1]、轨迹规划[2]等，这类系统的状态方程或量测方程为非线性的，或者是方程中的扰动项呈现非高斯分布。一般这类模型无法获得状态变量真实的后验分布，但可以构造近似的离散分布。随着计算机技术的不断进步，针对这类非线性、非高斯系统基于滤波理论的各种系统估算方法得到众多研究人员的青睐。

1960年，卡尔曼提出了最优化自回归数据处理算法[3]，即卡尔曼滤波（Kalman Filter，KF），并得到各个领域的关注。最初卡尔曼滤波主要针对线性动态系统问题的估算，应用于信号处理、目标跟踪、经济预测等方面[4]。之后，Busy、Sumhara等人提出了基于卡尔曼滤波的扩展卡尔曼滤波理论（Extended Kalman Filter，EKF），就是将卡尔曼滤波理论进一步应用到非线性领域[5]。近些年来，又有学者提出了无迹卡尔曼滤波（Unscented Kalman Filter，UKF）[6]并受到了广泛的关注。但是，卡尔曼滤波是假设系统的噪声满足高斯状态，实际工程应用中却很难满足，这将导致卡尔曼滤波失效。

另外，粒子滤波理论自20世纪60年代也被应用于自动控制领域。粒子滤波（Particle Filter，PF），也称之为序贯蒙特卡罗法（Sequential Monte Carlo，SMC）[7]，是一种结合了蒙特卡洛和贝叶斯的滤波理论。标准的粒子滤波易出现粒子贫化和权值退化的现象，为了解决这一问题提出了重要性重采样滤波（Sampling Importance Resampling particle filter，SIR）[8]。为进一步提高系统估计精度，学者们提出了辅助粒子滤波（Auxiliary Particle Filter，APF）[9]，这种滤波方式充分利用当前的观测值，从中产生样本引入重要性概率密度函数中，以解决粒子多样性匮乏问题。

本文旨在已知系统模型和观测量，每一步噪声未知的前提下，通过仿真实验针对非线性、非高斯状态系统运用EKF、UKF、SIR和APF等四种算法进行曲线估计，最后评估估算结果，探寻最优算法。

1 EKF和UKF

1.1 EKF

EKF是将非线性系统线性化，是一种伪非线性的卡尔曼滤波。与卡尔曼滤波相比，它利用局部线性化来估计非线性系统，再进行卡尔曼滤波，因而EKF属于次优滤波。

假设非线性系统的方程为

式中，k是时间步长，xk是系统的状态向量，zk是观测向量。

EKF先将系统状态空间向量xk和观测向量zk分别围绕第k-1时刻的滤波值(x^k-1,0)和预测值(x^k|k-1,0)，通过泰勒（Taylor）级数展开，只取其一阶量，舍弃高阶量，将非线性系统近似转换为线性系统，根据k-1时刻之前的参数对于k时刻状态进行最优估计，推出接近表达式，然后通过标准卡尔曼滤波计算状态向量的估计值。

EKF存在两个过程，分别为预测和校正。其一，预测过程是依托于系统状态转移方程的时间更新过程，是根据当前的系统状态转移方程获得观测值，并对应下一时刻的预先估计。校正过程是通过预测值进行估计的测量更新过程，是用观测值对上一步进行修正的过程。整个迭代循环示意图如图1所示。

图1 扩展卡尔曼滤波工作流程图

EKF算法是由针对线性系统的卡尔曼滤波衍生而来的，简单易行。当该非线性系统比较接近线性时，效果较好，EKF有可能接近真实值，否则将会产生很大的误差。

运用蒙特开罗仿真法针对EKF进行仿真。

图2分别是EKF在蒙特卡洛仿真时第0步、第1步、第50步、第100步时值的柱状图。可以看出，在第1、50和100时间步都出现了双峰。原因是卡尔曼滤波假设各个时刻的后验密度均是高斯分布。换句话说，一旦分布不是高斯分布，EKF可能会失效。因而，EKF在解决这种非线性和非高斯估算时可能会遇到问题。

图2 EKF的蒙特卡洛仿真

图3是EKF仿真结果（其中红色曲线为真实值，蓝色为EKF估计曲线），正如之前所预料的一样，EKF的仿真效果并不佳。

图3 EKF状态估计曲线

1.2 UKF

在非线性系统中EKF算法简单易行，但是EKF把一阶以上的高阶导忽略了，使得估算误差较大，进而人们又提出了UKF。

UKF的核心思想是摒弃了EKF对非线性系统近似线性化的处理思路，为了确保参数的均值与方差不变，通过无迹转换（Unscented Transformation,UT）模拟非线性函数的概率密度达到近似效果，利用扩展参数对非线性非高斯系统的进行估计[10]。

UKF的关键点就是无迹转换（也称UT转换），其步骤大致可以简述如下：

①假设x是一个nx维随机向量，非线性变换y=g(x)，x的均值为，方差为Px的随机变量，计算2nx+1个Sigma采样点ξi和对应的权值Wi:

上式中，κ是一个标量，表示分布的尺度函数；是矩阵方根的第i行或第i列。

②Sigma采样点ξi经非线性函数传播得到一组带有权值的新样本点

③用新样本点进行非线性变换，运用下列公式计算均值和协方差

将无迹转换应用到EKF中，就可以得到UKF。UKF利用了确定的采样来近似状态的后验概率密度，不需要对非线性模型线性化，且分布统计量的计算精度可以达到二阶[11]。因而，UKF可以有效地改善EKF当系统非线性强度增加导致的滤波发散问题。

图4是UKF的仿真结果（其中红色曲线为真实值，蓝色为UKF估计曲线），可以明显看出UKF滤波效果优于EKF。

图4 UKF状态估计曲线

2 粒子滤波

粒子滤波（Particle Filter,PF）是一种基于Monte Carlo仿真的最优回归Bayes估计的滤波算法。粒子指的是先通过状态的先验分布采样得到一组随机分布的样本，然后再运用Monte Carlo方法解决贝叶斯估计中的积分运算，以此计算出系统状态的最小均方差估计。粒子数量越多滤波结果越接近真实值，如当粒子数量趋向于无穷大时，估计结果可以逼近服从任意概率分布的系统状态。

2.1 序列重要性采样算法

Doucet等人于2000年提出序列重要性采样算法（Sequential importance sampling，SIS）。该算法可以改善贝叶斯重要性采样计算量大的问题[12]。

首先假设当前状态与将来的观测无关，由此得到重要性概率密度函数的分解式

假设给定系统状态遵循一阶马尔科夫性，且各观测值相互独立，由此得到

（17）（18）两式代入权值递推公式，得到重要性权值的分配公式

（19）式用来递归计算重要性权值，在给定重要性分布函数q(xk|xk-1,zk)的情况下，确保选取时重要性权值的方差最小。如果假设

可以使重要性分布函数满足要求，同时参考分布与真实分布相同，对于粒子设定权值均相等，此时

但是这种假设在实际应用中难以实现，一般都假设

并且代入重要性权值函数，得到

SIS算法通过采用序列的方式，以递归的形式计算重要性权重并进行状态估计。但也正因为如此，随着迭代计数次数的增加，导致样本权值的方差逐级增加，导致权值只体现在少数粒子上，其他粒子权值趋近于零，粒子的多样性严重丢失，这种现象称之为退化现象。避免退化的方式有两种：一是对粒子进行重采样，二是挑选合适的重要性函数。

2.2 SIR

SIR就是当系统检测到有效粒子数低于某个设定的阈值时，说明出现明显退化现象，通过对粒子权值的概率密度函数进行重新采样。换句话说，就是在确保粒子总数不变的前提下，抛弃权值小的非关键性粒子，集中力量计算权值大的重要性粒子，以减缓退化。

采样后的数值近似离散表示为

重采样是对每个粒子按照其权值生产Ni个样本，一旦Ni=0，该粒子被淘汰，但对于权值大的粒子进行有效复制。重采样过程会为当前的粒子带来一些额外的随机变量，从而增加粒子的多样性。

图5为SIR的仿真结果（其中红色曲线为真实值，蓝色为SIR估计曲线），可以看出估计曲线大多时刻可以准确地跟踪真实曲线，仿真效果比较理想。

图5 SIR状态估计曲线

2.3 APF

APF算法是一种基于重采样的粒子滤波的改良算法[13]。APF的核心思想是从SIS的重要性密度函数q(xk,i|z1:k)中生成粒子样本集,其中ij表示的是k-1时刻的标记。然后可以把(xk,i)中的i忽略掉，这样将从边缘密度函数q(xk|z1:k)产生样本集，其重要性密度函数为

假设

代入公式（25），得到

APF算法在原有粒子滤波算法的基础上加入了具有辅助变量的重要性函数q(xk,i|z1:k)，并以此采集样本点。因而，增加了估计的计算量，提高了计算成本。

图6是APF的仿真结果（其中红色曲线为真实值，蓝色为APF估计曲线），跟踪曲线与真实曲线只有较少数的时候存在偏差，估计效果佳。

图6 APF状态估计曲线

3 仿真结果分析

以每次独立仿真实验的均方差误差作为评估几种算法的评价标准

仿真每独立运行一次，每种算法的RMSE都会发生变化，但都在一定的范围内变化。

其中，EKF只能进行1次独立的蒙特卡洛的仿真，其余算法独立运行了100次。各种仿真的RMSE值的柱状图如图7所示。EKF不适合用于非线性程度较高的系统。EKF的值一般在10到40之间，UKF的值在10到25之间，SIR的值在4到12之间，APF的值在3到7之间。

图7 EKF、UKF、SIR和APF的RMSE柱状图

以几种算法的RSME平均值进行比较，详见表1所示。SIR和APF两种粒子算法的RMSE都降至6左右，效果明显优于上两种算法。单独对比这两种粒子算法如表2所示，可以看出不论粒子数量多少，APF的性能都是最好的。

表1 EKF、UKF、SIR和APF的RMSE平均值对比

表2 SIR和APF的RMSE的对比

4 总结

本文探讨了运用EKF、UKF、SIR和APF四种方法估算非线性、非高斯状态系统，并且利用蒙特卡洛仿真和RMSE均方根误差对估算方法的性能进行了评估，总体来说粒子滤波的效果更好。当无论粒子数是50、100还是200，APF都比其他方法有更好的性能。

本次仿真实验实现了大部分目标。但是在过程中发现了两个问题。第一，RMSE评价方法不够直观，因而采用RMSE平均值进行对比。第二，粒子滤波算法的RMSE值是波动的。即使我们使用100个独立的蒙特卡罗仿真并取RMSE平均值，RMSE平均值有时仍是波动的。

在很多的实际应用中，噪声不一定是高斯分布，我们可能不知道系统模型或模型中的参数。这时实际的测量值就很重要，传感器可能受到温度、距离、时间的影响，那么观测噪声可能也随距离、时间或温度的变化而变化。具体的现实问题的要求将指导我们选择更为适宜的方法。