EMD在GNSS时间序列周期项处理中的应用 Application of EMD to GNSS Time Series Periodic Term Processing 1 2 first-author 刘希康，博士生，主要从事地表形变与地下结构方面的研究工作。xikangliu@126.com 刘希康，博士生，主要从事地表形变与地下结构方面的研究工作。xikangliu@126.com xikangliu@126.com 1 corresp 丁志峰，博士，研究员。dingzf@cea-igp.ac.cn 丁志峰，博士，研究员。dingzf@cea-igp.ac.cn dingzf@cea-igp.ac.cn 2 2 中国地震局地球物理研究所，北京，100081 中国地震局地球物理研究所，北京，100081 Institute of Geophysics, China Earthquake Administration, Beijing 100081, China Institute of Geophysics, China Earthquake Administration, Beijing 100081, China 中国地震局第一监测中心，天津，300180 中国地震局第一监测中心，天津，300180 The First Crust Monitoring and Application Center, China Earthquake Administration, Tianjin 300180, China The First Crust Monitoring and Application Center, China Earthquake Administration, Tianjin 300180, China 全球导航卫星系统（global navigation satellite system，GNSS）位置时间序列蕴含有丰富的构造和非构造变形信息，具有成分复杂、建模困难、非构造信息难以有效分离等特点，利用自适应的经验模态分解（empirical mode decomposition，EMD）方法对川滇地区24个GNSS连续站时间序列作周期项修正。结果表明，周期项的修正十分必要，EMD方法能够根据每个台站信号的自身特性，自适应地提取不同频率、振幅的周期成分，这也更符合实际情况；相较于谐波模型，EMD方法对原始时间序列在N、E、U方向的改正均更加精确有效。使用修正后的连续站时间序列模拟流动观测，发现经过5~6 a/期观测即可得到相对可靠的台站运动速度，并通过距离较近的实际连续观测站对流动观测站周期项改正，验证了EMD方法的稳定性和可靠性，这也为流动GNSS观测实施、周期修正和资料使用提供了参考意见和理论依据。 Objectives 2 The position time series of global navigation satellite system (GNSS) contain rich tectonic and non-tectonic deformation information, and the complex components are difficult to accurately model and effectively separate non-tectonic information. It is very important to remove the non-tectonic deformation information for the accurate and effective application of the observation data. Methods 2 Empirical mode decomposition（EMD）is an adaptive time-frequency processing method to correct the period term of time series at 24 GNSS continuous stations and 2 mobile stations in Sichuan Province and Yunnan Province. Results 2 The results show that the correction of the period term is necessary, and the EMD method can extract the different frequency and amplitude period components adaptively according to its own characteristics of signal at each station. Compared with the original time series, the average root mean square error after correction is reduced significantly in N, E and U directions. And EMD method is more accurate and effective than harmonic model modification. We use the modified continuous station time series to simulate the mobile observations, and find that relatively reliable motion velocities can be obtained after an observation period of 5⁃6 years. The stability and reliability of EMD method are verified by correcting the period term of the mobile station with the actual continuous observation station at a close distance. Conclusions 2 This paper provides a reference and theoretical basis for the implementation of mobile GNSS observation and the correction of observation data. 全球导航卫星系统 周期运动 流动观测 谐波模型 经验模态分解 global navigation satellite system（GNSS） periodic motion mobile observation harmonic models empirical mode decomposition 全球导航卫星系统（global navigation satellite system，GNSS）技术能够提供高精度、大尺度、多维的地壳形变监测数据，在地壳形变监测领域应用已经十分广泛，也发挥着越来越重要的作用。1998年以来，中国地壳运动观测网络（简称网络工程）和中国大陆构造环境监测网络（简称陆态网络）建立了260多个连续观测台站和2 000多个流动观测站，其中90%以上流动观测站有5年及以上观测资料，部分超过20期。GNSS时间序列蕴含着丰富的构造和非构造变形信息，如海潮、大气、积雪、土壤水、海洋非潮汐等，如何去除其中的非构造变形信息对观测资料的准确、有效应用起着至关重要的作用。 诸多研究成果表明，中国的GNSS时间序列在水平方向以线性运动为主 ［ 1 1 - 2 2 1-2 ］ ，并有弱周期性变化，而高程方向则表现为复杂的周期振荡与非线性运动，且不同测站的周期振荡特性并不相同。在传统模型中，周期性信号一般通过周年或半周年正余弦谐波表示，可用最小二乘拟合方法分离，得到的振幅和相位是常数值。而实际上，不同年份间的GNSS测站环境变化及其响应具有不规则性 ［ 3 3 - 7 7 3-7 ］ ，用传统的谐波模型求解非构造周期性运动参数可能偏离甚至违背原时间序列的运动特性，进而导致用于噪声分析的残差不可避免地包含某种误差或错误信息 ［ 8 8 ］ 。 近年来，许多学者对振幅变化的周期项处理进行了研究，文献［ 9 9 ］提出一种基于卡尔曼滤波的方法来统计GNSS时间序列的周期性现象，但该模型仅适用于随机信号，而很多研究证实GNSS时间序列中含有闪烁噪声；文献［ 10 10 ］发现奇异谱分析（singular spectrum analysis，SSA）方法不需先验信息，比最小二乘拟合方法和卡尔曼滤波方法对变振幅周期项的提取效果更优；文献［ 11 11 ］利用多通道奇异谱分析（multichannel singular spectrum analysis，MSSA）方法提取了区域站点的共模误差；文献［ 12 12 ］利用MSSA估计了影响GNSS观测的几种常见环境效应因素。然而SSA方法和其改进的MSSA方法对滞后窗口的取值比较敏感，选取不合适的滞后窗口甚至会导致错误的结果 ［ 5 5 ］ 。文献［ 13 13 ］基于主成分分析（principal component analysis，PCA）方法发现全球大部分GNSS台站都存在非线性周期规律，周年项和半周年项占据主导地位；文献［ 7 7 ］对比分析了区域叠加法和PCA方法对高程时间序列共模分量的提取效果，并给出两种方法的使用建议。以上研究均验证了GNSS时间序列中变振幅周期信号的存在，并根据不同方法的性能对信号进行了分离，但如何准确地从包含了各种噪声、误差及周期性非构造运动信息的实际观测资料中有效提取出构造运动信息，还需不断探索不同的思路和方法。 经验模态分解（empirical mode decomposition，EMD）是近年发展起来的新型自适应信号时频处理方法，并成功应用在机械故障诊断、生物医学信号特征提取、地震数据分析等领域 ［ 14 14 ］ ，它不像MSSA或PCA方法需要多个信号一起分析提取共模分量，而仅依靠信号自身的时间尺度特性将信号分解为一定数量的本征模态函数（intrinsic mode function，IMF），IMF由观测信号唯一确定，各IMF分量包含了信号不同时间尺度的局部特征，具有清晰的物理意义，通过IMF分量可以提取信号不同周期的频率成分信息 ［ 15 15 ］ 。与传统谐波分析方法相比，EMD是直接的、后验的，不需要提前知晓信号的周期性成分及其对信号的影响程度。 目前，基于EMD发展的几种方法在GNSS数据处理方面已有许多探索研究和初步应用，主要包括震前异常特征信息提取 ［ 16 16 - 19 19 16-19 ］ 、时间序列噪声识别与信号提取 ［ 8 8 ， 20 20 ］ 和多路径模型改正 ［ 21 21 - 22 22 21-22 ］ 3个方面。本文将利用EMD方法对GNSS连续站时间序列进行周期项改正，并进一步利用改正后的时间序列模拟流动观测，查看积累多少期观测资料的情况下才能获得相对准确的运动速度，并用于实际数据周期项改正处理。 1　EMD方法原理 1 EMD方法就是基于信号本身的局部特征时间尺度，对原始信号进行平稳化处理，将复杂的信号分解成 m 个具有不同特征尺度的数据序列，计算式为： <math display="block" id="M1"><mi>x</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><mi>i</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>m</mi></mrow></munderover></mstyle><mrow><mi mathvariant="normal">I</mi><mi mathvariant="normal">M</mi><msub><mrow><mi mathvariant="normal">F</mi></mrow><mrow><mi>k</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo></mrow><mi mathvariant="normal">R</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">s</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117546&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117545&type= 44.36533356 8.55133343 式中， <math id="M2"><mi>x</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117548&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117547&type= 7.11199999 3.55599999 为原始信号； <math id="M3"><mi>t</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117550&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117549&type= 1.26999998 2.79399991 为原始信号时间； <math id="M4"><mi mathvariant="normal">I</mi><mi mathvariant="normal">M</mi><msub><mrow><mi mathvariant="normal">F</mi></mrow><mrow><mi>k</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117552&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117551&type= 13.20800018 4.48733330 为第 k 个IMF分量；Res <math id="M5"><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117554&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117553&type= 4.40266657 3.55599999 为趋势项。 对信号分解后，利用相关系数准则 ［ 8 8 ］ 确定 m 个IMF信号与噪声的分解本征模态函数 <math id="M6"><mi mathvariant="normal">I</mi><mi mathvariant="normal">M</mi><msub><mrow><mi mathvariant="normal">F</mi></mrow><mrow><mi>k</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117531&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117530&type= 8.21266651 3.80999994 。通过重构被分解序列获得原信号中的周期项、趋势项及噪声序列。 根据GNSS解算后时间序列可知，其采样间隔为1 d，即1/365 a，它的倒数即为最大采样频率，为365 Hz。由奈奎斯特采样定理可知，时频分析能识别出最大采样率的一半，即182.5 Hz，即可得到频谱的全部信息，换算得到的最小周期为1/182.5 a，也就是理论上能识别出周期最小为2 d的信号。 2　模拟信号分析 1 GNSS坐标时间序列一般由周期项、趋势项和噪声3个部分组成。由于实际连续GNSS坐标序列中周期性信号的振幅可能随时间变化，在剔除站点位置和阶跃式偏移前提下，模拟信号加入更接近实际的振幅变化因子 ［ 23 23 ］ ，计算式为： <math display="block" id="M7"><mtable><mtr><mtd><mi>x</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>=</mo><mi>a</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mi>b</mi><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mo stretchy="false">(</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mi>c</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mi>c</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mo stretchy="false">(</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mi>d</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mn mathvariant="normal">4</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>+</mo></mtd></mtr><mtr><mtd><mi>e</mi><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mo stretchy="false">(</mo><mn mathvariant="normal">4</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mi>c</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mn mathvariant="normal">4</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>+</mo></mtd></mtr><mtr><mtd><mi>c</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mo stretchy="false">(</mo><mn mathvariant="normal">4</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117533&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117532&type= 127.00000000 15.15533447 <math display="block" id="M8"><mo>+</mo><mi>ε</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117535&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117534&type= 10.32933331 4.48733330 式中， <math id="M9"><mi>x</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117537&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117536&type= 7.70466709 4.48733330 为添加时变周期性的模拟信号； <math id="M10"><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117569&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117568&type= 1.86266661 3.72533321 为第 i 年的时间； a 、 b 、 d 、 e 为常数，分别取值2、4、7、9； <math id="M11"><mi>ε</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117571&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117570&type= 6.85799980 4.48733330 为噪声； <math id="M12"><mi>c</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117555&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117572&type= 6.94266701 4.48733330 为振幅变化因子， <math id="M13"><mi>c</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>=</mo><mn mathvariant="normal">6</mn><msup><mrow><mi mathvariant="normal">e</mi></mrow><mrow><mn mathvariant="normal">0.3</mn><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117556&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117573&type= 22.77533340 4.74133301 。当 <math id="M14"><mi>c</mi><mo stretchy="false">(</mo><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117574&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=45117557&type= 6.94266701 4.48733330 为常数时， 式（2） 式（2） 则转化为传统的谐波周期模型。 图1 图1 为按 式（2） 式（2） 设计的时间长度为10 a的模拟信号， 图1 图1 （a）、 图1 图1 （b）分别为固定振幅和时变振幅的年周期、半年周期信号， 图1 图1 （c）为趋势项、周期信号+趋势项以及添加10 dB白噪声的合成信号。 图2 图2 为模拟信号的EMD分解。根据相关系数准则，将前5个分量（IMF1~IMF5）视为高频噪声，第6、第7个分量（IMF6、IMF7）为周期信号，第8个分量（Res）为趋势项。 图3 图3 为利用EMD方法和传统谐波模型提取的周期信号与模拟周期信号的对比图，其相关系数分别为0.991 4和0.969 1。由 图3 图3 可以看出，EMD方法提取的振幅随时间变化，且与模拟周期信号高度一致，其分解出的模态分量个数只与信号自身特性相关，而谐波模型需提前给出年周期和半年周期先验信息用于拟合，且振幅为固定值。此外，由 图2 图2 的IMF6和 图3 图3 可以看出，在两端位置由于端点效应出现幅值波动现象，在实际数据处理时可以通过数据延拓或舍弃两端部分数据来消除此问题的影响 ［ 24 24 ］ ，本文在实际数据处理时，两端各舍弃约30天数据。 3　GNSS时间序列周期项提取与应用 1 3.1　GNSS数据处理 2 本文选取陆态网络和网络工程的GNSS连续和流动观测数据，利用GAMIT/GLOBK ［ 25 25 ］ 软件解算。首先，利用GAMIT获得国际GNSS服务（international GNSS service，IGS）站与区域GNSS测站的单日松弛解，天线相位中心改正采用绝对天线相位中心模型，大气折射模型采用VMF1映射函数，全球固体潮、海潮模型改正遵循IERS2003协议，但未进行大气荷载模型改正。然后，将IGS站点作为固定框架站点，利用GLOBK进行平差，估算出区域站在ITRF2014参考框架下的点位坐标时间序列及误差。最后，修复时间序列因天线更换、强震同震位移等影响造成的阶跃，并剔除异常值。为了便于周期信号的提取，采用最小二乘法去除时间序列中的线性趋势项 ［ 12 12 ］ 。 3.2　周期项提取分析 2 为了检验谐波模型和EMD方法对周期项的改正效能，本文选取川滇地区数据记录完整率在95%以上的24个连续站和2个流动站（见 图4 图4 ），对其2011-01—2020-12期间的时间序列进行处理。因篇幅限制，以SCXD站N方向时间序列为例，对实际数据的EMD周期项提取进行说明，如 图5 图5 所示。 由 图5 图5 （a）可知，SCXD站N方向时间序列EMD分解的10个分量中，利用相关系数法确定IMF1~IMF5为高频噪声，IMF6~IMF9为周期项重构信号，Res为趋势项。其中，IMF6主要为0.5 a信号，此外还有弱0.4 a信号；IMF7主要为1.0 a信号和弱0.76 a信号；IMF8主要为2.67 a信号。由于信号总长度为10 a，由奈奎斯特采样定理可知，其最大只能表征周期为5 a的周期信号，故IMF9可能为10 a或更长趋势周期信号，快速 傅里叶变换（fast Fourier transform，FFT）给出其周期为10.67 a。 图5 图5 （b）为谐波模型周期项修正后的时间序列EMD分解后的IMF分量，由IMF6~IMF9分量可以发现，虽然经过谐波模型改正，但残余信号依然存在比较明显的不同成分周期信号。 图6 图6 为SCXD站N方向原始时间序列经EMD方法和谐波模型改正后的时间序列及其功率谱密度。 图6 图6 （b）为EMD改正后的残余时间序列利用FFT求取的功率谱密度，其不再含有明显的周期信号，主要为高频噪声。 图6 图6 （c）为谐波模型改正后残余时间序列的功率谱密度，结合 图5 图5 （b）可以发现，虽然幅值有所降低，但依然残留有显著的0.5 a、1.0 a和2.67 a周期信号，说明谐波模型改正并不能很好消除周期信号的影响。整体来看，原始时间序列经EMD改正后比谐波模型改正后更加平稳，经谐波模型和EMD改正后，均方根误差（root mean square error，RMSE）由1.90 mm分别降至1.69 mm和1.46 mm，表明EMD方法对周期信号提取更具优越性。 在实际资料利用EMD方法分解时，IMF分量一般小于13个，其中有效周期分量一般不多于7个，因此本文利用FFT探测其周期成分时主要选取前7种。利用EMD方法提取川滇地区24个台站N、E、U三方向数据的周期项，经FFT分别取前7种主要频率成分，统计结果如 图7 图7 所示。结果显示，N、E方向周期成分较为分散、复杂，主要周期为0.5 a和1.0 a，此外还有明显的1.1 a周 期项，但并不是每一个台站都具有年周期和半年周期成分，特别是E方向有6个台站不具有年周期；U方向优势周期较多且集中，24个台站均具有0.5 a和1.0 a周期项，此外还有显著的0.9 a、 1.1 a和1.2 a成分。值得注意的是，小于0.5 a周期并不是其主要的7个成分之一，虽然在单个台站资料处理时确实发现有小于半年周期项的存在，但在周期成分中占比较低。此外，部分台站三方向还具有2.0 a、2.5 a或3.5 a长周期成分，导致此现象的原因尚不明确，需要进一步研究。 表1 表1 为24个台站采用谐波模型和EMD方法提取周期项的相关系数，其中N、E、U方向平均相关系数分别为0.65、0.49和0.86，说明两种方法提取到的周期项具有较好一致性。综合 图7 图7 分析，整体上N方向振幅相对E方向大，周期信号明显，集中在0.5 a和1.0 a附近，相关系数较高。反之，E方向除优势周期外，其他信号成分分布较散，相关系数低，不同台站周期项形态差异也较大，其原因可能是水平方向的时间序列除了受区域整体年周期、半年周期因素影响外，还叠加了局部因素影响。U方向信号具有明显的周期性，且相关系数较高，说明两种方法在U方向对周期信号的提取均有较好的结果，但从形态和功率谱密度来看，EMD方法提取的周期信息更丰富、更全面。 2 站名 3 相关系数 2 站名 3 相关系数 N E U N E U SCMB 0.69 0.43 0.83 YNLJ 0.76 0.61 0.91 SCML 0.60 0.65 0.89 YNMJ 0.52 0.37 0.83 SCMN 0.67 0.36 0.91 YNMZ 0.36 0.60 0.82 SCPZ 0.70 0.54 0.90 YNRL 0.65 0.42 0.89 SCSM 0.74 0.53 0.91 YNSD 0.57 0.48 0.93 SCXC 0.61 0.79 0.91 YNTC 0.67 0.27 0.92 SCXD 0.73 0.42 0.90 YNTH 0.75 0.46 0.83 SCYY 0.63 0.43 0.90 YNYA 0.73 0.52 0.84 XIAG 0.52 0.56 0.80 YNYM 0.56 0.14 0.85 YNCX 0.65 0.58 0.88 YNYL 0.79 0.40 0.72 YNDC 0.62 0.55 0.79 YNYS 0.83 0.39 0.92 YNLC 0.53 0.51 0.85 YNZD 0.72 0.67 0.79 表2 表2 为24个台站时间序列的RMSE，结果表明，谐波模型和EMD方法对原始序列在N、E、U三方向均有明显的改进 ［ 26 26 ］ ，但二者改善程度不同，EMD方法比谐波模型对非构造周期信息的提取更充分，对周期项的修正效果更好。 2 站名 3 原始序列 3 谐波改正 3 EMD改正 N E U N E U N E U SCMB 2.26 2.58 6.89 2.10 2.53 5.40 1.92 2.28 4.90 SCML 2.17 2.93 9.71 2.08 2.77 6.28 1.82 2.61 5.29 SCMN 1.78 1.67 6.60 1.58 1.63 4.46 1.37 1.49 4.21 SCPZ 1.61 1.90 6.98 1.45 1.86 4.72 1.29 1.67 4.30 SCSM 1.95 1.93 7.82 1.78 1.85 5.23 1.65 1.70 4.72 SCXC 1.70 1.77 6.54 1.56 1.50 4.15 1.31 1.33 3.58 SCXD 1.90 1.91 6.80 1.69 1.88 4.90 1.46 1.75 4.30 SCYY 1.96 1.49 6.76 1.69 1.46 4.37 1.50 1.36 3.81 XIAG 3.15 4.71 9.45 3.05 4.55 8.14 2.95 4.37 7.04 YNCX 2.07 2.12 8.04 1.81 2.06 5.75 1.42 1.89 5.14 YNDC 2.01 2.44 8.91 1.89 2.36 7.19 1.75 2.13 6.12 YNLC 2.10 3.27 9.96 1.96 3.20 7.31 1.67 2.95 5.88 YNLJ 2.03 1.66 7.38 1.68 1.55 4.83 1.25 1.42 4.25 YNMJ 1.99 2.09 9.39 1.87 2.04 6.76 1.60 1.87 5.48 YNMZ 2.49 2.56 8.50 2.43 2.46 7.08 2.06 2.29 6.41 YNRL 2.32 2.35 8.62 2.10 2.24 5.99 1.74 1.87 5.11 YNSD 2.30 2.36 9.39 2.09 2.31 6.11 1.72 2.14 5.74 YNTC 1.80 1.97 10.27 1.62 1.94 6.67 1.39 1.62 6.46 YNTH 1.94 2.05 7.22 1.77 1.94 5.37 1.67 1.59 4.42 YNYA 1.77 1.90 8.20 1.60 1.81 5.62 1.42 1.72 4.41 YNYM 2.06 2.48 8.45 1.95 2.38 6.24 1.67 2.23 5.56 YNYL 2.09 2.35 11.57 1.79 2.31 9.08 1.58 2.10 5.89 YNYS 2.06 1.97 8.29 1.58 1.95 5.45 1.33 1.81 4.83 YNZD 3.53 3.66 7.20 3.38 3.54 5.70 3.30 3.42 4.51 平均值 2.13 2.34 8.29 1.94 2.26 5.95 1.70 2.07 5.10 3.3　周期信号对流动观测的影响 2 通过分析不同GNSS连续站时间序列，可以发现EMD方法能够根据信号的自身特性，自适应地对非构造变形的周期信号进行定量分离，刻画出更多细节变化，具有更高的应用价值。 由于流动GNSS观测采用漫游式观测，每期观测4 d，间隔1~3 a，观测时间不固定，需要耗费大量的人力物力，但取得的数据量较少，导致速度场求解精度不高，且部分区域N、U方向周期性明显，特别是U方向影响较大，所以周期项改正十分必要。目前文献［ 27 27 - 31 31 27-31 ］提出利用流动GNSS测站周边的多个连续站修正流动站的周期项影响，并取得了较好成果，但其采用的多为谐波模型。本文利用EMD方法对周期项进行修正，并进一步探讨利用改正后的连续站时间序列模拟流动观测，在积累多少期观测资料的情 况下才能获得比较准确的运动速度，以便为流动GNSS观测的周期修正和实施提供参考意见和理论依据。 选取研究区内8个连续站三分量时间序列，分别进行最小二乘拟合作为参考速度。假设每年随机选取连续4天作为一期流动观测资料，查看EMD方法需要多少期观测修正才能够达到与参考速度95%的差值在N、E方向小于1 mm/a，在U方向小于3 mm/a。按照上述要求，模拟了10 000次原始时间序列、谐波模型和EMD方法修正后拟合速度，统计结果见 表3 表3 。在N、E、U方向，谐波模型拟合速度的平均标准差减小比例分别为11.16%、5.63%、37.54%，EMD方法分别为24.33%、17.88%、47.15%，说明EMD方法比谐波模型对3个方向的拟合速度均具有更好的修正效果，特别是对U方向偏差大于3 mm/a的拟合速度。 3 站点 3 方向 4 原始数据 4 谐波改正 4 EMD改正 2 2 标准差减小比例/% 3 速度统计量 2 标准差/(mm·a -1 ) 3 速度统计量 2 标准差/(mm·a -1 ) 3 速度统计量 2 标准差/(mm·a -1 ) < 1 mm·a -1 [ 1 1 ， 3 3 ] mm·a -1 > 3 mm·a -1 < 1 mm·a -1 [ 1 1 ， 3 3 ] mm·a -1 > 3 mm·a -1 < 1 mm·a -1 [ 1 1 ， 3 3 ] mm·a -1 > 3 mm·a -1 谐波 EMD 3 SCMB N 9 529 471 0 0.49 9 823 177 0 0.43 9 838 162 0 0.42 13.09 13.89 E 9 125 875 0 0.57 9 250 750 0 0.55 9 739 261 0 0.46 3.48 18.90 U 5 048 4 590 362 1.45 6 835 3 128 37 0.99 7 756 2 235 9 0.83 31.62 42.60 3 SCXC N 9 226 774 0 0.59 9 490 510 0 0.55 9 880 120 0 0.47 7.42 20.12 E 9 469 531 0 0.52 9 826 174 0 0.42 9 956 44 0 0.35 20.27 32.76 U 3 803 4 865 1 332 1.97 6 385 3 571 44 1.07 7 631 2 364 5 0.82 45.68 58.52 3 SCXD N 9 245 755 0 0.56 9 626 374 0 0.49 9 919 81 0 0.39 13.15 30.56 E 9 093 907 0 0.61 9 187 813 0 0.59 9 510 490 0 0.54 2.30 10.59 U 3 853 4 974 1 173 1.95 5 738 4 097 165 1.25 6 501 3 462 37 1.06 35.73 45.75 3 SCYY N 8 786 1214 0 0.63 9 493 507 0 0.54 9 620 380 0 0.48 15.68 24.57 E 9 881 119 0 0.39 9 915 85 0 0.38 9 971 29 0 0.34 4.44 13.79 U 3 577 4 864 1 559 2.08 6 146 3 754 100 1.16 7 192 2 762 46 0.94 44.29 54.52 3 YNMJ N 9 387 613 0 0.53 9 573 427 0 0.49 9 889 111 0 0.41 7.98 22.33 E 9 306 694 0 0.55 9 410 590 0 0.53 9 682 318 0 0.46 4.17 16.53 U 2 599 4 423 2 978 2.84 3 972 4 722 1 306 1.98 4 181 5 251 568 1.61 30.40 43.26 3 YNSD N 8 878 1122 0 0.63 9 301 699 0 0.55 9 759 241 0 0.45 12.07 28.96 E 9 145 855 0 0.58 9 251 749 0 0.56 9 599 401 0 0.47 3.53 18.56 U 2 812 4 352 2 836 2.76 4 569 4 815 616 1.61 4 760 4 744 496 1.51 41.54 45.26 3 YNTH N 9 405 595 0 0.53 9 634 366 0 0.48 9 791 209 0 0.45 11.49 27.43 E 9 398 602 0 0.54 9 574 426 0 0.50 9 985 15 0 0.30 1.53 16.26 U 3 542 5 114 1 344 2.01 5 644 4 189 167 1.27 7 702 2 280 18 0.85 37.24 43.05 3 YNYM N 9 808 192 0 0.43 9 879 121 0 0.39 9 980 20 0 0.31 8.39 26.80 E 9 299 701 0 0.56 9 438 562 0 0.53 9 723 277 0 0.47 5.35 15.65 U 4 116 4 907 977 1.82 6 100 3 707 193 1.20 6 927 3 003 70 1.01 33.78 44.27 表4 表4 为利用原始时间序列和EMD方法改正两种方案进行周期项后，模拟流动观测满足上述条件要求时所需的观测期数。对比 表3 表3 可以发现，周期项对观测资料的影响不可忽略，有必要对其进行修正，周期项修正不仅可以提高数据的实际应用价值，在相同目标要求下还可以减少外业工作期数。EMD方法相比谐波模型对周期项的修正效能有较大的提升，周期性越明显，改进效果越好，这一点也可以通过3个方向结果进行对比验证，即U方向改正效果最好，N方向次之，E方向最差。 2 站点 2 观测期数 原始数据 EMD修正 SCMB 8 6 SCXC 6 5 SCXD 6 5 SCYY 6 5 YNMJ 6 5 YNSD 8 5 YNTH 7 5 YNYM 7 6 3.4　实际流动观测资料分析 2 如 图4 图4 所示，所选研究区内有2对相距较近的连续站和流动站，满足修正后理论上所需最少流动观测期数的观测数据，本文对其开展实际对比分析，分别利用SCMB和SCYY连续站资料对JB36和JB37流动站三分量进行周期项修正，连续站-流动站观测信息以及周期项修正后的拟合速度和RMSE见 表5 表5 。 2 连续站-流动站 2 方向 2 台站信息 4 拟合速度/(mm·a -1 ) 3 流动站RMSE/mm 距离/km 流动观测年份 连续站 流动站 谐波修正 EMD修正 原始序列 谐波修正 EMD修正 3 SCMB-JB36 N 3 <0.8 3 2011，2013—2015，2017，2018 -9.17 -10.09 -10.00 -9.87 21.31 21.07 20.98 E 34.61 35.32 35.35 35.05 75.27 75.31 75.00 U -0.04 -1.01 -1.36 -1.06 6.46 7.14 6.32 3 SCYY-JB37 N 3 <1.4 3 2011，2013—2018 -18.00 -17.82 -17.92 -17.86 48.01 48.29 48.02 E 38.11 37.60 37.63 37.51 103.50 103.59 103.29 U -0.25 2.73 2.46 2.74 10.58 9.03 8.61 图8 图8 为两对连续站-流动站周期项修正前后的时间序列及其线性拟合速度直线。由 图8 图8 可以看出，EMD方法修正后流动时间序列分布更加集中，但由于两个连续站和流动观测零值点起始年份分别为2010年和1998年，以及观测站址、环境条件、构造活动剧烈等因素综合影响，虽然两组台站相距很近，但拟合速度依然有明显的差别，特别是SCYY组垂直方向有着正负相反的变化趋势。从RMSE值看，谐波模型无论垂直方向还是水平方向对流动站的改正效果都不理想，修正效果不稳定，甚至越来越差，而EMD方法除JB37站N方向有0.01 mm升高外，其他值均有一定的改正效果，特别是垂直方向，说明该方法对非构造信息的去除更加稳定可靠。从整体上看，在水平方向，流动站和连续站均能很好地反映出本区域的构造运动趋势，但是由于其他构造和非构造因素的综合影响，即使相距很近的两个台站运动速度也可能有一定差异，数据长度足够时，周期信号对结果的影响并不太大；在垂直方向，周期项的影响很大，无论连续站还是流动站，都十分有必要对其进行修正，在同一区域不同台站的主要周期信号一般具有较好的一致性，但即使相距很近，其运动趋势也可能存在很大的差异，说明垂向运动不仅与区域整体运动有关，局部影响因素也很关键。 4　结　语 1 通过模拟实验与实际GNSS数据的坐标时间序列分析证明，EMD方法可以有效提取出非构造因素引起的周期项变化，其成分、幅值和周期变化更符合实际观测情况，蕴含丰富的非构造运动信息和物理意义，在时间序列分析中具有一定的优越性，相比于传统谐波改正模型也是一种方法上的新应用尝试。此外，本文探索利用EMD方法对连续观测资料进行周期项修正后模拟流动观测，查看需要多少期观测才能获取比较可靠的台站运动速度，以便为不同地区流动观测的实施和资料修正提供一些参考依据。得到的结论如下： 1）川滇地区不同台站模拟流动观测拟合速度经EMD方法修正后，最少5期、最多6期即可达到目标要求，如果只对水平方向修正，所需的期数会更少； 2）川滇地区GNSS连续站U方向运动周期性明显，大部分台站N方向周期性强于E方向，但也存在个别不同情况； 3）数据离散程度越大，所需观测期数越多，周期性越强，修正效果越好； 4）台站运动周期项的频率、幅值及其每年最大值出现的时间并不是固定不变的，可能与每年具体影响因素的强弱有关； 5）不同台站运动的时间序列不仅存在0.5 a、1.0 a，还可能还存在0.9 a、1.1 a、2.0 a等以往方法很难探测出的周期信号 ［ 9 9 ， 15 15 ， 24 24 ］ ； 6）本文10 000次统计模拟数据为每年随机选取，不排除包含相邻两期数据一个在年尾一个在年初的极端情况，如果近似均匀选取，所得结果可能会更好，在实际野外流动观测作业时，为了减弱周期项的影响，建议每期观测选择在相同的月份或季节； 7）在利用连续观测资料对流动观测资料进行周期项修正时，建议综合周围其他连续观测资料，选取最新时间序列变化平稳的时段资料进行周期项提取修正工作，以消除诸如大震、局部场地变化、仪器维修更换等引起的阶跃和趋势项转折等变化性影响； 8）EMD方法在周期项修正上明显优于谐波模型，并且不需给出先验拟合模型，它会根据信号的自身特性将信号分解为不同尺度的IMF分量，各分量所蕴含的具体物理意义还需进一步研究； 9）利用相关系数方法确定IMF分量噪声和信号分界点时，可能会存在临近一个高频分量划分不准的问题，这在后续工作中需要进一步完善，对于复杂信号趋势项和周期信号区分问题，可以通过FFT来确定； 10）EMD方法为周期项提取和修正提供了一种新的思路和途径，但也存在端点效应和模态混叠问题，在实际资料处理应用时还有很大改进空间。 参考文献 Jiang Weiping ， Li Zhao ， Liu Wanke ， et al . Some Thoughts on Establishment and Maintenance of Terrestrial Reference Frame Considering Non-Linear Variation ［J］. Geomatics and Information Science of Wuhan University ， 2010 ， 35 （ 6 ）： 665 - 669 . Some Thoughts on Establishment and Maintenance of Terrestrial Reference Frame Considering Non-Linear Variation 姜卫平 ， 李昭 ， 刘万科 ， 等 . 顾及非线性变化的地球参考框架建立与维持的思考 ［J］. 武汉大学学报（信息科学版） ， 2010 ， 35 （ 6 ）： 665 - 669 . 顾及非线性变化的地球参考框架建立与维持的思考 Wang Min ， Shen Zhengkang . Present-Day Tectonic Deformation in Continental China： Thirty Years of GPS Observation and Research ［J］. Earthquake Research in China ， 2020 ， 36 （ 4 ）： 660 - 683 . Present-Day Tectonic Deformation in Continental China： Thirty Years of GPS Observation and Research 王敏 ， 沈正康 . 中国大陆现今构造变形： 三十年的GPS观测与研究 ［J］. 中国地震 ， 2020 ， 36 （ 4 ）： 660 - 683 . 中国大陆现今构造变形： 三十年的GPS观测与研究 Dong D ， Fang P ， Bock Y ， et al . Anatomy of Apparent Seasonal Variations from GPS‐Derived Site Position Time Series ［J］. Journal of Geophysical Research ， 2001 ， 107 ： 2075 . Anatomy of Apparent Seasonal Variations from GPS‐Derived Site Position Time Series Jiang Weiping ， Li Zhao ， Liu Hongfei ， et al . Cause Analysis of the Non-Linear Variation of the IGS Reference Station Coordinate Time Series Inside China ［J］. Chinese Journal of Geophysics ， 2013 ， 56 （ 7 ）： 2228 - 2237 . Cause Analysis of the Non-Linear Variation of the IGS Reference Station Coordinate Time Series Inside China 姜卫平 ， 李昭 ， 刘鸿飞 ， 等 . 中国区域IGS基准站坐标时间序列非线性变化的成因分析 ［J］. 地球物理学报 ， 2013 ， 56 （ 7 ）： 2228 - 2237 . 中国区域IGS基准站坐标时间序列非线性变化的成因分析 Ming Feng ， Yang Yuanxi ， Zeng Anmin ， et al . Analysis of Seasonal Signals and Long-Term Trends in the Height Time Series of IGS Sites in China ［J］. Sience China Earth Sciences ， 2016 ， 46 （ 6 ）： 834 - 844 . Analysis of Seasonal Signals and Long-Term Trends in the Height Time Series of IGS Sites in China 明锋 ， 杨元喜 ， 曾安敏 ， 等 . 中国区域IGS站高程时间序列季节性信号及长期趋势分析 ［J］. 中国科学： 地球科学 ， 2016 ， 46 （ 6 ）： 834 - 844 . 中国区域IGS站高程时间序列季节性信号及长期趋势分析 Jiang W P ， Yuan P ， Chen H ， et al . Annual Variations of Monsoon and Drought Detected by GPS： A Case Study in Yunnan， China ［J］. Scientific Reports ， 2017 ， 7 （ 1 ）： 5874 Annual Variations of Monsoon and Drought Detected by GPS： A Case Study in Yunnan， China Zhan Wei ， Li Jingwei . Extraction and Analysis of the Common-Mode Component of GNSS Time Series in Yunnan ［J］. Journal of Seismological Research ， 2021 ， 44 （ 1 ）： 56 - 63 . Extraction and Analysis of the Common-Mode Component of GNSS Time Series in Yunnan 占伟 ， 李经纬 . 云南GNSS时间序列共模分量提取分析 ［J］. 地震研究 ， 2021 ， 44 （ 1 ）： 56 - 63 . 云南GNSS时间序列共模分量提取分析 Zhang Shuangcheng ， He Yuefan ， Li Zhenyu ， et al . EMD for Noise Reduction of GPS Time Series ［J］. Journal of Geodesy and Geodynamics ， 2017 ， 37 （ 12 ）： 1248 - 1252 . EMD for Noise Reduction of GPS Time Series 张双成 ， 何月帆 ， 李振宇 ， 等 . EMD用于GPS时间序列降噪分析 ［J］. 大地测量与地球动力学 ， 2017 ， 37 （ 12 ）： 1248 - 1252 . EMD用于GPS时间序列降噪分析 Davis J L ， Wernicke B P ， Tamisiea M E . On Seasonal Signals in Geodetic Time Series ［J］. Journal of Geophysical Research： Solid Earth ， 2012 ， 117 ： B01403 . On Seasonal Signals in Geodetic Time Series Chen Q ， Dam T V ， Sneeuw N ， et al . Singular Spectrum Analysis for Modeling Seasonal Signals from GPS Time Series ［J］. Journal of Geodynamics ， 2013 ， 72 ： 25 - 35 . Singular Spectrum Analysis for Modeling Seasonal Signals from GPS Time Series Zhou Maosheng ， Guo Jinyun ， Shen Yi ， et al . Extraction of Common Mode Errors of GNSS Coordinate Time Series Based on Multi-Channel Singular Spectrum Analysis ［J］. Chinese Journal of Geophysics ， 2018 ， 61 （ 11 ）： 4383 - 4395 . Extraction of Common Mode Errors of GNSS Coordinate Time Series Based on Multi-Channel Singular Spectrum Analysis 周茂盛 ， 郭金运 ， 沈毅 ， 等 . 基于多通道奇异谱分析的GNSS坐标时间序列共模误差的提取 ［J］. 地球物理学报 ， 2018 ， 61 （ 11 ）： 4383 - 4395 . 基于多通道奇异谱分析的GNSS坐标时间序列共模误差的提取 Gruszczynska M ， Rosat S ， Klos A ， et al . MultiChannel Singular Spectrum Analysis in the Estimates of Common Environmental Effects Affecting GPS Observations ［J］. Pure and Applied Geophysics ， 2018 ， 175 （ 5 ）： 1805 - 1822 . MultiChannel Singular Spectrum Analysis in the Estimates of Common Environmental Effects Affecting GPS Observations Zhang Jing ， Wang Xiaoya ， Hu Xiaogong . Analysis of GPS Stations' Time Series Based on PCA Method ［J］. Journal of Geodesy and Geodynamics ， 2019 ， 39 （ 6 ）： 613 - 619 . Analysis of GPS Stations' Time Series Based on PCA Method 张晶 ， 王小亚 ， 胡小工 . 基于PCA方法的GPS台站时间序列分析 ［J］. 大地测量与地球动力学 ， 2019 ， 39 （ 6 ）： 613 - 619 . 基于PCA方法的GPS台站时间序列分析 Huang N E ， Shen Z ， Long S R ， et al . The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis ［J］. Proceedings： Mathematical， Physical and Engineering Sciences ， 1998 ， 454 （ 1971 ）： 903 - 995 . The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis Cui Chengeng . Application of a Combined Method of Filtering to GPS Dynamic Observation Retraction ［J］. Journal of Gansu Sciences ， 2020 ， 32 （ 3 ）： 40 - 45 . Application of a Combined Method of Filtering to GPS Dynamic Observation Retraction 崔晨耕 . 一种组合滤波法在GPS动态观测中的应用 ［J］. 甘肃科学学报 ， 2020 ， 32 （ 3 ）： 40 - 45 . 一种组合滤波法在GPS动态观测中的应用 Chen C H ， Yeh T K ， Liu J Y ， et al . Surface Deformation and Seismic Rebound： Implications and Applications ［J］. Surveys in Geophysics ， 2011 ， 32 （ 3 ）： 291 - 313 . Surface Deformation and Seismic Rebound： Implications and Applications Chen C H ， Wen S ， Yeh T K ， et al . Observation of Surface Displacements from GPS Analyses Before and After the Jiashian Earthquake （M=6.4） in Taiwan， China ［J］. Journal of Asian Earth Sciences ， 2013 ， 62 ： 662 - 671 . Observation of Surface Displacements from GPS Analyses Before and After the Jiashian Earthquake （M=6.4） in Taiwan， China Chen C H ， Wang C H ， Wen S ， et al . Anomalous Frequency Characteristics of Groundwater Level Before Major Earthquakes in Taiwan， China ［J］. Hydrology and Earth System Sciences ， 2013 ， 17 （ 5 ）： 1693 - 1703 . Anomalous Frequency Characteristics of Groundwater Level Before Major Earthquakes in Taiwan， China Feng W ， Ren J W ， Jiang Z S . GPS Station Short-Term Dynamic Characteristics of Micro Displacement Before Menyuan M 6.4 Earthquake ［J］. Geo⁃desy and Geodynamics ， 2016 ， 7 （ 4 ）： 237 - 244 . GPS Station Short-Term Dynamic Characteristics of Micro Displacement Before Menyuan M 6.4 Earthquake Jia Ruisheng ， Zhao Tongbin ， Sun Hongmei ， et al . Micro-Seismic Signal Denoising Method Based on Empirical Mode Decomposition and Independent Component Analysis ［J］. Chinese Journal of Geophysics ， 2015 ， 58 （ 3 ）： 1013 - 1023 . Micro-Seismic Signal Denoising Method Based on Empirical Mode Decomposition and Independent Component Analysis 贾瑞生 ， 赵同彬 ， 孙红梅 ， 等 . 基于经验模态分解及独立成分分析的微震信号降噪方法 ［J］. 地球物理学报 ， 2015 ， 58 （ 3 ）： 1013 - 1023 . 基于经验模态分解及独立成分分析的微震信号降噪方法 Luo Feixue ， Dai Wujiao ， Wu Xixiu . EMD Filtering Based on Cross-Validation and Its Application in GPS Multipath ［J］. Geomatics and Information Science of Wuhan University ， 2012 ， 37 （ 4 ）： 450 - 453 . EMD Filtering Based on Cross-Validation and Its Application in GPS Multipath 罗飞雪 ， 戴吾蛟 ， 伍锡锈 . 基于交叉证认的EMD滤波及其在GPS多路径效应中的应用 ［J］. 武汉大学学报（信息科学版） ， 2012 ， 37 （ 4 ）： 450 - 453 . 基于交叉证认的EMD滤波及其在GPS多路径效应中的应用 Dai Wujiao ， Ding Xiaoli ， Zhu Jianjun ， et al . EMD Filter Method and Its Application in GPS Multipath ［J］. Acta Geodaetica et Cartographica Sinica ， 2006 ， 35 （ 4 ）： 321 - 327 . EMD Filter Method and Its Application in GPS Multipath 戴吾蛟 ， 丁晓利 ， 朱建军 ， 等 . 基于经验模式分解的滤波去噪法及其在GPS多路径效应中的应用 ［J］. 测绘学报 ， 2006 ， 35 （ 4 ）： 321 - 327 . 基于经验模式分解的滤波去噪法及其在GPS多路径效应中的应用 Lu Tieding ， Qian Wenlong ， He Xiaoxing ， et al . An Improved EMD Method for Determining Boundary IMF ［J］. Journal of Geodesy and Geodynamics ， 2020 ， 40 （ 7 ）： 720 - 725 . An Improved EMD Method for Determining Boundary IMF 鲁铁定 ， 钱文龙 ， 贺小星 ， 等 . 一种确定分界IMF分量的改进EMD方法 ［J］. 大地测量与地球动力学 ， 2020 ， 40 （ 7 ）： 720 - 725 . 一种确定分界IMF分量的改进EMD方法 Guo Xiang . Research on GPS Date Denoising Based on EMD ［D］. Fuzhou ： East China University of Technology ， 2016 . Research on GPS Date Denoising Based on EMD 郭翔 . EMD在GPS信号去噪中的应用研究 ［D］. 抚州 ： 东华理工大学 ， 2016 . EMD在GPS信号去噪中的应用研究 Herring T A ， King R W ， McClusky S C . GAMIT Reference Manual： GPS Analysis at MIT， Version 10.4 ［R］. Cambridge ： Massachusetts Institute of Technology ， 2010 . GAMIT Reference Manual： GPS Analysis at MIT， Version 10.4 van Dam T ， Wahr J ， Lavallée D . A Comparison of Annual Vertical Crustal Displacements from GPS and Gravity Recovery and Climate Experiment （GRACE） over Europe ［J］. Journal of Geophysical Research ， 2007 ， 112 （ B3 ）： B03404 . A Comparison of Annual Vertical Crustal Displacements from GPS and Gravity Recovery and Climate Experiment （GRACE） over Europe Liu Renli ， Li Jiancheng ， Jiang Weiping ， et al . Comparing Vertical Surface Displacements Using GRACE and GPS over Shanxi Province ［J］. Geomatics and Information Science of Wuhan University ， 2013 ， 38 （ 4 ）： 426 - 430 . Comparing Vertical Surface Displacements Using GRACE and GPS over Shanxi Province 刘任莉 ， 李建成 ， 姜卫平 ， 等 . 联合GRACE与GPS比较山西省垂向地表形变 ［J］. 武汉大学学报（信息科学版） ， 2013 ， 38 （ 4 ）： 426 - 430 . 联合GRACE与GPS比较山西省垂向地表形变 Zhan Wei ， Li Fei ， Zhu Shuang ， et al . Analysis and Test of Correction to Vertical Velocity Measured by GPS Mobile Observation with Continuous Observation ［J］. Geomatics and Information Science of Wuhan University ， 2016 ， 41 （ 7 ）： 963 - 968 . Analysis and Test of Correction to Vertical Velocity Measured by GPS Mobile Observation with Continuous Observation 占伟 ， 李斐 ， 朱爽 ， 等 . 应用GPS连续观测修正流动观测垂向速率的分析与实验 ［J］. 武汉大学学报（信息科学版） ， 2016 ， 41 （ 7 ）： 963 - 968 . 应用GPS连续观测修正流动观测垂向速率的分析与实验 Zhan W ， Li F ， Hao W F ， et al . Regional Characteristics and Influencing Factors of Seasonal Vertical Crustal Motions in Yunnan， China ［J］. Geophysical Journal International ， 2017 ， 210 （ 3 ）： 1295 - 1304 . Regional Characteristics and Influencing Factors of Seasonal Vertical Crustal Motions in Yunnan， China Li Jingwei ， Zhan Wei ， Liang Hongbao ， et al . Accurate Acquisition and Method Evaluation of GPS Annual Vertical Motion ［J］. Journal of Seismological Research ， 2020 ， 43 （ 1 ）： 95 - 100 . Accurate Acquisition and Method Evaluation of GPS Annual Vertical Motion 李经纬 ， 占伟 ， 梁洪宝 ， 等 . GPS垂向年周期运动的精确获取与方法评估 ［J］. 地震研究 ， 2020 ， 43 （ 1 ）： 95 - 100 . GPS垂向年周期运动的精确获取与方法评估 Liang H B ， Zhan W ， Li J W . Vertical Surface Displacement of Mainland China from GPS Using the Multisurface Function Method ［J］. Advances in Space Research ， 2021 ， 68 （ 12 ）： 4898 - 4915 . Vertical Surface Displacement of Mainland China from GPS Using the Multisurface Function Method 感谢中国地震局第一监测中心提供的GNSS数据资料。 10.13203/j.whugis20210029 P228 A 1671-8860（2023）01-0135-11 刘希康,丁志峰,李媛，等.EMD在GNSS时间序列周期项处理中的应用[J].武汉大学学报（信息科学版）,2023,48(1):135-145.DOI: 10.13203/j.whugis20210029 http://dx.doi.org/10.13203/j.whugis20210029 刘希康,丁志峰,李媛，等.EMD在GNSS时间序列周期项处理中的应用[J].武汉大学学报（信息科学版）,2023,48(1):135-145.DOI:10.13203/j.whugis20210029 LIU Xikang,DING Zhifeng,LI Yuan,et al.Application of EMD to GNSS Time Series Periodic Term Processing[J].Geomatics and Information Science of Wuhan University,2023,48(1):135-145.DOI: 10.13203/j.whugis20210029 http://dx.doi.org/10.13203/j.whugis20210029 LIU Xikang,DING Zhifeng,LI Yuan,et al.Application of EMD to GNSS Time Series Periodic Term Processing[J].Geomatics and Information Science of Wuhan University,2023,48(1):135-145.DOI:10.13203/j.whugis20210029 国家重点研发计划（2018YFC1503606）；中国地震局青年震情跟踪课题（2020010226，2021010511）；国家自然科学基金（41472180）。 2021-01-19 2023-01-05 2023-06-29 33226536 武汉大学学报（信息科学版） 1 48