导弹发射包线指数优化搜索仿真分析

2020-11-16胡东愿刘会亮岳龙飞杨任农左家亮

宇航学报 2020年10期

胡东愿，刘会亮，岳龙飞，杨任农,左家亮

(空军工程大学空管领航学院,西安 710051)

0 引言

导弹发射包线是评估导弹性能[1]和装备作战能力[2]的一个重要指标，也是飞行员控制战斗机隐蔽接敌[3-4]、机动决策[5]和目标拦截[6-7]的重要参考。在空战敌我态势快速变化的情况下，导弹发射包线也为飞行员及时掌握战场情况、预判空战结果、占据有利位置提供有效依据。因此对导弹发射包线解算的精度和实效性提出了更高的要求。

导弹发射包线即为导弹攻击区，是以敌机为中心，载机可以发射导弹的区域，主要由远边界和近边界构成。目前计算边界的主要方法有运动模拟法[8-10]、多项式拟合法[11]和查表插值法[7, 12]。

运动模拟法是其他攻击区解算方法的基础，通过弹道模型、目标运动模型和相对运动模型仿真计算得到。在边界值搜索时，经常会遇到计算量大、耗时长、精度低等问题。文献[13]将导弹模型简化为二维模型，利用二分法搜索包线值，计算速度有很大的提高，但与真实三维情况相差较大。文献[14]提出了近似式的快速生成方法，改进搜索初始值的选取，在一定程度上减小了搜索空间。但近似式中仅考虑起始方位角和初始高度，对其他影响攻击区的元素并没有考虑，使其搜索空间依然较大。文献[9, 15]提出自适应调整仿真计算步长，减少仿真计算次数，从而减小总体计算时间。但在每一个仿真步长中都使用简单的二分法或黄金分割法，每次搜索所需迭代次数和迭代时间并未改变。

为解决攻击区仿真解算中迭代次数多、计算量大、耗时长问题，在前人研究基础上，本文首先建立导弹和目标运动学模型，利用simulink搭建仿真模型；在边界值搜索时，以黄金分割法为基础，提出指数优化搜索算法，并与黄金分割法进行对比，分析算法性能。最后将算法运用到攻击区边界值搜索上，利用实验结果分析运动参数对攻击区的影响。实验表明，在给定的误差范围内，计算包线边界值时所需的迭代次数更少，在一定程度上减小了搜索时间。在相同的迭代计算次数下，能获取更高的精度，具有一定的工程价值。

1 运动学建模

1.1 坐标系的定义与转换

地面坐标系A-xyz：地面坐标系是与地球表面固连的惯性坐标系，随地球的自转而旋转。Ax轴为航迹面与水平面的交线，指向目标。Ay轴沿垂线向上，Az轴与其他两轴垂直并且构成右手坐标系。

导弹坐标系O-x2y2z2：弹道坐标系是原点取在导弹瞬时质心的动坐标系，随着导弹运动而运动。Ox2与导弹的速度矢量v重合，Oy2为包含v的铅锤面内并垂直于Ox2指向上，Oz2轴与其他两轴垂直并且构成右手坐标系。

两个坐标系可以通过旋转互换得到，如图1所示，旋转的角度θ为弹道倾角，即导弹的速度矢量(Ox2)与水平面Axz的夹角。ψV为弹道偏角，即导弹速度矢量在水平面内投影与Ax轴的夹角。

图1 地面坐标系与导弹坐标系示意图Fig.1 The earth and missile coordinate system

旋转矩阵如下：

1.2 导弹模型构建

导弹模型主要有三自由度[8, 16]和六自由度[17-18]两种，文献[16]表明三自由度模型解算的时间代价比六自由度模型小，但三自由度模型中，每一时刻的导弹都会有较大的偏差。文献[13]将快速计算法、三自由度和六自由度进行对比，结果表明六自由度模型在导弹命中概率上更为准确。考虑到目标机做大过载强机动飞行，导弹在跟踪过程中姿态角将会受到很大的影响，为尽可能减小命中误差，文中构建导弹六自由度模型。

构建导弹六自由度运动模型前，首先做出如下假设：

1)忽略导弹飞行时的次要因素对运动的影响，如弹体结构变形。导弹在高速飞行时，可以看作是薄翼的细长体的弹性结构，在气动力和结构弹性的作用下，会产生弹体微小的弹性和塑性变形。但对于弹道运动规律的研究来说，这些产生微小变形的力相比较导弹的重力和推力可以忽略不计，对运动轨迹不产生影响。

2)将导弹看作一个刚体时，可以将刚体的运动分解为质心的平移运动和绕质心转动。计算轨迹只需关注质心瞬时位置的三个自由度和刚体瞬时姿态的三个自由度。导弹的速度矢量方向与导弹轴近似重合，不考虑导弹的迎角和侧滑角。

3)飞行过程中，燃料不断消耗，导弹质量减小，假设在质量变化的过程中导弹的重心位置不变。根据“固化原理”，将导弹质量和燃气质量近似为一个整体，转化成一个常质量系。瞬时变质量系的导弹力学方程可以写成常质量的形式，因此任意时刻的质心可以用初始时刻的质心近似代替。

1.2.1导弹动力学模型

工程实践表明：以导弹坐标系为基础，对与研究导弹质心运动，方程的形式最为简单，便于分析运动特征。

导弹质心运动的动力学方程：

(1)

式中：α为攻角，β为侧滑角，ψV为弹道偏角，γV为速度倾斜角，X,Y,Z分别为阻力，升力和侧向力，P为推力。

导弹绕质心转动的动力学方程为：

(2)

式中：Jx1，Jy1，Jz1分别为导弹各轴的转动惯量；ωx1，ωy1，ωz1为转动角速度在弹体坐标系上各轴的分量；Mx1，My1，Mz1为作用在质心上的力对质心的力矩在弹体坐标系各轴上的分量。

1.2.2导弹运动学模型

以地面坐标系为基础，根据导弹质心运动轨迹，建立导弹质心运动学方程。

(3)

式中：(x,y,z)为导弹在地面坐标系中的位置坐标，V为导弹的速度，θ为弹道倾角，ψV为弹道偏角。

根据地面坐标系与弹体坐标系的转换关系可知：

(4)

式中：ψ为偏航角，ϑ为导弹俯仰角，γ为倾斜角，三者构成导弹的姿态角。

导弹绕质心转动方程为：

(5)

1.2.3质量变化方程

导弹在飞行过程中，由于发动机不断地消耗燃料，导弹质量不断减小，其变化可以简化如下：

(6)

(7)

式中：mc为燃料质量秒流量，即导弹单位时间内质量消耗量;m0为导弹的初始质量。

1.2.4导弹机动过载方程

导弹的机动性能是导弹飞行性能的重要参考指标，与导弹的运动有密切的联系，也是影响导弹攻击区的重要因素。导弹在飞行过程中的作用力以及加速度的变化情况可以用过载来衡量。

将过载矢量n向弹道坐标系上投影得各分量nx2，ny2，nz2如下：

(8)

由式(8)可推断在特殊情况下的过载，如导弹在铅锤平面内飞行时，nz2=0；在水平面内飞行时，ny2= 1；做直线飞行时，ny2= cosθ为常数，nz2=0；做等速直线飞行时，nx2= sinθ为常数，ny2= cosθ为常数，nz2=0；做水平直线飞行时，ny2= 1，nz2=0；做等速水平直线飞行时，nx2= 0，ny2=1，nz2=0。

1.3 目标运动方程

以地面坐标系为基础，目标运动方程如下：

(9)

(10)

式(9)中：(xt，yt，zt)为目标在地面坐标系中的坐标，vt，θt，φt为目标速度，航迹俯仰角，航迹偏转角。式(10)中:ntx，nty，ntz为目标纵向控制过载，偏航控制过载和转弯控制过载。

1.4 相对运动方程

比例导引法是指导弹速度矢量的旋转角速度与目标线的旋转角速度成比例的一种导引方法。弹道前段较弯曲，能充分利用导弹的机动性能；导弹后段较为平直，使导弹有充裕的机动能力。比例导引法是采用比较广泛的一种导引方式，其导引方程可表达为

(11)

式中：K为比例系数，如图2所示，式中q为目标线方位角，是目标线与基准线之间的夹角，从基准线逆时针转到目标线，则q为正。σ为导弹弹道角，为导弹速度矢量与基准线之间的夹角。η，ηT为导弹速度前置角和目标速度矢量前置角，分别为导弹、目标速度矢量与目标线之间的夹角。q，σ，η，ηT从基准线开始，逆时针转动到相应的夹角线上则为正值，反之顺时针转动则取负值。

导弹与目标之间的相对运动方程为：

(12)

图2 比例导引示意图Fig.2 Proportion guidance of missile

2 指数优化搜索

2.1 优化初始空间的黄金分割边界搜索法

黄金分割法适用于搜索区间上的任何函数求极值问题，对函数没有太多严格要求，适用范围较广。该算法通过选取试探点，使包含极值的空间不断缩短，直到区间的长度和所求值满足精度要求。

本文以黄金分割法为基础，仿真计算导弹轨迹，试探搜索攻击区的远边界。在算法的开始最主要是选择含有极值点的初始搜索空间，搜索空间距离太大会增加算法的计算次数，选择的空间太小则很容易错过极值点。利用导弹运动参数性能值，近似估计导弹包线边界的可能范围，为搜索的初始化提供启发式信息，可以有效避免搜索初始空间的随机性。

导弹能否追踪到目标并将其截获命中主要由两种因素决定，静态因素和动态因素。静态因素主要反映导弹的固有性能，如导弹的机动性能，携带燃料能支持的最大飞行时间，最大航程等。当目标发现被导弹跟踪时，会采取转弯机动的方式进行规避，机动性能决定了导弹能否及时调整前置角跟随目标转弯。目标转弯逃离时，最大飞行时间决定了导弹可跟踪目标的航程，在跟踪末端，因燃料的消耗而失去动能。若此时的弹目距离仍不在杀伤范围内，导弹将会丢失目标。导弹的固有性能随其生产而确定，变化不大。可用最大过载、导引头最大偏角、导引头最大跟踪角速度、导弹最大可控飞行时间来量化衡量。

动态因素主要反映当前两机所处的态势。文献[18]中主要根据两机的速度、相对距离、高度、进入角等参数将态势划分为4类：我机优势、敌机优势、双方中立、双方均势。双机的战场态势形成的主要原因在于敌方与我方导弹攻击区的空间关系，也可看作是导弹攻击区的一种表现形式。两机的速度和高度通过影响攻击区的包线从而决定了态势的变化。

静态因素和动态因素共同决定了导弹可攻击区的形状和大小。利用导弹最大过载、导引头最大跟踪角速度、导引头最大偏角、导弹最大可控飞行时间、进入角、初始速度、初始高度、两机相对速度、相对高度8个变量来近似拟合包线的搜索空间，初始化搜索值，减少盲目性。利用八元组Λ=[nmax,ωmax,tmax,ψ,v0,H0,ΔH,Δv]建立多项式：

R(pi)=(vmax·tmax·sinψ)p1+(H0·cosψ+

ΔH·nmax)p2+(ωmax·H0·tmax+Δv·tmax)p3

式中：nmax为导弹可用最大过载，ωmax为导引头最大偏角，tmax为导弹最大可控飞行时间，ψ为初始进入角，v0为导弹初始速度，H0为导弹初始飞行高度，ΔH为导弹目标的相对高度，Δv为导弹目标的相对速度。

根据每次搜索的结果对多项式系数进行调整，充分利用计算值，优化下一组初始的搜索空间。算法的步骤如下：

(1)利用八元组确定拟合多项式，确定初始搜索空间[a1,b1]。

(2)计算试探点λ1，μ1，在指定空间的0.382和0.618处取值，λ1=a1+0.382(b1-a1)，μ1=a1+0.618(b1-a1)，仿真计算两个点处的导弹轨迹。

(3)判断是否满足击中条件，若λ1满足，μ1也满足，转到(4);若λ1满足，μ1不满足，转到(5)；若λ1不满足，转到(6)。

(4)a1=μ1，b1不变，代替原搜索空间。

(5)a1=λ1，b1=μ1，代替原搜索空间。

(6)a1不变，b1=λ1，代替原搜索空间。

(7)判断该搜索空间是否满足精度要求，若满足，结束;若不满足，转到(2)。

2.2 指数优化边界值搜索

传统黄金分割法的区间长度缩短比率为常数，故其收敛速度较慢，需要计算的函数值也较多。考虑到实时性的要求，通过改变以指数收敛的区间长度缩短比率得到一种新的一维搜索指数优化算法。指数优化法和黄金分割法同样也只适用于单峰函数。在计算过程中，第一次迭代需要计算两个试探点，以后每次迭代只需要计算一点，另一点取自上次迭代，其与黄金分割法的主要区别在于区间的长度缩短比率不是以常数进行收敛，而是计算两次区间差的指数进行收敛，由于指数函数的引用，因此构造出一个具有更加光滑的收敛趋势并且更快收敛速度的指数优化函数。

指数优化搜索主要有以下步骤：

(1)给定初始区间[a1,b1]及精度要求ε>0，计算试探点λ1，μ1，有λ1=a1+0.382(b1-a1)，μ1=a1+0.618(b1-a1)。然后，计算λ1，μ1是否满足命中条件，置y1=b1-a1,e1=0，并令k=1。

(2)若bk-ak<ε，停止计算；否则当λ1满足时，转入(3)；当μ1满足时，转入(4)。

(3)置ak+1=λ，bk+1=bk，λk+1=μk，yk+1=bk+1-ak+1，ek+1=yk+1-yk，并判断ek+1大小，如果ek+1<0.694，有μk+1=ak+1+exp(-1.9ek+1)(bk+1-ak+1)，否则有μk+1=ak+1+(1-exp(-ek+1))(bk+1-ak+1)，并计算μk+1处的弹道轨迹，转(5)。

(4)置ak+1=ak，bk+1=μk，μk+1=λk，yk+1=bk+1-ak+1，ek+1=yk+1-yk，并判断ek+1大小，如果ek+1<0.694，有λk+1=ak+1+(1-exp(-1.9ek+1))(bk+1-ak+1);否则，有λk+1=ak+1+exp(-ek+1)(bk+1-ak+1)，然后计算μk+1处的弹道轨迹，转(5)。

(5)置k=k+1，返回(2)。

2.3 算法分析与对比

为验证两种算法的优劣性，首先以单峰函数求解极值为例，比较两算法的收敛速度和精度，构建函数如下：

(13)

式中：x∈[-1,1]。

图3 指数优化算法流程图Fig.3 Flow chart of exponential optimization search

图4 边界计算流程示意图Fig.4 Flow chart of boundary calculation

图5 求解极值函数示意图Fig.5 Schematic diagram of extreme value

图5为所求极值函数示意图，对两种算法的初始空间都取为[-1,1],给定初始精度，以迭代过程中的空间误差和函数值误差以及计算时间来比较两种算法的优劣性。空间误差即为相邻两次迭代搜索空间的距离变化量，函数值误差为相邻两次搜索极值的变化量。

图6显示指数优化算法比黄金分割法空间误差下降快，特别是在迭代前期，曲线下降速率较大，搜索区间的变化较明显。图7显示两种算法的函数值误差变化曲线，两者在迭代前期函数值误差变化差别不大，几乎相同，后期下降速度逐渐加快。

图6 空间误差变化曲线Fig.6 The lose change of zone with iteration

综合图5、图6可知，虽然指数优化算法在误差和迭代次数上有一定的优势，但两者的性能差距并未充分展现出来。主要原因在于测试的函数比较简单，函数斜率变化单一，函数值域区间较小，待搜索空间长度不大，因此两种算法的误差收敛值相似，而且整体迭代次数都较小。利用相同的方法，增加两种不同的函数(式14和式15)进行比较，两种算法的优劣如表1所示。

图7 函数值误差曲线Fig.7 The loss change of function value

|x|)cosx+4ln(15-|x|)sinx

(14)

式中：x∈[-10,10]。

|x|)cosx+3x0.1ln(15-|x|)sinx

(15)

式中：x∈[-10,10]。

表1 算法误差对比Table 1 Comparison of algorithm error

表1中函数2和函数3的搜索空间均为[-10,10]，后两个函数的斜率变化较函数1更为复杂。此时指数优化分割算法表现出较大的优势，迭代次数将近缩短一半，计算时间代价较小，同时函数值的误差也更小。

3 Simulink仿真试验

3.1 仿真结束条件判定

当导弹击中目标或确定导弹脱靶时都需要停止仿真，根据弹道脱靶量、工作时间和导弹末端与目标相对速度等参数确定仿真终止条件，攻击失败的检测条件可以简化如下：

1)距离限制

导弹脱靶量超过规定最小的相对距离值，dMD≥dMDmax目标丢失。

导弹飞行的高度不能高于25 km，不能低于200 m。

2)时间限制

导弹的弹道飞行时间大于可控时间，tD≥tmax，导弹失控。

导弹的飞行时间不能大于能源可支持的最大时间。

导弹与目标遭遇时，引信解锁时间必须在安全飞行距离之外，若引信未及时打开，目标丢失。

3)速度限制

导弹惯性飞行速度小于导弹飞行的最小速度时，无法提供足够的升力。

导弹接近目标的相对速度太大或太小，引信不能正常工作。

3.2 仿真流程

根据比例导引方式和导弹工作原理，构建简化的制导系统控制回路，如图8所示，主要由导引头回路，控制指令回路，自动驾驶仪三部分组成，通过弹体运动学方程和相对运动方程形成闭环回路。

导引头持续对目标进行跟踪，使天线瞄准目标，产生与目标线旋转角速度成正比的控制指令信号，自动驾驶仪对导弹的运动参数做出调整，使导弹以既定的程序导引规律接近目标。

在该平台下，利用Simulink对弹道运动学方程、相对运动方程、目标运动方程进行仿真求解，计算导弹运动轨迹，利用黄金分割法搜索边界值，根据仿真判定条件确定是否终止计算还是继续下一次循环。

图8 比例导引规律制导回路示意图Fig.8 Circuit of missile guidance and control system

4 结果分析及对比

4.1 算法对比

以仿真计算攻击区包线远边界值为例，在搜索时分别使用黄金分割法和指数优化算法，对比分析两种算法在迭代次数、计算时间和精度误差上的优劣性。假设目标和载机高度均为8000 m，初始速度为270 m/s，导弹比例导引系数K=4，两种算法的对比结果如图9所示。

图9 两种算法误差对比Fig.9 The comparison of loss in simutation between two algorithms

两种算法在迭代初期并未表现出较大的差距，随着迭代次数的增加，指数优化算法误差下降速度更快。指数优化算法在迭代23步时开始收敛，黄金分割法则需要40步。黄金分割法的误差在90 m左右，指数优化算法误差可控制在20 m之内。现在空空近距格斗弹弹片杀伤范围可允许误差在10～20 m之内，因此黄金分割对命中概率会带来更大的误差。无法满足实战要求。

4.2 导弹轨迹及参数变化分析

目标在相对导弹高度为5000 m下做水平匀速运动，速度为270 m/s,图10(a)为导弹以比例导引法追踪目标的轨迹，比例系数K=4，图11(a)表明在目标不做任何机动情况下弹目距离近似成线性变化。图12(a)显示在后期导弹前置角变化量减小，弹道在末端比较平直。

假设目标与载机高度差为4000 m，目标做蛇形机动运动，速度为270 m/s，导弹比例导引系数K=4，图10(b)为弹道轨迹三维曲线，图11(b)表明导弹在前期弹目距离缩短较快，后期弹目距离缩短较平缓，目标在过载机动，导弹也需要一定的机动时间，图12(b)表明导弹前置角变化量一直在增大，制导末端导弹迅速调整姿态角追踪目标。

图10 导弹-目标三维轨迹示意图Fig.10 Missile-target 3D trajectory

图11 导弹-目标距离曲线示意图Fig.11 The change of distance between missile and target

图12 比例导引前置角变化曲线示意图Fig.12 The change of lead angle with proportional guidance

4.3 高度对攻击区的影响分析

目标和载机飞行高度同时在8 km、10 km、12 km，速度均为270 m/s。图13表明目标和载机在同一高度时，高度越高，空气密度减小，空气阻力对导弹的影响变小，导弹攻击区作用距离会增大。图14显示高度差为2 km和4 km时攻击区的范围主要受载机高度的影响，当载机高度比目标机高度大时，导弹俯冲目标，载机会形成一定的高度优势。