基于刚体弹道模型的防空火控解算方法

何山，吴盘龙，李星秀，恽鹏，薄煜明

(1.南京理工大学 自动化学院, 江苏 南京 210094; 2.南京理工大学 理学院, 江苏 南京 210094)

0 引言

随着网络化、信息化技术的快速发展，防空反导作为现代战争的重要任务之一，发挥着越来越重要地位[1]。高炮以其灵活性高、射速快、机动性强和火力密集等优点成为了当前低空防御的主力[2]。

火控系统作为武器系统的神经中枢，其解算出的射击诸元精度直接影响武器系统的命中概率[3]。一般来说，火控系统解算的难点在于目标航迹和弹丸轨迹的准确预测[4]。而随着末端突防能力的日渐提升，目标强机动性使得在航迹的建模和状态估计上变得尤为困难；另一方面，复杂的战场环境和弹丸的一致性直接影响弹丸零升阻力、诱导阻力、升力以及马格努斯力等，从而给弹丸的弹道模型建立带来了极大的挑战。

现行防空火控系统中的弹道模型多以射表函数逼近为主，极大缩短了解算时间[5-6]。但在现今信息化作战的背景下，这一方法存在着精度有限、通用性差、设计过程繁琐等缺点，尤其是只能对各个弹种的具体弹道和各个修正量分别进行解算，极大地制约了火控系统的开发[7]。对此，基于弹道微分方程组的火控外弹道模型在该背景下应运而生，它具有通用性好、精度高等特点，但解算处理复杂，耗费时间较多，尤其是对动目标射击时，需以逆解法[6]为工具，通过试探法不断迭代求解出命中点的估计值，因此计算量更大。

目前工程应用中所选用的弹道微分方程组主要是质点弹道方程组和修正的质点弹道方程组[8]。质点弹道方程组是将弹丸作为质点处理，忽略了作用在弹丸上的升力、马氏力和攻角引起的诱导阻力，使弹道计算结果产生了一定的误差；修正质点弹道模型在质点弹道模型的基础上综合了这些因素，有效地提高了计算精度。然而对于高转速弹丸，修正质点弹道模型在偏流的计算上存在着一定的误差。对此有学者将弹丸在空中的运动视为一般刚体的运动，考虑弹丸飞行期间产生的进动、章动和自转等，极大地提高了对弹丸轨迹的预测精度[9]，但复杂的刚体弹道方程解算很难在实际火控系统的一个工作周期内完成。

鉴于此，本文将刚体弹道方程引入到防空火控系统的弹道模型中，极大地提升了火控系统整体的通用性与精确性。同时为了避免刚体弹道方程计算时给火控带来的实时性问题，本文从火控解算的方法入手，将无偏转换得到的噪声协方差矩阵做解耦，降低了目标状态估计时滤波的复杂度；并在割弦迭代法中，将虚拟脱靶量引入到射角和射向的迭代- 修正过程中，减少了火控解算中对弹道方程计算次数，加快射击诸元求解的收敛速度。从而在算法执行机制上，有效地克服了基于刚体弹道方程火控解算实时性不足的缺陷。

1 弹道建模

由于影响弹丸在空气中运动的因素极为错综复杂，本文针对某型防空高炮，根据各个因素对弹丸飞行影响的大小，在六自由度刚体弹道方程的基础上[10],忽略科式惯性力和赤道阻尼力矩；同时为了增大龙格- 库塔法求解时的步长，忽略弹轴高低角和方位角的2阶导项，建立如下弹道方程：

(1)

2 脱靶量

如图1所示，以武器回转中心O为坐标原点构建地面基准坐标系Oxyz，其中：Ox轴为射面与炮口水平面的交线，指向射击方向为正，Oy轴在射击面内，垂直Ox轴，向上为正；弹目偏差观测坐标系TxEyEzE以目标中心T为原点，沿着观目矢量OT为弹目偏差观测坐标系的纵深轴yE，过OT的铅锤面与迎光面(也叫Q平面)的交线为高低轴zE，向上为正，在迎光面内，垂直于zE轴的规定为方位轴xE，方向满足右手法则[11-12]；当弹丸以速度v发射时，记T*为目标T在Oxz面上的投影，Z*为弹迹穿过迎光面的穿越点，Δβ和Δε分别表示观目矢量OT与观炸矢量OZ*之间张角的水平分量和高低分量，βh为观目矢量在东、北、天地理坐标系下的方位角，εh为观目矢量在东、北、天地理坐标系下的高低角。则弹目偏差TZ*为

TZ*=(xE,0,yE,0,zE,0)=
(|OT|tan Δβ,0,|OT|tan Δε).

(2)

图1 脱靶量示意图Fig.1 Schematic diagram of miss distance

而在火控系统中，弹目偏差的数据处理都是建立在东、北、天地理坐标系下的[6]。因此，实际地面基准坐标系下的脱靶量Δβ′、Δε′分别为

(3)

3 快速防空火控解算方法

3.1 解耦无偏转换卡尔曼滤波

为了准确描述目标状态随着时间变化的过程，将目标位置、速度和加速度作为状态量，构建目标的运动方程

Xk=Fk-1Xk-1+Wk-1，

(4)

火控系统中对目标状态的估计离不开探测设备对目标的跟踪，在探测系统中，传感器的量测信息通常是建立在三维球坐标中，目标与传感器之间的相对位置如图2所示。图2中，rk、ηk、βk分别为径向距离、高低角、方位角。

图2 传感器测量坐标系Fig.2 Sensor measurement coordinate system

传感器的量测值Zk主要包括径向距离rk、高低角ηk和方位角βk，其测量方程[13-14]为

(5)

式中：Vrk、Vηk和Vβk是相互独立、均值为0且恒定方差为σrk、σηk和σβk的高斯白噪声。

球坐标可通过(6)式转化到笛卡尔坐标系下的量测：

(6)

因此，根据(5)式的量测值Zk，对真实均值和协方差矩阵求数学期望得到无偏量测偏差μk和协方差Rk[15]：

(7)

(8)

式中：

(9)

3.2 基于刚体弹道模型的快速射击诸元求解

3.2.1 命中方程建立

假设观瞄设备中心O′与武器回转中心O不在同一位置，记J=OO′为观炮基线，则实际防空武器的命中多边形如图3所示。

图3 命中多边形Fig.3 Hit polygon

相应的命中矢量方程为

Dh(t)-D(t)-J-Sh(t)=0，

(10)

式中：Dh(t)=D(t+tf)，tf为弹丸飞行时间。

为了计算方便，假设目标在弹丸飞行时间内做匀速直线运动，则

(11)

通过将(10)式中的矢量在东、北、天地理坐标系上进行投影，得到标量命中方程组

(12)

式中：(xoo′,yoo′,zoo′)和(xTh,yTh,zTh)分别为瞄具中心O和命中点Th的坐标；瞄准矢量D(t)=(T,T,T).

由(1)式可知，当射角φh和弹丸飞行时间tf已知时，利用4阶龙格- 库塔法可求得弹丸轨迹上任意一点坐标，即

(13)

再根据射向方位角βh，命中点坐标可写为

(14)

由(12)式、(13)式和(14)式即可求得射击诸元(φh,βh,tf)。

3.2.2 基于虚拟脱靶量的迭代——修正算法

考虑到火控解算的精度与计算速度等要求，可利用割弦迭代法计算出命中点坐标，在每次割弦过程中，需利用多次修正龙格- 库塔法求解弹道模型时的初值，才能计算出符合当前预测未来点的射击诸元。因此，本文引入虚拟脱靶量概念，即目标真实航迹与计算机仿真弹丸轨迹间偏差量，对初始射击射角和射向进行修正，极大地提高了计算当前预测未来点射击诸元的收敛速度。具体求解过程如下：

步骤2割弦迭代法初值的选取。假设在弹丸飞行时间内目标处于静止状态，即目标运动时间Tf=0 s，则命中点坐标为(xTh,k,yTh,k,zTh,k)，有

(15)

这时可将火控解算问题转化为一个外弹道学中的两点边值问题，即从上述边值条件出发，求出射击诸元。

步骤3弹道方程初值给定。根据所给的目标未来点坐标预估初始射角φh，0和射向方位角βh，0:

(16)

步骤4弹道方程求解。利用4阶龙格- 库塔法求解微分方程组(1)式，得到地面基准坐标系下弹道诸元的序列值：

{xi,yi,zi,vx,i,vy,i,vz,i,tf,i},i=0,1,2,…，

(17)

式中：i为龙格- 库塔法定步长求解时的序列号。

步骤5虚拟脱靶量计算。

1)将步骤4得到的弹道诸元序列通过(18)式转化到弹目偏差观测坐标系中，

(18)

2)利用拉格朗日插值，可计算出在Q平面内的二维弹目偏差(xE,0,zE,0)和弹丸飞行时间tf,i,0，利用(3)式将(xE,0,zE,0)转换到东、北、天地理坐标系下的虚拟脱靶量(Δβ′0,Δε′0)。

步骤6射角与射向迭代修正。判断(Δβ′0,Δε′0)是否满足系统所需的精度(一般取0.1 mil)。若满足，则射击诸元为φh=φh,0，βh=βh,0，tf=tf,i,0，并转到步骤7；否则对初始射角和射向进行修正，得到新的射角和射向方位角：

(19)

将修正后的射角和射向取代原来的变量，并返回步骤4.

至此，一个周期内的火控解算完成，相比传统火控解算的迭代修正算法，步骤5中利用虚拟脱靶量进行射角和射向的迭代修正，减少了火控解算中对弹道方程计算次数，极大提高了基于刚体弹道方程的防空火控解算方法求解速度。

4 仿真分析

4.1 仿真场景

本文以某型防空高炮为研究对象，构建三维仿真场景，对所提火控解算方法进行仿真验证。假设目标做匀速运动，起始位置为(3 500 m, 3 500 m, 1 000 m)，初始速度为(-300 m/s, -200 m/s, -25 m/s)；测量传感器采用雷达，采样周期为0.02 s，距离误差为5 m，方位和俯仰角误差为1 mrad，目标在雷达界面的显示如图4所示。

图4 目标在雷达界面显示Fig.4 Target on radar interface

弹丸初速v0=1 050 m/s，龙格- 库塔法的解算步长为0.01 s，在标准气象条件下，不同射角的弹丸轨迹如图5所示。

图5 弹丸轨迹Fig.5 Projectile trajectory

4.2 仿真分析

4.2.1 状态滤波

在实际工程中，通过传感器测量的目标运动参数包含众多干扰噪声，这些噪声不仅影响系统的计算精度，严重时系统将失稳。因此对目标进行状态滤波是准确预测目标未来点位置的关键技术。本文通过无偏转换将雷达对目标运动状态的非线性量测方程围绕目标运动状态进行展开，并将无偏转换得到的噪声协方差矩阵做解耦，降低了滤波算法计算的复杂度，其位置均方根误差(RMSE)、速度RMSE和滤波处理时间如图6、图7和图8所示。由图6～图8可以看出，解耦无偏转换卡尔曼滤波与无偏转换卡尔曼滤波算法[15]的误差差值不大于总误差的1%，而滤波时间减少了15%. 因此，该方法具有较高的精度和良好的性能。

图6 位置RMSEFig.6 RMSEs of position

图7 速度RMSEFig.7 RMSEs of velocity

图8 滤波的处理时间Fig.8 Processing time of filtering

4.2.2 火控解算误差

根据上述仿真场景，以2 s稳定建航后，分别用传统的火控解算方法[6]和本文所提方法对不同时间点进行火控解算，其解算点的射击诸元对比如表1所示，火控解算的误差如图9所示。

从图9仿真结果可以看出，本文所提方法的有效性和可行性。由表1可以看出：当目标在较远处，其解算误差很大，但随着目标越来越接近，弹丸飞行时间越来越短，误差呈现减小趋势，这是由于测量传感器所得目标数据存在着一定的量测误差，导致解算误差大部分来源于滤波误差；当目标距离3 000 m以内时，单轴方向上的解算误差不大于0.5 m，该方法可以很好地满足实际防空作战需求。

表1 射击诸元对比Tab.1 Comparison of firing data

图9 火控解算误差Fig.9 Errors of fire control calculation

4.2.3 解算时间

本文通过对不同距离目标的火控解算进行仿真实验来有效地说明所提出方法的实时性优势。在Matlab2014环境下(Windows 7, i5-3210M CPU, 2.5 GHz主频, 8 GB RAM)将传统的火控解算方法[6]与本文方法进行比较，解算时间对比如表2所示。

表2 Windows平台下的解算时间对比Tab.2 Comparison of calculation times based on Windows s

从表2中可以看出，本文所提方法在计算量上远远低于传统的火控解算方法[6]，这是由于火控诸元的解算主要包括状态滤波与命中点解算两大步骤。对于状态滤波，通过将无偏转换卡尔曼滤波的噪声协方差做解耦，一定程度上降低了滤波计算的复杂度；对于命中点解算，本文提出基于虚拟脱靶量的迭代——修正算法，利用虚拟脱靶量进行射角和射向的迭代修正，减少了火控解算中对弹道方程计算次数，大幅度提高火控求解速度。

为了有效说明该方法在工程时是可行的，将代码改编成C语言并移植到配置为i7-4790CPU、SylixOS实时操作系统的工控机上运行，解算时间对比如表3所示。从表3实验的数据中可以看出，本文所提方法在工程上具有可观的应用前景。

5 结论

为了提高防空火控解算的精确性与实时性，本文提出了基于刚体弹道方程的防空火控解算方法，仿真结果证明了该方法的有效性和可行性。然而，对于实际的防空系统，目标的强机动特性使得基于传统运动模型的滤波算法出现精度下降甚至发散现象，因此后续的工作是考虑在该类情况下，提高目标状态估计与轨迹预测的准确性，并在一定的基础上降低算法的复杂度，使得该算法可用于未来防空武器的火控系统中。

表3 SylixOS下的解算时间对比Tab.3 Comparison of calculation times based on SylixOS ms