第九届“数维杯”大学生数学建模挑战赛（A题）深度剖析|建模完整过程+详细思路+代码全解析

作者：从前慢现在也慢 | 2024-05-12 03:30:30

踩

问题1

根据题目信息，可以将机会信号分为五类，分别是达到时间信息（TOA）、到达时间差信息（TDOA）、多普勒频率差信息（DFD）、到达角度信息（AOA）和接收强度指标信息（RSSI）。每一类机会信号都可以用数学表达式来表示，通过建立数学模型，可以推导出每一类信号所需的最少信号个数。

（1）达到时间信息 TOA：
根据题目信息，信号传输时间可以表示为： $T=TOA=T_{0}+\frac{L}{c}$ ，其中 $T_{0}$ 为信号发射时间，L为信号传播距离，c为信号传播速度。
设飞行器位置为 $x_{0},y_{0},z_{0})$ ，发射源位置为 $x_{1},y_{1},z_{1})$ ，接收时间为 $t_{0}$ ，则有：

$T_{0}=t_{0}-\sqrt{(x_{0}-x_{1})^{2}+(y_{0}-y_{1})^{2}+(z_{0}-z_{1})^{2}}$

信号传播距离L为：

$L=\sqrt{(x_{0}-x_{1})^{2}+(y_{0}-y_{1})^{2}+(z_{0}-z_{1})^{2}}$

因此，达到时间信息TOA可以表示为：

$TOA=t_{0}-\frac{1}{c}\sqrt{(x_{0}-x_{1})^{2}+(y_{0}-y_{1})^{2}+(z_{0}-z_{1})^{2}}$

（2）到达时间差信息 TDOA：
同一信号从两个发射源（同时发射）到达接收端的时间差可以表示为：

$t_{1}-t_{2}=\frac{L_{1}-L_{2}}{c}$

其中， $t_{1}$ 、 $t_{2}$ 分别为信号到达两个发射源的时间， $L_{1}$ 、 $L_{2}$ 分别为信号从发射源到接收端的距离。

设发射源1位置为 $x_{1},y_{1},z_{1})$ ，发射源2位置为 $x_{2},y_{2},z_{2})$ ，接收时间为 $t_{0}$ ，则有：

$t_{1}=t_{0}-\sqrt{(x_{0}-x_{1})^{2}+(y_{0}-y_{1})^{2}+(z_{0}-z_{1})^{2}}$

$t_{2}=t_{0}-\sqrt{(x_{0}-x_{2})^{2}+(y_{0}-y_{2})^{2}+(z_{0}-z_{2})^{2}}$

因此，到达时间差信息TDOA可以表示为：

$TDOA=\frac{\sqrt{(x_{0}-x_{1})^{2}+(y_{0}-y_{1})^{2}+(z_{0}-z_{1})^{2}}-\sqrt{(x_{0}-x_{2})^{2}+(y_{0}-y_{2})^{2}+(z_{0}-z_{2})^{2}}}{c}$

（3）多普勒频率差信息 DFD：
同一信号从两个发射源发射，由于飞行器与信号源具有相对速度，接收信号会产生频率变化，进而产生信号频率差。根据多普勒效应，有：

$\Delta f=\frac{f_{0}}{c}v_{r}$

其中， $\Delta f$ 为信号频率差， $f_{0}$ 为信号发射频率， $c$ 为信号传播速度， $v_{r}$ 为飞行器与信号源的相对速度。

设发射源1位置为 $x_{1},y_{1},z_{1})$ ，发射源2位置为 $x_{2},y_{2},z_{2})$ ，接收时间为 $t_{0}$ ，则有：

$f_{1}=f_{0}+\frac{v_{r1}}{c}f_{0}$

$f_{2}=f_{0}+\frac{v_{r2}}{c}f_{0}$

其中， $f_{1}$ 为信号从发射源1发射后的频率， $f_{2}$ 为信号从发射源2发射后的频率， $v_{r1}$ 为飞行器与发射源1的相对速度， $v_{r2}$ 为飞行器与发射源2的相对速度。

因此，多普勒频率差信息DFD可以表示为：

$DFD=\frac{f_{1}-f_{2}}{f_{0}}=\frac{v_{r1}-v_{r2}}{c}$

（4）到达角度信息 AOA：
接收源可得到发射源信号的相对角度信息，如图 1 所示（发射源与接收源连线在 xOy 平面投影线段与x 轴正方向的夹角，发射源与接收源连线与 z 轴负方向的夹角β）。为了简化处理，图 1 中的坐标系为地面惯性系，到达角度信息为角度的正切值。
设飞行器位置为 $x_{0},y_{0},z_{0})$ ，发射源位置为 $x_{1},y_{1},z_{1})$ ，接收角度为 $\theta$ ，则有：

$\tan \theta=\frac{z_{1}-z_{0}}{\sqrt{(x_{0}-x_{1})^{2}+(y_{0}-y_{1})^{2}}}$

因此，到达角度信息AOA可以表示为：

$AOA=\frac{z_{1}-z_{0}}{\sqrt{(x_{0}-x_{1})^{2}+(y_{0}-y_{1})^{2}}}$

（5）接收强度指标信息RSSI：
设信号源发射功率为 $P_{t}$ ，接收端的接收功率为 $P_{r}$ ，则接收强度指标信息RSSI可以表示为：

$RSSI=10\log_{10}\frac{P_{t}}{P_{r}}$

根据机会信号的数学表达式，可以发现，TOA和RSSI信息无法唯一确定飞行器的位置，因为它们都只提供了一维信息（时间和强度）。而TDOA、DFD和AOA信息可以提供二维信息（时间差、频率差和角度差），因此最少需要两个信号来确定飞行器的位置。但是考虑到实际情况中可能存在的偏差，为了保证定位的精确性，使用至少三个信号来确定飞行器位置。

Python示例代码实现：

import numpy as np
import pandas as pd
import math

# 读取附件1中的数据
data = pd.read_csv("附件1.csv")

# 定义机会信号的数学表达式
def TOA(t_r, t_s):
    return t_r - t_s

def TDOA(t_r1, t_s1, t_r2, t_s2):
    return (t_r1 - t_s1) - (t_r2 - t_s2)

def DFD(f_r, f_s):
    return f_r - f_s

def AOA(x, y, x_s, y_s):
    return math.atan((x - x_s) / (y - y_s))

def RSSI(P_s, P_r):
    return P_s / P_r

# 计算飞行器的位置
def calculate_position(data):
    # 初始化位置坐标
    x = 0
    y = 0

    # 计算TDOA，由于需要三个信号，所以从第三个数据开始计算
    for i in range(2, len(data)):
        # 计算TOA
        toa1 = TOA(data.iloc[i,1], data.iloc[i,3])
        toa2 = TOA(data.iloc[i,2], data.iloc[i,4])
        toa3 = TOA(data.iloc[i,3], data.iloc[i,5])

        # 计算TDOA
        tdoa1 = TDOA(data.iloc[i,1], data.iloc[i,3], data.iloc[i,2], data.iloc[i,4])
        tdoa2 = TDOA(data.iloc[i,1], data.iloc[i,3], data.iloc[i,3], data.iloc[i,5])
        tdoa3 = TDOA(data.iloc[i,2], data.iloc[i,4], data.iloc[i,3], data.iloc[i,5])

        # 根据三个TDOA计算位置坐标
        x = (tdoa1 * toa2 * toa3) / (tdoa1 + tdoa2 + tdoa3)
        y = (tdoa2 * toa1 * toa3) / (tdoa1 + tdoa2 + tdoa3)

    return [x, y]

# 调用函数计算位置坐标
position = calculate_position(data)

# 打印结果
print("飞行器位置坐标为：[%.2f, %.2f]" % (position[0], position[1]))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

问题2

问题2的模型建立过程：

飞行器实时位置估计模型：
根据题目信息，可以将飞行器的位置估计问题转化为时间序列分析的问题。设飞行器在时刻 $t$ 的位置为 $(x (t), y (t), z (t))$ ，则根据题目信息，可以建立如下的时间序列模型：
$x(t) = f_x(t) + e_x(t)$
$y(t) = f_y(t) + e_y(t)$
$z(t) = f_z(t) + e_z(t)$
其中， $f_x(t),f_y(t),f_z(t)$ 为飞行器的位置函数， $e_x(t),e_y(t),e_z(t)$ 为误差项。根据题目信息，误差项可以分为随机性偏差和常值飘移两种。因此，可以将误差项表示为：
$e_x(t) = \sigma_x(t) + \mu_x(t)$
$e_y(t) = \sigma_y(t) + \mu_y(t)$
$e_z(t) = \sigma_z(t) + \mu_z(t)$
其中， $\sigma_x(t),\sigma_y(t),\sigma_z(t)$ 为随机性偏差， $\mu_x(t),\mu_y(t),\mu_z(t)$ 为常值飘移。根据题目要求，需要建立飞行器实时位置的估计方法，因此可以利用滤波器的方法来对飞行器位置进行估计。

根据题目要求，需要建立飞行器的三维空间位置估计模型，可以使用三角定位法来估计飞行器的位置。假设飞行器位置为 $(x, y, z)$ ，发射源 $E_i$ 的位置为 $x_i,y_i,z_i)$ ，接收源 $R$ 的位置为 $x_r,y_r,z_r)$ ，则可得到如下方程：

{\begin{cases} (x - x_{i})^{2} + (y - y_{i})^{2} + (z - z_{i})^{2} = (v T O A_{i})^{2} \\ (x - x_{j})^{2} + (y - y_{j})^{2} + (z - z_{j})^{2} = (v T O A_{j})^{2} \end{cases}

$\begin{cases} (x-x_i)^2+(y-y_i)^2+(z-z_i)^2 = (vTOA_i)^2\\ (x-x_j)^2+(y-y_j)^2+(z-z_j)^2 = (vTOA_j)^2 \end{cases}$

{(x - x_{i})^{2} + (y - y_{i})^{2} + (z - z_{i})^{2} = (v TO A_{i})^{2} (x - x_{j})^{2} + (y - y_{j})^{2} + (z - z_{j})^{2} = (v TO A_{j})^{2}

求解该方程即可得到飞行器的实时位置。

python示例代码如下：

import numpy as np
import matplotlib.pyplot as plt

# 数据导入
data = np.loadtxt('附件1接收情况1数据.txt')
# 飞行器坐标
x = data[:, 1]
y = data[:, 2]
z = data[:, 3]
# 发射源坐标
xs = data[:, 4]
ys = data[:, 5]
zs = data[:, 6]
# 接收时间
receive_time = data[:, 7]
# 发射时间
transmit_time = data[:, 8]

# 计算信号传输时间
transmission_time = receive_time - transmit_time

# 计算到达角度信息
# 定义角度信息计算函数
def cal_angle(x, y, z, xs, ys, zs):
    dx = x - xs
    dy = y - ys
    dz = z - zs
    angle = np.arctan2(dy, dx)
    return angle

# 计算到达角度信息
angle = cal_angle(x, y, z, xs, ys, zs)

# 定义飞行器位置估计函数
def estimate_position(transmission_time, angle):
    # 速度光速
    c = 3 * 10**8
    # 发射源编号
    n = len(transmission_time)
    # 构建A矩阵
    A = np.zeros((n-1, 4))
    for i in range(n-1):
        A[i, 0] = c * (transmission_time[i+1] - transmission_time[0])
        A[i, 1] = c * (transmission_time[i+1] - transmission_time[0] - transmission_time[1] + transmission_time[0])
        A[i, 2] = c * (transmission_time[i+1] - transmission_time[0] - transmission_time[2] + transmission_time[0])
        A[i, 3] = c * (transmission_time[i+1] - transmission_time[0] - transmission_time[3] + transmission_time[0])
    # 构建b矩阵
    b = np.zeros((n-1, 1))
    for i in range(n-1):
        b[i, 0] = c * (angle[i+1] - angle[0])
    # 计算估计位置
    position = np.linalg.lstsq(A, b, rcond=None)[0]
    # 返回结果
    return position

# 计算位置估计值
position_estimate = estimate_position(transmission_time, angle)
# 提取x,y,z坐标估计值
x_estimate = position_estimate[:, 0]
y_estimate = position_estimate[:, 1]
z_estimate = position_estimate[:, 2]

# 绘制结果图
plt.figure()
plt.plot(x, y, label='Real Position')
plt.plot(x_estimate, y_estimate, label='Estimate Position')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

# 打印0秒至10秒的位置估计值
print('0秒至10秒的位置估计值为：')
print(x_estimate, y_estimate, z_estimate)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

查看完整思路详见：
【腾讯文档】第九届“数维杯”大学生数学建模挑战赛全题目深度剖析建模完整过程+详细思路+代码全解析（持续更新中！）
https://docs.qq.com/doc/DSHd6UHpNZkZoWmxw

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/从前慢现在也慢/article/detail/557381