IGS基准站水平方向的运动以板块运动为主，而高程方向变化比较复杂，影响因素比较多^[1]。通过对连续参考站高程时间序列的研究分析，可以提取周期项和趋势项，获取地球振荡周期变化信息^[2]，进而研究造成测站非线性运动的物理机制，精化各种误差改正模型，获得测站准确的位置和线性速度，为建立和维持动态地球参考框架提供基础数据^[3]，为地球动力学方面的应用提供参考。

对GPS站坐标时间序列进行分析时，通常假定跟踪站存在周年和半周年等周期性运动，利用谐波函数进行拟合，求解各周期项的振幅和相位^[4]。但实际上，不同的GPS站高程时间序列的周期性变化非常复杂，周期可能并不严格等于周年和半周年，振幅也随时间发生变化。且已有研究表明，GPS站高程时间序列除存在半周年和周年运动，还存在其他周期信号^[5-6]。因此，常用的谐波拟合法并不能很好地适用于非平稳、非线性的GPS信号，而自适应性信号处理方法——经验模态分解(empirical mode decomposition，EMD)可以解决这些问题。EMD分解得到的IMF(intrinsic mode function)分量具有清晰的物理意义，揭示了信号里蕴含的多尺度振荡特征；对IMF分量进行Hilbert变换(HHT)，得到的频谱图可以反映信号的局部特征，以及能量随时间的分布。但EMD方法容易产生模式混叠现象，整体经验模态分解(ensemble empirical mode decomposition，EEMD)借助噪声可以减弱这种现象。为了研究GPS高程时间序列的周期振荡特性，本文利用HHT-EEMD方法获取IGS站的周期运动特征，并用小波变换对其进行比较和评价，同时结合环境负载数据和GRACE负载等初步探讨引起时间序列周期振荡的物理机制，以消除非构造运动的影响并提高GPS观测成果的应用价值。

1 HHT-EEMD算法的基本理论

Huang等^[7]于1998年提出一种分析非线性非平稳时间序列的新方法——经验模态分解(EMD)和希尔伯特(Hilbert)谱分析，统称为经验模态分析。经验模态分析的基本思路是：一个时间序列由多个时间尺度的振荡波构成，设法从经验资料中把这些固有的、内在的本征模态函数(IMF)分量逐级分离出来，并求其Hilbert谱，通过分析IMF分量及其相应的谱，即时-频分析，得到原序列的多尺度振荡特征^[8]。

1.1 经验模态分析(HHT-EMD)

经验模态分析方法实现的步骤如下：

1) 把物理系统的实测序列分解为频率逐渐降低的本征模态分量(IMF)和一个趋势项：首先逐级分离出IMF分量，分解完成后剩余的即为趋势项，这一过程称为EMD。

2) 对分解出的各个IMF分量作Hilbert变换，得到各自的瞬时振幅和瞬时频率，然后把瞬时振幅表示在时频表面上，就得到了Hilbert谱，这一过程称为Hilbert-Huang变换，简称HHT。

1.2 整体经验模态分解(EEMD)

信号出现中断时，EMD方法分解的IMF会产生模式混叠现象^{[7, 9]}。模式混叠现象指一个IMF分量包含着尺度互异的信号，或者尺度相似的信号存在于不同的IMF分量中，主要是由于信号中断引起的，不仅混淆了信号的时频分布特性，更重要的是使信号IMF分量的物理意义变得不清晰^[5]。模式混叠现象导致了EMD结果的不稳定性和不唯一性。为了消除模式混叠，使用一种借助噪声的资料分析方法，即整体经验模态分解EEMD。EEMD的想法是在进行每一次EMD分解之前，在原始序列信号x(t)中加入一定量的高斯白噪声w_k(t)，第k次待分解的信号x_k(t)可表示为:

$ {x_k}\left( t \right) = x\left( t \right) + {w_k}\left( t \right) $

(1)

经过EMD分解会生成一系列由高频至低频的IMF分量。每一次由系统随机生成的固定方差的白噪声序列都是不同的，通过重复分解过程会得到信号多次EMD分解的IMF分量，将相同尺度的IMF取平均值，即可得到EEMD最终分解的IMF分量。由于高斯白噪声是零均值的正态随机序列，当重复分解次数足够大时，加入白噪声的影响将被消除，从而分解得到的IMF全部来自于原始信号本身。

2 利用HHT-EEMD方法分析时间序列

采用IGS分析中心之一SOPAC(ftp://garner.ucsd.edu/pub/timeseries/)免费提供的全球坐标时间序列作为数据源。本文选用中国大陆的BJFS站，时间跨度为2000-01-01~2015-01-01，参考框架为ITRF2008。对BJFS站高程时间序列分别采用EMD和EEMD两种分解方法获取IMF分量。采用Wu等^[10]的建议，整体平均的次数选用100。为了研究不同比例白噪声对IMF分量的影响，加入的白噪声与原始信号的标准差之比分别取0.1、0.2、0.4，分别以EEMD01、EEMD02、EEMD04来表示。

2.1 EMD/EEMD比较

图 1显示了BJFS站2000~2015年高程方向的时间序列。从图 1可见，2010年以后数据缺失较为频繁(见粉色框)，典型数据异常如图中红色椭圆所示。

BJFS站时间序列经过EMD分解为12个频率逐渐减小的本征模态分量(IMF)和一个残差项。IMF分量如图 2所示，由图可见，每个分量均随时间呈周期性变化，从上到下周期依次增大，代表不同尺度的信号。因此EMD提供了一个很好的手段与方法来获得IGS站高程时间序列中的各种周期项。由于周期信号的成分比较复杂，某些周期信号并不表现为规则的正弦曲线形式，EMD分解的IMF分量振幅是随时间变化的，可以很好地刻画出信号的变化状态，因此EMD具有很强的从原始时间序列识别非线性信号成分的能力^[6]。但是图中部分IMF的周期不明，例如C9~C11分量均表现出模式混叠现象。为了克服模式混叠，利用EEMD方法对BJFS站高程时间序列进行分析，采用了3种方案(EEMD01，EEMD02，EEMD04)，与EMD的比对见图 3。

从图 3可见，3种EEMD分解结果比较一致，而EMD分解结果与EEMD有较大偏差。显然，C12应为序列中的长周期项，EEMD表明其周期大约为3 a，而EMD分解结果在2000~2009年区间上周期表现为1.5 a。C11分量EMD分解结果看不出明显的周期项，多种周期信号混叠在一起，从EEMD结果可以得出其应该为2 a项。可以看出，EMD分解时C12和C11分量均混入了周年信号，这是由于数据缺失和数据异常引起的模式混叠现象。同时，EEMD结果表明，C10与C9分量均为周年信号，而EMD仍然出现不同程度的混叠现象(图中紫色椭圆)。可以看出，EMD分解的C10分量混入了2 a信号，C9分量混入了0.5 a信号，说明EEMD方案可有效减弱模式混叠对分解后IMF分量的影响^[5]。然而，图中红色椭圆处表明，EEMD分解得到的2 a信号在2005~2009年区间上出现明显的周年信号，EEMD分解得到的周年信号也局部表现为0.5 a信号，这说明EEMD分解并没有完全消除模式混叠现象，只是在一定程度上予以减弱。EEMD结果表明，C8分量近似为0.5 a项，局部表现出模式混叠。图中只列出了C7~C12分量，另外C1~C4为高频分量，同时由表 1可得，C7为季节性信号，C6周期为2个月，C5周期为1个月。

经验模态分解只能大致看出IMF的周期项，并不能细致分析每个IMF的周期信息。将上面4种方法分解得到的各IMF作傅里叶频谱分析，图 4为EMD分解后IMF分量的Fourier频谱图，不同的颜色代表不同的IMF分量Fourier变换后的频谱。傅立叶变换是将图像从空间域转换到频率域，幅值最大的尖峰对应的频率就是IMF分量的周期。由图可见，不同颜色的频谱对称排列，这是由于FFT结果的对称性，通常只使用前半部分的结果，即小于采样频率一半的结果。部分IMF表现出的周期特性统计见表 1。

由表 1可得，EEMD01、EEMD02与EEMD04分解结果比较一致，只在频率较大的几个IMF分量上有微小差异，说明白噪声的比例对结果几乎没有影响。EMD分解结果则与EEMD结果存在差异，主要表现在周期较大的C11与C12分量上，这与上文分析得到的结论一致。总体来说，BJFS站存在近似2 a、1 a、0.5 a、0.25 a、2个月、1个月以及长周期变化。其中2 a项由于其信号强度很弱，很难通过传统的方法识别出来，说明EEMD方法具有从复杂序列中提取微弱信号的能力^[6]。由于几种EEMD方案的差异很小，下文提到的EEMD均指EEMD02分解方法。

图 5为BJFS站EMD/EEMD分解后的残差。结果表明，3种EEMD残差变化趋势比较一致，BJFS站在2000~2015年高程方向有逐渐升高的趋势，且变化速率逐渐减小。EEMD年线性变化速率为1.06 mm/a，小于EMD的趋势项变化速率为1.58 mm/a。EMD残差的趋势和EEMD比较一致，但残差曲线并不光滑，还存在着微小的波动，可能受到了序列中异常值的干扰。而EEMD结果由于随机噪声中和了异常值，残差曲线相对光滑。

综上所述，当GPS高程时间序列信号存在中断情形时，EMD方法产生明显的模式混叠现象，模式混叠导致了EMD的不稳定性和不唯一性。EEMD分解相比EMD可以明显减弱模式混叠对本征模态分量的影响，但并没有完全消除模式混叠现象。EEMD结果中还发现了潜在的近似2 a信号。此外，对比发现，加入白噪声的比例对EEMD的结果几乎没有影响，说明EEMD方法是稳定可靠的。总之，EEMD与EMD方法相比有较大的优势，但还需进一步改进。具体要从加入白噪声的振幅以及参与集合的数目方面深入研究。

2.2 HHT变换分析时间序列

进行EMD的主要目的是为了对每一个IMF进行Hilbert变换(HHT)并得到Hilbert谱，该谱能够精确地反映信号的能量在时间和频率上的分布规律。为了得到BJFS站高程运动的主要贡献项，本文对BJFS站EMD及EEMD分解的IMF分量作Hilbert-Huang变换，结果见图 6~7。

图 6显示了BJFS站EMD分解后12个IMF分量作Hilbert-Huang变换(HHT)以后的频谱图。它为时间-频率-振幅的谱图，颜色的变化表示能量的相对大小，而能量的大小则反映为IMF分量的振幅变化。可以看出，EMD分解得到的IMF分量混叠在一起，没有很好地分离，看不出每个分量的具体周期，反映出上文所述的模式混叠现象。

下面具体分析EEMD分解的IMF分量HHT后的结果。由图 7可见，随着频率的减小，4条频谱曲线依次对应C9~C12分量。其中C9与C10能量最大，为图中的红黄色曲线，频率在0.003 cpd附近波动，大致对应的是周年信号。其次为C11对应的曲线，为图中的蓝绿色曲线，频率在0.001 5 cpd附近波动，对应2 a信号。C12则为长周期信号，其余IMF分量的能量较小。这一结果表明，BJFS站周期运动主要以周年和2 a运动为主。

2.3 基于小波变换的比较和验证

为了验证HHT-EEMD方法在时间序列分析上的有效性，接下来对高程时间序列应用小波变换进行处理以对比分析。通过离散小波变换，一个信号可以分解为它的近似部分和细节部分。近似部分代表了信号的低频成分，对应较大的尺度；而细节部分代表了信号的高频成分，对应较小的尺度^[10]。以db4为基底对时间序列作小波分解的结果如下。

图 8为BJFS站时间序列经小波分解后得到的8层图像，分别对应8个高频分量。为了更详细地分析其周期特性，对部分高频分量作了傅里叶频谱分析，周期统计结果见表 2。小波分解结果表明，BJFS站的高频分量包括1 a、0.5 a、0.25 a、2个月以及1个月信号。此结果与EEMD的分解结果基本一致，验证了上述EEMD的分解结果。然而受小波基的限制，小波分解结果并没有得到2 a信号，EEMD分解得到的2 a信号是否真实存在还需要进一步证明。

小波变换不仅可以对序列进行分解，还可以对分解后的分量进行信号重构。经小波重构最终得到降噪数据。同样，EEMD也可以对信号进行降噪处理。将EEMD分解得到的C8~C12分量进行低通滤波:

$ {x_r}\left( t \right) = \sum\limits_{i = 8}^{12} {{C_i}\left( t \right)} + {r_n}\left( t \right) $

(2)

同时，选取db4和sym6两个小波基进行小波滤波，滤波结果见图 9。

由图 9可见，EEMD的重构信号与原始序列更加吻合，可以更好地捕获原始序列的波峰和波谷，尤其在数据缺失的地方。db4、sym6两种小波滤波曲线在时间序列异常值的位置均出现毛刺。相比小波滤波结果，EEMD滤波曲线较为光滑，没有毛刺，说明EEMD滤波结果受原始信号的影响较小，可以剔除瞬时强噪声^[11]。

用均方根(RMS)来定量比较3种滤波结果的优劣：

$ {\rm{RM}}{{\rm{S}}_u} = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{\left( {{x_i} - {l_i}} \right)}^2}} }}{n}} $

(3)

式中，n表示样本总数，x_i表示原始序列，l_i表示各种滤波序列。

表 3为3种滤波结果的RMS统计。由表 3可得，EEMD低通滤波的RMS小于db4和sym6两种小波滤波的RMS，表明由EEMD方法得到的重构信号更接近原始序列，验证了上文的结论。因此在滤波方面EEMD低通滤波效果更好。

通过对比可以看出，HHT-EEMD的优势在于：1)EEMD根据其自适应特性进行分解，不依赖于事先给定的周期拟合模型，即EEMD不用考虑基底函数的选择问题；2)Hilbert频谱图中的频率是瞬时频率，为时间的函数，可以表征信号的局部特征；3)HHT不是对原序列作变换，而是对分解以后的分量作变换，能更精确地反映出每个分量代表的具体物理意义，而且EEMD低通滤波也有更好的去噪重构效果。因此，HHT-EEMD在分析非平稳、非线性时间序列时具有一定的优越性。

3 GPS站周期性信号的来源分析

GPS位置时间序列中的周期信号主要由非构造形变引起^{[3, 12-13]}。引起非构造形变的地球物理因素主要包括两大类：一是潮汐形变；二是大气、非潮汐海洋、积雪和土壤水等质量负荷，也称为环境负载^[12]。为了验证HHT-EEMD方法得到的周期性信号确实存在于IGS站高程时间序列中，本文通过分析环境负载和GRACE负载来得到高程时间序列中存在的周期性振荡信号，并解释其形成原因。

3.1 环境负载的影响

本文采用EOST loading service提供的由环境负载引起的位移。主要考虑大气压负载、水文负载、非潮汐海洋负载，参考以往对环境负载的研究^{[3, 6]}，分别采用ATMIB、GLDAS、ECCO 3种全球负载模型，得到的测站位移属于CF框架。为了分析高程运动中各种周期项与负载的关系，将BJFS站3种环境负载引起的测站位移时间序列进行功率谱分析。功率谱密度(PSD)曲线见图 10。

由图可见，3种环境负载的功率谱密度均遵循逆幂律，是典型的随机过程^[6]。图中众多尖峰反映出了季节性，最强的周期分量出现在尖峰对应的频率上。3种负载引起的测站位移中均存在周年、半周年与季节性信号。能量最大的尖峰出现在周年频率附近，其中大气负载造成的测站位移最大，水文负载对测站位移的影响次之，非潮汐海洋负载造成的位移最小。由此可得，大气负载、水文负载、非潮汐海洋负载均会造成高程时间序列中的周年、半周年及季节性运动，并且对周年项的贡献最大。此外，图中的2 a信号并不明显，可能是被其他信号所掩盖。

3.2 GRACE负载的影响

为了证明2 a信号的存在，同样采用EOST loading service提供的基于GRACE系数的负载形变，时间跨度为2003~2016年。对GRACE负载形变进行功率谱分析，结果见图 11。

由图可知，GRACE负载形变中包含了周年、半周年与季节性信号，仍然是周年信号的能量最大。除此之外，GRACE负载中还发现了近似2 a信号以及长周期信号，说明BJFS站周期运动中确实存在2 a信号以及周期约为3 a的长周期信号，验证了EEMD的结果。因此EEMD方法在探测时间序列中的潜在信号方面具有重要作用。在计算测站垂直方向速度时2 a信号是应该被改正掉的，至于2 a变化具体由哪些因素引起，还需要进一步分析。

4 结语

1) 高程时间序列的周期信号。利用经验模态分解和整体经验模态分解方法，得到BJFS站高程时间序列不仅存在1 a、0.5 a、0.25 a、2个月、1个月以及长周期等周期性振荡，还存在以往方法很难探测出来的2 a信号。当数据存在异常时，经EMD分解后得到的某些IMF分量出现了模式混叠现象，EEMD可以明显减弱这种现象，但无法完全消除。Hilbert-Huang变换(HHT)得到的三维频谱图表明，周年和2 a振荡是BJFS站高程运动的主要贡献项。

2) 经验模态分析(HHT-EEMD)在数据资料分析中的有效性。对高程时间序列应用小波变换进行处理的结果与EEMD分解结果基本一致，但EEMD根据其自适应特性进行分解，不用考虑基底函数的选择问题，而且EEMD低通滤波也有更好的去噪重构效果。HHT的结果能够表征信号的局部特征，物理意义明确。因此，HHT-EEMD在高程时间序列分析中具有一定的优越性。

3) 非构造形变对时间序列的影响。通过环境负载和GRACE负载引起测站位移的功率谱分析得到，环境负载确实会造成IGS基准站的周年、半周年以及季节性运动，其中周年运动振幅最大。GRACE负载的功率谱密度也得到了同样的结论，而且还发现了2 a信号以及长周期信号的存在。至于2 a信号的形成原因，还需要进一步探索。

本文虽然只分析了BJFS站，但HHT-EEMD方法同样也适用于其他IGS站的高程时间序列分析。