基于滑动奇异谱分析方法的重力固体潮提取

王建波 王晋南 孙 洁 齐鑫敏 刘相威 高文宗

1 江苏海洋大学海洋技术与测绘学院，江苏省连云港市苍梧路59号，222005 2 菏泽市测绘研究院，山东省菏泽市珠江路1717号，274002 3 山东科技大学测绘与空间信息学院，青岛市前湾港路579号，266590

虽然CG-5重力仪自带软件可提供固体潮改正值，但理论模型是建立在理想条件和固定参数的基础上的，其计算得到的重力固体潮理论值与实际值之间存在一定差异。如何快速从相对重力观测数据中分离重力固体潮改正，对相对重力观测具有重要意义。

SSA方法可以对原先一维时间序列时延排序得到的时滞矩阵进行分解和重构，提取出原序列中代表不同成分的信号，如周期信号、趋势信号、噪声信号等。郭金运等[1]采用该方法提取重力固体潮，并进行周期为10 d的静态相对重力固体潮提取实验，取得很好的效果。虽然相对重力仪静态零漂线性度很好，但零漂总体呈下降趋势[2]。另外，相对重力仪在长时间静态重力观测时会受到观测时间长、零漂线性度改变、非线性误差[3]以及复杂环境噪声的影响，这些影响因素会降低SSA方法提取重力固体潮的精度。为此，本文提出滑动SSA方法，通过设置奇异谱分析数据长度时间窗口，解决长时间静态相对重力观测时重力固体潮的提取问题。

1 原理与方法

1.1 SSA方法、重力固体潮粗差识别及修复

SSA方法的原理、粗差识别及迭代插值方法可参见文献[1]。用SSA方法提取重力固体潮时，一般先利用ω-correlation方法[4]对重建成分(RCs)进行频谱分析和重构。而本文理论重力固体潮值选取与之不同，分为2种：1)对实测重力数据采用Eterna调和分析获得潮波参数，然后用TSoft软件[5]计算得到理论重力固体潮值；2)用CG-5重力仪自带软件进行计算。

粗差识别时，利用SSA方法获得原静态相对重力观测时间序列的主成分，计算原静态相对重力观测时间序列与主成分之间的残差，分析残差值时间序列，并以3σ阈值原则识别粗差。具体过程为：利用SSA方法处理原静态相对重力观测时间序列，利用ω-correlation方法分析各重建成分之间的相关性，找出趋势项、周期项以及第i阶次，当重建成分阶次大于i时，各个重建成分之间不能很好地分离；合并前i阶项得到主成分，计算原相对静态重力观测值序列和主成分序列之间的残差并识别粗差。

本文通过奇异谱迭代方法插补粗差处的重力值，与文献[1]不同之处在于：将残差时间序列中粗差位置处的残差值设置为0，然后通过SSA方法分解重构重建残差趋势项，用趋势项的值代替原缺失位置的0值。循环此过程，直至缺失位置处前后两次插补的数据残差RMS小于10-3μGal。这时，将重建的残差数据序列与之前提取的主成分数据序列相叠加，得到修正后的静态相对重力时间观测序列。

1.2 零漂和固体潮提取思路

针对长时间静态相对重力观测数据中零漂和固体潮的变化特点和SSA方法分析有限长度周期性信号的特点，本文提出采用滑动SSA方法提取长时间重力固体潮改正：

1)将原M天的静态相对重力观测值时间序列按L天为一个时间段进行滑动分组，滑动窗口长度为1 d，即任意相邻两组数据序列的起始时间相隔1 d，得到每个时间段M-L+1组数据序列。

2)对每个时间段内每组静态相对重力观测值序列，利用多项式拟合得到L天内的零漂值。从每组观测值中减去对应多项式拟合计算得到的零漂值，对其结果采用SSA方法提取固体潮，然后提取每组数据序列的固体潮改正。

3)对每个时间段内M-L+1组数据均通过多项式拟合提取零漂，将剔除零漂后的重力数据观测序列通过SSA方法提取重力固体潮改正值。重力固体潮(或零漂)按同一天取平均值获得当天重力固体潮改正(或零漂)的最终结果，即第1～M天分别在第1、2、…M-L+1组中找到对应值，并按1 d、2 d、…、M天取平均。为降低SSA中边界效应的影响，在精度对比统计中不统计每组前2 d和最后2 d的数据。

2 实验与分析

2.1 实验数据及数据预处理

相对静态重力观测时间序列中除包含重力固体潮外，还包括海潮负荷、零漂和噪声等。本次实验利用CG-5相对重力仪(序列号41221)在固定点(坐标120.10°E、36.00°N；大地高25.5 m)进行连续2个月(2016-03-05 00：00～05-04 23：59)的静态相对重力观测。相对重力仪采样频率为6 Hz，仪器前55 s进行数据采集，后5 s进行数据处理，因此1 min内可得到330个数据。相对重力仪自带地震滤波功能，可对获取的静态相对重力观测数据进行低频噪声过滤，舍弃高于6倍标准偏差的高频噪声。滤波之后对样本数据每分钟进行一次平均处理并且作为结果自动记录，61 d共获得87 840个数据。对相对重力观测数据进行预处理，采用本文提出的方法进行重力固体潮改正提取。作为对比，采用多项式拟合提取零漂，将剔除零漂后的重力数据采用SSA提取重力固体潮改正，技术路线如图1所示。

图1 滑动SSA提取重力固体潮技术路线

首先利用SSA方法进行粗差探测，将4 d共5 760个重力观测数据作为嵌入维数，对87 840个相对重力观测数据时间序列进行SSA分解和重构，获得RCs；通过ω-correlation方法分析各RC之间的相关系数，前30阶重构成分ω-correlation分析结果如图2所示。根据相邻两个RC之间的相关系数大于0.85可以判断，R2、R3、R4、R5、R7、R8、R9、R10、R16、R17、R20、R21、R26、R27、R28、R29均具有周期性。从重力数据序列SSA得到的前30阶重构部分FFT功率谱分析结果(图3)可以看出，从第15阶开始，不同频率混叠较严重，而且振幅非常小，各重构成分之间无法很好地分离，因此合并前14项重构成分作为原重力数据时间序列主成分。原静态相对重力观测与主成分之间残差时间序列如图4(a)所示，根据3σ准则选择阈值区间为(-0.009 8 mGal，0.009 8 mGal)，若超出此阈值区间则判定为粗差，共624个。根据奇异谱迭代方法修正残差处的粗差相对重力值，经过迭代后得到粗差处相邻两次迭代修复值残差的RMS为9.57×10-6μGal，修复后的残差值如图4(b)所示。可以看出，残差序列较为平稳，插补效果良好。修复后的残差叠加主成分可得到修复后的静态相对重力观测时间序列，结果如图5(a)所示。海潮和重力固体潮成因相同，因此具有相同的频率。本文基于NAO.99b海潮模型[6]，结合TSoft软件计算本次实验同历元下的重力海潮负荷值时间序列，结果如图5(b)所示，将海潮对相对重力的影响从原时间序列中去除，结果如图5(c)所示。

图2 静态相对重力数据序列SSA前30阶重建部分的ω-correlations

图3 静态相对重力数据序列SSA前30阶重构部分的FFT功率波谱分析

图4 原静态相对重力观测时间序列与主成分之间的残差

图5 静态相对重力观测时间序列

2.2 滑动SSA方法提取重力固体潮改正

为了验证滑动SSA方法提取重力固体潮改正的效果，对61 d静态相对重力观测时间序列分别以10 d、15 d、20 d、25 d、30 d、40 d、50 d为时间段，并以1 d为滑动窗口进行分组，得到每个时间段内包含的数据序列组数分别为52、47、42、37、32、22和12。

在提取重力固体潮之前，首先需要剔除每组静态相对重力观测数据中的零漂项。采用SSA方法识别CG-5相对重力观测数据中的趋势，将SSA窗口长度设为4 d，对每组数据进行分解和重构，得到趋势项RC1。然后采用多项式拟合提取零点漂移，拟合公式为：

y=ax2+bx+c

(1)

式中，a、b、c为拟合参数。利用SSA提取零漂时可剔除噪声，因此效果会更好。对每组数据序列剔除拟合后的零漂值，然后选择嵌入维数5 760进行 SSA，对重构成分进行ω-correlation相关性分析。个别分组中两个相邻RC之间的相关系数大于0.80时无数据，因此对比实验中均选择相关系数大于0.65(当相关系数为0.75时结果类似，相关系数为0.65时提取结果精度提高更明显)。从重构项i≤15中得到每组中同频率重构序列，对重构序列进行合并。对每组重构序列进行求和，计算每组提取的重力固体潮改正。SSA在分解重构过程中存在边界效应，因边界效应不是本文主要研究内容，为削弱其影响，将每组前2 d和最后2 d共5 760个重力固体潮提取结果舍去，再采用本文方法计算最终的重力固体潮改正值，得到的结果分别记为HSSA_10、HSSA_15、HSSA_20、HSSA_25、HSSA_30、HSSA_40和HSSA_50。

为了验证提取精度，将实测重力数据采用Eterna调和分析获得潮波参数后用TSoft软件计算得到重力固体潮值(记为ET_Data)和CG-5重力仪自带软件计算的重力固体潮值(记为CG5_Data)分别作为参考值，提取结果与同时刻ET_Data和CGS_Data之间的差值如图6所示，差值统计如表1(单位μGal)所示。

图6 不同方法提取的重力固体潮残差对比

表1 滑动SSA和SSA提取长时间重力固体潮残差统计

2.3 多项式拟合方法提取零漂和SSA提取重力固体潮

将61 d静态相对重力观测数据序列作为一个时间段长度，采用多项式拟合方法提取零漂值，设嵌入维数为5 760，利用SSA分解重构获得RC1趋势项，并利用式(1)进行多项式拟合。将零漂数据从静态相对重力序列中剔除，然后选择嵌入维数5 760进行SSA，对重构成分进行ω-correlation相关性分析，从结果中选出重构项i≤15和相关系数大于0.65的两个相邻RC项，得到同频率重构序列，并对重构序列进行合并。对重构序列进行求和，计算提取的重力固体潮改正，得到的结果记为SSA_Data。为削弱边界效应影响，同样将前2 d和最后2 d共5 760个重力固体潮提取结果舍去，其与同时刻ET_Data和CG5_Data的结果比较如图6所示，差值统计如表1所示。

2.4 不同方法提取重力固体潮残差分析

由表1可知，无论是滑动SSA还是SSA方法，提取重力潮改正结果与理论值差值的标准差(STD)均小于11 μGal，说明SSA方法提取重力固体潮改正结果较可靠，精度较高。对于同一种提取重力固体潮改正方法，选用CG5_Data作为理论参考值，其提取结果比选用ET_Data作为理论参考值时大，说明SSA或滑动SSA方法的提取结果更接近ET_Data理论值。采用滑动SSA提取重力固体潮改正时，每个连续时间段长度不同，其提取结果与理论值残差的STD和RMS不同。选择连续时间段长度为30 d和40 d时，滑动SSA提取重力固体潮改正值与理论值残差的STD和RMS比SSA方法更小，结果更优。选择滑动时间段长度为10 d、15 d、20 d、25 d和50 d时，滑动SSA方法提取重力固体潮改正值与理论值残差的STD和RMS大于SSA方法，表明SSA方法更优。选择滑动时间段长度为30 d时，滑动SSA方法提取重力固体潮改正值与ET_Data理论值残差的STD和RMS比SSA方法均小0.43 μGal；滑动SSA方法提取值与CG5_Data理论值残差的STD和RMS比SSA方法分别小0.201 μGal和0.219 μGal。选择滑动时间段长度为40 d时，滑动SSA方法提取重力固体潮改正值与ET_Data理论值残差的STD和RMS比SSA方法均小0.665 μGal；滑动SSA方法提取值与CG5_Data理论值残差的STD和RMS比SSA方法分别小0.614 μGal和0.602 μGal。选择滑动时间段长度为10 d时，采用滑动SSA方法提取重力固体潮时，提取结果与理论值差值的STD和RMS最大，即提取精度最低。残差最大值与最小值之间的差值可表明残差范围，范围越小说明残差的离散程度越小。当选用CG5_Data作为理论参考值时，除选择10 d进行滑动计算的结果外，滑动SSA方法提取结果的离散程度均小于SSA方法，提取结果残差离散程度最小的方法为选择20 d进行滑动计算；当选用ET_Data作为理论参考值时，除选择10 d、50 d滑动计算外，滑动SSA方法提取结果的离散程度均小于SSA方法，提取结果残差离散程度最小的方法为选择20 d进行滑动计算。结果表明，通过滑动SSA方法提取固体潮可降低提取结果残差的离散程度。

由图6可知，残差序列均在零值附近摆动，说明提取的零漂与参考信号拟合程度较高。残差时间序列仍然具有明显的周期性，因为相对重力序列中除零漂和固体潮外，还包括其他因素产生的周期信号，如天气变化、水文变化以及其他干扰噪声，这些因素在残差序列中会有所反映。由于ET_Data理论值是实测重力数据采用Eterna调和分析获得潮波参数，再利用TSoft软件计算得到的，以ET_Data为参考值时比以CG5_Data为参考值时残差振幅略小。选用不同连续时间段进行滑动计算的结果区别较大，其中选用连续时间段30 d、40 d计算结果较好。

3 结 语

本文将时间长度为2个月的静态相对重力观测数据序列按时间窗口10 d、15 d、20 d、25 d、30 d、40 d和50 d分别分为52组、47组、42组、37组、32组、22组和12组，采用滑动SSA方法进行重力固体潮改正提取实验。结果表明，SSA方法提取重力固体潮结果精度较高，采用滑动SSA和SSA方法提取的重力固体潮值与理论值(CG5_Data和ET_Data)残差的标准差均小于11 μGal。与SSA方法相比，选用合适时间段长度的滑动SSA方法提取固体潮可以降低提取结果残差的离散程度。选用CG5_Data值作为参考值、连续30 d和40 d数据作为一组时，滑动SSA方法提取的重力固体潮改正结果残差的RMS比未使用滑动SSA方法时分别小0.219 μGal和0.602 μGal；选用ET _Data值作为参考值、连续30 d和40 d数据作为一组时，滑动SSA方法提取的重力固体潮结果残差的RMS比未使用滑动SSA时分别小0.430 μGal和0.665 μGal。因此，选择合适的连续滑动时间段，滑动SSA方法比SSA方法提取重力固体潮改正结果精度略高。