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

问题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]))

问题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{array}{l}(x-x_i{)}^2+(y-y_i{)}^2+(z-z_i{)}^2=(vTO{A}_i{)}^2\\ (x-x_j{)}^2+(y-y_j{)}^2+(z-z_j{)}^2=(vTO{A}_j{)}^2\end{array}$$

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

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)

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