引言

随着全球地震观测台网的广泛布设，可用于研究的地震数据正在迅速增加，正确地利用这些数据可以极大地提高我们对很多地球科学问题的认识。然而，台站时间系统的故障和仪器布设错误等因素会造成观测数据出现的时钟错误和极性反转问题，对相关研究成果的可信度造成一定影响。波形数据中的时钟错误，会直接造成地震走时信息的拾取错误，影响地震波走时成像、震源定位等多方面研究（Tape et al，2009；Bondár，Storchak，2011）等。地震波形的极性是震源机制反演所需的关键数据，波形的初动极性对小震级地震的震源机制敏感，特别适用于余震序列研究（Tian et al，2013；Yukutake，Iio，2017），而极性反转问题的出现会直接影响到相关研究结果的可靠性。因此，从海量的地震数据中准确快速地检测和校正时钟错误和极性反转是地震数据预处理的重要步骤。

地震仪的内外时间系统不同步是造成时钟错误的主要原因。现代地震仪器的时间系统由内部时钟与外部同步装置组成（Schneider et al，1988），地震计内部的时钟存在线性、非线性时钟漂移，需要外部同步装置对它进行授时以校准时间偏差。二十世纪九十年代以来，大部分数字台站陆续开始使用全球定位系统（Global Positioning System，缩写为GPS）（Butler et al，2004）、全球导航卫星系统（Global Navigation Satellite System，缩写为GNSS）和网络进行授时与校时。但外部授时和校时系统有时会失效，例如：海底地震仪（ocean bottom seismometer，缩写为OBS）无法接收到GPS信号，只能使用内部时钟计时，从而导致数据存在时钟偏差，一般通过线性漂移公式对数据进行线性时钟漂移校正（Anchieta et al，2011）。此外，地震计硬件或软件故障也会造成数据较大的时间偏差（Gouédard et al，2014）。

地震台站波形数据的时钟核查是一件困难而耗时的工作（Sens-Schönfelder，2008）。目前主要使用两种方法对波形数据进行时钟检测和校正。Anchieta等（2011）利用海底地震仪与附近海岛上的陆地地震仪的P波到时差评估海底地震仪的时钟偏差。但是，这种方法需要研究区域有较高的地震活动性和清晰的震相记录，否则无法连续检测地震仪的时钟偏差。另一种方法是利用从台站间背景噪声提取的格林函数来检测时钟偏差，这种方法已经被用于陆地地震台站和OBS数据的时钟错误检查（Stehly et al，2007；Sens-Schönfelder，2008；王俊等，2013；Gouédard et al，2014；Abbas et al，2023）。但周期小于20 s的瑞雷面波在长距离传播后会大幅衰减，因此该方法只适用于检测中小孔径台阵（例如：OBS台阵、区域性台网或国家地震台网）的时钟偏差（Xia et al，2015）。针对早期全球地震台网（global seismographic network，缩写为GSN）布设稀疏的情况，Xia等（2015）利用在长距离传播衰减小的周期为26 s的噪声源实现了对洲际间超远台距的台站钟差检测，而对于南美洲的台站，由于周期为26 s的噪声源信号非常微弱，该方法的应用效果很差。全球面波成像研究需要使用数十年遍及全球的地震台站数据集（Ekström，2011；Ma，Masters，2014，Ma et al，2014；Priestley et al，2018），其中少数台站有时会出现较长时期的时钟错误，例如：台站GSC在1992年出现了持续数月的2 s时钟错误（Stehly et al，2007）、台站KOWA在2012年10—12月持续了约150 s的时钟错误（Xia et al，2015）。如果不对波形数据中的时钟错误进行检测和剔除，这些错误会严重影响相关地震学研究的结果。因此，针对长时期、全球范围的地震台站数据集，提出一种高效、准确的时钟错误检测方法，对地震学研究具有重要意义。

台站的地震计安装方向错误、电缆连接错误和数字化过程中的符号错误等因素都会导致波形的极性反转。近几十年来，地震台站数据极性反转的检测方法发展迅速。Hurst等（2002）和Roman等（2006）使用远震波形检测出混合台网中少量出现极性反转的台站。Ekström等（2006）通过对比理论模拟波形与实际观测的长周期地震图评估了大量永久性台站的波形数据质量，检测出少数发生极性反转的台站。Sens-Schönfelder （2008）利用从背景噪声中提取的格林函数的互相关异常检测出极性反转的波形记录。波形的极性反转现象已被多项研究发现（Ekström，Nettles，2010；Niu，Li，2011；Yang et al，2012；Young et al，2017；Petersen et al，2019）。Lentas （2021）提出了一种不依赖原始波形数据，仅使用国际地震中心（International Seismological Centre，缩写为ISC）公布的台站初动极性和相关震源模型的拟合优度概率检测波形数据的极性反转。

已有研究表明，虽然地震台站的时钟错误和极性反转问题的产生机制不同，但两者都会表现在波形记录的相位之中，进而导致测量得到的走时出现异常。基于这一共性，同时结合上述已有的成熟方法，本文拟提出一种利用地震面波相对走时测量值与P波到时相结合同时检测地震台站的时钟错误和极性反转的高效率方法，然后使用该方法检测全球2 010个地震台站出现时钟错误和波形反转的时段，以分析这种检测方法的可行性及检测结果的准确性。

1.   方法原理

全球M_W5.5以上地震平均每天发生1.2次（Pollitz et al，2012），地震激发的长周期面波信号可以在全球范围内被观测到（Ekström et al，2003）。这为实现全球台站以天为单位的时间尺度上的时钟错误和极性反转的连续检测提供了基础。

同一地震事件激发的面波，在不同台站所记录的相位差${\Delta \varphi }_{\mathrm{n}} $可以表示为

式中，$ \Delta {\varphi }_{\mathrm{s}} $为震源初始相位差，$ \Delta\varphi\mathrm{_{p\mathrm{\mathrm{ }}}} $为地震波从震源传播至台站产生的相位差，$ \Delta {\varphi }_{\mathrm{r}} $为地震台站的仪器响应不同而导致的相位差。在同一频率下，地震波传播时产生的相位和到时可以相互变换，即$ \Delta \varphi ＝ \omega \Delta t $，式中$ \Delta \varphi $，$ \omega $和$ \Delta t $分别为地震波相位变化量，角频率和时间变化量，因此可以先对波形进行窄带滤波，然后利用在时间域测量的台站相对走时换算得到台站的相位差（Jin，Gaherty，2015）。去除震源初始相位和仪器响应导致的相位差之后，不同台站之间的相对走时${t}_{\mathrm{r}} $可以表示为

式中：$ {\Delta t}_{2\mathrm{D}} $为面波从震源传播至不同台站所需时间差；$ {\sigma}_{t } $为测量误差，通过比较邻近台站的测量值，取为4 s （Ekström，2011；Ma et al，2014）；$ \Delta t $则为时钟错误导致的相对走时差，时钟正常台站的$ \Delta t $为0。为避免周波跳跃，先通过式（3）与式（4）去除根据已有相速度图预测的到时差。

式中，$ {{\tilde t}_{\mathrm{r}}} $为相对走时测量值，t_{2D_start}为已有相速度的预测到时，${\Delta{ \tilde t}_{2\mathrm{D}}}$为台站间实际到时差，$ c （ \theta \text{，} \phi ） $为某个频率下对应于经纬度$ （ \theta \text{，} \phi ） $的起始相速度（本文取10 mHz），$\varDelta $为震中距。去除已有相速度图预测的到时差后，绝大部分相对走时测量值$ {{\tilde t}_{\mathrm{r}}} $理论上应当处于一个较小的数值范围，因少部分地震波数据预测的$ {{\tilde t}_{2\mathrm{D}}} $受地球结构影响大，$ {{\tilde t}_{\mathrm{r}}} $超过理论范围，但该部分测量值可以在面波波速反演中得到较好的拟合，对应的拟合残差较小。

当台站在较长时间（数天乃至数月）内出现较大时钟偏差∆t （10 s以上，3.3节中有更详细的论述）时，相对走时测量值会连续出现较大异常，且其对应的拟合残差不会回落到较小的数值范围。当台站发生极性反转时，由互相关函数计算出的相对走时会包含半个周期的固定偏差，在反演中无法被拟合以致其残差仍为较大的异常值。因此，无论是时钟错误还是极性反转问题，仅在一个频率下对波形进行相位测量即可检测。而由于10 mHz频率波形的测量技术成熟，不易发生周期跳跃，且极性反转在周期为100 s的波形互相关测量中表现为半个周期50 s的偏差，有利于观察识别。因此，我们选择该频率波形进行检测。

然而，面波相对走时测量值会受一些因素的影响：① 面波相位测量之前需要利用震源位置参数去除震源初始相位$ {\varphi }_{\mathrm{s}} $，而震源水平定位存在的误差有时会达到20—50 km，使得相对走时的测量值出现偏差。事实上，震源定位的误差呈随机分布，由其产生的偏差也是随机分布。针对震源定位误差导致的偏差进行测试，结果显示，相较于时钟错误和极性反转表现出的连续且固定的偏差，震源定位导致的偏差明显不同（见附图1）。② 修正震源初始相位$ {\varphi }_{\mathrm{s}} $时利用了一维初始参考地球模型（preliminary reference earth model，缩写为PREM），因此需要考虑地下局部显著速度异常体的影响。此影响可以分为两种情况：一是台站正下方存在局部显著速度异常体异常，则该台站全部的相对走时测量值都会包含由其引起的固定偏差，该情况很容易从检测结果中排除；二是由于局部显著速度异常体不处于台站正下方，地震波的传播路径穿过该局部显著速度异常体是随机的，以致引起的偏差也随机出现。以上两种情况均无法对检测结果产生影响。③ M_W7.5以上地震的破裂尺度大且破裂过程复杂，会对相位校正造成影响（Ma，Masters，2015），但是本文所选用的数据中，2008—2012年间M_W7.5以上地震仅为26次，相对于M_W5.5以上地震的2 741次，数量非常少且不会集中在同一时期出现，因此对检测结果的影响可以忽略。

为了进一步排除利用面波相对走时检测可能导致的误判，同时区分时钟错误和极性反转，本文选择结合预测P波到时对初筛的异常数据进行二次检验。目前对全球P波波速的研究已具有相对高的精准度，远震P波到时相较于使用一维PREM模型计算的理论P波到时的偏差通常小于10 s （Xia et al，2015）。本文选用异常时段内具有清晰P波震相的全部地震事件并计算其理论P波到时，若台站出现时钟错误，其地震记录的P波震相较于预测到时会出现明显的提前或滞后；若台站记录发生极性反转，则地震记录P波震相与预测到时会基本一致。

2.   数据与应用

2.1   数据

根据美国地质调查局（United States Geological Survey，缩写为USGS）提供的地震目录，分别从美国地震学研究联合会（Incorporated Research Institutions for Seismology，缩写为IRIS）和德国波茨坦地学研究中心（Helmholtz Centre Potsdam-GFZ German Research Centre for Geosciences，缩写为GFZ）的GEOFON项目中收集了2008—2012年间M_W5.5以上地震的长周期面波波形数据，选定台站可供下载的全部长周期垂直分量（LHZ分量），涵盖全球30个台网的2 010个地震台站。目前从全球地震数据中心的两种主流数据下载方式分别为IRIS提供的BREQ_FAST和基于Python的ObsPy地震数据下载方式。因本研究使用的数据量大，对两种下载方式的效率进行测试。具体测试过程如下：选取台网II的6个台站2008年1月BHZ分量的mseed格式的波形数据；分别以BREQ_FAST和ObsPy两种方式在相同网络环境下交替进行三次下载，下载用时统计列于表1。结果显示基于ftp传输协议的BREQ_FAST实测下载效率比ObsPy高出162%。因此本文选用BREQ_FAST下载方式进行数据下载。实际应用的台网名称及其台站数量列于表2。

表  1  BREQ_FAST和ObsPy两种方式在相同网络环境中的下载用时统计

Table  1.  The download time of BREQ_FAST and ObsPy in the same network environment

下载方式	第一次用时/s	第二次用时/s	第三次用时/s	平均用时/s
BREQ_FAST	216.67	216.07	217.93	216.89
ObsPy	566.50	570.35	568.56	568.47

下载: 导出CSV 

| 显示表格

表  2  本文研究实际应用数据所属台网和该台网中包含的台站数量

Table  2.  The networks recorded the actual application data and the station number in networks

台网	台站数量		台网	台站数量		台网	台站数量		台网	台站数量		台网	台站数量		台网	台站数量
TA	1 439		CU	9		IU	79		NZ	10		G	32		XH	16
AF	1		CZ	2		KN	10		PM	5		GE	67		XN	11
9D	2		DK	16		MY	7		PS	5		IC	10		XV	2
BE	1		ER	1		US	62		SS	1		II	40		XY	58
C	12		MN	17		NR	13		TW	7		YT	45		ZM	30

下载: 导出CSV 

| 显示表格

2.2   测量相对走时

2.2.1   预处理

数据的预处理流程如下：① 利用傅里叶变换将原始波形从时间域变换到频率域，并采用中心频率为10 mHz的窄带通高斯滤波器对其进行滤波（Dziewonski et al，1969），在频率域去除仪器响应并进行归一化处理，再利用反傅里叶变换将其变换回时间域；② 修正震源初始相位、去除面波辐射花样及修正频散项相位，每个震源-台站对的震源初始相位$ {\varphi }_{\mathrm{s}} （ \omega ） $计算所使用的模型均为PREM，震源位置与震源深度参数使用全球矩心矩张量（Global Centroid-Moment-Tensor ，缩写为GCMT）公布的数据；③ 利用已有的全球10 mHz面波相速度图（Ekström，2011）计算理论走时，对同一地震事件的地震波数据进行初步时移。

2.2.2   测量

在波形测量过程中，为了避免小部分低信噪比记录干扰地震波形的相关性，首先选取信噪比大于4的数据（Houser et al，2008），信噪比定义为等长的信号段与非信号段上各点幅值平方和的比值。对不同台站记录的同一事件波形进行互相关计算，识别不同延迟时间对应的正峰值，之后利用VanDecar和Crosson （1990）提出的时域法计算同一地震事件不同台站的相对走时。

相对走时的测量使用聚类分析技术，即根据互相关系数衡量波形间的相似性，利用Hartigan （1975）的聚类算法生成类簇与聚类树。图1所示为不同台站记录的同一地震事件的波形聚类，其中左侧为利用相对走时进行时移的波形，右侧为相对应的聚类树。在聚类过程分析中，扭曲失真的波形记录和噪声记录可以被识别出并舍弃。利用聚类算法可以实现一次处理数百道地震波形数据，仅需要人工监督聚类时的分簇情况并决定保留或舍弃群簇，实现高效率高质量的相对走时测量。

图  1  一个波形聚类示例图

此示例中有827道波形，程序允许点击右侧的聚类树来选择合适的截止值。垂直黄线为此示例的合理截断值

Figure  1.  An example of waveform clustering

There are 827 waveforms in this example. The program allows to select the appropriate cut-off value by clicking on the cluster tree on the right side with the mouse. The vertical yellow line represents a reasonable cut-off for this example

下载: 全尺寸图片 幻灯片

2.3   面波相对走时筛选

本研究共进行了116万7 072道地震面波记录相对走时的测量，并统计分析了全部测量值。如图2所示，地震面波的整体测量平均值为−0.040 7，标准差为10.18。这表明大部分相对走时测量值介于−10—10 s，过小的时钟偏差无法识别。本文方法仅对大于10 s的时钟偏差能够进行有效检测，因此将−10 —10 s定义为正常范围。

图  2  相对走时测量值的分布直方图（灰色方块）与正态分布曲线（蓝色曲线）

Figure  2.  Distribution histogram （gray squares） of relative traveltime measurements and the normal distribution curve （blue curve）

下载: 全尺寸图片 幻灯片

使用本次测量得到的10 mHz面波数据集进行全球波速初步反演，反演的详细步骤参见文献（Ma et al，2014）。反演结果如图3所示，成像结果和拟合残差均表明本次测量的数据集整体质量良好。

图  3  利用10 mHz面波初步反演全球波速图及其拟合残差分布

图（a）−（c）分别为平滑因子$ \lambda $取2，10和100的结果；图（d）为$ \lambda $取10的相对走时测量值拟合残差直方图

Figure  3.  Global surface wave velocity map by 10 mHz surface wave inversion and the distribution of misfit residuals

Figs. （a）−（c） are the results when the smoothing factor λ is 2，10 and 100，respectively；Fig. （d） is histogram of the fitting residuals of relative traveltime measurements when λ is 10

下载: 全尺寸图片 幻灯片

对于单个台站，我们将相对走时测量值和拟合残差按其对应的地震事件发生时间由早到晚排列。如图4所示：时钟正常台站AAK的相对走时测量值绝大多数处于−10—10 s的正常范围内；由实际传播路径造成的较大相对走时测量值对应的拟合残差回落到正常范围内，如第50天的数据点所示；而受测量误差和震源定位误差等因素的干扰，个别相对走时测量值和拟合残差均异常大的情况随机地出现，如第110天的数据点。

图  4  2012年第1—170天台网II中台站AAK的相对走时测量值与拟合残差折线图

Figure  4.  Line chart of relative traveltime measurements （blue dots） and fitting residuals （green square dots） of the station AAK of Network II during the period from the 1st to the 170th day in 2012

下载: 全尺寸图片 幻灯片

然而，当台站在一段时期内发生时钟错误，则相对走时和拟合残差会在时钟错误期间连续超出正常范围。为了对比说明，本文选取了位于欧洲中部相距约1.4° （一个波长约4°）的两个台站进行对比。图5a与图5b分别展示了台网G的时钟错误台站ECH与台网GE的正常台站STU在2011年前80天的面波相对走时测量值与拟合残差的折线图。可见，台站ECH在2011年的第40—60天期间，相对走时测量值与拟合残差全部为远超10 s的异常值；对比台站STU在同一时期的相对走时测量值与拟合残差（图5b），可以看出台站ECH在此段时间发生了较大时钟偏差。为了证明检测结果是否准确可靠，本文将全部的钟错台站与距其1—2个波长内的正常台站进行了对比，具体内容见附表1和附图2，3，4。

图  5  台站ECH （a）和STU （b）在2011年第1—80天的相对走时测量值与拟合残差折线图

Figure  5.  Line chart of the relative traveltime measurements （blue dots） and fitting residuals （green square dots） of the stations ECH （a） and STU （b） from the 1st to the 80th day in 2011，respectively

下载: 全尺寸图片 幻灯片

图6显示出台网II的ERM台站在2011上半年的相对走时测量值和拟合残差均超出了正常范围，在50 s或−50 s附近。测量的10 Hmz面波数据周期为100 s，这种半个周期的固定偏差说明台站的波形记录发生了极性反转。下一步我们将利用预测P波到时给出更准确的判断依据。

图  6  台站ERM在2010年第1—180天期间的相对走时测量值与拟合残差折线图

Figure  6.  Line chart of the relative traveltime measurements （blue dots） and fitting residuals （green square dots） of the station ERM from the 1st to the 180th day in 2010

下载: 全尺寸图片 幻灯片

2.4   预测P波到时检验

利用面波相对走时测量值的异常，已经从大量地震数据初步检测出少量台站发生时钟错误和极性反转的时段。为了实现对异常时段的二次核查，本文从初步检测出的异常数据中挑选信噪比大于9且震中距在20°—140°之间的地震记录，以确保可以观察到清晰可辨的P波震相；然后将波形滤波为低频率波0.2 Hz，再根据PREM模型预测的理论到时将波形记录进行时移。时钟正常台站记录的P波震相会清晰地排列于预测P波到时线之上，而时钟错误台站记录的P波震相则会较理论P波到时线出现明显提前和滞后的情形。图7a绘出了由面波相对走时初步筛选出的一个时钟错误台站ULN在2010年第150—230天的相对走时和拟合残差，该台站时钟错误时段与时钟正常时段中的地震事件波形记录对比如图7b和图7c所示。将时钟正常时段的地震记录根据预测P波到时进行时移后，其P波震相与其它正常台站地震记录的P波震相对齐于预测P波到时线上。而时钟出错时段的地震记录时移后，其P波震相较预测P波到时线滞后约20 s。对于每个筛查出的异常时段，选取其中全部而非单个具有清晰P波震相的波形记录进行如上的检测，使最终的检测结果更加准确可靠。

图  7  台站ULN的相对走时和拟合残差折线图（a）及示例地震事件5个台站的波形记录（b，c）

（a） 2010年第150—230天台网IU的台站ULN的相对走时测量值与拟合残差折线图；（b） 2010年第167天M_W6.2地震的5个台站的波形记录；（c） 2010年第222天M_W7.3地震的5个台站的波形记录。红色波形为时钟错误台站记录，黑色波形为正常台站记录，蓝色虚线表示预测P波到时，下同

Figure  7.  Line chart of relative traveltime measurements and fitting residuals and waveforms of two example events

（a） Line chart of relative traveltime measurements and fitting residuals at the station ULN of network IU during the period from 150th to 230th day in 2010；（b） The waveforms of the M_W6.2 earthquake recorded by five stations on the 167th day in 2010；（c） The waveforms of the M_W7.3 earthquake recorded by the five stations on the 222nd day in 2010。The red waveform is the record of the ULN station （clock error），and the black waveform is the record of the normal station. The blue dotted line indicates the predicted arrival of P wave，the same below

下载: 全尺寸图片 幻灯片

此外，利用预测P波到时还可以很好地区分时钟错误与极性反转。如图8所示，由台网II的台站ERM在2012年第130—220天的地震记录与台网CU的台站GRGR在2011年第140—175天的地震记录测量得到的相对走时都为50 s左右。然而，在利用预测P波到时对波形时移后，台站GRGR的地震记录的P波震相无法与其它时钟正常地震记录的P波震相对齐于预测到时线上，存在近50 s的时间偏差；而台站ERM的地震记录的P波震相与其它时钟正常的地震记录的P波震相对齐于预测到时线。由此可知，台站GRGR的错误数据是因时钟错误所致，而台站ERM的错误数据则是由极性反转所致。

图  8  出现记录错误的台站的相对走时测量值、拟合残差折线图与示例地震事件波形图

图（a）和（b）分别为台站ERM与GRGR于2011年第130—230天期间相对走时测量值与拟合残差折线图；图（c）和（d）分别为台站ERM和GRGR与对比台站于2011年第175天记录的M_W7.3地震的波形图

Figure  8.  Line chart of relative traveltime measurements and fitting residuals and waveforms of the example event

Figs. （a） and （b） show the relative traveltime measurements （blue dots） and fitting residuals （green square dots） of the stations ERM and GRGR from the 130th to 230th day of 2011，respectively；Figs. （c） and （d） are the waveforms of the M_W7.3 earthquake recorded by five stations on the 175th day of 2011，respectively

下载: 全尺寸图片 幻灯片

3.   结果

表3和表4分别列出了台站时钟错误和波形极性反转出现时段的最终检测结果。结果表明，在检测的2 010个地震台站中，有12个台站的14个时段出现大于10 s的时钟错误，7个台站的8个时段发生极性反转。

表  3  本文检测到的台站出现时钟错误的时段

Table  3.  Periods of clock errors detected at the seismic stations in this study

台网	台站	时钟错误时段		台网	台站	时钟错误时段
MN	BNI	2 011年第97—2 09天		IU	KOWA	2 012年第261—350天
G	ECH	2 011年第42—59天		GE	MALT	2 009年第174—229天
GE	EIL	2 008年第2 02—252天		GE	PSZ	2 008年第199—2 07天
IU	FURI	2 009年第190—251天		G	SPB	2 010年第5—75天
GE	KMBO	2 012年第156—222天		IU	ULN	2 010年第160—2 00天
CU	GRGR	2 011年第140—177天 2 011年第286—301天		G	SCZ	2 010年第352天—2 011年第165天 2 011年第350天—2 012年第153天

下载: 导出CSV 

| 显示表格

表  4  本文检测到的台站出现极性反转的时段

Table  4.  Periods of polarity reversal at the seismic stations in this study

台网	台站	极性反转时段		台网	台站	极性反转时段
NZ	BKZ	2 008年第1天—2 012年第365天		II	ERM	2 008年第1天—2 012年第303天
NZ	KHZ	2 008年第1天—2 011年第54天		NZ	HIZ	2 008第1天—2 012年第365天
NZ	OUZ	2 008年第1天—2 010年第55天		NZ	ODZ	2 008年第1—190天
NZ	QRZ	2 008年第1—232天		NZ	QRZ	2 008年第290天—2 009年第22 0天

下载: 导出CSV 

| 显示表格

4.   讨论与结论

本研究中时钟错误时段的检测结果全部经过了P波到时的检验。台站时钟出错的时段短期的仅有数天，较长期的可达数月。本文提出的方法只能对偏差量较大（＞10 s）的时钟错误进行有效检测，即使对于使用大量数据的全球面波波速成像研究，这种程度的时钟偏差量也会在反演中产生不可忽略的影响。另一方面，在数据量使用较小的区域性地震成像研究中，此类较大时钟偏差可能会对结果产生严重影响。在预处理阶段去除时钟错误的数据可以避免反演结果中出现错误的成像结构。

本文方法利用了M_W≥5.5地震激发的全球范围内可观测的长周期信号，由于全球M_W≥5.5地震事件的发生频率约为1.2次/d，此方法可以给出以天数为分辨率的时段结果。同样是应用于超远台站间距的检测方法，Xia等（2015）利用几内亚湾26 s噪声源的互相关分析仅给出以月份为分辨率的检测结果。

尽管本文方法已经在互相关测量之前进行了初步时移以避免周波跳跃，但个别台站出现了数百秒的时钟偏差，仍导致了周波跳跃的出现。这种情况下，相对走时测量值无法反映出真实的时钟错误量。例如，台网IU的台站KOWA已被证实在2012年10—12月出现了约150 s的时钟错误（Xia et al，2015）。如图9a所示，该台站在2012年第270—350天的相对走时测量值与拟合残差在60 s上下波动，图9b显示钟错期间的P波到时与理论到时相差约为150 s。因此对于真实时钟错误量的评估，利用P波到时进行二次检验是十分必要的。

图  9  台站KOWA 2012年相对走时测量值、拟合残差折线图（a）与示例地震事件波形图（b）

（a） 台站KOWA 在2012年下半年的相对走时测量值与拟合残差折线图；（b） 5个台站记录的M_W7.8地震2012年第302天的波形记录

Figure  9.  Line chart of relative traveltime measurements and fitting residuals and waveforms of example event

（a） The relative traveltime measurement value （blue dots） and fitting residuals （green square dots） of the station ULN during the second half of 2012； （b） The waveforms of the M_W7.8 earthquake on the 302th day of 2012 recorded by five stations

下载: 全尺寸图片 幻灯片

对于波形极性反转的检测结果，我们使用远震初至P波波形进行检验。选取位于新西兰岛的两个正常台站SNZO和URZ和一个发生极性反转的台站BKZ。如图10所示，3个台站相对于2008年M_W7.9汶川地震震源方位角非常接近，台站BKZ出现明显的极性反转（本文所说的极性反转为垂直分量极性反转，并不涉及水平分量的方位角摆放错误）。出现极性反转的台站绝大部分归属台网NZ，我们就该问题联系了台网NZ的运营商GeoNet （Geological hazard monitoring system in New Zealand）。Geonet的回复证实了台网NZ部分台站的LHZ分量出现了极性反转的问题，并提供了一份已查明的发生极性反转的台站统计表单。与该表单的对比证明本文关于台网NZ发生极性反转台站的检测结果准确无误。对于垂直分量出现极性反转的原因，GeoNet认为相较于地震计的安装问题，布设时线缆连接错误的可能性更大。

图  10  2008年汶川M_S8.0地震事件及新西兰岛三个台站的位置分布（a）及其波形记录（b）（波形数据仅作0.01—0.02 Hz滤波处理）

Figure  10.  The 2008 Wenchuan M_S8.0 earthquake event and location of three stations at New Zealand （a） and corresponding waveform records （b） where waveform data were processed with a 0.01–0.02 Hz bandpass filter only

下载: 全尺寸图片 幻灯片

利用地震事件的波形数据，基于面波和P波技术检测时钟错误各有优劣。面波的能量较强，传播距离远，可以满足对全球范围内台站的检测，而由于面波测量时干扰因素的存在，难以保证结果的准确无误，而利用P波的检测技术虽然检测精度较高，但P波长距离传播后的低信噪比无法满足该方法对清晰P波震相的要求。故而本文将二者相结合，以达到更好的检测效果。然而，基于地震面波技术检测时钟偏差有其自身的局限性。受到测量误差、震源参数误差等因素的影响，该方法只能对较大的时钟偏差进行检测。此外，该方法适用于长时段的时钟偏差效应，对于单个事件的异常或者小于1天的钟错时段，无法做到有效检测，需结合其它资料进一步验证。针对目前方法的不足，未来考虑进一步改进。

本文提出利用面波相对走时测量值和P波到时相结合的方法实现对全球范围的地震台站同时检测出现时钟出错和极性反转的时段。该方法利用聚类分析进行面波相对走时的高效测量，进而实现了对大量地震台站数据准确快速地检测，且P波到时的二次检验保证了最终检测结果的准确性。之后使用该方法对全球30个台网的2 010个台站在2008—2012年期间记录到的M_W≥5.5地震的长周期面波波形数据进行了检测应用。结果表明，其中12个台站的14个时段出现了大于10 s的时钟错误，7个台站的8个时段出现了极性反转。对实际波形数据的应用表明本文的方法能够快速准确地检测出地震台站的时钟错误和极性反转时段，在提高全球面波地震成像质量的同时也为使用相关数据的震源机制反演等研究提供了参考。而该方法的不足之处在于，难以处理单个事件或小于1天的时钟错误和极性反转问题，无法对数秒或更小的时钟偏差进行有效检测。

附录

通过对去除初始相位时使用的震源位置参数添加50 km以内的随机水平误差，以定量地评估震源水平定位误差对面波相对走时测量的影响。此外，我们对震源深度参数添加±10 km随机误差以检测震源深度误差对面波相对走时测量的影响。具体结果如附图1所示，可见无论是水平定位误差还是深度误差，对面波相对走时的测量值的影响都很小，少量较大异常值随机出现，不会对最终结果造成影响。

  1  正常台站AAK与时钟错误台站ECH 的测量对比

（a） 正常测量；（b） 震源参数添加50 km内的随机水平定位误差；（c） 震源参数添加$ \pm 10 $ km内的随机深度定位误差

  1.  Measurement comparison between normal station AAK and clock error station ECH

（a） Normal measurements； （b） Random horizontal location errors within 50 km added to the source parameters；（c） Random depth location errors within ±10 km added to the source parameters

下载: 全尺寸图片 幻灯片

我们将台站间距与对应的相对走时测量值之差进行了分析（附图2），台间距小于10°的相对走时测量值之差大致在4 s之内。将全部的时钟错误台站与距其1$ — $2个波长之内台站的测量值与拟合残差进行对比。其中台站KOWA与SPB附近没有可用来对比的正常台站（附图3），但文中已经将KOWA与前人关于时钟错误的研究结果进行了对比分析，其余台站对列于附表1，对比如附图4所示。

  2  台站间距与相对走时测量值差的对应关系

误差棒取每个直条的比例中位数绝对偏差

  2.  Correspondence between station spacing and the difference in relative travel time measurements

The error bars take scaled median absolute deviation （SMAD） of each bin

下载: 全尺寸图片 幻灯片

  3  时钟错误台站（红色三角形）与附近正常台站（蓝色三角形）的分布

  3.  Distribution of clock error stations （red triangles） and nearby normal stations （blue triangles）

下载: 全尺寸图片 幻灯片

  1  对比台站对及其台间距

  1.  Comparison of station pairs and their spacing

对比台站	台站间距/°		对比台站	台站间距/°		对比台站	台站间距/°
BNI与SSB	1.52		FURI与DAMY	7.95		ECH与STU	1.46
GRGR与BBGH	2.28		SCZ与TPNV	4.15		ULN与TLY	4.41
EIL与KSDI	3.56		MALT与CSS	5.28		KMBO与MBAR	6.54
PSZ与MORC	2.42

下载: 导出CSV 

| 显示表格

  4  时钟错误台站的时钟错误时段（左）和附近同时段正常台站（右）的相对走时与拟合残差对比

  4.  Comparison of relative travel-time measurements with fitting residuals in the periods of clock errors at affected stations （left） and corresponding periods at nearby normal stations （right）

下载: 全尺寸图片 幻灯片

  4  时钟错误台站的时钟错误时段（左）和附近同时段正常台站（右）的相对走时与拟合残差对比

  4.  Comparison of relative travel-time measurements with fitting residuals in the periods of clock errors at affected stations （left） and corresponding periods at nearby normal stations （right）

下载: 全尺寸图片 幻灯片