中国指挥与控制学会会刊 
RAM(Rocket, Artillery and Mortar)类武器是火箭弹、榴弹、迫击炮弹等武器的统称,其威力大,火力猛,机动性能好,射程远,火力覆盖能力强,操作方便,是当前战场主要的火力压制武器,所以,用于防御这类武器的C-RAM系统是新形势下部队防护的迫切需求。RAM类武器作为弹道目标,其实际运动规律与高炮火控系统在解相遇问题的过程中常用的匀速、 匀加速直线目标运动规律不符[1]。因此,有必要研究面向RAM类目标的跟踪与预测算法。
RAM类目标被动段会沿着一个由目标类型以及相应的性能特征事先确定好的飞行路径运动,具有显著的动力学特性,对该类目标建模的主要依据是目标的动力学因素,其中,最典型的参数就是质阻比[2]。对于合作目标,其质阻比已知,可直接根据目标运动方程进行跟踪,并对后续轨迹进行预测。
然而,在实际作战场景下,来袭的RAM类目标往往是非合作目标,质阻比未知,需要对其进行估计[3]。质阻比的估计精度和收敛速度与目标状态的估计相互耦合,当前估计质阻比主要有两种方法。1)利用滤波算法实时估计目标的距离、速度、加速度,根据质阻比物理意义与目标状态及雷达量测的关系推导出质阻比计算公式[4],该方法在建模过程中不考虑质阻比的影响,往往模型与运动特性不符,估计精度较差。2)基于目标的运动规律,将质阻比与选取的各状态分量一起组成扩维的状态向量,利用非线性滤波实时估计质阻比,该方法充分将目标运动规律与质阻比相结合,估计精度较高。文献[5]将质阻比建模为常数, 采用高斯分布对其进行描述,但因为对质阻比建模不满足其慢时变特性,会影响滤波器的收敛速度和精度。文献[6]采用指数函数在每一个小的时间间隔内对质阻比建模, 该方法计算量大且收敛速度不理想,达不到实际作战场景下对来袭目标快速跟踪的要求。
基于此,本文采用一阶马尔科夫过程对质阻比建模,通过CKF算法对质阻比进行估计,并加入强跟踪思想来应对环境干扰导致的模型失配问题。
1 RAM类目标动力学模型
RAM类目标为典型的弹道目标,主动段结束后,被动段一般保持零攻角[7]。对RAM类目标建模如图1 所示,不考虑目标的外形和姿态变化,目标在运动过程中受到重力和空气阻力作用。空气阻力与速度方向相反,重力指向地心。
忽略地球自转,在雷达东北天坐标系下,对RAM类目标被动段的运动状态进行描述,其状态模型可以表示为:
式中, 为零均值高斯白噪声, 为RAM类目标在雷达东北天坐标系下的位置, 为目标到地心的距离, 为地球半径, 。 为目标在雷达东北天坐标系下的速度分量, 为目标合速度,h为海拔高度,ρ(h)为海拔高度h处的空气密度。 为地球重力常数。β(t)为质阻比,是弹头质量和外形的组合参数[8],其数学表达式为 ,式中,m为炮弹质量,S为目标有效气动面积,CD(t)为空气阻力系数。
为了对质阻比进行实时估计,需要对前面建立的模型进行状态扩维作为滤波器的状态模型。常见方法是直接将质阻比作为一个常量来扩维,但这样不能保证质阻比的非负性。为了防止滤波发散,这里采用了式(2)的指数模型:
式中, 为质阻比初始值,可选取典型RAM类武器的质阻比参数。考虑质阻比变化为慢时变过程,q(t)应为一个慢时变参数,这里使用能够较好地模拟慢时变的一阶马尔科夫过程对q(t)建模,即
式中, 为一阶马尔科夫过程的时间常数,常取为一个较大值, 为零均值高斯白噪声。增广后的RAM类目标状态向量为
采用四阶龙格-库塔法对状态方程进行离散,最终得到目标的离散非线性系统状态方程为
2 雷达量测模型
3 强跟踪容积卡尔曼滤波(STCKF)
CKF适用于任何非线性模型,高维的情况下也能够有效跟踪。算法具体流程如下:
步骤1:初始化。
给定滤波器初始条件 ${{\hat{x}}_{0}}{{P}_{0}}{{Q}_{0}}{{R}_{0}}$。循环执行步骤2~步骤3。
步骤2:预测。
根据k-1时刻估计均值 ${{\hat{x}}_{k-1}}$和协方差矩阵,计算容积点:
表示n维空间的点集,即:
$\chi_{\text{k}-1}^{i}={{\hat{x}}_{k-1}}+{{S}_{\text{k}-1}}{{\xi }_{i}}$
状态方程传递容积点:
$\hat{\chi}_{\text{k}|\text{k}-1}^{i}=f(\chi_{\text{k}-1}^{i})$
计算状态预测均值和预测协方差:
${{\hat{x}}_{\text{k}|\text{k}-1}}=\frac{1}{m}\overset{m}{\mathop{\underset{i=1}{\mathop \sum }\,}}\,\hat{\chi}_{\text{k}|\text{k}-1}^{i}$
${{P}_{\text{k}|\text{k}-1}}=\frac{1}{m}\overset{m}{\mathop{\underset{i=1}{\mathop \sum }\,}}\,\left({\hat{\chi}}_{\text{k}|\text{k}-1}^{i}-{{{\hat{x}}}_{\text{k}|\text{k}-1}} \right){{(\hat{\chi}_{\text{k}|\text{k}-1}^{i}-{{\hat{x}}_{\text{k}|\text{k}-1}})}^{\text{T}}}+{{Q}_{\text{k}-1}}$
步骤3:更新。
根据预测均值 ${{\hat{x}}_{\text{k}|\text{k}-1}}$和协方差,计算容积点:
$\chi_{\text{k}|\text{k}-1}^{i}={{\hat{x}}_{\text{k}|\text{k}-1}}+{{S}_{\text{k}|\text{k}-1}}{{\xi }_{i}}$
量测方程传递容积点:
$\hat{Z}_{\text{k}|\text{k}-1}^{i}=h\left( \chi _{\text{k}|\text{k}-1}^{i} \right),i=0,\ldots,m=2n$
计算量测预测均值、新息协方差、状态与量测间的互协方差和滤波增益:
${{{\hat{z}}}_{\text{k}|\text{k}-1}}=\frac{1}{m}\overset{m}{\mathop{\underset{i=1}{\mathop \sum }\,}}\,\hat{Z}_{\text{k}|\text{k}-1}^{i}$
${{S}_{k}}=\frac{1}{m}\overset{m}{\mathop{\underset{i=1}{\mathop \sum }\,}}\,\left(\hat{Z}_{\text{k}|\text{k}-1}^{i}-{{{\hat{z}}}_{\text{k}|\text{k}-1}} \right){{(\hat{Z}_{\text{k}|\text{k}-1}^{i}-{{\hat{z}}_{\text{k}|\text{k}-1}})}^{\text{T}}}+{{R}_{k}}$
${{C}_{k}}=\frac{1}{m}\overset{m}{\mathop{\underset{i=1}{\mathop \sum }\,}}\,\left(\hat{\chi}_{\text{k}|\text{k}-1}^{i}-{{{\hat{x}}}_{\text{k}|\text{k}-1}} \right){{(\hat{Z}_{\text{k}|\text{k}-1}^{i}-{{\hat{z}}_{\text{k}|\text{k}-1}})}^{\text{T}}} $
计算k时刻后验状态估计均值和协方差:
${{\hat{x}}_{k}}={{\hat{x}}_{\text{k}|\text{k}-1}}+{{K}_{k}}\left( {{z}_{k}}-{{{\hat{z}}}_{\text{k}|\text{k}-1}} \right)$
CKF算法实现简单且估计精度高,但当目标受到外界干扰模型失配时,估计精度将会变差,甚至发散。强跟踪滤波器更加强调量测信息的作用,可以更好地应对目标的状态突变,其原理是选择一个适当的时变增益矩阵 ,使得状态估计残差方差最小,强迫输出残差序列保持相互正交,即:
$E[\left( {{x}_{k}}-{{{\hat{x}}}_{k}} \right)( {{x}_{k}}-{{{\hat{x}}}_{k}}{{)}^{\text{T}}} ]=\text{min}$
式中,ε为残差。
在式(19)的新息协方差中引入时变渐消因子 ,即
${{S}_{k}}={{\lambda }_{k}}\left[ \frac{1}{m}\overset{m}{\mathop{\underset{i=1}{\mathop \sum }\,}}\,(\hat{Z}_{\text{k}|\text{k}-1}^{i}-{{{\hat{z}}}_{\text{k}|\text{k}-1}} ){{(\hat{Z}_{\text{k}|\text{k}-1}^{i}-{{{\hat{z}}}_{\text{k}|\text{k}-1}})}^{\text{T}}} \right]+{{R}_{k}}$
的计算公式为
${{\lambda }_{k}}=\text{max}\left\{ 1,\frac{\text{tr}( {{C}_{0,k}}-{{R}_{k}})}{\text{tr}(\frac{1}{m}\overset{m}{\mathop{\mathop{\sum }_{i=1}}}\,\left(\hat{Z}_{\text{k}|\text{k}-1}^{i}-{{{\hat{z}}}_{\text{k}|\text{k}-1}} \right)\left(\hat{Z}_{\text{k}|\text{k}-1}^{i}-{{{\hat{z}}}_{\text{k}|\text{k}-1}}{{)}^{\text{T}}} \right)} \right\}$
其中, ,0<ρ<1为遗忘因子,通常取0.95^[13]。
4 仿真验证
为验证本文方法,设计地面雷达对RAM类目标跟踪的仿真场景。以某口径的火箭弹为例,根据炮弹三自由度运动方程,在雷达站坐标系下生成目标的真实轨迹,发射高低角为45°,落点在雷达附近。雷达扫描周期为20 ms,测距精度为 10 m,测角精度为 1 mrad(包括高低角和方位角)。节选目标飞行中段及末段25 s轨迹作为跟踪对象,并加入10 m/s的横风模拟外界环境的干扰,使轨迹产生横偏,运动轨迹如图3 所示。
为验证所建立模型的优势,采用三种方法对比仿真,分别对有风和无风状态下的目标进行跟踪。
1)把质阻比建模为常数的CKF方法,状态噪声协方差为
2)基于改进模型的CKF方法,马尔科夫模型的时间常数设为 ,状态噪声协方差为
3)基于改进模型的STCKF方法,参数设置与方法2一致。
本文共进行100次蒙特卡洛仿真,质阻比初值均为β0=4 000 kg/m2。 无风条件下的跟踪结果如图4 、5 、6 和表1 所示。从图4 可看出,利用传统模型估计质阻比可以收敛至真值附近,但收敛速度缓慢。相比之下,本文模型收敛速度大幅提高,比普通方法快了40%,且精度较好。因为基于一阶马尔科夫过程的质阻比模型可以更好地表征RAM类目标质阻比的慢时变过程,从而加快了滤波的收敛速度。从图6 和图7 可以看出,目标状态估计的收敛速度和精度也有了明显提高,精度提高了30%。但STCKF方法和CKF方法在精度上并没有太大的区别,因为无风条件下本文建立的模型与目标较匹配,引入强跟踪算法并不会对跟踪效果有明显的提升。
有风条件下的跟踪结果如图7 、8 、9 和表2 所示,从图7 可以明显看出,因为存在横风的扰动,模型失配,基于改进模型的CKF方法对质阻比的估计效果会明显变差。对比之下,引入了强跟踪的STCKF算法对质阻比的估计依然能够收敛到真值,位置和速度跟踪误差对比普通CKF方法也有一定改善,因为强跟踪算法更加强调量测信息的作用,对系统模型不准确具有一定的鲁棒性,增加了估计的稳定性。
5 结束语
本文以RAM类目标为对象,根据其运动特性,提出了一种状态与质阻比联合估计的目标跟踪方法。考虑传统建模中将质阻比视为常量,与质阻比的慢时变过程不符,本文采用一阶马尔科夫过程对质阻比建模。为提高对弹道类目标跟踪的可靠性,引入强跟踪容积卡尔曼滤波器。经仿真验证,该模型相比于常用的常数模型,跟踪的收敛速度提高了40%,精度提高了30%, 强跟踪滤波器的引入也提高了算法的鲁棒性。