基于射击误差状态方程的高炮毁伤概率建模

王向民，王 军，谢杰涛，郭 治

(1.南京理工大学 自动化学院，南京 210094;2.中国人民解放军 32200 部队，辽宁 锦州 121000)

为应对来自空中的威胁，防空反导系统作为最直接有效的反制手段，受到了各国的高度重视.以高炮为主的近程末端防空武器系统，其具有较高的费效比，常作为防御体系的最后屏障，因而它对来袭目标的毁伤关系到整个防空作战的成败.为获取高炮武器系统的毁伤概率：一种方法是通过大量的实弹射击，统计出射弹数目和实际毁伤目标的数量，二者比值即为武器的毁伤概率,该方法可信度高，但检验费用昂贵，且不能适用于武器系统的设计与论证阶段；另一种方法是建立基于高炮射击误差的毁伤概率数学模型，利用数值计算的方法获得高炮武器系统的毁伤概率，该方法通过对高炮的射击误差进行适当分解，利用误差模型转换法构建毁伤概率的计算模型.在我国现行的国家军用标准(简称国军标)中，把误差转换为重复误差和不可重复误差，虽然降低了计算复杂度，但存在较大的模型近似[1-4].射击误差中的弱相关误差源对于毁伤概率计算具有显著影响，陶德敬等[5]提出将射击误差分解为共有分量和非共有分量，运用误差序列预测理论给出了确定强相关误差和弱相关误差下射击过程的平均毁伤概率表达式，但该模型只能满足误差序列为2阶状态方程以下的情况；姚志军等[6]利用弱相关射击误差状态方程给出了递推形式的毁伤概率计算方法，但同样对于高阶自相关误差序列的相关性计算没有明确给出解析式.随着现代高炮武器系统射频的提高，射击过程中前一发弹药的毁伤概率对其后各发都有影响，因而每次发射的弹药对毁伤概率的贡献都是有差别的，而将误差序列简单地近似为1阶或2阶系统的模型简化方法[7]，有悖于高炮射击的真实模型.高炮的误差序列从本质上来说，可以视为武器身管的随机振动、跟踪误差以及火控解算误差等综合而成的平稳时间序列，有关时间序列相关性的分析方法，包括线性和非线性回归方法进行预测[8]、基于数据驱动的时间序列预测等[9-10]，需要有足够的实验数据才能获得较高的预测精度[11-12]，并且随着时间序列相关性阶数的增加，建模的复杂性和计算量也会大幅增加[13].

为了建立更精确的毁伤概率计算模型，本文以射击误差序列所构建的状态方程为出发点，运用卡尔曼滤波的思想，以递推方式给出连续脱靶条件下射击诸元误差的密度函数，将包括由射击冲击载荷影响在内的射击误差的弱相关分量分解为完全可预测与完全不可预测的分量，再分别与完全可预测的强相关分量和完全不可预测的不相关分量合并，构成2个随机分量，根据高炮进行点射时的时空特性，构建射击误差的随机状态模型，进而推导出相应的毁伤概率递推计算模型.

1 射击误差序列状态方程构建

x(k)=xb(k)+xr(k)+xq(k)

(1)

式中：xr(k)与xq(k)之和称为射击诸元误差xs(k).射击误差序列的状态方程可进一步表示为

(2)

作为随机状态方程，还要已知它的初始状态方差才能是完备的.设它的初始状态方差为

(3)

考虑到射击误差为平稳高斯过程，因而可得[5]

(4)

式中：Rn为n阶预测系数；w(k)为白噪声.此时随机误差序列可分解成可预测分量与不可预测分量.

2 射击误差序列状态方程构建

射击误差序列的一步预测特性是影响毁伤概率的一个决定性因素.对一个正常的射击过程而言，其射击误差是由两部分构成：射击诸元误差，它是经历了足够长的时间后，暂态分量可以忽略了的相关平稳序列；射弹散布误差，它虽然也是平稳序列，但是它是在击发后才加入到射击误差序列之中，并且在有限次数的发射后，即行终止的平稳序列.

2.1 射击误差序列的统计特性

对射击误差而言，它的不相关、弱相关和强相关3个平稳随机序列中的每一个都是由更多的独立的子序列之和构成.对误差序列的一步预测方程,可表示为

(5)

(6)

并且有

(7)

(8)

两个独立序列和的均值函数是其k步预测值，而协方差函数是其k步预测方差.两个随机序列之和存在一个过渡过程.当0<|r|<1且k→∞时，此过渡过程才能结束，两个独立序列的和的均值函数与方差函数才能如同随机变量一样求和与方差.显然，只有两种情况是例外，即两个独立序列均为强相关序列或均为不相关序列，这是因为在这两种情况下，其任何一步的预测均值与方差都是不变的，且与其初值相同.

如果射击误差序列中的弱相关分量的过渡过程可以在射击准备时间内衰减到可以忽略的条件，射击过程中的所有弱相关分量与强相关分量才可以依随机变量求和规则予以处理.例如：射击诸元误差xs(k)有3个独立序列构成，即

xs(k)=xg(k)+xw(k)+xq(k)

(9)

(10)

且有xs(k)的一步预测系数表达式为

(11)

2.2 连续脱靶条件下射击诸元误差的密度函数

f(xs(M),xs(M-1),…,xs(n),…,xs(1))=

(12)

Xs(k-1)=

[xs(k-1)xs(k-2) …xs(k-n)]T

(13)

若射击诸元误差的最高阶次为n，根据式(12)给出的处理方法，xs(k)在当k>n时，则预测样本空间为

{xs(k)|X(k-1)}∈

(14)

此时xs(k)为一种双重正态分布，它作为射击误差的动态均值，与前面的各个动态均值都有关.

由于k∈[1,n]期间，高炮武器控制系统完成了振动的过渡过程，当k>n，xs(k)进入平稳序列，表明高炮随时可以实施射击.

对首次击发的第一枚弹药，由于射弹散布误差xb(1)使x(1)=xb(1)+xs(1)的方差产生突变，此时，x(1)的条件密度函数可表述为

f(x(1))=f(x(1)|xs(1)f(xs(1)))

(15)

因而在首发弹药脱靶条件下，射击诸元误差的密度函数为

g1(xs(1))=

(16)

式中:s1为当k=1时刻靶标的迎弹面，其补集为R1.因此首发弹药脱靶的概率为

(17)

根据条件概率密度的定义，连续两次发射均脱靶条件下，射击诸元误差的密度函数为

g1(xs(1))dxs(1)

(18)

而连续两发均脱靶的概率为

P(x(1)∉s1,x(2)∉s2)=

(19)

依据上述推导方式，可以得到前N发弹药均脱靶条件下，射击诸元误差的密度函数为

gN-1(xs(N-1))dxs(N-1)=

gN-1[xs(N-1)]dxs(N-1)

(20)

上述密度函数的积分即为连续发射N发弹药均脱靶的概率.对于满足0-1毁伤律的高炮武器系统，在连续发射N发弹药的条件下，至少命中一发的概率可表示为

(21)

如果将毁伤定义为命中，则式(21)为射击诸元误差存在相关性时，高炮武器进行N发连续射击的毁伤概率.

3 单管高炮武器系统的毁伤概率

单管高炮武器系统是指在同一瞬时仅能发射一发弹药.对于转管火炮，虽然在外形上具有多个发射管，但每个火炮身管只能在转到特定位置后，依单管规则顺序发射，因此在计算毁伤概率时，仍然可按照单管方式加以计算.

(22)

再考虑到射弹散布误差序列是在开始时与射击准备误差序列求和的，故必须以已知射击诸元Zs(k)条件下的分布特性来表述.基于上述分析可知第一发弹药射击诸元误差的密度函数为

(23)

从式(23)出发，可递推出前k发均不毁伤条件下，射击诸元误差的密度函数为

gk(Zs(k))=

gk-1(Zs(k-1))dZs(k-1)

(24)

则单管高炮武器一次N发的点射的毁伤概率为

(25)

式中：Σb和Σs为协方差矩阵；det为行列式；I为单位矩阵.

4 数值仿真与分析

当射击误差中的弱相关分量为1阶平稳序列模型时，分别采用国军标、Monte Carlo(MC)模拟、文献[6]递推法以及本文方法进行计算，仿真计算结果如图1所示.当射击误差中的弱相关分量采用3阶平稳序列模型时，误差序列的各阶自相关系数为：ri,i=0.81,ri,i+1=0.72,ri,i+2=0.66 (i=1,2,…,21)，仿真计算结果如表1所示.当射击误差中的弱相关分量采用4阶平稳序列模型时，误差序列的各阶自相关系数为：ri,i=0.81,ri,i+1=0.72,ri,i+2=0.66，ri,i+3=0.60 (i=1,2,…,20)，仿真计算结果如表2所示.

图1 相关系数与点射毁伤概率关系对比图Fig.1 Comparison of bust firing damage probability with correlation coefficient

表1 具有3阶自相关性的射击误差序列毁伤概率Tab.1 Damage probabilities of firing error sequence with third-order autocorrelation

表2 具有4阶自相关性的射击误差序列毁伤概率Tab.2 Damage probabilities of firing error sequence with fourth-order autocorrelation

上述仿真结果表明，当射击误差序列考虑为1阶平稳序列模型时，4种计算方法的结果毁伤概率结果差异性很小.但随着误差相关性阶数的增加，国军标的计算结果与其他3种存在较为明显的差距，这表明将高阶自相关性误差序列模型近似为1阶模型，会产生一定的模型计算误差.而运用本文给出的预测系数计算方法可以较好地适应各类高阶模型，计算结果更加符合射击误差的实际情况.由于仿真中的各弱相关误差序列的随机特性已经给定，采用MC方法模拟射击过程中的各弱相关误差源，可以忽略模型误差，通过大量的模拟抽样所得到的数值计算结果是可信的.因此可以将MC方法的计算结果作为比较基准进行准确性验证，通过表1和表2计算结果对比可看出，国军标方法有20%以上的偏差，文献[6]方法的偏差为10%以上，而本文计算结果的偏差小于3%.

对于计算过于复杂而难以求得解析解的随机过程问题，虽然MC模拟可以通过构造符合一定规则的随机数来进行数值求解，但正如前文所述，考虑到高炮武器系统验证的费效比问题，高炮武器系统从设计到部队列装，不可能通过大量的试验进行毁伤概率的验证，在试验数据较少的情况下，使用MC模拟射击打靶过程，存在模型误差的问题，有可能产生较大的计算误差.为进一步说明此问题，利用GPML Toolbox version 4.1工具箱在MATLAB R2013a下，生成一个均值为0，方差为1的高斯随机过程，并产生 10 000 组样本，每个样本均匀采样20个点，统计采样点落在某一区间的占比，用于模拟高炮打靶发射弹药数与命中目标数的比例，即为命中概率.根据MC模拟的原则，分别抽取不同数量的样本m1，每组总采样点个数为20m1，记录采样点落在[-0.5，0.5]之间的个数为n1，记p1=n1/(20m1)，称为抽样落点占比.而总的采样点个数为 200 000 个，统计得到全体采样点落在[-0.5，0.5]之间的个数为 93 662 个，记p2=93 662/200 000=0.468，称为全体落点占比.不同采样数量下的落点占比与全体落点占比的统计结果如图2所示.

图2 MC算法抽样数量与落点占比的关系Fig.2 Number of samples versus proportion of dropping points in MC algorithm

从图2中可知，当抽样数量较少时，MC模拟的结果具有较大的不确定性.模拟过程中还发现，即使抽样数量相同的不同样本，MC模拟的数值统计结果也存在明显差异，这对于试验数据较少，且要求毁伤概率计算具有严格数学推理的武器研制与验证单位是不能接受的.而本文的递推估计模型具有严格的数学推导，只需确定弱相关误差的均值、方差和相关系数，就可以准确计算出毁伤概率.

5 结语

提出了一种利用高炮武器射击误差的状态方程估计毁伤概率的方法，在将射击误差序列中的弱相关部分，分解为可预测误差和不可预测误差的，给出了可预测系数的计算方法，并利用递推方法建立了高炮武器连续点射的毁伤概率计算模型，可以较好地适用于射击误差序列具有多阶自相关性时高炮武器系统的毁伤概率计算.本方法也可以作为各类速射火炮进行连续射击时毁伤概率计算的一种参考.