再入弹道目标跟踪的球面单纯形-径向容积卡尔曼滤波算法

李春月，廖育荣，倪淑燕，陈 帅

0 引 言

弹道导弹具有射程远、速度快、精度高、突防能力强、杀伤威力大、效费比高等优点，已成为现代战争非常重要的进攻性武器。为了遏制日益加剧的弹道导弹威胁，世界各国都在大力发展导弹防御系统。雷达是导弹防御系统中的核心探测器，对弹道导弹的实时跟踪精度直接影响到导弹防御系统对导弹拦截的成功率。因此，高精度弹道导弹实时跟踪算法是目前研究的热点问题。

弹道目标实时跟踪问题本质上为高维非线性系统的最优状态估计问题。卡尔曼滤波是目前最优状态估计中应用最为广泛的一种算法，经典的非线性卡尔曼滤波算法主要有扩展卡尔曼滤波（Extended KalmanFilter，EKF）算法[1]和无迹卡尔曼滤波（Unscented Kalman Filter，UKF）算法[2,3]。其中，EKF算法的核心思想是将非线性的系统状态函数做线性化处理，具体方法是利用泰勒级数展开法对非线性状态方程进行一阶展开，然后再用离散卡尔曼滤波进行最优状态估计，这种算法对非线性较强的系统估计精度不够高；UKF算法的出发点是基于“对概率分布进行近似要比对非线性函数近似容易很多”的思想，采用Sigma点的分布近似表示非线性函数的分布，有效提高了估计精度，但是对于高维非线性系统（维数n≥4），UKF算法中的自由调节参数κ＜ 0，使得中心采样点的权值κ＜ 0，从而使 UKF 算法在滤波过程中可能会出现协方差非正定情况，导致最终收敛结果不稳定甚至发散，并且随着系统维数的增加，采样点与中心点距离不断增大，导致非局部效应，严重影响 UKF 算法的估计精度。2009年，Arasaratnam等人[4,5]采用Spherical -Radial规则近似非线性函数传递的后验均值和协方差的方法，依据高斯滤波框架提出容积卡尔曼滤波（Cubature Kalman Filter，CKF）算法，相比于UKF算法，CKF算法有严格的数学推导过程，降低了运算量，提升了估计精度，得到广泛应用[6～9]。文献[10]对UKF算法和CKF算法的估计精度和数值稳定性做了比较，指出当系统状态维数小于3维时，UKF算法估计精度高于CKF算法，当系统状态维数等于3维时，二者估计精度相当，当系统状态维数大于3维时，无论估计精度和数值稳定性，CKF算法都高于 UKF算法；受到Spherical -Radial规则的启发，文献[11]提出了更多容积点的CKF算法，以提高估计精度，但是也相应增加了计算量；文献[12]为了进一步提高 CKF 算法的估计精度，提出迭代 CKF 算法，并将其应用于再入弹道目标状态估计，其效果优于 CKF 算法。

针对如何进一步提高对再入弹道目标的实时跟踪精度，本文提出将另一种基于Spherical Simplex-Radial准则的SSRCKF算法用于再入弹道目标实时跟踪中。该算法将Spherical Simplex准则引入CKF算法，以一种新的容积点选择方法来计算球面积分，进而近似计算高斯域下的贝叶斯积分，有效改善了CKF算法的估计精度，将其用于再入弹道目标跟踪中，提高弹道目标的实时跟踪精度，最后通过仿真分析，验证了算法的有效性。

1 弹道目标跟踪数学模型

1.1 再入段动力学模型

对再入弹道目标进行动力学建模，首先，对再入弹道目标进行受力分析。设弹道目标一般保持零度攻角再入，受到的空气动力表现为大气阻力，大气阻力加速度方向与再入速度方向相反。定义以弹道测量雷达为原点的北天东（North-Up-East）坐标系，O-xyz，Ox指向正北方向，Oy垂直于地球表面指向正上方，Oz指向正东方向，由于再入弹道导弹速度快，再入时间短，忽略地球自转角速度的影响。对再入目标建模如图 1所示，假设地球为标准球体，EO为地球质心，ag和af分别为引力加速度和大气阻力加速度，弹道目标的位置、速度分量和弹道系数构成七维状态向量，给出再入弹道目标的系统状态方程。

图1 雷达观测再入弹道目标示意Fig.1 Radar Observation Reentry Ballistic Target Signal

式 中 [ω1( t ) … ω7(t )]Τ为 系 统 状 态 噪 声 ， 满 足ω～0(,)Q的统计特性 ；µ为地球重力常数，µ= G M = 3 .986× 1 014m3/s2；为 地 球 半 径 ，Re= 6 378137m ； R(t)为目标到地心距离，()Vt为目标再入速度，ρ(t)为大气密度模型，ρ(t ) = ρ0exp ( −(R (t) − Re)H0)，H0为大气密度标高，H0=6700km，ρ0为海平面的大气密度，ρ0= 1 .225kg/m3；β( t)为弹道系数模型，有[13]：

式中DC为气动阻力系数；m为弹道目标质量；A为弹道目标有效截面积。为了保证弹道系数的非负性，防止滤波器发散采用式（2）的指数模型；()tγ的变化率可以用零均值高斯白噪声表示，即：

由此可以得到再入弹道目标的7维状态方程，状态量，为了保证滤波器的精度，采用四阶龙格-库塔方法对状态方程进行离散，最终得到离散非线性系统状态方程：

1.2 量测模型

以观测雷达为坐标系原点，在北天东坐标系下建立再入弹道目标量测模型。如图1所示，雷达对弹道目标的距离为 r，俯仰角为η，方位角为ε，有z = [r,ε,η]Τ。由观测值相对弹道目标位置的几何关系建立非线性量测方程：

式中 vk为雷达测量噪声向量，，满足v～0(,)R的统计特性，且kv，kω和kx互不相关。由此得到系统的量测方程：

2 球面单径容积卡尔曼滤波跟踪算法

2.1 Spherical Simplex-Radial准则

高斯域下的贝叶斯滤波可以一般化为下面的积分公式：

的求解一般可以采用一系列函数通过点集加权求和的近似方法，即：

式中 N为总点数； xi， wj分别为正交点集和权重。对式（7）采用Spherical Radial变换： x =ry且 yΤy = 1；xΤx =r2，则式（7）可重写为

式中 σ（·）为球面区域= { y ∈ RnyΤy = 1}的面积微元。通过式（9）可以看出，高斯域下的贝叶斯滤波公式通过Spherical-Radial准则变换，被分解为一个球面积分和一个径向积分。这两个积分同样无法直接求得，此时再利用式（8），球面积分用包含sN个点的球面积分准则计算[4]：

径向积分可由包含rN个点的高斯求积准则求得：

下面以这两个积分为基础介绍 Spherical Simplex准则和Radial准则。

2.1.1 Spherical Simplex准则

首先计算球面积分，球面积分 ()Sr的一种高效率的计算方法是：选取一系列包含 n个正则单行顶点的点集 aj=[aj,1, aj,2, … ,aj,n]Τ, j = 1 ,2,… ,n +1作为容积点[14,15]（n为状态维数），采用中心对称的容积准则来近似球面面积。其中容积点的每个元素,jia 选取规则为

利用ja构造Spherical Simplex准则形式如下：

式中为单位球的表面积，，且

计算径向积分S(r)rn−1e−r2dr ，与式（13）相对应，通过与式（11）进行匹配可以得到 Nr=1的Radial准则[12]：

2.1.2 Radial准则

由 Γ ( n + 1 ) = n Γ ( n)对第 2个等式进行化简，解出，权重 wr,1=Γ( n/2)/2。

非线性函数与高斯概率密度的乘积的积分是高斯域 下 贝 叶 斯 滤 波 器 的 核 心 ， 当(x ) = e−xΤx，w2(x ) = N (x ; 0, I)时，通过式（7）可得：

由式（13）、式（14）和式（20）可以求得：

把式（13）、式（16）代入式（18），即可得到 Spherical Simplex-Radial准则，用 ()Qf表示：

权重为 N (x;µ,P)时，µ，P分别为x的均值与协方差，对式（19）进行线性变换，得到一般形式Spherical Simplex-Radial 准则：

2.2 SSRCKF算法

基于Spherical Simplex-Radial 准则的容积卡尔曼滤波采用一系列等权值的容积点对高斯域下的贝叶斯滤波进行非线性逼近。其中容积点的选取规则为

式中 2m n= ，（n为状态空间维数）；为权值；,j = 1,2,… ,n +1，由式（12）得到；的第i列。

算法具体实现步骤如下：

a）步骤1：滤波器初始化。

SSRCKF滤波器初始化如下：

循环，完成以下步骤。

b）步骤2：时间更新。

首先对进行 Cholesky分解，有公式，取下三角阵，计算容积点：

状态方程传递容积点：

估计k时刻状态预测值：

估计k时刻的先验估计误差协方差：

c）步骤3：量测更新。

对先验估计误差协方差k−P进行Cholesky分解：

计算容积点：

量测方程传递容积点：

估计k时刻量测预测值：

估计k时刻的量测误差协方差：

估计k时刻的一步预测互相关误差协方差：

计算k时刻滤波增益：

计算k时刻的状态估计值：

计算k时刻的后验估计误差协方差：

通过完成上述3个步骤，完成了由k-1时刻到k时刻的后验估计值与后验估计误差协方差的更新。

3 仿真实验与分析

在Matlab（R2010b）环境下对雷达观测载入弹道目标系统进行建模仿真，对提出的基于 Spherical Simplex-Radial 准则的容积卡尔曼滤波算法进行验证。首先基于状态模型，在雷达站坐标系下生成带有状态噪声的载入弹道目标的真实轨迹，用四阶龙格-库塔方法对状态方程进行离散时，离散步长h=1；然后根据量测模型生成带有量测噪声的测量信息用于滤波。再入目标初始状态以及弹道系数常值为

雷达每秒对再入弹道目标进行一次测量，状态噪声协方差Q与量测噪声协方差R分别为

通过仿真，得到再入弹道目标的运动轨迹，如图2所示。

图2 再入弹道目标的运动轨迹Fig.2 The Trajectory of the Ballistic Trajectory

给定滤波初值如下：

分别以UKF，CKF和SSRCKF算法对再入弹道目标进行实时跟踪，设蒙特卡洛打靶次数为500。根据速度和位置均方根误差对比 3种算法的性能。位置均方根误差定义为

式中 N为蒙特卡洛打靶次数；，分别为第 n次打靶时 k时刻再入弹道目标的真实值与滤波值，速度均方根误差与式（37）方法相同。

仿真结果如图3、图4所示。图3为弹道目标速度均方根误差，图4为弹道目标位置均方根误差。从图3、图4中，可以看出，SSRCKF算法相比于UKF算法与CKF算法有更高的再入弹道目标实时跟踪精度。

图3 弹道目标速度均方根误差Fig.3 Velocity Root-mean-square Error of the Ballistic Target

图4 弹道目标位置均方根误差Fig.4 Root Mean Square Error of Trajectory Target

统计 200～250 s，UKF、CKF、SSRCKF 3 种算法对再入弹道目标实时跟踪的速度与位置的均方根误差，求出其平均值，如表1所示。由表1可以看出：相同条件下对再入弹道目标进行实时跟踪，SSRCKF算法在定位精度上比CKF算法提高了约4.5 m，比UKF算法提高了5 m；在定速精度上比CKF算法提高了约0.6 m/s，比UKF算法提高了0.7 m/s，充分证明了算法的有效性。

表1 200～250 s速度与位置均方根误差平均值Tab. 1 Mean Square Root Mean Square Error of 200～250 s

4 结 论

本文首先对再入弹道目标进行动力学建模，然后以推导的基于Spherical Simplex-Radial准则的容积卡尔曼滤波算法对再入弹道目标进行实时跟踪，该算法在没有明显提高计算量的前提下，有效提高了对再入弹道目标的实时跟踪精度；最后通过仿真验证，证明所提出的方法较UKF与经典的CKF算法，有更好的性能。