基于外弹道的防空火控改进算法及模型仿真

李 智,吴盘龙

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

在现代计算机仿真技术的加持下,火控解算已经基本可以摆脱射表进行[1]。文献[2]给出了使用龙格-库塔算法对外弹道微分方程组进行实时积分解算的流程,并通过Matlab仿真验证。文献[3]解决了给定位置解射击诸元与解命中问题,并与射表数据对比,得到较好结果。文献[4]针对弹道解算的边值问题提出了迭代-修正法,使弹道问题可以脱离射表进行。文献[5]基于二分法探讨了一种射击诸元的解算方法,缩短了弹道解算时间。传统的射击诸元解算流程繁复、低效、误差显著,对此,本文提出了改进的射击诸元解算的迭代-修正流程,并基于Simulink模块搭建了防空火控解算系统。

1 防空火控解算系统设计

防空火控解算系统作为末端防御系统综合控制技术研发环境的重要组成部分,主要用于打击典型目标的射击诸元解算,通过对典型目标航路生成以及武器弹道解算的研究,为火控防御系统的开发提供有利环境,其不仅有效地降低了系统的研制成本,而且提高了拦截武器的研制成功率。

防空火控解算系统主要包括坐标转换模块、外弹道模拟模块、目标信息模拟模块、火控解算模块等,其整个框架如图1所示。

图1 防空火控解算系统

坐标转换模块用来统一不同坐标系下弹体运动方程与目标运动方程,方便后续解算;外弹道模拟模块通过弹道方程、气象模拟和弹道解算来模拟一个或多个完整的弹体运动轨迹;目标信息模拟模块用来产生火控解算所需要的典型运动目标的位置、速度、加速度等信息,并可以将其与外弹道模拟结合来判断是否命中目标;火控解算模块根据目标实时运动状态,通过解命中方程组,不断计算未来时态弹目相会点和当前时态火箭弹射击诸元。

2 基于外弹道的火控算法设计

2.1 坐标系定义与坐标转换模块

在火箭弹外弹道理论中,弹体的数学模型比较复杂,必须采用合理的坐标系,以保证获取准确的弹道轨迹和目标航迹。建立火箭弹运动方程组所涉及的坐标系繁多,为了方便理解,对本文火控解算中主要的坐标系定义如下。

1)发射坐标系O-XVYVZV

O:武器载体中心为原点;XV:沿着水平线指向射击方向为正;YV:垂直于水平面方向铅直向上为正;ZV:垂直于射击面且按右手法则确定其正方向。

2)东北天坐标系O-XEYEZE

O:武器载体中心为原点;XE:水平面内指向正东方向为正;YE:水平面内指向正北方向为正;ZE:垂直于水平面方向,铅直向上为正。

弹道微分方程组主要基于发射坐标系,该坐标系可以用于确定火箭弹质心的空间坐标位置和火箭弹在空间的姿态等。目标运动和命中方程组主要基于东北天坐标系,该坐标系用于描述目标运动的空间坐标位置和火箭弹射击诸元等。所以,需要将弹道解算数据和观测目标数据统一到东北天坐标系上,发射坐标系与东北天坐标系关系如图2所示。

图2 发射坐标系与东北天坐标系转换关系

其中,火箭弹发射高低角(射击方向与XV轴的夹角)为θ,火箭弹发射方位角(射击方向与YE轴的夹角)为φ。东北天坐标系可由发射坐标系经两次旋转而成。第一次是发射坐标系先绕YV顺时针转动90°-φ,第二次是再绕旋转后的XV轴顺时针转动90°到东北天坐标系[6]。发射坐标系与东北天坐标系转换矩阵如下:

(1)

2.2 外弹道模拟模块

2.2.1 弹道方程

弹道方程方面,本文中的防空武器采用某型70 mm防空火箭弹,假设火箭弹仅在纵向平面内飞行,以弹体飞行时间t为自变量,弹体运动微分方程组为[7]:

(2)

其中,x、y为弹体质心坐标,v为弹体速度,m为弹体质量,g为重力加速度;α为弹体攻角,θ为弹道倾角,ϑ为俯仰角,且有ϑ=α+θ,ωz为俯仰角速度,Jz为绕z轴转动惯量。包含的力和力矩方程为:

(3)

其中,X为弹体受到的阻力,是总空气动力在飞行速度方向上的分力,指向与速度v相反;Y为受到的升力,在弹轴和速度v所确定的阻力面内,垂直于速度v且与v同侧;Mz为静力矩,能使攻角增大的称为翻转力矩,能使攻角减小的称为稳定力矩,统称为静力矩,其方向垂直于攻角平面;P为发动机推力,Isp为发动机比冲,mp为推进剂质量,tp为推进时间;ρ为空气密度函数,S为横截面积,L为特征长度;Cx、Cy为阻力系数、升力系数,依据弹体马赫数,可以采用分段线性插值算法得到:

(4)

其中Ma是弹体在当前时刻的马赫数,Ma1、Ma2为弹体阻力系数表中Ma所在区间范围的下界和上界。Cy同理可得。

2.2.2 弹道方程解算

外弹道运动方程组建立完成后,需要对其进行求解。龙格-库塔法(Runge-Kutta)是求解常微分方程的常用迭代算法[8],弹道方程的数值解法采用定步长四阶龙格-库塔法,四阶龙格-库塔法求解微分方程过程如下。

为方便描述,将弹道微分方程和初值记为

(5)

其中,Y=[v,θ,x,y,ϑ,m]T,Y0为方程解算的初始值。四阶龙格-库塔法迭代公式如式(6)所示。

(6)

对于四阶微分方程,通过引入新的变量k1、k2、k3、k4,将原方程转化为较为简单的一阶微分方程,然后进行求解。通过调整龙格-库塔法步长的大小来调整解算的时间。调整时首先要注意满足时间要求,然后尽可能地增大精度。通过该算法可以解算任意射角下和任意弹体飞行时间下的弹道诸元。

2.3 火控解算模块

2.3.1 命中方程组

在火控解算中,命中点指的是火箭弹和目标的相遇点,通过求解弹道轨迹与目标轨迹的交点,得到命中点坐标和火箭弹飞行时间[9]。一般情况下,我们假设观测瞄具中心与武器载体中心同处于一点,在这种情况下,目标运动时,瞄准矢量、命中矢量与提前矢量将构成一个三角形,被称为命中三角形[10],命中点问题如图3所示。

图3 命中三角形示意图

由于现代火控计算机的计算形式是数字而非矢量的,所以,需要在标量的情况下进行解算,并做出以下基本假设:

1)以武器载体中心为原点建立坐标系;

2)发射时刻目标所在位置(即目标现在点)为M0,其坐标为M0(x0,y0,z0);

3)目标从初始位置运动Tm时间到命中点,弹体从发射到运动到该点的飞行时间为Td,命中时,Tm=Td,命中点即相遇点M,其坐标为M(xm,ym,zm)。以此建立目标命中方程为:

(7)

当发射高低角θ0和弹飞时间Td已知时,利用四阶龙格-库塔法可求得弹体运动轨迹上任意一点坐标:

(8)

联立得命中方程为:

(9)

式(9)即为在直角坐标系下命中点坐标与弹飞时间的非线性方程组,对该方程组求解即可解决命中问题。

2.3.2 射击诸元解算

在命中三角形中,M点为命中点,武器载体中心O′为原点,射击诸元解算是两点边值问题。从该边值条件出发,解算出射击时所需要的高低角θ、方位角φ和火箭弹飞行时间Td,是本章防空火控解算系统的任务。为了获得较快的收敛速度,将边值问题化为初值问题求解。这里选择现在时刻弹目视线角(目标现在点M0和武器载体中心原点O′连线与水平面夹角)为高低角初值进行解算。主流的求解方法中,迭代-修正法具有较快的收敛速度。

根据文献[3]的传统迭代算法,射击诸元求解的流程首先视目标位置为固定点对高低角进行迭代,当火箭弹高度zd与目标现在点z0的差值小于一个设定好的数值ez,则视为命中固定位置目标;然后在目标运动后,根据运动时间对目标位置进行迭代,当弹目距离D小于一个设定好的数值ed,则完成对运动目标的命中问题求解。即使其对时间的迭代方式有所改进,但总体上来说该算法依然效率较低,还易将高低角迭代过程产生的误差引入位置的迭代过程中,在复杂环境下易使结果发散。根据文献内容完善其对运动目标射击诸元求解算法的流程如图4所示。

图4 基于传统的迭代-修正法射击诸元解算过程

由于传统流程在求解针对动目标的射击诸元时,将目标先视为静止状态,基于当前状态设定发射初值,迭代解出修正高低角;再考虑运动状态下目标飞行时间,显然此时目标已经飞离原位置;继而,再次将目标视为静止状态,基于当前的目标位置设定更新后的发射初值,再次迭代解出修正高低角,如此往复循环。可以看到,为了解算射击诸元,传统流程中包含了多重循环,众所周知,循环层数对算法速度有着至关重要的影响,并且可能会因多重循环的误差积累导致结果发散。

图5 基于改进的迭代-修正法射击诸元解算过程

在流程图中,弹目高度偏差值ez、弹目终点距离偏差值ed根据步长与解算精度的要求设定,两者的设定对解算结果是否合理以及解算是否快速有影响。迭代初值通常设置为

(10)

龙格-库塔法解弹道的结束判断条件为:静目标下火箭弹横坐标xd=x0,动目标下火箭弹横坐标xd=xm,当满足此条件时结束弹道解算过程。微分方程组解算积分步长为T,到弹道终点的积分步数为n,则此时火箭弹飞行时间Td和目标运动时间Tm为

Tm=Td=(n-1)*T

(11)

目标终点的坐标为(xme,yme,zme),目标终点的高低角为

(12)

火箭弹终点的坐标为(xde,yde,zde),火箭弹终点的高低角为

(13)

则高低角修正为

θi+1=θi+θme-θde,i=0,1,2,…

(14)

迭代过程结束的条件有两个,当满足任意一个,就结束射击诸元解算过程。

结束条件1:设置迭代精度ez、ed,当满足下式:

|zd-z0|

(15)

D

(16)

说明解算的结果达到设定的精度需求,结束对应的迭代解算过程。

结束条件2:为了避免因精度设置过小或结果发散等导致的迭代-修正次数过多,保证火控解算在合理有效时间内完成,设置了迭代次数上限,当迭代次数达到设定次数上限时,结束迭代解算过程。

3 模块搭建与系统仿真

针对防空火控命中问题,为了有效说明经过改进的迭代-修正法进行射击诸元解算在实时性方面具有明显的提升,在仿真系统中使用改进的射击诸元解算方法与传统方法进行火控解算和仿真,对结果进行比较,解算情况如下。

根据目标模型、外弹道模型及火控解算算法,得到Simulink仿真模型[11]如图6所示。“火控解算”模块具体如图7所示。“弹体动力学模块”具体如图8所示。“坐标系转换”具体如图9所示。

图6 火控解算Simulink模型

图7 火控解算模块

图8 弹体动力学模块

图9 坐标系转换模块

采用龙格-库塔法在仿真系统中解算某型70 mm火箭弹外弹道方程组,积分步长为0.01 s,仿真初始参数如表1所示。仿真结果如图10所示。

表1 火箭弹火控解算初始参数

图10 弹道轨迹和目标运动轨迹三维视图

1)基于传统算法进行解算

雷达采样周期为0.02 s,解算步长为0.01 s,总迭代次数i=9,在最终脱靶量D=0.940 m情况下,解算得到射击诸元:高低角226.700 mil,方位角422.621 mil,命中时间Tm=Td=7.44 s。Matlab平台下解算所需时间为14. 195 392 s。具体如表2所示。

表2 传统算法迭代过程数据

可以看到,由于该算法迭代分两步进行,只有当修正值Δθ达到一定精度时,方可对弹飞时间及弹目距离进行计算,因而其迭代效率低下,并且由于算法精度有限,第9次迭代时数据出现重复,增大迭代次数上限也不会使结果更加精准,而此时命中精度只有1 m左右,无法满足高精度火控需求。

2) 基于改进算法进行解算

雷达采样周期为0.02 s,解算步长为0.01 s,总迭代次数i=4,在最终脱靶量D=0.220 m下,解算得到射击诸元:高低角226.987 mil,方位角422.620 mil,命中时间Tm=Td=7.44 s。Matlab平台下解算所需时间为5. 791 897 s。具体如表3所示。

表3 改进算法迭代过程数据

可以看到,由于该算法同步进行迭代,无论当修正值Δθ达到精度与否,均对弹目关系进行计算,因而其迭代效率较高,并且算法精度相对较好,第4次迭代后结果便达到了可观的精度,两者弹飞时间在一定精度内可以看作基本相同,说明了本文改进算法的精确性。

进一步,为了说明改进算法在火控解算中实时性较好,针对不同起始位置的运动目标分别在Matlab环境与SylixOS实时操作系统工控机上进行了火控解算和仿真,多次对比两种算法所需解算时间。仿真结果如表4与表5所示。

表4 Windows平台下两种算法耗时对比

表5 SylixOS平台下两种算法耗时对比

从表中可以看出,本文引入的改进算法在多种情况下解算时间均远低于传统的解算方法,且算法的实时性较好,在实际工程中有一定的应用价值。

4 结束语

本文设计了防空火控解算系统,搭建了坐标系转换模块、目标信息模拟模块、外弹道解算模块、火控解算模块的内容。然后针对传统射击诸元解算流程繁复、低效、误差显著的特点,改进了解算流程。最后对改进算法模型进行仿真,分析证明了改进算法有效性和可行性。

对于存在目标探测误差的防空火控系统,目标轨迹偏差会使得基于典型运动模型的算法精度下降,下一步是针对目标航迹的状态滤波研究,并在此基础上继续研究火控解算系统。