EMD信号分析方法是Huang等^[1]于1998年提出，主要用于非线性、非平稳时间序列的处理，利用EMD方法对GNSS坐标时间序列进行降噪处理是当前研究的热点。张双成等^[2-3]利用EMD方法将线性趋势项、周期项、噪声从原始序列中有效分离出来，并准确地从GNSS时间序列中提取出周期项。罗飞雪等^[4]将交叉证认的方法应用于EMD滤波去噪，自适应地选择低频IMF分量的个数，可实现噪声与信号的分离。贾瑞生等^[5]将EMD方法与独立成分分析方法相结合，提出一种新的微震信号降噪方法。王文波等^[6]提出一种基于主成分分析(PCA)的EMD降噪方法，可进一步提取分界IMF分量中的信号。孙伟峰等^[7]根据连续均方根误差(CMSE)准则来判定信号与噪声的分界点，提出一种改进的EMD降噪方法。陈凤林^[8]将每个IMF分量与原始信号p值第1次取得极小值时作为信号与噪声的分界点，基于此提出一种可靠的指标来确定模态分界点。张恒璟等^[9-10]利用EMD方法从GPS高程时间序列中提取振荡周期，采用平均周期与能量密度乘积指标来改进CMSE准则，识别信号与噪声的分界IMF分量。

然而，大多数研究采用相关系数准则确定分界IMF分量，取互相关系数达到第1个局部极小值时作为分界IMF分量，但有时所确定的模态分界点并不准确，且无法直接给出分界IMF分量的K值。CMSE准则的本质为各IMF分量的能量密度或均方根误差，在实际过程中，若连续的IMF分量的均方根误差比较接近，则很难判定信号与噪声的分界点。本文基于EMD方法提出一种降噪的新思路，综合考虑r值及RMSE值，采用复合评价指标来自动确定分界IMF分量的K值，并与相关系数准则所确定的分界IMF分量进行对比，验证本文方法的准确性。

1 EMD降噪理论及其改进 1.1 EMD降噪基本原理

EMD方法将信号分解为若干频率由高到低的IMF分量和1个残余项，将残余项与部分低频IMF分量重构可以实现降噪。每个IMF分量必须满足2个条件^[11]：1)在整个序列中，极值点的数量与过零点的数量最多相差1个；2)由局部极大值所构成的上包络线及局部极小值所构成的下包络线的平均值为0。然而，在实际信号分解过程中，IMF分量很难严格满足第2个条件，因此给出每个IMF分量停止筛选的阈值表达式：

${\rm{SD}} = \mathop \sum \limits_{t = 0}^{N - 1} \left[ {\frac{{{{\left( {{c_k}\left( t \right) - {c_{k - 1}}\left( t \right)} \right)}^2}}}{{c_k^2\left( t \right)}}} \right]$

(1)

式中，c_k(t)、c_k-1(t)为每个IMF分量筛选过程中2个相邻的数据序列；SD为每个IMF分量停止筛选的阈值，通常情况下取0.2~0.3。

信号x(t)进行EMD降噪的过程为：1)分别找出原始数据序列x(t)中所有极大值点和极小值点，计算上、下包络线的均值m₁，将原始数据序列减去该均值，从而得到新的数据序列c₁(t)；2)重复步骤1)，直到满足IMF阈值条件，则认为得到IMF分量；3)将原始数据序列x(t)减去第1个IMF分量，形成新的数据序列x₂(t)，然后重复步骤1)与步骤2)，得到m个IMF分量和1个残余项w(t)，当残余项满足单调条件时才符合要求。因此，原始数据序列可表示为：

$x\left( t \right) = \mathop \sum \limits_{k = 1}^m {\rm{IM}}{{\rm{F}}_k} + w\left( t \right)$

(2)

原始数据序列经EMD后，需要确定噪声与真实信号的分界IMF分量。通过相关系数准则选取分界IMF分量，相关系数第1次取得极小值时所对应的IMF分量为分界IMF分量。EMD方法将分界IMF分量并入噪声部分，分界之后对低频IMF分量及残差项进行重构，得到降噪信号。降噪后的信号可表示为：

$\hat x\left( t \right) = \mathop \sum \limits_{k = K + 1}^m {\rm{IM}}{{\rm{F}}_k} + w\left( t \right)$

(3)

式中，$\hat x\left( t \right) $为降噪后的信号；IMF_k为分界IMF分量，即k之前的IMF分量均为纯噪声IMF分量。

1.2 EMD降噪方法改进

考虑到利用相关系数准则确定分界IMF分量的复杂性及利用单一指标确定分界IMF分量的不准确性，本文采用复合评价指标T^[12]，同时考虑信号的细节信息和逼近信息，计算曲线的RMSE值与r值，给出一种自动确定分界IMF分量的算法，直接完成分界IMF分量K值的确定。为方便公式表达，将残余项视为最后1个IMF分量。

均方根误差计算公式为：

${\rm{RMSE}} = \sqrt {\frac{1}{N}\mathop \sum \limits_{t = 0}^{N - 1} {{\left( {{\rm{IM}}{{\rm{F}}_k}\left( t \right) - x\left( t \right)} \right)}^2}} {\rm{}}$

(4)

平滑度计算公式为：

$r = \frac{{{{\mathop {\sum\limits_{t = 0}^{n - 2} {{{\left( {{\rm{IM}}{{\rm{F}}_k}\left( {t + 1} \right) - IM{F_k}\left( t \right)} \right)}^2}} }\nolimits_{}^{} }^{}}}}{{\mathop {\sum\limits_{t = 0}^{n - 2} {{{\left( {x\left( {t + 1} \right) - x\left( t \right)} \right)}^2}} }\nolimits_{}^{} }}$

(5)

式中，IMF_k(t)为第k个IMF分量，k取值为1，2，…，m+1；x(t)为带噪的原始数据序列。

将RMSE与r两个指标进行归一化处理：

${P_{{\rm{RMSE}}}} = \frac{{{\rm{RMSE}}}}{{\mathop {\sum {\left( {{\rm{RMSE}}} \right)} }\nolimits^{} }}$

(6)

${P_{\rm{r}}} = \frac{r}{{\mathop {\sum {\left( r \right)} }\nolimits^{} }}$

(7)

采用变异系数定权法对归一化后的2个指标进行赋权处理，定权过程为：

${\rm{C}}{{\rm{V}}_{{P_{{\rm{RMSE}}}}}} = \frac{{{\sigma _{{P_{{\rm{RMSE}}}}}}}}{{{\mu _{{P_{{\rm{RMSE}}}}}}}}{\rm{}}$

(8)

${\rm{C}}{{\rm{V}}_{{P_r}}} = \frac{{{\sigma _{{P_r}}}}}{{{\mu _{{P_r}}}}}{\rm{}}$

(9)

${W_{{P_{{\rm{RMSE}}}}}} = \frac{{{\rm{C}}{{\rm{V}}_{{P_{{\rm{RMSE}}}}}}}}{{{\rm{C}}{{\rm{V}}_{{P_{{\rm{RMSE}}}}}} + {\rm{C}}{{\rm{V}}_{{P_r}}}}}{\rm{}}$

(10)

${W_{_{{P_r}}}} = \frac{{{\rm{C}}{{\rm{V}}_{{P_r}}}}}{{{\rm{C}}{{\rm{V}}_{{P_{{\rm{RMSE}}}}}} + {\rm{C}}{{\rm{V}}_{{P_r}}}}}$

(11)

式中，σ、μ分别为标准差及均值，CV为变异系数，W为基于变异系数定权的权值。

复合评价指标T的表达式为：

$ T = {W_{_{{P_{{\rm{RMSE}}}}}}} \times {P_{{\rm{RMSE}}}} + {W_{{P_r}}} \times {P_r} $

(12)

信号与噪声分界阈值表达式为：

$ {R_{k - 1}} = \left| {\frac{{{T_{k - 1}}}}{{{T_k}}}} \right| $

(13)

式中，k≥2，当阈值首次取1≤R_k-1≤3时，则判定k-1为分界点，即分界IMF分量的K值为k-1。将前k个IMF分量视为噪声，k之后的IMF分量视为信号，将信号主导的IMF分量进行重构，得到降噪后的数据序列。

利用p(相关系数)、RMSE及SNR等传统评价指标分析去噪效果。p值越接近于1，降噪信号与原始数据序列的相似度越高，拟合效果越好，去噪效果也越好；RMSE可体现降噪信号与原始信号之间的偏差程度，数值越小表明去噪效果越好；SNR主要体现噪声信号在整体信号中的比重，其值越大，去噪效果越好。

2 模拟数据降噪分析

GPS坐标时间序列一般由季节项、趋势项及噪声3个部分组成。先剔除模拟数据Ⅰ、Ⅱ、Ⅲ中站点位置、趋势项及阶跃式偏移，考虑到3个恒定振幅的周期项和噪声，设置采样频率为1 Hz，采样点数为1 024个，表达式为：

$ \left\{ \begin{array}{l} {y_1} = 5{\rm{sin}}\left( {\frac{{2{\rm{\pi }}t}}{{600}}} \right){\rm{sin}}\left( {\frac{{2{\rm{\pi }}t}}{{350}}} \right)\\ {y_2} = 7\sin \left( {\frac{{2{\rm{\pi }}t}}{{500}}} \right)\\ {y_3} = 2{\rm{sin}}\left( {\frac{{2{\rm{\pi }}t}}{{50}}} \right)\\ \varepsilon = {\rm{noise}}\\ x = {y_1} + {y_2} + {y_3} + \varepsilon \end{array} \right. $

(14)

模拟数据Ⅰ、Ⅱ中分别加入信噪比为4 dB、6 dB的高斯白噪声；模拟数据Ⅲ中加入白噪声与幂律噪声的组合噪声，白噪声振幅为5 mm，幂律噪声振幅为0.02 mm，幂律噪声的谱指数为-1.2。

由于连续GPS坐标序列中含有振幅时变季节性信号，其他6个模拟数据在剔除站点位置、趋势项及阶跃式偏移的前提下，加入振幅变化因子^[13]。通过下式分别模拟10 a的坐标序列时变季节性信号数据：

$ \begin{array}{l} {\rm{y}}\left( {{t_i}} \right) = {\rm{asin}}\left( {2\pi {t_i}} \right) + bcos\left( {2\pi {t_i}} \right) + c\left( {{t_i}} \right)sin\left( {2\pi {t_i}} \right) + \\ \;\;\;\;\;c\left( {{t_i}} \right)cos\left( {2\pi {t_i}} \right) + dsin\left( {4\pi {t_i}} \right) + ecos\left( {4\pi {t_i}} \right) + \\ \;\;\;\;\;\;\;\;c\left( {{t_i}} \right)sin\left( {4\pi {t_i}} \right) + c\left( {{t_i}} \right)cos\left( {4\pi {t_i}} \right) + \varepsilon \left( {{t_i}} \right) \end{array} $

(15)

式中，y(t_i)为模拟的时变季节性信号数据；a、b、d、e为常数；t_i为GPS年积日，c(t_i)=2e^0.3sin(t_i)为振幅变化因子，ε(t_i)为噪声。

在模拟数据Ⅳ、Ⅴ、Ⅵ中，a、b、d、e取值分别为2 mm、3 mm、4 mm、5 mm。模拟数据Ⅳ、Ⅴ中加入信噪比为4 dB、6 dB的高斯白噪声；模拟数据Ⅵ中加入白噪声与幂律噪声的组合噪声，白噪声振幅为5 mm，幂律噪声振幅为0.02 mm，幂律噪声的谱指数为-1.2。

为验证本文方法与所取振幅的大小无关，在模拟数据Ⅶ、Ⅷ、Ⅸ中，a、b、d、e取值分别为6 mm、7 mm、8 mm、9 mm。模拟数据Ⅶ、Ⅷ中分别加入信噪比为4 dB、6 dB的高斯白噪声；模拟数据Ⅸ中加入白噪声与幂律噪声的组合噪声，白噪声振幅为5 mm，幂律噪声振幅为0.02 mm，幂律噪声的谱指数为-1.2。

表 1、2为9个模拟数据计算得到的各IMF分量与原始信号之间的相关系数值(p)及利用本文方法计算得到的复合评价指标值(T)。由表 1、2可知，在9个模拟数据中，除模拟数据Ⅲ、Ⅴ外，其他模拟数据计算的相关系数值与复合评价指标值所确定的分界IMF分量相同，即所确定的K值相同，信号与噪声的分界点一致。在模拟数据Ⅲ、Ⅴ中，相关系数准则确定的IMF分量K值为4，复合评价指标确定的IMF分量K值为3。

图 1为模拟数据Ⅲ、Ⅴ的加噪信号波形，图 2为基于模拟数据Ⅲ采用相关系数与复合评价指标所确定的重构信号与真实信号波形。从图 2可以明显看出，复合评价指标确定的重构信号与真实信号的波形非常接近，而相关系数准则确定的重构信号与真实信号的波形存在一定差异。图 3为基于模拟数据Ⅴ采用相关系数与复合评价指标确定的重构信号与真实信号波形图。从图 3可以看出，在历元端点处，相较于相关系数准则确定的重构信号，复合评价指标所确定的重构信号与真实信号波形更为接近。在其他历元处，复合评价指标与相关系数准则确定的重构信号与真实信号的波形无明显差异。上述分析表明，复合评价指标准则比相关系数准则所确定的K值更加准确，改进的EMD方法比传统EMD方法提取真实信号的效果更佳，去噪效果更优。

为定量说明2种准则的去噪精度，表 3为模拟数据Ⅲ、Ⅴ利用这2种准则得到的重构信号与真实信号的3种降噪效果评价指标值。由表 3可知，在模拟数据Ⅲ中，与互相关系数准则得到的各项评价指标值相比，复合评价指标得到的重构序列与真实信号的p值提高0.048 4，SNR值提高11.554 8 dB，RMSE值降低1.091 1；在模拟数据Ⅴ中，与互相关系数准则得到的各项评价指标值相比，复合评价指标得到的重构序列与真实信号的p值提高0.024 5，SNR值提高8.719 dB，RMSE值降低1.127 4。综上所述，在9个模拟数据中，模拟数据Ⅲ、Ⅴ采用本文复合评价指标确定的K值相较于相关系数准则确定的K值更准确；在其余7个模拟数据中，2种指标准则确定的K值相同，从而验证本文方法确定分界IMF分量的可靠性。

3 实测GPS高程时间序列

本文利用乌鲁木齐(URUM)站与武汉(WUHN)站的实测高程(U)时间序列进行降噪研究与分析。数据来源于中国地震局GNSS数据产品服务平台(http://www.cgps.ac.cn/)，URUM站所取数据为2009~2016年跨度为7 a的高程时间序列数据，WUHN站所取数据为2006~2016年共10 a的高程时间序列数据，2个站点的采样间隔均为1/365.25 a，最大采样频率为365.25 Hz。在进行EMD方法去噪之前，利用Hector软件剔除序列中较大的粗差，将去除粗差后的时间序列作为本文降噪分析的原始数据序列，所得结果见图 4、5。从图 4、5可以看出，在某些时刻，原始数据序列的波动明显小于实测信号，表明Hector软件能很好地剔除实测数据中的粗差。

表 4为利用URUM站与WUHN站实测高程数据计算得到的各IMF分量与原始信号之间的p值及利用本文方法计算得到的T值。分析表 4可知，基于URUM站实测高程时间序列数据，相关系数第1次取得极小值为0.197 1，确定的分界IMF分量K值为3；R值第1次取得符合条件值为2.97，此时确定的IMF分量K值也为3。基于WUHN站实测高程时间序列数据，相关系数第1次取得极小值为0.188 5，R值第1次取得符合条件值为1.56，2个准则确定的IMF分量K值均为4。综上分析可知，本文给出的复合评价指标准则能够准确地确定分界IMF分量K值，即可以准确找出信号与噪声模态的分界点，表明本文降噪改进方法具有可靠性。

4 结语

本文针对EMD降噪过程中采用相关系数单一指标确定分界IMF分量K值存在不准确性的问题，综合考虑信号r值及RMSE值，提出一种新的复合评价指标准则来确定分界IMF分量K值，从而改进传统的EMD降噪方法。该复合评价指标准则给出信号与噪声分界的阈值表达式，当阈值R_k-1首次处于^{[1, 3]}时，判定k-1为分界IMF分量K值。

本文通过对9个模拟数据及2个实测数据进行实验，证实该方法能够准确确定真实信号与噪声的分界点。模拟数据Ⅰ、Ⅱ、Ⅳ、Ⅵ、Ⅶ、Ⅷ、Ⅸ与URUM站及WUHN站的实测GPS高程时间序列数据结果表明，相关系数准则与本文复合评价指标准则所确定的分界IMF分量K值相等，从而验证本文方法确定分界IMF分量的可行性。模拟数据Ⅲ、Ⅴ的分析结果表明，本文给出的复合评价指标准则确定的分界IMF分量K值比相关系数准则确定的分界IMF分量K值更准确，表明本文方法的可靠性。由于本文方法确定的信号与噪声分界点更准确，因此基于新准则改进的EMD方法去噪精度优于基于相关系数准则的EMD方法。但本文仅利用2个GPS实测站高程时间序列进行验证，存在一定局限性，下一步需要对更多站点的时间序列数据进行研究与分析。另外，考虑到实测数据真值未知的情况，后期研究需确定一种去噪评价指标来说明本文改进的EMD方法降噪效果的优劣。