基于射流微分方程组的水炮命中问题求解*

张兴江，王 军

（南京理工大学先进发射协同创新中心，南京 210094）

0 引言

随着海上领土纠纷的升级和维权活动的增多，消防水炮作为一种非杀伤性武器得到了越来越多的关注和应用。消防水炮可以用来驱逐外来骚扰和自卫反击，同时又可以避免伤亡和引起战争，是海上出巡和海上捕鱼作业时理想的武器装备。当前，我国海洋维权局势日益严峻，维权任务十分艰巨，中国海监海上维权执法力量亟须加强。对此有分析人士认为，在中国现有海监船无武器配备情况下，高压水炮是巡航执法最有利的装备，可有效震慑违法船只[1]。

智能消防水炮在灭火系统中有着大量的应用与研究，但其较多地集中在射流方程的建立。目前命中问题求解主要有基于数学模型和基于计算机视觉两种方法。基于计算机视觉的方法[1]，主要是通过图像识别对射流落点和目标进行闭环控制，但其在调整速度上仍不及基于数学的方法，在此不进行分析。水炮中基于数学模型的计算方法的研究大多为针对静止目标、固定高度，文献[2]虽然针对不同高度的目标进行了研究，但其根本仍然为事先固定好高度，不能跟随目标高度变化实时修正。

对于运动目标的命中问题求解，在火炮控制中有着大量的研究[3]，但因为求解时间限制，火炮解算中使用更多的是射表法。而对于水炮来说，其解算周期长，精度明显低于火炮，所以可以使用求解射流微分方程组的方法来求解水炮的命中问题。

1 船用水炮射击解命中问题描述

解命中问题的实质就是确定射流与运动目标在空间相遇点即提前点的坐标。将目标航迹的数学模型与射流的数学模型联立求解，可以求得命中点坐标，进而换算出水炮的射击诸元。

图1 解命中问题描述

假设该时刻对应的命中点为Mq。已知经过滤波得到射击瞬时目标在大地坐标系下的位置坐标和速度，目标运动TF时间后到达命中点Mq，Mq点在水炮大地坐标系下的直角坐标为，计算得到Mq点在水炮大地坐标系下的极坐标为，且对应的水平距离和高度分别为dq1和hq1。目标在水炮大地坐标系下方位角β，高低角α。

1.1 目标航迹数学模型

考虑到船和目标都是慢速目标，在一个解算周期内往往没有明显的机动，所以解算时将目标与船的运动看作匀速运动。因此，每个解算周期内目标航迹的数学模型便是以目标相对于水炮的坐标为起点，以目标相对于水炮的速度做匀速直线运动，即目标在水炮大地坐标系内做匀速直线运动。航迹的数学模型如下：

其中，Vs、Xs分别为目标在水炮大地坐标系的速度、位置矢量，d为目标距炮口的水平距离，hs为目标距炮口的垂直距离。

1.2 射流数学模型

目前大多数水射流运动模型都基于牛顿第二定律建立。无论是考虑水射流中的一段射流微团或是考虑射流中的水滴，一般都只考虑空气阻力和重力[2]。针对各种情况下的射流方程也有着许多的研究，但对于水炮不同型号、不同流量目前并没有一种通用的数学模型，常用的方法是利用射流运动学公式进行推导，通过射流实测数据确定射流方程中的待定系数，从而获得适用于此种水炮的射流模型。这种方法虽然通用性差，但针对特定的水炮也可以获得较高的精度。

显然射流的运动可以分解为垂直面和水平面两种运动，二者互不影响，求解射流方程时，为便于分析，将水炮大地坐标系绕Zp轴旋转得到射流坐标系Xh-Yh-Zh，使射流沿Yh-Zh平面射出，则射流垂直面内的运动将完全在Yh-Zh平面内，水平面内的运动将完全在Xh-Yh平面内，即Yh轴反映射流的水平距离，Zh轴反映射流的垂直距离，Xh轴反映射流的横向偏移量，如图2所示。

图2 射流在射流坐标系的运动

射流运动模型的不同关键在于空气阻力的选取不同，不同的空气阻力公式可以构造不同的射流方程，其能达到的精度也各有不同。HATTON等[4-5]提出了一种空气阻力模型，证明其可以模拟大流量消防炮的射流轨迹，而船用水炮正是一种大流量消防炮，故本文用此种模型建立射流方程，其空气阻力模型为：

式中，S为射流轨迹长，满足

其中，Fyz为Yh-Zh平面内单位质量的空气阻力，Fx为Xh-Yh平面内单位质量的空气阻力，待定系数k、b、A、B可以参考文献[4-6]由实验数据获得。

按运动学规律，射流的数学模型最终建立如下：

2 命中方程求解

命中方程的求解即联立射流微分方程组式（5）和目标航迹数学模型式（1）进行数值求解，可以同时获得射流与目标相遇时间tf和水炮射击诸元。其中求解微分方程组的方法很多，如欧拉矩形法、梯形法、龙塔库塔法等，在此，选用常用的四阶龙塔库塔法进行数值求解[7-8]。

2.1 初始条件的求取

要求解射流微分方程组，初始条件至关重要。在水炮的射流命中问题中，已知条件有：

这是典型的两点边值问题，如何从边值条件（6）出发，求解出水炮射击诸元和射流飞行时间tf，正是水炮命中问题求解的任务。在此，以射流落点为依据，使用变步长的四阶龙塔库塔法，将边值问题变为初值问题迭代求解。

2.1.1 初始高低角αt和轨迹长S0

轨迹长S0=0；高低角αt初值选取为目标坐标对应的高低角，即

船用水炮在射击时往往采用固定的流量和压力，因此，射流相对炮口初速vp0便也确定[9]，后面也不考虑vp0变化的情况。射流初速主要影响因素有：vp0、风速、船的运动、水炮的转动。对于船的姿态角在解算周期内变化很小，造成的影响可以忽略。

图3 合成速度示意图

其中，vw为风速大小，βw为风向，vb为船速大小，βw为船航向，βq为目标点相对于水炮的方位角。由于风的垂直分量比较小，因此，风只在水平面上分解为横风和纵风，同时船基本没有垂直方向的运动。

假设水炮在每个解算周期内的运动为匀速转动，忽略加减速过程，那么由于水炮转动产生的射流初速度便不能忽略。考虑水炮在持续射击的过程中，每个运动周期水炮高低转动的角度很小，所以水炮转动对vy、vz产生的影响可以忽略。

解算周期开始获得的水炮方位角β0，对方位角作线性预测，假设解算周期结束对应架位β，可得水炮的方位速度前馈为

所以射流在炮口处的线速度：

其中，T为解算周期，Vβmax为水炮方位转动最大角速度，L为水炮身管长。

所以射流初速为

2.2 求解流程

2）使用四阶龙塔库塔法求解射流微分方程组式（5），设定初始步长为0.01 s。每次求解过程中有关射流诸元序列值：

i为以步长h的计算序列号。求解终止条件为：

当求解终止时，此时i对应的序列值则为射流在射流坐标系的落点坐标和速度，飞行时间：

3）利用求得的tf和目标航迹模型式（1）获得目标未来点坐标，误差分别为：

4）对高低角修正，新的射流高低角为

αt为上一次龙塔库塔法求解的高低角初值。

利用新的目标未来点和修正量确定初值。

3 实验仿真分析

针对上述数学模型和命中求解流程，使用matlab2014进行实验仿真，验证求解方法的可行性，并进行相关的分析。

仿真时相关参数如下：

各空气阻力系数：k=0.000 8，b=0.023，A=20，B=0.04，g=9.8（在此省略单位），水炮部分参数：L=0.5 m，vp0=50 m/s，解算周期 T=0.2 s，解算精度0.1 m。

因为解算过程纵风和船速最终会分解到对应轴，本质相当于船静止。所以设定固定横风vwx=2m/s，船速vb=0，纵风vyz=0。

水炮高低角为30°时，射流轨迹如图4所示。

图4 射流轨迹

3.1 静止目标、匀速运动和匀加速运动目标的对比

各种目标仿真参数如下：

静止目标航迹：位置（10，50，-3），匀速运动目标航迹：起始点（10，50，-3），速度（3，2，0.1），匀加速运动目标航迹：起始点（10，50，-3），速度（3，2，0.1），加速度（1，-1，0）。

3.1.1 单点的解算误差

对于解算点来说，解算误差和射击诸元仿真如图5～图7所示。

图5 静止目标解算点误差和射击诸元

图6 匀速运动目标解算点误差和射击诸元

图7 匀加速目标解算点误差和射击诸元

从以上仿真结果可以看出，因设定迭代误差为0.1 m，假定目标求解过程做匀速运动，所以当目标实际静止或做匀速运动时，解算误差小于0.1 m；而当目标做匀加速运动时，可以看到解算误差大于0.1 m，且具有随时间误差增大的趋势，这是由于目标航迹模型假定不一致造成的，增大的误差主要是加速度产生的位移，符合公式，增大趋势是因为随着目标越来越远，射流飞行时间越来越长，加速度产生的影响也越来越大。

3.1.2 整个解算周期内匀速转动的解算误差

对于解算周期内匀速转动水炮来说，以目标做匀速运动为例，解算误差如图8所示。

图8 匀速运动目标解算周期误差和射击诸元

与解算点误差进行对比，解算周期内实际射流落点和目标误差变化趋势与解算点误差变化趋势相同，只是误差变大。原因在于水炮在每个解算周期内匀速转动，因水炮高低角不变，所以射流的水平距离不变，随着目标的运动，误差便会越来越大，到达下一个周期，误差便会重新变小，然后逐渐变大，所以呈现的结果就是锯齿状误差。而对于静止目标来说，因没有转动，所以误差不变。

3.2 不同解算周期仿真对比

以目标匀速运动为例，不同的解算周期对于解算点的误差影响不大，都保持在0.1 m以内，而对于每个解算周期内匀速转动的情况有明显影响，仿真结果如图9所示。

图9 匀速运动目标不同解算周期误差

其中，按从左往右、从上往下顺序，解算周期分别为：0.1 s、0.2 s、0.5 s、1 s。

从图中可以看出，随着解算周期的增大，解算误差越来越大，原因在于解算周期越大，周期内水炮匀速转动的时间越长，积累的误差便会越大。同理，可以知道对于固定解算周期，目标速度越快，误差越大。

3.3 不同解算精度的对比

求解过程中迭代误差设定直接影响解算精度，当设定0.1 m时，解算点精度在0.1 m内，当设定0.5 m时，解算点解算精度在0.5 m内，而在整个解算周期内，解算误差仿真结果如图10所示。

图10 解算精度0.5时解算误差

与图8对比，迭代误差0.1 m和0.5 m时，解算误差变化不大。所以可以根据解算误差适当的选择迭代误差，以减少迭代次数和时间。

3.4 解算时间

按目标静止，设定距离，解算时间如表1所示。

表1 不同距离解算时间

从表1可以看出，目标距离越远，解算时间越长。原因在于步长和高低角的调整不理想。但总体来说，选择适当的解算周期和距离，解算时间是可以满足要求的。同时降低解算精度可以明显减少解算时间。

4 船姿态角变化的影响

船在海上行驶时，势必会因为海面的运动造成船姿态角的变化，相应的水炮姿态角也会发生变化。由前面所述，通过命中方程求解可以求得水炮在水炮大地坐标系下的射击诸元，但因为水炮的姿态运动，水炮最终输出的射击诸元不能不发生变化。如果水炮驱动控制系统自带有处理水炮姿态变化的程序，那么只需输出水炮在大地坐标系下的射击诸元即可；但如果其没有，在命中问题求解后输出的射击诸元必须进行相应的处理，以抵消水炮姿态角变化带来的影响。

通过船体测量系统和相应的姿态滤波等处理程序，可以获得水炮射击时的水炮姿态角，然后只需将水炮的射击诸元由水炮大地坐标系进行姿态坐标变换，转换到水炮的载体坐标系即可[3]，此时水炮在水炮载体坐标系下的射击诸元即为输出到水炮驱动控制系统的射击诸元。

5 结论

提出了基于射流微分方程组水炮的实时命中问题求解方法，介绍了目标航迹数学模型和射流数学模型的建立，对于影响射流的因素：风、船的运动、水炮转动等统一处理，使用四阶龙塔库塔法求解微分方程组，并通过修改步长和高低角不断迭代，直至满足精度要求，获得水炮射击诸元。通过matlab仿真验证，其解算误差在0.1 m内，证明了方法的可行性，为水炮自动控制系统的开发提供了一种可行思路。

通过实验对不同目标航迹、不同解算周期和解算精度等情况的仿真，得出了不同情况下的解算误差，为水炮的实际应用，如数学模型的建立、控制策略的选择、解算精度和解算时间的合理优化等，提供了一定的参考价值。

但在命中求解过程中仍有许多可以改进的地方[10-14]，比如龙塔库塔法步长的选择、高低角的修正、初值的确定等，并且射流微分方程可以进行简化，上述的微分方程组具有7个变量，可以将其减少或分离，以减少解算的时间。同时在每个解算周期内水炮匀速转动的控制策略会造成一定误差，可以尝试其他控制策略以优化误差。