2024年数维杯国际大学生数学建模挑战赛
C题 时间信号脉冲定时噪声抑制与大气时延抑制模型
原题再现：
脉冲星是一种快速旋转的中子星，具有连续稳定的旋转，因此被称为“宇宙灯塔”。脉冲星的空间观测在深空航天器导航和时间标准维护中发挥着至关重要的作用。
脉冲星时间在原子计时中的应用有望提高本地原子钟的稳定性和可靠性，代表着未来计时发展的一个长期方向。脉冲星时间研究的关键挑战之一是如何解决脉冲星定时噪声导致的精度和稳定性下降的问题。
脉冲星定时噪声是脉冲星旋转参数在长时间尺度（通常是数月或数年）内发生的连续扰动。它表现为预测脉冲到达时间（PT）和实际到达时间（PT-TT）之间的差异，该差异从不等于零。定时噪声通常是“红噪声”，几乎存在于所有脉冲星中，包括毫秒脉冲星。一些表现出随机变化，而另一些表现出准周期性，如图1所示。
脉冲星定时噪声主要包括脉冲星旋转产生的红噪声、色散测量（DM）变化引起的噪声、观测设备噪声和特定频带的噪声。IPTA发布的数据使用功率谱模型描述了红噪声强度，如下所示：
其中P0表示红噪声的强度，f是傅里叶频率，fc是角频率，q是频谱指数。除了功率谱估计方法（Matsakisd，1997）外，用于脉冲星定时噪声估计的方法包括Δ8 模型（Arzoumannian，1994；杨廷高，2014），指数模型（Shannon，2010），σz(i) 估计（Reardon，2016）和经验模态分解（EMD）方法（Huang，1998；Gao Feng，2018）。然而，这些方法只能在一定程度上提高脉冲星模型的精度。脉冲星定时噪声的来源广泛而复杂，因此寻找有效的处理方法是一项挑战。这仍然是精对苯二甲酸研究的一项重要任务。一些学者甚至尝试使用人工智能方法来提取和建模时序噪声的特征（梁洪涛，2023），旨在更好地解决时序噪声的去除和预测。
脉冲星时间可以从原子时间中创建独立的时间尺度，并为低地球轨道、地球同步轨道、高椭圆地球轨道、月球轨道、星际导航和深空导航中的航天器提供丰富的导航信息，如位置、速度和时间。脉冲星时间的主要考虑因素是脉冲到达时间（TOA），它受到各种延迟效应的影响。因此，延迟消除是决定脉冲星时间精度的关键因素。通常，观测到的脉冲到达时间（∆t）校正为太阳系重心。这种修正取决于脉冲星的位置、速度、质量和太阳系的天体。该校正的方程式可以总结为：
式中，∆c表示时钟延迟，∆A表示大气延迟，ΔEo表示爱因斯坦延迟，△Ro表示Romer延迟，△So表示Shapiro延迟，D/f^2表示色散延迟，∆·vp表示视差运动延迟，∆-B表示二元轨道运动延迟。
对于大气延迟，电磁波在大气中的传播速度比在真空中慢。在典型的观测频率下，总电子含量的显著变化可以导致10ns到几百ns的OA波动（Liu，2020）。
图3是大气延迟的图示，其中O代表地球中心，h₀ 是地球的表面高度。电磁波的路径来自h₀ to h₂,通过不同的大气层：平流层（h₁ and h₂) 和对流层（h₀ 和h1）。P-是电磁波的路径。角ε是通过电离层、平流层和对流层的电磁波的仰角，直到它到达地面望远镜（Liu，2020）。常用的Saastamoinen（1972）模型仅考虑平流层和对流层的折射时间延迟。折射延迟的方程为：
折射时间延迟仅对低于20GHz的无线电频率有效。这有助于减少色散效应，并提高超宽带和高频观测的脉冲星时间精度。大气时间延迟需要更好的建模，特别是对于仰角较小（10度或更小）的观测，其中映射函数的不准确可能会显著影响TOA（到达时间）精度（Liu，2020）。
请开发一个模型来解决以下问题：
问题（1）：考虑使用功能模型模拟图2中的脉冲星定时噪声，目标是95%或更高的模型拟合。建模所需的数据见附件1。可以参考和不使用的数据关系包括：脉冲星的观测频率在带宽为320 MHz的无线电频带中为1540 MHz，MJD 52473到56081的RMS值为75268.376µs，而MJD 52 473到56 646的RMS价值为78502.322µs。通常假设红噪声的强度与RMS值成比例，尽管不相等。
问题（2）：考虑对图2中脉冲星定时噪声的未来趋势进行短期（从几天到一个月）和长期（从几个月到几年）预测。预测验证所需的数据见附件1。
问题（3）：考虑为20 GHz以上的无线电观测频率建模折射时间延迟，确保天顶延迟小于或等于7.69 ns。相关参数和计算过程也可参考《时空参考系》第6.1章。
问题（4）：考虑为具有较小仰角（10度或更小）的观测建模大气延迟，以提高TOA精度。请提供您的模型，并描述您的考虑因素和可实现的目标。
整体求解过程概述(摘要)
本研究旨在解决脉冲星定时中噪声去除和大气延迟校正的挑战，以提高时间信号的精度。该研究为脉冲星在精密计时和深空导航中的未来应用奠定了科学基础。
对于问题1，任务涉及对脉冲星的红噪声进行建模。作为响应，我们的团队开发了一种基于功率谱密度（PSD）方法的红噪声产生和参数拟合模型。首先，使用FFT将PT-TT时间序列数据从时域转换到频域，然后计算PSD。然后，我们构建了一个基于PSD的红噪声模型，并应用最小二乘法来最小化模型预测和实际观测之间的误差。拟合过程强调了低频范围，这最好地捕捉了红噪声的特征。采用随机模拟方法产生红噪声，并用两个时间段的观测数据验证了结果。R²值为0.9507896654和0.9587555241的结果满足95%以上的目标精度。
对于问题2，该任务需要预测脉冲星定时噪声的短期和长期未来趋势。我们的团队开发了SARIMA模型和LSTM模型来解决这些趋势。根据文献和数据观测，脉冲星信号表现出周期性波动，使得SARIMA模型适合于短期预测，而长期趋势的存在需要使用LSTM模型。为了评估这些模型的性能，我们比较了SARIMA和ARIMA的MSE和RMSE度量，以及其他模型，如LSTM、GRU、MLP和随机森林。我们的比较表明，SARIMA模型和LSTM模型表现最好，R²值分别为0.96751486047和0.9936738159777331。此外，LSTM超参数的敏感性分析证实了该模型的鲁棒性，在各种设置下，R²保持在95%以上。
对于问题3，任务是对20 GHz以上无线电频率的折射时间延迟进行建模。为了解决这个问题，我们的团队开发了一个高频大气延迟校正模型。对于电离层延迟，我们使用了基于电子密度和信号频率的色散模型，而对于对流层延迟，则使用了改进的Saastamoinen模型。通过改变TEC值和频率范围来进行灵敏度分析，以评估天顶延迟变化。结果表明，该模型在不同条件下保持了良好的稳定性和准确性。
对于问题4，重点是低仰角下的精确大气延迟测量。该团队开发了一个用于小仰角的延迟映射模型。电离层延迟使用问题3中的相同色散模型进行建模，而对流层延迟分为干分量和湿分量，使用改进的Saastamoinen方法进行建模。我们为低仰角优化了Herring映射函数。通过改变仰角和TEC值进行敏感性分析。结果表明，在高频下，大气延迟效应减小，增加仰角显著减少延迟，表明使用更高频率的信号和增加仰角可以有效地提高测量精度和信号稳定性。
模型假设：
为了简化模型，做出了以下假设：
1.假设大气层是均匀的，折射率随高度连续变化。
2.高频无线电信号的延迟主要是由对流层效应引起的，没有考虑电离层延迟。
3.假设大气压力、温度和湿度遵循标准大气模型。
4.假设大气中的无线电信号传播路径是平滑和连续的，而不考虑湍流。
5.大气折射率被建模为温度、湿度和压力的线性函数。
6.在高频（20 GHz以上）下，电离层的影响可以忽略不计。
问题分析：
问题一的分析
问题一涉及建立一个模型来模拟脉冲星定时噪声的特性，目标是实现超过95%的模型拟合精度。然后，该模型将用于预测未来的噪声趋势，并支持脉冲星时间信号的改进。首先，我们的团队对附件1.xlsx中的数据进行了预处理，提取了MJD和PT-TT数据，以确保均匀采样。接下来，使用傅里叶变换计算功率谱密度（PSD），揭示噪声在频域中的分布特征。随后，基于PSD建立了红噪声模型，并通过最小化模型预测值和观测值之间的误差来应用最小二乘拟合。特别注意低频PSD特性，因为该区域的波动最有效地反映了红噪声的特性。最后，利用模型结果，我们的团队应用随机模拟来产生红噪声，并进行了几次验证实验，以确保模型在模拟脉冲星噪声时的可靠性。
问题二的分析
问题二旨在预测脉冲星定时噪声的短期和长期。脉冲星定时噪声具有周期性、长期趋势和随机波动，这就需要使用多种建模技术来准确捕获这些特征。为了有效地预测短期波动和长期趋势，我们的团队考虑了传统的统计特性和更复杂的非线性特征。因此，我们选择SARIMA时间序列模型来捕获周期性和趋势变化，并采用LSTM深度学习模型来增强对复杂非线性特征的预测能力，利用递归神经网络（RNN）来捕获长期相关性和复杂模式。通过结合统计模型和深度学习模型，我们全面地解决了噪声特性，并实现了更准确的预测。
问题三的分析
问题三需要对20 GHz以上频率的大气折射延迟进行建模，约束是折射延迟不得超过7.69纳秒。大气折射延迟主要由电离层和对流层引起。电离层延迟取决于射频和总电子含量（TEC），而对流层延迟受气压、温度、湿度和天顶角等因素的影响。为了降低模型复杂度，我们使用简单的色散延迟模型来模拟电离层延迟，而对流层延迟被建模为计算总折射延迟的基础。对于对流层折射延迟，我们使用改进的Saastamoinen模型，结合气压、温度和湿度，以确保延迟保持在指定的范围内。
问题四的分析
问题四侧重于在小仰角（10度或更低）下建模大气延迟，以提高脉冲星到达时间（TOA）精度。在这些角度下，信号通过大气的路径增加，导致折射诱导延迟显著增加。准确的建模需要考虑大气折射如何随高程变化。在标准折射延迟模型中，映射函数将天顶延迟转换为任何高程的延迟。我们的团队使用Herring映射函数，但在小角度下，它可能会由于大值而导致数值不稳定。为了解决这个问题，我们通过在分母中添加校正因子来修改函数，确保稳定的计算。此外，通过分离干延迟（电离层）和湿延迟（对流层），我们为这两个层开发了不同的大气延迟模型。
模型的建立与求解整体论文缩略图
全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# Read the Excel data
file_path = 'Attachment 1.xlsx'
data = pd.read_excel(file_path, header=None, names=["MJD", "PT-TT"])
# Print the structure of the data to ensure correct reading
print(data.head())
# Ensure the PT-TT column is of numeric type and remove any missing or anomalous
values
data['PT-TT'] = pd.to_numeric(data['PT-TT'], errors='coerce') # Convert PT-TT
column to numeric, invalid values become NaN
data = data.dropna(subset=['PT-TT']) # Drop rows containing NaN values
# Check the data again
print(data.head())
# Select the PT-TT data column for smoothing
pt_tt = data['PT-TT'].values
# Set the window size, using 10 days as the window size
window_size = 10
# Calculate the simple moving average
sma = np.convolve(pt_tt, np.ones(window_size)/window_size, mode='valid')
# Create a new DataFrame to save the results
result = pd.DataFrame({ "MJD": data['MJD'][window_size-1:].values, # Align the MJD with the
moving average "PT-TT (Moving Average)": sma # Moving average results
})
# Save the results to an Excel file
output_file_path = 'Attachment 2.xlsx'
result.to_excel(output_file_path, index=False)
print(f"Moving average results have been saved to the file: {output_file_path}")
# Visualize the results
plt.figure(figsize=(10, 6))
plt.plot(data['MJD'], pt_tt, label='Original PT-TT (s)', color='blue', alpha=0.6)
plt.plot(data['MJD'][window_size-1:], sma, label=f'Moving Average (Window
size={window_size})', color='red', linewidth=2)
plt.xlabel('MJD (days)')
plt.ylabel('PT-TT (s)')
plt.title('Pulsar Timing Noise with Moving Average')
plt.legend()
plt.grid(True)
plt.show()
import pandas as pd
import numpy as np
from scipy.fft import fft, fftfreq
import matplotlib.pyplot as plt
# Step 1: Load data and extract relevant columns
data = pd.read_excel('Attachment 1.xlsx') # Adjust based on the data file path
# Clean the data
data_cleaned = data.iloc[1:, 1:3] # Skip the header row, select MJD and PT-TT
columns
data_cleaned.columns = ['MJD(days)', 'PT-TT(s)'] # Rename columns
data_cleaned['MJD(days)'] = pd.to_numeric(data_cleaned['MJD(days)'],
errors='coerce')
data_cleaned['PT-TT(s)'] = pd.to_numeric(data_cleaned['PT-TT(s)'],
errors='coerce')
data_cleaned.dropna(inplace=True) # Drop missing values
# Extract time and signal data
time = data_cleaned['MJD(days)'].values
signal = data_cleaned['PT-TT(s)'].values
# Step 2: Calculate time interval
# Since the Fourier transform requires a uniformly sampled signal, we
approximate the time interval by the average difference between adjacent time points
delta_t = np.mean(np.diff(time))
# Step 3: Perform Fourier Transform
# Use the fft function to perform a Fast Fourier Transform, converting the
time-domain signal to the frequency domain
N = len(signal) # Number of data points in the signal
signal_fft = fft(signal) # Perform Fourier transform on the signal
frequencies = fftfreq(N, delta_t) # Calculate corresponding frequencies
# Step 4: Calculate Power Spectral Density
# PSD (Power Spectral Density) is the energy intensity at each frequency,
calculated from the Fourier transform results
psd = np.abs(signal_fft)**2 / N # Take the square of the absolute value of the
Fourier transform result, then divide by the number of data points
# Retain only the positive frequencies
positive_frequencies = frequencies[frequencies > 0]
positive_psd = psd[frequencies > 0]
# Step 5: Plot Power Spectral Density
plt.figure(figsize=(10, 6))
plt.loglog(positive_frequencies, positive_psd, color='orange')
plt.xlabel('Frequency (1/days)')
plt.ylabel('Power Spectral Density')
plt.title('Power Spectral Density of Pulsar Timing Noise')
plt.grid(True)
plt.show()
