image
Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 
Online first  /  Current issue  /  All issues  /  Top access  /  RSS

Progress in Geophysics >

2024 , Vol. 39 >Issue 6: 2137 - 2152

DOI: https://doi.org/10.6038/pg2024HH0487

Extraction of interannual signals from superconducting gravimeter data and their characterization

  • XiZhi HU , 1, 2 ,
  • GenYou LIU , 1, *
Expand
  • 1 State Key Laboratory of Geodesy and Earth's Dynamics, Innovation Academy for Precision Measurement Science and Technology, CAS, Wuhan 430077, China
  • 2 University of Chinese Academy of Sciences, Beijing 100049, China

Received date: 2023-12-26

  Online published: 2025-01-14

Copyright

Copyright ©2024 Progress in Geophysics. All rights reserved.
Fold

Abstract

The gravity variation of the Earth is a comprehensive reflection of the global mass and its motion, and the superconducting gravimeter is capable of observing small changes in the gravity field. In order to characterize the cyclic variations of gravity on interannual scales and the physical mechanisms behind them, in this paper, for the observation data of seven superconducting gravimeter stations, after deducting the instrument drift, synthetic tide, and atmospheric loading effects to obtain the gravity residuals, the normal Morlet wavelet transform method is utilized to identify and extract the signals with the periods of 1 a and 1.2 a, respectively, and the gravimetric factors of the pole tide associated with the Chandler wobble are estimated by the least-squares method. Then residuals are deducted from gravity pole tide, obvious interannual signals common to the stations are extracted, which have periods of 2.3 a, 2.7 a, 3.3 a, 3.7 a, 4.8 a and 6 a, and their corresponding mean amplitudes are obtained in the time domain. The sources of these periodic signals are preliminarily analyzed, in which the variations on the 2 a to 5 a scales have a high correlation with hydrological loads, and the 6 a signals in the residuals are highly correlated with the sequence of gravity variations computed from the GNSS-observed model of the vertical displacements of the Earth's surface 6 a.

Key words： Time-varying gravity field; Superconducting gravimeter; Normal Morlet wavelet transform; Gravity pole tide; Interannual signal

Cite this article

XiZhi HU , GenYou LIU . Extraction of interannual signals from superconducting gravimeter data and their characterization[J]. Progress in Geophysics, 2024 , 39(6) : 2137 -2152 . DOI: 10.6038/pg2024HH0487

0 引言

重力时变信号反映了地球内部物质的分布及运动规律，对重力时变特征的分析有助于分析地球各圈层之间相互作用.重力、极移与日长、地表位移或形变等不同观测数据的年际变化及其背后可能蕴含了地球物理激发源的地球内部的动力学过程(如Mound and Buffett, 2006; Gillet et al., 2010; Holme and de Viron, 2013; Chao，2017；Ding et al., 2021).地球内部质量分布变化导致的重力变化可能发生在与之相同的时间尺度上，并存在于地球各圈层导致的重力变化中.由于观测技术多在地表和地球外部的原因，可能来自地球内部的长周期重力信号一般较微弱，然而超导重力仪观测具有高精度(最高能够达到10-10 m/s2)、低漂移率、高的稳定性和灵敏度特性(孙和平等，2021)，使得其观测数据成为我们目前提取长周期微小信号、认识地球各圈层内部活动及圈层间耦合等现象的重要资料.
从超导重力仪观测数据中提取信号的研究多集中于从秒到数月的频段.徐华君等(2009)利用标准Morlet小波变换(Normal Morlet Wavelet Transform，NMWT)从超导重力仪观测序列中提取出振幅为1 nm/s2的径向球型振荡(0S0) 信号(周期约为20.47 min)，将Q值计算由传统的频域转到时域，简化了Q值计算过程，同时也表明NMWT是一个分析自由振荡的理想工具；Li等(2020)采用自相关方法，联合STS-1地震仪，从超导重力仪长期观测资料中提取出周期介于133.3~500 s的长周期面波和地球自由振荡；徐建桥等(2005)利用超导重力仪观测序列识别了周期在10~17 h的4个核模信号；孙和平等(2003)基于全球超导重力仪观测数据测定出地球液核近周日共振参数.随着超导重力仪观测精度的提高以及资料处理方法的改进，考虑了一些更精细的误差与改正，比如Hu等(2005)基于Daubechies高消失矩小波构建带通滤波器，得到4~40天周期的时-频气压导纳值，提高了超导重力仪大气压在相应周期的改正精度.
由于目前全球超导重力仪持续观测数据的时间跨度相比极移、日长变化等数据较短(最长约25年)，因此针对其年际变化的时频分析相对较少，且多集中于重力极潮尤其重力极潮因子估计的研究.极移的周年项和钱德勒摆动的能量谱峰在超导重力仪观测序列的功率谱中被检测到，并根据所用台站结果的叠加准确地计算出钱德勒周期的极潮潮汐因子(杨学峰等, 2004；Xu et al., 2004)；Ducarme等(2006)利用POLAR程序对全球9个扣除固体潮和海潮负荷以及大气负荷的超导重力数据进行回归分析，并顾及平衡海洋极潮效应获得了极潮潮汐因子的精确估值；Hu等(2007)对超导重力残差利用Daubichies小波滤波方法得到其256~512天频带的数据，并对其进行最小二乘估计得到钱德勒周期的重力极潮潮汐因子；Wang等(2015)利用集合经验模态分解方法(Ensemble Empirical Mode Decomposition，EEMD)从长度为7年的超导重力仪观测资料中提取了包含重力极潮在内的长周期信号；Ziegler等(2016)使用堆叠法算得重力极潮潮汐因子，并用加权组合的方法显著减少了局部水文负荷的影响；Ding和Chao(2017)对GPS三维位移场和超导重力观测序列利用最优序列估计(Optimal Sequence Estimation, OSE)方法获得了钱德勒周期的复Love数，为地幔弹性和地球内部的衰减模型提供了更精确的约束；孙和平等(2018)基于EEMD从5个超导重力仪台站的观测资料中算得钱德勒周期的重力极潮潮汐因子；申文斌和刘任莉(2009)采用小波分析方法从超导重力仪数据中探测到与内核超速旋转有关的重力变化.黄乘利和金文敬(1999)从地球自转轴方向变化引起的地表点经纬度和自转速度变化对离心力位的影响入手，结合极移与日长观测数据对日长和极移的重力效应分别作了定量估计；对海洋极潮进行自洽平衡响应建模的结果显示，自洽平衡海洋极潮对超导重力仪站重力极潮的负荷效应可达nm/s2量级，对重力极潮振幅因子δ的影响达1%, 而其引起的相位延迟可忽略不计(Chen et al., 2008).
本文利用NMWT对观测连续性较强的7台超导重力仪观测数据的残差进行处理，对提取的年际信号的时域变化特征进行了初步分析，为进一步完善超导重力仪观测数据的各项物理改正模型以及研究地球各圈层质量变迁和运动提供一定的参考.

1 数据来源及预处理说明

从国际地球动力学与固体潮服务组织(International Geodynamics and Earth Tide Service IGETS)管理的43个超导重力仪台站的观测序列(ISDC: IGETS Data Base(gfz-potsdam.de))中选择了观测质量较好、记录持续较长的7个台站，如表 1所示.其中，CB、MB、MC、MO、ST和WE站均采用采样间隔为1 min的Level 3数据，为IGETS的工作人员处理的结果，可直接使用；Level 3数据结果考虑了固体潮、海潮负荷、大气负荷、地球自转变化(包括极移与日长)改正以及趋势项影响，处理方法和流程如图 1所示(以Membach为例).用于时频分析的数据均为降采样至1天的重力残差数据.
表1 所选超导重力数据情况

Table 1 The case of selected superconducting gravity data

台站名称 仪器名称 经度和纬度 高程/m 观测时段
Canberra(CB) GWR C031 35.3206°S, 149.0077°E 762.75 1997-07-01—2021-12-31
Medicina(MC) GWR C023 44.5219°N, 11.6450°E 28.00 1998-01-01—2019-04-30
Membach(MB) GWR C021 50.6085°N, 6.0095°E 250.00 1995-08-05—2021-03-31
Moxa(MO) CD034_L 50.6450°N, 11.6160°E 455.00 2000-01-02—2019-12-31
Strasbourg(ST) GWR C026 48.6217°N, 7.6838°E 180.00 1997-02-27—2018-11-12
Wettzell(WE) GWR CD029_U 49.1440°N, 12.8780°E 613.70 1998-11-05—2010-08-31, 2013-01-01—2018-03-31
GWR CD030_U 2010-09-01—2012-12-31
Wuhan(WU) GWR T004 TT70 C032 30.5159°N, 114.4898°E 80.00 1997-12-21—2012-07-28
OSG 065 30.5158°N, 114.4898°E 89.30 2013-03-28—2021-03-31
图1 Membach站重力残差的导出

Fig 1 Output process of gravity residuals at Membach station

图 1为比利时Membach台站的GWR C021仪器观测的Level 3数据处理过程，观测周期为1995年8月5日至2021年3月31日，图中重力数据的单位均为nm/s2(1 nm/s2=0.1 μGal)，被扣除的部分为红色，剩余部分为蓝色，图 1中各子图表示的分别是：
(a) Membach台站的原始重力观测值；
(b) 随时间变化的固体潮合成潮与海潮负荷.在高频下使用由最小二乘法调整的局部潮汐模型；而在低频下使用DDW99重力潮汐因子对潮汐信号进行建模和固体地球潮汐的HW95引潮位，以及FES2014c的海洋潮汐负荷；
(c) 由原始重力观测值扣除合成潮负荷获得的残差；
(d) 采用ERA5地表气压数据基于反变气压计假设计算的大气负荷效应，且对于与站之间的角度距离小于0.10°的记录，ERA5大气由1 min的本地大气记录代替；
(e) 蓝色实线为扣除固体潮合成潮与海潮负荷和大气负荷效应获得的残差，红色实线为趋势项；
(f) 蓝色实线为扣除固体潮合成潮与海潮负荷、趋势项和大气负荷效应获得的残差，红色实线为理论重力极潮和日长引起的重力变化，这使用了IERS EOPC04每日序列(http://hpiers.obspm.fr/IERS/eop/EOPC04/)并假设潮汐重力因子δ为1.16进行建模，日长变化引起的重力变化通常比重力极潮小两个数量级，在本研究中可忽略(Wahr, 1985; Xu et al., 2004)；此外对海洋极潮进行自洽平衡响应建模；
(g) 扣除固体潮合成潮与海潮负荷、趋势项、大气负荷效应和重力极潮效应获得的残差.该残差用于提取其他年际信号.
WU站提供的Level 3数据较短，因此从采用率为1 h间隔的Level 2数据(未经固体潮、海潮和大气负荷改正)进行了重新处理获得Level3数据，具体处理流程和结果如图 2所示.针对Level 2残差中2012年7月29日至2013年3月27日因仪器原因缺失的数据，进行了三次样条插值来补充.图 2中，(a)为1 h采样率的Level2观测数据; (b)为gotic2软件计算的Wuhan本地合成潮负荷；(c)是扣除固体潮和海潮负荷获得的残差；(d)是来自EOST Loading Service(EOST/ITES Loading Service (u-strasbg.fr))的Wuhan大气ERA5(ECMWF reanalysis)重力负荷模型；(e)中蓝色实线是扣除固体潮、海潮负荷和大气负荷的残差，红色分段直线为拟合的趋势项；(f)中蓝色实线是经趋势项扣除、突跳修正、插值后的1 h采样残差，未扣除理论重力极潮和日长变化的影响；(g)是在(f)蓝色实线的基础上扣除理论重力极潮和日长变化影响得到的残差；WU站的重力残差与其他台站统一为1天采用，用于年际时频分析.
图2 Wuhan站重力残差的导出

Fig 2 Output process of gravity residual at Wuhan station

2 信号提取方法与仿真实验

2.1 方法描述

本文主要采用NMWT、正交Daubechies小波、傅里叶变换以及“极值点延拓”策略(Duan and Huang, 2020; Duan et al., 2017)相结合的信号提取流程来准确识别与提取重力年际变化中的长周期(准)谐波信号.
NMWT由Liu等(2007)提出，该方法具备频率分辨能力高的特性，可以直接无偏地将谐波信号从数据中提取出来而无需进行逆变换，是进行谐波分析的有力工具，在地球物理与大地测量等领域的定量研究有着成功的应用(Su et al., 2014; Cai et al., 2018; Wang et al., 2018;Duan and Huang, 2020).
NMWT方法的定义公式为：
$W_{\mathrm{g}} f(a, b)=\frac{1}{|a|} \int_{-\infty}^{+\infty} f(t) g^*\left(\frac{t-b}{a}\right) \mathrm{d} t, $
其中，a是尺度因子，b是时移因子，a, b∈R(a≠0)，R表示实数域；上角标“*”表示共轭运算符；核函数g(t)称为标准Morlet小波基，其表达式为
$g(t)=\frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{t^2}{2 \sigma^2}+\mathrm{i} 2 \pi t\right)$
其中，σ是窗宽因子，(∀σ∈R)≥1，控制初始的时频分辨率，该因子是一个经验值，经后面的仿真实验验证，当σ=3时，在提取重力长周期信号时结果相对最优.
NMWT具有小波变换共有的边缘效应问题，会导致提取的信号在振幅确定上存在失真现象，如图 3g中的蓝色实线，实际处理中周期越长的频带其振幅失真程度越大.虽然傅里叶变换方法的频率分辨率相对于NMWT较低、无法展现信号振幅随时间的变化，但振幅的确定更直观准确，因此傅里叶变换计算的振幅可以作为NMWT提取信号提供初始参考.
3 极值点对称延拓与直接对称延拓的比较

Comparison of symmetric extreme points extension and direct symmetric extension

实际序列往往频率范围较宽，相近的周期之间会形成相互干扰，单个年际信号的周期通常较为稳定、所在的频率范围较窄.采用正交Daubechies小波可以从原始序列中获取含有目标信号的频段，缩小频率范围、减少周边干扰.将Daubechies小波选取的频段和原始序列的傅里叶频谱进行对比，以两者重叠的频率区域作为的确存在的有效频段，如图 3e中的两傅里叶频谱中2.8~5.6 a频段.
NMWT存在边缘效应，影响区间长度的可由 $ R_{\mathrm{g}}(a)=1.643 \delta|a|$ 计算，其中δ为窗宽因子，a为时间尺度即目标信号的周期(Liu et al., 2007).本文采取极值点延拓策略处理边界效应：利用NMWT方法的相位无偏特性确定观测序列两端的极点位置，保留两极值点之间的序列，然后将保留的序列在两端作对称延拓.
为说明极值点延拓在减少边缘效应方面的优势，假设原始数据是一个采样间隔1 d、长度10000 d、周期为4.8 a的谐波信号，如图 3a中白色矩形区域内的绿色实线(真实值)，对其分别进行了直接延拓和极值点延拓.采用直接双向对称延拓后得到的新序列为图 3a中绿色实线，利用NMWT获得的信号如图 3中所有的蓝色实线，在原始数据的两端出现了周期、相位和振幅的较大偏差；采用极值点对称延拓后进行NMWT获得的信号为红色虚线，在白色矩形区域内与真实值吻合得很好.因此，在实测数据的处理中，我们采用极值点延拓策略去除边缘效应.

2.2 仿真实验与信号提取流程

基于2.1中的方法和策略，我们结合仿真实验来说明信号提取的整体流程.
仿真数据构建了采样间隔为一天、长度为10000天的超导重力残差，观测噪声为±3 nm/s2的白噪声，其中包含了7个谐波信号，周期(振幅) 分别为2.3 a(5 nm/s2)、2.7 a(5 nm/s2)、3.3 a(8 nm/s2)、3.7 a(3 nm/s2)、4.8 a(12 nm/s2)、6 a(5 nm/s2)和8.5 a(4 nm/s2)，如式(3)所示：
$f(t)=\sum_{n=1}^7 f_n(t)+\operatorname{noise}(t), $
其中: t为时间变量，$ f_1(t)=5 \cos \left(\frac{2 \pi}{2.3 \times 365.25} t\right)$、$f_2(t)=5 \cos \left(\frac{2 \pi}{2.7 \times 365.25}(t-180)\right)$、$\begin{aligned} f_3(t)=8 \cos\left(\frac{2 \pi}{3.3 \times 365.25}(t-100)\right)\end{aligned}$、$ f_4(t)=3 \cos\left(\frac{2 \pi}{3.7 \times 365.25}(t+100)\right)$、$\begin{aligned}f_5(t)=12 \cos \left(\frac{2 \pi}{4.8 \times 365.25}(t-100)\right)\end{aligned}$、$\begin{aligned}f_6(t)=5 \cos \left(\frac{2 \pi}{6 \times 365.25}(t+100)\right)\end{aligned}$、$f_7(t)=4 \cos \left(\frac{2 \pi}{8.5 \times 365.25} t\right)$，noise(t)是复合信号(图 4a).
图4 仿真实验中4.8 a信号的识别与提取

Fig 4 Identification and extraction of 4.8 a signal in simulation experiment

当所选频段中仍存在多个信号时，为避免其他信号对目标信号提取造成干扰，采用信号强度由高到低的顺序依次处理，先识别和提取出能量较高的信号，从残差中扣除这一信号获得新的残差，再从新的残差识别和提取能量较低的信号，以此类推.信号提取流程如下：
(1) 明确目标：从原始序列的傅里叶频谱中获悉其中有哪些明显的周期信号(见图 4e蓝色虚线)，注意到周期在4~5 a的谱峰最高，因此首先对其进行提取.
(2) 对目标信号所在频段进行定位：使用正交Daubechies小波滤波方法从原始序列中获取包含目标信号的频段，这里运用正交Daubechies小波(N=14)将复合信号按频率划分为趋势项a1和细节项d1~d14，其中d10为2.8~5.6 a频段的信号，d10与复合信号的傅里叶频谱对比如图 4e所示，可见正交Daubechies小波对复合信号频段的划分准确.
(3) 确定目标信号两端极值点的位置：使用NMWT方法从目标频段中提取相位无偏但振幅受到边缘效应影响的目标信号，以确定极值点的位置，如图 4g中蓝色实线.
(4) 提取目标信号：利用(3)确定的极值点位置，进行极值点延拓获得新序列，使用NMWT方法从新序列中提取目标信号，最后将提取结果截取至原始时间区间，如图 4g中红色虚线.
根据以上步骤获得的信号(图 4g中红色虚线)，可以确定其周期为4.8 a，振幅为12.0 nm/s2，与设计的仿真信号一致.
为评定4.8 a信号结果的精度，我们分别在时、频两域采用每一历元提取结果与仿真值之差的均方根误差(Root Mean Square Error)来衡量，即：
$\left\{\begin{array}{l}\delta_{\mathrm{A}}=\sqrt{\frac{\sum_{t=1}^n\left(A_{\mathrm{obs}, t}-A_{\mathrm{model}, t}\right)^2}{n}} \\\delta_{\mathrm{P}}=\sqrt{\frac{\sum_{t=1}^n\left(P_{\mathrm{obs}, t}-P_{\mathrm{model}, t}\right)^2}{n}}, \end{array}\right.$
其中，δA、δP分别为提取结果的振幅和周期误差，n为历元总数，t为历元数，下角标obs, t、model, t分别对应NMWT提取结果和仿真信号，A、P分别对应振幅与周期.
根据式(4)计算出4.8 a信号结果的周期与振幅精度分别为0.07 a、0.5 nm/s2，可表示为：4.8±0.07 a、12.0±0.5 nm/s2；其他信号的提取结果及精度分别为：2.3±0.04 a、5.0±0.2 nm/s2；2.7± 0.05 a、5.0±0.2 nm/s2；3.3±0.05 a、8.0±0.2 nm/s2；3.7±0.05 a、3.0±0.2 nm/s2；6.0±0.06 a、5.0±0.3 nm/s2；8.5±0.18 a、4.0±0.4 nm/s2，提取结果与仿真数据吻合得很好.

3 实测数据的处理结果与分析

3.1 结果概述

各台站观测在扣除固体潮、海潮负荷、大气负荷以及趋势项的基础上，扣除理论重力极潮前、后的重力残差(分别对应图 1f、g)的傅里叶频谱如图 5所示.为便于区别和描述，我们将扣除重力极潮前的残差称为残差1，扣除重力极潮后的残差称为残差2.
图5 各台站扣除理论重力极潮前、后的重力残差的傅里叶频谱

Fig 5 Fourier spectrum of gravity residuals before and after subtracting the theoretical gravity pole tide at each station

从各台站残差的傅里叶频谱中可以看出，各台站的残差1都存在较为明显的周期约为1 a和1.2 a的谱峰，分别对应重力极潮的周年项和未能扣除的周年影响、重力极潮的钱德勒项；残差2存在较多其他谱峰，这说明残差中蕴含着丰富的年际信号，且各台站具有区域差异.因所用数据的长度有限，傅里叶频谱在长周期频段的频率分辨率低，显示的低频信号周期可能不够准确，另外傅里叶频谱无法展示各周期信号的时变情况，所以仅作为信号提取的参考.
在初步确定各台站重力残差1、2中可能存在的年际信号后，使用标准Morlet小波变换和“极值点延拓策略”对各台站重力残差1进行处理得到1 a和1.2 a信号，对各台站重力残差2进行处理得到2.3 a、2.7 a、3.3 a、3.7 a、4.8 a和6 a信号.各周期信号的平均振幅如表 2所示，时域提取结果(左列)和对应的NMWT时频振幅谱(右列)如图 6所示.
表2 各台站重力残差中年际周期信号的周期和平均振幅

Table 2 Period and average amplitude of the interannual cycle signal in the gravity residuals at each station

CB MB MC MO ST WE WU
1 19.7 20.1 8.2 18.2 28.3 18.2 17.3
1.2 14.0 15.3 17.9 12.3 15.2 14.9 11.0
2.3 2.19 7.20 2.57 1.53 4.37 3.82 7.71
2.7 3.3 4.4 2.3 2.5 7.8 8.0 4.3
3.3 1.4 1.1 3.9 4.3 2.6 2.3 11.3
3.7 5.3 4.2 2.7 2.6 2.5 2.3 3.2
4.8 3.2 9.0 8.3 12.9 4.9 5.8 2.1
6 4.6 3.7 3.9 4.3 5.4 4.3 5.2
图6 各站时域提取的年际信号的周期、振幅和NMWT时频系数谱

(a)1 a与1.2 a信号；(b)2.3 a信号；(c)2.7 a信号；(d)3.3 a信号；(e)3.7 a信号；(f)4.8 a信号；(g)6 a信号.图(a)左列子图每两行为同一台站的1 a和1.2 a信号，其中奇数行为1.2 a信号，偶数行为1 a信号.

Fig 6 Period, amplitude and NMWT time-frequency coefficient spectra of the interannual signals extracted in the time domain from each station

(a)1 a and 1.2 a signals; (b)2.3 a signal; (c)2.7 a signal; (d)3.3 a signal; (e)3.7 a signal; (f)4.8 a signal; (g)6 a signal. In the subgraph of the left column of (a), every two rows are 1 a and 1.2 a signals of the same station, where the odd rows are 1.2 a signals and the even rows are 1 a signals.

3.2 对结果的分析

3.2.1 1 a与1.2 a信号

对于残差1中的1 a信号和1.2 a信号，以Membach站为例，将1 a和1.2 a信号在时域求和，得到的结果与考虑地球弹性响应(Wahr, 1985; 黄乘利和金文敬, 1999)与自洽海洋平衡极潮(Chen et al., 2008)的理论重力极潮理论值符合得很好(图 7).相比杨学峰等(2004)利用自回归模型估计残差的功率谱密度和积谱密度的频域结果，NMWT和极值点对称延拓策略提取的结果不仅在频域上与理论重力极潮吻合，在时域上也得到了较为准确的极移重力变化结果，表明了NMWT方法在准谐波探测上的优势.
图7 利用NMWT方法从超导重力数据中探测极移导致的重力变化

Fig 7 Detection of gravity changes due to polar motion from residuals of superconducting gravity data using the NMWT method

极潮潮汐因子的确定可为了解地球形变、地球内部构造和物理参数提供约束，也可以为台站进行重力极潮改正提供参考.为估计重力极潮潮汐因子(包括振幅因子δ和相位滞后κ)，则需利用最小二乘匹配法对1 a和1.2 a信号的提取结果进行估算.考虑到重力极潮的周年项受到目前难以精确建模的水文负荷等环境因素的影响较大，因此本文仅计算钱德勒周期的极潮潮汐因子.
对刚体地球，理论重力极潮的计算公式(Wahr, 1985)为：
$g_{\mathrm{PM}}=R \mathit{\Omega}^2 \sin 2 \theta[x(t) \cos \lambda+y(t) \sin \lambda], $
其中，R=6.371×106 m为地球平均半径；Ω=7.292115×105 rad/s为地球平均自转角速度；θ和λ分别为台站所在地的余纬和经度；x(t)和y(t)分别为t时刻该点的极坐标，这里采用的是国际地球自转服务中心(International Earth Rotation Service，IERS)提供的EOPC04极坐标序列代入式(5)算得各台站对应的采样率为1 d的理论重力极潮.
虽然NMWT方法和极值点对称延拓策略可以将重力极潮的两个主要周其成分依频率分离，但极移是一个整体，所以通常在计算重力极潮潮汐因子时需要对这两个主要成分之和(图 7)进行整体最小二乘估计.分别将理论重力极潮的滤波结果、重力残差1的滤波结果看作两个周期分别为365.25天和432天的正(余)弦信号之和：
$\left\{\begin{array}{l}\Delta g_{\mathrm{th}}=A_{\mathrm{a}}^{\mathrm{th}} \cos \left(\omega_{\mathrm{a}} t+\varphi_{\mathrm{a}}^{\mathrm{th}}\right)+A_{\mathrm{c}}^{\mathrm{th}} \cos \left(\omega_{\mathrm{c}} t+\varphi_{\mathrm{c}}^{\mathrm{th}}\right) \\\Delta g_{\mathrm{re}}=A_{\mathrm{a}}^{\mathrm{re}} \cos \left(\omega_{\mathrm{a}} t+\varphi_{\mathrm{a}}^{\mathrm{re}}\right)+A_{\mathrm{c}}^{\mathrm{re}} \cos \left(\omega_{\mathrm{c}} t+\varphi_{\mathrm{c}}^{\mathrm{re}}\right)\end{array}, \right.$
其中，Δgth、Δgre分别是理论重力极潮的滤波结果和重力残差1的滤波结果，上角标th、re分别对应滤波结果来自理论重力极潮和残差1，下角标a、c分别对应重力极潮的周年项和钱德勒项，$\omega_{\mathrm{a}}=\frac{2 \pi}{365.25}\left(\mathrm{rad} \cdot \text { day }^{-1}\right)$、$\omega_c=\frac{2 \pi}{432}\left(\mathrm{rad} \cdot \text { day }^{-1}\right)$.利用最小二乘方法，可以估算出未知参数Aath、Acth、Aare、Acre、φath、φcth、φare和φcre，进而代入式(7)估算出如表 3所示的钱德勒周期的极潮潮汐因子：
$\left\{\begin{array}{l}\delta=A_{\mathrm{c}}^{\mathrm{re}} / A_{\mathrm{c}}^{\mathrm{th}} \\\kappa=\varphi_{\mathrm{c}}^{\mathrm{re}}-\varphi_{\mathrm{c}}^{\mathrm{th}}\end{array} .\right.$
表3 基于NMWT和极值点延拓策略、Daubechies小波滤波和EEMD得到的各台站Chandler周期的极潮潮汐因子

Table 3 Pole tide gravimetric factors at the Chandler period at each station derived from NMWT and extreme point extension strategy, Daubechies wavelet filtering and EEMD

台站 方法 δ κ/(°)
CB NMWT+EPE 1.1556±0.0092 0.0351±1.5009
Daubechies 1.0974±0.0144 11.1854±1.2708
EEMD 1.2677±0.0122 16.7882±1.0961
MB NMWT+EPE 1.2295±0.0103 8.4993±0.6281
Daubechies 1.1400±0.0152 7.7133±1.0466
EEMD 1.3165±0.0221 -2.7277±1.8523
MC NMWT+EPE 1.1191±0.0099 0.9178±0.8698
Daubechies 1.1103±0.0127 4.5842±1.0548
EEMD 1.3232±0.0216 -2.2010±2.0141
MO NMWT+EPE 1.1648±0.0123 0.0087±0.9882
Daubechies 1.0599±0.0141 5.3561±1.0155
EEMD 1.2874±0.0201 0.8317±1.3013
ST NMWT+EPE 1.2628±0.0338 2.7441±0.4809
Daubechies 1.2852±0.0424 19.3586±0.7759
EEMD 1.1202±0.0831 -3.0710±1.0014
WE NMWT+EPE 1.2010±0.0121 16.4063±0.5642
Daubechies 1.0890±0.0182 19.4593±0.8266
EEMD 1.2602±0.0402 13.1109±1.4165
WU NMWT+EPE 1.2028±0.0137 12.4857±0.6156
Daubechies 1.0756±0.0173 4.6767±1.1960
EEMD 1.2128±0.0637 9.1753±1.1708
此外，为说明NMWT结合极值点对称延拓策略所得结果的准确度，表 3还特别添加了分别对Daubechies小波滤波结果(256~512 d频段，详见Hu等(2007)和徐华君等(2008))和EEMD滤波的结果(对式(5)模型的滤波结果取其IMF6、7两个频段之和，对残差1的滤波结果取其IMF7、8两个频段之和，详见孙和平等，2018)进行的上述最小二乘估计.
从表 3中的结果可以看出，三种方法的结果大同小异.NMWT和极值点对称延拓策略对钱德勒周期极潮潮汐因子的估计精度普遍略高于Daubechies小波滤波和EEMD的结果，这有可能是Daubechies小波滤波的边缘效应影响了振幅因子的估计，此外，由于Daubechies小波和EEMD对重力残差2的滤波结果在频带上比NMWT宽(Daubechies小波滤波结果在256~512 d频带，EEMD的滤波结果取IMF6、7或IMF7、8之和，而NMWT得到的结果在频域上更集中于周年和钱德勒周期)，因此前二者包含了该频带的其他微小的干扰信号，这对振幅与相位的估计也产生了一定影响.

3.2.2 2~5 a信号

对于周期为2.3 a、2.7 a、3.3 a、3.7 a和4.8 a的5个信号，它们的周期接近于日长变化序列中识别到的2.3 a(2.4 a)、3.3 a、3.7 a和4.8 a信号，大气角动量数据(风向部分)被认为是这些日长信号的主要激发源(孔昭洋等，2021；张欣峰和刘根友，2020；廖德春，2000).Chen等(2019)在日长变化序列中提取出了3.65 a和4.9 a信号，并认为这在很大程度上可以由本地或非本地的大气、海洋和水文效应来解释，本研究发现的3.7 a和4.8 a重力信号在周期上与之非常相近，因此极有可能与之对应；此外，Chen等(2019)在日长变化序列里也检测到2.5 a信号，本研究获得的2.3 a和2.7 a信号可能与之相对应.
上述2~5 a尺度的重力信号因为相似的周期而在大气、海洋和水文效应中能找到对应，而超导重力仪观测数据的残差是由重力观测数据扣除各物理改正项的模型所得，因此残差中2~5 a尺度的信号可能由大气、海洋的重力模型误差和未扣除的水文负荷引起，从而在重力残差中得以体现.
以Canberra站为例，在2~5 a尺度上探究大气负荷、海潮负荷与海洋非潮汐负荷之和(后面简称海洋负荷)、水文负荷三者与残差的关系.其中，大气负荷是超导重力仪Level3数据中的大气改正项，海潮负荷由gotic2软件计算，水文负荷是EOST Loading Service提供的ERA5(ECMWF再分析)数据(loading.u-strasbg.fr).
首先，通过对比大气负荷、海洋负荷和水文负荷三者的傅里叶频谱(图 8)，可知在2~5 a尺度上，水文负荷的谱峰非常明显，而大气与海洋负荷的影响极小(海洋影响小于1 nm/s2，大气的影响小于2 nm/s2)，因此大气与海潮负荷对这一尺度的贡献不显著，水文负荷的贡献较大.虽然不能排除海洋负荷与大气负荷模型误差的影响，尤其是由观测不完备造成的模型误差，但这种模型误差的估计非常困难，因此目前只能研究量级较大的水文负荷对重力残差的影响.
图8 Canberra站大气、海潮与水文负荷数据及其傅里叶频谱

Fig 8 Atmospheric, oceanic tidal and hydrographic loading data and their Fourier spectra at Canberra station

为了从残差中选取出2~5 a频带，使其避免受其他相邻频带干扰，遂用正交Daubechies小波(消失矩阶数为14)将本台站的水文负荷(图 9a)分为14层，取d9和d10(dN所对应的频带为2N~2N+1天，这两层之和表示1.4~5.6 a频带的信号)，并用NMWT方法剔除1.8 a、5.3 a和4.2 a的干扰信号后(图 9b、c)，可在傅里叶频谱上看到非目标频带信号很好的剔除效果(红色虚线为水文负荷的全频傅里叶频谱，蓝色实线为目标频带的傅里叶频谱)，并在NMWT时频谱中发现在2~5 a的时间尺度上有着4个明显的周期分别为2.3 a、2.7 a、3.3 a、4.8 a的信号和微弱的3.7 a信号(如图 9d所示，为了同时展现NMWT频谱中的多个显著信号，这里采用未经过振幅改正的NMWT频谱).
图9 Canberra站水文负荷2~5 a频带信号和NMWT时频系数谱

Fig 9 Signal and NMWT time-frequency coefficient spectra for hydrologic loading in bands 2~5 a at Canberra station

对水文负荷中的上述信号进行提取，结果如图 10所示：红色实线为水文负荷中的信号，蓝色实线为残差中的相应信号，水文负荷中的2.3 a、2.7 a、3.3 a、3.7 a和4.8 a信号在相位上与残差中相应的信号近似一致或近似相反，水文负荷与残差中2.3 a、2.7 a、3.3 a、3.7 a和4.8 a的相关系数分别为-0.83，0.85，0.90，0.72，-0.83(标注在图 10各子图右侧)，这说明了该站残差中的2~5 a信号与水文负荷有着较高的相关性；至于为什么残差中不同的年际信号与水文负荷的相应年际信号的相关系数的符号不同，这可能涉及复杂的水文局部效应，有待进一步的研究.
图10 Canberra站2~5 a频带水文负荷信号(红)和残差信号(蓝)对比

Fig 10 Comparison of hydrologic loading signal (red) and residual signal (blue) in band 2~5 a at Canberra station

其他台站的重力残差2~5 a信号与相应的水文负荷信号在振幅上处于同一数量级，按照Canberra站的操作，计算残差2~5 a信号与相应的水文负荷信号的相关系数，结果如表 4所示.
表4 各台站重力残差2与水文负荷中2~5年周期信号的相关系数

Table 4 Correlation coefficients between the 2~5-year cycle signals extracted from the gravity residuals 2 and those extracted from the hydrologic loadings at each station

CB MB MC MO ST WE WU
2.3 -0.83 -0.74 0.63 -0.86 0.69 0.78 -0.18
2.7 0.85 0.82 0.89 -0.12 0.73 0.67 0.94
3.3 0.90 -0.82 -0.73 0.70 0.86 -0.54 -0.95
3.7 0.72 0.70 0.61 0.75 0.61 0.43 0.68
4.8 -0.83 -0.85 0.93 -0.95 0.57 0.62 -0.92
由表 4可看出，残差2中的2~5 a各信号与水文负荷的相关系数表现出较大区域性差异，除了个别台站在某些周期上(WU站的2.3 a信号、MO站的2.7 a信号)的相关系数接近0，总体上呈现出较高相关性.对于为何不同台站在相同周期信号上与水文负荷表现出相反的相关性(比如某两个台站在某一周期信号上表现出符号相反的相关系数)，原因可能是与各台站所处的海拔与当地海拔或地下水位的相对关系有关，且表现出各频段、区域性的差异性；对于同一台站的不同周期信号与水文负荷表现出相关系数差异过大甚至符号相反，需要进一步的研究.

3.2.3 6 a信号

对于本文提取出的6 a重力信号，除WE站为13.1 nm/s2，其余振幅在5 nm/s2左右，周期均处于5.7~5.95 a；其中，位于欧洲的相邻5个台站MB、MC、MO、ST和WE在相位上较为相近(图 6g).
6 a信号在地表位移观测中也有所体现：Ding和Chao(2018)从扣除了固体潮、海潮、极移与一些季节性和年际影响后的GNSS观测地表位移中检测到一个空间分布符合2阶2次扇形球谐相位模式(Y22)的6年周期性变化，并为此提供了一个如公式(8)所示的经验模型(Ding et al., 2020)：
$u_{\mathrm{U}}=\frac{a}{4} \sqrt{\frac{15}{2 \pi}} \sin ^2 \theta \cos \left(2 \pi f t+\varphi_1+2 \lambda\right), $
其中，uU为地表垂直位移，参数a=4.37 mm和φ1=1.119 rad分别是垂直位移序列的振幅估计参数和初始相位，θ和λ分别是所用超导重力仪台站的余纬和经度.
为探究两类观测的6 a信号的联系，利用此经验模型，结合自转地球上一阶垂直形变和重力位的关系，以及地表重力变化与重力位之间的关系(公式(9))(Wahr, 1985)，得到式(10)所示的6 a垂直位移伴随的地表重力变化序列.并将两类重力信号进行对比和相关性分析.公式(9)、(10)为：
$\left\{\begin{array}{l}u_{\mathrm{U}}=\Delta V \cdot h / g \\\Delta g=-2 \delta \Delta V / R\end{array}\right.,$
$\Delta g=-\sqrt{\frac{15}{8 \pi}} \frac{g a \delta}{h R} \sin ^2 \theta \cos \left(\frac{2 \pi t}{T}+\varphi_1+2 \lambda\right),$
其中h和δ分别是垂直位移和重力勒夫数，分别等于0.61和1.16(Dehant et al., 1999)，Δg是地表重力变化，ΔV是重力位变化，g=9.80 m/s2是重力加速度，R=6371.393 km是地球平均半径.
图 11中红色实线是根据公式(10)得到6 a垂直位移所致的地表重力变化序列，蓝色实线是NMWT从残差2中检测出的6 a信号.
图11 各台站残差2中6 a信号与GNSS观测U向6 a位移所致的地表重力变化模型对比

Fig 11 Comparison of the 6 a gravity signal in the residuals 2 of observations from each station and the surface gravity change due to GNSS observations of U-direction 6 a displacement model

在振幅方面，6 a重力模型决定了纬度相近的台站具有相近的振幅，而本研究所用的所有欧洲的台站(MB、MC、MO和ST台站)的提取结果振幅都在5 nm/s2左右，WE台站的振幅略大，为13.1 nm/s2，考虑这可能与数据预处理得到Level 3数据的过程相关，因为WE站Level 1数据存在较多间断、阶跃等，并且阶跃和间断往往伴随发生，而对阶跃的处理相比其他预处理过程而言，主观性较强，因此可能会引入新的阶跃，从而影响了6 a乃至所有年际信号提取结果的准确性.在相位方面，NMWT方法提取出的6 a信号的结果与6 a垂直位移的相位较为相近(MC的两类重力信号相位差最大，约为0.78 a=0.27π)，即NMWT方法从残差中检测的6 a信号在相位上与地表垂直位移的2阶2次扇形球谐相位模式吻合得较好.对各站的两类重力信号进行相关性分析，可得CB、MB、MC、MO、ST和WU对应的相关系数分别是0.94、0.90、0.81、0.98、0.97与0.92，表明重力6 a信号与地表垂直位移可能同源.

4 结论

本文采用了Canberra等7个超导重力仪观测数据，扣除了固体潮、海潮负荷、大气负荷影响和仪器漂移后获得的超导重力仪重力残差，使用NMWT方法识别和提取出了重力极潮，并基于此利用最小二乘方法计算出钱德勒周期的极潮潮汐因子，精度比Daubichies小波滤波和EEMD方法所得的结果略高.对上述残差扣除重力极潮后，用NMWT方法识别和提取出了若干较为稳定的年际信号，这些信号的周期集中出现在2.3 a、2.7 a、3.3 a、3.7 a、4.8 a和6 a，通过上述信号的时域提取结果计算获得了相应的平均振幅.其中2~5 a尺度的变化与水文负荷有着较高的相关性，但也不能排除大气和海潮负荷模型误差的影响；6 a信号与GNSS观测地表垂直位移模型的相关性较高，与垂直位移可能同源.重力变化是地球整体运动的综合反映，物理模型的精确计算与扣除仍然是一个有待完善的问题，更长的观测数据有助于年际变化的精确提取和物理解释.

感谢陈晓东副研究员提供了Wuhan站1997年12月至2021年3月的h1数据，IGETS提供了超导重力仪Level3数据，IGETS的Dr. Christoph Förste将超导重力仪数据库中Wuhan站部分错误的Level2数据予以及时修正，EOST Loading Service提供了地表重力载荷数据.

References

Publishing order | Descend order by publishing year | Descend order by cited within
Cai S , Liu L T , Wang G C . Short-term tidal level prediction using normal time-frequency transform. Ocean Engineering. 2018, 156 489-499

DOI

Chao B F . Dynamics of the inner core wobble under mantle-inner core gravitational interactions. Journal of Geophysical Research: Solid Earth. 2017, 122(9 7437 7448

DOI

Chen J L , Wilson C R , Kuang W J , et al. Interannual oscillations in earth rotation. Journal of Geophysical Research: Solid Earth. 2019, 124(12): 13404-13414

DOI

Chen X D , Ducarme B , Sun H P , et al. Loading effect of a self-consistent equilibrium ocean pole tide on the gravimetric parameters of the gravity pole tides at superconducting gravimeter stations. Journal of Geodynamics. 2008, 45(4-5): 201-207

DOI

Dehant V , Defraigne P , Wahr J M . Tides for a convective earth. Journal of Geophysical Research: Solid Earth. 1999, 104(B1): 1035-1058

DOI

Ding H , An Y C , Shen W B . New evidence for the fluctuation characteristics of intradecadal periodic signals in length-of-day variation. Journal of Geophysical Research: Solid Earth. 2021, 126(2

DOI

Ding H , Chao B F . Solid pole tide in global GPS and superconducting gravimeter observations: signal retrieval and inference for mantle anelasticity. Earth and Planetary Science Letters. 2017, 459 244-251

DOI

Ding H , Chao B F . A 6-year westward rotary motion in the Earth: Detection and possible MICG coupling mechanism. Earth and Planetary Science Letters. 2018, 295 50-55

Ding H , Xu X Y , Pan Y J , et al. A time-varying 3-D displacement model of the 5.9-year westward motion and its applications for the global navigation satellite system positions and velocities. Journal of Geophysical Research: Solid Earth. 2020, 125(4

DOI

Duan P S , Huang C L . Intradecadal variations in length of day and their correspondence with geomagnetic jerks. Nature Communications. 2020, 11(1): 2273

DOI

Duan P S , Liu G Y , Hu X G , et al. Possible damping model of the 6 year oscillation signal in length of day. Physics of the Earth and Planetary Interiors. 2017, 265 35-42

DOI

Ducarme B , Venedikov A P , Arnoso J , et al. Global analysis of the GGP superconducting gravimeters network for the estimation of the pole tide gravimetric amplitude factor. Journal of Geodynamics. 2006, 41(1-3): 334-344

DOI

Gillet N , Jault D , Canet E , et al. Fast torsional waves and strong magnetic field within the Earth's core. Nature. 2010, 465(7294): 74-77

DOI

Holme R , de Viron O . Characterization and implications of intradecadal variations in length of day. Nature. 2013, 499(7457): 202-204

DOI

Hu X G , Liu L T , Ducarme B , et al. Estimation of the pole tide gravimetric factor at the chandler period through wavelet filtering. Geophysical Journal International. 2007, 169(3): 821-829

DOI

Hu X G , Liu L T , Hinderer J , et al. Wavelet filter analysis of local atmospheric pressure effects on gravity variations. Journal of Geodesy. 2005, 79(8): 447-459

DOI

Huang C L , Jin W J . The effect of the variation of earth rotation on gravimetry. Astronomical Techniques and Instruments. 1999 3): 41-43

DOI

Kong Z Y , Zhou Y H , Xu X Q , et al. Correlation analyses among ΔLOD, AAM and ENSO, and the 2020—2021 La Nina Event. Progress in Astronomy. 2021, 39(4): 532-543

DOI

Li H , Xu J Q , Chen X D , et al. Extracting long-period surface waves and free oscillations using ambient noise recorded by global distributed superconducting gravimeters. Seismological Research Letters. 2020, 91(4): 2234-2246

DOI

Liao D C . Main excitation sources of ΔLOD on interannual time scale. Acta Astronomica Sinica. 2000, 41(2): 139-147

DOI

Liu L T , Hsu H , Grafarend E W . Normal Morlet wavelet transform and its application to the Earth's polar motion. Journal of Geophysical Research: Solid Earth. 2007, 112(B8): B08401

DOI

Mound J E , Buffett B A . Detection of a gravitational oscillation in length-of-day. Earth and Planetary Science Letters. 2006, 243(3-4): 383-389

DOI

Shen W B , Liu R L . Study of detecting the inner core super rotation using SG data. Geomatics and Information Science of Wuhan University. 2009, 34(1): 72-76

Su X Q , Liu L T , Houtse H , et al. Long-term polar motion prediction using normal time-frequency transform. Journal of Geodesy. 2014, 88(2): 145-155

DOI

Sun H P , Xu J Q , Ducarme B . Experimental earth tidal models in considering nearly diurnal free wobble of the Earth's liquid core. Chin. Sci. Bull.. 2003, 48(9): 935-940

Sun H P , Zhang M M , Xu J Q , et al. Etraction and study of the gravity pole tide using EEMD. Chin. J. Geophys.. 2018, 61(2): 521-530

DOI

Sun H P , Zhou J C , Xu J Q , et al. High-precision superconducting gravimetric observations and investigations provide national surveying and mapping and earth's dynamics with theoretical fundamentals. Bulletin of Chinese Academy of Sciences. 2021, 36(2): 216-223

DOI

Wahr J M . Deformation induced by polar motion. Journal of Geophysical Research: Solid Earth. 1985, 90(B11): 9363-9368

DOI

Wang G C , Liu L T , Xu A G , et al. On the capabilities of the inaction method for extracting the periodic components from GPS clock data. GPS Solutions. 2018, 22(4): 92

DOI

Wang L S , Chen C , Du J S , et al. Noise reduction, atmospheric pressure admittance estimation and long-period component extraction in time-varying gravity signals using ensemble empirical mode decomposition. Terrestrial Atmospheric and Oceanic Sciences. 2015, 26(2 Part 1): 111-120

Xu H J , Liu L T , Luo X W . Normal Morlet wavelet transform and its application to detection of the earth's free oscillations. Acta Geodaetica et Cartographica Sinica. 2009, 38(1): 16-21

DOI

Xu H J , Liu L T , Xu H Z , et al. Wavelet approach to study gravity pole tide. Geomatics and Information Science of Wuhan University. 2008, 33(11): 1114-1117

Xu J Q , Sun H P , Fu R S . Detection of long-period core modes by using the data from global superconducting gravimeters. Chin. J. Geophys. 2005, 48(1): 69-77

DOI

Xu J Q , Sun H P , Yang X F . A study of gravity variations caused by polar motion using superconducting gravimeter data from the GGP network. Journal of Geodesy. 2004, 78(3): 201-209

DOI

Yang X F , Sun H P , Xu J Q . Detection of gravitational effect of polar motion by use of auto-regression model. Journal of Geodesy and Geodynamics. 2004, 24(2): 110-114

DOI

Zhang X F , Liu G Y . Extraction of interannual signals in the length-of-day variation and correlation analysis with the atmosphere. Journal of Geodesy and Geodynamics. 2020, 40(11): 1188-1193

DOI

Ziegler Y , Hinderer J , Rogister Y , et al. Estimation of the gravimetric pole tide by stacking long time-series of GGP superconducting gravimeters. Geophysical Journal International. 2016, 205(1): 77-88

黄 乘利 , 金 文敬 . 地球自转变化的重力效应. 云南天文台台刊. 1999, 3 41-43

DOI

孔 昭洋 , 周 永宏 , 许 雪晴 , et al. 日长变化、大气角动量和ENSO年际信号的相关分析及2020—2021年拉尼娜事件. 天文学进展. 2021, 39(4): 532-543

DOI

廖 德春 . 日长年际变化的主要激发源. 天文学报. 2000, 41(2): 139-147

DOI

申 文斌 , 刘 任莉 . 利用超导重力数据探测内核超速旋转的研究. 武汉大学学报·信息科学版. 2009, 34(1): 72-76

孙 和平 , 徐 建桥 , Ducarme B . 基于全球超导重力仪观测资料考虑液核近周日共振效应的固体潮实验模型. 科学通报. 2003, 48(6): 610-614

DOI

孙 和平 , 张 苗苗 , 徐 建桥 , et al. 基于集合经验模态分解法的重力极潮提取与研究. 地球物理学报. 2018, 61(2): 521-530

DOI

孙 和平 , 周 江存 , 徐 建桥 , et al. 高精度超导重力观测与研究为国家精密测绘和全球地球动力学提供理论基础. 中国科学院院刊. 2021, 36(2): 216-223

DOI

徐 华君 , 柳 林涛 , 罗 孝文 . 标准Morlet小波变换检测地球自由振荡. 测绘学报. 2009, 38(1): 16-21

DOI

徐 华君 , 柳 林涛 , 许 厚泽 , et al. 重力极潮的小波分析. 武汉大学学报·信息科学版. 2008, 33(11): 1114-1117

徐 建桥 , 孙 和平 , 傅 容珊 . 利用全球超导重力仪数据检测长周期核模. 地球物理学报. 2005, 48(1): 69-77

DOI

杨 学峰 , 孙 和平 , 徐 建桥 . 利用自回归模型检测地球极移重力效应. 大地测量与地球动力学. 2004, 24(2): 110-114

DOI

张 欣峰 , 刘 根友 . 日长年际信号提取及其与大气相关性分析. 大地测量与地球动力学. 2020, 40(11): 1188-1193

DOI

Options
Outlines

/

〈 〉