Accuracy estimation of global tide models using continuous gravity observation data in northern China
-
摘要: 为评估全球潮汐模型在我国潮汐改正中的适用性，本文首先对10个重力站2016—2018年的观测数据进行了精度评定，而后基于均方根、和方根、纬度依赖关系以及重力残差等指标对7个全球潮汐模型进行了精度评定。结果表明：10个重力站的一些评价指标达到甚至超越了早期超导重力仪，例如M2波潮汐因子的中误差普遍小于0.000 70，其中最高精度约为0.000 14，5个主要潮波的稳定度均≤0.001 5。在10个观测模型和7个全球潮汐模型中，DDW-NHi和M2001模型考虑了地球扁率的影响，基于这两个模型计算的和方根较其它模型所得的和方根均小，约为0.288×10−8 m/s2。基于最高精度的乌什站数据对Molodensky， DDW-NHi，M2001与观测模型的改正精度的对比显示，DDW-NHi模型改正计算的重力残差（±0.4×10−8—±1.0×10−8 m/s2）不及观测模型（±0.1×10−8—±0.5×10−8 m/s2），但依然优于M2001模型（±0.7×10−8—±1.4×10−8 m/s2），且DDW-NHi模型改正获得的残差比传统的Molodensky模型所得残差（±0.5×10−8—±1.5×10−8 m/s2）小1×10−8—2×10−8 m/s2.
-
关键词:
- 全球潮汐模型 /
- 海潮负荷改正 /
- 重力残差 /
- 和方根 /
- 潮汐因子纬度依赖关系
Abstract: To evaluate the applicability of the global tide models for tide correction in China, this paper firstly assessed the data accuracy of ten gravity observation from 2016 to 2018, and then made the accuracy assessment for the seven global tide model by the evaluation indicators such as the root mean square （RMS）, the root sum square （RSS）, latitude dependence of the gravimetric amplitude factor, and residual gravity. The result shows that some of the evaluation indicators have achieved or even surpassed the accuracy of the early superconducting gravimeters in ten gravity observations. For example, the RMS of the M2 gravimetric amplitude factor was less than 0.000 70, withthe highest RMS about 0.000 14, and the stability of five main tide waves was less than 0.001 5. As for the ten observed and the seven global tide models, the RSS of DDW-NHi and M2001 models with the Earth’s oblateness influence is only about 0.288×10−8 m/s2, which is smaller than others. Comparison of the tide corrected accuracies by Molodensky, DDW-NHi, M2001, and the observed tide models suggests that residual gravity corrected by DDW-NHi （±0.4×10−8—±1.0×10−8 m/s2） is larger than that by observed one （±0.1×10−8—±0.5×10−8 m/s2） for Wushi gravity observatory with the highest tide accuracy, but still smaller than that by M2001 （±0.7×10−8—±1.4×10−8 m/s2）. Furthermore, the residual gravity corrected by DDW-NHi is less about 1×10−8—2×10−8 m/s2 than by the traditional Molodensky model. -
引言
随着我国地震观测台站数量的不断增加，区域地震监测能力的逐渐增强，可监测震级的下限不断降低。密集地震观测台网产出高质量地震观测数据的同时，也对现有业务系统提出了更多的新需求，其中从实时产出的地震观测数据中准确判别和拾取不同类型的震相是最基础的、最迫切的需求之一。这些震相信息为更好地认识地壳介质结构、了解发震断层、获取震源信息、进行地震预警服务等提供数据支持。
为了更快更准确地从连续地震记录波形中判别和拾取出地震震相，许多学者都开展了系列研究，并提出了很多实用的算法。这些震相拾取算法可归纳为能量准则（Allen，1982；Baer，Kradolfer，1987；Lomax et al，2012；马强等，2013）、极化分析（Vidale，1986；Ruud，Husebye，1992；Amoroso et al，2012）、模糊分析（Chu，Mendel，1994）、人工神经网络（Dai，MacBeth，1995）、高阶统计量（Saragiotis et al，2002；Küperkoch et al，2012）、小波变换（Anant，Dowla，1997；刘希强等，1998，2000）及基于机器学习（Kong et al，2019；Perol et al，2018；Ross et al，2018a，b；Zhu et al，2018；李安等，2020）等几大类相关算法。
相较于初动清晰、信噪比较高的P波震相，S波震相由于受到P波尾波和（或）其它转换波震相的干扰，信噪比通常较P波要低得多，因此S波震相拾取结果的准确性普遍低于P波震相。即使是有经验的分析人员，仍无法保证S波震相拾取的准确性。现有的S波震相自动拾取算法大都是基于极化特征分析的方法（Flinn，1965；Vidale，1986；Jurkevics，1988；Cichowicz，1993），这些算法利用三分向地震波形记录，基于P波震相与S波震相的不同偏振特征，通过计算质点运动的偏振度、线性度等参数，寻找特征函数突变发生的点，进而判定S波震相到时的位置。此类方法大都适用于单台记录，算法有效性依赖于其有效剥离多种震相的能力，如将剪切波与面波分离等。Lois等（2013）即利用一个基于特征值分析的特征函数，提出了一种S波自动拾取时域算法，该算法对于P波、S波的极化方向不作任何假设，是一种近乎“零参数”的时域实时算法。本文则在该算法的基础上，进一步优化相关参数，并利用福建地震台网2015—2018年期间的9 855条三分向波形记录分析进行验证。
1. 方法
在实时地震数据处理中，通常采用滑动窗的方式将连续观测数据离散为固定长度的数据，而后再采用相应离线数据处理算法进行分析。同样，对于一段已知P波到时位置的波形记录，本文在P波到时后采用滑动窗将记录划分，并根据一定的滑动步长进行滑动处理，然后对各小窗的三分向波形记录求取特征值，进而可获得特征值时程，本文将该特征值时程作为S波识别特征函数。其中，由三分向波形求取特征值可采用数学分析中协方差矩阵分解的方式得到，特征值时程的物理含义是确定质点运动偏振方向。假定由三分向波形数据计算得到的协方差矩阵为
$${\boldsymbol{C}} {\text{＝}} \left[ {\begin{array}{*{20}{c}} {{\rm{cov}} {\text{（}}x{\text{，}}x{\text{）}}}&{{\rm{cov}} {\text{（}}x{\text{，}}y{\text{）}}}&{{\rm{cov}} {\text{（}}x{\text{，}}z{\text{）}}} \\ {{\rm{cov}} {\text{（}}y{\text{，}}x{\text{）}}}&{{\rm{cov}} {\text{（}}y{\text{，}}y{\text{）}}}&{{\rm{cov}} {\text{（}}y{\text{，}}z{\text{）}}} \\ {{\rm{cov}} {\text{（}}z{\text{，}}x{\text{）}}}&{{\rm{cov}} {\text{（}}z{\text{，}}y{\text{）}}}&{{\rm{cov}} {\text{（}}z{\text{，}}z{\text{）}}} \end{array}} \right]{\text{，}}$$ (1) 式中，
${\rm{cov}} {\text{（}}x{\text{，}}y{\text{）}} = \Big[\sum\limits_{i {\text{＝}} 1}^L {{\text{（}}{x_i} {\text{－}} \overline x{\text{）}}{\text{（}}{y_i} {\text{－}} \overline y{\text{）}}}\Big]\Big/ L$ ，L为滑动窗长度，且${\rm{cov}} {\text{（}}x{\text{，}}y{\text{）}} {\text{＝}} {\rm{cov}} {\text{（}}y{\text{，}}x{\text{）}}$ ，${\rm{cov}} {\text{（}}x{\text{，}}z{\text{）}} {\text{＝}} $ $ {\rm{cov}} {\text{（}}z{\text{，}}x{\text{）}}$ ，${\rm{cov}} {\text{（}}y{\text{，}}z{\text{）}} {\text{＝}} {\rm{cov}} {\text{（}}z{\text{，}}y{\text{）}}$ 。将式（1）分解即可分别得到三个特征值${\lambda _1}{\text{，}}{\lambda _2}{\text{，}}{\lambda _3}$ （${\lambda _1} {\text{＞}} {\lambda _2} {\text{＞}} {\lambda _3}$ ）及相应的三个特征向量${ {\boldsymbol{u}}_{{1}}}{\text{，}}{ {\boldsymbol{u}}_{{2}}}{\text{，}}{ {\boldsymbol{u}}_{{3}}}$ 。对于连续波形记录，即可得到特征值时程${\lambda _1}{\text{（}}t{\text{）}}{\text{，}}{\lambda _2}{\text{（}}t{\text{）}}{\text{，}} $ $ {\lambda _3}{\text{（}}t{\text{）}}$ 。据此也可分别计算得到信号偏振线性度、质点运动平面度、P波传播方位角、线性运动入射角等参数（Jurkevics，1988；马强，2008）。由于
${\lambda _1}{\text{（}}t{\text{）}}$ 对信号传播方向上的能量突变更为敏感（图1d），因此本文仅选择其作为特征函数。为了更好地体现信号中能量的微弱变化并突显信号突变位置，本文求取${\lambda _1}{\text{（}}t{\text{）}}$ 的平方根$f{\text{（}}t{\text{）}} {\text{＝}} \sqrt {{\lambda _1}{\text{（}}t{\text{）}}}$ （图1e），放大${\lambda _1}{\text{（}}t{\text{）}}$ 中低值成分的同时压低${\lambda _1}{\text{（}}t{\text{）}}$ 中高值成分进而减小其动态范围，实现对特征函数的压缩。尤其对于信噪比较低的事件波形，采用此方式的特征函数凸显信号中的能量突变位置更为有效。随后，选取P波初至到S尾波段对$f{\text{（}}t{\text{）}}$ 求取峰度系数$K{\text{（}}t{\text{）}} {\text{＝}} $ $ kur[f{\text{（}}t{\text{）}}]$ ，即可从特征函数中将S波到时位置进一步凸显出来（图1f）。应用该时间窗是为了避免将后续事件的初至震相错误识别为本次事件S波震相，但对于密集事件序列，该时间窗的选取比较困难。由图1f可见，在S波到时附近$K{\text{（}}t{\text{）}}$ 呈现一个显著的陡坎，但S波到时点位于陡坎起始最低处，因此在实际应用中很难准确拾取，需要进一步处理，即对$K{\text{（}}t{\text{）}}$ 求取一阶差分$\Delta K{\text{（}}t{\text{）}}$ ，此时S波到时与$\Delta K{\text{（}}t{\text{）}}$ 峰值位置相一致，S波到时位置将更容易判别，据此即可对S波粗到时进行拾取。最后，应用自回归赤池信息准则（Akaike information criteria，缩写为AIC）（Akaike，1974）在所拾取S波粗到时位置附近精确判别（前后各取0.3 s），即可获得准确可靠的S波到时信息。![]()
图 1 本文S波拾取方法算例示意图事件发震时刻为2015-01-02 00:33:56.89，震中距为3.8 km的GTSK台站记录，窗长为0.2 s。图（a−c）为三分向波形记录；图（d）为特征值时程；图（e）为特征函数f （t）；图（f）为峰度系数时程K（t）；图（g）为求取差分后特征函数ΔK（t）Figure 1. An example of S phase picking by using the algorithm of this studyThe origin time of the event is 00:33:56.89 on 2 January 2015. The wavforms were recorded by the station GTSK with epicentral distance 3.8 km，and window length is taken as 0.2 s. Figs. （a） to （c） are three-component seismic records；Fig. （d） shows the history of three eigenvalues；Fig. （e） is the characteristic function f （t）；Fig. （f） shows time history of Kurtosis coefficient K（t），and Fig. （g） shows the differential characteristic function ΔK（t）滑动时窗长度是本文算法中需要确定的重要参数。若选择较短的时窗，该算法对于特征函数中微小的改变过于敏感，拾取到的S波到时位置将超前；若选择过长的时窗，算法将缺乏必要的敏感度，拾取到的S波到时位置将显著滞后。实时地震数据处理中，算法的鲁棒性是一个重要的指标，选择某一固定时间窗长度显然无法满足所有事件的需求。为此，本文采用多窗口综合加权方法，窗口长度为0.2—1.4 s，间隔0.2 s，共7个。首先在各窗口中拾取S波到时Si，并计算S震相信噪比Ri，然后依据信噪比综合加权给出S波到时位置，即
$${S_{\rm{f}}} {\text{＝}} \dfrac{{\displaystyle\sum\limits_{i {\text{＝}} 1}^7 {{S_i}{R_i}} }}{{\displaystyle\sum\limits_{i {\text{＝}} 1}^7 {{R_i}} }}{\text{.}}$$ (2) 2. 数据
本文研究中收集了福建台网2015—2018年记录的1 862次地方震事件，震级范围为ML−0.5—4.2，仅挑选其中震中距小于100 km的9 855条三分向地震动记录计算分析，在此震中距范围内，初至S波震相均为Sg。本文以人工地震编目提供的Sg震相拾取结果为参考，对比采用本文方法的拾取震相到时拾取位置，以验证本文算法的准确性及实用性。本研究中使用的台站观测记录以仙游震群和台湾海峡南部震群为主，这两个震群的观测数据量约占所有分析记录1/3左右，各台站记录数量空间分布如图2所示。本文所用台站记录在震中距100 km内基本呈均匀分布，震中距分布如图3所示。
![]()
图 2 本文所用台站记录随震中距和震级分布Figure 2. Epicentral distance and magnitude distribution for earthquake records used in this study3. 结果与讨论
应用本文方法得到的S波到时拾取结果与人工编目S震相拾取结果间的偏差统计情况如图4所示。结果显示，相较于人工拾取结果，本文方法的S波震相到时拾取平均偏差为（−0.003±1.34） s，其中拾取偏差小于0.2 s的记录所占比例为63.9%，拾取偏差小于0.5 s的记录所占比例为79.6%，即应用本文方法能够较准确地对大部分S波震相到时进行拾取。但从图3中也可以明显看出，应用本文方法的S波震相到时拾取结果有超前于人工拾取结果的趋势，因而造成图4b中统计直方图的偏态分布。本文将测试数据进一步划分为仙游震群、台湾海峡震群及其它事件三类，每类分别包括1 276条、2 330条和6 249条记录，并分析不同类别的S波到时拾取结果精度，结果如图5所示。由该图可见，不同的类别事件应用本文方法的S波拾取结果偏差分布较为一致，均表现出超前于人工拾取结果的趋势，三类事件的S波拾取偏差分别为（−0.018±1.68） s，（−0.053±1.56） s和（0.011±1.17） s。综上可见，本文所研究的相关算法具有较好的适用性，应用于不同区域的事件时均表现稳定。
![]()
图 4 应用本文方法的S波震相到时拾取偏差统计图（a）为各记录S震相拾取偏差散点图，图（b）和（c）为S震相拾取偏差统计直方图Figure 4. Statistics of S phase arrival time pick errors by using our methodFig. （a） gives S phase picking error for each record，Figs. （b） and （c） show the histogram of the phase picking error![]()
图 5 分区域S波震相到时拾取偏差结果统计（a） 仙游震群；（b） 台湾海峡南部震群；（c） 其它事件Figure 5. Statistics of S phase arrival time picking error in different regions（a） Xianyou sequence；（b） Southern Taiwan Strait sequence；（c） Other events此外，由图4还可以看出，仍有约11.2%的记录的S波震相到时拾取偏差超过±1.0 s，其中约4.1%的记录的S波震相到时拾取偏差超过±2.0 s，图6分别列举出三个拾取偏差较大的记录情况。进一步分析后，可将引起拾取偏差过大的原因归纳为三个方面。首先，信噪比较低是造成S波震相拾取偏差较大的显著原因。如图6a中所示波形记录，由于S波震相受到前序P波尾波的污染，信噪比较低，即使对于有经验的分析人员仍较难准确拾取，这也是通用S波震相到时拾取算法准确度不高的主要困难所在。本文采用式（3）截取人工拾取S波到时前后各2 s窗长内的波形记录并计算均方根进而得到各记录的信噪比，对于采样率为100 的记录：
![]()
图 6 S波震相到时拾取偏差较大的波形记录（a） 低信噪比记录结果；（b） 图（a）中记录1—20 Hz带通滤波后的结果；（c） 多事件叠加记录结果；（d） EW向异常记录结果Figure 6. Waveforms with large S phase arrival time picking error（a） The result for a low SNR record；（b） The result for the record in Fig. （a） after 1−20 Hz band-pass filtering；（c） The result for a multi-event record；（d） The result for an abnormal record in EW component$${\rm{SNR}} {\text{＝}} {\dfrac{{\sqrt {\dfrac{1}{n}\displaystyle\sum\limits_{i {\text{＝}} s}^{s {\text{＋}} n} {{x_i}} } } }{ {\sqrt {\dfrac{1}{n}\displaystyle\sum\limits_{i {\text{＝}} s {\text{－}} n}^s {{x_i}} } }}}{\text{，}}$$ (3) 式中s为人工拾取S波到时位置，n为数据长度（取200）.
随后，对406条S波到时拾取偏差超过±2.0 s的记录的信噪比进行统计，并与其它记录的信噪比情况进行对比，如图7所示。可以直观地看出低信噪比（＜2.0）记录在拾取偏差超过±2.0 s的记录中占有最大比例（图7a），与之相应，对于其余拾取误差较小的记录，高信噪比记录所占的比例较高（图7b）。图7以信噪比2.0为界，分别给出了低信噪比和高信噪比情况下S波震相拾取误差分布情况。由该结果可见：对于低信噪比波形记录，S波震相拾取误差超过1.0 s的比例为19%，显著多于高信噪比波形记录中的比例（9%）；且高信噪比波形记录中S波震相拾取误差小于0.2 s的比例（70%）更明显高于低信噪比波形记录中的比例（52%）。由此，作者认为观测记录信噪比是显著影响S波震相拾取精度的最重要因素。表1分别列出了针对不同信噪比记录，应用本文方法统计得到的S波震相拾取精度结果。可见，随着记录信噪比的提高，S波拾取残差呈平均值减小、标准差降低的趋势，与图7、图8所示结果一致。此外，从应用本文方法拾取的S－P到时差与人工拾取S－P到时差的对比（图9a）中也可以看出，偏差较大（如＞±2.0 s）记录的信噪比大多较低，而且随着震中距增加（S－P到时差增大），S波拾取结果的偏差分布也随之更为离散，显然这也与记录信噪比逐渐降低有关。同时，本文也对比了不同震群事件的S波拾取精度（图9b-d），本文算法在不同震群中的表现也较为一致，其中：仙游震群（图9b）中震中距较小的记录占比较高，算法对于此部分记录（＜5.0 s）也显示了更佳的性能，显然其信噪比也更高；同理海峡震群（图9c）和其它事件（图9d）中大震中距记录占比也相对更高，S波拾取残差随信噪比降低也变得更为离散，再次说明了记录的信噪比在应用本文S波拾取算法中的重要性。
![]()
图 7 地震记录信噪比的分布（a） 拾取偏差大于2.0 s的406条记录；（b） 拾取偏差小于2.0 s的9449条记录Figure 7. SNR distribution of seismic records（a） 406 records with S phase picking errors larger than ±2.0 s；（b） 9449 records with S phase picking errors less than ±2.0 s表 1 不同信噪比记录的S波到时拾取偏差统计Table 1. Statistic on S phase picking error for different SNR recordsSNR 记录数量 偏差均值/s 偏差中值/s 偏差标准差/s ＜2 2 953 0.135 0.040 1.92 2—5 6 416 −0.051 0.050 0.99 ≥5 486 −0.202 0.040 0.96 ![]()
图 8 不同信噪比记录S波震相拾取偏差分布（a） 信噪比小于2.0的记录；（b） 信噪比大于2.0的记录Figure 8. Pi chart of S phase detection error for different SNR records（a） Records with SNR less than 2.0；（b） Records with SNR larger than 2.0![]()
图 9 应用本文方法所得的S－P到时差与人工拾取结果对比及其信噪比分布（a） 所有记录；（b） 仙游震群；（c） 台湾海峡震群；（d） 其它事件Figure 9. Comparison of the arrival time differences for S−P by our method with those by manual picking and SNR distributions（a） All records；（b） Xianyou sequence；（c） Taiwan Strait sequence；（d） Other events若对事件波形进行带通滤波（如1—20 Hz），则可在一定程度上提高记录质量，也将有利于提高S波震相的拾取精度。图6b展示了图6a所示波形记录经1—20 Hz零相移带通滤波后再次应用本文方法的S波震相拾取结果。图10则分别为该记录滤波前后应用本文方法的具体识别效果，其中各子图所展示内容与图1一致。如前文已提及，不同滑动窗长下，S波到时拾取位置也会有所差异，因而本文采用式（2）对不同窗长下的S波到到时拾取位置进行加权综合，而图1和图10仅为某一窗长（分别为0.2 s和1.0 s）下的拾取结果，与图6中所示S波到时综合加权判别结果（虚线）会有所差异。对比可见，由于波形记录质量显著提高，S波震相拾取也更为准确，相较于人工拾取结果的偏差仅为0.36 s。需要说明的是，本文图4中所示震相拾取偏差均基于原始波形，未经带通滤波。鉴于上述分析，可推测出，若能够通过带通滤波提高部分低信噪比记录的信噪比水平，则震相拾取的偏差将得以改善。作者也将收集相关资料，继续深化该部分的研究。
![]()
图 10 图6a中信噪比较低记录带通滤波前（a）、后（b）的S波到时拾取各子图意思同图1，事件发震时间为2015-01-08 05:48:44.56，台站ZPCH的震中距为99.8 km，窗长为1.0 sFigure 10. S phase picking for a low SNR record in Fig. 6a before （a） and after （b） band-pass filteringEach subfigure has the same meaning as Fig. 1. The event occurred at 05:48:44.56 on 8 January 2015，which was recorded by the station ZPCH with epicentral distance 99.8 km. And time window length is taken as 1.0 s针对地震序列（尤其是大震后的密集余震序列）的处理能力也是考验震相拾取算法实用性的一个重要指标。鉴于序列震事件记录的复杂性及多样性，现有常规处理算法（如STA/LTA，AR-AIC，Z-detect等）均不能较为有效地处理和应对，本研究所涉及相关算法本质上也属于这类算法范畴。如前文所述，该算法不仅依赖于P波到时位置，而且仅能在P波到时之后的一段时间窗内识别唯一一个S波到时位置。显然，该算法基本不具备处理更为复杂的情况下序列地震的能力。图6c展示了一个密集序列的波形记录，在30 s内记录了多个小震事件，由于记录过于复杂，人工分析也仅标注了其中主要的S波震相到时（图6c中所标识Sg震相），而本文方法则只拾取到了序列中的一个可能的S波震相。由于本文方法尚不具备连续拾取后续震相的功能，因此本文方法在此类序列地震记录中的处理结果欠佳。在实际业务系统应用中，可适当提高S波尾波段阈值（即缩短事件持续时间），部分程度上可弥补在序列事件处理中的不足。
此外，由于波形记录本身存在的质量缺陷也是造成S波震相拾取偏差过大的显著影响因素。本文中所讨论相关算法的有效性和可靠性均依赖特征值分析结果，而一旦所分析的波形存在异常，则将直接影响特征值的分解结果，进而造成S波到时拾取结果的较大偏差。由图6d可见，由于EW向分量波形显著异常（该分向无输出），特征值计算错误，最终导致S波拾取结果出现较大偏差。尽管此类问题在实际业务应用中无法完全避免，但通过对记录波形质量的分析监控，能够有效地避免将异常波形纳入处理，从而提高S波到时拾取结果的准确性。
4. 结论
本文基于特征值分解算法，研究了一种可用于地方震S波震相实时提取的实用化算法。该算法对P、S波偏振方向不作任何假设，通过应用多个滑动时间窗即可较准确地拾取出S波震相到时，避免由于窗口长度选择不合理而造成的拾取偏差，满足地震信号的实时处理需求。应用福建地震台网地震观测的数据测试结果表明，该算法S波震相拾取精度高，适用于实时地震信息处理系统，但尚不具备密集序列连续处理能力。此外，由于应用该算法前需要首先已知P波到时信息，因此S波震相到时拾取的可靠性一定程度上也依赖于P波震相到时拾取的准确程度。但相较于S波震相，P波震相初动更为明显，因而更易于拾取，拾取精度更有保障。本文仅以福建地震台网100 km范围内的地方震事件为主，所拾取震相也主要为Sg震相，重点对该算法的适用性展开分析讨论，因而所得相关结论也可能不够全面。作者将继续收集相关波形数据，以更充分地论证该算法的适用性及其在实际数据处理系统中应用的可靠性。
随着地震观测台站数量的不断增加及观测数据质量的逐步提升，近年来深度学习等人工智能算法在地震数据处理领域内的应用也越来越多。尤其是在震相拾取方面，已经有多个学者提出了高性能、高可信的相关震相拾取算法，如GPD算法（Ross et al，2018a）、U-net算法（赵明等，2019）、PhaseNet算法（Zhu，Beroza，2018）等。作者也将应用这类算法开展相关研究，并将其与传统方法进行整合，以期更快、更准、更全地产出震相到时信息。
审稿专家为本文提出了建设性的修改意见，作者在此表示感谢。
-
![]()
图 1 我国10个重力观测站空间分布示意图
Figure 1. The spatial distribution of ten gravity observatories in China
![]()
图 2 2016—2018年主要潮波的潮汐因子稳定度变化
（a） O1波；（b） M2波；（c） P1波；（d） K1波
Figure 2. Variation of the stability of gravimetric amplitude factor of the main tidal waves from 2016 to 2018
（a） O1 wave；（b） M2 wave；（c） P1 wave；（d） K1 wave
![]()
图 3 O1波，M2波和K1波的潮汐振幅（a）和潮汐因子（b）的纬度依赖关系
Figure 3. The latitude dependence of the gravimetric amplitude （a） and gravimetric amplitude factors （b） of O1，M2 and K1 waves
![]()
图 4 潮汐模型潮汐改正后的重力残差时间序列及其功率谱密度
（a） 乌什重力站固体潮预处理数据；（b） Molodensky模型潮汐改正的重力残差；（c） DDW-NHi模型潮汐改正的重力残差；（d） M2001模型潮汐改正的重力残差；（e） 观测模型潮汐改正的重力残差；（f） 重力残差的功率谱密度
Figure 4. Residual gravity and power spectral density corrected by several tidal models
（a） The preprocessing gravity data from Wushi observatory；（b） The residual gravity corrected by Molodensky model；（c） The residual gravity corrected by DDW-NHi model；（d） The residual gravity corrected by M2001 model；（e） The residual gravity corrected by observed model；（f） The power spectral density of the residual gravity
![]()
图 5 周日波（a）和半日波（b）的重力残差功率谱密度频段特征
Figure 5. The power spectral density of the residual gravity for diural （a） and semi-diurnal （b） waves
![]()
图 6 海潮负荷改正前后观测的与模型的重力残差振幅比较
图中P代表进行了气压改正；＋N代表进行了Nao.99b海潮负荷改正；＋NC代表进行了Nao.99b和osu.chinasea.2010海潮负荷改正；＋NN代表进行了Nao.99b和naoregional.1999海潮负荷改正
Figure 6. Comparison of the residual gravity amplitude before and after corrected by ocean loading models and global body tide model
P represents the gravimetric amplitude factors corrected by atmospheric pressure；＋N represents those corrected by global ocean model of Nao.99b；＋NC represent the those corrected by the composite ocean model of Nao. 99b and osu.chinasea.2010；＋NN represent those corrected by the composite ocean model of Nao.99b and naoregional.1999
![]()
图 7 海潮改正对潮汐因子的影响
图中P代表进行了气压改正的潮汐因子；DUT10，EOT11a，FES2004，GOT4.7，HAM11a，Nao.99b和TPXO7.2分别代表进行了气压和对应全球海潮模型改正的潮汐因子；DDW-NHi代表了DDW-NHi理论潮汐因子。横坐标台站按经度由东向西排列
Figure 7. Impaction of the gravimetric amplitude factors corrected by seven global ocean models
P represents the gravimetric amplitude factors corrected by atmospheric pressure；DUT10，EOT11a，FES2004，GOT4.7，HAM11a，Nao.99b and TPXO7.2 represent those corrected by atmospheric pressure and corresponding global ocean tide model respectively；DDW-NHi represents those in DDW-NHi. The stations are arranged from east to west according to their longitude
![]()
图 8 格尔木观测站潮汐因子和调和分析中误差随格值系数的变化
（a） M2波；（b） O1波；（c） K1波；（d） P1波；（e） Q1波；（f） 中误差随格值系数的变化
Figure 8. Varaition of gravimetric amplitude factor and RMS of the harmonic analysis with the scale factor for Golmud observatory
（a） M2 wave；（b） O1 wave；（c） K1 wave；（d） P1 wave；（e） Q1 wave；（f） The RMS of the residual gravity
表 1 全球潮汐模型的特征
Table 1 Feature statistics of global body tide models
来源 长周期波 周日波 半日波 1/3日波 地球模型 Molodensky （1961） - 4 1 - 圆球 Warh （1981a，b，c） 1 19 1 0 1066A，PREM-C，C2 Dehant et al （1999） 1 9 1 1 流体静力弹性地球模型 11 9 1 1 非流体非弹性地球模型 Mathews （2001） - 16 - - PREM 
下载: 导出CSV
表 2 7个全球潮汐模型潮汐因子的常数项和纬度依赖项
Table 2 Constant term and latitude-dependent term of the gravimetric amplitude factor in seven global body tide models
全球潮汐模型 O1波 P1波 K1波 Q1波 M2波 ${G_0} $ 
$G_ { \pm }'/10^{-5}$ 
${G_0} $ 
$G_ { \pm }'/10^{-5}$ 
${G_0} $ 
$G_ { \pm }'/10^{-5}$ 
${G_0} $ 
$G_ { \pm }'/10^{-5}$ 
${G_0} $ 
$G_ { \pm }'/10^{-5}$ 
Molodensky 圆球 1.16 0 1.153 0 1.137 0 1.137 0 1.16 0 1066A 1.152 −600 1.147 −600 1.132 −600 1.152 −600 1.16 −500 Warh PREM-C 1.152 −700 1.147 −700 1.132 −600 1.152 −700 1.16 −500 C2 1.151 −700 1.147 −600 1.132 −600 1.151 −700 1.159 −500 Dehant DDW-He 1.154 24 8 1.147 77 −13 1.147 77 −13 1.132 83 −60 1.160 30 7 DDW-NHi 1.154 24 8 1.149 15 −10 1.134 89 −57 1.154 03 9 1.161 72 10 Mathews M2001 1.154 07 5 1.148 91 5 1.136 1 5 1.154 1 5 - - 注： ${G_0} $ 为潮汐因子的全球常数项，
$G_ { \pm }'$ 为纬度依赖项系数，其中M2001模型的
${G_0} $ 为原模型常数项的模，
$G_ { \pm }'$ 为原模型的周期项。
下载: 导出CSV
表 3 2016年重力潮汐观测站基本信息
Table 3 Basic information of the gravity observatories in 2016
序号 台站 省份 台站
类型观测环境
等级台站简介 中误差
/（10−8 m·s−2）1 高台 甘肃 地下室 优秀 位于黑河以北不足2 km处 ±0.637 2 格尔木 青海 山洞 优秀 位于昆仑布尔汉达山北麓山洞中 ±0.706 3 海拉尔 内蒙 山洞 优秀 位于海拉尔盆地与大兴安岭北坡海拉尔河与伊敏河交汇处 ±0.601 4 兰州 甘肃 山洞 优秀 位于兰州东盆地盐场堡以北、白塔山以东、黄河北岸Ⅲ级阶级地的后缘，
距黄河最近距离约1.5 km±0.740 5 牡丹江 黑龙江 山洞 良好 距市区约1.6 km，近年在距离台站500 m处建有成片住宅小区 ±0.718 6 沈阳 辽宁 山洞 良好 位于沈阳市东郊东陵区天柱山西南麓，南距沈抚公路北线约0.4 km ±0.747 7 泰安 山东 山洞 良好 位于泰山南麓正南，0.5 km处有小型水库 ±0.989 8 炭山 宁夏 山洞 良好 位于海原县炭山乡丘陵村，距海原县新区约20 km，台站供电系统时常故障 ±0.891 9 乌加河 内蒙 山洞 优秀 位于内蒙古自治区巴彦淖尔市乌拉特中旗乌加河镇正北1 km处 ±0.838 10 乌什 新疆 山洞 优秀 阿克苏地区乌什县西郊的七女坟旁 ±0.605
下载: 导出CSV
表 4 连续重力观测站的潮汐因子δ和气压导纳值的统计结果
Table 4 The result of gravimetric amplitude factor δ and barometric admittance of gravity observatories
重力站 O1波 P1波 K1波 M2波 S2波 气压导纳/（10−8 m·s−2·hPa−1） δ σ/10−4 δ σ/10−4 δ σ/10−4 δ σ/10−4 δ σ/10−4 数值 中误差 格尔木 1.152 94 4.0 1.145 22 14.3 1.131 84 9.9 1.165 22 1.5 1.165 10 3.4 −0.269 0.015 兰州 1.157 16 4.7 1.148 77 9.0 1.136 50 3.4 1.158 77 2.4 1.160 37 5.4 −0.319 0.011 高台 1.156 37 3.8 1.150 97 6.7 1.137 49 2.5 1.161 62 1.2 1.163 21 2.7 −0.391 0.008 海拉尔 1.155 19 5.9 1.157 16 11.2 1.138 65 4.1 1.167 74 4.2 1.171 67 9.4 −0.371 0.010 牡丹江 1.153 84 4.7 1.151 92 8.4 1.137 44 3.2 1.163 77 3.4 1.167 07 7.5 −0.347 0.009 沈阳 1.153 70 7.3 1.150 00 13.1 1.137 17 4.9 1.159 45 2.8 1.159 77 6.2 −0.356 0.010 泰安 1.158 32 11.3 1.153 09 20.6 1.139 41 7.7 1.155 73 2.4 1.156 20 5.6 −0.288 0.018 炭山 1.154 74 7.8 1.145 53 16.9 1.133 36 6.4 1.161 86 2.4 1.161 45 5.9 −0.312 0.024 乌加河 1.155 54 3.1 1.150 56 5.6 1.137 52 2.1 1.157 09 1.3 1.159 06 3.1 −0.299 0.010 乌什 1.155 14 3.9 1.151 40 7.0 1.137 26 2.6 1.160 97 1.4 1.161 00 3.2 −0.430 0.010 注：δ为潮波的潮汐因子，σ为中误差。
下载: 导出CSV
表 5 乌什重力观测站的潮汐分析结果
Table 5 The tidal analysis results of Wushi gravity observatory
起始频率
/cpd结束频率
/cpd潮波名称 潮汐振幅
/（10−9 m·s−2）潮汐因子 相位滞后/° 数值 中误差 数值 中误差 0.501 37 0.911 39 Q1 58.938 1 1.154 66 0.001 80 −0.097 0.09 0.911 39 0.947 99 O1 307.830 6 1.155 14 0.000 39 0.020 0.019 0.947 99 0.981 85 M1 24.209 7 1.504 93 0.004 20 −0.427 0.160 0.981 86 0.998 63 P1 143.232 9 1.151 40 0.000 70 −0.060 0.035 0.998 63 1.001 37 S1 3.387 2 1.209 77 0.042 58 16.412 2.319 1.001 37 1.004 11 K1 432.929 1 1.137 26 0.000 26 0.007 0.014 1.004 11 1.006 85 ψ1 3.387 1 1.260 70 0.030 07 −0.923 1.368 1.006 85 1.023 62 φ1 6.164 9 1.169 38 0.016 14 −0.622 0.791 1.023 62 1.057 49 J1 24.208 8 1.153 00 0.004 61 −0.174 0.229 1.057 49 1.470 24 OO1 13.246 6 1.162 69 0.011 38 0.117 0.561 1.470 24 1.880 26 2N2 13.001 1 1.161 73 0.003 27 0.201 0.161 1.880 27 1.914 13 N2 81.411 6 1.159 33 0.000 71 −0.046 0.035 1.914 13 1.950 42 M2 425.214 3 1.160 97 0.000 14 0.013 0.007 1.950 42 1.984 28 L2 12.018 8 1.157 81 0.006 06 −0.428 0.300 1.984 28 2.002 74 S2 197.832 1 1.161 00 0.000 32 −0.289 0.020 2.002 74 2.451 94 K2 53.780 9 1.158 76 0.001 54 −0.137 0.076 2.451 94 7.000 00 M3M6 6.287 4 1.073 14 0.006 61 0.109 0.353
下载: 导出CSV
表 6 各型重力仪潮汐观测精度和稳定度分析
Table 6 Accuracy and stability analysis of tidal observatory for different gravimeters
重力仪
类型数据时段 M2波潮汐因子
中误差中误差
/（10−8 m·s−2）导纳中误差
/（10−8 m·s−2·hPa−1）稳定度 来源 LRC 1973—1990 0.002 00—0.005 00 - - - Melchior和de Becher （1983） GEO 1971—1990 0.001 00—0.002 00 - - - GWR 1975—1988 0.000 08 - - - LRC 1983—1986 ≤0.005 00 - - 0.005 0 毛慧琴等（1989） GS 1989—1993 0.000 30—0.003 00 - - - 魏望生和喻节林（1995） DZW 0.000 80—0.004 00 - - - LRC 0.000 40—0.002 00 - - - GEO 0.001 00—0.005 00 - - - T，CT，CD 1989—2000 - ±0.05—±0.70 ±0.000 4—0.003 0 0.001 5 Ducarme et al （2002） ASK - - ±0.010 0 - LRC 2000—2001 0.000 51 ±0.20 ±0.022 0 - 孙和平等（2002） OSG 0.000 06 ±0.02 ±0.002 0 - T，CT，CD 1997—2001 0.000 70 - ±0.000 4—0.003 0 0.001 0 Xu et al （2004） OSG 2009—2010 0.000 04 ±1.10—±1.40 ±0.000 5 - 韦进等（2012）
下载: 导出CSV
表 7 7个全球潮汐模型计算所得重力残差的均方根RMS及和方根RSS
Table 7 RMS and RSS of the residual gravity corrected by seven global body tide models
潮汐模型 地球模型 RMS/（10−8 m·s−2） RSS/（10−8 m·s−2） M2波 K1波 O1波 P1波 Q1波 Tamura
潮波表郗钦文
潮波表Molodensky 均质圆球 0.143 9 0.168 7 0.173 4 0.064 8 0.047 8 0.292 8 0.291 2 Warh 1066A 0.190 9 0.179 8 0.169 3 0.068 9 0.039 9 0.322 1 0.322 0 PREM-C 0.190 9 0.179 8 0.169 4 0.069 3 0.040 0 0.322 3 0.320 6 C2 0.199 9 0.179 8 0.171 2 0.068 9 0.039 9 0.328 5 0.329 0 DDW DDW-He 0.143 6 0.227 5 0.167 1 0.066 6 0.053 7 0.328 0 0.328 1 DDW-NHi 0.144 5 0.171 3 0.167 1 0.065 2 0.040 2 0.289 8 0.289 6 Mathews M2001 0.144 5 0.168 3 0.167 2 0.065 2 0.040 2 0.288 1 0.287 9
下载: 导出CSV
表 8 两个全球潮汐模型与观测站的潮汐因子的相对误差
Table 8 The relative errors of the gravimetric amplitude factors between DDW-NHi and M2001 models and the observatories
台站 M2波潮汐因子相对误差 O1波潮汐因子相对误差 P1波潮汐因子相对误差 K1波潮汐因子相对误差 RMS/（10−8 m·s−2） DDW-NHi M2001 DDW-NHi M2001 DDW-NHi M2001 DDW-NHi M2001 DDW-NHi M2001 格尔木 0.251 3% 0.251 3% 0.178 2% 0.150 4% 0.253 3% 0.247 4% 0.259 5% 0.459 4% 0.120 6 0.141 0 高台 0.015 2% 0.015 2% 0.224 4% 0.253 1% 0.187 8% 0.189 1% 0.312 9% 0.094 1% 0.087 5 0.068 3 海拉尔 0.015 3% 0.495 8% 0.062 1% 0.093 9% 0.815 2% 0.800 9% 0.574 1% 0.290 9% 0.223 8 0.198 0 兰州 0.015 4% 0.211 4% 0.354 1% 0.381 8% 0.269 4% 0.275 8% 0.291 5% 0.093 8% 0.114 6 0.105 1 牡丹江 0.015 5% 0.200 8% 0.188 0% 0.157 6% 0.582 2% 0.575 4% 0.244 1% 0.008 6% 0.227 8 0.220 6 乌什 0.015 6% 0.236 3% 0.025 6% 0.003 9% 0.168 3% 0.165 9% 0.246 3% 0.011 9% 0.107 8 0.095 5
下载: 导出CSV
-
李大炜,李建成,金涛勇,胡敏章. 2012. 利用验潮站资料评估全球海潮模型的精度[J]. 大地测量与地球动力学,32(4):106–110. Li D W,Li J C,Jin T Y,Hu M Z. 2012. Accuracy estimation of recent global ocean tide models using tide gauge data[J]. Journal of Geodesy and Geodynamics,32(4):106–110 (in Chinese).
刘子维,李辉,韦进,郝洪涛,吴云龙. 2011. 利用M2潮波振幅因子精密测定gPhone弹簧重力仪的标定因子[J]. 大地测量与地球动力学,31(5):146–150. Liu Z W,Li H,Wei J,Hao H T,Wu Y L. 2011. Accurate determination of calibration factor of gPhone spring gravimeters by using M2 tidal wave amplitude factor[J]. Journal of Geodesy and Geodynamics,31(5):146–150 (in Chinese).
毛慧琴,许厚泽,宋兴黎,陈振邦. 1989. 中国东西重力潮汐剖面[J]. 地球物理学报,32(1):62–69. doi: 10.3321/j.issn:0001-5733.1989.01.006 Mao H Q,Hus H T,Song X L,Chen Z B. 1989. East-west gravity tidal profile of China[J]. Acta Geophysica Sinica,32(1):62–69 (in Chinese).
孙和平,许厚泽,徐建桥,柳林涛. 2000. 重力场的潮汐变化观测及其研究[J]. 地球科学进展,15(1):53–57. doi: 10.3321/j.issn:1001-8166.2000.01.008 Sun H P,Xu H T,Xu J Q,Liu L T. 2000. Observation and study of the tidal variation in gravity field[J]. Advance in Earth Sciences,15(1):53–57 (in Chinese).
孙和平,陈晓东,刘明,周百力. 2002. LCR-ET20弹簧型潮汐重力仪观测结果的分析和比较[J]. 地震学报,24(5):510–515. doi: 10.3321/j.issn:0253-3782.2002.05.008 Sun H P,Chen X D,Liu M,Zhou B L. 2002. Analysis and comparison of the tidal gravity observations obtained with LCR-ET20 spring gravimeter[J]. Acta Seismologica Sinica,24(5):510–515 (in Chinese).
韦进,李辉,刘子维,康开轩,郝洪涛. 2012. 武汉九峰地震台超导重力仪观测分析研究[J]. 地球物理学报,55(6):1894–1902. doi: 10.6038/j.issn.0001-5733.2012.06.010 Wei J,Li H,Liu Z W,Kang K X,Hao H T. 2012. Observation of superconducting gravimeter at Jiufeng seismic station[J]. Chinese Journal of Geophysics,55(6):1894–1902 (in Chinese).
魏望生,喻节林. 1995. 中国大陆重力潮汐参数及其特征研究[J]. 地壳形变与地震,15(4):77–85. Wei W S,Yu J L. 1995. Study on the gravity tide parameters and their characteristics in the mainland of China[J]. Crustal Deformation and Earthquake,15(4):77–85 (in Chinese).
郗钦文,侯天航. 1987. 新的引潮位完全展开[J]. 地球物理学报,30(4):349–362. Xi Q W,Hou T H. 1987. A new complete development of the tide-generating potential for the epoch J2000.0[J]. Acta Geophysica Sinica,30(4):349–362 (in Chinese).
张锐,韦进,刘子维,李辉,郝洪涛. 2011. 用SGC053超导重力仪观测资料对gPh058重力仪格值的精密测定[J]. 大地测量与地球动力学,31(5):151–155. Zhang R,Wei J,Liu Z W,Li H,Hao H T. 2011. Accurate determination of scale value of gPh058 gravimeter by use of observations with SGC053 superconducting gravimeter[J]. Journal of Geodesy and Geodynamics,31(5):151–155 (in Chinese).
周江存,徐建桥,孙和平. 2009. 中国大陆精密重力潮汐改正模型[J]. 地球物理学报,52(6):1474–1482. doi: 10.3969/j.issn.0001-5733.2009.06.008 Zhou J C,Xu J Q,Sun H P. 2009. Accurate correction models for tidal gravity in Chinese continent[J]. Chinese Journal of Geophysics,52(6):1474–1482 (in Chinese).
Dehant V,Ducarme B. 1987. Comparison between the theoretical and observed tidal gravimetric factors[J]. Phys Earth Planet Inter,49(3/4):192–212.
Dehant V,Defraigne P,Wahr J M. 1999. Tides for a convective Earth[J]. J Geophys Res,104(B1):1035–1058. doi: 10.1029/1998JB900051
Ducarme B,Sun H P,Xu J Q. 2002. New investigation of tidal gravity results from the GGP network[J]. Bull Inf Marees Terrestres,136:10761–10776.
Francis O,Niebauer T M,Sasagawa G,Klopping F,Gschwind J. 1998. Calibration of a superconducting gravimeter by comparison with an absolute gravimeter FG5 in Boulder[J]. Geophys Res Lett,25(7):1075–1078. doi: 10.1029/98GL00712
Hábel B,Meurers B. 2014. A new tidal analysis of superconducting gravity observations in western and central Europe[J]. Contribut Geophys Geodesy,44(1):1–24. doi: 10.2478/congeo-2014-0001
Love A E H. 1909. The yielding of the Earth to disturbing forces[J]. Mon Not R astr Soc,69(551):476.
Mathews P M. 2001. Love numbers and gravimetric factor for diurnal tides[J]. J Geodetic Soc Jpn,46(4):231–236.
Melchior P,de Becker M. 1983. A discussion of world-wide measurements of tidal gravity with respect to oceanic interactions,lithosphere heterogeneities,Earth’s flattening and inertial forces[J]. Phys Earth Planet Inter,31(1):27–53. doi: 10.1016/0031-9201(83)90064-X
Melchior P. 1994. A new data bank for tidal gravity measurements (DB 92)[J]. Phys Earth Planet Inter,82(2):125–155. doi: 10.1016/0031-9201(94)90085-X
Molodensky M S. 1961. The theory of nutation and diurnal Earth tides[J]. Communs Obs Roy Belg,288:25–56.
Tamura Y. 1987. A harmonic development of the tide-generating potential[J]. Bull Inf Marees Terrestres,99:6817–6855.
Wahr J M. 1981a. Body tides on an elliptical,rotating,elastic and oceanless earth[J]. Geophys J R astr Soc,64(3):677–703. doi: 10.1111/j.1365-246X.1981.tb02690.x
Wahr J M. 1981b. A normal mode expansion for the forced response of a rotating earth[J]. Geophys J R astr Soc,64(3):651–675. doi: 10.1111/j.1365-246X.1981.tb02689.x
Wahr J M. 1981c. The forced nutations of an elliptical,rotating,elastic and oceanless earth[J]. Geophys J R astr Soc,64(3):705–727. doi: 10.1111/j.1365-246X.1981.tb02691.x
Xu J Q,Sun H P,Ducarme B. 2004. A global experimental model for gravity tides of the Earth[J]. J Geodyn,38(3/5):293–306.
-
期刊类型引用(6)
1. 刘文玉，程正璞，年秀清，陈闫，胡钰铃，覃祖建，邵明正. 基于三维剩余密度结构的松原地震成因. 地震地质. 2024(02): 462-476 . 
百度学术
2. 张仁鹏，雷焕珍，陈聪，陈琳荣，冯靖乔. 2019年松原M_S5.1地震烈度评定及震害特征简析. 防灾减灾学报. 2023(01): 34-40 . 百度学术
3. 于晨，卢军，解滔，刘长生. 2017—2019年松原地区4次中强地震前绥化台地电阻率异常分析. 中国地震. 2023(02): 314-324 . 百度学术
4. 张仁鹏，冯靖乔，雷焕珍，于宏伟，马飞，王佳蕾. 松原市5.7级和5.1级地震房屋震害特征分析与启示. 防灾科技学院学报. 2022(03): 80-87 . 百度学术
5. 李梦莹，张志宏，焦明若，方禹心. 2018-05-28松原M_S5.7地震地电场变化特征研究. 大地测量与地球动力学. 2022(12): 1276-1280+1299 . 百度学术
6. 李梦莹，杨士超，张志宏，孔祥瑞，孙庆山. 2018年松原M_S 5.7地震地电场变化特征. 地震地磁观测与研究. 2022(S1): 120-122 . 百度学术
其他类型引用(2)
计量
- 文章访问数: 947
- HTML全文浏览量: 420
- PDF下载量: 74
- 被引次数: 8
