大多数基于EMD降噪的研究采用相关系数准则确定分界IMF函数^[1-6]，取互相关系数达到第1个局部极小值时的IMF作为分界IMF，而不能直接给出分界IMF函数的K值，且这些研究都是将低频有用分量进行累加重构实现降噪目的，未考虑高频噪声与低频信号IMF分量的个数差异对计算量的影响。本文基于EMD方法，提出一种降噪的新思路，采用平均周期与能量密度乘积指标的方法来自动地确定分界IMF的K值，并与相关系数准则所确定的分界IMF进行对比，验证该方法的准确性。针对IMF低频信号分量个数多于高频噪声分量个数的情况，将噪声分量进行重构，然后从原始信号中扣除重构的噪声，从而达到降噪的目的。

1 EMD降噪理论的改进

根据EMD降噪的基本原理^[7]可知，将低频IMF分量进行重构，能够达到降噪的目的。然而，当分界IMF之下的低频模态分量个数多于高频IMF噪声分量的个数时，直接重构会增加计算量，甚至会带入少量的系统误差。针对该问题，本文提出一种新的思路来获得降噪后的“干净”信号，可以减少计算量。原始数据经EMD分解后，从原始数据序列中减去重构后的高频噪声分量，得到“干净”信号。降噪信号的表达式为：

$ \hat{x}(t)=x(t)-\sum\limits_{k=1}^{K} \operatorname{IMF}_{k} $

(1)

式中，x(t)为原始序列，$\hat{x}(t)$为降噪后信号，K为噪声与真实信号的分界IMF所对应顺序数。

此外，考虑到利用相关系数准则确定分界IMF的复杂性，以及不能直接给出分界IMF函数的K值，本文在张恒璟等^[8]的基础上，利用平均周期(T_k)与能量密度(E_k)的乘积(ET_k)作为指标，给出一种自动确定分界IMF的算法，能够直接完成分界IMF函数K值的确定^[8]。

信号与噪声分界点确定的阈值表达式为：

$ R_{k-1}=\left|\frac{E T_{k}}{\left(\frac{1}{k-1} \sum\limits_{i=1}^{k-1} E T_{i}\right)}\right| $

(2)

式中，k≥2，当R_k-1≥C(本文中C取2)时，则判定k为分界点，文中N为各个模态的数据长度。将前k-1个噪声主导的模态分量进行重构，用原始数据序列减去重构的噪声，得到降噪后的信号。

利用相关系数(ρ)、均方根误差(RMSE)、信噪比(R_sn)等传统评价指标来衡量去噪的效果^[9]。相关系数越接近1，降噪信号与原始数据序列的相似度越高，拟合效果越好，去噪效果也就越好；均方根误差体现了降噪信号与原始信号之间的偏差程度，值越小，去噪效果越好；信噪比主要是体现噪声信号在整体信号中的比重，其值越大，去噪效果越好。

2 模拟数据降噪分析

由于实测数据含有误差，首先采用真值已知的模拟信号对本文的方法进行验证。模拟信号的采样频率为1 Hz，采样点数为1 024个，noise是信噪比为4 dB的高斯白噪声，其构成分量信号波形图如图 1所示，加噪后所得的模拟原始信号及其频谱图如图 2所示。模拟数据可用式(3)表示：

$ \left\{\begin{array}{l} y_{1}=4 \sin \left(\frac{2 \pi t}{800}\right) \sin \left(\frac{2 \pi t}{250}\right) \\ y_{2}=2 \sin \left(\frac{2 \pi t}{600}\right) \\ y_{3}=\sin \left(\frac{2 \pi t}{50}\right) \\ \varepsilon=\text { noise } \\ x=y_{1}+y_{2}+y_{3}+\varepsilon \end{array}\right. $

(3)

利用相关系数准则确定模拟信号的噪声与信号的分界IMF，计算的各IMF分量与模拟加噪数据序列之间的相关系数如图 3所示。由图 3可知，第4个IMF分量第1次取得局部极小值，所以将前4个分量视为高频噪声，从原始加噪数据序列中剔除高频噪声，获得降噪后的纯净数据序列。

采用平均周期与能量密度的乘积指标(ET_k)确定分界IMF，所得到的平均周期与能量密度以及ET_k如表 1所示。由表 1可知，IMF₁~IMF₃ 3个分量平均周期与能量密度的乘积指标的均值为5.817，而IMF₄的指标为19.159，明显大于均值的2倍，将前3个分量视为高频噪声分量。本文方法确定分界IMF的K值时，所计算的R值分别为0.600，0.728，3.293，此时的K值为4，即前3个IMF分量为高频噪声，剩余的分量为低频真实信号。本文给出部分IMF分量的波形图及其对应的频谱，如图 4所示。图 5(a)为利用本文方法所得的重构去噪信号与原始不加噪的模拟信号对比图，图 5(b)为重构噪声与真实噪声的对比图。从图 5(a)中可看出，本文方法去噪后信号波形图与真实未加噪信号波形图拟合效果较好，图 5(b)中重构噪声与真实模拟噪声图基本吻合，验证了本文方法用于降噪的可靠性。

由于本文方法确定分界IMF函数的K值与相关系数准则所得结果不同，分别计算2种方法所得“干净”信号与真实未加噪信号的相关系数、均方根误差及信噪比，所得结果如表 2所示。表 2中，RPN₁表示真实信号与原始信号的各评价指标，DPN₁表示降噪信号与原始信号的各评价指标，DPR₁表示本文方法所得降噪信号与真实信号的各评价指标，DPR₂表示采用相关系数准则所得降噪信号与真实信号的各评价指标。由表 2可知，相较于DPN的各评价指标值，DPN₁中ρ值与R_sn值更大，而RMSE值更小，说明本文方法降噪效果较优；与DPR₂中的各指标值相比，DPR₁中ρ值与R_sn值更大，而RMSE值更小，表明本文方法确定的分界IMF函数降噪效果优于相关系数准则。综合表 2中的各项指标值，验证了本文方法降噪效果的可靠性。

根据本文方法所确定的分界IMF函数的K值可知，本文所模拟的加噪原始数据序列的低频真实信号分量的个数明显大于高频噪声分量的个数，采用传统方法直接利用低频信号进行去噪信号重构需要累加6N₁(N₁表示模拟数据的长度)次；而本文方法只需将前3个高频噪声分量进行累加，需2N₁次运算，加上N₁次减法运算，总共只需3N₁次运算，明显比直接重构得到“干净”信号的计算量要少。综上所述，在模拟数据降噪分析中，本文方法能够自动确定分界IMF函数的K值，较之相关系数准则法降噪效果更精确，且能够减少计算量。

3 GPS实测高程时间序列降噪分析

为了进一步验证本文方法的可靠性，使用北京房山(BJFS)站的实测高程(U)时间序列来进行降噪研究与分析。数据来自中国地震局GNSS数据产品服务平台(http://www.cgps.ac.cn/)，时间跨度为2002~2018年，采样间隔为1/365.25年，最大采样频率为365.25 Hz。由于实测数据中含有一些粗差，进行分析之前利用Hector软件对序列中较大的粗差进行剔除，同时对剔除粗差之后的序列进行趋势项估计，将剔除粗差之后的坐标高程时间序列作为本文降噪分析的原始数据序列，所得结果如图 6所示。从图中可以看出，剔除粗差后的时间序列图的波动明显小于实测信号图，所估计的趋势项与U方向时间序列拟合效果也较好，且高程时间序列呈现出明显的周期性，说明确实达到了较好的剔除粗差的效果。

用EMD方法对BJFS站原始高程时间序列进行处理，得到11个本征模态分量(IMF)，计算各IMF分量与剔除粗差之后的高程时间序列的相关系数，所得结果如图 7所示。从图中可以明显看出，第5个IMF与原数据序列的相关系数第一次取得局部极小值，即前5个高频IMF为噪声，剩余的为真实信号。

利用各IMF分量的平均周期与能量密度的乘积指标确定分界IMF，所得结果如表 3所示。由表 3可知，前3个IMF分量的平均周期与能量密度的乘积指标(ET_k)的平均值为29.842，IMF₄的ET_k指标为69.604，明显大于前3个IMF分量均值的2倍，即第4个IMF分量为分界IMF，前3个为高频噪声分量，剩余的IMF分量为真实信号。根据本文自动确定分界IMF分量K值的办法，在满足R阈值的情况下，R的取值依次分别为1.303 5，1.648 6，2.332 4，当R为2.332 4时停止筛选，此时自适应K值为4，即前3个IMF分量为噪声，后8个IMF分量为真实信号，此结果与利用相关系数所求得的结果不一致。为了比较本文方法与平均周期与能量密度指标方法确定分界IMF的准确性，并且对最后的降噪效果进行评价，分别利用2种不同结果重构后信号与原始信号的相关系数、均方根误差以及信噪比进行对比分析，结果如表 4所示。由表 4可知，本文算法在实测数据中能够直接给出分界IMF的K值，且相比于相关系数准则所得降噪信号，本文方法所得降噪信号与原始信号的ρ值与R_sn值更大，RMSE值更小，表明本文方法降噪效果比相关系数准则更佳。

图 8给出了部分IMF分量的波形，前3个IMF分量为高频噪声，第4个IMF分量为分界IMF函数。根据本文提出的降噪思路，将前3个高频噪声分量进行重构，然后用原始高程时间序列减去重构后的高频噪声，得到降噪的高程时间序列，如图 9所示。由图 9可知，降噪后“干净”信号较之原始高程时间序列明显更清晰，验证了本文方法应用于GPS高程时间序列降噪的有效性。本文方法得到“干净”信号只需3N₂(N₂表示实测数据的长度)次运算，传统的EMD重构方法需要进行7N₂次累加运算，减少了计算量。

4 结语

本文在EMD降噪分解过程中，结合平均周期与能量密度乘积的指标，给出一种自动确定分界IMF分量K值的算法，并且与相关系数准则确定分界IMF分量的方法进行对比。模拟数据结果表明，本文算法所得的降噪信号较之相关系数准则所得降噪信号，相关系数提高0.017，RMSE减小0.156，信噪比提高5.142 dB。BJFS站的实测GPS高程时间序列结果表明，本文所提出的改进EMD方法较之传统EMD方法，相关系数提高0.04，均方根误差减小0.51 mm，信噪比提高3.14 dB。模拟数据和实测GPS高程时间序列数据的结果都表明，改进的EMD方法用于降噪分析精度优于传统EMD方法，是一种更为精确的降噪方法。

同时，本文针对EMD分解过程中高频噪声IMF分量的个数少于低频真实信号IMF分量的个数的情况，提出一种新的EMD降噪思路。模拟数据结果表明，本文的降噪方法相对于传统的直接重构算法，减少了4N₁次累加运算。GPS实测高程时间序列结果表明，本文提出的新方法较之直接重构方法降噪，减少了4N₂次累加运算。综上所述，本文所提出的基于EMD的GPS高程时间序列降噪新方法，能够自动确定分界IMF函数的K值，减少计算量。但是，当高频噪声IMF分量个数多于低频噪声分量个数时，该方法并不适用，且信号与噪声分界的阈值为经验值，需要对该阈值的选取展开进一步的研究。