火炮射击中的数学建模问题

火炮射击中的数学建模问题

第三题

考虑更实际的情况，火炮的射击的准确性受到诸多因素的综合影响，从而导致炮弹的落点与实际目标点产生偏差，称为射击误差。射击误差通常主要包含两类，一类称之为诸元误差，与火炮射击时的准备误差、气象数据误差等有关，可理解为系统误差，在一次射击任务中诸元误差是不变的，在多次任务中该误差并不相同；另一类称之为散布误差，是每一次射击过程中不可观不可控因素引起的随机偏差，可理解为随机误差。诸元误差和散布误差均服从二维正态分布，二维正态分布的两个维度分别为距离向和方位向，其中距离向是炮弹在射击方向上与目标点产生的偏差，方位向是炮弹在目标点的左右偏离产生的偏差。以诸元误差为例，设其分布参数为 ${\sigma }_f^2,{\sigma }_d^2$，则火炮在某次射击中诸元误差$x_e,y_e$ 满足：

$$f_e(xe,ye)=\frac{1}{2\pi {\sigma }_f{\sigma }_d}{e}^{-\frac{1}{2}(\frac{x_e^2}{{\sigma }_f^2}+\frac{y_e^2}{{\sigma }_d^2})}$$

如图 1 所示，O 点为火炮的目标点 $(x_0,y_0)$，Y 轴为距离向，X 轴为方位向，实际火炮落点存在一个诸元误差 $(x_e,y_e)$ 和散布误差 $(x_b,z_b)$ 叠加在其中，因此，最终火炮落点满足：
$$\{\begin{array}{ll}& \text{x=x0+xb+xe}\\ & \text{y=y0+yb+ye}\end{array}$$

图1：火炮射击误差示意图

对火炮的射击误差的测定，主要是测定其诸元误差及其分布参数以及散布误差的分布参数，以便对火炮的打击效果有更好的评估，进而根据作战目的和火炮的性能能够对兵力弹药有更好的测算。
（1）附件 2 中给出了某次任务中采集到的某门火炮射击 60 枚炮弹的落点位置（仅考虑二维情况），假设该门火炮的目标点为 (40, 60) ，请根据这些数据估计火炮的在这次射击中的诸元误差和散布误差的分布参数。
（2）附件 3 给出了某门火炮前后两次任务，每次任务各射击 60 枚炮弹的落点位置，第一次任务的目标点为 (40, 60) ，第二次任务重新对火炮对行了调整，将目标点调整为 (43, 70) ，请根据这些数据估计火炮的散布误差参数和两次任务中的诸元误差。

（1）

思路

可以根据附件给出的射击数据，通过最小二乘法来估计误差分布的参数。具体步骤如下：

1. 首先计算每个炮弹的距离向误差$\Delta x=x_i-x_0$ 和方位向误差 ​$\Delta y=y_i-y_0$。其中，$(x_0,y_0)$是目标点坐标，$(x_i,y_i)$分别是第$i$枚炮弹的实际落点坐标， $(x_e,y_e)$是诸元误差， $(x_b,z_b)$是散布误差。
2. 对距离向误差和方位向误差分别进行拟合：
   根据题目所给出的公式：
   $$f_e(x_e,y_e)=\frac{1}{2\pi {\sigma }_f{\sigma }_d}\exp \left(-\frac{1}{2}\left(\frac{x_e^2}{{\sigma }_f^2}+\frac{y_e^2}{{\sigma }_d^2}\right)\right)$$
   我们可以看出 $x_e$ 和 $y_e$ 的分布是独立的，因此可以将其拆分为 $x_e$ 和 $y_e$ 两个单变量的二元高斯分布：
   $$f(x_e)=\frac{1}{\sqrt{2\pi }{\sigma }_f}\exp \left(-\frac{x_e^2}{2{\sigma }_f^2}\right)$$
   $$f(y_e)=\frac{1}{\sqrt{2\pi }{\sigma }_d}\exp \left(-\frac{y_e^2}{2{\sigma }_d^2}\right)$$
   其中，${\sigma }_f$ 和 ${\sigma }_d$ 分别是距离向误差和方位向误差的标准差，分别对应诸元误差。

每个单变量高斯分布的拟合函数都可以采用最小二乘法来求得，即通过最小化垂直方向上的数据点到拟合函数的距离来确定模型参数。
考虑 $x_e$ 拟合函数，将样本数据 $x_1,x_2,\dots ,x_n$ 看做是从一个期望为 $\mu$、方差为 ${\sigma }_f^2$ 的高斯分布中随机产生的，于是可以构建出目标函数：
$$f(\mu ,{\sigma }_f)=\sum _{i=1}^n{\left(x_i-\frac{1}{\sqrt{2\pi }{\sigma }_f}\exp \left(-\frac{(x_i-\mu {)}^2}{2{\sigma }_f^2}\right)\right)}^2$$
其中 $x_i$ 是样本数据中第 $i$ 个数的值。通过最小化目标函数来求得参数 $\mu$ 和 ${\sigma }_f$ 的最优解，最终的拟合函数为：
$$f(x_e)=\frac{1}{\sqrt{2\pi }\widehat{{\sigma }_f}}\exp \left(-\frac{(x_e-\hat{\mu }{)}^2}{2{\widehat{{\sigma }_f}}^2}\right)$$
其中 $\hat{\mu }$ 和 $\widehat{{\sigma }_f}$ 是最优解。$y_e$ 的拟合函数可以类似进行求解。

最优解的求解

首先，我们需要选取一个初始值（${\mu }_0$ 和 ${\sigma }_{f0}$）来开始迭代。我们可以自行选择一个数值或者通过一些启发式方法（如根据数据的均值和方差来选择）来确定初始值。对于每次迭代，根据当前的参数 $\mu$ 和 ${\sigma }_f$ 的值来计算目标函数 $f(\mu ,{\sigma }_f)$ 的梯度，然后根据梯度方向更新参数：

$$\begin{array}{rl}{\mu }_{t+1}& ={\mu }_t-\eta \frac{\partial f({\mu }_t,{\sigma }_{ft})}{\partial \mu }\\ {\sigma }_{ft+1}& ={\sigma }_{ft}-\eta \frac{\partial f({\mu }_t,{\sigma }_{ft})}{\partial {\sigma }_f}\end{array}$$

其中 $t$ 是迭代次数，$\eta$ 是学习率，通常需要通过试错或者其他方法来确定。

现在，我们需要求出目标函数 $f(\mu ,{\sigma }_f)$ 对 $\mu$ 和 ${\sigma }_f$ 的偏导数。根据题目所给出的目标函数的表达式：

$$f(\mu ,{\sigma }_f)=\sum _{i=1}^n{\left(y_i-\frac{1}{\sqrt{2\pi }{\sigma }_f}\exp \left(-\frac{(x_i-\mu {)}^2}{2{\sigma }_f^2}\right)\right)}^2$$

我们可以分别对 $\mu$ 和 ${\sigma }_f$ 求偏导数。先是对 $\mu$ 求偏导数：

$$\begin{array}{rl}\frac{\partial f(\mu ,{\sigma }_f)}{\partial \mu }& =\sum _{i=1}^n-2\left(y_i-\frac{1}{\sqrt{2\pi }{\sigma }_f}\exp \left(-\frac{(x_i-\mu {)}^2}{2{\sigma }_f^2}\right)\right)\frac{1}{\sqrt{2\pi }{\sigma }_f}\exp \left(-\frac{(x_i-\mu {)}^2}{2{\sigma }_f^2}\right)\frac{x_i-\mu }{{\sigma }_f^2}\\ & =\sum _{i=1}^n-2\left(y_i-f(x_i,\mu ,{\sigma }_f)\right)\frac{x_i-\mu }{{\sigma }_f^2}\end{array}$$
其中 $f(x_i,\mu ,{\sigma }_f)=\frac{1}{\sqrt{2\pi }{\sigma }_f}\exp \left(-\frac{(x_i-\mu {)}^2}{2{\sigma }_f^2}\right)$。这里的偏导数表达式是基于链式法则得到的，借助于复合函数求导的方法，同理我们可以得到对 ${\sigma }_f$ 的偏导数：
$$\begin{array}{rl}\frac{\partial f(\mu ,{\sigma }_f)}{\partial {\sigma }_f}& =\sum _{i=1}^n-2\left(y_i-\frac{1}{\sqrt{2\pi }{\sigma }_f}\exp \left(-\frac{(x_i-\mu {)}^2}{2{\sigma }_f^2}\right)\right)\frac{1}{\sqrt{2\pi }{\sigma }_f^2}\exp \left(-\frac{(x_i-\mu {)}^2}{2{\sigma }_f^2}\right)\frac{(x_i-\mu {)}^2}{{\sigma }_f^3}\\ & =\sum _{i=1}^n-2\left(y_i-f(x_i,\mu ,{\sigma }_f)\right)\frac{(x_i-\mu {)}^2}{{\sigma }_f^3}-\left(y_i-f(x_i,\mu ,{\sigma }_f)\right)\frac{1}{{\sigma }_f^2}\end{array}$$
我们现在已经得到了梯度的表达式，可以使用梯度下降方法或者其他优化算法来迭代求解最优解。

3.对拟合函数分别进行最小二乘拟合，求出标准差 ${\sigma }_f,{\sigma }_d$;

根据这些数据，可以估计出这门火炮在这次射击中的诸元误差和散布误差的分布参数。

import math
import numpy as np

def f_e(xe, ye, sigma_f, sigma_d):
    exponent = -0.5 * ((xe ** 2) / (sigma_f ** 2) + (ye ** 2) / (sigma_d ** 2))
    coefficient = 1 / (2 * math.pi * sigma_f * sigma_d)
    return coefficient * np.exp(exponent)

import numpy as np

def calculate_expectation_and_variance(data):
    # 计算期望
    expectation = np.mean(data)
    
    # 计算方差
    variance = np.var(data)
    
    return expectation, variance

# 测试数据
data = [4.2, 5.5, 6.7, 7.2, 8.4, 9.3, 10.1, 11.0, 12.4, 13.8]

# 计算期望和方差
expectation, variance = calculate_expectation_and_variance(data)

# 打印结果
print("期望：", expectation)
print("方差：", variance)

第四题

问题转化：

我们需要通过目标点和落点数据来判定炮弹落点可能来自于哪一门火炮。由于无法获取火炮的诸元误差和散布误差参数，我们需要建立一个模型来进行判定。我们可以将该问题视为一个分类问题，即分类一个炮弹落点所属的火炮编号。

模型建立：

我们可以使用K近邻分类器来解决这个问题。K近邻分类器是一种有监督学习算法，能够利用已知类别的训练数据集对新的数据进行分类。它的基本思想是：如果一个样本在特征空间中的K个最相邻的样本中的大多数属于某一个类别，则该样本也属于该类别。

模型实现：

1. 数据预处理
   将目标点和落点坐标合并为一个18维的特征向量，即将数据的特征维度从2维增加到18维，我们首先需要将原始数据进行处理，生成18维特征向量。可以使用numpy库中的concatenate()函数将两个数组按列合并起来：

import numpy as np
# 加载数据
target = np.loadtxt('target.txt')
impact = np.loadtxt('impact.txt')
# 将目标点和落点坐标按列合并为一个数组
data = np.concatenate([target, impact], axis=1)

2. 数据划分
   将数据集划分为训练集和测试集。我们可以使用sklearn库中的train_test_split()函数将数据随机分为训练集和测试集：

from sklearn.model_selection import train_test_split
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(data, label, test_size=0.2, random_state=42)

3. 模型训练
   使用训练数据对K近邻分类器进行训练，并对测试集进行验证。可以使用sklearn库中的KNeighborsClassifier类进行训练和验证：

from sklearn.neighbors import KNeighborsClassifier
# 创建K近邻分类器模型
knn = KNeighborsClassifier(n_neighbors=3)
# 在训练集上拟合K近邻分类器
knn.fit(X_train, y_train)
# 在测试集上进行预测
y_pred = knn.predict(X_test)

这里我们设置K值为3（即只考虑最近的3个邻居），可以根据实际情况调整K值。在训练集上拟合K近邻分类器后，我们可以使用predict()函数对测试集进行预测，并得到预测结果。
4. 模型评估
使用准确率作为模型的评估指标。准确率是指分类器预测结果正确的样本数占总样本数的比例。可以使用sklearn库中的accuracy_score()函数来计算准确率：

from sklearn.metrics import accuracy_score
# 计算准确率
accuracy = accuracy_score(y_test, y_pred)
print('准确率：{:.2f}%'.format(accuracy*100))

5. 自行生成数据进行验证：
   我们可以自行生成一批数据来验证模型的准确性。生成数据的方法可以如下：
   首先，我们从目标点中随机选择一门火炮，然后在选中的火炮附近随机生成100个落点，模拟100个炮弹的轨迹。根据这100个炮弹的落点数据，将这些数据输入我们训练好的K近邻分类器模型，预测炮弹落点可能来自于哪一门火炮，最后将预测结果与该火炮的编号进行比较，判断模型的准确性。
   代码如下：

import random

# 随机选择一门火炮
selected_gun = random.randint(0, 17)

# 在选中的火炮附近随机生成100个落点
test_impact = []
for i in range(100):
    target = target[selected_gun]
    impact_x = np.random.normal(target[0], 100)
    impact_y = np.random.normal(target[1], 100)
    test_impact.append([impact_x, impact_y])
test_impact = np.array(test_impact)

# 将目标点和落点坐标按列合并为一个数组
test_data = np.concatenate([np.tile(target[selected_gun], (100, 1)), test_impact], axis=1)

# 对测试数据进行预测
test_label_pred = knn.predict(test_data)

# 统计预测结果
pred_count = np.bincount(test_label_pred)

# 输出预测结果
print('预测结果：', pred_count.argmax())
print('实际结果：', selected_gun)

这里我们随机选择一门火炮，然后在该火炮的目标点附近随机生成100个落点，模拟100个炮弹的轨迹。根据这100个炮弹的落点数据，将这些数据输入我们训练好的K近邻分类器模型，预测炮弹落点可能来自于哪一门火炮。最后将预测结果与该火炮的编号进行比较，输出预测结果和实际结果。
如果模型准确性较高，则预测结果和实际结果应该接近。如果模型准确性较低，则预测结果和实际结果可能相差较大。可以根据多次的预测结果来评估模型的准确性。