几种非线性滤波算法的研究内附程序

1、2017 年 秋 季学期研究生课程考核考核科目:雷达系统导论 学生所在（系）:电子与信息工程学院学生所在学科:电子与同学工程学 生 姓 名:学 号:学 生 类 别:考核结果阅卷人（读书报告、研究报告） 第 1 页 （共 页）几种非线性滤波算法的介绍与性能分析作者姓名： 学号：专业院系：电信学院电子工程系 电子邮件：摘要非线性滤波算法在雷达目标跟踪中有着重要的应用，对雷达的跟踪性能有着至关重要的影响。好的滤波算法有利于目标航迹的建立及保持，能够得到较精确的目标位置，为发现目标后的后续工作提供可靠的数据依据。本文重点介绍了雷达数据处理中的几种非线性滤波算法：扩展卡尔曼滤波（EKF）、不敏卡尔曼滤波

2、（UKF）、粒子滤波（PF），并且给出了一个利用这三种算法进行数据处理的一个实例，通过这个实例对比分析了这三种算法的性能以及优劣。关键字非线性滤波算法；扩展卡尔曼滤波；不敏卡尔曼滤波；粒子滤波； I. 概述 （一级表题格式）在雷达对目标进行跟踪前要先对目标进行检测。对于满足检测条件的目标就需要进行跟踪，在跟踪的过程中可以利用新获得的数据完成对目标的进一步检测比如去除虚假目标等，同时利用跟踪获得数据可以进一步完成对目标动态特性的检测和识别。因此对目标进行准确的跟踪是雷达性能的一个重要指标。在检测到满足条件的目标后，根据目标运动状态建立目标运动模型，然后对目标跟踪算法进行设计，这是雷达目标跟踪中的

3、核心部分。目前主要的跟踪算法包括线性自回归滤波，两点外推滤波、维纳滤波、滤波、加权最小二乘滤波、维纳滤波和卡尔曼滤波1。对于线性系统而言最优滤波的方法就是卡尔曼滤波，卡尔曼滤波是线性高斯模型下的最优状态估计算法。但是实际问题中目标的运动模型往往不是线性的，因此卡尔曼滤波具有很大的局限性。目前主要用的非线性滤波算法可以分为高斯滤波和粒子滤波2。不敏卡尔曼滤波和扩展卡尔曼滤波就是高斯滤波中的典型代表，也是应用相对较为广泛的。粒子滤波的应用范围比高斯滤波的适用范围要广，对于系统状态非线性，观测模型非高斯等问题都有很好的适用性。本文具体分析阐述了扩展卡尔曼滤波算法，不敏卡尔曼滤波算法，粒子滤波算法，并

4、且通过一个实例利用仿真的方法分析了这三种算法在滤波性能上的优劣，最后对这三种算法做了一定的总结。我本科毕业设计题目为基于历史数据的路径生成算法研究，由于我是跨专业保研到电信学院，该课题所研究内容不属于雷达系统研究范围，是一种城市路网最快路径生成算法。II. 几种非线性滤波算法A. 扩展卡尔曼滤波扩展卡尔曼滤波是将非线性系统转换为近似的线性系统的一种方法，其核心思想是围绕滤波值将非线性函数展开成泰勒级数并略去二阶及以上的项，得到一个近似的线性化模型，然后应用卡尔曼滤波完成状态估计。扩展卡尔曼滤波状态空间模型： 状态方程 观测方程其中和为非线性函数在扩展卡尔曼滤波中，状态的预测以及观测值的预测由非

5、线性函数计算得出，线性卡尔曼滤波中的状态转移矩阵A阵和观测矩阵H阵由f和h函数的雅克比矩阵代替。对和Taylor展开，只保留一次项有：其中：为对求导的雅克比矩阵为对求导的雅克比矩阵，于是可以得出：通过以上变换，将非线性问题线性化。接下来EKF滤波过程同线性卡尔曼滤波相同，公式如下:通过EKF算法线性化状态转移矩阵和观测矩阵后，剩下的滤波过程与普通的卡尔曼滤波无异，滤波过程简单且容易进行。正因EKF简单易于实现的特性，使得该算法一直以来都应用广泛，但是它的局限性也是非常明显的。在这种滤波方法中非线性因子的存在对滤波稳定性和状态估计精度都用很大的影响，其滤波结果的好坏与量测噪声和状态噪声的统计特性

6、也有很大的关系，对于高斯噪声，该算法有很好的适用性，但是对于非高斯噪声，该算法的滤波精度会受到很大的影响。在滤波过程中由于需要预先估计过程噪声协方差矩阵和量测噪声协方差矩阵，如果这两个矩阵的值的估计出现较大的误差，将使得滤波结果出现很大的偏差，容易产生误差积累从而导致滤波结果发散。在用该算法滤波进的初始时刻要预先假设状态的初始值和初始协方差，如果这两个值的估计出现较大偏差，滤波结果也会出现发散的现象。总的来说只有当系统的动态模型和观测模型都接近线性时，利用扩展卡尔曼滤波算法跟踪目标会取得较好的效果，滤波结果较接近真实值。B. 不敏卡尔曼滤波算法不敏滤波也叫无损卡尔曼滤波，它的核心是不敏变换，摒

7、弃了对非线性函数进行线性化的传统做法，采用卡尔曼线性滤波框架，依据不敏变换计算非线性变换的随机变量的统计特性，对于一步预测方程，使用不敏变换来处理均值和协方差的非线性传递。不敏变化不需要像EKF那样对非线性状态方程和量测方程线性化，它是对状态向量的概率密度函数进行近似化，用一系列确定样本来逼近状态的后验概率密度，表现为一些列选取好的采样点，不需要求导计算雅克比矩阵。这些采样点完全体现了高斯密度的真实均值和协方差。近似化后的概率密度函数仍然是高斯的。当这些点经过任何非线性系统的传递后，得到的后验均值和协方差都能够精确到二阶3。由于不需要对非线性系统进行线性化，可以很容易的应用于非线性系统的状态估

8、。对于不敏变换可以做如下阐述：假设X为一个维随机向量，g：为一非线性函数并且。X的均值和协方差分别为。计算UT变化的步骤如下：首先计算式中，是一个尺度参数，可以为任何数值，只要。是均方根矩阵的第i行或第i列，为状态向量的维数，其中。每个采样点通过非线性函数传播，得到计算y的估计均值和协方差，滤波模型如下：由状态方程可以计算得到点的一步预测：状态预测估计和状态预测协方差：式中由量测方程可得点量测的预测值：量测的预测值和协方差为：式中量测和状态向量的互协方差为：状态更新和协方差更新表示为：不敏卡尔曼滤波不必计算雅克比矩阵，不必对非线性系统函数进行任何形式的逼近；在预测阶段只是标准的线性代数运算；对

9、于系统函数来说可以不连续。不敏卡尔曼滤波算法的计算量一般扩展卡尔曼滤波算法，这是由于扩展卡尔曼滤波算法通过线性化处理来实现非线性滤波估计，而不敏卡尔曼滤波是利用样本来逼近状态的概率密度函数，计算量主要发生在选取点时的方根分解运算。在计算速度上扩展卡尔曼滤波算法拥有明显优势，但他的性能随着非线性强度变大而明显下降。不敏卡尔曼滤波算法因不用线性化处理而很好的解决了这一问题。但是不敏卡尔曼滤波算法是用高斯分布来逼近系统状态的后验概率密度，如果系统状态的后验概率密度函数时非高斯的，那么滤波结果将产生极大的误差。C. 粒子滤波算法粒子滤波是一种非线性滤波算法，是一种基于Monte Carlo仿真的最优回

10、归贝叶斯滤波算法4。这种滤波方法将所关心的状态矢量表示为一组带有相关权值的随机样本，并基于这些样本和权值可以计算出状态估值。这种方法不受线性化误差或高斯噪声假定的限制，适用于任何环境下的任何状态转化或量测模型。粒子滤波算法的核心思想便是利用一系列随机样本的加权和表示后验概率密度，通过求和来近似积分操作。滤波模型表述如下：假定k时刻，一组随机样本是根据后验概率密度所获得的采样。其中表示为0到k时刻的第i个样本集合，即粒子集合；为相关权值，并且权值满足，Ns为样本采样数，即粒子数；代表传感器k时刻的量测集合；表示为0到k时刻的所有状态向量集合。则在k时刻，后验概率密度可近似表示为：由于很难直接从抽

11、取样本，通常可利用一个重要性概率密度来获得样本值。从而，权值可以按序贯重点抽样的方法获得。如果是从获得的样本，未归一化的权值可以定义为：如果所选择的重要性概率密度满足：由以上两式可得：为了能够方便的采用回归贝叶斯滤波算法，我们希望重要性概率密度只与前一时刻的测量和状态有关，即：根据以上两式有：在粒子滤波算法中，经过几个迭代周期后，大多数的粒子权值会趋近于零，即粒子衰减现象。由于粒子权值的协方差随着时间的增长而不断变大，这种现象是无法避免的。为了减弱这种影响，一种最直接的方法就是使用大量的粒子数目。当然，这经常是不实际的，因此，目前采用的两种方法：（1）选择最优的重要性概率密度；（2）进行重抽样

12、。最优的重要性概率密度为：最优重要性概率密度可以使采样点权值的协方差最小。目前有两种情况经常采用最优重要性概率密度。第一种情况是Xk为有限集合，例如用于跟踪机动目标的跳变马尔科夫线性系统；第二种情况是状态方程为非线性的，量测方程为线性的系统。重抽样的基本思想是削减小权值的粒子，并集中较大权值的粒子。目前经常采用的集中重抽样方法有：分层抽样和残差抽样、系统重抽样。为了得到正确的状态估计，通常希望粒子权值的方差尽可能趋近于零。然而，蒙特卡洛模拟方法一般都存在权值退化问题。在实际计算中，经过数次迭代，只有少数粒子的权值较大，其余粒子的权值可忽略不计。粒子权值的方差随着时间增大，状态空间中的有效粒子数

13、较少。对于大多数系统来说最优重要性概率密度往往是无法实现的因此经常利用线性化的技术对最优重要性概率密度进行次优近似。粒子滤波算本质上是一种状态搜索算法5。我们估计得到一系列状态，然后将这些状态与我们的观测进行比较，依据似然函数，我们判断哪一些状态是最有可能生成当前观测的状态，然后保留他们，并且给予他们较大的权重，其他的状态责备削弱或者抛弃。但是如果似然函数有若干个窄的峰，那么粒子将有可能聚集到其中的某些峰处，对其他的峰视而不见，也就是说粒子滤波不能保证对任意复杂形状的概率分布都实现准确的描述和跟踪，尤其当空间维数较高时，这一缺陷会更加突出，这归根结底是重采样破坏了粒子的多样性。一个好的重采样算

14、法应该在增加粒子多样性和减少权值较小的粒子数目之间进行有效折衷。III. 仿真分析在本文中我们分别利用扩展卡尔曼滤波算法，不敏卡尔曼滤波算法，粒子滤波算法对一个非线性例子进行仿真分析。系统的状态方程和量测方程都是非线性的。模型如下：实验参数设定6：系统噪声，量测噪声，初始状态，初始方差，时间序列长度T=50,。不敏卡尔曼滤波中的参数为：，。粒子滤波中的粒子数N=300.本文对各种滤波方法进行了M=50次的Monte Carlo仿真，三种算法的滤波结果如图1所示。图1. 三种算法的状态估计曲线图1给出了三种非线性滤波算法的经过50次Monte Carlo仿真后的状态估计。从图中可以看出扩展卡尔曼

15、滤波算法的滤波效果很差，滤波结果大部分时间都偏离真实状态较大。不敏卡尔曼滤波效果相对来说较好，大部分时间滤波结果和真实状态之间的偏差不大。粒子滤波的效果最好，滤波结果和真实状态之间的差距很小，二者具有很好的一致性。三种滤波算法的均方根误差图像如图2所示。图2. 三种算法的均方根误差曲线由图二可以看出扩展卡尔曼滤波算法的均方根误差最大，而且误差数值相对真实值偏离较大，这主要是由于线性化的过程中忽略了高阶项，引入了线性化误差。不敏卡尔曼滤波算法的均方根误差相对较小，但是粒子滤波的均方根误差最小，跟随真实值的效果最好。因此可以看出在非线性条件下，粒子滤波的效果好于不敏卡尔曼滤波好于扩展卡尔曼滤波。表

16、1 各种滤波算法的实验结果比较滤波 算法运算复 杂度运行时间(s) (归一化)存储量非线性 强度EKFO(n3)0.142低弱非线性UKFO(n4)0.220中无限制PFO(Nn3)0.855高无限制表1给出了各种滤波算法的运算复杂度，运算时间，运算所需存储量和滤波对非线性条件的要求。从运算速度上看，扩展卡尔曼滤波算法具有绝对优势，这是由于他将系统方程和量测方程线性化以后的结果。但是根据图1和图2可以看出在非线性较强的条件下这种滤波算法的效果很差，滤波精度较低，甚至出现发散。从滤波精度上看不敏卡尔曼滤波算法和粒子滤波算法差距不大，但是粒子滤波算法的计算量要比不敏卡尔曼滤波算法的大的多，而且算法

17、复杂度是随之粒子数目的增长而线性增加的。对于这两种算法，在一般的高斯白噪声环境下宜使用不敏卡尔曼滤波算法，但是在更复杂的非高斯情况下粒子滤波具有更好的适用性。IV. 结论本文详细的介绍了三种非线性滤波算法：扩展卡尔曼滤波算法，不敏卡尔曼滤波算法，粒子滤波算法。在给出三种滤波算法的详细的理论基础上对每一种算法的优缺点都进行了一定的分析和总结。在文章的最后通过一个仿真实验分析比较了这三种算法在具体滤波过程中的性能。扩展卡尔曼滤波算法在线性化状态方程和量测方程的过程中引入了线性化误差，滤波的精度出现严重下降，滤波结果甚至出现发散；但是线性化后的模型让它的计算量大大地减小了。不敏卡尔曼滤波算法和粒子滤

18、波算法在非线性问题的目标跟踪中表现的都很好，不敏卡尔曼滤波算法的计算量相对来说是比较小的，不过这种算法只适用于高斯白噪声的环境中，如果系统状态的后验概率密度函数时非高斯的，那么滤波结果将产生极大的误差。粒子滤波最复杂的环境有极高的适应性，而且也不要求噪声服从高斯分布，但是该算法的计算量较大，计算量是随着粒子数目线性增加的，而粒子数目往往较大，这就对计算机数据处理性能提出了较高的要求。致谢在本论文的撰写过程中袁子寅师兄给予了很多的建设性意见，宋河晏同学在编程方面提供了一定的帮助，特此表示感谢。参考文献1 秦勤. 雷达目标跟踪的卡尔曼滤波方法的研究D.大连海事大学,2006.2 殷俊丽. 天波超视

19、距雷达机动目标检测与跟踪D.西安电子科技大学,2009.3 夏建涛,任震,陈立,景占荣.极坐标下卡尔曼滤波算法的研究J.西北工业大学学报,2000(03):396-399.4 龚亚信. 基于粒子滤波的弱目标检测前跟踪算法研究D.国防科学技术大学,2009.5 吴兆平. 雷达微弱目标检测和跟踪方法研究D.西安电子科技大学,2012.6 刘丽丽.几种非线性滤波算法的性能分析J.价值工程,2010,29(34):190-191% EKF UKF PF三种算法对比clcclose allclear;% tic;x = 0.1; % 初始状态x_estimate = 1;%状态的估计e_x_estima

20、te = x_estimate; %EKF的初始估计u_x_estimate = x_estimate; %UKF的初始估计p_x_estimate = x_estimate; %PF 的初始估计Q = 10; % 过程状态协方差R = 1; % 测量噪声协方差P =5;%初始估计方差e_P = P; %UKF方差u_P = P;%UKF方差pf_P = P;%PF方差tf = 50; % 模拟长度x_array = x;%真实值数组e_x_estimate_array = e_x_estimate;%EKF 最优估计值数组u_x_estimate_array = u_x_estimate;%

21、UKF 最优估计值数组p_x_estimate_array = p_x_estimate;%PF 最优估计值数组u_k = 1; %微调参数u_symmetry_number = 4; % 对称的点的个数u_total_number = 2 * u_symmetry_number + 1; % 总的采样点的个数linear = 0.5;N = 500; %粒子滤波的粒子数close all;%粒子滤波初始N 个粒子for i = 1 : Np_xpart(i) = p_x_estimate + sqrt(pf_P) * randn;endfor k = 1 : tf% 模拟系统x = line

22、ar * x + (25 * x / (1 + x2) + 8 * cos(1.2*(k-1) + sqrt(Q) * randn; % 状态值y = (x2 / 20) + sqrt(R) * randn; % 观测值%扩展卡尔曼滤波器%进行估计第一阶段的估计e_x_estimate_1 = linear * e_x_estimate + 25 * e_x_estimate /(1+e_x_estimate2) + 8 * cos(1.2*(k-1);e_y_estimate = (e_x_estimate_1)2/20; % 这是根据k=1 时估计值为1 得到的观测值；%相关矩阵e_A =

23、 linear + 25 * (1-e_x_estimate2)/(1+e_x_estimate2)2);% 传递矩阵e_H = e_x_estimate_1/10; % 观测矩阵%估计的误差e_p_estimate = e_A * e_P * e_A + Q;%扩展卡尔曼增益e_K = e_p_estimate * e_H/(e_H * e_p_estimate * e_H + R);%进行估计值的更新第二阶段e_x_estimate_2 = e_x_estimate_1 + e_K * (y - e_y_estimate);%更新后的估计值的误差e_p_estimate_update =

24、e_p_estimate - e_K * e_H * e_p_estimate;%进入下一次迭代的参数变化e_P = e_p_estimate_update;e_x_estimate = e_x_estimate_2;% 粒子滤波器% 粒子滤波器for i = 1 : Np_xpartminus(i) = 0.5 * p_xpart(i) + 25 * p_xpart(i) / (1 + p_xpart(i)2) + 8 * cos(1.2*(k-1) +sqrt(Q) * randn;p_ypart = p_xpartminus(i)2 / 20; % 预测值p_vhat = y - p_y

25、part;% 观测和预测的差p_q(i) = (1 / sqrt(R) / sqrt(2*pi) * exp(-p_vhat2 / 2 / R); % 各个粒子的权值end% 平均每一个估计的可能性p_qsum = sum(p_q);for i = 1 : Np_q(i) = p_q(i) / p_qsum;% 各个粒子进行权值归一化end% 重采样权重大的粒子多采点，权重小的粒子少采点, 相当于每一次都进行重采样；for i = 1 : Np_u = rand;p_qtempsum = 0;for j = 1 : Np_qtempsum = p_qtempsum + p_q(j);if p_

26、qtempsum = p_up_xpart(i) = p_xpartminus(j); % 在这里xpart(i) 实现循环赋值；break;endendendp_x_estimate = mean(p_xpart);u_x_par = u_x_estimate;for i = 2 : (u_symmetry_number+1)u_x_par(i,:) = u_x_estimate + sqrt(u_symmetry_number+u_k) * u_P);endfor i = (u_symmetry_number+2) : u_total_numberu_x_par(i,:) = u_x_es

27、timate - sqrt(u_symmetry_number+u_k) * u_P);end%计算权值u_w_1 = u_k/(u_symmetry_number+u_k);u_w_N1 = 1/(2 * (u_symmetry_number+u_k);%把这些粒子通过传递方程得到下一个状态for i = 1: u_total_numberu_x_par(i) = 0.5 * u_x_par(i) + 25 * u_x_par(i)/(1+u_x_par(i)2) + 8 * cos(1.2*(k-1);end%传递后的均值和方差u_x_next = u_w_1 * u_x_par(1);f

28、or i = 2 : u_total_numberu_x_next = u_x_next + u_w_N1 * u_x_par(i);endu_p_next = Q + u_w_1 * (u_x_par(1)-u_x_next) * (u_x_par(1)-u_x_next);for i = 2 : u_total_numberu_p_next = u_p_next + u_w_N1 * (u_x_par(i)-u_x_next) * (u_x_par(i)-u_x_next);endfor i = 1 :u_total_numberu_y_2obser(i,:) = u_x_par(i);e

29、nd%通过观测方程得到一系列的粒子for i = 1: u_total_numberu_y_2obser(i) = u_y_2obser(i)2/20;end%通过观测方程后的均值y_obseu_y_obse = u_w_1 * u_y_2obser(1);for i = 2 : u_total_numberu_y_obse = u_y_obse + u_w_N1 * u_y_2obser(i);end%Pzz测量方差矩阵u_pzz = R + u_w_1 * (u_y_2obser(1)-u_y_obse)*(u_y_2obser(1)-u_y_obse);for i = 2 : u_tot

30、al_numberu_pzz = u_pzz + u_w_N1 * (u_y_2obser(i) - u_y_obse)*(u_y_2obser(i) - u_y_obse);end%Pxz状态向量与测量值的协方差矩阵u_pxz = u_w_1 * (u_x_par(1) - u_x_next)* (u_y_2obser(1)-u_y_obse);for i = 2 : u_total_numberu_pxz = u_pxz + u_w_N1 * (u_x_par(i) - u_x_next) * (u_y_2obser(i)- u_y_obse);end%卡尔曼增益u_K = u_pxz/u_pzz;%估计量的更新u_x_next_optimal = u_x_next + u_K * (y - u_y_obse);% 第一步的估计值+ 修正值；u_x_estimate = u_x_next_optimal;%方差的更新u_p_next_update = u_p_next - u_K * u_pzz * u_K;u_P = u_p_next_update;%进行画图程序x_array = x_array,x;e_x_estimate_array = e_x_estimate_array,e_x_estimate;p_x_estimate_array