弹道再入目标轨迹跟踪中非线性滤波算法研究＊

田亚菲 莫 骅 于 飞

（1．兰州大学信息学院 兰州 730000）（2．总参谋部通信训练基地 张家口 075100）

1 引言

弹道再入目标的轨迹跟踪问题一直以来受到人们的普遍关注，它是弹道导弹防御系统的核心技术，直接影响对来袭导弹拦截的成功率。针对跟踪目标的动力学模型和雷达的观测模型具有典型的非线性特性，国内外相关文献提出许多非线性滤波算法，在算法选择上存在明显不同的观点［1～3］。常见的算法中，大致可分为两大类［4］：一类是基于点的滤波算法，主要包括扩展卡尔曼滤波器（EKF）和无迹卡尔曼滤波器（UKF）及对应的改进算法；另一类是基于密度的滤波算法，包括贝叶斯粒子滤波器（PF）及对应的改进算法。

本文依据弹道再入目标动力学模型和雷达观测模型特点，从算法基础理论出发对比分析EKF和UKF算法的异同点，通过仿真实验全面比较了EKF、UKF和PF三种基本算法的性能，得出到在高斯噪声条件下，EKF要优于UKF和PF，在闪烁噪声条件下，PF要优于EKF和UKF的结论。

2 基于点的滤波算法理论

2．1 线性递归估计算法的基本思想

本质上讲，EKF和UKF在状态更新阶段是一个基于最小均方误差（MMSE）条件下的线性递归估计算法。一个非线性动态系统的状态方程和观测方程可用下面两式来描述。

使用线性回归估计方法可以从当前观测值zk中估计出k时刻的状态变量xk，其中经典的算法是线性最小均方误差估计法（LMMSE）。观测值序列的概率密度函数未知的平稳过程在前二阶矩已知的条件下，可以通过递归迭代的算法完成估计，此时基于观测值的估计是观测值的线性变换。若符合高斯约束条件时，LMMSE算法是基于MMSE约束条件下的最优估计。

用zk表示从初始时刻到k时刻所得到的一组观测数据，LMMSE 估计中E＊xk｜zk具有如下解析形式［5］：

如果以上式中状态一步预测、一步预测误差协方差，观测值预测、观测预测估计协方差Sk，状态估计值与观测估计值互协方差都能通过给定的k-1时刻的状态估计值和误差协方差Pk-1｜k-1（即初值）求出，此LMMSE算法就是真正意义上的递归运算。若系统的动力学方程和观测方程能满足线性和高斯条件，以上计算公式就是著名的kalman滤波算法。

对于具有观测方程（2）的非线性系统

和都依赖于k-1时刻以前的观测值zk-1和相应时刻LMMSE 估计值。为实现真正的递归线性估计，要求k时刻的估计值仅从和Pk-1｜k-1就能求得，即：

以上两式中，Pred f（·），，Pk-1｜k-1表示｛，Pk-1｜k-1｝通过非线性函数f（·）的传递作用来近似得到一步状态预测值和对应的协方差Pk｜k-1。Pred h（·），，Pk｜k-1表示 通过状态预测值来近似得到观测预测值和相应的误差矩阵。近似运算由于舍弃了部分信息，从而导致信息传递的不完整，使非线性滤波器的性能发生退化。通常采用两种方法来克服近似运算带来的弊端：一是尽可能地使用最优的近似，这样就能避免过多的舍弃信息，保证预测的准确性，例如采用高阶的泰勒级数展开方法（TSE）；另外一种方法是用一种适当的分布来逼近非线性函数而不进行舍弃，这样的预测中包含了所有非线性系统的统计信息，例如无迹变换方法（UT）［6］。

无论观测方程是否满足线性条件，只要能实现非线性系统的近似线性化，得到预测值的一阶矩和二阶矩信息，LMMSE估计器就可以非常理想地直接应用在非线性系统的状态更新上，这就是EKF和UKF本质上相同之处。两者的区别仅在于在式（7）～式（8）中对非线性系统的近似处理上采用不同的策略。

2．2 EKF的线性化策略

EKF的近似策略是对非线性函数的Taylor展开式进行一阶或二阶的线性化截断，忽略高阶项。假设过程噪声和量测噪声是0均值且互不相关的的高斯噪声，在一阶EKF中，非线性系统的状态方程（1）和观测方程（2）通过以下的步骤完成线性化：

F和H分别是函数f（·）和h（·）对状态变量x求偏导的雅克比矩阵，即和W和V分别是函数f（·）和h（·）分别对过程噪声w和量测噪声v求偏导数的雅克比矩阵，即W［i，j］＝根据以上线性化变换，预测及状态更新过程描述如下

式（11）～式（15）与LMMSE 的递归线性运算式共同构成了一阶EKF算法公式。

EKF的估计精度严重依赖于系统的非线性程度。为提高近似精度，可采用二阶EKF。虽然二阶EKF可以比一阶EKF获得更好的滤波性能，但是计算量较大，同时要求非线性函数的海森矩阵在数学上存在，而在计算海森矩阵中会产生计算误差。因此，针对弹道再入目标观测方程的强非线性特性，在EKF线性化近似之前，可以采用坐标系转换的方式来尽可能地减弱其非线性特性［7］。对于极坐标下所描述的观测方程，可以通过以下两式转化成直角坐标系的方程［8］。

转化后的观测变量成为伪线性的形式，减弱了雷达观测变量在初始极坐标系下的非线性程度。

2．3 UKF的线性化策略

在基于近似非线性函数的概率分布比近似其本身更容易［9～10］的思想指导下，Julier等提出了基于unscented transformation的Kalman滤波算法，即UKF。x为nx维的随机变量，y＝g（x）是一个非线性变换，假设x的均值和方差分别表示为和Cx，算法描述如下：

第一步：计算sigma点的值χi及其权重Wi：

其中κ是缩放比例系数，是矩阵（nx＋κ），Cx的方根的第i行（列）。

第二步：通过非线性函数传播sigma点：

第三步：计算y的均值和方差，以及x和y的协方差：

非线性系统的状态方程（1）和观测方程（2）通过UT 完成线性化变换：

和分别表示在k-1 时刻预测值和更新值，和是在k-1时刻过程噪声和观测噪声。状态更新阶段同EKF一样使用LMMSE估计算法。

尽管采用有限个采样点简化了近似过程，但是计算采样点的方根运算增加了计算量，且根式中的矩阵一旦失去正定势必导致滤波器发散。

3 基于密度的PF算法理论

Gordon和Salmond在1993年提出了一种基于序列重要抽样滤波思想的Bootstrap非线性滤波方法［11］，从此奠定了粒子滤波（PF）算法的基础。该算法的基本原理是对于一个平稳的时变动态系统，设k-1时刻系统的后验概率密度为p（xk-1｜zk-1），根据一定的原则选择n个随机样本点（粒子），k时刻获得观测信息后，经过预测更新和状态更新，n个粒子的后验概率密度近似为p（xk｜zk），随着粒子数目的增加，粒子的概率密度函数就逐渐逼近状态变量的概率密度函数，此时粒子滤波器的估计就达到了最优贝叶斯估计的效果［12］。

基本的粒子滤波算法流程如下：

第一步初始化：取k＝0，按p（x0）抽取N个样本点，i＝1，2，…，N；

第四步归一化权重：

第五步重采样：根据各自归一化权重的大小复制或者舍弃样本，得到N个近似服从p分布的样本。

进博会食品及农产品展区6万平方米，分为乳制品、蔬果农产品、肉制品和水产品、休闲食品甜食调味品、酒类和饮料等五大专业展区。徜徉其间，从参观者身上你会感受到中国对美好生活的追求和向往是那么强烈！中国当自强！中国农业当自强！

第七步k＝k＋1，重复第二至六步进行迭代运算直到结束。

PF算法通过寻找一组在状态空间传播的随机样本来对概率密度函数完成近似，以样本均值代替积分运算，从而获得状态变量最小方差分布。它与EKF、UKF不同，对状态变量没有高斯约束，因此近年来成为了非线性非高斯系统状态估计的“最优”估计器。但是PF算法滤波精度与粒子数量、重要性函数、重采样方法的选取有很大关系［13］。近年来国内外出现了大量的改进算法，但计算量巨大、实时性差的缺点始终未能克服，在一定程度上限制了其在弹道再入目标轨迹跟踪问题上的工程应用。

4 仿真研究

4．1 弹道再入目标的数学模型及仿真环境描述

在不影响模型精确度的前提下，假定地球表面为平面，目标为质点且不考虑自旋。相关文献对弹道再入目标的数学模型进行了较为系统的分析［7，14］。本文用如下两式描述动力学方程和观测方程。

xk＝［xk，˙xk，yk，˙yk］′是状态变量，各元素分别表示目标在直角坐标系下x轴和y轴方向上的位置和速度分量。zk＝［rk，θk］′是观测变量，各元素分别表示雷达对目标的测量距离和角度。wk是过程噪声，代表观测噪声，各元素分别是距离和角度的观测误差。假设过程噪声和观测噪声互不相关，且服从不依赖于初值x0的0均值高斯分布，它们的协方差矩阵通过以下三个矩阵式描述，其中Q是过程噪声协方差矩阵，q是过程噪声的方差，R是观测噪声协方差矩阵。

式（18）中非线性函数Ψ（·）定义为

式（20）右端后两项描述的是空气阻力和重力对目标状态变量的影响作用。g是重力加速度，f（xk，β）是空气动力阻力函数，β为弹道系数，一般的，该系数是时间变量t的函数，即β（t）。文献［15］讨论了参数β（t）的建模方式对非线性滤波算法的性能影响。为简化模型，本文采用中的弹道系数为确知参数。式（20）中的状态转移矩阵F和驱动矩阵G分别表示为

空气动力阻力函数f（xk，β）定义为式（21）。

ρ（·）是随海拔高度呈指数衰落的空气密度函数。观测方程中的非线性测量函数h（·）表示为

仿真中雷达扫描周期T＝2s，仿真跟踪步数60步，目标弹道系数β＝40000kg／ms2，蒙特卡洛仿真次数为1000次，目标初始状态［332000m，2290cos（190°）m／s，98000m，2290sin（190°）m／s］，高斯噪声环境下观测噪声的标准差分别是σr＝100m 和σε＝0．05rad，过程噪声方差q＝2。

在高斯噪声条件和闪烁噪声条件下，分别使用EKF、UKF及PF进行滤波，通过比较不同方向上位置、速度的估计精度和运算时间来对比三种算法的基本性能。EKF采用坐标变换的方法减弱观测方程的非线性特性，使用一阶近似，PF 中粒子数为500，采用残差重采样方法。

4．2 仿真结果及分析

图1 三种算法在X 轴上位置的RMSE

图2 三种算法在X 轴上速度的RMSE

图3 三种算法在Y 轴上位置的RMSE

图4 三种算法在Y 轴上速度的RMSE

图1比较了高斯条件下三种滤波算法在X和Y轴方向上位置与速度的均方根误差（RMSE），由于各滤波器的初值已给定，从初始时刻开始三种滤波器就能得到良好的跟踪效果。图1中，EKF 和PF对X轴方向上位置变量的跟踪精度几乎一致，要略优于UKF。从理论上讲，如果适当增加PF的粒子数量和尝试不同的重采样算法，会提高其跟踪精度，但是会带来运算量的剧增。在仿真的第60步，三种滤波器的结果趋于稳定状态，分别收敛于60m、70m 和120m 的误差水平。图2表明，EKF和PF对X方向上速度估计的性能要优于UKF，误差仅是后者的1／3左右。图3中尽管UKF在第40步左右出现190m 的误差，但是三种滤波器在第50步后，误差稳定于160m，可以认为三种滤波器的性能相当。图4可以得到和图2基本相同的结论。

通过对仿真结果的分析可看出PF 由于在粒子数量和重采样方法选择策略两个条件限制下，并未能达到理想的跟踪精度，但仍体现出良好的性能，基本与EKF 的性能相当。相同条件的仿真中EKF则要优于UKF，这并非与文献［10］的论述相悖，只是由于在弹道再入目标轨迹跟踪问题上，若弹道系数预先确知，对观测方程采用适当的伪线性变换，削弱其非线性程度，EKF 是可以达到了较为精确的滤波效果的。理论上讲，UKF 能够得到比EKF更高的近似精度，而并非在任何情形下都要优于EKF。UT 作为一种线性化的技术，只是为了获得滤波器状态更新时必要的信息，而并不能完全取代其它的线性化技术。

表1给出了MATLAB 中测量的三种滤波器相对于EKF在一次蒙特卡洛仿真中的运算时间。

表1 三种算法的相对运行时间

三种滤波算法中，UKF的运算时间为EKF 的1．85 倍，运行时间与EKF 相当，估计精度不及EKF。PF尽管在运行精度上和EKF相当，运行时间却是其91．05倍。

表2比较了闪烁噪声环境下三种滤波算法的X 和Y 轴方向上位置和速度的均方根误差的均值（ARMSE）。闪烁噪声采用高斯噪声和具有“厚尾”特性的拉普拉斯分布噪声加权和来实现。在严重闪烁噪声干扰条件下，仿真程序中闪烁噪声权重Weight＝0．9。

表2中，闪烁噪声条件下EKF 出现了极大偏差，UKF发散，滤波器无法正常工作，因此基于点的滤波器在此条件下失去了对目标的跟踪性能。PF的ARMSE尽管相对于高斯环境下略有升高，但仍然收敛。

表2 三种算法在闪烁噪声环境下的ARMSE

5 结语

本文基于弹道再入目标的轨迹跟踪问题，比较了EKF、UKF 和PF 的性能。通过分析LMMSE滤波器的基本原理，可知所有基于点的滤波算法在状态更新阶段的本质是相同的，仅在线性化策略上有所区别。尽管在高斯噪声条件下非线性滤波问题中，UKF要比EKF 具有更高阶的信息，其性能优于EKF，但在弹道再入目标轨迹跟踪问题中，如果目标的弹道系数确知，EKF 具有更高的估计精度和更短的运算时间，能够接近优化后的PF 算法的性能。因此高斯条件下在滤波器的选择上本文更倾向于EKF 算法。闪烁噪声条件下，基于高斯假设的基于点的滤波算法失去了理论基础，滤波器性能出现偏差甚至是发散，此时基于密度的PF 算法就成为了“最优”的估计器，同时对弹道再入目标滞空时间短暂与PF 算法运算时间长的矛盾应该予以重视。

应当指出的是，初始状态对非线性滤波器的稳定性非常关键，弹道再入目标的初始状态的确定方法对滤波算法的性能会产生很大影响，这个问题有待进一步研究。当弹道系数作为状态变量的一个待估参量时，其建模方法对非线性滤波器的性能会产生影响。下一步将继续研究这些因素对弹道目标轨迹跟踪问题中非线性滤波器选择策略的影响。

［1］Farina A，Ristic B，Benvenuti D．Tracking a Ballistic Target：Comparison of Several Nonlinear Filters［J］．IEEE Transactions on Aerospace and Electronic Systems，2002，38（3）：854-867．

［2］江晓东，谢京稳，郭军海．三种非线性滤波在再入弹道估计中的分析研究［J］．战术导弹技术，2012（1）：43-49．

［3］刘鑫，郭云飞，薛安克．几种滤波器跟踪性能的比较［J］．计算机与数字工程，2013（3）：391-393．

［4］Zhao Z L，Chen H M，Chen G S，et al．Comparison of Several Ballistic Target Tracking Filters［M］．New York：IEEE，2006：2197-2202．

［5］Zhao Z L，Li X R，Jilkov V P．Best Linear Unbiased Filtering with Nonlinear Measurements forTarget Tracking［J］．Aerospace and Electronic Systems，IEEE Transactions on，2004，40（4）：1324-1336．

［6］Li X R，Jilkov V P．A Survey of Maneuvering Target Tracking：Approximation Techniques for Nonlinear Filtering［C］／／Pefense and Security．Intemational Society for Optics and Photonics，2004：537-550．

［7］Li X R，Jilkov V P．A Survey of Maneuvering Target Tracking—PartⅢ：Measurement Models［C］／／Internation Symposium on Optical Science and Technology．International Society for Optics and Photonics，2001：423-446．

［8］Yaakov B S，Li X R．Estimation with Applications to Tracking and Navigation：Theory Algorithms and Software［M］．New York：John Wiley ＆Sons，Inc．，2001．

［9］Julier S J，Uhlmann J K．Unscented Filtering and Non-linear Estimation［J］．Proceedings of the IEEE，2004，92（3）：401-422．

［10］Julier S J，Uhlmann J K，Whyte H F．A New Method for the Nonlinear Transformation of Means and Covariances inFilters and Estimators［J］．Automatic Control，IEEE Transactions on，2000，45（3）：477-482．

［11］Gordon N J，Salmond D J，Smith A F M．Novel Approach to Nonlinear／non-Gaussian BayesianState Estimation［J］．Iee Proceedings F Radar and Signal Processing，1993，140（2）：107-113．

［12］Doucet A，Freitas N，Gordon N．Sequential Monte Carlo Methods in Practice［M］．New York：Springer-Verlag，2001．

［13］冯驰，王萌，汲清波．粒子滤波器重采样算法的分析与比较［J］．系统仿真学报，2009（4）：1101-1105．

［14］Li X R，Jilkov V P．A Survey of Maneuvering Target Tracking-partⅡ：Ballistic Target Models［C］／／Proceedings of 2001SPIE Conference on Signal and Data Processing of Small Targets，2001，4473：559-581．

［15］Ristic B，Farina A，Benvenuti D．Tracking a Ballistic Object on Reentry：Performance Bounds and Comparison of NonlinearFilters［C］／／Radar，Sonar and Navigatio，IEE Proceedings，2003，150（2）：65-70．