2024年数维杯数学建模A题多源机会信号建模与导航分析解题全过程论文及程序

2024年数维杯 数学建模

A题多源机会信号建模与导航分析

原题再现：

（一）问题背景
尽管全球卫星定位系统下的定位导航技术已成熟，但考虑到室内、隧道、建筑密集区等复杂环境或全球卫星定位系统被毁失灵等突发场景，会发生全球卫星定位系统拒止情况，无法有效定位导航。因此，需要发展基于新型信号的自主定位导航方法，机会信号导航是目前一种可行的自主导航技术。机会信号是指存在于空间域中的各类无线电信号，具有不同的频段。在多个发射源发射多种机会信号（机会信号中均带有发射源位置信息、发射时间信息）的场景下，考虑接收飞行器的自主导航解算问题，即：飞行器根据接收到的机会信号信息，通过飞行器机载设备，实时计算飞行器自身的三维空间位置。
根据机会信号蕴含的信息来分类，可以将机会信号分为五类：达到时间信息（TOA）、到达时间差信息（TDOA）、多普勒频率差信息（DFD）、到达角度信息（AOA）、接收强度指标信息（RSSI）。五类机会信号的具体信息含义如下：
（1）达到时间信息TOA：信号传输时间，即飞行器接收信号时间与信号发射时间的差。
（2）到达时间差信息TDOA：同一信号从两个发射源（同时发射）到达接收端的时间差。
（3）多普勒频率差信息DFD：同一信号从两个发射源发射，由于飞行器与信号源具有相对速度，接收信号会产生频率变化，进而产生信号频率差，具体计算方式如下
多普勒频率差信息=发射信号的频率/信号传播速度∙ (发射源1对应的频率变化比率−发射源2对应的频率变化比率)
其中发射源对应的频率变化比率=发射源相对接收源的速度向量⋅发射源相对接收源的位移向量/发射源与接收源的距离
（4）到达角度信息AOA：接收源可得到发射源信号的相对角度信息，如图1所示（发射源与接收源连线在xOy平面投影线段与x轴正方向的夹角，发射源与接收源连线与z轴负方向的夹角β）。为了简化处理，图1中的坐标系为地面惯性系，到达角度信息为角度的正切值，即tanα/tanβ。
在这里插入图片描述
（5）接收强度指标信息RSSI：通过对比标称距离下的标称信号强度，可以获得接收信号的强度指标信息，具体计算方式如下接收强度指标信息=标称信号强度−10∙信道衰减系数∙ lg(发射源与接收源相对距离/标称距离)现有一例飞行器在空中做动态飞行，以0.01秒的采样间隔时间接收到不同种类的机会信号，发射源数据、机会信号参数和接收到的机会信号具体数据见附件1。
（二）题目数据参见附件1。
（三）需要解决的问题
1. 结合题目信息与数据，建立机会信号的数学表达式。同时，针对每一类机会信号，讨论能够唯一确定飞行器位置的最少的机会信号个数。
2. 根据附件1的接收情况1数据，不考虑数据偏差的情况下（认为所有数据可信），请设计飞行器实时位置的估计方法，并给出飞行器0秒至10秒的导航定位结果（实时的位置估计值）。
3. 在附件1的接收情况1数据中，某些机会信号可能有较大的偏差，请建立机会信号的实时筛选方法，筛选出偏差较大的机会信号。根据建立的机会信号筛选方法，给出此时飞行器0秒至10秒的导航定位情况。
4. 机会信号的偏差可以分为两种，一种是随机性偏差（零均值），一种是常值飘移。请建立评价判断方法，并依据所提出方法，判断接收情况2中的机会信号的随机性偏差程度以及常值飘移量。设计合理的机会信号筛选方法，给出接收情况2下的飞行器0秒至10秒的定位结果。

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

针对全球卫星定位系统（GNSS）在室内、隧道或城市密集环境中的局限性，提出了一种创新的飞行器定位系统，该系统采用牛顿迭代算法融合多种机会信号以实现飞行器的自主定位。系统设计的核心在于利用环境中的无线电信号，通过解析信号中蕴含的位置和时间信息，从而在GNSS信号受限的情况下提供精确的导航定位。
在问题一中，建立了机会信号的数学模型，并确定了最少信号个数以唯一确定飞行器的位置。模型基于信号处理和定位理论，综合考虑了TOA、TDOA、AOA、DFD和RSSI 信号，通过几何和代数方法分析，得出了在三维空间中定位所需的最少信号数量。其中，TOA、TDOA和DFD信号至少需要4个发射源，而AOA信号至少需要2个发射源，RSSI信号则至少需要4个发射源。在特殊情况下，如两发射源之间的距离为2D时TOA信号和RSSI信号最少可以用2个发射源确定飞行器的三维空间位置。
问题二中，设计了一种实时位置估计方法，通过数据预处理将AOA转换为角度值，RSSI 转换为信号强度值，TOA转换为距离估计，TDOA转换为距离差估计，DFD转换为频率差估计。随后，建立了基于三维坐标系的飞行器位置参数模型，采用牛顿迭代法求解非线性方程，以实现飞行器位置的实时估计。该方法利用发射源的初始位置坐标和飞行器的初始位置坐标，通过迭代逼近过程，快速找到满足精度要求的近似解，从而在无数据偏差的假设下，给出了飞行器在指定时间段内的导航定位结果。
问题三中，整合了TOA、TDOA、DFD、AOA和RSSI等多种信号信息，定位误差分析采用了统计定位误差和克拉美罗界（CRLB）两种方法，以评估定位系统的性能。利用卡尔曼滤波器算法处理数据并进行实时定位，通过预测和更新两个步骤，结合实时观测数据，提供连续且准确的定位结果。此外，设计了一种机会信号筛选模型，使用CRLB作为筛选标准，移除偏差较大的信号，以提高定位的准确性，并不断更新飞行器的位置估计。
问题四的分析中，建立了评价判断方法，以判断信号的随机性偏差和常值飘移，并设计了合理的机会信号筛选方法。通过假设检验和置信区间的概念来评估信号偏差的统计特性，设计了一种综合考虑偏差特性的信号筛选方法，以优化飞行器的定位结果。将牛顿迭代算法应用于机会信号的实时筛选和定位问题，提出了一种新的信号融合模型，并通过卡尔曼滤波器算法实现了实时定位。此外，还对定位误差进行了深入分析，包括统计定位误差和克拉美罗界，为飞行器定位系统的优化提供了理论依据。通过仿真实验，验证了所提方法的有效性，为提高飞行器的自主导航能力提供了一种新的解决方案。

问题分析：

问题一的分析
问题一要求我们建立机会信号的数学模型，并确定最少信号个数以唯一确定飞行器的位置。这一问题将基于信号处理和定位理论，考虑到每种信号类型提供的几何和物理信息。例如，TOA和TDOA信号可以用来确定飞行器与发射源之间的距离或距离差，而AOA信号提供了角度信息，DFD信号包含了频率变化信息，RSSI信号则提供了信号强度信息。我们将利用这些信息构建数学模型，并通过几何和代数方法分析最少信号个数，以实现三维空间中的定位。
问题二的分析
问题二要求设计一种实时位置估计方法，并在无数据偏差的假设下，给出飞行器在指定时间段内的导航定位结果。通过预处理信号数据并转换为可用于定位的数值，进而利用牛顿迭代法对基于不同信号类型的非线性方程进行求解，以实现对飞行器在特定时间段内连续位置变化的精确估计。
问题三的分析
问题三要求开发一种实时筛选方法，以识别并排除偏差较大的信号，提高飞行器的导航定位精度。关键点在于构建一个信号质量评估机制，并通过该机制实时监测和筛选信号。解决思路涉及信号融合模型的建立，定位误差的统计分析，以及卡尔曼滤波器算法的应用，通过这些步骤实现对偏差信号的有效识别和排除，优化定位结果。
问题四的分析
问题四要求建立评价判断方法，以判断信号的随机性偏差和常值飘移，并设计合理的机会信号筛选方法。这一问题的分析将涉及到概率论和数理统计的知识。我们将使用假设检验和置信区间的概念来评估信号偏差的统计特性。对于随机性偏差，可能会采用基于最小二乘法的线性回归模型；对于常值飘移，则可能采用时间序列分析中的自适应滤波技术。完成偏差评估后，将设计一种综合考虑偏差特性的信号筛选方法，以优化飞行器的定位结果。

模型假设：

1．信号覆盖假设：假设飞行器在其飞行范围内始终能够接收到至少三个发射源的机会信号。
2．信号接收假设：飞行器的接收设备能够无误差地接收到所有机会信号，并且接收时间与发射时间的同步是精确的。
3．数据处理能力假设：飞行器的机载设备具有足够的计算能力，能够实时处理所有接收到的信号，并进行快速的定位解算。
4．算法稳定性假设：所设计的算法在面对不同发射源布局和信号特性时，均能保持稳定性和有效性。
5．模型普适性假设：建立的数学模型和算法适用于所有类型的飞行器，包括但不限于无人机、飞机、卫星等。

论文缩略图：

在这里插入图片描述

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

部分程序代码：

import numpy as np
import pandas as pd
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#Constants
C=3e8 #Speedoflight in m/s
f0 = 3e8 #Frequency in Hz
k =20 #Channelattenuation coefficient
s =100 #Nominalsignal strength
nominal_distance = 1000 # Corresponding to nominal signal strength
#Initial positions and velocity vectors of the sources
source_positions = np.array([[100, 300, 50],[1000, 100, 150],
[-1000,-100, 450],
[-200,-100, 1050]])
source_velocities = np.array([[0, 0, 0],[0, 0, 0],[0, 0, 8],[20, 0, 0]])
#Load the data from Excel
file_path = 'D:/python/LDcode/database/飞行器接受情况 1.xlsx'
data = pd.read_excel(file_path)
#Extracting columns for different signals
times = data['时间 t（s）'].values
toa = data[['TOA(发射源 1)', 'TOA(发射源 4)']].values
tdoa = data[['TDOA(发射源 1,发射源 2)', 'TDOA(发射源 3,发射源4)']].values
dfd = data[['DFD(发射源 1,发射源 2)']].values
aoa = {'alpha': data['AOA[tan(alpha)](发射源 4)'].values,'beta': data['AOA[tan(beta)](发射源 4)'].values}
rssi = data[['RSSI(发射源 1)', 'RSSI(发射源 2)', 'RSSI(发射源 3)']].values
#Initialize aircraft position (starting guess)
initial_guess = np.array([0, 0, 0])
#Function to calculate the residuals for Newton's method
def residuals(position, toa, tdoa, dfd, rssi, aoa, times):x, y, z = positionres = []for t_idx, t in enumerate(times):toai = toa[t_idx]tdoaij = tdoa[t_idx]dfdij = dfd[t_idx]rssi_i = rssi[t_idx]aoa_alpha = aoa['alpha'][t_idx]aoa_beta = aoa['beta'][t_idx]for i in range(len(source_positions)):xi = source_positions[i, 0] + t * source_velocities[i, 0]yi = source_positions[i, 1] + t * source_velocities[i, 1]zi = source_positions[i, 2] + t * source_velocities[i, 2]di = np.sqrt((x- xi)**2 + (y- yi)**2 + (z- zi)**2)if i < len(toai) and not np.isnan(toai[i]):res.append(toai[i] * C- di)if i < len(rssi_i) and not np.isnan(rssi_i[i]):res.append(1000 * np.exp((s- rssi_i[i]) / (10 * k))- di)if i == 3:res.append(y- aoa_alpha * (x- xi))res.append(np.sqrt(di**2- (zi- z)**2)- aoa_beta * (zi- z))for j in range(i + 1, len(source_positions)):if j < len(tdoaij) and not np.isnan(tdoaij[j- 1]):xj = source_positions[j, 0] + t * source_velocities[j, 0]yj = source_positions[j, 1] + t * source_velocities[j, 1]zj = source_positions[j, 2] + t * source_velocities[j, 2]dj = np.sqrt((x- xj)**2 + (y- yj)**2 + (z- zj)**2)res.append(tdoaij[j- 1] * C- (di- dj))if j < len(dfdij) and not np.isnan(dfdij[j- 1]):vi_ri = np.dot(source_velocities[i], [x- xi, y- yi, z- zi]) / divj_rj = np.dot(source_velocities[j], [x- xj, y- yj, z- zj]) / djres.append(dfdij[j- 1]- (f0 / C) * (vi_ri- vj_rj))return res# Generate time intervals from 0 to 10 seconds with a step of 0.01 secondstime_intervals = np.arange(0, 10.01, 0.01)# Store the resultspositions = []# Perform the optimization for each time stepfor t in time_intervals:res_lsq = least_squares(residuals, initial_guess, args=(toa, tdoa, dfd, rssi, aoa, [t]))positions.append(res_lsq.x)initial_guess = res_lsq.x # Update initial guess for next iteration# Convert to numpy arraypositions = np.array(positions)# Print the positions at each time stepfor idx, (t, pos) in enumerate(zip(time_intervals, positions)):print(f"Time {t:.2f} s: Position (x, y, z) = {pos}")# Visualizationfig = plt.figure()ax = fig.add_subplot(111, projection='3d')ax.plot(positions[:, 0], positions[:, 1], positions[:, 2], label='Aircraft Trajectory')ax.set_xlabel('X')ax.set_ylabel('Y')ax.set_zlabel('Z')ax.legend()plt.show()