CN103926599A - 基于emd迭代阈值滤波的gnss多径效应抑制方法 - Google Patents

Abstract

本发明公开一种基于EMD迭代阈值滤波的GNSS多径效应抑制方法，对测量误差EMD分解得到N个本征模态函数分量，多径误差常集中在高阶IMF和余项中，一方面，采用奇异谱分析对EMD分解的N个IMF进行噪声估计，根据估计结果选择前M个IMF用作噪声辅助数据分析，得到测量误差的EMD分解结果；另一方面，采用随机采样第1个IMF数据的方式，构造与原信号具有相同信噪比的多个序列，采用模态单元阈值滤波对序列除噪，平均除噪后的多组结果分离得到GNSS多径误差，将其用于补偿GNSS测量结果即实现其测量过程中多径效应的抑制。本发明提高EMD分解精度，减小噪声的影响。

Description

基于EMD迭代阈值滤波的GNSS多径效应抑制方法

技术领域

本发明涉及全球卫星导航***，具体涉及一种基于EMD迭代阈值滤波的GNSS多径效应抑制方法。

背景技术

GNSS泛指卫星导航***，包括全球性、区域性和增强的卫星导航，例如美国的GPS、俄罗斯的Glonass、欧洲的Galileo以及中国的北斗卫星导航***。GNSS软件接收机通过接收卫星的广播星历信号，测量信号的传输延时，并结合信号的传播速度计算出接收机与卫星的距离来实现定位。由于接收机周围可能含有反射物，使得其接收到的信号是卫星信号与其反射信号的混合信号，信号反射增加了卫星信号的传输延时，所以产生距离计算误差，并最终影响定位结果，这种现象称为多径效应。影响卫星定位精度的误差有很多，近年来研究发现，基于短基线的双差分模型能有效的去除***误差中的具有相关特性的误差，如卫星钟差、接收机钟差以及电离层、对流层传输延迟等等，但是由于多径误差和接收机内部器件产生的热噪声不具有显著的相关性，所以无法通过双差分模型进行消除。

目前抑制多径误差的方法主要包括：空间选择法，即选择合适的接收机工作环境；采用先进的天线技术和接收机硬件技术以及各种运行于接收机内部的自适应滤波方法。基于硬件技术的多径抑制方法，不仅增加了***的体积和成本（例如增加扼流圈，虽然能抑制来自天线下方的大部分多径误差，但是对来自天线上方的多径误差效果有限），更重要的是其使用效果并不理想，进而限制了其在低成本特别是民用接收机中的应用。由于基于软件的实时滤波方法具有成本低，实施简单的特点，近年来受到了广泛关注。基于多径误差的时变特性，即其统计特征是随时间变化的，传统的固定结构的滤波器无法使用，所以许多自适应滤波方法被引入到多径误差的抑制过程中，例如基于小波变换的自适应滤波、基于变结构的有限冲激响应(FIR)滤波以及基于EMD的常规自适应滤波，上述方法均在一定的条件下取得了很好的效果，但是由于多径误差与接收机环境密切相关，动态环境下自适应滤波设计复杂，其适用性有待进一步研究。EMD是一种数据驱动的自适应分解方法，它基于任意信号均由一组本征模态函数组成的假设，将信号按不同的波动和趋势进行逐级分解，最终得到频率逐级降低的一组本征模态函数(IMF)。如果接收机工作于动态模式下其多径误差的变化具有一定的随机性，但是其幅值往往大于热燥声，频率小于热噪声。

发明内容

发明目的：为解决现有技术中存在的不足，本发明提供一种基于EMD迭代阈值滤波的GNSS多径效应抑制方法。

技术方案：本发明的一种基于EMD迭代阈值滤波的GNSS多径效应抑制方法，包括以下步骤：

（1）对GNSS短基线的双差分残差信号e(t)进行EMD分解，得到一组 $\mathrm{IMF}{\mathrm{\Σ}}_{i=1}^N{\mathrm{imf}}_i$ (t)和余项r(t)，即 $e\left(t\right)={\mathrm{\Σ}}_{I=1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$

其中，IMF是指本征模态函数，i为模态分量阶数，imf_i(t)为第i阶模态分量序列，N为分解得到的IMF总阶数，t为残差信号的时间变量，其取值与样本数据的长度有关系；

（2）利用奇异谱分析对步骤（1）中的IMF进行噪声估计，确定出噪声占优的前M个IMF分量；

（3）设Y(t)为噪声辅助分析的序列，初始时Y(t)=e(t)，求解imf(t)_i时，将Y(t)逐次加入辅助噪声经EMD分别得到1阶IMF，记结果集为IMFs，设共加入I次噪声，对得到的IMFs取平均，则得到所求IMF分量：

且1≤q≤I

Y(t)=Y(t)-imf_i(t)作为分解对象逐次加入噪声继续上述过程得到1阶IMF，依次得到imf₁～imf_M，记R(t)=Y(t)-imf_M为余项；其中，q是指第q次加入辅助噪声得到的IMF编号；

（4）对R(t)进行EMD分解得到信号的低频IMF分量imf_M+1～imf_N，从而实现含噪声信号的完整EMD分解，N为应用噪声辅助分析后EMD分解得到的IMF总阶数；

（5）基于步骤（4）中的IMF分量进行阈值滤波，滤波后重构的信号表示为 $\mathrm{er}\left(t\right)={\mathrm{\Σ}}_{i=n}^m{f}_{\mathrm{thr}}\left({\mathrm{imf}}_i\left(t\right)\right)+{\mathrm{\Σ}}_{i=m+1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$ 其中f_thr(imf_i(t))为阈值滤波函数，其中1≤n≤m-1，m、n分别为阈值滤波的上界和下界，m的选择为滤波提供了灵活性；

（6）将er(t)用于GNSS双差分测量结果X(t)的实时补偿数据，X(t)-er(t)即为去掉多径误差后的GNSS测量结果。

进一步的，所述步骤（2）中的噪声的具体估计步骤为：

（I）对任意长度为L的序列imf_i(t)，1≤t≤L，选择长度为m的时延窗口，构造一个m×n的轨迹矩阵X=[X₁,X₂,…,X_n]；

其中，X_k=(imf_i(k),imf_i(k+1),…,imf_i(k+m-1))^T，(·)^T表示矩阵的转置，n=L-m+1，1≤k≤n，2≤m≤L/2，上述过程可以视为在R^L空间中将时间序列imf_i(t)线性映射到向量空间构造出m×n的Hankel矩阵；奇异值分解(SVD)可以将复杂的矩阵分解成几个简单矩阵相乘的形式，具有明显的物理意义,其不要求待分析数据为方阵形式，较特征值分解的形式适用性更广；

（II）设λ₁≥λ₂≥…≥λ_m为XX^T的特征值，其对应的特征向量为u₁,u₂,…u_m，定义 $X=X_1+X_2+\·\·\·+X_m,X_i={\mathrm{\λ}}_i^{\frac{1}{2}}{u}_i{v}_i^T;$

其中 $v_i=X^T{u}_i/{\mathrm{\λ}}_i^{\frac{1}{2}},X^T$ 为X的转置矩阵；由于 ${\left|\right|X\left|\right|}^2=\mathrm{trace}\left\{{\mathrm{XX}}^T\right\}={\mathrm{\Σ}}_{i=1}^m{\mathrm{\λ}}_i,$ trace{·}表示矩阵的迹，X_i ²=λ_i(1≤i≤m)；

（III）定义为分量X_i对X能量的贡献，Ψ越小说明贡献越小，以含噪声信号的IMF分量为分析对象即可绘制多个分量序列的奇异谱。由于白噪声奇异值相差不大，所以其奇异向量对信号能量变化的贡献是稳定的，即白噪声的奇异谱是平坦的，而传感器输出的非平稳信号在延迟窗内其奇异值变化较大，能量变化多集中到前几个奇异值，其奇异谱有陡峭的下降趋势，基于上述现象可以对IMF分量进行噪声水平的定性评判。

进一步的，所述步骤（3）中加入的辅助噪声为白噪声，其满足（0，σ²）分布，σ的取值范围为0.01～0.5。

进一步的，所述步骤（3）得到第1～M阶IMF的具体步骤如下：

（a）对Y(t)进行EMD分解得出第1个IMF就终止筛选过程，以减少运算量；

（b）Y(t)每次叠加随机产生的辅助噪声,其均值为零，方差不变；

（c）重复步骤（a）和步骤（b），直到分解得到M个IMF为止。

进一步的，上述步骤（a）中的筛选过程的具体步骤为：

步骤1：初始时定义q=1、Y(t)=e(t)；

步骤2：加入辅助白噪声wn_q(t)，构造h(t)=Y(t)+wn_q(t)；

步骤3：找出序列h(t)的全部极大值与极小值；

步骤4：用三次样条插值构造序列的上下包络时间序列：u(t)和d(t)，计算包络线均值m(t)=(u(t)+d(t))2；

步骤5：求取h(t)=h(t)-m(t)，如h(t)满足以下条件：①任意时刻由局部极大值和局部极小值定义的包络均值为零；②序列过零点的个数与极值点的个数相等或者最多相差一个，则定义h(t)为1个IMF，结束本次筛选，记h(t)为imf₁ ^q(t)，更新q=q+1，返回步骤2直到q=I+1终止全部筛选过程，如h(t)不满足筛选停止条件（①②），返回步骤3继续执行直至满足一次筛选停止条件；

步骤6： ${\mathrm{imf}}_1\left(t\right)=(1/I){\mathrm{\Σ}}_{q=1}^I{\mathrm{imf}}_1^q\left(t\right),$ 其中I为噪声辅助的次数。

其中，加入辅助噪声的次数I理论上是越大越好，但是考虑到实际应用中的实时性要求，应当适当选择噪声加入次数，例如加入噪声的标准差σ=0.01此时可取I=200次会取得较好的效果。依次得到imf₁～imf_M后，再直接对运算的余项做常规EMD得到高阶的IMF分量，这样就得到了双差分测量残差信号完整的IMF分量imf₁～imf_N，N为总的IMF个数。

进一步的，步骤（5）所述的阈值滤波的基础上进行迭代，设共迭代K次，其步骤包括：

A、构造与原信号有相同信噪比的序列 $e^p\left(t\right)=\mathrm{alter}\left({\mathrm{imf}}_1\left(t\right)\right)+{\mathrm{\Σ}}_{i=2}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$ 其中alter(imf₁(t))表示imf₁(t)循环移位或者随机打乱得到的信号序列，p为迭代次数；

B、对ep(t)进行步骤（5）中的阈值滤波得到

${\tilde{e}}^p\left(t\right)={\mathrm{\Σ}}_{i=n}^m{f}_{\mathrm{thr}}\left({\mathrm{imf}}_i\left(t\right)\right)+{\mathrm{\Σ}}_{i=m+1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right);$

C、重复步骤A和步骤B K-1次后，得到最终滤波结果：

$\mathrm{er}\left(t\right)=(1/K){\mathrm{\Σ}}_{P=1}^K{\tilde{e}}^p\left(t\right).$

有益效果：本发明将信号的噪声压缩到低阶IMF中，同时消除噪声位置对EMD分解结果的影响，随机采样第1阶IMF，构造与原信号有相同信噪比的多个序列，对上述序列进行阈值滤波得到多种滤波结果，最后对所有滤波结果取平均得到最终的降噪结果，本发明具有以下的优点：

（1）减小含噪声信号EMD分解中噪声使信号的小幅度波形淹没到噪声IMF中的概率。

（2）降低含噪声信号EMD分解的阶次，由于噪声会改变信号的极值点分布，增加第一阶IMF含噪量，进而提高后续的EMD分解精度，即使信号EMD分解受噪声影响减少。

（3）增加第1阶IMF的含噪量，使得后续的迭代区间阈值滤波中基于第1阶IMF构造与信号具有相同信噪比的新序列更准确。

（4）由于EMD分解结果对随机信号位置具有敏感性，所以多次分解相同信噪比的信号，最后平均处理能减少EMD分解位置敏感性带来的误差。

附图说明

图1为本发明的原理图；

图2为本发明的实施例中的多接收机短基线测量***图；

图3为本发明的实施例中的多径误差提取效果图。

具体实施方式

下面对本发明技术方案结合附图和实施例进行详细说明。

本发明的一种基于EMD迭代阈值滤波的GNSS多径效应抑制方法，包括以下步骤：

（1）对GNSS短基线的双差分残差信号e(t)进行EMD分解，得到一组 ${\mathrm{\Σ}}_{i=1}^N{\mathrm{imf}}_i\left(t\right)$ 和余项r(t)，即 $e\left(t\right)={\mathrm{\Σ}}_{i=1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$

其中，IMF是指本征模态函数，i为模态分量阶数，imf_i(t)为第i阶模态分量序列，N为分解得到的IMF总阶数；

（2）利用奇异谱分析对步骤（1）中的IMF进行噪声估计，确定出噪声占优的前M个IMF分量；

（3）设Y(t)为噪声辅助分析的序列，初始时Y(t)=e(t)，求解imf(t)_i时，将Y(t)逐次加入辅助噪声经EMD分别得到1阶IMF，记结果集为IMFs，设共加入I次噪声，对得到的IMFs取平均，则得到所求IMF分量：

${\mathrm{imf}}_i\left(t\right)=(1/I){\mathrm{\Σ}}_{q=1}^I{\mathrm{imf}}_i^q\left(t\right),$ 且1≤q≤I

Y(t)=Y(t)-imf_i(t)作为分解对象逐次加入噪声继续上述过程得到1阶IMF，依次得到imf₁～imf_M,记R(t)=Y(t)-imf_M为余项；

（4）对R(t)进行EMD分解得到信号的低频IMF分量imf_M+1～imf_N，从而实现含噪声信号的完整EMD分解，N为应用噪声辅助分析后EMD分解得到的IMF总阶数；

（5）基于步骤（4）中的IMF分量进行阈值滤波，滤波后重构的信号表示为 $\mathrm{er}\left(t\right)={\mathrm{\Σ}}_{i=n}^m{f}_{\mathrm{thr}}\left({\mathrm{imf}}_i\left(t\right)\right)+{\mathrm{\Σ}}_{i=m+1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$ 其中f_thr(·)为阈值滤波函数,其中1≤n≤m-1，m、n分别为阈值滤波的上界和下界，m的选择为滤波提供了灵活性；

（6）将er(t)用于GNSS双差分测量结果X(t)的实时补偿数据，X(t)-er(t)即为去掉多径误差后的GNSS测量结果。

进一步的，所述步骤（2）中的噪声的具体估计步骤为：

（I）对任意长度为N的序列x(t)，1≤t≤N，选择长度为m的时延窗口，构造一个m×n的轨迹矩阵X=[X₁,X₂,…,X_n]；

其中X_k=(x(k),x(k+1),…,x(k+m-1))^T，(·)^T表示矩阵的转置，n=N-m+1，2≤m≤N/2，上述过程可以视为在R^N空间中将时间序列x(t)线性映射到向量空间构造出m×n的Hankel矩阵；奇异值分解(SVD)可以将复杂的矩阵分解成几个简单矩阵相乘的形式，具有明显的物理意义，其不要求待分析数据为方阵形式，较特征值分解的形式适用性更广；

（II）设λ₁≥λ₂≥…≥λ_m为XX^T的特征值，其对应的特征向量为u₁,u₂,…u_m，定义 $X=X_1+X_2+\·\·\·+X_m,X_i={\mathrm{\λ}}_i^{\frac{1}{2}}{u}_i{v}_i^T;$

其中 $v_i=X^T{u}_i/{\mathrm{\λ}}_i^{\frac{1}{2}};$ 由于 ${\left|\right|X\left|\right|}^2=\mathrm{trace}\left\{{\mathrm{XX}}^T\right\}={\mathrm{\Σ}}_{i=1}^m{\mathrm{\λ}}_i,$ trace{·}表示矩阵的迹，Xi²=λ_i(1≤i≤m)，

（III）定义为分量X_i对X能量的贡献，Ψ越小说明贡献越小，以含噪声信号的IMF分量为分析对象即可绘制多个分量序列的奇异谱。由于白噪声奇异值相差不大，所以其奇异向量对信号能量变化的贡献是稳定的，即白噪声的奇异谱是平坦的，而传感器输出的非平稳信号在延迟窗内其奇异值变化较大，能量变化多集中到前几个奇异值，其奇异谱有陡峭的下降趋势，基于上述现象可以对IMF分量进行噪声水平的定性评判。

进一步的，所述步骤（3）中加入的辅助噪声为白噪声，其满足（0，σ²）分布，σ的取值范围为0.01～0.5。

进一步的，所述步骤（3）得到第1～M阶IMF的具体步骤如下：

（a）对Y(t)进行EMD分解得出第1个IMF就终止筛选过程，以减少运算量；

（b）Y(t)每次叠加随机产生的辅助噪声，其均值为零，方差不变；

（c）重复步骤（a）和步骤（b），直到分解得到M个IMF为止。

进一步的，上述步骤（a）中的筛选过程的具体步骤为：

步骤1：初始时定义j=1、Y(t)=e(t)：加入辅助白噪声wn_j(t)，构造h(t)=Y(t)+wn_j(t)；

步骤2：找出序列h(t)的全部极大值与极小值；

步骤3：用三次样条插值构造序列的上下包络：u(t)和d(t)，计算包络均值m(t)=(u(t)+d(t))2；

步骤4：求取h(t)=h(t)-m(t)，如h(t)满足以下条件：①任意时刻由局部极大值和局部极小值定义的包络均值为零；②序列过零点的个数与极值点的个数相等或者最多相差一个，则定义h(t)为1个IMF，结束本次筛选，返回步骤1，记h(t)为imf₁ ^j(t)，更新j=j+1，重复上述步骤直到j=I终止筛选，如h(t)不满足筛选停止条件（①②），返回步骤2继续执行直至满足筛选停止条件；

步骤5： ${\mathrm{imf}}_1\left(t\right)=(1/I){\mathrm{\Σ}}_{j=1}^I{\mathrm{imf}}_1^j\left(t\right),$ 其中I为噪声辅助的次数。

其中，加入辅助噪声的次数I理论上是越大越好，但是考虑到实际应用中的实时性要求，应当适当选择噪声加入次数，例如加入噪声的标准差σ=0.01此时可取I=200次会取得较好的效果。依次得到imf₁～imf_M后，再直接对运算的余项做常规EMD得到高阶的IMF分量，这样就得到了双差分测量残差信号完整的IMF分量imf₁～imf_N，N为总的IMF个数。

进一步的，步骤（5）所述的阈值滤波的基础上进行迭代，设共迭代K次，其步骤包括：

A、构造与原信号有相同信噪比的序列，由于实际应用中只含有噪声的IMF分量是很难存在的，第1阶IMF中很可能混有小幅度的波动信号，本发明将噪声进一步压缩至低阶IMF中，故第1阶IMF为噪声的可信度进一步提高，随机采样imf1(t)的数据可以选择的方法有随机打乱和循环移位两种，alter(·)表示随机采样过程，构造的新序列如下所示：

$e^p\left(t\right)\mathrm{alter}\left({\mathrm{imf}}_1\left(t\right)\right)+{\mathrm{\Σ}}_{i=2}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$ 其中alter(imf₁(t))表示imf₁(t)循环移位或者随机打乱得到的信号序列，p为迭代次数；

B、基于双差分残差信号只包含热噪声和多径误差的特征，多径误差主要存在于残差信号的高阶IMF中，考虑到动态接收机应用中的多径效应可能存在高频成分，所以对低阶IMF进行阈值滤波，对高阶IMF和余项直接重构，f_thr(·)为阈值处理函数，为每次对imf₁(t)随机取样构造新序列的除噪结果，其形式如下对e^p(t)进行步骤（5）中的阈值滤波得到：

${\tilde{e}}^p\left(t\right)={\mathrm{\Σ}}_{i=n}^m{f}_{\mathrm{thr}}\left({\mathrm{imf}}_i\left(t\right)\right)+{\mathrm{\Σ}}_{i=m+1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right);$

C、由于随机信号的位置对EMD分解的结果有很大的影响，这种位置敏感性特征会导致分解误差，并最终影响除噪结果，构造与原含噪声信号具有相同信噪比的序列，多次滤波后对多个结果求平均能降低这种误差，即：重复步骤A和步骤B K-1次后得到最终滤波结果：

$\mathrm{er}\left(t\right)=(1/K){\mathrm{\Σ}}_{P=1}^K{\tilde{e}}^p\left(t\right).$

得到er(t)后，即可用于补偿GNSS观测信号，接收机运行于静态模式时直接补偿同一时刻同一地点的GNSS测量信号，运行于动态模式下的接收机当多径误差变化不剧烈时亦可以采用此种方法。

实施例：

本实施例主要包括以下步骤：

（1）建立GNSS短基线测量模型，接收机的分布如图（2）所示，包含1个主接收机、3个从接收机的差分测量***，各接收机间的距离均已知，记主接收机与从接收机1的距离为d1=21.140m，利用卫星定位的测量结果，获取到另一组主接收机与从接收机1的距离数据d2，记e(t)=d1-d2为包含多径误差的短基线测量误差；

（2）对e(t)进行常规EMD分解，确定需要噪声辅助分析的IMF个数，本例取M=2，选择辅助噪声为高斯白噪声wn(0,σ²)，其中σ=0.01,迭代次数为I=200次，采用本发明提出的EMD改进分解方法得到imf₁～imf₇；

（3）选择n=2,m=6（基于多径误差多为低频并集中在最后1阶IMF和余项r(t)中的经验知识），进行阈值滤波处理，得到滤波后的imf₂～imf₆，一次阈值滤波迭代的结果为： ${\tilde{e}}^p\left(t\right)={\mathrm{\Σ}}_{i=n}^m{f}_{\mathrm{thr}}\left({\mathrm{imf}}_i\left(t\right)\right)+\mathrm{\Σ}{\mathrm{imf}}_7\left(t\right)+r\left(t\right);$ alter(·)选择随机打乱方法，取迭代次数为10次，最终的滤波结果为 $\mathrm{er}\left(t\right)=(1/10){\mathrm{\Σ}}_{P=1}^{10}{\tilde{e}}^p\left(t\right).$ 其滤波效果如图（3）所示，图左侧从上到下依次为包含多径误差和接收机热噪声的基线误差信号，右侧分别为其对应的多径误差；

（4）将er(t)作为图（2）所示的多接收机短基线测量***的补偿数据，得到处理前后的基线测量误差如下所示：

测量数据	基线误差1	基线误差2	基线误差3
				滤波前RMS（m）	0.1384	0.1946	0.1863
滤波后RMS（m）	0.0223	0.0236	0.0232
				改善度(%)	83.87	87.85	87.54

表中RMS为均方误差，使用本发明方法提取出多径误差后，发现三种基线测量误差均显著降低。

Claims (4)

1.一种基于EMD迭代阈值滤波的GNSS多径效应抑制方法，其特征在于包括以下步骤：

（1）对GNSS短基线的双差分残差信号e(t)进行EMD分解，得到一组IMF和余项，分别为 ${\mathrm{\Σ}}_{I=1}^N{\mathrm{imf}}_i\left(t\right)$ 和r（t），即 $e\left(t\right)={\mathrm{\Σ}}_{I=1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$

其中，IMF是指本征模态函数，i为模态分量阶数，imf_i(t)为第i阶模态分量序列，N为分解得到的IMF总阶数，t为残差信号的时间变量；

（2）利用奇异谱分析对步骤（1）中的IMF进行噪声估计，确定出噪声占优的前M个IMF分量；

（3）设Y(t)为噪声辅助分析的序列，初始时设Y(t)=e(t)，求解imf(t)_i时，将Y(t)逐次加入辅助噪声经EMD分别得到1阶IMF，记结果集为IMFs，设共加入I次噪声，对得到的IMFs取平均，则得到所求IMF分量：

且1≤q≤I

将Y(t)=Y(t)-imf_i(t)作为分解对象逐次加入噪声继续上述过程得到1阶IMF，依次为imf₁～imf_M，记R(t)=Y(t)-imf_M为余项；其中，q是指第q次加入辅助噪声得到的IMF编号；

（4）对R(t)进行EMD分解得到信号的低频IMF分量，依次为imf_M+1～imf_N，进而完成含噪声信号的完整EMD分解，N为应用噪声辅助分析后EMD分解得到的IMF总阶数；

（5）基于步骤（4）中的IMF分量进行阈值滤波，滤波后重构的信号表示为 $\mathrm{er}\left(t\right)={\mathrm{\Σ}}_{i=n}^m{f}_{\mathrm{thr}}\left({\mathrm{imf}}_i\left(t\right)\right)+{\mathrm{\Σ}}_{i=m+1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$ 其中f_thr(imf_i(t))为阈值滤波函数，其中1≤n≤m-1，m、n分别为阈值滤波的上界和下界；

（6）将er(t)用于GNSS双差分测量结果X(t)的实时补偿数据，X(t)-er(t)即为去掉多径误差后的GNSS测量结果。

2.根据权利要求1所述的基于EMD迭代阈值滤波的GNSS多径效应抑制方法，其特征在于：所述步骤（3）中加入的辅助噪声为白噪声，该白噪声满足（0，σ²）分布，σ的取值范围为0.01～0.5。

3.根据权利要求1所述的基于EMD迭代阈值滤波的GNSS多径效应抑制方法，其特征在于：所述步骤（3）得到第1～M阶IMF的具体步骤如下：

（a）对Y(t)进行EMD分解得出第1个IMF就终止筛选过程；

（b）Y(t)每次叠加随机产生的辅助噪声，其均值为零，方差不变；

（c）重复步骤（a）和步骤（b），直到分解得到M个IMF为止。

4.根据权利要求1所述的基于EMD迭代阈值滤波的GNSS多径效应抑制方法，其特征在于：将所述步骤（5）的阈值滤波进行迭代，设共迭代K次，其步骤包括：

A、构造与原信号有相同信噪比的序列 $e^p\left(t\right)=\mathrm{alter}\left({\mathrm{imf}}_1\left(t\right)\right)+{\mathrm{\Σ}}_{i=2}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right),$ 其中，alter(imf₁(t))表示imf₁(t)循环移位或者随机打乱得到的信号序列，p为迭代次数，N为IMF总阶数；

B、对e^p(t)进行步骤（5）中的阈值滤波得到：

${\tilde{e}}^p\left(t\right)={\mathrm{\Σ}}_{i=n}^m{f}_{\mathrm{thr}}\left({\mathrm{imf}}_i\left(t\right)\right)+{\mathrm{\Σ}}_{i=m+1}^N{\mathrm{imf}}_i\left(t\right)+r\left(t\right);$ m、n分别为阈值滤波的上界和下界；

C、重复步骤A和步骤BK-1次后，得到最终滤波结果：

$\mathrm{er}\left(t\right)=(1/K){\mathrm{\Σ}}_{P=1}^K{\tilde{e}}^p\left(t\right).$