中国地震  2022, Vol. 38 Issue (1): 102-111
琼中台重力非潮汐变化气压与海潮负荷改正
李盛, 胡久常, 周雯, 张慧, 王慧琳, 谢小玲     
海南省地震局,海口 570203
摘要:本文对琼中台连续重力观测数据进行收集整理并处理,基于处理后的数据,进行了潮汐分析和非潮汐分析。潮汐分析采用VAV调和分析方法;非潮汐分析则分别进行了零漂改正、固体潮改正、气压改正和海潮改正。其中,零漂改正采用一般多项式拟合零漂的方法;气压改正采用VAV软件;海潮改正运用SPOTL程序,以NAO.99b潮汐模型计算了琼中台海潮负荷值。最终获得了改正后的琼中台重力非潮汐变化,结果表明琼中台的重力气压导纳值为-0.34×10-8m/s2/mbar,气压改正幅度约为10×10-8m/s2,海潮改正幅度约为5×10-8m/s2。改正后,琼中台重力非潮汐变化数据,比仅进行零漂固体潮改正后的重力非潮汐变化数据中的潮汐信号更加微弱,说明进行海潮改正后的效果是明显的,该方法可进一步去除其中的潮汐信号。
关键词连续重力    潮汐分析    零漂    重力气压导纳值    海潮改正    
The Air Pressure and Ocean Tide Correction of Gravity Non-tidal Changes in Qiongzhong Seismic Station
Li Sheng, Hu Jiuchang, Zhou Wen, Zhang Hui, Wang Huilin, Xie Xiaoling     
Hainan Earthquake Agency, Haikou 570203, China
Abstract: In this paper, we collected the continuous observed gravity data of Qiongzhong seismic station, and performed tidal and non-tidal analysis based on the processed data. Tide analysis was done by using VAV harmonic method and Non-tidal analysis was done by performing zero drift correction, solid tide correction, barometric pressure correction, and tide correction. Among them, the zero drift correction was done by the general polynomial fitting zero drift method, the air pressure correction was performed by VAV software, and the tide correction was done by using the SPOTL program to calculate the tidal load of Qiongzhong seismic by NAO.99b tide model. Finally, we obtained the non-tidal changes about the gravity of Qiongzhong seismic station. Our results show that the gravity pressure admittance of Qiongzhong Platform is -0.34×10-8m/s2/mbar, the pressure correction range is about 10×10-8m/s2, and the tide correction range is about 5×10-8m/s2. After the above correction, the gravity non-tidal change data of Qiongzhong seismic station is weaker than the tidal signal in the gravity non-tidal change data after only zero drift solid tide correction. It shows that the effect of tide correction is obvious, and the tide signal can be further removed.
Key words: Continuous gravity     Tide analysis     Zero drift     Gravity admittance value     Tide correction    
0 引言

重力随时间的变化可分为潮汐变化和非潮汐变化,前者如地球潮汐与地球自转变化,后者一般是地球物质迁移引起的重力随时间的变化(王谦身,2003)。连续重力观测是在台站固定点上进行的连续观测,获得的台站重力变化是台站高程以及台站和其周围地壳物性特征随时间变化的结果。连续重力观测能够精确探测到地球系统各圈层物质迁移引起的重力变化效应,包括叠加在一起的长期(多年的)或短期(季节性的)的变化(Sun et al,200120092011;周江林等,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 琼中地震台位置
1 原始资料分析

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

图 2 琼中台预处理后的重力潮汐变化数据(a)、拟合零漂(b)与零漂改正后残差(c)时序图
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波潮汐因子的变化。

表 1 琼中台连续重力各潮汐波参数
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.3 气压改正

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

表 2 不同机构的大气重力导纳值

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

表 3 琼中台2009年1月—2019年9月每年各频段大气重力导纳值(单位: ×10-8m/s2/mbar)

基于上述结果,取2009年1月至2019年9月的大气重力导纳值,即-0.34×10-8m/s2/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,20112013)。

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

图 4 琼中台重力潮汐变化零漂固体潮气压海潮改正时序图

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

表 4 琼中台重力潮汐变化、零漂固体潮改正、气压改正与海潮改正数据各主要潮波振幅(单位: ×10-8m/s2)
3.5 小结

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

图 5 琼中台重力气压负荷和海潮负荷改正值

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

图 6 琼中台海潮改正后的重力非潮汐年变序列
4 结论

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

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

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

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

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

参考文献
李辉、李瑞浩, 1994, 固体潮观测数据的滤波预处理方法, 地壳形变与地震, 14(1): 23-31.
周江林、沈萍、田鑫, 2015, 北京地震台gphone重力仪同震响应特征分析, 中国地震, 31(3): 553-561. DOI:10.3969/j.issn.1001-4683.2015.03.010
翟丽娜、孔祥瑞、杨锦玲等, 2020, 大连台重力固体潮受海潮负荷影响特征研究, 中国地震, 36(4): 953-960.
李盛、孙三健、张慧等, 2016, 热带气旋对琼中台重力非潮汐变化影响研究, 地震工程学报, 38(增刊Ⅰ): 101-105.
李盛、廖桂金, 2020, 雷琼地区重力场动态变化特征研究, 大地测量与地球动力学, 40(11): 1118-1125.
李英冰, 2003. 固体地球的环境变化响应. 博士学位论文. 武汉: 武汉大学.
罗少聪、孙和平, 2000, 大气压力变化对武汉台站重力场观测的影响, 测绘学报, 29(增刊Ⅰ): 75-79.
孙和平、许厚泽、周江存等, 2005, 武汉超导重力仪观测最新结果和海潮模型研究, 地球物理学报, 48(2): 299-307. DOI:10.3321/j.issn:0001-5733.2005.02.010
孙和平、周江存, 2002, 中国地壳运动观测网络基准站重力场变化的海潮负荷信号改正问题, 地球科学进展, 17(1): 39-43. DOI:10.3321/j.issn:1001-8166.2002.01.007
王迪晋、申文斌, 2009, 利用超导重力数据研究若干地球科学问题的进展, 测绘科学, 34(增刊Ⅰ): 10~13, 152.
王谦身, 2003, 重力学, 6-7, 北京: 地震出版社.
韦进、李辉、刘子维等, 2011, 利用积谱研究气压对SGC053超导重力仪的影响, 大地测量与地球动力学, 31(4): 47-51.
杨锦玲、关玉梅、钟继茂等, 2016, 厦门重力固体潮海潮负荷改正研究, 地球物理学进展, 31(3): 992-998.
周江存, 李辉, 孙和平, 等, 2004. 地震网络重力固体潮台站观测的海潮负荷影响. 见: 重力学与固体潮学术研讨会暨祝贺许厚泽院士70寿辰研讨会会议论文集. 武汉: 中国地球物理学会.
Agnew D C, 1997, NLOADF: a program for computing ocean-tide loading, J Geophys Res, 102(B3): 5109-5110. DOI:10.1029/96JB03458
Farrell W E, 1972, Deformation of the Earth by surface loads, Rev Geophys, 10(3): 761-797. DOI:10.1029/RG010i003p00761
Hu X G, Liu L T, Hinderer J, et al, 2006, Wavelet filter analysis of atmospheric pressure effects in the long-period seismic mode band, Phys Earth Planet Int, 154(1): 70-84. DOI:10.1016/j.pepi.2005.09.003
Kim T H, Shibuya K, 2013, Verification of the ellipsoidal earth model with an inelastic and convective mantle using tidal gravity variations revisited, Geophys J Int, 194(1): 230-248. DOI:10.1093/gji/ggt108
Kim T H, Shibuya K, Doi K, et al, 2011, Validation of global ocean tide models using the superconducting gravimeter data at Syowa Station, Antarctica, and in situ tide gauge and bottom-pressure observations, Polar Sci, 5(1): 21-39. DOI:10.1016/j.polar.2010.11.001
Sun H P, 1992, Comprehensive researches for the effect of the ocean loading on gravity observations in the Western Pacific Area, Marées Terrestres, 113: 8271-8292.
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), J Geodesy, 74(10): 720-729. DOI:10.1007/s001900000139
Sun W K, Wang Q, Li H, et al, 2009, Gravity and GPS measurements reveal mass loss beneath the Tibetan Plateau: geodetic evidence of increasing crustal thickness, Geophys Res Lett, 36(2): L02303.
Sun W K, Wang Q, Li H, et al, 2011, A reinvestigation of crustal thickness in the Tibetan Plateau using absolute gravity, GPS and GRACE Data, Terr Atmos Oceanic Sci, 22(2): 109-119. DOI:10.3319/TAO.2010.06.07.01(TibXS)
van Camp M, Vauterin P, 2005, Tsoft: graphical and interactive software for the analysis of time series and Earth tides, Comput Geosci, 31(5): 631-640. DOI:10.1016/j.cageo.2004.11.015
Venedikov A P, Arnoso J, Vieira R, 2003, VAV: a program for tidal data processing, Comput Geosci, 29(4): 487-502. DOI:10.1016/S0098-3004(03)00019-0
Tamura Y, Sato T, Ooe M, et al, 1991, A procedure for tidal analysis with a Bayesian information criterion, Geophys J Int, 104(3): 507-516.
Yu N H, 2006, Ocean tide models' performance in coastal regions, Annals of Shanghai Astronomical Observatory Chinese Academy of Sciences, 17(27): 17-25.