2018, Vol. 61
2. 山东科技大学矿山灾害预防控制省部共建国家重点实验室, 青岛 266590;
3. 同济大学测绘与地理信息学院, 上海 200092;
4. 中国科学院测量与地球物理研究所, 武汉 430077
2. State Key Laboratory of Mining Disaster Prevention and Control Co-founded by Shandong Province and the Ministry of Science and Technology, Shandong University of Science and Technology, Qingdao 266590, China;
3. College of Surveying and Geo-informatics, Tongji University, Shanghai 200092, China;
4. Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China
由固体潮观测和研究得到的地球物理信息是了解地球内部结构的重要依据，研究表明高精度重力固体潮观测是区别于地震技术用于检测地球内部动力学和反演地球内部构造的有效途径，为研究地球自由振荡、液核近周日共振、固态内核的平动振荡和地球自转等提供有效约束，是地震监测的重要技术，是解释诸如大气、海洋和陆地水等地表质量迁移以及各种区域环境效应的重要手段(管泽霖和宁津生, 1981; Crossley et al., 1999; Sun et al., 2001, 2002, 2004; 孙和平等, 2006).固体潮一直是大地测量学和地球物理学的基础研究，其研究方法主要包括仪器观测、理论模拟以及实际观测值与模拟值之间的对比分析.
19世纪末，固体潮主要采用水平摆来进行观测，但由于受到当时制造技术水平所限，其稳定性和观测精度较差，未对固体潮的研究起到关键性的推动作用.随着重力仪技术的兴起和发展，利用相对重力仪进行相对重力测量成为观测固体潮的主要手段(Ducarme and Sun, 2001; Sun et al., 2002).一个物体(质量为m)在地球重力场内的重量mg随着重力加速度g的变化而变化，如果用另一种力(如弹力)来平衡这种重量的变化，则通过对物体平衡状态的观测，就可以测量出两点间的重力差值.在进行相对重力测量时，相对重力仪除了受到重力的影响外，还会受到其他因素的影响，因此在计算时要加入相应的改正，主要的改正项包括:重力固体潮、零点漂移(陆仲连, 1996; Yu et al., 2015, 2018; 于红娟等, 2017)、温度改正、气压改正、倾斜改正(Torge, 1989)和海潮改正(Gokhan and Jin, 2016; Lei et al., 2017). CG-5相对重力仪是由加拿大的Scintrex公司继CG-3型重力仪之后设计生产的一种改进型自动电子读数相对重力仪，该仪器基于熔凝石英弹性系统，采用微处理器装置实现自动测量，在测量过程中，可以自动实现温度补偿、仪器倾斜改正、地震滤波、线性漂移改正等功能(Scintrex Limited, 2012)，读数分辨率达到1×10-8m·s-2.在相对重力测量中，重力固体潮也反映在观测结果中，把重力固体潮成分从相对重力观测数据中分离对相对重力观测具有重要意义.
大量相关学者从重力固体潮的物理机制出发，根据已有的重力观测资料，建立了具有较高精度的固体潮模型(Lefèvre et al., 2000; Sun et al., 2003; Penna et al., 2008; Taguchi et al., 2014).然而，理论模型是在一些理想条件和固定参数的基础上而建立的，由理论模型计算得到的重力固体潮值是理想化的，不能完全反映出实际观测中的变化异常，因此由重力固体潮模型计算的重力固体潮值存在一定的局限性.
奇异谱分析(Singular Spectrum Analysis，SSA)(Vautard et al., 1992)作为一种数字信号处理技术，在气候学(Wyatt et al., 2012)、测量学(Chen et al., 2013; Zabalza et al., 2014)、海洋学(Kondrashov and Berloff, 2015)等领域都有应用.对于事先未知物理本质的系统，SSA可以从它包含噪声的有限长度的观测序列中提取尽可能多的可靠信息，并基于这些信息对数据进行趋势识别、周期判断、降噪处理和数据重构等.SSA作为经验正交函数展开的一种扩展形式，主要用于单一时间序列的分析，它从时间序列的动力重构出发，不受正弦波假定的约束，无需先验信息，具有稳定识别和强化周期信号的优点，特别适合于分析有周期震荡的时间序列数据(吴洪宝和吴蕾, 2005).根据SSA能够提取观测序列中主成分的特性，利用SSA可以对时间序列进行粗差探测、插值和外推等操作(Shen et al., 2017, 2018).然而，SSA在进行信号提取时会产生边界效应，边界效应的存在会使得时间序列边缘处的数据处理精度较差，从而对整个时间序列的处理结果造成一定的影响.因此，需要用一定方法削弱SSA的边界效应.
小波变换是一种时间和频率的局域变换方法，能够从信号中提取有用的信息，并通过伸缩平移等功能对信号进行多尺度的细化分析.小波变换的核心思想就是按照尺度来分析信号，通过小波伸缩和平移来研究信号与小波之间的相关性.小波变换在科学研究和工程技术中得到广泛应用，其在大地测量学和地球物理学中也是进行信号去噪和信号分析的一种常用方法(杨文采等, 2001; 柳林涛和许厚泽, 2004; Hu et al., 2006, 2007).通过SSA方法与小波变换方法的对比可以对SSA方法的可靠性进行评价.
为了将重力固体潮从相对重力观测数据中分离出来，在不需要测站位置、理论潮波信息等先验信息的基础上，本文提出一种利用SSA提取相对重力观测数据中重力固体潮影响量的新思路.
1 原理与方法 1.1 奇异谱分析设x1, x2, …, xN为观测得到的一维相对重力时间序列，其时间长度为N，SSA可以把该时间序列分解成诸如趋势、周期或准周期、噪声等一系列具有各自特征的时序，以达到对原始时序的趋势提取、周期识别、数据平滑和消噪等目的(Elsner, 2002; Shen et al., 2018).下面将根据Golyandina和Zhigljavsky(2013)等人的研究方法，给出SSA的主要计算步骤：
(1) 构建时滞矩阵.设整数
|
(1) |
其中Xi=(xi, …, xi+M-1)T, (1≤i≤K). X为M×K阶时滞矩阵，它的每一行每一列都是初始时间序列的子序列，且其为Hankel矩阵.
(2) 奇异值分解.令矩阵S = XXT，求得S的特征值并按照大小降序排列(λ1≥…≥λM≥0)，对应的特征向量分别为U1, …, UM.设d=rank(X)，Vk= 
|
(2) |
其中初等矩阵
(3) 对角平均化.对角平均化的目的就是把上步中分解得到的初等矩阵X k 重新转换成为长度为N的新时间序列，称之为重建成分(Reconstruction component, RC)，所有RC之和等于原序列.定义矩阵Z = Xk，z1, z2, …, zN为Z经过对角平均化得到的时间序列.设M*=min(M, K), K*=max(M, K)，如果M＜K，则让zij*=zij，否则就让zij*=zji，对角平均化的公式可表示为：
|
(3) |
长期重力监测数据中通常包含着年周期、半年周期、长期趋势变化和噪声等，利用w-correlation方法(Hassani, 2007)对各RC之间的相关性进行分析，将具有相同信号特征的RC进行分组.假设得到的时间序列为Yi，则任意两个重建时间序列的w-correlation可表示为：
|
(4) |
其中
奇异谱粗差探测的基本原理是利用SSA提取原时间序列中的包括趋势项、周期项在内的主成分，通过判断原时间序列与主成分序列之间的残差值大小来辨别粗差.其具体步骤如下：
(1) 将原时间序列进行SSA分解，利用w-correlation方法分析重建成分之间的相关性，设当重建成分阶次大于i时各个重建成分之间不能很好地分离；
(2) 合并前i阶的重建成分作为原时间序列的主成分；
(3) 对原时间序列与主成分序列的残差序列进行分析，依照3σ准则确定阈值，判断粗差并剔除.
1.3 奇异谱迭代插值对于因剔除粗差造成的数据缺失，或原时间序列中本就存在的数据缺失，利用SSA迭代插值来补全数据.SSA迭代插值(Schoellhamer, 2001; Beckers and Rixen, 2003; Kondrashov et al., 2010; Kondrashov and Ghil, 2006; Shen et al., 2017)的基本原理是将时间序列中有数据缺失的位置用0填充，再进行SSA分解，用分解得到的重建成分的值代替原来缺失位置的0值，反复循环此过程，以达到最好的效果.在计算过程中，窗口长度M取时间序列已知周期的公倍数，重构阶次P取i.选取原时间序列中的一段数据作为交叉验证数据，求取交叉验证数据与SSA重构序列之间的残差RMS值，将此RMS值的大小作为评价插值质量的标准.具体步骤为：
(1) 对原时间序列中的缺失数据用0补全并做标记，对补全后的时间序列做中心化处理；
(2) 对补全后的时间序列做奇异谱分解，确定前i阶RC为原序列的主要成分，原序列中缺失位置处的数值用第一个RC(记为R1)中的值代替，反复循环此过程，直至前后两次插补的数据残差RMS值小于0.1 μGal；
(3) 添加R2重建缺失数据，即由R1和R2线性叠加得到缺失数据，重复步骤(2)，直至R1…Ri均线性叠加到缺失数据中去.
1.4 SSA边界效应的削弱奇异谱分析通常存在较为明显的边界效应，为了减弱其影响，本文提出利用SSA在长度为N的原始相对重力时间序列两端分别外推n个数据.SSA迭代外推的基本原理与SSA插值的原理相同，SSA插值是将有数据缺失的时间序列用0补全，而SSA迭代外推是将要外推数据的时间序列两端分别补n个0，其余步骤与奇异谱迭代插值相似，不同的是SSA外推数据时重构阶次P参照重构成分的贡献率或其信号特征进行选择；在进行外推数据质量评价时，将整个原时间序列作为交叉验证数据进行质量评价.
将数据外推后得到的长度为N+2n的时间序列进行SSA分解、重构，得到包含重力固体潮成分的重构时间序列.对此重构时间序列的两端分别截掉n个数据得到最终的SSA重力固体潮提取值.
1.5 小波变换原理小波分析融合多时间尺度的理念，可以使时间序列中不易被发现的多种周期性变化被挖掘出来.小波变换具有良好的时频局部化性能，即在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，因而对信号在时频域都具有很强的表征信号局部特性的自适应性.
在小波变换的信号分析中，如果在每个可能的尺度下都计算小波系数，计算量将会非常的大.因此，为了数值计算的简化和理论分析的简便，需要对小波变换进行离散化处理(Stephane, 1998).离散小波不需要知道尺度函数和小波函数的具体结构，根据系数即可以实现信号的分解和重构，从而将原信号分解为具有不同信号特征的信号.具体步骤为：给定长度为N的相对重力序列s，第一步使用给定的小波产生两个系数集：低频(逼近)系数cA1和高频(细节)系数cD1；下一步分解就是使用同样的框架将低频系数cA1分解成两部分，即cA1替换为cA2和cD2，依次类推，信号在第j层的分解结构为[cAj, cD1, …, cDj]；最后通过给定的小波重构各个系数向量得到小波变换后的重构序列.
2 模拟数据分析由于重力固体潮实际值与相应理论值之间是有一定差异的，因此采用SSA从实测数据中提取得到的重力固体潮与理论值进行对比，并不能很好地说明SSA提取重力固体潮的可靠性.为了证明SSA方法提取重力固体潮的可靠性，先采用模拟的相对重力数据进行SSA分析并与小波分析方法进行对比.
2.1 模拟数据CG-5相对重力仪测得的相对重力数据可以分为三个部分：重力固体潮、海潮负荷、零漂和噪声，其中海潮负荷效应对重力的影响在实际观测中可以通过模型值进行改正，因此本实验模拟数据中只包含重力固体潮、零漂和噪声，分别采用以下方法进行模拟：①重力固体潮，应用国际地潮中心推荐的TSoft软件(Van Camp and Vauterin, 2005)采用Tamura潮波表中的Q1、O1、M1、P1、K1、J1、N2、M2、L2、S2、K2、M3主潮波计算实验地点(经纬度：36.00070°N，120.11930°E)2017-01-30—2017-02-08(UTC时间)时间长度为10 d的相对重力时间序列，采样间隔为60 s，共14400个数据；②零漂，CG-5相对重力仪的静态零点漂移线性度很好(于红娟等, 2017)，因此以斜率为-200 μGal/d的趋势项来模拟零漂；③噪声，加入标准差为1.3 μGal的高斯白噪声来模拟实测数据中的随机噪声(马玄龙等, 2010).
2.2 模拟数据SSA分析对模拟相对重力数据进行最小二乘拟合，得到10 d内的模拟数据的零漂值为-200.448 μGal/d，与实际值相差0.448 μGal/d，根据最小二乘拟合结果将零点漂移从模拟数据中去掉.对去掉零漂的模拟序列按照1.1节中介绍的SSA方法(选择窗口长度M=5760)进行分解和重构.前30阶重构成分w-correlation分析结果如图 1所示，根据SSA原理，当序列中存在一个周期成分时，SSA将其分解得到一对接近相等的特征值，即当两个相邻RC之间的相关系数大于一定的阈值，则判断属于同一周期成分.若阈值选择过大，会导致一些周期成分被过滤掉，反之则不能较好地排除一些周期性不明显的重构成分.经实验，当ρijw>0.86时，可以认定模拟数据中两个重构时间序列属于同一周期成分，并能够较好地将非周期成分排除，据此认为RC1、RC2、…、RC8、RC10、RC11、RC17、…、RC30均具有周期性.对重构成分进行快速傅里叶变换(Fast Fourier Transform，FFT)频谱分析后发现，第17阶及其之后的重构成分均为高频信号，且其振幅很小，可视为高频噪声进行剔除，部分高频信号及其频谱信息如图 2所示.
|
图 1 模拟数据SSA前30重构成分的w-correlation Fig. 1 The w-correlations of the first 30 reconstructed components of the simulated signal with SSA |
|
图 2 部分高频重构成分及其FFT功率波谱 Fig. 2 Part high frequency reconstructed components and the FFT power spectrums |
将属于同一周期成分的重构序列进行合并，利用FFT谱分析探寻合并项中存在的周期特性，如图 3所示.在利用SSA提取重力固体潮时，不能将各个潮波信号完全分离，具有相似频率的潮波会被划分成同一组信号，将FFT谱分析得到的各个合成项的频率与各主潮波的理论频率进行对比，进行大致的潮波划分.如表 1所示，RC1与RC2的合并项(简称RC1+2，下同)的频率为1.9336 cpd(cycles per day)，但其中也混杂着频率为1.0547 cpd的潮波信号；RC3+4的主频率为0.9668 cpd，其中也混杂了频率为1.9336 cpd的半日周期信号；而主频率为0.8789 cpd的RC5+6项则混入了更多的其他信号，包括频率为2.0215 cpd的周期信号；RC7+8和RC10+11项分离情况较好，基本只包含一个周期信号，频率分别为2.0215 cpd、2.9004 cpd.
|
图 3 模拟数据中由SSA提取的主要周期项及其FFT功率波谱 Fig. 3 Main periodic components extracted from simulated signal with SSA and the FFT power spectrums |
|
表 1 模拟数据中主要周期项的潮波划分 Table 1 Tide division of main periodic components of the simulated signal |
利用TSoft软件采用Tamura潮波表计算得到重力固体潮理论值，与SSA得到的重力固体潮提取值进行对比，结果如图 4所示，其中红色实线(图中SSA)为利用SSA提取的重力固体潮，黑色虚线(图中TSoft)为由TSoft计算的重力固体潮理论值.提取值与理论值的差异在对比图中不能清晰辨别，因此图 4中给出了二者的残差序列图，残差统计结果见表 2.图 4显示残差序列在左右边界处的波动较大，残差最小值-1.1 μGal与最大值0.9 μGal也分别出现在残差序列左右边界附近，说明SSA方法存在较为明显的边界效应.SSA提取值与理论值残差的RMS与STD值均为0.3 μGal，说明利用SSA方法获取相对重力信号中的重力固体潮时，虽然不能将重力固体潮的各个潮波进行很好的分离，但是可以将大部分重力固体潮信号从相对重力信号中很好地提取出来，从而可以对相对重力仪观测数据进行重力固体潮改正.
|
图 4 模拟数据SSA重力固体潮提取值与理论值对比与残差 Fig. 4 Residuals between the gravity tides extracted from the simulated signal with SSA and the theoretical values |
|
|
表 2 SSA和小波变换提取的重力固体潮残差统计(mGal) Table 2 Statistics of the gravity tides residuals obtained by SSA and wavelet transform (mGal) |
利用Myer小波对去除零漂的模拟相对重力信号s进行9层小波分解和重构，得到一个逼近信号A9和九个细节信号D1、D2、D3、…、D9，原信号s与重构得到的信号满足s=A9+D9+D8+D7+…+ D1.经过FFT频谱分析发现，D1—D7这七个细节信号属于高频信号，作为噪声进行剔除，A9、D9、D8则具有较为明显的周期特性，如图 5所示，A9具有1.0547 cpd的日周期，D9具有1.9336 cpd的半日周期，而D8中混入了一些噪声信号，但其明显具有2.9004 cpd的1/3日主周期.将A9、D9、D8合并项作为重力固体潮提取值与TSoft计算得到的理论值进行对比，如图 6所示，残差序列具有一定的线性趋势，并且在残差序列右边界处存在较大的边界效应，使得残差序列最大值达到12.7 μGal.由于残差中存在线性趋势和边界效应的影响，其STD与RMS值大于1 μGal量级.
|
图 5 模拟数据中由小波变换获取的主要周期项及其FFT功率波谱 Fig. 5 Main periodic components obtained from the simulated signal with wavelet transform and the FFT power spectrums |
|
图 6 模拟数据小波变换重力固体潮提取值与理论值残差 Fig. 6 Residuals between the gravity tides obtained from the simulated signal with wavelet transform and the theoretical values |
经以上实验可以发现，SSA方法与小波方法在获取相对重力信号中的重力固体潮时均受到边界效应的影响，尤其以小波分析更为明显，由于边界效应的影响使得其重力固体潮残差最大值达到了12.7 μGal.从表 2的残差统计数据中可以看出，由SSA方法和小波分析得到的重力固体潮残差最大值与最小值之间差值分别为2.0 μGal、17.6 μGal，可见利用SSA方法获取的残差序列波动更小.小波变换得到的残差结果还存在一定的线性趋势，而SSA得到的相对重力残差序列较为平稳.两种方法得到的重力固体潮与理论值的残差RMS值分别为0.3 μGal、1.6 μGal，所以在本实验中SSA方法比Myer小波变换具有更强的可靠性.
3 实例分析为了进一步研究利用SSA提取重力固体潮的可行性和可靠性，以及SSA边界效应对重力固体潮提取的影响，本文按照图 7中的技术路线进行实验.
|
图 7 SSA提取重力固体潮技术路线 Fig. 7 Sketch of extracting gravity tides with SSA |
将序列号为140541221的CG-5相对重力仪置于山东科技大学内一固定点进行连续的静态相对重力观测.读数自动记录时间间隔为1 min，采样频率为6 Hz，1 min内(其中55 s进行数据采集，5 s进行数据处理)可获得330个样本数据，利用仪器自身具有的地震滤波功能对这些样本数据中的低频噪声进行过滤，并舍弃高于6倍标准偏差的高频噪声，然后对样本数据进行平均得到1 min观测时间内的最终读数.观测过程中仪器还进行了温度补偿和倾斜改正(Scintrex Limited, 2012).通过17000 min的静态相对重力观测共得到17000个相对重力观测数据，中间的14400个数据是本文研究的主要对象.为了消除边界效应，要对原始序列进行SSA迭代外推，观测开始和结束各1300个数据作寻求最佳的数据外推个数之用.
根据1.2节中的方法对由14400个相对重力数据组成的时间序列进行粗差探测.对时间序列进行SSA并通过w-correlation方法分析各RC之间的相关系数，发现当RC阶次大于10的时候各RC之间就不能很好地分离，合并前10阶RC作为原时间序列的主成分，其与原相对重力时间序列的残差见图 8，根据3σ准则选择阈值区间为(-3 μGal, 3 μGal)，残差值超出此区间的相应相对重力数据被认为是粗差，据此共剔除了117个粗差.对于剔除粗差后的相对重力时间序列，使用1.3节中的SSA迭代插值方法进行数据插补.选取窗口长度为1440，由于已知前10阶RC为原始相对重力时间序列的主成分，因此选择重构阶次为10，此时交叉验证数据与重构时间序列之间的残差RMS为0.86 μGal，插补后的相对重力时间序列与主成分序列之间的残差如图 9所示.残差序列较为平稳，说明粗差剔除与数据插补的效果良好.
|
图 8 原始相对重力时间序列与主成分序列残差 (虚线表示阈值) Fig. 8 Residuals between original relative gravity time series and principal component time series (The dashed lines represent the thresholds) |
|
图 9 SSA粗差剔除和插补后的相对重力时间序列与主成分序列残差 Fig. 9 Residuals of the relative gravity time series after gross errors were eliminated and the missing data are interpolated with SSA relative to the principal component time series |
由于海潮和重力固体潮都是由日月等引力变化引起的，二者具有相同的频率，利用SSA不能将二者区分开来(许厚泽等, 1982)，只能采用表面负荷理论(Farrell, 1972)和已有的海潮模型(Matsumoto et al., 1995; Lefèvre et al., 2002)，通过表面负荷Green函数与海潮潮高的全球褶积从理论上计算出海潮的负荷影响(Agnew, 1997).因此在利用SSA进行重力固体潮提取之前，本文基于NAO99b海潮模型(Matsumoto et al., 2000)，应用国际地潮中心推荐的TSoft软件(Van Camp and Vauterin, 2005)计算获取本实验地点相应历元下的重力海潮负荷值时间序列.将海潮对相对重力的影响从原时间序列中扣除再利用SSA提取重力固体潮.
3.2 静态零漂在利用SSA提取相对重力数据中的重力固体潮之前，需要将其中的零漂项进行剔除，基于SSA可以将信号中的趋势项进行识别的特性，采用SSA的方法进行零漂剔除.CG-5相对重力仪的静态零点漂移线性度很好(于红娟等，2017)，采用SSA方法需要识别相对重力数据中的线性趋势，因此选择的窗口长度不宜过长，否则识别的趋势项会带有一定的非线性趋势，这与相对重力仪的零漂特性不符.按照1.1节的方法采用窗口长度M=1440对相对重力时间序列进行分解和重构，得到的RC1对应的方差贡献率接近于1，代表相对重力时间序列的总体线性趋势，将其进行最小二乘拟合作为提取的零点漂移.为了说明SSA提取的零点漂移的可靠性，对相对重力时间序列直接进行最小二乘拟合求出10 d的零漂值与SSA求得的估计值作比较.由SSA估算得到的零漂值为229.450 μGal/d，对相对重力时间序列进行最小二乘拟合求得的零漂值为229.061 μGal/d，二者之间的差值仅为0.389 μGal/d，具有非常好的一致性.由于利用SSA提取零漂时剔除了噪声，所以效果更好.
3.3 重力固体潮的提取选择窗口长度M=5760对外推后的相对重力时间序列进行SSA分解和重构，对重建成分进行w-correlation分析，图 10显示了前30项由Yi得到的重建成分之间的ρijw值，当i>14时，各RC之间就不能很好的相互分离，这说明它们中噪声占较大部分.若ρijw>0.85，则认为两个重构时间序列属于同一周期成分，据此将RC1与RC2，RC3与RC4，RC5与RC6，RC7与RC8、RC12与RC13分别合并为一组.利用FFT确定分组后的时间序列的周期性，如图 11所示，每个合并项的频率信息见表 3.
|
图 10 前30重建部分的w-correlation Fig. 10 The w-correlations of the first 30 reconstructed components |
|
图 11 主要周期项及其FFT功率波谱 Fig. 11 Main periodic components and the FFT power spectrums |
|
|
表 3 主要周期项的潮波划分 Table 3 Tide division of main periodic components |
将SSA提取的周期项进行合并，合并的结果视为重力固体潮的提取值.利用TSoft软件采用Tamura潮波表中的Q1、O1、M1、P1、S1、K1、J1、OO1、2N2、N2、M2、L2、S2、K2、M3、M4潮波计算得到重力固体潮理论值，与SSA得到的重力固体潮提取值进行对比，残差如图 12中虚线(图中SSA)所示，统计结果见表 4.SSA提取值与理论值残差的RMS与STD值均不超过2.5 μGal，说明由SSA提取的重力固体潮是比较可靠的.
|
图 12 残差对比 Fig. 12 Comparison of different residuals |
|
|
表 4 不同处理方式下的重力固体潮残差统计(mGal) Table 4 Statistics of the gravity tides residuals obtained by different processes (mGal) |
由2.2节模拟数据的SSA实验可知，SSA在提取相对重力信号中的重力固体潮时受边界效应影响较大，因此采用1.4节中介绍的方法对SSA的边界效应影响进行削弱.首先需要确定SSA外推数据时合适的窗口长度M和重构阶次P，经实验，当窗口长度M=5760(时间长度4 d)时能够较好地将重力固体潮提取出来，此时重力固体潮只包含在前13阶的重构成分中，因此采用窗口长度M=5760、重构阶次P=13进行SSA数据外推.此外，外推数据个数的选择也是一个值得讨论的问题，外推数据个数过少会使得边界效应的削弱效果不明显，过多则会增加不必要的运算.为了确定最佳的外推数据个数，从多余观测的2600个观测值中分别取不同个数的观测值作为外推值进行边界效应的削弱，将提取得到的重力固体潮值与理论值的残差RMS和STD作为判断最佳外推数据个数的标准，当外推数据个数为1200(原始数据量的8.33%)时，残差RMS和STD均不再有明显的减小，因此确定1200为最佳的数据外推个数.按照以上的参数进行SSA数据外推，获得序列长度为15600的相对重力时间序列.
将外推得到的时间序列按照1.1节中的方法进行SSA分解，提取并截取对应历元下的重力固体潮成分，提取值与理论值的残差如图 12中实线(图中Extrapolated Data by SSA)所示，统计结果见表 4.从图 12中可以看出，用这两种方法获取的重力固体潮与理论值的残差在边界处显示出一定的差异，而在中间段则差异较小，这说明边界效应现象在SSA的解算过程中是明显存在的.与直接利用SSA从原始序列中提取重力固体潮相比，利用SSA外推数据削弱边界效应后的残差最大值有所减小，最小值有所增大，说明残差序列的震荡幅度有所减小.从其他统计结果来看，标准差和均方根值更小，表明这种削弱边界效应的方法是有效的，可以提高重力固体潮提取的效果.
由图 12发现残差时间序列依然具有一定的周期性，这是因为在采用TSoft计算理论重力固体潮时只采用了潮汐表中主要的16个潮波，所以会有其他的重力固体潮潮波反应在残差序列中.此外，本次相对重力观测实验是将CG-5相对重力仪置于山东科技大学实验楼内某一固定点进行的静态相对观测，实验期间由于人们上下班造成的周期震动不能被CG-5相对重力仪全部滤除，剩余的震动会使得重力观测数据中含有周期性的噪声并反映到残差时间序列中.
3.4 SSA拟合序列与原序列残差分析相对重力观测数据中包含的主要信号为零漂和重力固体潮，将SSA提取的零漂与重力固体潮组合，可视为拟合的相对重力数据.除了零漂和重力固体潮信号，相对重力观测数据中还包含了大气变化、水文变化等对重力产生的影响，这些信号也反映在了残差序列中.图 13是原数据与SSA拟合的相对重力数据之间的残差，在右侧边界处有明显的向下漂移的现象，造成这种现象的原因有两个：(1)SSA提取的零漂是观测历元内的相对重力的总体线性趋势，而在原始观测数据中，第10天的零漂突然减小，而这就导致了SSA拟合序列与原序列之间的残差序列在右边界处存在向下漂移的现象；(2)本文虽然对SSA的边界效应进行了削弱，但是边界效应还是存在的，序列边界处受到剩余的边界效应影响.二者残差均值非常小，为-0.1 μGal，说明残差序列在0值附近分布均匀，STD与RMS均为1.1 μGal，可见由SSA提取的零漂与重力固体潮组合信号与原相对重力信号拟合程度较高.
|
图 13 SSA拟合序列与原序列残差 Fig. 13 Residuals between the observed time series and the fitted time series obtained by SSA |
基于相对重力观测数据中包含着重力固体潮、海潮、零漂等信号及噪声的特点，本文提出利用SSA将相对重力时间序列中的重力固体潮信号进行识别和提取.利用SSA从模拟数据中提取的重力固体潮与理论值的残差RMS小于1 μGal，与小波方法相比，采用SSA方法从相对重力信号中提取的重力固体潮更为可靠.CG-5实验结果表明，SSA提取的重力固体潮与采用潮波表计算的理论值残差RMS小于3 μGal，提取的零点漂移与最小二乘线性拟合求得的零漂值之间的差值不到0.4 μGal/d.利用SSA外推数据处理边界效应后残差序列的标准差、均方根值和振幅均有所减小，说明这种方法可在一定程度上抑制边界效应的影响.综上所述，利用SSA方法提取重力观测数据中的重力固体潮是可行、可靠的.但在本实验中，这种方法在提取相对重力数据中的重力固体潮时也存在一些不足：①窗口长度M和相关系数阈值ρijw的选择具有经验性，不同的参数选择所获得结果不同，但总体差异不会过大；②SSA方法无法将具有相似频率的潮波信号进行分离；③SSA方法不能很好地提取相对重力信号中的长周期信号，在进行长周期信号提取时，选择不同的窗口长度提取出的长周期信号存在较大差异，而日周期和半日周期信号则不会因窗口长度的变化而产生较大变化，所以在本实验中未能对长周期的重力固体潮信号进行精确提取.然而，本方法主要针对于弹簧相对重力仪观测数据中的重力固体潮提取，相对于超导重力仪观测，弹簧相对重力仪的观测长度一般较短，且精度较低，所以长周期的重力固体潮潮波对于相对重力仪观测数据的重力固体潮提取和改正无较大影响.
致谢 感谢匿名审稿专家的意见和建议.感谢Scintrex公司的R. Lachapelle先生提供的关于CG-5重力仪的技术支持.感谢瑞典Chalmers理工大学H.G.Scherneck先生对于计算海潮模型潮波参数方面的帮助.
Agnew D C. 1997. NLOADF:A program for computing ocean-tide loading. Journal of Geophysical Research, 102(B3): 5109-5110. DOI:10.1029/96jb03458 |
Beckers J M, Rixen M. 2003. EOF calculations and data filling from incomplete oceanographic datasets. Journal of Atmospheric and Oceanic Technology, 20(12): 1839-1856. DOI:10.1175/15200426(2003)020〈1839:ecadff〉2.0.co;2 |
Chen Q, Van Dam V, Sneeuw N, et al. 2013. Singular spectrum analysis for modeling seasonal signals from GPS time series. Journal of Geodynamics, 72: 25-35. DOI:10.1016/j.jog.2013.05.005 |
Crossley D, Hinderer J, Casula G, et al. 1999. Network of superconducting gravimeters benefits a number of disciplines. Eos, Transactions American Geophysical Union, 80(11): 121-126. DOI:10.1029/99eo00079 |
Ducarme B, Sun H P. 2001. Tidal gravity results from GGP network in connection with tidal loading and Earth response. Journal of the Geodetic Society of Japan, 47(1): 308-315. |
Elsner J B. 2002. Analysis of time series structure:SSA and related techniques. Journal of the American Statistical Association, 97(460): 1207-1208. DOI:10.1198/jasa.2002.s239 |
Farrell W E. 1972. Deformation of the Earth by surface loads. Reviews of Geophysics and Space Physics, 10(3): 761-797. DOI:10.1029/rg010i003p00761 |
Gokhan G, Jin S G. 2016. Evaluation of ocean tide loading effects on GPS-estimated precipitable water vapour in Turkey. Geodesy and Geodynamics, 7(1): 32-38. DOI:10.1016/j.geog.2015.12.008 |
Golyandina N, Zhigljavsky A. 2013. Singular Spectrum Analysis for Time Series. Berlin: Springer. DOI:10.1007/978-3-642-34913-3
|
Guan Z L, Ning J S. 1981. Earth Shape and External Gravity Field (in Chinese). Beijing: Surveying and Mapping Press.
|
Hassani H. 2007. Singular spectrum analysis:methodology and comparison. Journal of Data Science, 5(2): 239-257. |
Hsu H T, Chen Z B, Yang H B. 1982. The effects of ocean tides on gravitational tidal observations. Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese), 25(2): 120-129. |
Hu X G, Liu L T, Ducarme B, et al. 2006. Wavelet filter analysis of local atmospheric pressure effects in the long-period tidal bands. Physics of the Earth and Planetary Interiors, 159(1-2): 59-70. DOI:10.1016/j.pepi.2006.06.001 |
Hu X G, Liu L T, Ducarme B, et al. 2007. Estimation of the pole tide gravimetric factor at the Chandler period through wavelet filtering. Geophysical Journal International, 169(3): 821-829. DOI:10.1111/j.1365-246X.2007.03330.x |
Kondrashov D, Ghil M. 2006. Spatio-temporal filling of missing points in geophysical data sets. Nonlinear Processes in Geophysics, 13(2): 151-159. DOI:10.5194/npg-13-151-2006 |
Kondrashov D, Shprits Y, Ghil M. 2010. Gap filling of solar wind data by singular spectrum analysis. Geophysical Research Letters, 37(15): L15101. DOI:10.1029/2010GL044138 |
Kondrashov D, Berloff P. 2015. Stochastic modeling of decadal variability in ocean gyres. Geophysical Research Letters, 42(5): 1543-1553. DOI:10.1002/2014gl062871 |
Lefèvre F, Lyard F H, Le Provost C. 2000. FES98:A new global tide finite element solution independent of altimetry. Geophysical Research Letters, 27(17): 2717-2720. DOI:10.1029/1999gl011315 |
Lefèvre F, Lyard F H, Le Provost C, et al. 2002. FES99:A global tide finite element solution assimilating tide gauge and altimetric information. Journal of Atmospheric and Oceanic Technology, 19(9): 1345-1356. DOI:10.1175/1520-0426(2002)019〈1345:fagtfe〉2.0.co;2 |
Lei M F, Wang Q J, Liu X L, et al. 2017. Influence of ocean tidal loading on InSAR offshore areas deformation monitoring. Geodesy and Geodynamics, 8(1): 70-76. DOI:10.1016/j.geog.2016.09.004 |
Liu L T, Xu H Z. 2004. Wavelets in airborne gravimetry. Chinese Journal of Geophysics (in Chinese), 47(3): 490-494. DOI:10.3321/j.issn:0001-5733.2004.03.019 |
Longman I M. 1959. Formulas for computing the tidal accelerations due to the moon and the sun. Journal of Geophysical Research, 64(12): 2351-2355. DOI:10.1029/jz064i012p02351 |
Lu Z L. 1996. Theory and Method of Earth Gravity Field (in Chinese). Beijing: PLA Press.
|
Ma X L, Xiao Y, Liu L. 2010. Analysis on residual errors by CG-5 gravimeter's static observation. Resources Environment & Engineering (in Chinese), 24(6): 699-700. DOI:10.3969/j.issn.1671-1211.2010.06.015 |
Matsumoto K, Ooe M, Sato T, et al. 1995. Ocean tide model obtained from TOPEX/Poseidon altimetry data. Journal of Geophysical Research, 100(C12): 25319-25330. DOI:10.1029/95jc02777 |
Matsumoto K, Takanezawa T, Ooe M, et al. 2000. Ocean tide models developed by assimilating TOPEX/Poseidon altimeter data into hydrodynamical model:A global model and a regional model around Japan. Journal of Oceanography, 56(5): 567-581. DOI:10.1023/A:1011157212596 |
Penna N T, Bos M S, Baker T F, et al. 2008. Assessing the accuracy of predicted ocean tide loading displacement values. Journal of Geodesy, 82(12): 893-907. DOI:10.1007/s00190-008-0220-2 |
Schoellhamer D H. 2001. Singular spectrum analysis for time series with missing data. Geophysical Research Letters, 28(16): 3187-3190. DOI:10.1029/2000gl012698 |
Scintrex Limited. 2012. CG-5 Scintrex autograv system operation manual V8.0. Concord: Scintrex Limited.
|
Shen Y, Guo J Y, Liu X, et al. 2017. One hybrid model combining singular spectrum analysis and LS+ARMA for polar motion prediction. Advances in Space Research, 59(2): 513-523. DOI:10.1016/j.asr.2016.10.023 |
Shen Y, Guo J Y, Liu X, et al. 2018. Long-term prediction of polar motion using a combined SSA and ARMA model. Journal of Geodesy, 92(3): 333-343. DOI:10.1007/s00190-017-1065-3 |
Stephane M. 1998. Wavelet Tour of Signal Processing. San Diego: Academic Press.
|
Sun H P, Takemoto S, Hsu H T, et al. 2001. Precise tidal gravity recorded with superconducting gravimeters at stations Wuhan (China) and Kyoto (Japan). Journal of Geodesy, 74(10): 720-729. DOI:10.1007/s001900000139 |
Sun H P, Hsu H T, Jentzsch G, et al. 2002. Tidal gravity observations obtained with a superconducting gravimeter at Wuhan/China and its application to geodynamics. Journal of Geodynamics, 33(1-2): 187-198. DOI:10.1016/s0264-3707(01)00063-1 |
Sun H P, Xu J Q, Ducarme B. 2003. Experimental earth tidal models in considering nearly diurnal free wobble of the Earth's liquid core. Chinese Science Bulletin, 48(9): 935-940. DOI:10.1007/bf03325679 |
Sun H P, Xu J Q, Ducarme B. 2004. Detection of the translational oscillations of the Earth's solid inner core based on the international superconducting gravimeter observations. Chinese Science Bulletin, 49(11): 1165-1176. DOI:10.1360/03wd0242 |
Sun H P, Hsu H, Chen W, et al. 2006. Study of Earth's gravity tide and oceanic loading characteristics in Hong Kong area. Chinese Journal of Geophysics (in Chinese), 49(3): 724-733. DOI:10.1002/cjg2.v49.3 |
Taguchi E, Stammer D, Zahel W. 2014. Inferring deep ocean tidal energy dissipation from the global high-resolution data-assimilative HAMTIDE model. Journal of Geophysical Research Oceans, 119(7): 4573-4592. DOI:10.1002/2013jc009766 |
Torge W. 1989. Gravimetry. Berlin, New York: Walter de Gruyter.
|
VanCamp V, Vauterin P. 2005. Tsoft:graphical and interactive software for the analysis of time series and Earth tides. Computers & Geosciences, 31(5): 631-640. DOI:10.1016/j.cageo.2004.11.015 |
Vautard R, Yiou P, Ghil M. 1992. Singular-spectrum analysis:a toolkit for short, noisy chaotic signals. Physica D:Nonlinear Phenomena, 58(1-4): 95-126. DOI:10.1016/0167-2789(92)90103-t |
Wu H B, Wu L. 2005. Methods for Diagnosing and Forecasting Climate Variability (in Chinese). 2nd ed. Beijing: China Meteorological Press.
|
Wyatt M G, Kravtsov S, Tsonis A A. 2012. Atlantic multidecadal oscillation and northern hemisphere's climate variability. Climate Dynamics, 38(5): 929-949. DOI:10.1007/s00382-011-1071-8 |
Yang W C, Shi Z Q, Hou Z Z, et al. 2001. Discrete wavelet transform for multiple decomposition of gravity anomalies. Chinese Journal of Geophysics (in Chinese), 44(4): 534-541. |
Yu H J, Guo J Y, Li J L, et al. 2015. Zero drift and solid Earth tide extracted from relative gravimetric data with principal component analysis. Geodesy and Geodynamics, 6(2): 143-150. DOI:10.1016/j.geog.2015.01.006 |
Yu H J, Guo J Y, Liu Y, et al. 2017. Precision analysis on CG-5 relative gravimeter through field experiments. Science of Surveying and Mapping (in Chinese), 42(3): 155-160. |
Yu H, Guo J, Kong Q, et al. 2018. Gravity tides extracted from relative gravimeter data by combining empirical mode decomposition and independent component analysis. Pure and Applied Geophysics, 175(5): 1683-1697. DOI:10.1007/s00024-018-1864-3 |
Zabalza J, Ren J C, Wang Z, et al. 2014. Singular spectrum analysis for effective feature extraction in hyperspectral imaging. IEEE Geoscience and Remote Sensing Letters, 11(11): 1886-1890. DOI:10.1109/lgrs.2014.2312754 |
管泽霖, 宁津生. 1981. 地球形状及外部重力场. 北京: 测绘出版社.
|
柳林涛, 许厚泽. 2004. 航空重力测量数据的小波滤波处理. 地球物理学报, 47(3): 490-494. DOI:10.3321/j.issn:0001-5733.2004.03.019 |
陆仲连. 1996. 地球重力场理论与方法. 北京: 解放军出版社.
|
马玄龙, 肖毅, 刘磊. 2010. CG-5重力仪静态观测残差分析. 资源环境与工程, 24(6): 699-700. DOI:10.3969/j.issn.1671-1211.2010.06.015 |
孙和平, 许厚泽, 陈武, 等. 2006. 香港地区重力固体潮和海潮负荷特征研究. 地球物理学报, 49(3): 724-733. DOI:10.3321/j.issn:0001-5733.2006.03.016 |
吴洪宝, 吴蕾. 2005. 气候变率诊断和预测方法. 2版. 北京: 气象出版社.
|
许厚泽, 陈振邦, 杨怀冰. 1982. 海洋潮汐对重力潮汐观测的影响. 地球物理学报, 25(2): 120-129. DOI:10.3321/j.issn:0001-5733.1982.02.003 |
杨文采, 施志群, 侯遵泽, 等. 2001. 离散小波变换与重力异常多重分解. 地球物理学报, 44(4): 534-541. DOI:10.3321/j.issn:0001-5733.2001.04.012 |
于红娟, 郭金运, 刘扬, 等. 2017. CG-5相对重力仪野外实验精度分析. 测绘科学, 42(3): 155-160. |