2024年数维杯国际大学生数学建模挑战赛A题飞行器激光测速中的频率估计问题解题全过程论文及程序
2024年数维杯国际大学生数学建模挑战赛
A题 复合直升机的建模与优化控制问题
原题再现：
（一） 问题的背景
空速，即飞机相对于空气的速度，是飞行期间需要监控的关键参数。空速与飞行状态密切相关，如迎角和侧滑角。如果空速数据异常，很容易导致失速等事故。因此，准确测量空速具有重要意义。
图1：激光测速空速原理图
激光测速是一种可行的空速测量方法。如图1所示，该原理涉及以固定频率发射激光，然后通过空气中气溶胶粒子的Mie散射效应获得多普勒频移的信号光。利用相干干涉原理，获得了包含多普勒频移信息的信号，并估计了该信号的频率。利用估计的多普勒频移信息，可以计算空速。空速测量的关键步骤是估计时间序列信号的频率信息，其中信号满足以下表达式：
其中Are表示信号的幅度，f0是信号的频率，φ是信号的相位，z（t）表示噪声信息。由于来自空气压力和温度等因素的干扰，飞机接收的信号与大量噪声信息z（t）混合，使得信号x（t）的频率估计问题非常具有挑战性。
目前，存在一个采样间隔为Ts=2×10^（-9）秒的飞机在太空中飞行的示例，接收多个多普勒频移信号，实际接收的数据提供在附件1中。为确保信息安全，不同的飞行周期可能导致飞机接收到的激光信号的幅度、频率和相位的变化，并且环境噪声可能不同。因此，接收到的激光信号的特征，包括频率和相位，在不同的飞行周期中可能不相同。此外，环境噪声也不同，即附件1中的噪声特性在不同的子表中不相同。
（二） 问题的数据
请参阅附件1。
（三） 需要解决的问题
1.分析实际接收信号的噪声特性可以帮助设计信号频率估计算法。在附件1的飞行周期1中，接收信号的非噪声部分的已知幅度为4，频率为30×10^6 Hz，相位为45°。请分析飞行周期1中接收数据的噪声z（t）特性。
2.在实际场景中，接收信号的非噪声部分的频率未知，需要估计。在附件1的飞行周期2中，已知实际接收信号的幅度为2，并且相位为0°。请设计一种方法来估计飞行周期2中接收信号的非噪声部分的频率。（注：飞行周期2的噪声特性可能与飞行周期1的噪声特性不同。）
3.在实践中，通常不可能预先知道接收信号的非噪声部分的幅度和相位信息，但仍然需要估计频率。请根据附件1中来自飞行周期3的数据，设计一种方法来估计飞行周期3中接收信号的频率。（注：飞行周期3中的噪声特性可能与飞行周期1和2中的噪声特征不同。）
4.为了避免信号之间的干扰，在实践中使用间歇接收，限制可用信息的量。参考附件1中飞行周期4的数据。基于该数据，请分析间歇接收的模式，并设计一种方法来估计飞行周期4中接收信号的频率。
整体求解过程概述(摘要)
对于问题1，区分混合变量数据与建模统计特征时使用的计算公式在数值上的差异。计算数据的基本统计量，如平均值、最大值和最小值、方差、偏度和峰度，以描述数据的分布和形状，计算傅里叶变换以实现频谱分析，并通过计算异常幅度波动和时间不规则性来进一步分析数据。发现噪声序列是一个独立的序列，服从正态分布，平均值为0.00584，方差为3.9904。此外，噪声序列在很大程度上逼近白噪声数据，并且具有白噪声数据的特征。
对于问题二，首先处理包含接收信号的时间和值的时间序列数据。使用孤立森林算法来识别数据中的离群点。为了深入分析数据并检测可能的小位移或趋势变化，进一步引入滑动窗口技术和累积和（CUSUM）方法来检测频率周期。最后，引入隐马尔可夫模型来估计非噪声部分的频率，发现FFT结果中最大幅度的频率为41×106Hz。
对于问题三，本文采用FFT、PSD、MUSIC算法和自相关函数法对周期3信号数据进行频率估计。通过去除直流分量、平滑和滤波，提高了信号质量，信噪比达到20.69dB。分析表明，该信号的主频约为35MHz，其中FFT、PSD和MUSIC算法估计分别为34.99986MHz、35.16MHz和35.4MHz。所有方法的结果基本一致，验证了频率估计的可靠性，为高精度信号分析提供了参考。
对于问题四，本文提出了两种估计第四飞行周期内接收到的间歇信号频率的方法。第一种方法是卡尔曼滤波器，它利用递归算法来预测和更新信号，最终确定频率为56.29MHz。利用FFT谱图进一步验证了该结果，证明了其可靠性。第二种方法是小波变换，它从中心频率计算主频，以获得间歇信号的频率为55.56MHz。这两种方法都表明频率估计是可靠的。
模型假设：
1.过程噪声符合零均值高斯分布，其中在状态转变期间引入的噪声是随机的，并且其分布遵循平均值为零的高斯（正态）分布。
2.观测噪声符合零均值高斯分布。类似地，在观测过程中引入的噪声也是随机的，其分布也遵循均值为零的高斯分布。
3.通过对一组小波基函数进行尺度和位移变换，得到第四个飞行信号。小波变换通过对原始信号进行连续细分，并选择小波基进行多尺度和位移变换来实现这一点，从而获得一组不同的函数。
问题分析：
问题一的分析
在附录1的飞行循环1中，已知接收信号的非噪声部分具有4的幅度、30×106 Hz的频率和45的相位◦. 请分析飞行循环1接收数据的噪声z（t）特性。针对这一问题，我们首先对数据进行预处理，将数据序列分解为噪声序列和纯信号序列，分别对原始序列、噪声序列和纯粹信号序列进行分析后，对这三个变量的值进行统计表征，并通过概率密度函数获得噪声序列的特征。
问题二的分析
在附件1的飞行循环2中，已知实际接收的信号具有振幅2和相位0◦ . 噪声序列可以被假设为白噪声和非白噪声两种情况，并且通过识别离群点、检测频率周期和构建隐马尔可夫模型来估计非噪声部分的频率来估计这两种情况中的每一种情况的频率。
问题三分析
本文采用各种频率估计方法，包括FFT、PSD、MUSIC和自相关函数，来分析周期3信号数据。预处理步骤（如DC去除、平滑和滤波）增强了信号的质量，其特征是高采样率（500MHz）和强SNR（20.69dB）。频率分析揭示了35MHz左右的一致主频估计（例如，FFT:34.99986MHz，PSD:35.16MHz，MUSIC:35.4MHz），尽管存在微小差异，但确认了方法的可靠性。该研究强调了这些技术在准确识别信号主频方面的有效性。
问题四分析
第四个飞行周期采用间歇接收模式。该模式的挑战是接收的信号数据可能包含间隙，导致信息丢失并影响频率估计的精度。为了解决这个问题，有必要分析信号的采样间隔，识别数据间隙，并设计一种合理的频率估计方法。
模型的建立与求解整体论文缩略图
全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：
import pandas as pdimport matplotlib.pyplot as pltfrom scipy.statsi mport gaussian_kde,kurtosis,skew,shapirofrom scipy.fft import fft,fftfreqfrom stats models.tsa.stattools import acf,pacf,adfuller,q_statfrom sklearn.ensemble import IsolationForestfrom hmmlearn import hmmplt.rcParams[’font.sans-serif’] = ’SimHei’plt.rcParams[’axes.unicode_minus’]=False#Filepathfile_path=’C:/Users/25833/Desktop/noise_data.xlsx’data=pd.read_excel(file_path)#Extractingcolumnstimes=data[’Time’]noises=data[’Noise’]xt=data[’X(t)’]signal=data[’Signal’]#%%Plotnoiseovertimefig=plt.figure(figsize=(10,6), facecolor=’white’,dpi=150)ax=fig.add_subplot(111,facecolor=’gray’)ax.plot(times,noises,linestyle=’-’,color=’white’)ax.set_xlabel(’Time’,color=’black’,fontsize=24)ax.set_ylabel(’Noise’,color=’black’,fontsize=24)ax.spines[’bottom’].set_color(’black’)ax.spines[’left’].set_color(’black’)forspineinax.spines.values():spine.set_linewidth(2)plt.margins(x=0,y=0)ax.tick_params(axis=’both’,which=’major’,labelcolor=’black’,colors=’black’,labelsize=15)ax.grid(True,color=’lightgray’,linestyle=’--’,linewidth=0.5)plt.show()
#%%KDEforX(t)noise_values=xtkde=gaussian_kde(noise_values,bw_method=’scott’)x=np.linspace(min(noise_values),max(noise_values),1000)pdf_values =kde(x)plt.figure(figsize=(10,6),dpi=150)plt.plot(x,pdf_values,color=’green’,lw=2,label=’X(t)PDF’)plt.fill_between(x,pdf_values, color=’green’,alpha=0.3)plt.title(’ProbabilityDensityFunctionofX(t)’)plt.xlabel(’X(t)Value’,fontsize=24)plt.ylabel(’Density’,fontsize=24)plt.legend()plt.show()#%%Signal statisticsdf=pd.DataFrame(data)mean_value =df[’Signal’].mean()median_value=df[’Signal’].median()std_dev=df[’Signal’].std()range_value=df[’Signal’].max()-df[’Signal’].min()cv=std_dev/mean_valueskewness= skew(df[’Signal’])kurtosis_value=kurtosis(df[’Signal’])-3print(f"Mean:{mean_value}")print(f"Median:{median_value}")print(f"StandardDeviation:{std_dev}")print(f"Range:{range_value}")print(f"CoefficientofVariation: {cv}")print(f"Skewness:{skewness}")print(f"Kurtosis:{kurtosis_value}")#%%ACFandPACFplotsfornoiselag_acf=acf(noises,nlags=20)lag_pacf= pacf(noises,nlags=20, method=’ols’)plt.figure(figsize=(12,6),dpi=150)plt.subplot(211)plt.plot(lag_acf)plt.title(’AutocorrelationFunction(ACF)ofNoise’)plt.subplot(212)plt.plot(lag_pacf)plt.title(’PartialAutocorrelationFunction(PACF)ofNoise’)plt.tight_layout()plt.show()#%%FFTof noisefft_noise=fft(noises)frequencies=fftfreq(len(noises))plt.figure(figsize=(12,6),dpi=150)plt.plot(frequencies,np.abs(fft_noise))plt.title(’FrequencySpectrumofNoise’)plt.xlabel(’Frequency(Hz)’)plt.ylabel(’Amplitude’)plt.tight_layout()plt.show()#%%Correlationheatmapcorrelation_matrix=np.corrcoef(signal,noises)plt.figure(figsize=(6,6),dpi=150)plt.imshow(correlation_matrix,cmap=’hot’,interpolation=’nearest’)plt.colorbar()plt.title(’CorrelationHeatmapbetweenSignalandNoise’)plt.xticks([])plt.yticks([])plt.show()#%%Ljung-Boxtestlb_value,p_value=q_stat(lag_acf,len(noises))plt.figure(figsize=(12,6),dpi=150)plt.plot(lb_value[1:],p_value[1:],’o-’)plt.title(’Ljung-BoxTestforNoise’)plt.xlabel(’Lag’)plt.ylabel(’p-value’)plt.axhline(y=0.05,color=’red’,linestyle=’--’)plt.show()#%%ADFtestresult=adfuller(noises)