地球重力场变化与各种物理过程(如潮汐、构造活动、地球自由振荡等)紧密联系，国内外许多学者都试图通过重力潮汐的观测来探索地球内部的奥秘并解释与之相关的地球物理现象^{[1, 2, 3]}。1997年7月启动的“全球地球动力学计划”(Global Geodynamics Projects，GGP)为研究全球动力学现象，如地球固体潮、重力场与海洋大气耦合、地球自由核章动、地球自转变化等导致的重力效应，提供了全球超导重力观测资料^{[4, 5, 6]}。

超导重力仪被认为是精度最高，稳定性和连续性最好的重力仪，对重力变化的灵敏度极高 (在频段低于0.8 mHz时超过当前最好的宽频地震仪^[7])，适合长周期信号的观测，在自由振荡等研究中优势明显^{[8, 9, 10]}。美国GWR公司生产的iGrav超导重力仪是最新型的便携式相对重力仪，其基本原理与台站型超导重力仪(OSG)解释相同，体积更小更为便携。其标称精度在频域中为1 nanoGal，在时域中1 min平均值为0.05 μGal，零漂小于0.5 μGal/月^[11]。

利用超导重力仪获取高精度地球重力信号，需进行各项必要的改正及残差分析。本文基于Tsoft软件^[12]处理武汉九峰台站连续5个月的超导重力观测数据，进行气压负荷、极移重力改正及尖峰、阶越、间断等干扰信号的修正；进一步分析重力残差及仪器的长期漂移特征，评估观测精度。采用功率谱密度估计(Spanos-Tratskas方法)^[13]，检测2013年Mw6.6级芦山地震激发的地球自由振荡球型模态₀S_m，并分析各频段上的信噪比。

武汉大学的iGrav-007超导重力仪引进于2012年，次年3月安装于中科院测地所武汉九峰台站。由于观测显示安装前3个月仪器尚未完全稳定，因此本文取2013年6月25日~11月25日连续5个月的原始秒采样观测数据，并采用最小二乘滤波获得分采样数据。

利用Tsoft软件进行数据预处理的过程包括：(1) 根据台站潮汐参数构制重力合成潮，得到重力残差；(2) 利用Tsoft多种修正方法去除尖峰、阶跃和间断等干扰信号；(3) 由回归分析法计算大气重力导纳值，改正大气负荷效应；(4) 根据IERS提供的地球自转参数，改正极移重力效应。对经大气负荷和极移重力改正后的重力残差分别进行多项式拟合和滤波，去除仪器漂移和残余潮汐信号，得到最终重力残差 ε(t)(单位为nm·s^-2)：

式中，t为采样历元；y_obs(t)为iGrav-007的重力潮汐观测值；y_mod(t)为台站重力合成潮；dg_P(t)为气压引起的台站重力变化；dg_pole(t)为极移引起的台站重力变化； 表示拟合仪器漂移和残余潮汐信号的p次多项式，a_i是多项式系数。

台站重力合成潮y_mod(t)可表示为：

式中，n为波群数；N为波群总数；δ_n和Δφ_n为第n个波群数的模型振幅因子和相位滞后；m为第n个波群的潮汐分量数，M_n为第n个波群的潮汐分量总数；H_nm、ω_nm和φ_nm为第n个波群第m个潮汐分量的模型振幅、角频率和初始相位^[2]。

本文采用台站提供的重力合成潮波群(共17个波群)参数，并置测站点的大地经、纬度λ=114.489 80°E，φ=30.515 87°N，仪器高程H=89.3 m，构制y_mod(t)，如图1(b)所示。

图 1 重力残差中的台站大气负荷和极移重力改正 Fig. 1 Correction of Air Pression and Polar Motion at Station on the Gravity Residuals

图选项

大气负荷改正的常用方法是将重力潮汐观测值减去重力合成潮获得重力残差y_resd(t)，如图1(c)所示，再将y_resd(t)与台站气压观测值P(t)进行最小二乘回归求得大气重力导纳值C，乘积CP(t)即为大气负荷改正量。

在时域内获得大气重力导纳值，一元线性回归模型为：

气压引起的台站重力变化dg_P(t)可表示为：

采用最小二乘线性迭代回归法，以残差标准差σ为中误差，3σ为限差，每次迭代时均去掉原始重力残差和气压观测值中的二次多项式拟合部分，当C的变化小于0.000 1 nm·s^-2/mBar时停止迭代，实际迭代4次即可，得到台站大气重力导纳值C=(-2.480 7± 0.003 5) nm·s^-2/mBar，斜距C₁= 4 361.244 8 nm·s^-2。由图1(d)可知，台站大气负荷的影响振幅约100 nm·s^-2。经大气负荷改正后，重力残差得到明显改善，如图1(e)所示，可见大气负荷效应是重力残差中的主要成分。

极移引起的地表任一点重力变化可根据该点坐标和极坐标的变化(从IERS提供的地球自转参数中得到)算得，故极移重力效应可作为已知量扣除^[14]。极移引起的台站重力变化dg_pole可表示为：

式中，地球平均半径 a=6 378 137.0 m；地球平均自转角速度 ω=7.292 115×10^-5 rad/s；(x，y)为测站点的极坐标；大地经度λ、大地纬度φ与上文给出的值相同。

根据式(5)算得武汉九峰台站的极移重力效应，并取与图1(e)中残差相同的时间段。因IERS提供的数据采样率为1 d，采用三次样条插值获得分采样的极移重力，如图1(f)所示，极移影响的振幅约20 nm·s^-2，较明显。

对最终重力残差线性拟合得iGrav-007的漂移率为0.40 μGal/月，如图2(a)所示，较其标称值小。采用20阶多项式拟合重力残差，如图2(b)所示，得到拟合残差均方根RMS=0.423 μGal/min (经实验，采用更高阶多项式RMS变化均小于0.01 μGal)。该值大于其标称值0.050 μGal/min。究其原因，一是因台站背景噪声较大；二是大气负荷改正不彻底，未考虑中远和远区的气压影响^[15]以及对重力潮汐观测影响的频率依赖，仅采用线性回归计算大气重力导纳值^{[15, 16, 17]}。

图 2 重力残差的线性拟合、20次多项式拟合及多项式拟合残差 Fig. 2 Linear Fitting, 20 Times Polynomial Fitting and Polynomial Residual Error on the Gravity Residuals

图选项

对最终重力残差分别进行最小二乘高通、带通、低通滤波，截止频率分别以合成潮M3波的最高频率3.081 254 cpd和MM波的最低频率0.020 885 cpd为界。此处认为比M3波频更高的部分由仪器的高频噪声引起，比MM波频更低的部 分为仪器漂移，中间部分主要为残余潮汐。故将重力残差分为短周期、潮汐频段、长周期三部分。由图3(a)可知，重力残差中存在微弱的残余潮汐信号，能量约在±3 nm·s^-2之间，相对较小。以短周期部分扣除平均值后的均方根0.040 μGal/min衡量其观测精度，表明iGrav-007观测性能良好。

图 3 重力残差的最小二乘滤波 Fig. 3 Least Squares Filtering on the Gravity Residuals

图选项

2013年芦山地震在2.3~5.3 mHz频段激 发了微弱的球型简正模信号，为研究iGrav-007的 高频特性提供了条件。采用功率谱密度估计，检测该地震激发的地球自由振荡球型模态₀S_m，分析各频段上的信噪比，并与PREM模型^[18]理论预测值比较，评价iGrav-007的残差水平。芦山地震发生于2013年4月20日0:2:48(UTC)，选取了地震后18 h的iGrav-007观测数据，如图4所示。

图 4 iGrav-007记录的芦山地震后18 h的原始观测数据/V Fig. 4 Hours Original Recording After Lushan Earthquake by iGrav-007 /V

图选项

在0.28~1.80 mHz频段上，自由振荡信号受噪声干扰严重，₀S₀~₀S₁₀基型振荡无法有效检测。分析原因，一是这些频率的低阶振型受地震强度较小的限制，难以激发；二是地球自转和椭率可导致自由振荡信号的分裂和耦合^{[19, 20]}。

图5(a)、图5(b)分别为iGrav-007在1.8~3.9 mHz和3.9~5.5 mHz频段检测到的芦山地震激发的自由振荡振型。数据长度为18 h，傅里叶变换采用Hanning窗。由图5可知iGrav-007清晰地检测到了₀S₁₁~₀S₄₈之间除₀S₁₆振型之外的振型基频系列。为对比观测到的自由振荡周期与PREM理论值，图5中用虚线指示PREM模型的理论频率值。₀S₁₆振型未形成突出的谱峰，主要是受附近的噪声影响较大所致。1.80~3.00 mHz频段信噪比较低，但功率谱密度能较好地反映基频球型振荡频率。在3.00~4.75 mHz频段检测到₀S₂₁~₀S₄₀的所有振型，信噪比均大于3，谱峰清晰，大部分基频振型观测值与PREM理论值偏差较小。其中₀S₃₀振型虽形成清晰的谱峰，但观测值与理论值频率偏差相对其他频段较大，达到20.0×10^-6 Hz(约8倍的傅氏谱分析频率分辨率)。究其原因，一是由于该振型附近噪声较大，二是这一基频振型可能受到地球上地幔剪切波的干扰^[21]，且₀S₃₁振型处可能出现谱线分裂现象。在4.75~5.50 mHz频段检测到₀S₄₁~₀S₄₈振型，该频段附近噪声较低，谱峰清晰可见，与PREM理论值较吻合。

图 5 iGrav-007记录到的芦山地震激发的自由振荡信号傅里叶变换频谱 Fig. 5 FFT Transform Spectrum of Earth’s Free Oscillation Excited by Lushan Earthquake Recorded by iGrav-007

图选项

由表1可知，1.80~5.50 mHz频段iGrav-007检测到的基频自由振荡信号和PREM理论值偏差仅3个基频振型偏差绝对值超过0.30%，偏差绝对值小于0.15%的检测振型数占总振型数的65.7%。

本文基于iGrav-007在武汉九峰台站连续5个月的观测数据，应用Tsoft软件进行大气负荷和极移重力效应改正，采用最小二乘线性迭代回归，在时域内得到武汉九峰台站大气重力导纳值为(-2.480 7±0.003 5) nm·s^-2/mBar。利用线性拟合最终重力残差，得到iGrav-007的漂移率为0.40 μGal/月；利用多项式拟合重力残差，去除仪器长期漂移量和残余潮汐信号，得到最终重力残差均方根为0.423 μGal/ min；利用滤波，以重力残差短周期部分扣除平均值后的均方根0.040 μGal/min衡量观测精度，表明iGrav-007观测性能良好。

iGrav-007很好地记录了2013年Mw6.6芦山地震，检测出1.80~5.50 mHz频段内₀S₁₁~₀S₄₀振型中除₀S₁₆振型外的基频振型自由振荡系列。观测值与PREM理论值偏差绝对值超过0.30%的有3个基频振型，小于0.15%的检测振型数占总振型数的65.7%。

下一步将对iGrav-007固体潮观测数据作调和分析，综合考虑固体潮、气压、极移、日长等变化的影响，获得更准确的潮汐参数和频域内的大气重力导纳值。

致谢: 感谢中国科学院测量与地球物理研究所提供的武汉九峰台站潮汐参数。