Detection of clock error and polarity reversal based on surface wave phase measurement and P-wave arrival time
-
摘要:
地震观测中会遇到台站出现时钟错误和波形极性反转等问题,需要在地震数据处理过程中对其进行检测和校正。基于此,本文提出了一种基于地震面波相位测量和P波到时检测台站时钟错误和极性反转的方法,并使用该方法检测了全球30个台网的2 010个台站2008—2012年间MW≥5.5的长周期地震波形数据,对使用本文方法处理结果显示,有12个台站的14个时段出现大于10 s的时钟错误;7个台站的8个时段发生了波形极性反转。表明此方法可以快速准确地检测出地震台站出的现较大时钟错误和极性反转时段,在提高全球面波地震成像质量的同时也为使用相关数据的震源机制反演等研究提供参考。但本文方法难以处理单个事件或时间长度小于1天的时钟错误和极性反转问题,无法对数秒或更小的时钟偏差进行有效检测。
Abstract:With the extensive deployment of seismic observatory networks around the world, the amount of seismic data available for research is rapidly increasing, and the proper utilization of these data can greatly improve our understanding of the dynamics of Earth. However, factors such as failures in station time systems and deployment errors can cause observational data containing clock error and polarity reversal, thereby reducing credibility of the relevant research results. Clock error in waveform data can directly cause time errors in the picking of major phases, which can affect results of source location studies, seismic travel-time tomography etc. Therefore, accurate, fast detection and subsequent correction of clock error and polarity reversal of terabyte sized seismic data, made possible by recent deployments of national and international networks, are important steps in seismic data preprocessing.
The unsynchronized internal and external time systems of station seismometers lead to clock error. Currently, two main methods are used for clock detection and correction of waveform data. The first method uses the difference in P-wave arrival-times of nearby stations to evaluate the clock error of seismometers. However, this method requires high seismicity and clear phase records in the study area, otherwise continuous detection of seismograph clock error cannot be achieved. Another more common method is to use the Green’s function extracted from the inter-station ambient noise to detect the clock error. However, Rayleigh surface waves with a period of less than 20 s attenuate greatly after long-distance propagation, so this method is only suitable for detecting clock error in small- and medium-aperture arrays. Although longer period noise signal can enable the detection of station clock error for ultra-long intercontinental station distances, the study of global surface wave velocity tomography that requires decades of seismic station datasets spread across the globe, both of these detection methods are not applicable. Therefore, it is important to propose new methods for efficient and accurate clock error detection for a decades and worldwide seismic station dataset.
Factors such as, incorrect orientation of the seismometers installed at the station, incorrect cable connections, and symbol errors during digitization, can lead to polarity reversal of the waveform. Like clock error, detection and correction of polarity reversal can increase accuracy of results e.g. in seismic source mechanism inversion.
Although the clock error and polarity reversal problems of seismic stations have different mechanisms, both of them are manifested in the phase of the waveform recordings, which results in anomalies in the measured travel-times. Therefore, it is feasible and meaningful to propose a method that can detect clock error and polarity reversal simultaneously.
The method in this paper utilizes the combination of seismic surface wave relative travel-time measurements and P-wave arrivals to achieve an efficient method of simultaneously detecting clock error and polarity reversal at stations. Firstly, globally observable long-period surface waves excited by events with MW≥5.5 are selected. The relative travel-times are obtained by measuring the same seismic event recordings from different stations. The anomalies of the relative travel-time measurements at a certain time period of one station can reflect the clock error or polarity reversal problems of the station. In order to further exclude the possible misjudgment caused in the detection using surface wave relative travel-time measurements and to distinguish the clock error and polarity reversal, the proposed method combines the prediction of P-wave arrivals as a secondary check on the anomalous data.
The proposed method was applied to detect long-period surface waveform data from 2010 stations of 30 networks during 2008−2012. The results show that 12 stations have clock error greater than 10 s in 14 periods, and 7 stations occur polarity reversal in 8 periods.
Seismic wave data clock error detection based on surface wave or P-wave techniques has advantages and disadvantages. Surface waves have stronger energy and usually propagate over long distances, which can satisfy the detection of stations on a global scale. Due to the existence of interference factors during the measurement of surface waves, it is difficult to ensure absolute accuracy of the results. Although the detection accuracy of P-wave method is higher than the surface wave method, the lower signal-to-noise ratio of P-wave propagation over long distances could not satisfy the requirements of the method for a clear P-wave phase. The proposed method combines the two techniques to achieve better detection results.
The proposed method utilizes the cluster analysis for the efficient measurement of the relative travel-times of surface waves and it contains two checks. The first check by surface waves realizes the accurate and fast detection of a large number of data, and the secondary check by the P-wave arrival time ensures the accuracy of the final detection results. The application of the actual waveform data proves that the method can quickly and accurately detect the clock error and polarity reversal of the stations, which improves the quality of global surface wave velocity tomography and also provides a reference for the research on the inversion of source mechanism using relevant data. However, the proposed method based on seismic surface wave technology has its own limitations. Due to measurement errors and source parameter errors, the method can only detect large clock error. Additionally, the method is applicable to the clock error over long time and cannot effectively detect the anomalies of single events or clock error of less than 1 day, which need to be further verified by combining with other data. Considering the limitations of the current method, further improvements may be worth considering in the future.
-
Keywords:
- seismic station /
- clock error /
- polarity reversal /
- cross-correlation analysis /
- P wave arrival time
-
1. 引言
随着全球地震观测台网的广泛布设,可用于研究的地震数据正在迅速增加,正确地利用这些数据可以极大地提高我们对很多地球科学问题的认识。然而,台站时间系统的故障和仪器布设错误等因素会造成观测数据中存在时钟错误和极性反转问题,从而降低了相关研究成果的可信度。波形数据中的时钟错误,会直接造成地震走时信息的拾取错误,影响地震学多方面的研究,例如:地震波走时成像、震源定位研究(Tape et al,2009;Bondár,Storchak,2011)等。地震波形的极性是震源机制反演所需的数据,地震波形的极性是震源机制反演所需的关键数据,波形的初动极性对小震级地震的震源机制敏感,特别适用于研究余震序列(Tian et al,2013;Yukutake,Iio,2017),而极性反转问题的出现会直接影响到相关研究的结果。因此,从海量的地震数据中准确快速地检测和校正时钟错误和极性反转是地震数据预处理的重要步骤。
地震仪的内外时间系统不同步是造成时钟错误的主要原因。现代地震仪器的时间系统由内部时钟与外部同步装置组成(Schneider et al,1988),安置于地震计内部的时钟存在线性和非线性的时钟漂移,需要外部同步装置对它进行授时以校准时间偏差。20世纪90年代以来,大部分数字台站开始逐步使用全球定位系统(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个地震台站出现时钟错误和波形反转的时段,以期分析了这种检测方法的可行性与检测结果的准确性。
2. 方法原理
全球MW5.5以上地震平均每天发生1.2次(Pollitz et al,2012),地震激发的长周期面波信号可以在全球范围内被观测到(Ekström et al,2003)。这为实现对全球台站以天为单位的时间尺度上连续检测时钟错误和极性反转提供了基础。
同一地震事件激发的面波,在不同台站所记录的相位差${\Delta \varphi }_{\mathrm{n}} $可以表示为
$$ {\Delta \varphi }_{\mathrm{n}}=\Delta {\varphi }_{\mathrm{s}} + \Delta {\varphi }_{\mathrm{l}} + \Delta {\varphi }_{\mathrm{r}} \text{,} $$ (1) 式中,$ \Delta {\varphi }_{\mathrm{s}} $为震源初始相位差,$ \Delta {\varphi }_{\mathrm{l}} $为地震波从震源传播至台站产生的相位差,$ \Delta {\varphi }_{\mathrm{r}} $为地震台站的仪器响应不同而导致的相位差。在同一频率下,地震波传播时产生的相位和到时可以相互转化($ \Delta \varphi = \omega \Delta t $,其中$ \Delta \varphi $,$ \omega $和$ \Delta t $分别为地震波相位变化量,角频率和时间变化量),因此可以先对波形进行窄带滤波,然后利用在时间域测量的台站相对走时换算得到台站的相位差(Jin,Gaherty,2015)。去除了震源初始相位和仪器响应导致的相位差后,不同台站之间的相对走时${t}_{\mathrm{r}} $可以表示为
$$ {t}_{\mathrm{r}}={t}_{2\mathrm{D}} + {t}_{\sigma } + \Delta t \text{,} $$ (2) 式中:$ {t}_{2\mathrm{D}} $为面波从震源传播至不同台站所需时间差;$ {t}_{\sigma } $为测量误差,通过比较临近台站的测量值,该项约为4 s (Ekström,2011;Ma et al,2014);$ \Delta t $则为时钟错误导致的相对走时差,时钟正常台站的$ \Delta t $为0。为避免周波跳跃,我们先通过式(3)与式(4)去除根据已有相速度图预测的到时差
$$ \widetilde {{t}_{\mathrm{r}}}={t}_{\mathrm{r}}-{t}_{2{\mathrm{D}}\_{\mathrm{start}}}= \widetilde {{t}_{2\mathrm{D}}} + {t}_{\sigma } + \Delta {t} \text{,} $$ (3) $$ {t}_{2\mathrm{D}\_\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}={\int }_{0}^{\varDelta }\frac{\mathrm{d}l}{c ( \theta \text{,} \phi ) } \text{,} $$ (4) 式中,$\widetilde {{t}_{2\mathrm{D}}} $为台站间实际到时差,$ c ( \theta \text{,} \phi ) $为在某个频率下对应于经纬度$ ( \theta \text{,} \phi ) $的起始相速度(本文取10 mHz),$\varDelta $为震中距。去除已有相速度图后的绝大部分相对走时测量值$ \widetilde {{t}_{\mathrm{r}}} $,理论上其应当处于一个较小的数值范围,少部分地震波数据预测的$ \widetilde {{t}_{2\mathrm{D}}} $受地球结构影响大,导致$ \widetilde {{t}_{\mathrm{r}}} $超过理论范围,但该部分测量值可以在面波波速反演中得到较好的拟合,对应的拟合残差较小。
当台站在较长时间(数天乃至数月)内出现较大时钟偏差量∆t 时(10 s以上,3.3节中有更详细的论述),相对走时测量值会连续出现较大异常值,且其对应的拟合残差不会回落到较小的数值范围。当台站发生极性反转时,由互相关函数计算出的相对走时会包含半个周期的固定偏差,在反演中无法被拟合以致其残差仍为较大的异常值。因此,无论是时钟错误还是极性反转问题,仅在一个频率下对波形进行相位测量即可检测。而由于10 mHz频率波形的测量技术成熟,不易发生周期跳跃,且极性反转在周期为100 s的波形互相关测量中表现为半个周期50 s的偏差,利于进行观察识别。因此,我们选择该频率波形进行检测。
然而,面波相对走时测量值会受一些因素的影响:① 面波相位测量之前需要利用震源位置参数去除震源初始相位$ {\varphi }_{\mathrm{s}} $,而震源水平定位存在的误差,有时会达到20—50 km,进而造成相对走时的测量值出现偏差。事实上,震源定位的误差呈随机分布,由其产生的偏差也是随机出现。针对震源定位误差导致的偏差进行测试结果显示,相较于时钟错误和极性反转表现出的连续且固定的偏差,震源定位导致的偏差明显不同(附录1)。② 修正震源初始相位$ {\varphi }_{\mathrm{s}} $时利用了1D初始参考地球模型(preliminary reference earth model,缩写为PREM),因此需要考虑地下局部强速度结构异常的影响。此影响可以分为两种情况:一是台站正下方存在局部强速度结构异常,则该台站全部的相对走时测量值都会包含由其引起的固定偏差,很容易从检测结果中排除掉该情况;二是由于局部强速度结构不处于台站正下方,地震波的传播路径穿过该局部强速度结构是随机的,以致引起的偏差也随机出现。这两种情况都无法对检测结果产生影响。③ MW7.5以上地震的破裂尺度大且过程复杂,会对相位校正造成影响(Ma,Masters,2015),但是本文所选用的数据中,2008—2012年间MW7.5以上地震仅为26次,相对于MW5.5以上地震的2 741次,数量非常少且不会集中在同一时期出现,因此对检测结果的影响可以忽略。
为了进一步排除利用面波相对走时检测可能导致的误判,同时区分时钟错误和极性反转,本文选择结合预测P波到时对初筛的异常数据进行二次检验。目前对全球P波波速的研究已经具有相对高的精准度,远震P波到时相较于使用1D PREM模型计算的理论P波到时的偏差通常小于10 s (Xia et al,2015)。本文选用异常时段内具有清晰P波震相的全部地震事件并计算其理论P波到时,若台站出现时钟错误,其地震记录的P波震相较于预测到时会出现明显的提前或延后;若台站记录发生极性反转,则地震记录P波震相与预测到时线会基本一致。
3. 数据与应用
3.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年间MW5.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下载方式 第1次用时/s 第2次用时/s 第3次用时/s 平均用时/s BREQ_FAST 216.67 216.07 217.93 216.89 ObsPy 566.50 570.35 568.56 568.47 表 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 3.2 测量相对走时
3.2.1 预处理
数据的预处理流程如下:① 将原始波形利用傅里叶变换从时间域变换到频率域,利用窄带通高斯滤波器将其滤波到10 mHz频率(Dziewonski et al,1969),并在频率域去除仪器响应并进行归一化处理,再利用反傅里叶变换将其转换回时间域;② 进行震源初始相位的修正,面波辐射花样去除与频散项相位修正,每个震源−台站对的震源初始相位$ {\varphi }_{\mathrm{s}} ( \omega ) $计算所使用的模型均为PREM,震源位置与震源深度参数使用全球矩心矩张量(Global Centroid-Moment-Tensor ,缩写为GCMT)公布的数据;③ 利用已有的全球10 mHz面波相速度图(Ekström,2011)计算出理论走时,对同一地震事件的地震波数据进行初步时移。
3.2.2 测量
在波形测量过程中,为了避免小部分低信噪比记录干扰地震波形的相关性,我们首先选取信噪比大于4的数据(Houser et al,2008),信噪比定义为等长的信号段与非信号段各点幅值平方和的比值。对不同台站记录的同一事件波形进行互相关计算,识别对应着不同延迟时间的正峰值,之后利用VanDecar和Crosson (1990)提出的时域法计算同一地震事件不同台站的相对走时。
相对走时的测量使用聚类分析技术,即根据互相关系数衡量波形间的相似性,利用Hartigan (1975)的聚类算法生成类簇与聚类树。如图1所示,不同台站记录的同一地震事件波形聚类。利用相对走时将波形进行时移并将其绘在图形界面的左边,相对应的聚类树绘于图形界面的右边。在此阶段中,扭曲失真的波形记录和噪音记录可以被识别出并舍弃。利用聚类算法可以实现一次处理数百道地震波形数据,仅需要人工监督聚类时的分簇情况并做出保留和舍弃群簇的决策,实现了高效率高质量的测量。
图 1 一个波形聚类示例图此示例中有827道波形,程序允许点击右侧的聚类树来选择合适的截止值。垂直黄线为此示例的合理截止值Figure 1. An example of waveform clusteringThere are 827 waveforms in this example. The program allows to select the appropriate cut-off value by clicking oon the cluster tree on the right side with the mouse. The vertical yellow line represents a reasonable cutoff for this example3.3 面波相对走时筛选
本研究共进行了116万7 072道地震面波记录相对走时的测量,并统计分析了全部测量值。如图2所示,地震面波的整体测量值的平均值为−0.040 7,标准差为10.18。这表明大部分相对走时测量值介于−10 —10 s,过小的时钟偏差无法识别。本文的方法仅能对大于10 s的时钟偏差进行有效检测,因此将−10 —10 s定义为正常值界。
使用本次测量得到的10 mHz面波数据集进行全球波速初步反演,反演的详细步骤参见文献(Ma et al,2014)。反演结果如图3所示,成像结果和拟合残差均表明本次测量的数据集整体质量良好。
对于单个台站,我们将相对走时测量值和拟合残差按其对应发生地震事件的时间顺序进行排列。如图4所示,时钟正常台站AAK的相对走时测量值绝大多数处于−10—10 s的正常范围内;由实际传播路径造成的较大相对走时测量值对应的拟合残差回落到正常范围内,如第50天的数据点所示;而受测量误差和震源定位误差等因素的干扰,随机地出现个别相对走时测量值和拟合残差均异常大的情况,如第110天的数据点。
然而,当台站在一段时期内发生时钟错误,则相对走时和拟合残差会在时钟错误期间连续超出正常范围。为了对比说明,本文选取了位于欧洲中部相距约1.4° (一个波长约4°)的两个台站进行对比。图5a与5b分别展示了台网G的时钟错误台站ECH与台网GE的正常台站STU在2011年前80天的面波相对走时测量值与拟合残差的折线图。由图可见,台站ECH在2011年的第40—60天期间,相对走时测量值与拟合残差全部为远超10 s的异常值;对比台站STU在同一时期的相对走时测量值与拟合残差(图5b),可以看出台站ECH在此段时间发生了较大时钟偏差。为了证明检测结果是否准确可靠,本文将全部的钟错台站与距其1—2个波长内的正常台站进行了对比,具体内容见附录2。
图6显示出台网II的ERM台站在2011上半年的相对走时测量值和拟合残差均超出了正常范围,在50 s或−50 s附近。测量的10 Hmz面波数据周期为100 s,这种半个周期的固定偏差说明台站的波形记录发生了极性反转。下一步我们将利用预测P波到时给出更准确的判断依据。
3.4 预测P波到时检验
利用面波相对走时测量值的异常,已经从大量地震数据初步检测出少量台站发生时钟错误和极性反转的时段。从初步检测出的异常数据中挑选信噪比大于9且震中距在20°—140°之间的地震记录,以保证可以观察到清晰可辨的P波震相。将波形滤波为低频率波0.2 Hz,之后根据PREM模型预测的理论到时将波形记录进行时移。时钟正常台站记录的P波震相会清晰地排列于预测P波到时线之上;而时钟错误台站记录的P波震相则会较理论P波到时线出现明显提前和延后的情形。如图7所示:由面波相对走时初步筛选出的一个时钟错误的台站ULN(图7a),并对比其时钟出错时段与时钟正常时段中的地震事件波形记录(图7b、c)。将时钟正常时段的地震记录根据预测P波到时进行时移后,其P波震相与其它正常台站地震记录的P波震相对齐于预测P波到时线上。而时钟出错时段的地震记录时移后,其P波震相较预测P波到时线延后约20 s。对于每个筛查出的异常时段,选取其中全部而非单个具有清晰P波震相的波形记录进行如上的检测,使最终的检测结果更加准确可靠。
图 7 台站ULN的相对走时和拟合残差折线图(a)与示例地震事件5个台站的波形记录(b,c)(a) 2010年第150—230天台网IU的台站ULN的相对走时测量值与拟合残差折线图;(b) 2010年第167天MW6.2地震的5个台站的波形记录;(c) 2010年第222天MW7.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 station ULN of network IU during the period from 150th to 230th day in 2010;(b) The waveforms of the MW6.2 earthquake recorded by five stations on the 167th day in 2010;(c) The waveforms of the MW7.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天MW7.3地震的波形图Figure 8. Line chart of relative traveltime measurements and fitting residuals and waveforms of the example eventFigs.(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 MW7.3 earthquake recorded by five stations on the 175th day of 2011,respectively4. 结果
表3和表4分别列出了台站时钟错误和波形极性反转出现时段的最终检测结果。结果表明,在检测的2 010个地震台站中,有12个台站的14个时段出现大于10 s的时钟错误,7个台站的8个时段发生极性反转。
表 3 本文检测到的台站出现时钟错误的时段Table 3. Periods of clock errors stations台网 台站 钟错时段 台网 台站 钟错时段 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天表 4 本文检测到的台站出现极性反转的时段Table 4. Periods of polarity reversal stations台网 台站 极性反转转时段 台网 台站 极性反转时段 NZ BKZ 2 008年—2 012年 II ERM 2 008年—2 012年第303天 NZ KHZ 2 008年—2 011第54天 NZ HIZ 2 008年—2 012年 NZ OUZ 2 008年—2 010第55天 NZ ODZ 2 008年第1—190天 NZ QRZ 2 008年第1—第232天 NZ QRZ 2 008年第290天—2 009年第22 0天 5. 讨论与结论
钟错时段的检测结果全部经过了P波到时的检验。台站时钟出错的时段短期的仅有数天,较长期的可达数月。本文提出的方法只对偏差量较大(>10 s)的时钟错误进行有效检测,即使对于使用大量数据的全球面波波速成像研究,这种程度的时钟偏差量也会在反演中产生不可忽略的影响。另一方面,在数据量使用较小的区域性地震成像研究中,此类较大时钟偏差可能会对结果产生严重影响。在预处理阶段去除时钟错误的数据可以避免反演结果中出现错误的成像结构。
本文的方法利用了MW5.5以上地震激发的在全球范围内可观测的长周期信号,由于全球MW≥5.5地震事件的发生频率约为1.2次/天,此方法可以给出以天数为分辨率的时段结果。同样是应用于超远台站间距的检测方法,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) 台站KOWA2012年下半年的相对走时测量值与拟合残差折线图;(b) 5个台站记录的2012年第302天MW7.8地震的波形记录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 MW7.8 earthquake on the 302th day of 2012 recorded by five stations对于波形极性反转的检测结果,我们使用远震初至P波波形进行了检验。选取位于新西兰岛的两个正常台站SNZO和URZ和一个发生极性反转的台站BKZ。如图10所示,3个台站相对于2008年MW7.9汶川地震震源方位角非常接近,台站BKZ出现明显的极性反转(本文所说的极性反转为垂直分量极性反转,并不涉及水平分量的方位角摆放错误)出现极性反转的台站绝大部分归属台网NZ,我们就该问题联系了台网NZ的运营商GeoNet (Geological hazard monitoring system in New Zealand)。Geonet的回复证实了台网NZ部分台站的LHZ分量出现了极性反转的问题,并提供了一份已查明的发生极性反转的台站统计表单。与该表单的对比证明本文关于台网NZ发生极性反转台站的检测结果准确无误。对于垂直分量出现极性反转的原因,GeoNet认为相较于地震计的安装问题,布设时线缆连接错误的可能性更大。
图 10 2008年汶川MS8.0地震事件及台站位置(a)与新西兰岛三个台站记录到的波形(b)(波形数据仅作0.01—0.02 Hz滤波处理)Figure 10. Example event and stations location map and corresponding waveforms diagram(a) The focal mechanism and location of the selected earthquakes,and the distribution locations of the three stations. Figure (b) shows the first-arrival of the 2008 MW7.9 earthquake recorded by three stations,and the azimuths of the stations relative to the hypocentral location (Waveform data is only filtered 0.01—0.02 Hz)利用地震事件的波形数据,基于面波和P波技术检测时钟错误各有优劣:面波的能量较强,传播距离远,可以满足对全球范围内台站的检测。而由于面波测量时干扰因素的存在,难以保证结果的准确无误;而利用P波的检测技术虽然检测精度较高,但P波长距离传播后的低信噪比无法满足该方法对清晰P波震相的要求。故而本文选择将二者相结合,以达到更好的检测效果。然而,基于地震面波技术检测时钟偏差有其自身的局限性。受到测量误差、震源参数误差等因素的影响,该方法只能对较大的时钟偏差进行检测;此外,该方法适用于长时段的时钟偏差效应,对于单个事件的异常或者小于1天的钟错时段,无法做到有效检测,需结合其它资料进一步验证。针对目前方法的不足,未来可以考虑近一步的改进。
本文提出利用面波相对走时测量值和P波到时相结合的方法实现对全球范围的地震台站同时检测出现时钟出错和极性反转的时段。该方法利用聚类分析进行面波相对走时的高效测量,进而实现了对大量地震台站数据准确快速地检测,而P波到时的二次检验保证了最终检测结果的准确性。之后使用该方法对全球30个台网的2 010个台站在2008—2012年间≥MW5.5地震的长周期面波波形数据进行了检测应用。结果表明,其中12个台站的14个时段出现了大于10 s的时钟错误;7个台站的8个时段出现了极性反转。对实际波形数据的应用表明本文的方法能够快速准确地检测出地震台站的时钟错误和极性反转时段,在提高全球面波地震成像质量的同时也为使用相关数据的震源机制反演等研究提供了参考。而该方法的不足之处在于,难以处理单个事件或小于1天的时钟错误和极性反转问题,无法对数秒或更小的时钟偏差进行有效检测。
附录1
通过对去除初始相位时使用的震源位置参数添加50 km以内的水平随机误差,以定量地评估震源水平定位误差对面波相对走时测量的影响。此外,我们对震源深度参数添加±10 km随机误差以检测震源深度误差对面波相对走时测量的影响。无论是水平定位误差还是深度误差,对面波相对走时的测量值的影响都很小,少量较大异常值也呈随机出现,不会对最终结果造成影响。具体结果如图1所示。
图 1 正常台站AAK (上行)与钟错台站ECH (下行)的测量对比(a)正常测量;(b)震源参数添加50 km内的随机水平定位误差;(c)震源参数添加$ \pm 10 $ km内的随机深度定位误差Figure 1. Measurement comparison between normal station AAK (upper) and clock error station ECH (bottom)(a) are normal measurements;(b) are source parameters with random horizontal localization errors within 50 km added;(c) are source parameters with random depth localization errors within ±10 km added附录2
我们将台站间距与对应的相对走时测量值之差进行了分析(图2),台间距小于10°的相对走时测量值之差大致在4 s之内。将全部的钟错台站与距其1$ — $2个波长之内台站(表1的测量值和拟合残差进行了对比(图4)。台站KOWA与SPB附近没有可用来对比的正常台站(图3),但文中已经将KOWA与前人关于时钟错误的研究结果进行了对比分析。
表 5 对比台站对及其台间距Table 5. 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° -
图 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 oon the cluster tree on the right side with the mouse. The vertical yellow line represents a reasonable cutoff for this example
图 7 台站ULN的相对走时和拟合残差折线图(a)与示例地震事件5个台站的波形记录(b,c)
(a) 2010年第150—230天台网IU的台站ULN的相对走时测量值与拟合残差折线图;(b) 2010年第167天MW6.2地震的5个台站的波形记录;(c) 2010年第222天MW7.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 station ULN of network IU during the period from 150th to 230th day in 2010;(b) The waveforms of the MW6.2 earthquake recorded by five stations on the 167th day in 2010;(c) The waveforms of the MW7.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
图 8 出现记录错误的台站的相对走时测量值、拟合残差折线图与示例地震事件波形图
图(a)和(b)分别为台站ERM与GRGR在2011年第130—第230天期间相对走时测量值与拟合残差折线图;图(c)和(d)分别为台站ERM和GRGR与对比台站记录的2011年第175天MW7.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 MW7.3 earthquake recorded by five stations on the 175th day of 2011,respectively
图 9 台站KOWA 2012年相对走时测量值、拟合残差折线图(a)与示例地震事件波形图(b)
(a) 台站KOWA2012年下半年的相对走时测量值与拟合残差折线图;(b) 5个台站记录的2012年第302天MW7.8地震的波形记录
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 MW7.8 earthquake on the 302th day of 2012 recorded by five stations
图 10 2008年汶川MS8.0地震事件及台站位置(a)与新西兰岛三个台站记录到的波形(b)(波形数据仅作0.01—0.02 Hz滤波处理)
Figure 10. Example event and stations location map and corresponding waveforms diagram
(a) The focal mechanism and location of the selected earthquakes,and the distribution locations of the three stations. Figure (b) shows the first-arrival of the 2008 MW7.9 earthquake recorded by three stations,and the azimuths of the stations relative to the hypocentral location (Waveform data is only filtered 0.01—0.02 Hz)
图 1 正常台站AAK (上行)与钟错台站ECH (下行)的测量对比
(a)正常测量;(b)震源参数添加50 km内的随机水平定位误差;(c)震源参数添加$ \pm 10 $ km内的随机深度定位误差
Figure 1. Measurement comparison between normal station AAK (upper) and clock error station ECH (bottom)
(a) are normal measurements;(b) are source parameters with random horizontal localization errors within 50 km added;(c) are source parameters with random depth localization errors within ±10 km added
表 1 BREQ_FAST和ObsPy两种方式在相同网络环境中的下载用时统计
Table 1 The download time of BREQ_FAST and ObsPy in the same network environment
下载方式 第1次用时/s 第2次用时/s 第3次用时/s 平均用时/s BREQ_FAST 216.67 216.07 217.93 216.89 ObsPy 566.50 570.35 568.56 568.47 表 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 表 3 本文检测到的台站出现时钟错误的时段
Table 3 Periods of clock errors stations
台网 台站 钟错时段 台网 台站 钟错时段 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天表 4 本文检测到的台站出现极性反转的时段
Table 4 Periods of polarity reversal stations
台网 台站 极性反转转时段 台网 台站 极性反转时段 NZ BKZ 2 008年—2 012年 II ERM 2 008年—2 012年第303天 NZ KHZ 2 008年—2 011第54天 NZ HIZ 2 008年—2 012年 NZ OUZ 2 008年—2 010第55天 NZ ODZ 2 008年第1—190天 NZ QRZ 2 008年第1—第232天 NZ QRZ 2 008年第290天—2 009年第22 0天 表 5 对比台站对及其台间距
Table 5 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° -
王俊,郑定昌,詹小艳,江昊琳,缪发军,高景春. 2013. 基于背景噪声互相关格林函数的单台时间误差估计[J]. 地震学报,35(6):888–901. doi: 10.3969/j.issn.0253-3782.2013.06.012 Wang J,Zheng D C,Zhan X Y,Jiang H L,Miao F J,Gao J C. 2013. Estimation of time error for a single station based on Green’s function by ambient noise cross-correlation[J]. Acta Seismologica Sinica,35(6):888–901 (in Chinese).
Abbas A,Zhu G H,Zi J P,Chen H,Yang H F. 2023. Evaluating and correcting short-term clock drift in data from temporary seismic deployments[J]. Earthquake Research Advance,3(2):100199. doi: 10.1016/j.eqrea.2022.100199
Anchieta M C,Wolfe C J,Pavlis G L,Vernon F L,Eakins J A,Solomon S C,Laske G,Collins J A. 2011. Seismicity around the Hawaiian islands recorded by the PLUME seismometer networks:Insight into faulting near Maui,Molokai,and Oahu[J]. Bull Seismol Soc Am,101(4):1742–1758. doi: 10.1785/0120100271
Bondár I,Storchak D. 2011. Improved location procedures at the international seismological centre[J]. Geophys J Int,186(3):1220–1244. doi: 10.1111/j.1365-246X.2011.05107.x
Butler R,Lay T,Creager K,Earl P,Fischer K,Gaherty J,Laske G,Leith B,Park J,Ritzwolle M,Tromp J,Wen L X. 2004. The global seismographic network surpasses its design goal[J]. Eos,Trans,Am Geophys Union,85(23):225-229.
Dziewonski A,Bloch S,Landisman M. 1969. A technique for the analysis of transient seismic signals[J]. Bull Seismol Soc Am,59(1):427–444. doi: 10.1785/BSSA0590010427
Ekström G,Nettles M,Abers G A. 2003. Glacial earthquakes[J]. Science,302(5645):622–624. doi: 10.1126/science.1088057
Ekström G,Dalton C A,Nettles M. 2006. Observations of time-dependent errors in long-period instrument gain at global seismic stations[J]. Seismol Res Lett,77(1):12–22. doi: 10.1785/gssrl.77.1.12
Ekström G,Nettles M. 2010. Performance of the GSN station CASY-IU,1996-2009[R]. New York:Waveform Quality Center.
Ekström G. 2011. A global model of Love and Rayleigh surface wave dispersion and anisotropy,25-250 s[J]. Geophys J Int,187(3):1668–1686. doi: 10.1111/j.1365-246X.2011.05225.x
Gouédard P,Seher T,McGuire J J,Collins J A,van der Hilst R D. 2014. Correction of ocean-bottom seismometer instrumental clock errors using ambient seismic noise[J]. Bull Seismol Soc Am,104(3):1276–1288. doi: 10.1785/0120130157
Hartigan J A. 1975. Clustering Algorithms[M]. New York:John Wiley and Sons:351.
Houser C,Masters G,Shearer P,Laske G. 2008. Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms[J]. Geophys J Int,174(1):195–212. doi: 10.1111/j.1365-246X.2008.03763.x
Hurst A W,Bibby H M,Robinson R R. 2002. Earthquake focal mechanisms in the central Taupo volcanic zone and their relation to faulting and deformation[J]. N Z J Geol Geophys,45(4):527–536. doi: 10.1080/00288306.2002.9514989
Jin G,Gaherty J B. 2015. Surface wave phase-velocity tomography based on multichannel cross-correlation[J]. Geophys J Int,201(3):1383–1398. doi: 10.1093/gji/ggv079
Lentas K. 2021. Indications of seismic station phase reversals detected from parametric data in the ISC bulletin[J]. J Seismol,25(1):1–23. doi: 10.1007/s10950-020-09960-1
Ma Z,Masters G. 2014. A new global Rayleigh- and love-wave group velocity dataset for constraining lithosphere properties[J]. Bull Seismol Soc Am,104(4):2007–2026. doi: 10.1785/0120130320
Ma Z T,Masters G. 2015. Effect of earthquake locations on Rayleigh wave azimuthal anisotropy models[J]. Geophys J Int,203(2):1319–1333. doi: 10.1093/gji/ggv369
Ma Z T,Masters G,Laske G,Pasyanos M. 2014. A comprehensive dispersion model of surface wave phase and group velocity for the globe[J]. Geophys J Int,199(1):113–135. doi: 10.1093/gji/ggu246
Niu F L,Li J. 2011. Component azimuths of the CEArray stations estimated from P-wave particle motion[J]. Earthq Sci,24(1):3–13. doi: 10.1007/s11589-011-0764-8
Petersen G M,Cesca S,Kriegerowski M,the AlpArray Working Group. 2019. Automated quality control for large seismic networks:Implementation and application to the AlpArray seismic network[J]. Seismol Res Lett,90(3):1177–1190. doi: 10.1785/0220180342
Pollitz F F,Stein R S,Sevilgen V,Bürgmann R. 2012. The 11 April 2012 east Indian Ocean earthquake triggered large aftershocks worldwide[J]. Nature,490(7419):250–253. doi: 10.1038/nature11504
Priestley K,McKenzie D,Ho T. 2018. A lithosphere-asthenosphere boundary:A global model derived from multimode surface‐wave tomography and petrology[M]//Yuan H Y,Romanowicz B. Lithospheric Discontinuities. Hoboken:John Wiley & Sons:111-123.
Roman D C,Neuberg J,Luckett R R. 2006. Assessing the likelihood of volcanic eruption through analysis of volcanotectonic earthquake fault-plane solutions[J]. Earth Planet Sci Lett,248(1/2):244–252.
Schneider J F,Aster R C,Powell L A,Meyer R P. 1988. Timing of portable seismographs from omega navigation signals[J]. Bull Seismol Soc Am,78(2):1035–1035. doi: 10.1785/BSSA0780021035
Sens-Schönfelder C. 2008. Synchronizing seismic networks with ambient noise[J]. Geophys J Int,174(3):966–970. doi: 10.1111/j.1365-246X.2008.03842.x
Stehly L,Campillo M,Shapiro N M. 2007. Traveltime measurements from noise correlation:Stability and detection of instrumental time-shifts[J]. Geophys J Int,171(1):223–230. doi: 10.1111/j.1365-246X.2007.03492.x
Tape C,Liu Q Y,Maggi A,Tromp J. 2009. Adjoint tomography of the Southern California crust[J]. Science,325(5943):988–992. doi: 10.1126/science.1175298
Tian Y,Ning J Y,Yu C Q,Cai C,Tao K. 2013. Focal mechanism solutions of the 2008 Wenchuan earthquake sequence from p-wave polarities and SH/P amplitude ratios:New results and implications[J]. Earthquake Science,26(6):357–372. doi: 10.1007/s11589-014-0067-y
VanDecar J C,Crosson R S. 1990. Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least squares[J]. Bull Seismol Soc Am,80(1):150–169.
Xia Y J,Ni S D,Zeng X F,Xie J,Wang B S,Yuan S Y. 2015. Synchronizing intercontinental seismic networks using the 26 s persistent localized microseismic source[J]. Bull Seismol Soc Am,105(4):2101–2108. doi: 10.1785/0120140252
Yang W Z,Hauksson E,Shearer P M. 2012. Computing a large refined catalog of focal mechanisms for Southern California (1981—2010):Temporal stability of the style of faulting[J]. Bull Seismol Soc Am,102(3):1179–1194. doi: 10.1785/0120110311
Young B A,Chen K C,Huang B S,Chiu J M. 2017. Identification and elimination of data peculiarities in the strong‐motion downhole array in Taipei Basin[J]. Seismol Res Lett,88(1):82–95. doi: 10.1785/0220160164
Yukutake Y,Iio Y. 2017. Why do aftershocks occur? relationship between mainshock rupture and aftershock sequence based on highly resolved hypocenter and focal mechanism distributions[J]. Earth Planets Space,69(1):68. doi: 10.1186/s40623-017-0650-2