2020, Vol. 40
钻孔应变仪一般埋深几十米以上，既可排除部分地表噪声，又能记录到地震仪可能遗漏的慢地震、蠕滑等物理事件。应变仪在国内外地球物理事件研究中具有大量应用，邱泽华等[1-3]对姑咱台YRY-4型应变仪记录到的汶川地震数据进行分析，认为其可记录到震前的应变信号；杨选辉等[4]对倾斜应变高频观测进行研究，发现倾斜应变仪的地震记录与附近摆式地震仪的地震记录频谱非常一致；Barbour等[5]基于板块观测计划(PBO)的应变仪数据对南加利福尼亚地区的地震进行同震分析；Johnston等[6]结合连续钻孔应变数据对Parkfield地震进行近场分析，同时对震颤及断层响应进行研究。距离震中较近的应变仪可记录到震前应变信号的观点虽然已经得到很多实例证明，但仍然需要更具说服力的例子。高采样应变仪的记录中包含更多的低频成分，但应变仪记录的波形振幅较小，会对研究造成一定困难。
近年来，钻孔应变观测已经成为热门课题，国内外学者也都在积极研究如何才能使钻孔应变观测数据发挥更大的价值。本文旨在通过部分台站的钻孔应变仪记录的同震信号和QSSP计算得到的理论应变数据，分析两者的体波段互相关情况并计算两者的互相关系数，探究钻孔应变观测在震源反演和地球结构模型检验方面的可靠性，为今后钻孔应变观测在该方面的研究提供参考。
1 印尼苏门答腊海域7.8级地震同震信号分析2016-03-02印度尼西亚苏门答腊岛西部海域发生7.8级地震，Global CMT公布的发震时刻为12:50:5.4 GMT(北京时间2016-03-02 20:50:5.4)，震中坐标为4.75°S、94.22°E，震源深度为37.2 km。
1.1 数据预处理本文采用9个台站的钻孔应变仪记录的同震信号数据，为便于研究和描述，计算各台站的震中距并对台站进行编号，见表 1，台站和震中位置分布见图 1。
|
表 1 台站信息表 Tab. 1 Stations information |
|
图 1 钻孔应变仪台站及震中位置分布 Fig. 1 Borehole strainmeter stations and epicenter location |
钻孔应变仪记录的原始数据包含长周期的潮汐和非潮汐信号，为提取有效的目标地震信号，需对数据进行处理。本次选用东西分量和北南分量秒采样数据进行研究，时间范围为北京时间2016-03-02 00:00:00~03-03 23:59:59。针对存在突跳的分量，需要去除突跳并进行有效的插值，以减小突跳对后续研究的影响。各台站各分量均存在明显的潮汐变化，需要对信号进行高通滤波处理，以去除低频潮汐和非潮汐长周期信号。本文设计FIR汉宁高通滤波器的阻带截止频率为0.02 Hz，最小阶数为320[7]。
1.2 同震信号震相识别与到时分析监测信号的可靠性不仅与监测仪器的精度和灵敏度有关，还与监测仪器的授时系统是否准确有很大关系。为判断本文研究数据的可靠性，需要对钻孔应变仪记录的数据进行应变地震波到时分析，对可能存在的时间差进行讨论和修正。同时，为选取互相关分析数据，需要对各分量数据进行P波震相识别。
采用TauP软件[8]对8个台站的P波到时进行计算，结果见表 2。若该台站2个方向上均可识别P波到时，则表中数据为两者平均值；若只有1个方向可识别，则表中数据为该方向P波到时。将8个台站按震中距由小到大排列，把钻孔应变仪记录到的2个分量的数值分别绘制成图，横坐标零点为地震发震时刻。由于嘉祥台2个分量都未记录到此次地震，故不可用于后续研究，需要剔除。在北南分量中(图 2)，除荣成台外，P波到时均可清晰识别，营口台实际到时与理论到时接近，其他台站实际到时与理论到时差别较大；在东西分量中(图 3)，榆树沟台和荣成台无法清晰识别P波到时，营口台实际到时与理论到时接近，其他台站实际到时与理论到时差别较大。各台站P波到时时差见表 2。
|
图 2 北南分量 Fig. 2 NS component |
|
图 3 东西分量 Fig. 3 EW component |
|
|
表 2 各台站P波理论与实际到时 Tab. 2 Theoretical and actual arrival time of P wave of each station |
从理论上分析，P波到时应随震中距的增大而增大，而实际得到的结果与理论分析差别较大，且不存在规律性。除营口台实际到时与理论到时接近外，其他台站均有较大程度的延后，可能是因为应变仪的授时系统出现故障，但8个台站中7个同时出现故障的概率较小，所以最有可能的原因为台站的授时系统长时间未校准。
2 基于QSSP的理论应变数值计算 2.1 QSSP简介QSSP程序由GFZ中心汪荣江教授编写，使用一种基于球对称和自重的地球模型计算完整的合成地震图的混合方法，该地球模型具有大气、海洋、地幔、液态外核和固态内核多层结构[9]。QSSP可以计算包含应变在内的多种参数，应变计算为某一点上从地震发生至指定时间范围内的应变值变化。
2.2 理论应变地震波到时分析和震相识别将QSSP计算所得的应变数据进行P波震相识别，在北南分量中(图 4)各波均可清楚识别，在东西分量中(图 5)灵山台、信宜台、库米什台和榆树沟台未能识别到P波。QSSP计算的理论数据所显示的各波到时符合表 2中的理论到时。
|
图 4 北南理论分量 Fig. 4 NS theoretical component |
|
图 5 东西理论分量 Fig. 5 EW theoretical component |
表 3为2个信号相关程度类型其中ρxy为2个信号的互相关程度，ρxy值越大，表明2个信号波形越相似。ρxy=1时，2个信号完全正相关；ρxy=－1时，2个信号完全负相关；ρxy=0时，2个信号完全不相关。
|
|
表 3 2个信号相关程度类型 Tab. 3 Correlation types of two signals |
根据本文监测数据和理论数值进行到时分析和震相识别的结果，对可以清晰识别P波的监测数据和理论数据进行相关性分析，各台站互相关分析所选取的分量见表 4。
|
|
表 4 各台站选取分量汇总 Tab. 4 Components selected from each station |
信宜台、库米什台和榆树沟台东西分量理论数据未识别到P波，荣成台东西分量和北南分量的监测数据未识别到P波，故无法对以上台站的相关分量数据进行互相关分析。
根据监测数据到时分析可知，除营口台实际到时与理论到时接近外，其他台站的实际到时与理论到时存在很大差别。为最大程度地减少互相关分析的干扰因素，得到有效的互相关分析结果，采取对齐P波到时的方法对监测数据进行截取，使得截取后的监测数据P波到时与理论数据P波到时保持一致。对截取后的数据进行互相关分析，结果见图 6、7(按震中距从小到大排列)。
|
图 6 北南分量比较 Fig. 6 Comparison of NS component |
|
图 7 东西分量比较 Fig. 7 Comparison of EW component |
对比图 6、7可以发现，监测数据中存在许多高频干扰，这些干扰可能会对互相关分析造成较大影响。由于已经对数据进行高通滤波，因此可使用FIR汉宁低通滤波器对监测数据进行二次滤波[7]，而各台站受到的高频干扰不同，所需滤波器的截止频率(Wc)也不同。将滤波后的数据再与理论数据进行互相关分析，结果见图 8、9(按震中距从小到大排列)。
|
图 8 北南分量比较 Fig. 8 Comparison of NS component |
|
图 9 东西分量比较 Fig. 9 Comparison of EW component |
根据互相关分析结果可以发现，在本文研究数据范围内，监测数据与理论数据的互相关类型有显著相关、高度相关和微相关3类，其中6个分量显著相关，3个分量高度相关，1个分量微相关(表 5)，互相关系数越大，表明实际应变地震波的路径与QSSP所采用的模型越接近。
|
|
表 5 低通滤波前后监测数据与理论数据互相关系数汇总 Tab. 5 Summary of correlation coefficient between monitoring data and theoretical data before and after low-pass filtering |
在北南分量中，除灵山台、信宜台和榆树沟台外，其他台站的互相关系数都接近高度相关。信宜台互相关系数较小的原因不明，需要继续结合其他资料进行研究；而榆树沟台互相关系数较小的原因可能是受到的干扰较多，这一点可从图 2中得到验证。在东西分量中，除营口台外，汕头台和会昌台都为高度相关，营口台东西分量是唯一一个微相关的分量，原因不明，还需要结合其他资料继续进行研究。
对比低通滤波前后互相关系数发现，除营口台东西分量外，其余各分量均有明显提高，这也验证了此前高频干扰会对互相关系数产生影响的猜测。在此次互相关分析所采用的6个北南分量中，除信宜台和榆树沟台由于干扰较大造成互相关系数较小外，可以发现，在震中距36°~52°范围内，各台站互相关系数约为0.8；在3个东西分量中，除营口台外，各分量的互相关系数都大于0.8。
综上所述，震中距在36°~52°范围内的钻孔应变观测数据与理论应变数据具有较好的一致性。测震观测地震波是地震反演中经典且成熟的应用资料，一般来说，其他高数采率仪器记录到的观测地震波与测震观测地震波相关程度高，则说明其他仪器记录到的观测地震波可用于地震反演。考虑到测震仪与钻孔应变仪所观测的物理量不同，仪器原理和整个观测系统的频率特性也不同，而由QSSP正演所得到的理论应变地震波属于应变物理量，与钻孔应变仪观测的物理量相同，且不涉及仪器原理和整个观测系统的频率特性，所以认为观测应变地震波可用于地震反演，同时应变观测地震波与理论应变地震波的相关程度比应变观测地震波与测震观测地震波的相关程度更高。
4 结语本文对印尼苏门答腊海域7.8级地震秒采样钻孔应变仪的监测数据与QSSP计算所得的理论应变数据的互相关情况进行分析，结果表明，在排除已知干扰的情况下，震中距在36°~52°范围内的钻孔应变观测数据可用于震源反演和新地球结构模型的检验，且结果可靠。另外，本文研究中也存在以下问题需要改进：
1) 本文选取的台站监测数据受到的干扰较大，虽然已进行滤波等一系列处理，但不能完全排除这些干扰对研究结果的影响。在后期的研究中，可选用受干扰较小的台站的监测数据，进一步减小对研究结果的影响，得到更具说服力的结论。
2) 本文选用的数据过少，且缺少高纬度和低纬度地区台站的监测数据，在后期的研究中可以选用纬度分布更加全面的台站数据。
3) 由于钻孔应变仪对体波波段的记录能力明显不如其对面波的响应，该特性对本次研究中的震相识别造成一定困难。同时，在震相识别过程中可以发现，大部分台站可以清晰识别各震相，P波震相较S波震相识别难度大，P波振幅较小，与噪声振幅相近，S波振幅与噪声振幅区别明显。在后期的研究中可选用对体波波段响应较灵敏的仪器的监测数据进行研究。
4) 在到时分析过程中可以发现，实际到时与理论到时存在较大差别，在今后的研究中应该选取授时系统准确的台站数据。
致谢: 湖北省地震局吕品姬副研究员为本文研究提供思路和应变资料，GFZ中心汪荣江教授为本研究提供QSSP程序，在此一并表示感谢！
| [1] |
邱泽华, 唐磊, 周龙寿, 等. 四分量钻孔应变台网汶川地震前的观测应变变化[J]. 大地测量与地球动力学, 2009, 29(1): 1-5 (Qiu Zehua, Tang Lei, Zhou Longshou, et al. Observed Strain Changes from 4-Component Borehole Strainmeter Network before 2008 Wenchuan Earthquake[J]. Journal of Geodesy and Geodynamics, 2009, 29(1): 1-5)
(  0) |
| [2] |
邱泽华, 唐磊, 阚宝祥, 等. 钻孔应变台网记录的汶川地震前异常变化[J]. 国际地震动态, 2008(11): 102 (Qiu Zehua, Tang Lei, Kan Baoxiang, et al. Abnormal Changes Prior to the Wenchuan MS8[J]. Recent Developments in World Seismology, 2008(11): 102 DOI:10.3969/j.issn.0253-4975.2008.11.103)
( 0) |
| [3] |
邱泽华, 张宝红, 池顺良, 等. 汶川地震前姑咱台观测的异常应变变化[J]. 中国科学:地球科学, 2010, 40(8): 1-1 039 (Qiu Zehua, Zhang Baohong, Chi Shunliang, et al. Abnormal Strain Changes Observed at Guza before the Wenchuan Earthquake[J]. Science China: Earth Sciences, 2010, 40(8): 1-1 039)
( 0) |
| [4] |
杨选辉, 杨树新, 张彬, 等. 测震数采仪记录的钻孔应变[J]. 中国地震, 2011, 27(1): 65-71 (Yang Xuanhui, Yang Shuxin, Zhang Bin, et al. Strains Recorded by Seismic Acquisition Unit[J]. Earthquake Research in China, 2011, 27(1): 65-71 DOI:10.3969/j.issn.1001-4683.2011.01.007)
( 0) |
| [5] |
Barbour A J, Agnew D C, Wyatt F K. Coseismic Strains on Plate Boundary Observatory Borehole Strainmeters in Southern California[J]. Bulletin of the Seismological Society of America, 2015, 105(1): 431-444 DOI:10.1785/0120140199
( 0) |
| [6] |
Johnston M J S, Borcherdt R D, Linde A T, et al. Continuous Borehole Strain and Pore Pressure in the near Field of the 28 September 2004 M6.0 Parkfield, California, Earthquake: Implications for Nucleation, Fault Response, Earthquake Prediction, and Tremor[J]. Bulletin of the Seismological Society of America, 2006, 96(4B): 56-72 DOI:10.1785/0120050822
( 0) |
| [7] |
万永革. 数字信号处理的MATLAB实现[M]. 北京: 科学出版社, 2012 (Wan Yongge. MATLAB Implementation of Digital Signal Processing[M]. Beijing: Science Press, 2012)
( 0) |
| [8] |
Crotwell H P, Owens T J, Ritsema J. The TauP Toolkit: Flexible Seismic Travel-Time and Ray-Path Utilities[J]. Seismological Research Letters, 1999, 70(2): 154-160 DOI:10.1785/gssrl.70.2.154
( 0) |
| [9] |
Wang R J, Heimann S, Zhang Y, et al. Complete Synthetic Seismograms Based on a Spherical Self-Gravitating Earth Model with an Atmosphere-Ocean-Mantle-Core Structure[J]. Geophysical Journal International, 2017, 210(3): 1 739-1 764 DOI:10.1093/gji/ggx259
( 0) |