2024年国赛高教杯数学建模D题反潜航空深弹命中概率问题解题全过程文档及程序

2024年国赛高教杯数学建模

D题反潜航空深弹命中概率问题

原题再现

应用深水炸弹（简称深弹）反潜，曾是二战时期反潜的重要手段，而随着现代军事技术的发展，鱼雷已成为现代反潜作战的主要武器。但是，在海峡或浅海等海底地形较为复杂的海域，由于价格低、抗干扰能力强，仍有一些国家在研究和发展深水炸弹反潜技术。反潜飞机攻击水下目标前，先由侦察飞机通过电子侦察设备发现水下潜艇目标的大致位置，然后召唤反潜飞机前来进行攻击。当潜艇发现被侦察飞机电子设备跟踪时，通常会立即关闭电子设备及发动机，采取静默方式就地隐蔽。
本问题采用目标坐标系：潜艇中心位置的定位值在海平面上的投影为原点 𝑂，正东方向为 𝑋 轴正向，正南方向为 𝑌 轴正向，垂直于海平面向下方向为 𝑍 轴正向。正北方向顺时针旋转到潜艇航向的方位角记为 𝛽，假定在一定条件下反潜攻击方可获知该航向（见图1）。
在这里插入图片描述
由于存在定位误差，潜艇中心实际位置的3个坐标是相互独立的随机变量，其中 𝑋，𝑌均服从正态分布 𝑁(0,𝜎2)，𝑍 服从单边截尾正态分布 𝑁(ℎ0,𝜎𝑧2,𝑙)，其密度函数为

这里 ℎ0 是潜艇中心位置深度的定位值，𝑙 是潜艇中心位置实际深度的最小值，𝜙 和 𝛷 分别是标准正态分布的密度函数与分布函数。将潜艇主体部分简化为长方体，深弹在水中垂直下降。假定深弹采用双引信（触发引信+定深引信）引爆，定深引信事先设定引爆深度，深弹在海水中的最大杀伤距离称为杀伤半径。深弹满足以下情形之一，视为命中潜艇：
(1) 航空深弹落点在目标平面尺度范围内，且引爆深度位于潜艇上表面的下方，由触发引信引爆；
(2) 航空深弹落点在目标平面尺度范围内，且引爆深度位于潜艇上表面的上方，同时潜艇在深弹的杀伤范围内，由定深引信引爆；
(3) 航空深弹落点在目标平面尺度范围外，则到达引爆深度时，由定深引信引爆，且此时潜艇在深弹的杀伤范围内。
请建立数学模型，解决以下问题：
问题1 投射一枚深弹，潜艇中心位置的深度定位没有误差，两个水平坐标定位均服从正态分布。分析投弹最大命中概率与投弹落点平面坐标及定深引信引爆深度之间的关系，并给出使得投弹命中概率最大的投弹方案，及相应的最大命中概率表达式。针对以下参数值给出最大命中概率：潜艇长100 m，宽20 m，高25 m，潜艇航向方位角为 90∘，深弹杀伤半径为20 m，潜艇中心位置的水平定位标准差 𝜎=120 m，潜艇中心位置的深度定位值为150 m.
问题2 仍投射一枚深弹，潜艇中心位置各方向的定位均有误差。请给出投弹命中概率的表达式。针对以下参数，设计定深引信引爆深度，使得投弹命中概率最大：潜艇中心位置的深度定位值为150 m，标准差 𝜎𝑧=40 m，潜艇中心位置实际深度的最小值为 120 m，其他参数同问题1。
问题3 由于单枚深弹命中率较低，为了增强杀伤效果，通常需要投掷多枚深弹。若一架反潜飞机可携带9枚航空深弹，所有深弹的定深引信引爆深度均相同，投弹落点在平面上呈阵列形状（见图2）。在问题2的参数下，请设计投弹方案（包括定深引信引爆深度，以及投弹落点之间的平面间隔），使得投弹命中（指至少一枚深弹命中潜艇）的概率最大。
在这里插入图片描述

整体求解过程概述(摘要)

随着潜艇技术的不断发展，反潜作战面临着日益严峻的挑战。深弹作为重要的反潜武器，其投放策略的优化对于提高反潜作战效率至关重要。本文旨在通过数学建模和数值分析，找出在不同条件下深弹投放的最大命中概率策略。
潜艇深度已知的单深弹投放模型，在潜艇深度位置信息确定的情况下，本文首先在二维平面上建立了深弹的毁伤概率模型。该模型考虑了深弹爆炸的威力范围、潜艇的尺寸以及投放角度等因素，将问题转化为炸弹最大毁伤区域的构建与最优概率积分问题。根据深弹爆炸深度与潜艇深度的位置关系，本文详细分析了五种可能的投放情况，包括深弹在潜艇上方、下方、侧面以及直接命中等。通过动态讨论，确定了当深弹的落点位于潜艇的上表面和下表面之间时，投弹可达最大投弹概率。
潜艇深度未知的单深弹投放模型，考虑到实际作战中潜艇深度往往未知，本文在原有模型基础上增加了深度参数，并假设其服从单边截尾正态分布。这使得模型更加贴近实际，也增加了问题的复杂性。为了处理深度不确定性，本文建立了三维立体的命中概率模型。通过数值分析和简化处理，将三维模型分解为二维命中概率模型与一维深度概率模型的最大化积分问题。在二维命中概率模型的基础上，本文利用数值积分方法求解了潜艇深度定位在最大可能区间的命中概率，为实际作战中的深弹投放提供了决策支持。
潜艇深度未知的多深弹投放模型，首先通过数学分析确定了深弹的最佳爆炸深度，以确保在给定深度范围内达到最大的毁伤效果。接着，本文将最大命中概率的决定因素分解为深弹二维毁伤区域与潜艇的平面分布概率两部分。通过综合考虑深弹的投放位置、爆炸范围以及潜艇的可能位置，构建了多深弹投放的命中概率优化模型。利用数学解析和数值解析与网格算法相结合的方法，本文求解了多深弹投放模型下的最大命中概率。当深弹间距d=150时，命中概率达到最大值0.333。这一结果对于指导实际作战中的多深弹投放具有重要意义。

模型假设：

1.潜艇以及深弹默认为单一质点；
2.假设深弹之间不会互相影响；
3.假设潜艇时刻处在动态移动中。

问题分析：

首先，我们分析了深弹引爆深度与命中概率之间的内在联系。在概率密度函数已知且非负的前提下，求解最大投射命中概率被转化为求解最大积分区域的问题。通过数学推导，我们证明了只有当引爆深度位于潜艇上下表面之间时，积分区域才能达到最大，进而确定了固定投弹点处的命中概率。接着，我们考虑了引爆深度变化时积分区域虽然大小相同但命中概率差异显著的情况。通过将该问题转化为二维正态分布函数的最值求解问题，我们得到了最大命中概率的表达式。在已知潜艇长宽高等参数（深度定位参数无误差）且深弹杀伤半径为20m的条件下，我们计算出了最大的命中概率。
在潜艇深度未知的情况下，我们引入了深度定位误差，并假设其服从单边截尾正态分布。为了简化问题，我们将问题转化为在潜艇最有可能出现的深度范围内选择最佳投弹深度的问题。基于问题一的研究成果，我们已知潜艇深度一定时的最佳投弹深度和投弹点。因此，在已知潜艇中心位置深度的定位误差服从的分布参数后，我们以潜艇的高度为度量区间，找出了该区间内最有可能出现的深度范围，并取其中点作为深弹的定深引信引爆深度。这样的选择使得投弹命中概率达到最大。
在问题二的基础上，我们进一步考虑了多深弹的投放布局问题。给定深弹的定深引信引爆深度，我们根据潜艇深度参数服从的单边截尾正态分布，计算了潜艇中心落在不同区域内的概率。通过参数调整，我们计算了每个区域内9枚航空深弹至少一次命中潜艇的概率。然后，将这两组概率对应相乘并求和，得出了命中概率的数学期望。最后，通过数值优化方法，我们找出了使命中概率数学期望最大的投弹方案，即最佳投弹布局。

模型的建立与求解整体论文缩略图

在这里插入图片描述

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码(完整论文以及代码请联系博主)：

sigma = 120; % 标准差，单位：米
L = 100; % 潜艇长度，单位：米
R = 20; % 杀伤半径，单位：米
W = 20; % 潜艇宽度，单位：米

% 定义被积函数
f = @(x, y) exp(-(x.^2 + y.^2) / (2 * sigma^2));

% 计算第一个积分项（左侧部分）
I1 = integral2(f, -R-L/2, -L/2, @(x)-W/2-sqrt(R^2 - (x + L/2).^2), @(x)W/2+sqrt(R^2 - (x + L/2).^2));

% 计算第二个积分项（中间部分）
I2 = integral2(f, -L/2, L/2, -W/2-R, W/2+R);

% 计算第三个积分项（右侧部分）
I3 = integral2(f, L/2, R+L/2, @(x)-W/2-sqrt(R^2 - (x - L/2).^2), @(x)W/2+sqrt(R^2 - (x - L/2).^2));

% 计算总积分
p00 = (1 / (2 * pi * sigma^2)) * (I1 + I2 + I3);

% 显示结果
disp(['The probability p(0,0) is: ', num2str(p00)]);

function [d,I]=x2
sigma = 120; % 标准差
L = 100; % 潜艇长度
R = 20; % 杀伤半径
W = 20; % 潜艇宽度
H = 25; % 高度
sigma_z = 40; % Z 的标准差
l1 = 120;
h0 = 150;

f = @(x, y) (1 / (2 * pi * sigma^2)) * exp(-(x.^2 + y.^2) / (2 * sigma^2));
Phi = @(x) normcdf(x, 0, 1);
dm=1/(1 - Phi((l1 - h0) / sigma_z));
% 定义函数 g(z)
g_z = @(z) (1/sigma_z)*dm * (1 / sqrt(2 * pi)) * exp(-((z - h0).^2) / (2 * sigma_z^2));
%test=integral(@(z) g_z(z),120,200);
fun = @(x,y,z) f(x,y).*g_z(z);

d = 152.5:1:180;

I1 = arrayfun(@(d) integral3(@(x, y, z) f(x, y) .* g_z(z), -L/2, L/2, -W/2, W/2, l1, d-R-H/2), d);
I2=[];
I3=[];
I4=[];
I5=arrayfun(@(d) integral(@(z) g_z(z), d-H/2,d+H/2), d);
I5=0.083734*I5;
I6=[];
I7=[];
I8=[];

for i=1:length(d)
    dx=0.5;
    dy=0.5;
    dz=0.5;
    %%%以下计算 I2
    % 初始化黎曼和
    sum=0;
    % 计算黎曼和
    zmin = d(i) - R - 0.5 * H;
    zmax = d(i) - 0.5 * H;
    xmin = @(z) -L/2-sqrt(R^2 - (d(i) - z - H/2).^2);
    xmax = @(z) -L/2;
    ymin = @(x,z) -W/2-sqrt(R^2 - (d(i) - z - H/2).^2-(x+L/2).^2);
    ymax = @(x,z) W/2+sqrt(R^2 - (d(i) - z - H/2).^2-(x+L/2).^2);
    for z = zmin:dz:zmax
        x1=xmin(z);
        xu=xmax(z);
        for x =x1:dx:xu
            y1=ymin(x,z);
            yu=ymax(x,z);
            for y = y1:dy:yu
                sum = sum + fun(x,y,z) * dx * dy * dz;
            end
        end
    end
    I2=[I2 sum];
end

%%%以下计算 I3
% 初始化黎曼和
sum=0;
% 计算黎曼和
zmin = d(i) - R - 0.5 * H;
zmax = d(i) - 0.5 * H;
xmin = -0.5*L;
xmax = 0.5*L;
ymin = @(z) -W/2-sqrt(R^2 - (d(i) - z - H/2).^2);
ymax = @(z) W/2+sqrt(R^2 - (d(i) - z - H/2).^2);
sum = 0;
for z = zmin:dz:zmax
    xl=xmin;
    xu=xmax(z);
for x =x1:dx:xu
    yl=ymin(z);
    yu=ymax(z);
    for y = y1:dy:yu
        sum = sum + fun(x,y,z) * dx * dy * dz;
    end
end
end
I3=[I3 sum];

%%%以下计算 I4
% 初始化黎曼和
sum=0;
% 计算黎曼和
zmin = d(i) - R - 0.5 * H;
zmax = d(i) - 0.5 * H;
xmin = @(z) L/2;
xmax = @(z) L/2+sqrt(R^2 - (d(i) - z - H/2).^2);
ymin = @(x,z) -W/2-sqrt(R^2 - (d(i) - z - H/2).^2-(x-L/2).^2);
ymax = @(x,z) W/2+sqrt(R^2 - (d(i) - z - H/2).^2-(x-L/2).^2);
sum = 0;
for z = zmin:dz:zmax
    xl=xmin(z);
    xu=xmax(z);
    for x =x1:dx:xu
        yl=ymin(x,z);
        yu=ymax(x,z);
        for y = y1:dy:yu
            sum = sum + fun(x,y,z) * dx * dy * dz;
        end
    end
end
I4=[I4 sum];

%%%以下计算 I6
% 初始化黎曼和
sum=0;
% 计算黎曼和
zmin = d(i) + 0.5 * H;
zmax = d(i) +R+ 0.5 * H;
xmin = @(z) -L/2-sqrt(R^2 - (d(i) - z + H/2).^2);
xmax = @(z) -L/2;
ymin = @(x,z) -W/2-sqrt(R^2 - (d(i) - z + H/2).^2-(x+L/2).^2);
ymax = @(x,z) W/2+sqrt(R^2 - (d(i) - z + H/2).^2-(x+L/2).^2);
for z = zmin:dz:zmax
    xl=xmin(z);
    xu=xmax(z);
    for x =x1:dx:xu
        yl=ymin(x,z);
        yu=ymax(x,z);
        for y = y1:dy:yu
            sum = sum + fun(x,y,z) * dx * dy * dz;
        end
    end
end

1. import numpy as np
 2. from scipy.stats import truncnorm
 3. from scipy.integrate import quad
 4. import pandas as pd
 5. import matplotlib.pyplot as plt
 6. import matplotlib
13. h0 =200
14. sigma_z =50 
15. l= 150
18. def truncated_normal_density(z):
 19. a, b= (l-h0) /sigma_z, float('inf') 
 20. scale = sigma_z
 21. loc =h0
 22. return truncnorm.pdf(z, a,b, loc=loc, scale=scale)
25. def calculate_integral(z0):
 26. lower_bound =z0-10
 27. upper_bound =z0+10
 28. result, _= quad(truncated_normal_density, lower_bound, upper_bound)
 29. return result
32. def main():
 33. z0_values = np.arange(100,200, 0.0001)
 34. integral_values =[calculate_integral(z0) for z0 in z0_values]
 37. df =pd.DataFrame({
 38. 'z0 (meters)': z0_values,
 39. 'Integral Value':integral_values
 40. })
 41. df.to_excel('result.xlsx',index=False)
44. plt.figure(figsize=(30, 20))
 45. plt.plot(z0_values, integral_values, marker='o', linestyle='*',color='y')
 46. plt.xlabel()
 47. plt.ylabel()
 49. plt.grid(True)
 50. plt.savefig('2.jpg')
 51. plt.show()
 53. if __name__ == "__main__":
 54. main()