库塔默森法的解命中方法进行仿真分析，总结了诸元解算过程中积分造成的误差与积分误差估计值之间的统计经验，基于该经验提出了积分实际误差与估计值的换算关系，基于该换算关系应用有最优步长控制策略的龙格库塔默森法进行C-RAM解命中计算。将该方法与采用4阶龙格库塔法以及龙格库塔夏普法的C-RAM解命中方法进行对比验证。仿真结果表明：通过最优步长控制，该方法可在不同斜距离下保证解算精度，实现对解算时间的最优控制，且默森法也是最适合用于C-RAM解命中计算中变步长控制的方法。摘要：为解决传统高炮火控系统无法实现C-RAM解命中计算的问题，针对RAM目标的运动模式，提出了一种基于变步长龙格-库塔法的高炮C-RAM火控解命中方法。通过对应用定步长龙格库塔默森法的解命中方法进行仿真分析，总结了诸元解算过程中积分造成的误差与积分误差估计值之间的统计经验，基于该经验提出了积分实际误差与估计值的换算关系，基于该换算关系应用有最优步长控制策略的龙格库塔默森法进行C-RAM解命中计算。将该方法与采用4阶龙格库塔法以及龙格库塔夏普法的C-RAM解命中方法进行对比验证。仿真结果表明：通过最优步长控制，该方法可在不同斜距离下保证解算精度，实现对解算时间的最优控制，且默森法也是最适合用于C-RAM解命中计算中变步长控制的方法。
C-RAM(counter-rocket artillery mortar)是一种使用我方防空火力主动拦截来袭的火箭弹、榴弹和迫击炮弹，以达到防御敌方地面曲射火力打击的主动防护系统[1]。高炮C-RAM火控系统的任务是根据雷达对来袭RAM目标的跟踪数据估计出其实时位置、速度并辨识其弹道系数，根据以上信息通过解命中计算得出目标的未来命中点并射击，从而实现对RAM目标的拦截。但在火控处理时，由于RAM目标的运动模式的非线性性[2]，传统的高炮解命中方法因假定模型以及计算方式的问题无法得到正确的解命中结果。由此需要设计一种基于RAM目标运动方程假定，通过积分方法计算未来点的高炮C-RAM解命中方法，以实现高炮C-RAM作战的需要。
为解决相关问题，文献[3]提出了以质点弹道方程为基础，采用龙格库塔法积分弹道方程，解算RAM目标未来点的解命中方式。但由于此方法采用固定步长进行积分，在不同情况下积分误差存在一定差距，导致解算精度会根据目标斜距离变化。在不同斜距离下，定步长解命中方法的积分次数会根据弹飞时间的增长而增加，多次的积分会累积更大的积分误差，导致随着目标斜距离的增加解算误差会成倍增大。
本文就高炮C-RAM火控解命中计算的问题，在使用定步长龙格库塔法的C-RAM解命中方法的基础上，总结了使用定步长的龙格库塔默森法进行解命中时，诸元解算实际积分误差与积分误差估计值之间的统计经验。以此统计经验为基础提出了两者的换算关系，基于该换算关系应用最优步长控制策略[4]的龙格库塔默森法进行C-RAM解命中计算的外弹道外推。仿真证明该方法相比于使用定步长法及龙格库塔夏普法的C-RAM解命中方法，能够保证不同斜距离上的解命中精度，控制积分步长为误差允许范围内的最大值，减少积分次数，降低解命中计算时间。
由于在实际进行C-RAM作战时目标弹丸转速、攻角、弹翼升力等状态无法测得，这些状态对于弹丸运动轨迹在短时间内的影响较小，在进行短时间的外弹道外推时相对于弹道系数估计误差造成的未来点偏移其影响基本可以忽略不计，故本文采用弹丸质心运动模型作为外推RAM目标外弹道轨迹的公式，将除去空气阻力与重力外的其他力对弹丸运动的影响看作是弹道系数误差所造成的影响的一部分或是将其忽略。基于C-RAM作战时一般在RAM目标的外弹道末端进行拦截，此时目标速度通常处于亚音速状态，弹道系数值的变化幅度较小，故可以将弹道系数视为常值进行处理，并假设气象条件是标准的，无风雨干扰。弹丸质心运动模型示意图如图1。
图1 弹丸质心运动模型示意图
Fig.1 Schematicdiagram of projectile centroid motion model
弹丸质心运动方程组如式(1)所示:
(1)
式中: x、y、z、vx、vy、vz分别为弹丸的位置坐标与速度分量;c为弹道系数；g为重力加速度；cs为声速；H(z)为空气密度特征函数；G(v,cs)为空气阻力函数[5-7]。
库塔默森法由于RAM目标运动模式的非线性性，应用代数方法外推其外弹道会产生较大的误差[8]，故应使用积分法外推外弹道[9]，为保证诸元解算的精度，本文使用了龙格库塔默森法这一数值积分方法外推目标外弹道。
龙格库塔默森法是一种可估计积分误差的龙格库塔积分法，其计算公式为式(2)所示的四阶五级公式，其积分结果为yk+1，另外采用了一个三阶四级公式求了一个低阶积分结果
利用
来估计积分误差[10]。此误差估计值实质为积分过程中的截断误差，只能作为实际误差的参考值使用，如需得到较为准确地误差估计值，则还需对此值做进一步处理。
龙格库塔默森法计算公式为式(2)所示[10-11]:
(2)
其中：
(3)
误差估计值计算公式如式(4)所示：
(4)
此误差估计值的精度为3阶，而此方法的积分精度为4阶，对于大部分使用情况已经能够满足精度要求。
在命中方程求解过程中，由于目标外弹道需要使用龙格库塔法进行外推，计算量较大。为了尽可能减少积分次数，本方法在普通的命中方程迭代求解方法[12-13]的基础上，做出了一些针对性的改进。结合龙格库塔积分法的C-RAM命中方程求解方法的具体步骤如下：
1) 初值选取。迭代的初始值的选取直接影响收敛迭代次数，且还需要确保初始值不超过实际命中点过多，因此首次解命中要优化初始值选取，弹飞时间初值可取
其中Tf(db,zb)为弹飞时间计算公式，C为一个小于1大于0的实数，本文中此处取C=0.8。否则，可取上次解命中的弹飞时间的0.9倍为本次迭代的初值
同时，设置i=1。
2) 命中点坐标求解。目标外弹道在使用龙格库塔法进行外推时，若积分步长过大会导致截断误差较大从而影响诸元解算精度[14-15]，因此外推时需要进行多步积分完成一次命中点坐标的求解，为了减少总的积分步数，需要在一次命中点求解时舍去最后一段小于积分步长的时间。命中点迭代坐标数学表达示如式(5)所示：
(5)
其中各命中点迭代坐标
为目标实际位置经过j次积分的结果，并如式(6)所示，
为弹飞时间的迭代值，hk为第k次积分时的积分步长，
为舍去最后一段小于积分步长的j次积分的总时间。
(6)
同时，若
则
此时只需要积分一次，故也不再需要对外推时间长度进行控制。
以上命中点积分方法使得在进行弹道外推时预测命中点不会超过实际命中点过多，从而最大程度上的减少积分次数。
3) 命中点水平距离和高度求解:
(7)
(8)
其中:
与
分别为命中点水平距离与高度的第i次迭代计算值；
为水平距离修正；
为高度修正。
4) 弹丸飞行时间求解:
(9)
5) 循环判断。若
为命中时间允许误差，则令i=i+1，转到第2步循环，否则结束循环，得到命中点空间位置及弹丸命中时间
龙格库塔积分法的C-RAM命中方程求解方法解算过程如图2所示。
图2 结合龙格库塔积分方法的迭代法收敛曲线
Fig.2 Convergence diagram of iterative method combined with Runge-Kutta integral method
库塔默森法下积分误差估计值及诸元解算误差的统计经验基于式(4)计算出的的库塔默森法的积分误差估计值并不准确，若是要在保证命中精度的情况下对积分步长进行最优控制，就必须要找到积分误差估计值与实际解算误差值之间的统计经验。
本文采用以0.001 s作为积分步长的4阶龙格库塔法生成原始弹道作为误差分析时的命中诸元标准值，并以本文中的迭代法以不同积分步长对此弹道进行诸元解算并计算解算误差。对比不同积分步长下的实际解算误差与积分误差估计值总结出了其之间的统计经验。
总结出以下结论：
1) 解命中误差与积分步长成正比；
2) 对于拥有同一初速度、速度方向、高度、弹道系数但高炮与弹丸轨迹水平投影斜距不同的弹道，解算误差除以弹丸飞行时间的值与积分步长为线性关系，且斜率相对一致；
3) 不同斜距的外弹道轨迹下，x轴与y轴位置坐标估计误差与积分步长的3次方的比值与积分步长为斜率相对一致的线性关系；
4) 不同斜距的外弹道轨迹下，z轴位置坐标估计误差与积分步长的2次方的比值与积分步长为斜率相对一致的线性关系；
5) 同一条弹道在进行误差估计时在在积分步长一致的情况下估计误差的变化较小。
各项结论所总结的关系如图3—图6所示。
图3 x轴解算误差与弹飞时间的比值曲线
Fig.3 Ratio curve between x-axis solution error and missile flying time
图4 x轴估计误差与步长3次方的比值曲线
Fig.4 Ratio curve of x-axis estimation error to the third power of step size
图5 z轴估计误差除以步长2次方值的曲线
Fig.5 Ratio curve of z--axis estimation error to the quadratic of step size
图6 同一弹道在定步长下不同次积分的估计误差变化
Fig.6 Estimation error variation of different times of integration for the same trajectory under fixed step size
对不同弹道进行仿真后确认了以上统计经验具有广泛性，故可以将实际误差、误差估计值、积分步长之间的关系总结为如式(10)所示的公式。
(10)
式(10)中： Cx、Cy、Cz为各轴的误差换算比值常数；Ex、Ey、Ez为各轴位置估计误差；ht为积分步长；Ex_ture、Ey_ture、Ez_ture为实际位置解算误差； tf为弹丸飞行时间。
在对多种弹丸的弹道进行大量仿真验证后发现，在不同弹道系数及弹丸速度下，误差换算比值常数Cx、Cy、Cz是不同的。对一飞行轨迹起点固定的弹丸进行定步长解命中仿真，设定弹道系数为0.5，改变每次计算时x轴的初始速度，x轴误差换算比值常数随积分步长与x轴速度变化曲面如图7。积分步长为0.05 s时，x轴误差换算比值常数在不同弹道系数下随x轴速度变化的曲面如图8。
图7 不同速度、积分步长下的误差换算比值
常数变化曲面
Fig.7 Variation curve of error conversion ratio constant under
different speeds and integration steps
图8 不同速度、弹道系数下的误差换算比值
常数变化曲面
Fig.8 Variation curve of error conversion ratio constant under
different speeds and ballistic coefficient
从图7、图8的结果中可以看出该比值常数在不同速度以及弹道系数下是不同的，变化规律也呈非线性，但不随积分步长变化。在使用式(10)估计实际解算误差时，比值常数的取值将决定估计出的实际解算误差准确与否。
对于比值常数的取值问题，在实际应用时的解决方式有3种：
1) 通过大量的实验来制作对应不同弹道系数、速度的比值常数表，在计算时根据目标的弹道系数以及速度从表中取值；
2) 将比值变化曲线近似线性化，在计算时根据近似函数得出比值用于计算；
3) 取一固定值，保证大部分高概率情况的解算效果，牺牲小部分低概率情况的解算精度和速度。
此外基于库塔默森法的统计经验，还对库塔夏普法进行了统计经验的总结，并计算出了夏普法误差换算比值常数变化曲线。
龙格库塔夏普法计算公式为式(11)所示[10-11]:
(11)
其中：
(12)
其误差估计值计算公式如式(13)所示：
(13)
夏普法中K5正好是下一次积分计算时的K1，因此只需要在第一步多算一次K1，之后每步积分只需要计算4次右端函数。其误差换算公式如式(14)所示，形式与默森法基本相同，但积分步长次方值不同。
(14)
根据式(14)计算出的夏普法误差换算比值常数变化曲面如图9、图10。可以看出相比于默森法的比值常数变化曲线，夏普法的比值常数在速度和弹道系数较大时，同一速度、弹道系数下其不同积分步长对应比值常数变化较大。故在速度和弹道系数较大的情况下，若要进行最优步长控制，夏普法的步长控制难度要大于默森法，故相比于默森法，夏普法并不适用于C-RAM的解命中计算。
图9 弹道系数=0.2，不同速度、积分步长下
夏普法的误差换算比值常数变化曲面
Fig.9 Variation curve of error conversion ratio constant
under different speeds and integration steps
图10 弹道系数=0.5，不同速度、积分步长下
夏普法的误差换算比值常数变化曲面
Fig.10 Variation curve of error conversion ratio constant
under different speeds and integration steps
库塔默森法以式(10)所描述的估计误差与实际误差之间的关系为基础，便可以实现应用最优步长控制的库塔默森积分法进行解命中计算，此处认为使解命中误差保持在可接受范围内的最大积分步长为最优步长。对于比值常数的取值问题，本文采用了取一固定值作为比值常数的方法。
具体的最优步长控制策略为：
① 首先以一较为合适的固定步长hm进行第一次解命中，求得弹丸飞行时间tf。
② 给定误差限制ε0，以固定步长hm作为第一次积分步长，以上一次解命中的弹飞时间tf作为本次解算参考弹飞时间，求若以本次积分步长进行解命中时各轴诸元解算的误差值εx、εy、εz，公式为式(15)所示。
(15)
若εx≤ε0且εy≤ε0且εz≤ε0，则本次积分成功，通过式(16)确定下一步步长hm+1。
(16)
若εx>ε0或εy>ε0或εz>ε0，则本次积分误差过大，同样通过式(16)求出一个积分步长hm+1。为了尽量不重复进行积分运算，减小计算时间，在一次诸元解算中除第一次积分外均无需重新积分。在同一次解算中同一积分步长造成的误差相差不大，除第一次积分外若出现计算误差值大于误差限制值的情况，只需要减小下一次积分步长保证下一次积分的精度即可。若对解算精度的要求较高，则只需要将误差限制值设置的比理想误差值稍小就能控制每步积分的误差在范围内。
③设置最小步长hmin，一旦积分步长小于此步长值则不再适用步长控制，此举是为了避免计算机在步长过小时由于舍入误差造成数值偏差比过大，从而计算出一个过大的误差值，引起步长振荡导致解算失败。
为验证基于统计特征的最优步长控制策略是否能够有效控制积分步长使解命中误差保持在可接受范围，对某2个RAM目标进行解命中计算，假设目标一被探测到时的弹道起点为x=500 m、y=1 000 m、z=1 000 m,高炮位于坐标原点，弹丸速度向量为x轴方向-250 m/s、y轴方向-250 m/s、z轴方向-90 m/s，弹道系数为0.2。
目标二被探测到时的弹道起点为x=500 m、y=1 000 m、z=1 500 m,高炮位于坐标原点，弹丸速度向量为x轴方向-200 m/s、y轴方向-150 m/s、z轴方向-150 m/s，弹道系数为0.3。
假设解命中过程中所得到的目标运动状态以及弹道系数均为理想值，每次诸元解算间隔为0.01 s，共解算600次。基于变步长默森法的解命中法为方法一，基于变步长夏普法的解命中法为方法2，文献[2]中的定步长解命中法为方法3，变步长法设置误差限制为各轴位置误差小于0.1 m，方法1比值常数取值为Cx=Cy=14 000、Cz=110；方法2比值常数取值为Cx=Cy=2 200、Cz=20；定步长法设置积分步长为0.05 s。
分别使用3种解命中方法对目标一及目标二进行仿真。诸元解算后进行逆解得到的解算误差如图11—图13所示。
图11 x轴位置解算误差图
Fig.11 x-axis position calculation error
图12 y轴位置解算误差图
Fig.12 y-axis position calculation error
图13 z轴位置解算误差图
Fig.13 z-axis position calculation error
从仿真结果中可以看出，相对于定步长的C-RAM解命中的方法，变步长法能够使各轴方向上的解算误差基本保持在设定的误差范围内。但由于比值常数使用了取固定值的方法，导致对误差的控制精度较差，解算误差超出设置范围的情况，若是采取使用比值常数表或者拟合比值变化曲线的方法，将能进一步的提升误差控制的精度。
在变步长法中，相比默森法的误差控制效果，夏普法在对目标一解命中时效果与默森法接近，对目标二解命中时误差小于默森法。但误差较小并不意味着控制效果好，夏普法与默森法的积分精度同为四阶且误差估计精度也同为三阶，在相同的精度条件下，更小的解算误差代表消耗了更多的解算时间，证明了在积分步长控制效果上夏普法要弱于默森法。且仿真中发现对夏普法应用本文中的最优步长控制时不能将上一次计算出的K5代入下一次计算中，否则将会因为步长的变化导致误差估计值的计算出现错误，进而无法有效控制积分步长，故在每步计算量上夏普法与默森法也基本相同，使其不具备计算量的优势。
变步长默森法积分步长的控制效果如图14、图15所示。
图14 变步长法解算目标一时积分步长控制云图
Fig.14 The time integral step control of using variable step
method to solve the target one
图15 变步长法解算目标二时积分步长控制云图
Fig.15 The time integral step control of using variable step
method to solve the target two
由图14、图15所展示的积分步长控制情况中可以看出，本文中的变步长C-RAM解命中方法能够实现对积分步长的最优控制。其能够实时的调节积分步长的大小，在积分误差较大的情况下减小步长保证精度，积分误差较小时增加步长减少积分次数。实现在保证解命中精度的情况下尽可能增加积分步长的长度，减少计算时间。
本文提出了一种使用最优步长控制的龙格库塔默森法的高炮C-RAM解命中方法，并将其与使用定步长龙格库塔法以及应用同样步长控制方法的龙格库塔夏普法的高炮C-RAM解命中方法进行了对比仿真。仿真结果证明：
1) 在基于弹丸质心运动模型模拟理想RAM目标外弹道轨迹的情况下，基于变步长龙格库塔法的高炮C-RAM解命中方法可以完成高炮C-RAM解命中任务。
2) 基于C-RAM解命中误差与积分估计误差规律应用的最优步长控制策略能够完成对积分步长的最优控制，相比使用定步长积分法进行C-RAM解命中计算，该方法可控制目标在不同距离、速度、弹道系数情况下的诸元解算误差。
3) 在同样应用最优步长控制方法同样精度的变步长积分法中，默森法相比于夏普法步长控制效果更好，更适合用于C-RAM解命中计算。
若能够在工程应用方面解决误差换算比值常数表的编制问题，则该方法便可以工程化应用于高炮防空系统中，可提高高炮的C-RAM作战能力。
[1] 毛征,孟灿,杨俊强,等.近程防空、反导与C-RAM武器系统现状与发展[C]//系统仿真技术及其应用学术论文集(第15卷).中国科学技术大学出版社,2014:25-34.
Mao Z,Meng C,Yang J Q,et al.Development of short-range air-defense,anti-missile and C-RAM system[C]//Collection of Academic Papers on System Simulation Technology and its Application.China University of Science and Technology Press,2014:25-34.
[2] 李文才.一种弹道目标的轨迹预测方法[J].四川兵工学报，2012,33(03):25-27.
Li W C.A trajectory prediction method of ballistic target[J].Journal of Sichuan Ordnance Industry,2012,33(03):25-27.
[3] 刘剑威.基于外弹道算法的C-RAM火控处理技术[J].指挥控制与仿真，2012,34(03):47-49.
Liu J W.Method for C-RAM FCS data processing based on ballistic equation[J].Command Control & Simulation,2012,34(03):47-49.
[4] 张华军,谢呈茜,苏义鑫,等.船舶操纵运动仿真中改进变步长龙格库塔算法[J].华中科技大学学报(自然科学版)，2017,45(07):122-126.
Zhang H J,Xie C X,Su Y X,et al.Improved variable stepsize Runge-Kuttaalgorithm in the ship manoeuvringmotion simulation[J].J.Huazhong Univ.of Sci.& Tech.(Natural Science Edition),2017 45(07):122-126.
[5] 李强,欧阳攀.基于外弹道的高炮火控算法与仿真[J].火力与指挥控制，2014,39(05):131-134.
Li Q,Ouyang P.Anti-aircraft gun fire control algorithm and simulation based on ballistic[J].Fire Control & Command Control,2014,39(05):131-134.
[6] 韩子鹏.弹箭外弹道学[M].北京：北京理工大学出版社，2014:76-78.
Han Z P.Exterior ballistics of projectiles and rockets[M].Beijing：Beijing Institute of Technology Press,2014:76-78.
[7] 王海宁,刘序旻,张瑶瑶,等.防空高炮对RAM弹类目标跟踪射击可行性分析[J].兵器装备工程学报,2017,38(07):94-97.
Wang H N,Liu X M,Zhang Y Y,et al.Tracking feasibility analysis of antiaircraft gun against the projectile of RAM[J].Journal of Ordnance and Equipment Engineering,2017,38(07):94-97.
[8] 李智，吴盘龙.基于外弹道的防空火控改进算法及模型仿真[J].指挥控制与仿真，2021,43(03):18-24.
Li Z,Wu P L.Improved algorithm and model simulation for anti-aircraft fire control based on external ballistics[J].Command Control & Simulation,2021,43(03):18-24.
[9] 徐亚军，王钰，宋贤龙.基于四阶龙格-库塔公式的火控系统瞄准角修正模型[J].火力与指挥控制，2020,45(05):111-124.
Xu Y J,Wang Y,Song X L.The aimed angle modified model fir fire control system based on fourth order Runge-Kuttaformula[J].Fire Control & Command Control,2020,45(05):111-124.
[10] 黄柯棣.系统仿真技术[M].长沙:国防科技大学出版社，1998:102-106.
Huang K D.System simulation technology[M].Changsha:National Defense University Press,1998:102-106.
[11] 肖田元.系统仿真导论[M].北京:清华大学出版社，2000:38-40.
Xiao T Y.Introduction to system simulation[M].Beijing:Tsinghua University Press,2000:38-40.
[12] 陈群斋.快速解命中的一种有效方法[J].火力与指挥控制，1994,20(01):20-24.
Cheng Q Z.An efficient method of quick solving hit problem[J].Fire Control & Command Control,1994,20(01):20-24.
[13] 徐国亮，王勇.舰炮反导火控原理[M].北京:北京理工大学出版社，2018:130-134.
Xiu G L,Wang Y.Anti missilefire control principle of naval gun[M].Beijing:Beijing University of Technology Press,2018:130-134.
[14] 冯勇，黄燕，黄兵.拟合弹道系数法进行弹道外推的研究[J].弹箭与制导学报，2014,34(02):117-119.
Feng Y,Huang Y,Huang B.Trajectory extrapolation based on fitted ballistic coefficient[J].Journal of Projectiles,Rockets,Missiles and Guidance,2014,34(02):117-119.
[15] 王干，熊风，欧能杰，等.基于反向扩展卡尔曼滤波的弹道外推算法[J].电光与控制，2020,27(12):49-52,89.
Wang G,Xiong F,Ou N J,et al.A ballistic extrapolation algorithm based on inverse extended kalmanfilter[J].Electronics Optics & Control,2020,27(12):49-52,89.
Citation format:LIU Xinyu, SHU Lipeng, LIN Zhiwei, et al.Hit solving method of antiaircraft gun C-RAM based on variable step size Runge-Kutta method[J].Journal of Ordnance Equipment Engineering，2022，43(11)：297-304.