在线弹道参数滤波与辨识方法分析

李小元,王中原,常思江

(1.海军研究院,北京 100161;2.南京理工大学 能源与动力工程学院,江苏 南京 210094)

炮弹在发射飞行过程中,存在火炮初始扰动、初速或然误差、炮弹结构参数误差、气象条件误差等,造成射弹散布。

弹道修正技术是指对发射出去的炮弹的一段飞行弹道参数进行实时测量,同预定弹道参数进行比较,进行弹道偏差计算,形成控制指令,通过弹上控制机构作用调节弹道向预定弹道逼近,从而大幅度减小炮弹散布的一种技术。如果仅对纵向弹道进行调节,称为一维弹道修正弹。当弹上控制机构采用阻力环装置来调节阻力,则为打远修近、阻力环体制的一维弹道修正弹,它是目前广泛采用、结构简洁的一种一维弹道修正技术[1-3]。

一维弹道修正技术在实际中要取得好的应用效果,一个重要问题是如何根据实测的一段飞行弹道参数快速准确地辨识出该弹的弹道特征参数,进而准确地预报后续弹道、落点,并准确地知道这发弹同预定弹道的偏差,从而实施弹道修正。目前,采用卡尔曼滤波、特征参数辨识方法是一种较为有效的方法。文献[4-5]分别采用卡尔曼滤波、速度-时间序列匹配等方法对该问题进行了研究。文献[6]通过建立弹道滤波方程,采用卡尔曼滤波法外推弹道,对弹道修正弹的落点进行了预测。为了进一步提高弹箭密集度,文献[7]比较了线性卡尔曼模型和扩展卡尔曼模型在GPS量测数据下的差异,并分析了2种模型在不同条件下的优缺点。但是上述文献在数据测量精度、基础阻力系数变化等因素对预报弹道落点的影响方面研究较少。

在卡尔曼滤波、特征参数辨识方法中,在对应整个可提供的实时测量弹道段上,分为预报和校正2个环节交替实施,核心目的是快速、准确辨识出广义气动力系数。在预报环节,弹的阻力系数状况对预报结果极为敏感,影响很大;在校正环节,雷达实时测量的弹道参数的精度对校正结果影响极大。它们对最终辨识出的广义气动力系数影响很大,直接影响在线弹道预报精度。本文对阻力环修正体制下的这些问题进行分析、讨论,给出相关建议。

1 一维弹道修正弹的飞行弹道控制方程

针对本文研究的一维弹道修正弹,根据其工作原理,可建立一维修正弹道方程:

(1)

控制方程为

(2)

式中:Fx2,Fy2,Fz2为作用在弹丸上的合力F在弹道坐标系x2y2z2中的投影,在计算Fx2,Fy2时需要使用阻力系数符合系数kcx和升力系数符合系数kcy,上述2个系数为待辨识的参数;Mξ,Mη,Mζ为作用在弹丸质心处的合力矩M在弹轴坐标系ξηζ中的投影;MD为阻力环附加的极阻尼力矩;FD为阻力环张开后增加的阻力;ΔCx为阻力系数增量;v为速度;θa为弹道倾角;ψ2为弹道偏角;ωξ,ωη,ωζ为弹丸角速度在弹轴坐标系ξηζ中的分量;φa为弹轴高低摆动角;φ2为弹轴侧向摆动角;γ为滚转角;x,y,z为弹丸在地面坐标内的位置分量;JA,JC分别为极转动惯量和赤道转动惯量;m为弹丸质量;ρ为大气密度;S为弹丸特征面积。各坐标系的定义参见文献[8]。

2 卡尔曼滤波与特征参数辨识的弹道预报模型和方法

2.1 卡尔曼滤波方法

如前所述,如何根据实测的一段飞行弹道参数快速、准确地辨识出该弹的弹道特征参数是一个核心问题。本文选取阻力系数符合系数kcx和升力系数符合系数kcy作为弹道特征参数,分别用于纵向射程与横向距离的弹道、落点预报。为了能够利用一段实测弹道数据尽可能准确地估计出这2个弹道特征参数,本文采用卡尔曼滤波方法,将kcx和kcy作为状态变量进行估计。卡尔曼滤波方法的优点是考虑了物理系统本身的动态模型和测量系统的量测误差,利用从某时刻k-1以前直到现在的测量数据Y(j个)逐次递推,可估计出k时刻最优的系统状态X。卡尔曼滤波基本方程简介如下。

预测方程:

(3)

(4)

(5)

滤波方程:

(6)

卡尔曼滤波增益:

(7)

协方差矩阵:

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

(8)

式中:Δt为卡尔曼滤波的采样间隔;Φk为基本矩阵,H为量测矩阵,Qk为过程噪声方差矩阵,Rk为量测噪声方差矩阵,其他各符号定义参见文献[9]。

下面以上述卡尔曼滤波方程为基础,结合具体的弹道模型,构建一维修正弹的弹道滤波模型。

2.2 弹道滤波模型

由于在利用卡尔曼滤波辨识弹道特征参数的过程中,不涉及阻力环机构作用。因此,在构建一维修正弹弹道滤波模型时,令方程(1)中FD=0,并将参数kcx和kcy分别置于阻力系数和升力系数之前作为乘子,并补充如下微分方程:

(9)

则状态变量取为

X=(vθaψ2xyzωξωηωζφaφ2γkcxkcy)T

(10)

将方程(1)改写成状态空间形式为

(11)

式中:fi(i=1,2,…,14)为方程(1)和方程(9)对应的右端函数。

对方程(11)进行线性化,可得:

(12)

式中:F为动力矩阵,其元素为fi(i=1,2,…,14)对每个状态变量的偏导数,由于具体形式冗长,为节约篇幅,这里略去。

在求出了动力矩阵F的具体表达式后,可求得弹道滤波系统的基本矩阵Φk为

Φk≈I14×14+F14×14·Δt

(13)

式中:I14×14为单位矩阵。

假设可直接或间接测得的弹道参数为

Z=(vθaψ2xyz)T

(14)

则弹道滤波系统的量测方程可以表达为

Z=h(X)+n=H·X+n

(15)

式中:n为量测噪声。

由于状态变量X为14×1维矩阵,Z为6×1维矩阵,故量测矩阵H为6×14维矩阵,n为6×1维矩阵。

(16)

量测噪声n的方差矩阵为

(17)

式中:σ为测量值的标准差。

另外,与状态方程相对应的还有过程噪声方差矩阵Qk,由于很难描述弹道模型与实际弹丸飞行状态的差异,故一般取Qk=O14×14。

数值迭代求解如下形式的Ricatti方程:

(18)

式中:Pk/k-1为k时刻数据更新前的状态估计值协方差矩阵;Pk为数据更新后的状态估计值协方差矩阵,相应地,Pk-1为上一时刻的状态估计值协方差矩阵;Kk为卡尔曼增益矩阵。

为了使卡尔曼滤波系统正常工作,状态估计值协方差矩阵的初始值P0根据量测噪声的标准偏差来选取,有:

P0=diag(P0,11,P0,22,…,P0,1414)

(19)

卡尔曼滤波的基本方程为

(20)

2.3 弹道预报方法

根据卡尔曼滤波方程(20)可以实时估计出弹道方程组对应的状态变量,将测量数据段最后一点的状态变量最优估计值作为弹道预报的起点,对弹道方程组进行数值积分直至落点,可预报出阻力环机构全弹道不作用条件下的该发弹的无控弹道及落点。将事先获取的阻力环阻力系数ΔCx代入弹道方程组中(即FD≠0),并预设一个阻力环启动作用时刻,可预报出该发弹在某一时刻张开阻力环后的修正弹道及落点。如果将无控弹道预测落点与目标弹的射程相减,取得该发弹的射程修正量,利用该修正量进行迭代计算,可很容易求出该射程修正量对应的阻力环机构启动时刻。

3 在线弹道预报分析

为了验证上述模型与方法的正确性和有效性,本节以某大口径一维弹道修正弹的实际数据为例,开展在线弹道预报分析。在该一维弹道修正弹的某次试验中,连续射击6发无控弹,采用两部弹道跟踪雷达同时对每发弹丸实时跟踪,分别测得对应的飞行弹道数据。这两部雷达的跟踪精度有差异,其中一部为高精度弹道跟踪雷达(记为HA雷达),而另一部的跟踪精度相对较低(记为LA雷达)。采用不同精度雷达同时跟踪的目的是为了研究测量数据精度对弹道预报的影响。需要说明的是,由于一维弹道修正弹的弹道预报主要关注纵向距离x,故下面分析均对应预报纵向距离x。

具体分析包括几个方面:研究不同阻力系数对弹道预报精度的影响;分析测量飞行弹道数据对弹道预报的影响。

3.1 在线弹道预报结果

利用HA雷达实测5 s数据,包括弹丸速度和位置,进行在线弹道预报,所得结果如图1～图3和表1所示。

其中,图1～图3分别为其中1发弹丸的位置、速度及阻力系数符合系数kcx滤波曲线,表1为6发弹丸预报射程与实测射程之差。从图1～图3可以看出,采用卡尔曼滤波方法得到了弹丸位置、速度及阻力系数符合系数的最优估计值。由表1所示的在线射程预报结果看,6发预报射程与实测射程之差出现了正负号,即对于这一组射弹,有的射程预报值大于射程实测值,有的则小于实测值。这表明由于各种随机误差因素的存在,在线弹道预报误差本身也存在随机性,而一维弹道修正弹的主要功能是提高射程密集度,因此应尽可能减小弹道预报误差的跳差范围,即尽量控制极差。从数值上看,试验弹的平均射程为45 700 m,而6发弹的射程预报极差为73 m(约为平均射程的0.16%),射程预报误差平均值为21 m(不到平均射程的0.05%),能够满足实际工程应用的要求。

图1 弹丸位置实测值与最优估计值

图2 弹丸速度实测值与最优估计值

图3 参数kcx的最优估计结果

表1 6发弹丸预报射程与实测射程之差

3.2 基础阻力系数变化对弹道预报的影响

在弹丸阻力系数标准值基础上变化±2.5%和±5%,将变化后的阻力系数代入前面建立的在线弹道预报模型中,所得预报射程与实测射程之差ΔXc列于表2中。

表2 不同基础阻力系数对应的预报射程与实测射程之差

表2中结果说明,弹丸阻力系数标准值变化对预报结果有一定影响。按照相对于标准值的变化百分比,预报射程误差呈线性变化。从一组数据来看,预报射程误差相对于标准值条件下的预报误差产生了一个系统偏差,6发弹丸预报射程误差的极差变化很小。从这个意义上讲,弹丸阻力系数标准值发生变化时,主要影响在线弹道预报结果的准确度,而不影响反映弹丸发与发之间差异的散布,即对一维修正落点的密集度无影响。从武器使用角度,应兼顾准确度和密集度,故有必要采用尽可能准确的阻力系数。由于本文通过大量试验取得了较为精确的弹丸阻力系数,因此预报准确度也相对较高。

3.3 测量数据对弹道预报的影响

试验中LA雷达斜距离标准差比HA雷达斜距离标准差大2 m,LA雷达高低角标准差比HA雷达高低角标准差大0.65 mrad,LA雷达方位角标准差比HA雷达方位角标准差大0.5 mrad。表3所示为采用HA和LA两部雷达跟踪数据预报出的射程值,其他条件参数均取为相同。Xm为实测射程,Xc,HA为HA雷达预报射程,Xc,LA为LA雷达预报射程,ΔXc,HA=Xc,HA-Xm,ΔXc,LA=Xc,LA-Xm。

表3 利用两部雷达跟踪数据预报的射程

从平均值上看,基于LA雷达跟踪数据的预报误差要小于基于HA雷达跟踪数据的预报误差,特别是对于2号弹、3号弹以及6号弹。但对于一维弹道修正弹,其原理主要是减小弹丸的纵向散布,故对于在线弹道预报而言,更重要的是尽可能减小一组射弹射程预报误差的极差。这一极差主要受跟踪数据精度的影响,由表3可知,HA雷达对应的射程预报误差极差为73 m,而LA雷达对应的射程预报误差极差为116 m。因此,应尽可能采用跟踪精度高的雷达。

测量数据的长度对预报结果也有影响。前面计算中,都是采用了5 s跟踪数据。下面针对HA雷达数据,以1号弹为例,研究不同数据长度对预报结果的影响,如图4所示。

图4 不同数据长度对应的射程预报误差

由图4可以看出跟踪数据长度对预报结果的影响。当数据时长大于5 s,预报误差略有波动,但基本上维持在25 m左右,数据时长越长对预报误差的影响较小,而数据越多,在线计算量也越大;当数据时长小于5 s,预报误差随着数据时长的减小而迅速增大,3 s时长对应的预报误差已接近70 m。由于图4所示是1号弹对应的结果,如果是对于6号弹,3 s时长对应的预报误差接近150 m,对修正弹最终密集度指标的实现极为不利。上述结果表明,当综合考虑在线弹道预报的精度和计算量(实时性),选取5 s数据是比较合适的。

3.4 关键参数对弹道预报的影响

前面所建基于卡尔曼滤波的在线弹道预报模型中,参数kcx也是作为一个状态变量。实际研究中发现,该状态变量对应的协方差矩阵元素初始值Pk,0对参数kcx的最优估计结果有影响,并且较为敏感,进而影响射程预报的结果。本节仍采用HA雷达跟踪数据,研究不同Pk,0值对预报结果的影响,结果如表4所示。表中,ΔXc为预报射程与实测射程之差。

表4 不同Pk,0值对应的预报射程与实测射程之差

图5为不同Pk,0值对应的射程预报误差极差(ΔXc,m)曲线图。

图5 不同Pk,0值对应的射程预报误差极差

从表4和图5中不难看出,尽管Pk,0的量值不大,但对预报结果却具有较大影响。前面研究表明,对于一维弹道修正弹,一组射程预报误差的极差大小可在一定程度上反映在线弹道预报的性能。当Pk,0>0.01,随着Pk,0值的增大,射程预报误差的平均值和极差都增大;当Pk,0<0.01,随着Pk,0值的减小,射程预报误差的平均值和极差也增大。图5所示极差曲线中,显然Pk,0=0.01是一个最优值。需要指出的是,对于不同弹丸、不同测量设备或不同模型,该最优值也不尽相同,可采用试验手段绘制类似图5的极差曲线,从而确定Pk,0的最优值。工程实践表明,这一方法是比较有效的。

4 结束语

本文针对阻力环体制的一维弹道修正弹,建立了飞行控制弹道模型,并基于卡尔曼滤波理论建立了弹道滤波模型,给出了弹道预报的方法,并通过算例对本文研究结果的正确性和有效性进行了分析和验证。通过对实际工程应用结果的分析, 提出以下建议:为了提高无控弹道和修正弹道的预报精度,应尽可能准确地获取全弹的阻力系数,系数值应覆盖全弹道马赫数。在确定某一实时测量弹道参数的基础上,进行若干发炮射试验,通过反复调整预报-校正分段、测量参数方差同实测弹道的对比,确定出较佳匹配(协方差、分段等)关系,来提高“广义阻力系数”的辨识效果,提高后续弹道预报精度。