基于MonteCarlo法的经典射击问题中的杀伤概率估计问题

• 问题的建模

  

​ 若假设一个火力单位有$m$件同种类型的武器，同步发射，一次发射时每件武器各发射出$n$发射弹。在这样的条件下，击毁一个目标$V$的概率，一般称为杀伤概率。若我们假设每发射弹弹着点的射击误差均位二维随机向量。假设第$i$件武器第$j$发射弹弹着点的射击误差向量是：
$$X_{ij}=Y_{ij}+Z+A_i+A$$
​ 其中，$X_{ij}$取值$x=(x_1,x_2)$;对不同件的武器，不同的发射弹是非重复误差$Y_{ij}$的取值$y=(y_1,y_2)$,且服从如下概率密度的二维正态分布:
$${\varphi }_1(y_1,y_2)=\frac{1}{2\pi |B{|}^{\frac{1}{2}}}{e}^{-\frac{1}{2}(y_1,y_2)B^{-1}(y_1,y_2{)}^T}$$
​ 其中$B$为二阶协方差矩阵。其中对于不同的武器，误差$Z$是重复误差取值为$z=(z_1,z_2)$,服从如下概率密度的二维正态分布：
$${\varphi }_2(z_1,z_2)=\frac{1}{2\pi |C{|}^{\frac{1}{2}}}{e}^{-\frac{1}{2}(z_1,z_2)B^{-1}(z_1,z_2{)}^T}$$
​ 其中$C$为二阶协方差矩阵。且$A_i$是第$i$件武器的系统误差，取值为$a_i=(a_{i1},a_{i2})$。$A$是火力单位的系统误差，取常数值$a=(a_1,a_2)$。比如说我们假设一发命中弹的杀伤概率为:
$$g=\frac{1}{{\omega }_0}$$
​ ${\omega }_0:$击毁目标所需要的平均命中弹数目。如果我们设计的目标$V$是边长为$2l$的正方形。因此具有$m$件

武器组成的火力单位，一次发射$n$发射弹，对目标$V$的杀伤概率:
$$R={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }Q(z_1,z_2){\varphi }_3(z_1,z_2)d{z}_{1}d{z}_2$$
​ 其中：$Q(z_1,z_2)$表示一次射击中$m$件武器$n$发导弹至少一发导弹击中$V(z_1,z_2)$的概率：
$$Q(z_1,z_2)=1-\prod _{i=1}^m[1-P_i(z_1,z_2){]}^n$$
​ $P_i(z_1,z_2)$表示第$i$种武器击中$V(z_1,z_2)$的概率：
$$P_i(z_1,z_2)={\iint }_{D}g{\varphi }_1(x_1-z_1-a_{i1}-a_1,x_2-z_2-a_{i2}-a_2)d\sigma$$
​ 其中$D$为边长为$2l$的正方形区域。其中${\varphi }_3(z_1,z_2)$表示击种目标$V(z_1,z_2)$的总的误差概率密度函数，是$X_{ij}$的概率密度函数。

求解以上模型理论上是计算其二重积分比较方便，详细的计算过程就不赘述，有想了解具体的计算过程的同学请参考我之前写的博客《通俗易懂的MonteCarlo积分方法(六)》

• 问题的仿真步骤

  在这里我们不用以前那种比较复杂的积分计算方法。我们将直接采用投点法对这整个射击过程进行模拟。先产生服从给定分布的随机变量，再用随机投点法求积。

  *STEP1:*产生二维正态分布随机向量$y_{ij}\sim {\varphi }_1(y_1,y_2),z_i\sim {\varphi }_2(z_1,z_2)$

  *STEP2:*计算代表一发导弹对的射击误差的随机向量的$x_{ij}=(x_{ij1},x_{ij2})$,其中：
  $$\begin{gathered} x_{ij1}=y_{ij1}+z_1+a_{i1}+a_1 \\ {x}_{ij2}=y_{ij2}+z_2+a_{i2}+a_2 \end{gathered}$$
  *STEP3:*计算$F_{ij}=[X^{\ast }-x_{ij};X^{\ast }+x_{ij}]$为实际打击区域。

  *STEP4:*判断$F_{ij}\in V?$($V$为被打击对象所在区域)如果$F_{ij}\in V,S_{ij}=1,$否则$S_{ij}=0$。

  *STEP5:*计算一个火力单位一次射击命中的总次数$S:$
  $$S=\sum _{i=1}^m\sum _{j=1}^n{S}_{ij}$$
  ​ 那么这个火力单位的命中率$H:$
  $$H=\frac{S}{mn}$$
  ​ 那么现在开始计算杀伤概率，只要这么多次尝试有一次击中$V$,便可以彻底摧毁$V$,那么总的 杀伤概率为:
  $$k=1-\prod _{i=1}^m[1-P_i(V)]$$
  ​ 其中$P_i(V)$为武器$i$能击中$V$的概率，在这里我们做$M$次仿真模拟,求得：
  $$P_i(V)=\frac{1}{Mn}\sum _{k=1}^M\sum _{j=1}^n{S}_{ijk}$$
  *STEP6:*重复实验$N$次，得到$k_i(i=1,2,3...)$ 其杀伤概率的均值为$R:$
  $$R=\frac{1}{N}\sum _{i=1}^N{k}_i$$

• 求解代码：

  我们不妨设有$m=4$,$n=10$,被打击的点为原点坐标为$(0,0)$,二维随机变量$Y_{ij}$的均值为$(0,0)$(炮口永远是对准目标的)。$m\times n$个单位的方差矩阵$Y_{{\sigma }_1},Y_{{\sigma }_2}$分别用随机数产生:

  >> clear
  >> Ysi1 = rand(4,10);
  >> Ysi2 = rand(4,10);
  

  $Y_{{\sigma }_1}$具体的值为：

  

  $Y_{{\sigma }_2}$具体的值为：

  

  >>  Zsi1 = rand(4,1);
  >> Zsi2 = rand(4,1);
  

  $Z_{{\sigma }_1}$的值为：

  

  $Z_{{\sigma }_2}$的值为：

  

协方差矩阵$O$为:
$$O=\left(\begin{array}{cc}{\sigma }_1^2& \rho {\sigma }_1{\sigma }_2\\ \rho {\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right)$$
我们设$\rho =0.8$。

我们分别设置$a_1,a_2=0$,即假设火力系统摆放位置等因素考虑确定，不存在系统误差。

我们设置$a_{i1},a_{i2}$分别由在$[-0.1,0.1]$上均匀分布随机数产生:

>> a_i1 = unifrnd(-0.1,0.1,[4 1]);
>> a_i2 = unifrnd(-0.1,0.1,[4 1]);

$a_{i1}$的值为:

$a_{i2}$的值为:

被打击目标的$l=0.5m$,即要打击一个$1m^2$的正方体$V$。

以下附上仿真代码：

function Attack
intialValue;%赋初值
m = 4;
n = 10;
l = 0.5;%一边长0.5m
rho = 0.8;
X_mean = 0;
Y_mean = 0;%被打击物体的坐标为(0,0);
N = 100;%仿真100次
M = 100;%子仿真100次
a1 = 0;
a2 = 0;
S = zeros(m,n,M);
P = zeros(1,m);
k = zeros(1,N);
Propablity = zeros(1,N);
for i = 1:N
    %模拟打击过程
   for j = 1:m
       z = mvnrnd([X_mean,Y_mean],XFC(rho,Zsi1(j),Zsi2(j)));
       for times = 1:M
         for k = 1:n
           y = mvnrnd([X_mean,Y_mean],XFC(rho,Ysi1(j,k),Ysi2(j,k)));
           x(1) = y(1)+z(1)+a_i1(j)+a1;
           x(2) = y(2)+z(2)+a_i2(j)+a2;
           F(1) = X_mean + x(1);%实际打击的一个坐标
           F(2) = Y_mean + x(2);%实际打击的另一个坐标
           if (-l<=F(1)<= l)&&(-l<=F(2)<= l)
               S(j,k,times) = 1;
           else
               S(j,k,times) = 0;
           end
         end
       end
   end
   %计算杀伤率
   Pro = 1;
   for j = 1:m
      PV = 0;
      for k = 1:n
          for times = 1:M
               PV = PV + S(j,k,times);
          end
      end
      P(j) =  PV/(M*n);
      Pro = Pro * (1-P(j));
   end
   k(i) = 1 - Pro;
   Propablity(i) = sum(k(1:i))/i;
end
   plot(1:N,Propablity,'g'); 
   title('仿真迭代曲线');
   xlabel('仿真次数');
   ylabel('杀伤率');
   disp(['杀伤率是：',num2str(Propablity(N))]);
end


function OS = XFC(rho,sigma1,sigma2)
OS = [sigma1^2,rho*sigma1*sigma2;rho*sigma1*sigma2,sigma2^2];
end

function intialValue 
    Ysi1 = rand(4,10);
    Ysi2 = rand(4,10);
    Zsi1 = rand(4,1);
    Zsi2 = rand(4,1);
    a_i1 = unifrnd(-0.1,0.1,[4 1]);
    a_i2 = unifrnd(-0.1,0.1,[4 1]);
end

然后就：

???合着我这半天白搞了???