中国指挥与控制学会会刊 
天基红外导弹预警系统探测弹道导弹,主要是通过探测发动机尾焰辐射信号的角度信息来估计目标运动状态,其典型代表是美国的天基红外系统(Space Based Infrared System,SBIRS)[1]。由于探测器存在系统误差,且长距离探测受到大气层较强的干扰,探测数据中同时含有随机误差。为了实时估计目标的运动状态,需要对导弹进行运动建模,并设计滤波算法对导弹运动进行滤波。助推段弹道导弹机动性强,当前统计(Current Statistic, CS)模型[2]将导弹机动加速度看作一阶时间相关过程,能较好地匹配导弹在助推段的运动情况。目前,常用的非线性滤波算法主要有扩展卡尔曼滤波(Extended Kalman Filtering,EKF)[3]、无迹卡尔曼滤波(Uncented Kalman Filtering,UKF)[4]和容积卡尔曼滤波(Cubature Kalman Filtering,CKF)[5]等。SBIRS探测模型的非线性较强,采用EKF进行估计可能会产生较大的截断误差[6]。导弹状态矢量维数较高,采用UKF进行估计可能出现状态协方差矩阵不正定的情况,导致滤波中断[7]。CKF算法是利用容积点非线性传播的后验统计特性来近似高斯积分,可以克服上述问题,但时间更新和测量更新两个阶段均需构造2n个容积点(n为机动目标状态矢量的维数),当n较大时运算量较大。由于导弹在助推段运动时间短,弹道估计对滤波算法时效性要求比较高,有必要从运算效率上对CKF算法进行改进。
本文提出一种简化CS-CKF的助推段弹道估计方法,测量更新采用构造容积点非线性传播的方式实现,而状态和协方差一步预测按照线性状态方程计算。
下面依次介绍基于CS-CKF的弹道估计方法和基于简化CS-CKF的弹道估计方法,通过仿真实验来验证简化CS-CKF弹道估计方法的优越性。
1 简化CS-CKF的弹道估计方法
1.1 CS模型
用CS模型描述导弹在助推段的运动状态,离散形式的状态方程为
其中,Xk、Φk、Gk、wk、 分别为k时刻导弹状态矢量、状态转移矩阵、输入控制矩阵、过程噪声矢量、导弹机动加速度的均值。
当探测时间间隔较小时,一个时间间隔内导弹的加速度可认为是常数,即输入控制矩阵Gk为0。因此,CS模型可以简化为分段常加速度模型为
Xk+1=ΦkXk+wk
其中,Φk= ,wk为过程噪声矢量,I3×3为三阶单位矩阵,03×3为三阶零矩阵,T为探测时间间隔。
过程噪声矩阵
其中,α为机动频率, 为机动加速度的方差,
机动频率α为机动时间常数的倒数,假设导弹在整个助推段都在机动,则机动时间常数就等于助推段飞行时间。CS模型采用修正的瑞利分布来描述机动加速度均值 ,机动加速度方差 的表达式为
式中amax为机动加速度的最大值。
1.2 滤波算法初始化
按照文献[8]提出的方法,求出k=0~2时刻导弹等效探测位置r0、r1、r2,使用r0、r1、r2对滤波器进行初始化。具体方法如下:
假设导弹前三个时刻做常加速运动,X2为k=2时刻的九维状态矢量,则k=0~2时刻的探测方程可写成
式中,νS为k=0~2时刻等效位置探测噪声组成的矢量,相应的协方差矩阵RS=E(νS )。令权重矩阵W= ,利用加权最小二乘算法计算k=2时刻X2的估计值 =[ WhS]-1 WYS,则相应的状态协方差矩阵P2/2=[ WhS]-1 WRSWhS[ WhS]-1,式中的下标2均表示k=2时刻。
1.3 CS-CKF算法
在获得k时刻的状态估计值 以及相应的协方差矩阵Pk后,一个完整的CS-CKF运算周期分为两个阶段,共8个步骤。
时间更新阶段:
1)构造容积点
对Pk矩阵作Cholesky分解,即Pk=Sk ,构造容积点
$\boldsymbol{X}_{k}(i)=\boldsymbol{S}_{k} \xi_{i}+\hat{\boldsymbol{X}}_{k}, i=1,2, \cdots 2 n $
其中,n为导弹运动状态矢量 的维数,ξi= [I3×3 -I3×3]i,[I3×3 -I3×3]i代表[I3×3 -I3×3]矩阵的第i列,I3×3为3阶单位矩阵。
2)容积点基于状态方程的非线性传播
χk+1/k (i)=f(χk (i))+wk
3)计算一步状态预测均值 和一步状态预测误差协方差矩阵Pk+1/k
$\left\{\begin{array}{l} \hat{\boldsymbol{X}}_{k+1 / k}=\frac{1}{m} \sum_{i=1}^{m} \chi_{k+1 / k}(i) \\ \boldsymbol{P}_{k+1 / k}=\frac{1}{m} \sum_{i=0}^{m}\left[\chi_{k+1 / k}(i)-\hat{\boldsymbol{X}}_{k+1 / k}\right] \\ {\left[\chi_{k+1 / k}(i)-\hat{\boldsymbol{X}}_{k+1 / k}\right]^{\mathrm{T}}+\boldsymbol{Q}_{k}} \end{array}\right.$
其中,Qk为CS模型构造的过程噪声矩阵。
测量更新阶段:
4)基于一步状态预测重新构造容积点
Pk+1/k=Sk+1/k ,则重新构造的容积点
$\chi_{k+1 / k}(i)=S_{k+1 / k} \xi_{i}+\hat{X}_{k+1 / k},(i=0,1, \cdots, 2 n) $
5)容积点基于探测方程的非线性传播
$\hat{\boldsymbol{Y}}_{k+1 / k}(i)=\boldsymbol{h}\left(\chi_{k+1 / k}(i)\right), i=1,2, \cdots 2 n $
6)计算探测量一步预测均值
$\hat{\boldsymbol{Y}}_{k+1 / k}=\frac{1}{m} \sum_{i=1}^{m} \hat{\boldsymbol{Y}}_{k+1 / k}(i)$
7)计算滤波增益
首先,计算k+1时刻探测误差自协方差矩阵 和一步预测互相关协方差矩阵 。
$\left\{\begin{array}{l} \boldsymbol{P}_{k+1}^{Y Y}=\frac{1}{m} \sum_{i=1}^{m}\left[\boldsymbol{Y}_{k+1}(i)-\hat{\boldsymbol{Y}}_{k+1 / k}\right]\left[\boldsymbol{Y}_{k+1}(i)-\hat{\boldsymbol{Y}}_{k+1 / k}\right]^{\mathrm{T}}+\boldsymbol{R}_{k} \\ \boldsymbol{P}_{k+1 / k}^{X Y}=\frac{1}{m} \sum_{i=1}^{m}\left[\boldsymbol{X}_{k+1}(i)-\hat{\boldsymbol{X}}_{k+1 / k}\right]\left[\boldsymbol{Y}_{k+1}(i)-\hat{\boldsymbol{Y}}_{k+1 / k}\right]^{\mathrm{T}} \end{array}\right. $
其中Rk为探测噪声矩阵。而后,计算k+1时刻滤波器的滤波增益Kk+1= ( )-1。
8)计算k+1时刻状态估计值和状态误差协方差矩阵
$\left\{\begin{array}{l} \hat{\boldsymbol{X}}_{k+1}=\hat{\boldsymbol{X}}_{k+1 / k}+\boldsymbol{K}_{k+1}\left(\boldsymbol{Y}_{k+1}-\hat{\boldsymbol{Y}}_{k+1}\right) \\ \boldsymbol{P}_{k+1}=\boldsymbol{P}_{k+1 / k}-\boldsymbol{K}_{k+1} \boldsymbol{P}_{k+1}^{\mathrm{T}} \boldsymbol{K}_{k+1}^{\mathrm{T}} \end{array}\right. $
1.4 简化CS-CKF算法
在获得k时刻的状态估计值 以及相应的协方差矩阵Pk/k后,简化CS-CKF算法的递推处理流程如图1 所示。
算法的一个周期分为两个阶段。
时间更新阶段:
1)状态一步预测
$\hat{\boldsymbol{X}}_{k+1 / k}=\boldsymbol{\Phi}_{k} \hat{\boldsymbol{X}}_{k / k} $
2)状态协方差一步预测
$\boldsymbol{P}_{k+1 / k}=\boldsymbol{\Phi}_{k} \boldsymbol{P}_{k} \boldsymbol{\Phi}_{k}^{\mathrm{T}}+\boldsymbol{Q}_{k} $
测量更新阶段与CS-CKF算法相同,不再赘述。
从以上步骤可以看出,单个周期内,简化CS-CKF算法的时间更新阶段仅需计算1次,而CS-CKF算法时间更新阶段需计算2n次。
2 仿真实验
2.1 仿真流程
仿真流程如图2 所示,主要分为6个步骤。
1)仿真一条助推段弹道,设置双星位置、探测时间间隔和视线测量误差,利用双星探测模型构造一组含噪声的探测数据。
2)利用双星位置,将双星探测数据转化为导弹等效探测位置。
3)利用前3个时刻等效探测位置,对滤波算法进行初始化。
4)设置CS模型相关参数,构造过程噪声协方差矩阵。
5)利用简化CS-CKF算法进行滤波递推处理,计算导弹各个时刻的估计值,进而得到估计的弹道数据。
6)将估计的弹道数据与仿真的弹道数据进行对比分析,验证算法的有效性。
2.2 仿真场景
以某三级弹道导弹为例,假设导弹发射点的地理坐标为160°W、30°N,海拔为40 m,发射方位角为302.28°,导弹最大负攻角为6.57°。通过计算,得到导弹助推段三维弹道曲线,如图3 所示。
分别定点于162°W和110°E的两颗GEO发现导弹目标,首次探测到导弹的时间为20 s,探测间隔为1 s,最后探测到导弹的时间为196 s,视线测量误差为50 μ rad。
2.3 参数设置
CS模型中,导弹机动频率为0.005 Hz,在地球固联坐标系中描述,导弹x轴、y轴、z轴方向加速度最大值均为95 m/s2,卫星探测时间间隔为1 s。
2.4 仿真结果及分析
1)估计误差
对CS-CKF和简化CS-CKF算法蒙特卡洛仿真1000次,以x方向为例,状态估计误差曲线如图4 所示。用实线描述CS-CKF算法的估计均方根误差,用点虚线描述简化CS-CKF算法的估计均方根误差。从图4 可以看出,两种算法在x方向状态估计均方根误差差别非常小。y方向、z方向与x方向情况相同,说明两种算法的状态估计误差相当。
2)一步预测标准差
对CS-CKF和简化CS-CKF算法进行单次仿真,对每个时刻一步状态协方差预测矩阵Pk+1/k的元素Pk+1/k(1,1)、Pk+1/k(4,4)、Pk+1/k(7,7)进行开平方,分别得到一步预测x方向位置、速度、加速度的标准差。对各个时刻两种算法的一步预测标准差进行相减,得到两种算法x方向一步预测标准差的差值情况,如图5 所示,与状态估计值相比,两种算法的一步预测标准差差值非常小,可以忽略。其他方向一步预测标准差的差别情况与x方向相同,由此可知,两种算法一步预测的一致性比较强,本文算法是合理的。
3)运算时间
实验平台计算机CPU为Core(TM) i7-4790(3.60 GHz),内存为8 GB,显卡为AMD Radeon R7350;软件编程环境为Matlab2017b。对CS-CKF算法和简化CS-CKF算法进行1000次蒙特卡洛仿真,简化CS-CKF滤波函数运算时间为34.76 s,CS-CKF滤波函数运算时间为48.09 s,简化CS-CKF滤波函数运算时间比CS-CKF滤波函数运算时间缩短38%。
进一步对两种滤波函数的代码进行比较,耗费时间相差较大的步骤如表1 所示。由表1 可以看出,两种滤波函数运算时间的主要差别在时间更新阶段。CS-CKF滤波函数时间更新阶段耗费的总时间为13.19 s,其中构造容积点、基于容积点进行状态一步预测和协方差一步预测耗费的时间分别3.06 s、3.97 s和6.16 s。简化CS-CKF滤波算法的时间更新是直接利用线性状态方程实现的,因此运算时间较短。
表1 两种滤波函数耗费时间差别较大的步骤运算时间表 |
| 耗费时间差别 较大的步骤 | CS-CKF 滤波函数 | 简化CS-CKF 滤波函数 |
|---|---|---|
| 构造运动模型的 容积点 | 3.06 s | 无 |
| 基于容积点求 状态一步预测 | 3.97 s | 无 |
| 基于容积点求 协方差一步预测 | 6.16 s | 无 |
3 结束语
针对CS-CKF的弹道估计运算量较大的问题,提出一种简化CS-CKF的助推段弹道估计方法。依次介绍运动模型、估计方法、滤波初始化方法和递推处理步骤等,并进行仿真实验。实验结果表明,与传统CS-CKF的弹道估计算法相比,本文算法合理有效,在估计误差相当的条件下,运算时间大约缩短38%。对于实时性要求很高的助推段弹道估计来说,简化CS-CKF的算法优势比较明显。