0 引言

重力随时间的变化可分为潮汐变化和非潮汐变化，前者如地球潮汐与地球自转变化，后者一般是地球物质迁移引起的重力随时间的变化(王谦身，2003)。连续重力观测是在台站固定点上进行的连续观测，获得的台站重力变化是台站高程以及台站和其周围地壳物性特征随时间变化的结果。连续重力观测能够精确探测到地球系统各圈层物质迁移引起的重力变化效应，包括叠加在一起的长期(多年的)或短期(季节性的)的变化(Sun et al，2001、2009、2011;周江林等，2015;翟丽娜等，2020)。连续重力观测获得的数据由重力潮汐变化、重力非潮汐变化和误差构成，即

$ g=g_{\text {tide }}+g_{\text {non-tide }}+\mathit{\epsilon}(t) $

(1)

重力潮汐变化主要是日月的引潮力，而重力非潮汐变化主要由观测仪器本身的零漂、自然环境干扰、人为因素和残差所构成，即

$ g_{\text {non-tide }}=g_{\text {linear-drift }}+g_{\text {environment }}+g_{\text {human }}+\mathit{\epsilon}(t) $

(2)

自然环境干扰引起的重力变化主要是由气压负荷、地下水、极移、海潮等引起，即

$ g_{\text {environment }}=g_{\text {air-pressure }}+g_{\text {polar-shift }}+g_{\text {groundwater }}+g_{\text {ocean-tide }}+\mathit{\epsilon}(t) $

(3)

为获得重力非潮汐变化，需对重力潮汐变化数据进行零漂、固体潮、气压和海潮等的改正，本文所要研究的就是分别定量分析气压和海潮对琼中台重力非潮汐变化的影响，获得气压和海潮负荷改正后的琼中台重力变化。海南岛位于中国大陆南端，四面环海，琼中地震台(109.8°E，19.0°N，属于热带地区，地理位置见图 1)位于海南岛中部，2008年3月开始连续重力观测，运行仪器为PET重力仪(美国Micro-g LaCoste，Inc公司生产)，该重力仪是一种全自动的金属弹簧(零长弹簧)连续重力观测仪，是在LCR-G型基础上改进而成，工作原理与LCR-G相同。琼中台距离最近的海岸线约60km，台站受热带风暴、台风和强降雨气候影响频繁，属于热带季风气候，台风季节集中在夏末秋初，8至10月是台风多发期，气压变化具有典型的热带特征(李盛等，2016;李盛等，2020)。基于琼中台这一鲜明的地理环境和气候特征，对琼中台重力非潮汐变化进行气压与海潮负荷改正很有必要，可准确认识琼中台重力变化，充分发挥琼中台重力基准作用。

1 原始资料分析

对原始数据进行分析，可掌握仪器的运行情况，为后续的非潮汐分析打下基础。对琼中台2009年1月至2019年9月连续重力原始数据进行台阶、突跳等处理，获得该段时间琼中台重力潮汐变化数据(图 2(a))，由图可知琼中台重力潮汐变化总体比较稳定，固体潮清晰，但漂移存在非线性的情况(数据采样率为分钟值)。

2 潮汐分析

1966年保加利亚的Venedikov教授提出了利用48阶奇偶滤波器采用滤波方式分别计算周日、半日和1/3日波的潮波分量，并利用上述分量进行潮汐参数的估计(Venedikov et al，2003)。目前该方法被国家重力台网中心^①在结合Nakai检验方法和数字滤波方法(李辉等，1994)基础上用于进行连续重力观测资料的评价工作。2003年，Venedikov进一步改进了程序，引入日本京都大学Tamura等(1991) 1200个分波的潮波表，将气压等辅助观测序列、分段多项式等方法引入潮汐分析数学模型中，使得程序能够对任意潮汐观测数据进行分析，即为VAV潮汐分析软件。

① http://www.gncc.ac.cn

运用VAV潮汐分析软件，对2009年1月至2019年9月琼中台连续重力观测数据进行潮汐分析，得到琼中台重力潮汐变化观测数据的各主要潮波参数(表 1)。其中，M2波(月亮的主半日波)的振幅最大，S2波(太阳的主半日波)次之。实际上，在各潮波中，M2波振幅最大，也最为稳定，因此在日常分析中主要跟踪M2波潮汐因子的变化。

3 重力非潮汐变化改正 3.1 零漂改正

弹簧重力仪的零漂与弹性材料以及仪器的设计方式有着直接的关系。如果考虑弹性系统的弹性滞后、蠕变变形、弹性后效、弹簧年龄、内部传感器温度和弹簧的应力加载等因素影响，可以模拟弹性系统的零漂。但由于上述因素的差异性以及建模的困难，本文选择一般多项式拟合方法模拟弹黃重力仪的零漂。

根据琼中台连续重力观测数据的处理情况，将2009年1月至2019年9月重力潮汐变化数据分段，然后对分段数据的零漂进行一般多项式分段曲线拟合(图 2(b))，即可获得琼中台重力潮汐变化数据的零漂拟合值以及零漂改正值(图 2(c))。

3.2 固体潮改正

采用国际地潮中心推荐的Tsoft固体潮数据预处理软件作为观测数据预处理的分析软件(van Camp et al，2005)。该软件可计算基于WDD模型(非流体静力平衡的地球模型)的理论固体潮值。运用该软件，计算可得琼中台理论固体潮，对零漂改正后琼中台重力残差(图 2(b))进行理论固体潮改正，得到零漂固体潮改正后的重力残差曲线(图 3)。

3.3 气压改正

气压是量级仅次于固体潮的重力影响因素。重力固体潮数据中的气压影响主要是由大气质量变化引起的负荷效应(直接负荷效应)和由大气负荷形变引起的地表重力场变化(间接负荷效应)组合而成(罗少聪等，2000)。其中直接负荷效应的影响和本地的气压资料密切相关，可通过数值模拟的方法进行模拟(王迪晋，2009)，是本地重力影响的主要因素。气压对重力的影响还具有频率依赖性。根据目前的研究，气压对重力可影响到频段7cpd的重力信号(韦进等，2011)。因此，为了尽可能地消除气压对重力的影响，建立了既依赖于时间变化又依赖于频率变化的重力气压导纳值，以消除重力观测值中的气压噪声(Hu et al，2006)。气压和重力之间的变化规律可用大气重力导纳值(Admittance Factor)来描述，其基本数值在(-0.25~-0.40)×10^-8m/s²/mbar，表 2列举了不同机构给出的气压导纳值结果(李英冰，2003;孙和平等，2002)。

为了对琼中台重力非潮汐变化进行气压改正，运用VAV软件，计算了琼中台2009年1月至2019年9月每年的大气重力导纳值，计算结果包含了各频段(周期1~7cpd)的大气重力导纳值(表 3)。由表 3可知，2009年1月以来，琼中台每年各频段的气压导纳值介于(-0.15~-0.49)×10^-8m/s²/mbar之间，11年中有5年的气压导纳值位于(-0.30~-0.40)×10^-8m/s²/mbar之间。而若整体计算2009年1月至2019年9月的气压导纳值，则该值为-0.34×10^-8m/s²/mbar，与表 2中的结果接近。说明短期来看，琼中台的大气重力导纳值有所差异，但长期而言，琼中台的大气重力导纳值与其他学者研究不同区域的结果较为一致。

基于上述结果，取2009年1月至2019年9月的大气重力导纳值，即-0.34×10^-8m/s²/mbar，对每年的重力非潮汐变化进行改正，最终获得了气压改正后的琼中台重力非潮汐变化(以下简称气压改正)。

3.4 海潮改正

海洋潮汐是海水在太阳和月亮的引潮力作用下周期性涨落的自然现象。海洋潮汐是重力固体潮观测的重要干扰源。研究表明，在沿海地区重力固体潮中的海潮负荷一般可达几十微伽量级。为有效利用重力资料研究地球物理学和地球动力学问题，对重力观测资料的海潮负荷信号进行改正就显得特别重要(孙和平等，2002)。Farrell(1972)给出了计算海潮负荷的负荷格林函数，并指出海潮负荷即海潮潮高和负荷格林函数的褶积积分

$ L(\theta, \lambda, t)=R^{2} \iint \rho H(\varphi, A, t) G(\varphi) \sin \varphi \mathrm{d} \varphi \mathrm{d} A $

(4)

式中，R为地球半径，ρ为海水密度，θ为计算点的余纬，λ为计算点的经度，(φ，A)为模板坐标，t为时间，G(φ)为重力负荷格林函数。通常，负荷结果均以余弦波的振幅Amp和相位Φ给出，即任一时刻的负荷值为

$ L(\theta, \lambda, t)=A \operatorname{Amp} \cdot \cos (\varPhi+\varphi) $

(5)

因此，总结海潮负荷包括三方面内容: ①海潮模型;②格林函数;③计算方法。

由于存在多种海潮模型，因此需要对其进行选择。海潮模型对重力负荷计算的影响具有明显的地域性特点(孙和平等，2005)，全球区域没有绝对优越的海潮模型。研究表明，近海潮汐对沿海台站的影响非常大(Sun，1992;周江存等，2004)。目前应用较为广泛的有CSR4.0、NAO.99b、FES2004、GOT04、TPXO7.2、DTU10、EOT11a、HAMTIDE11a等8个海潮模型(Kim et al，2011、2013)。

针对众多海潮模型，Yu(2006)认为NAO.99b为最好的潮汐模型，该模型是根据日本天文台应用T/P卫星数据，结合日本验潮站资料构建的区域模型，精度为0.25°×0.25°。因此，本文基于SPOTL(Agnew，1997)程序，以NAO.99b潮汐模型计算了琼中台海潮负荷值，并运用该负荷值对琼中台重力非潮汐零漂、固体潮、气压改正后的数据进行海潮改正(以下简称海潮改正)。将零漂固体潮改正与海潮改正后的琼中台重力非潮汐变化数据进行对比(图 4)，认为海潮改正后的重力残差值变化幅度约为35×10^-8m/s²，大于仅进行潮汐、零漂改正后的重力残差值变化幅度。

运行VAV软件，分别对零漂固体潮改正、气压改正、海潮改正后的琼中台重力非潮汐变化数据进行Venedikov调和分析，获得以上数据的各主要潮汐波参数，并将其与所有改正前的琼中台重力潮汐变化数据的主要波群参数进行比较，结果表明经海潮改正后的各潮波振幅明显小于仅进行零漂固体潮改正的振幅，其中变化幅度最大的为M2波，其振幅由1.70438×10^-8m/s²降至0.44651×10^-8m/s²。因此，进行海潮负荷改正可进一步消除重力潮汐变化数据中的潮汐信号，提高提取重力非潮汐变化的精度。

3.5 小结

根据以上计算的结果，将琼中台重力非潮汐变化气压与海潮负荷改正值绘制在同一图中(图 5)，可知气压改正值约为10×10^-8m/s²，海潮改正值约为5×10^-8m/s²。

从2012年以来琼中台海潮改正后的重力非潮汐年变化来看(图 6)，2015年以来年变规律基本呈春冬低、夏秋高的变化趋势，其中6—7月为一年中重力最高值，12—1月为重力最低值，2012年以来琼中台海潮改正后重力非潮汐变化未出现趋势变化的异常情况。

4 结论

通过对2009年1月至2019年9月琼中台重力潮汐变化数据的潮汐分析以及非潮汐分析，得到如下认识:

(1) 琼中台的大气重力导纳值总体上与不同学者对不同区域研究所得的结果较为一致，未来对琼中台重力非潮汐变化数据进行气压改正时，可直接采用-0.34×10^-8m/s²/mbar作为大气重力导纳值。

(2) 杨锦玲等(2016)运用8个海潮模型计算厦门台海潮负荷，经对比，显示经8个全球海潮模型改正后，主要潮波的海潮负荷振幅差异变小。因此，对琼中台重力非潮汐变化的海潮改正选取NAO.99b潮汐模型是合适的。计算结果表明琼中台受到的海潮负荷变幅约5×10^-8m/s²，小于厦门的海潮负荷变幅7.3×10^-8m/s²，这应与琼中台距离海洋稍远有关(厦门台距离海岸线不到1km)。

(3) 琼中台重力非潮汐变化气压改正变幅约为10×10^-8m/s²，海潮改正变幅约为5×10^-8m/s²，气压改正变幅大于海潮改正变幅。经海潮改正后的琼中台重力非潮汐变化数据中的潮汐信号比仅进行零漂固体潮改正更加微弱，说明进行海潮改正具有一定效果，进一步消除了重力残差中的潮汐信号，有助于更准确地提取地球内部引起的重力变化，以获得地球内部动力学信息，为地震预测服务。海潮改正后仍然残留潮汐信号，可能是因为在进行固体潮改正时使用的是理论值，且海潮改正基于的是海潮模型。因此计算的固体潮、海潮负荷与实际真实值仍有误差，这些均可能导致海潮改正后的重力非潮汐变化数据仍残留潮汐信号。

需要指出的是，以上重力非潮汐变化改正过程中，零漂改正具有一定的不确定性，主要是因为受PET/gphone相对重力仪自身结构特性的限制，很难完全真实地拟合仪器的零漂。尤其当仪器出现诸如停电或仪器故障等原因导致停测一段时间，待仪器恢复观测后，仪器的零漂非线性零漂特征明显，故真实拟合仪器零漂更为困难，以致在数据处理过程中，可能会把真实的重力变化异常当做零漂进行改正。