范 旭，陈国光,王光志，赵军强，杨智杰

(1 中北大学 机电工程学院， 太原 030051； 2.豫西工业集团有限公司， 河南 南阳 473000)

【弹道工程和火药工程】

卡尔曼弹道滤波器抗野值方法研究

范 旭1，陈国光1,王光志2，赵军强2，杨智杰1

(1 中北大学 机电工程学院， 太原 030051； 2.豫西工业集团有限公司， 河南 南阳 473000)

在弹道诸元参数最优估计获取过程中,为减少量测过程中可能会形成的野值对卡尔曼滤波过程所造成的不利影响，提出了一种基于多项式拟合的野值去除方法。该方法通过对滤波前期弹道数据的选取，进行拟合估计，获得三方向位置拟合多项式。在此多项式的基础上，以残差作为判别依据，以三倍标准差作为判别标准，将病态点去除。仿真结果表明，该方法行之有效，可将野值去除，保证滤波快速收敛。

弹道滤波；卡尔曼滤波；野值去除；多项式拟合

对于弹道测量，其测量数据中包含着干扰或噪声引起的误差，直接用测量数据进行分析计算将会带来很大的误差，甚至是错误，而必须从量测数据中滤取所需要的信息，这种数据处理方法称为滤波。目前工程技术中使用较多的是最小二乘滤波、最大似然估计和卡尔曼滤波，其中卡尔曼滤波在控制领域，在火箭、炮弹、导弹弹道和飞行状态以及卫星轨道观测的数据处理上得到广泛的应用。卡尔曼最初提出的滤波基本理论只适用于线性系统，并且要求量测也必须是线性的。如果模型是非线性的，通常在推导滤波方程时，增加线性化步骤。在状态估计时，对系统方程在前一状态值处做实时的线性泰勒近似；在预测步骤中，对量测方程在相应的预测位置也进行线性泰勒近似；所得到的卡尔曼滤波被称为扩展卡尔曼滤波。处理非线性模型的这一思想是很自然的，且滤波过程简单有效[1]。

而量测系统的噪声是真实存在且不可忽视的。在常规的滤波过程中，我们视系统噪声为正态分布的高斯噪声，虽然正态分布的规律使大部分的随机数位于三倍标准差之内。但在实际情况中，量测信号要经过传输系统、数据采集转换系统等环节，才能转化为数字形式的目标信号，这样不但会带有大量的随机噪声而且有可能受到外界无规律的干扰产生一些与真实数据相差较大的野值，以至于影响滤波过程。为避免这一现象，对滤波过程的抗野值方法的研究势在必行[2]。

1 外弹道滤波模型

弹道滤波即是根据在一段弹道上测得的弹箭飞行弹道数据(坐标，速度、加速度、姿态角和姿态变化率)，利用数学方法从这些数据中提取当前飞行状态及弹道数学模型中的参数的最优估计。而弹道预测则是利用这种最优估计和弹道数学模型计算弹道，求取在所需点的弹道参数。弹道滤波是弹道预测的基础，而弹道预测是弹道修正的重要依据，所以滤波的收敛速度及准确度也是弹道修正的重要研究方向。

1.1 滤波状态模型

地面炮位侦查雷达对空中飞行的弹箭进行跟踪和探测，地面计算机对雷达实测到的弹道参数进行数据处理，获得较准确的弹道参数，与目标坐标比较，得到弹道偏差，并由雷达向弹载计算机传输指令，使修正控制机构实现一维或二维弹道修正。根据上述弹道修正弹的执行过程，在雷达量测数据的滤波的过程中，需要用到弹道模型，考虑弹道参数获取的快速性、实时性等，本文采用质点弹道模型作为弹道滤波的状态方程[3]：

(1)

则有系统状态方程：

(2)

式中W(k)是均值为零的高斯白噪声，且服从方差为Q的正态分布，W～N(0,Q)。

1.2 滤波量测方程

设雷达测量值为斜距r、方位角β和高低角ε，雷达坐标系为球坐标系，可得雷达测量值与地面坐标系的关系[4]：

β=arctan(z/x)

(3)

令量测变量为Z，即Z=(r,β,ε)T，则得量测方程为:

(4)

式中，V是雷达的量测噪声，假定为零均值高斯白噪声，且服从方差为R的正态分布，即V～N(0,R)，h(x)为三维矢量函数。

2 卡尔曼弹道滤波方程

卡尔曼滤波方程组分为两大部分，第一部分为卡尔曼滤波方程，这一部分负责向下一时刻推算轨迹状态；第二部分为卡尔曼滤波器的增益矩阵递推公式，这一部分用于反馈先验估计，并对预测进行修正[4]。

设k时刻的被估计状态Xk受系统噪声序列w驱动，驱动机理由下述状态方程描述

Xk=Φk,k-1Xk-1+Γk-1Wk-1

(5)

量测方程：

Zk=HkXk+Vk

(6)

式中:Φk,k-1为tk-1时刻的一步转移矩阵；Γk为系统噪声驱动阵；Hk为量测阵；Vk为量测噪声序列；Wk为系统激励噪声序列。

如果被估计状态Xk满足式(5)，对Xk的量测Zk满足式(6)，系统噪声Wk和量测噪声Vk满足式(6)，系统噪声Wk方差阵Qk非负定，量测噪声Vk方差阵Rk正定，k时刻的量测为Zk，则Xk的估计Xk按下述方程求解：

状态的一步预测

(7)

量测预测方程：

(8)

(9)

滤波增益：

(10)

一步预测均方误差矩阵：

(11)

估计均方误差

Pk=(I-KkHk)Pk/k-1

(12)

3 野值的判别与处理

3.1 三方向量测位置分量的方差计算

设弹道时间为t，将量测数据rt,βt,εt由下式转化为x,y,z三方向数据即xt,yt,zt。

xt=rt·cosβt·cosεt

yt=rt·sinεt

zt=rt·cosεt·sinβt

(13)

(14)

同理可求得δy，δz。

3.2 弹道轨迹的多项式拟合

在实际情况下，弹体的在空中的运动受到重力，空气阻力及横风等外界因素影响，其运动轨迹为曲线[8]。在曲线拟合的过程中，可选择线性回归模型，指数类或幂指数类等，及其组合形式作为拟合算法的基础，其中线性回归模型即多项式拟合方法计算复杂度低，计算速度快。同时，分析真实的弹道曲线，可发现弹道曲线形状与变化规律与抛物线类似，从实际情况与模型的相似性出发，可选择多项式进行曲线拟合。从上述两点分析，将采用线性回归模型进行弹道拟合并分析。在曲线拟合的计算过程中，常用到最小二乘法，高斯消去法，三次样条法和追赶法等方法，来对系数矩阵进行求解，进而得到多项式系数。基于弹道曲线的形状和拟合弹道多项式计算的快速实时性，选取最小二乘法进行进行弹道多项式拟合。

设拟合多项式的为n次幂多项式，其表达形式为:

f(x)=A1·xn+A2·xn-1+…+An·x+An+1

(15)

取m次量测数据(xi,yi)(m>n)，根据最小二乘法估计准则可得目标函数：

(16)

因目标函数为实际问题，其极小值必然存在，由极小值存在条件：

(17)

将上述条件方程转化为矩阵表达形式如下：

(18)

由数值分析方法可知，当系数矩阵A非奇异时，则上式有唯一解向量x。x向量中的参数即A1、…、An+1，就是多项式各次幂的系数。在已知完整多项式及量测时间间隔T的情况下， 将xm=t+T代入，便可获得预估位置分量xm+1,ym+1,zm+1。

3.3 野值的判别与处理

在绿波起始时刻t之前，分别取基准弹道数据和量测数据，数据长度为N。将量测数据段代入按上式(14)，求得三方向标准差δx，δy，δz。将量测数据与基准弹道数据比对，若三方向上任意方向上差值大于3δ，则将该数据点去除，即：

用剩余数据点进行多项式拟合，进而得到多项式全系数。

滤波开始后，利用基准多项式估计下一采样时刻数据，以上述原则为判据。若量测数据在判据下非野值，则将其用于滤波，并替代数据段中最初时刻数据，重新拟合多项式用以下次估计。

若量测数据被判别为野值，则不将此点代入滤波过程，以多项式估计值作为新量测数据，

代替原始量测，进入滤波过程且多项式无需改变。

4 仿真校验

为验证野值去除方法的可行性和实用性，需进行计算机仿真计算。以某型火箭弹为例，通过6D弹道方程求解出弹道数据，将弹道数据转化为极坐标的雷达量测数据。并通过随机发生器产生均值为零的高斯白噪声作为要引入的雷达噪声，加入到原始量测数据，不定时间间隔加入较大误差作为野值，产生受噪声影响的雷达量测数据。

在仿真计算中，取雷达数据刷新间隔为0.1 s，径向距离上的方差σr为10 m，方向角的噪声方差σβ为0.001 5 rad，俯仰角的噪声方差σε为0.001 5 rad，即量测噪声矩阵R的主对角元素值。模拟在火箭弹发射3.5 s后雷达开始跟踪目标测量，在量测数据在40 s，开始对弹道数据进行滤波处理，持续监控20 s左右至60 s时停止滤波处理。

多项式的拟合精度受拟合数据点个数及拟合模型与实际轨迹契合度等方面影响。所以将从采用数据点的时间长度，多项式阶数两方面来分析其变化对多项式拟合精度的影响。

从图1可看出所估位置的残差随着拟合多项式阶数的增加先减少而后略有增加。虽然拟合数据时间长度不同，但一次多项式基本对应残差最大值，造成这种情况的原因在于直线模型与真实三方向轨迹情况存在很大差异，具有较大的模型误差，拟合效果差。

图1 多项式阶数与拟合位置残差的关系

由一次多项式至二次多项式时，拟合残差迅速减小，证明拟合模型与真实轨迹状况的匹配度较高。图2表示了二次多项式拟合残差与拟合数据长度的关系。

由二次多项式至五次多项式时，拟合残差却略有增长，原因在于拟合模型与真实轨迹匹配度逐渐减弱。基于精度考虑，将选取二阶多项式作为拟合模型。

图2 二次多项式拟合残差与拟合数据长度的关系

从图2可看出随数据点数的增加，基于二次多项式的拟合残差逐渐减小，而到使用滤波前五秒数据时残差却有所增长。其原因在于，当数据时间长度在一秒和四点五秒之间时，随着数据点的增多，真实弹道与抛物线的相似度逐渐提高，使模型误差逐渐减小，残差随之下降。但真实弹道虽与抛物线相似，但其轨迹为非对称性，并不是二次多项式完全契合的抛物线，所以当数据时间长度再次增加时，这种差异性会逐渐凸显，造成拟合残差增大。所以拟合数据时间长度选取为4.5 s。

从图3、图4和图5可看出存在野值时，若不去除，滤波过程波动较大，滤波收敛速度受到影响，大幅度了延长了滤波收敛时间。可见野值的存在对滤波过程的影响持续时间较长且波动较大。证明去除野值应是弹道滤波过程必须考虑到的因素。

图3 X方向是否去野值对滤波的影响

图4 Y方向是否去野值对滤波的影响

图5 Z方向是否去野值对滤波的影响

5 结论

本文针对雷达量测的火箭弹弹道数据，通过对滤波前期的量测数据进行坐标转换，对所量测的三方向数据进行多项式拟合，进而估计滤波下一时刻的位置分量。以滤波前量测数据作为样本，求取三方向位置分量方差，以此来作为量测系统的总体参数。选取拟合数据长度和拟合多项式阶数作为变量，通过计算机仿真分析其对弹道参数拟合精度的影响，最后选取最优组合确定其为弹道拟合模板。通过滤波仿真，验证此种野值去除方法的可行性。由仿真结果可看出，文中所提出的去野值方法有效，简单，计算速度快，符合实用性要求，可有效去除野值对滤波过程的影响，保证滤波过程平缓，无较大波动，使滤波快速收敛。

[1] 秦永元,张洪钺,汪叔华.卡尔曼滤波与组合导航原理[M].西安：西北工业大学出版社，1998:33-41.

[2] 闫小龙,陈国光,杨东.抗野值卡尔曼滤波在火箭弹落点估计中的应用[J].弹箭与制导学报,2016(3):94-98.

[3] 韩子鹏.弹箭外弹道学[M].北京：北京理工大学出版社，2014:74-77.

[4] 李岩,任睿,王旭刚.两种卡尔曼滤波模型在修正弹弹道数据处理中的应用比较[J].弹道学报,2011(1):27-30.

[5] 吴汉洲,宋卫东,徐敬青.基于多项式拟合的扩展卡尔曼滤波算法[J].计算机应用,2016(5):1455-1457,1463.

[6] 闫小龙,陈国光,白敦卓.自适应卡尔曼滤波器在火箭弹落点估计中的应用(英文)[J].Journal of Measurement Science and Instrumentation,2015(3):212-217.

[7] 余小琴,沈文苗.扩展卡尔曼滤波算法初值选取方法[J].声学与电子工程,2012(1):12-13,17.

[8] 秦伟,田晓丽,刘超.基于卡尔曼滤波的导弹姿态最优控制法[J].机械工程与自动化,2016(1):183-185.

[9] 战杰. 一种外弹道测量数据的斑点型野值剔除方法[J].航天控制,2016,34(1):75-77+83.

[10] 卓宁. 靶场外弹道数据处理中野值点剔除方法[J].测试技术学报,2008(4):313-317.

[11] 李景熹,王宇,王树宗,等.观测值中野值的判别与处理方法仿真研究[J].微计算机信息,2006(13):140-142.

TheStudyofKalmanBallisticFilterforAdaptingtheOutliers

FAN Xu1, CHEN Guoguang1, WANG Guangzhi2, ZHAO Junqiang2, YANG Zhijie1

(1.The North University of China, Taiyuan 030051, China; 2.Yuxi Industries Group Co.,Ltd., Nanyang 473000, China)

In the acquisition process of the optimal ballistic parameter estimation,Kalman filter is an algorithm that the application and usability of project is strong.To reduse the adverse effects of the outliers that might be formed during the measurement on Kalman filtering process, a kind of wild value removed method based on polynomial fitting is proposed. The method use the early trajectory data of the filter to polynomial fit.On the basis of this polynomial, the residual as the criterion,triple std as the judge standard, wipe off the outliers. Simulation results show that,the method is effective,the outliers can be removed, the process of filtering can be fast convergence.

trajectory filter; Kalman filter; outliers removed; polynomial fit

2017-07-20；

2017-08-11

范旭，男，硕士研究生，主要从事智能弹药、弹箭制导与控制研究; 陈国光，男，博士，教授，博士生导师，主要从事信息化武器系统研究。

10.11809/scbgxb2017.12.022

本文引用格式：范旭，陈国光,王光志，等.卡尔曼弹道滤波器抗野值方法研究[J].兵器装备工程学报，2017(12):93-97.

formatFAN Xu, CHEN Guoguang, WANG Guangzhi,et al.The Study of Kalman Ballistic Filter for Adapting the Outliers[J].Journal of Ordnance Equipment Engineering,2017(12):93-97.

TJ765

A

2096-2304(2017)12-0093-05

(责任编辑唐定国)