地球物理学报  2021, Vol. 64 Issue (6): 1937-1948   PDF    
首次直接观测到与理论预测一致的同震静态应力偏量变化——2016年4月7日山西原平ML4.7地震的钻孔应变观测
石耀霖1, 尹迪1, 任天翔2, 邱泽华3, 池顺良4     
1. 中国科学院计算地球动力学重点实验室, 中国科学院大学, 北京 100049;
2. 中国地质科学院地质研究所, 北京 100037;
3. 中国地震局地壳应力研究所, 北京 100085;
4. 鹤壁市地震局, 河南 鹤壁 458000
摘要:地震会引起同震位移和应力场变化,造成现存各断层上的库仑应力变化量的增加或减小并影响后续地震风险,对于地震预测具有重要意义.现有的大量地震实例中,理论计算的同震位移场与GPS、InSAR等实测的同震位移吻合,间接证实了基于均匀理想弹性假定下可以使用半无限空间Okada公式计算与位移微分相关的同震应力变化.然而,由于地壳并非严格的均匀理想弹性介质,理论计算和直接测量的同震应力变化进行比较仍是一个重要的研究课题.中国对四分量钻孔应变仪的数百个地震同震阶变的研究,发现记录绝大多数阶变不能满足自检条件;美国板块边界观测(PBO)中数百个同震钻孔应变记录同样发现观测幅度与理论预期值相差甚多,可达一两个数量级,这成为国际钻孔应变观测中的一个未解难题.2016年4月7日山西原平发生ML4.7地震,该地震虽然不大,但原平钻孔应力台站距离震中仅19 km,得到了优质的同震阶变记录.本文以此展开研究,结果表明台站观测到的主应力方向和应力偏量大小,均与基于震源机制解的破裂模型的理论预计值很好吻合,这是首次实际观测到与理论计算预测基本一致的同震水平应力偏量变化,为利用库仑应力概念估计后续地震活动性提供了观测基础.然而,观测到的应力变化与理论计算值相差109 Pa的围压,本文也对其可能原因进行了探索,认为在探头、固结水泥和围岩三个可能影响观测的因素中,固结水泥影响可能性最大,并对今后工作提出了建议.
关键词: 地应力      钻孔应变      同震应力变化      库仑应力     
The variation of coseismic static stress deviation consistent with theoretical prediction was observed for the first time —Observation of borehole strain of the Yuanping ML4.7 earthquake in Shanxi on April 7, 2016
SHI YaoLin1, YIN Di1, REN TianXiang2, QIU ZeHua3, CHI ShunLiang4     
1. Key Laboratory of Computational Geodynamics, University of Chinese Academy of Sciences, Beijing 100049, China;
2. Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037, China;
3. Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China;
4. Hebi Seismological Bureau, Henan Hebi 458000, China
Abstract: Earthquakes can cause coseismic displacement and stress changes, increasing or decreasing the Coulomb stress on existing faults and affect subsequent earthquake risks, which is of great significance for earthquake prediction. Observed coseismic displacements by GPS, InSAR, etc. are in good agreement with the calculations, indirectly confirming that using Okada formula for semi-infinite space, calculated coseismic stress changes should also be true, since the uniform half space elastic model works well for displacements. However, due to the earth's crust is not an ideal uniform elastic medium, the comparison between theoretical calculations and directly measurements of coseismic stress changes is still an important research topic. Hundreds records of coseismic off-sets from four-component borehole strain gauges have been obtained in China, but most of them cannot pass the self-check conditions; hundreds of coseismic off-sets in borehole strain records have also been obtained by the PBO (plate boundary observation) observations in the United States, the observed amplitude are significantly different from the theoretical expectations, with huge deviation up to one or two orders of magnitude. It has become an international unsolved problem in borehole strain observation. An earthquake of magnitude ML4.7 occurred in Yuanping, Shanxi on April 7, 2016. Although the earthquake was not large, the Yuanping borehole stress station was only 19 km away from the epicenter, and a high-quality coseismic off-set was recorded. Our study shows that the direction of principal stress and the magnitude of the deviatoric stress observed at the station are in good agreement with the theoretical predictions. This is the first time to show such consistency, which provides an observation basis for the use of the Coulomb stress concept to estimate the subsequent seismic activity. However, there is still a 109 Pa difference in confining pressure between the observed and the theoretical values. We also try to explore its possible causes. There are three affecting factors: the sensor, the cement and the surrounding rocks. It is indicated that, among them, the cement for consolidation is the most possible source for erroneous measurement, and improvement measures for future work are suggested.
Keywords: Crustal stress    Borehole strain    Coseismic stress variation    Coulomb stress    
0 引言

地应力超过岩体强度会造成岩体破坏、断层错断从而发生地震.地震时从震源辐射出地震波造成动态应力扰动,地震波传播经过衰减平息之后,地震断层错断留下了永久变形,造成了同震位移场、应变场和应力场的变化,被称为同震静态变化,这种现象得到了固体地球物理学者们的高度重视(Gomberg et al., 1998; Kilb et al., 2002; Ma et al., 2005).由于同震应力变化,可能会促使其他断层上的某些地震更加容易发生、甚至触发新的地震;也可能形成“应力影”,即改变的应力状态降低了某些断层的地震危险性.这对于研究大地震对后续地震活动性的影响,开展物理地震预报研究具有十分重要的意义.

地震断层造成怎样的位移、应变和应力变化,早在20世纪60年代开始就有了不少研究(Chinnery, 1961, 1963; 陈运泰等, 1975; 黄福明和王廷韫, 1980),特别是Okada(1980, 1985, 1992; Okada and Yamamoto, 1991)的一系列工作已经形成了一个十分成熟的研究体系,理论预测的位移场也被大量地震后GPS、InSAR等观测所证实(Tsuji et al., 1995; Simons et al., 2002).因此,人们广泛开展了大震的同震应力变化的计算,以及应力变化对断层破裂的影响,提出库仑应力变化量ΔCFS的概念:

其中τn为断层面上剪应力,σn为断层面上正应力,P为孔隙流体压力,μ为断层摩擦系数(King et al., 1994; Cocco and Rice, 2002).断层在原有错动趋势方向上的剪应力增加会增大地震危险性,减小则相反;断层上的正应力如果是压应力减小(张应力增加)会使断层更加容易错动,增加地震危险性,压性正应力增加则相反;为负值的孔隙流体压力的增加会造成断层上正应力有效压力的减小,增加地震危险性,孔隙流体压力的减少则相反.这些影响都反映在库仑应力变化量上,即ΔCFS变化为正值表示更危险,ΔCFS为负值表示更安全.

库仑应力变化的概念被普遍认可和广泛应用,国内外已经发表了很多相关的文章用于参考借鉴.理论预测的位移场已被GPS、InSAR等观测资料证实,如果把地壳视为均匀理想弹性介质,应变场和应力场变化是根据位移的空间微分计算出来的,因此库仑应力的存在也得到了间接的证实.然而,如同山的高度虽然决定了山体的平均坡度,但山路上具体的每一点坡度作为高度的空间微分仍会有大小起伏偏离平均坡度一样,地壳岩石广泛存在垂向和横向不均匀性,并充斥着断层裂隙和孔隙流体.因此,把直接观测的同震应力变化和基于均匀弹性介质简化的理论计算值进行对比研究分析,仍是一项重要的研究课题.毕竟,迄今只有间接证据,而缺失同震静态应力变化直接观察的确凿证据,是地球物理研究者们共同面临的难题.

应力是一个张量,完全描述三维应力状态需要六个独立分量,平面应力状态则需要三个独立分量.人们通过垂直钻孔内孔壁三个独立方向上各自变形的测量,依据圆孔在均匀应力场内弹性变形的规律,计算水平面内平面应力状态(潘立宙,1980).应力张量可以分解为球应力张量和偏应力张量两部分的和,前者影响弹性介质体积的涨缩,后者影响其形状的变化.作为平面应力张量有三个独立分量,因此至少要在三个独立方位上进行测量得到孔径变化.我国首次提出和实行了安置45°角的4个元件测量的方法(欧阳祖熙和张宗润, 1988),根据应力偏张量对圆孔变形影响规律,方位角互相垂直的一号和三号元件测值变化之和应该等于二号和四号元件测值变化之和,即所谓1+3=2+4的准则,判断记录是否可靠,进行自检.国外最开始采用三分量,之后也采纳了我国四分量的做法(Gladwin and Michael, 1984).

20世纪60年代至70年代,我国分量式钻孔应变仪灵敏度尚不足以观测到固体潮,20世纪80年代,我国自主研制了RZB分量式钻孔应变仪相继投入使用.“十五”计划期间,在中国地震局“数字地震观测网络”计划支持下,全国各地布设了数十套YRY-4型四分量钻孔应变仪,多年实际观测表明它们取得了较为可靠的观测数据(池顺良等,2007),可以记录到高频地震信息及低频的固体潮、年周期变化和长年变化的信息(邱泽华和池顺良,2013),并满足自检条件(吴立恒等, 2012).在营口地震台同时安装有伸缩仪,虽然与钻孔应变仪基线长度和安装方式不同,但两种仪器观测结果具有可比性(池顺良等,2008).

在多年多台大量观测中,四分量钻孔应变观测数据在从地震到长年变化超宽频段内均满足1+3=2+4的自检条件,唯独在记录同震阶变时,大部分观测数据却出现了相关性很差的局面.全国34个钻孔四分量应变仪记录的汶川地震同震应变阶跃变化(唐磊和邱泽华, 2011),昌平台站钻孔四分量应变仪观测到的1998年和1999年张北2次地震阶跃变化(邱泽华和石耀霖,2003),不仅1+3和2+4不相等,而且观测应变阶跃值比断层错断理论计算值有的高出一两个数量级,主应力方向与理论值也没有可比性.近年我国对一些钻孔应变资料加强了分析和利用(Kong et al., 2018; Gong et al., 2019),但展示的不同地震、不同台站的阶变记录仍不能满足自检条件.

国外虽然很早也注意到了阶跃现象,但当时没能做出定量的分析(Gwyther et al., 1992).最近Barbour等(2015)对9个PBO钻孔应变台站的34个地震记录中记到的311个阶变进行了分析,发现阶变幅度与理论预期的幅度量级相差甚多,变化范围在低一个数量级到高三个数量级之间.而且观测应变绝对阶变量与反映地震能量密度的峰值动态应变(PDS)有关,他们觉得“这种效应的确切机制尚不清楚,并且很难从手头的数据中分析出来”.

关于库仑应力的理论计算蓬勃开展,而且GPS对位移场的观测也证明了其结果不容置疑,应力观测提供同震应力变化的直接观测验证成为一个重要而紧迫的科研课题,也是世界上钻孔应变测量观测中共同面临的未解难题.本文将利用原平地应力站对原平地震的同震钻孔应变观测数据,进行探索、分析和讨论.

1 原平地震

据中国地震台网测定,2016年4月7日山西原平发生ML4.7地震.山西地震台网测定震源深度为11 km,中国地震台网给定震源深度为16 km.李斌等(2019)利用Hpyo2000定位方法对地震重新定位得到震中位置为38.86°N, 112.92°E,震源深度为14.7 km,发生在忻定平原内部,五台山北麓断裂与五台山西麓断裂附近(吴昊昱等, 2018)(图 1).同时李斌等(2019)吴昊昱等(2018)采用不同方法反演得到原平地震的震源机制解,具体见表 1,根据实际地质断层考察,五台山北麓断裂走向N50°~70°E,倾角60°~70°(王秀文等,2018), 与震源机制解节面Ⅱ更接近,因此选取节面Ⅱ为原平地震的实际震源机制解参数.

图 1 山西原平地震震源机制解 其中蓝色震源球为Dreger与Snoke方法得到的震源机制解,红色震源球为CAP方法得到的震源机制解. Fig. 1 Focal mechanism solution of Shanxi Yuanping earthquake The blue beach balls are the focal mechanism solution obtained by the Dreger and Snoke method, and the red beach ball is the focal mechanism solution obtained by the CAP method.
表 1 原平ML4.7地震震源机制参数 Table 1 Parameters of the Yuanping ML4.7 earthquake focal mechanism solution
2 原平地震造成的同震应力变化——理论计算

理论计算过程中,我们将根据震级推断地震震源尺度和错动量,根据震源机制,推断断层错动模型,计算原平地震在原平地应力台站造成的应力变化,然后与台站实际观测到的应力变化进行比较和分析.

推断震源尺度和错动量需要了解地震的矩震级,虽然地震目录提供的是近震震级ML,但在地震震级转换关系方面前人做了许多工作, 李莹甄等(2014)进行了总结,给出了矩震级与各种震级标度之间的对应关系.原平地震震级为ML4.2,根据我国常规震级标度与矩震级之间的对应关系,该区MLMSMW,根据震级与破裂长度及破裂面积的经验公式(龙峰等,2006)计算断层的破裂长度L及破裂面积A.

(1)

式中,G为剪切模量,slip为滑动量.

参照前人对山西原平P波速度结构的研究成果,我们取地壳P波平均速度为6.3 km·s-1,略低于华北地壳平均速度(王长在等,2018),该区泊松比取为0.25(武岩等,2018).根据纵横波速度与拉梅常数关系计算得到研究区域杨氏模量:E=73 GPa,剪切模量G=29 GPa.最终计算出断层破裂长度为1.60 km,破裂面积为1.17 km2,破裂宽度为0.73 km.根据地震矩M0和矩震级MW的关系(陈培善和白彤霞,1991),求得M0=2.493×1015N·m,与李斌等(2019)结果一致.根据M0=GAS,其中A为断层面积,S为平均错动量(Kanamori and Anderson, 1975),计算得到断层错动量为0.073 m.

根据地震目录给出的震中位置和震源深度、震源机制解得到的断层走向、倾角和滑动方向,以及推导得到的断层尺度和错动大小,利用汪荣江开源程序(edcmp/edgrn)计算三种震源机制下(李斌等,2019;吴昊昱等,2016)原平地震对原平台造成的应力应变张量及最大最小主应力,具体理论计算结果见表 2,规定压应力为负.同时计算了三种不同震源深度和机制情况下最大最小主应力及方位角均值和标准偏差,其中σ1的标准偏差为11.8 Pa(相对误差11.7%),σ2的标准偏差为2.44 Pa(相对误差17.3%),θ的标准偏差为12.4°,这种量级的偏差在地球物理观测中一般是可以接受的.鉴于山西省从地表到35~40 km Moho边界的地壳速度变化不大(宋美琴等,2012),而震源深度只有14 km,因此计算中没有分层,退化为弹性半无限空间模型.

表 2 原平台站处三种震源机制下的应力分量及最大最小主应力理论解 Table 2 Theoretical solution of stress components and maximum and minimum principal stresses at Yuanping station for the three focal mechanisms
3 原平地震造成的同震应力变化——实际观测

原平地应力观测台站(38.72°N, 112.8°E)位于原平地震震中(38.86°N, 112.92°E)西南方向19 km.该台站高程857 m,处于山西盆地边缘向山区过渡地带,岩性为花岗岩,岩石完整.2015年10月13日台上安装了YRY-4钻孔应变计,分辨率为5×10-11,采样频率为10 Hz.孔深45 m,探头安装在接近孔底处,孔内水位距离井口40 m.虽然该套仪器不带温度计,但根据其他台站经验这种深度探头处日温变化小于0.001 ℃.仪器从安装至今数据完整率高于99%,固体潮记录清晰,M2波潮汐因子相对精度达到0.003.地震时记录到了地震波动,不仅波动满足自检条件(图 2),而且阶跃也满足1+3=2+4的自检要求(图 3).这是迄今钻孔应变仪记录中不多见的优质记录.

图 2 山西原平台高通滤波后自检曲线图 Fig. 2 Self-test curve after high-pass filtering in Yuanping platform, Shanxi
图 3 山西原平四分量钻孔应变台站观测到的2016年4月7日地震记录(10 Hz采样率) Fig. 3 The earthquake record of April 7, 2016 (10 Hz sampling rate) observed by the four-component borehole strain station in Yuanping, Shanxi

原平地应力观测台站的分量式钻孔应变仪4个径向位移传感器方位角(相对磁南北)分别为:Ⅰ:-56°;Ⅱ:-11°;Ⅲ:34°;Ⅳ:79°. 经过磁偏角校正(-4.5°)后原平台四路元件的方位角(相对磁南北)为:Ⅰ:-60.5°;Ⅱ:-15.5°;Ⅲ:29.5°;Ⅳ:74.5°.计算过程中取东向为x轴,角度为0°;北向为y轴,角度为90°,对四路元件进行坐标变换得到四路元件与x轴的夹角如图 4b所示.

图 4 原平台四路元件相对磁南北的方位角,其中黑色虚线为磁偏角校正前数值,黑色实线为校正后数值;(b) 原平台四路元件校正后用于计算中相对x轴的夹角 Fig. 4 (a) The azimuth of the four component sensor relative to the magnetic north in Yuanping platform. The black dotted line is the value before magnetic declination correction, and the black solid line is the value after correction. (b) The angle of the four components relative to the x-axis after correction used for calculation in Yuanping platform

其中Ⅰ, Ⅱ, Ⅲ, Ⅳ号元件阶跃值分别为-9.011×10-10,-0.912×10-10,-9.396×10-10,-16.662×10-10,阶跃满足1+3=2+4的自检要求(图 3).按照石耀霖和范桃园(2000),知道钻孔附近岩石力学性质时,每三个元件一组,可以计算共四组水平应变分量,及最小二乘平均值,并计算它们的标准方差,结果如表 3.其中根据钻孔应变仪的四个元件的观测值计算最小二乘法水平应变张量的三个分量的公式如下,其余各组公式不再一一列出:

(2)

表 3 四种原件分组下的水平应变分量及最小二乘均值 Table 3 The horizontal strain component and the least square mean value of the four original groups

其中K4O4K3O3K2O2K1O1分别代表Ⅰ、Ⅱ、Ⅲ、Ⅳ四路元件记录到的同震应变阶跃值.

根据应变张量计算最大主应变以及方位角,根据弹性介质胡克定律由应变计算应力,具体计算结果见表 4,规定压应力为负,其中σ1σ2的标准偏差分别为2.91和2.90 Pa.

表 4 由观测四路应变计算得到的平面应变/应力分量及平面主应变/主应力 Table 4 Strain/stress components and principal strains/stresses calculated from the observation of four-component strain data
4 讨论

原平地震震中靠近五台山北麓断裂(图 1),五台山北麓断裂走向N70°E, 与CAP方法反演得到的震源机制解250°更接近,震源深度14.7 km结果更为可靠,因此在对比分析时Okada解选用震源机制解(250°/63°/-120°)下震源深度14.7 km计算得到的理论结果与观测结果进行对比分析,具体Okada理论值与观测值计算结果对比见表 5.

表 5 原平台平面应力变化的观测值与理论解对比 Table 5 Comparison of observational and theoretical stresses at Yuanping Station

表观上,观测结果与Okada理论结果似乎存在显著差异:理论计算值σ1为张应力101 Pa,σ2为压应力-14 Pa,最大主应力与x轴夹角为-6.7°;实际观测值σ1为-8 Pa,σ2为-123 Pa,σ1x轴夹角为16.2°.邱泽华(2017)的计算结果也发现了类似问题.但应该注意到,不仅观测值和理论值的σ1的方位角接近,而且最大剪应力(σ1-σ2)/2的观测值为57.5 Pa, 理论值为57.5 Pa,相差接近于0.实际上,注意到(σ1+σ2)/2的观测值(-65.5 Pa)与理论值(+43.5 Pa)的差值约为-109 Pa,将观测结果加上109 Pa的均匀围压后,σ1为101 Pa,σ2为-14 Pa,与Okada理论解在数值上基本一致,角度差异为23°(图 5).考虑到国际上PBO对同震应变300多次观测结果,90%以上的观测值与理论值在数值上存在一个数量级以上的差异,而方向上未能提供比较(Barbour et al., 2015).也考虑到震源深度和机制本身测量存在很大的误差,例如本次地震三个震源深度和机制解最大压应力的相互差异可以达到20 Pa、最大剪应力相互差异可以达到12 Pa,主应力方位角相互差异可以达到22°(表 2).我们认为本次测量到的同震应力偏量变化已经能较好的与理论计算预期值吻合.唯一的欠缺是,四个观测探头还显示受到理论值预期之外均匀挤压,折算成水平均匀压应力相当于109 Pa的压应力.

图 5 原平台平应力张量变化观测值与理论预测值对比示意图 (a) 最大最小主应力大小和方向理论值与观测值直接对比,其中蓝色箭头为理论值,红色箭头为观测值;(b) 理论值与校正后观测值对比;(c) 观测值校正原理示意图,实际观测值来自两部分贡献:观测值可视为地震产生的理论预测张量值与109 Pa均匀围压二者之和. Fig. 5 A schematic show for comparison between the observed and the theoretical predicted changes of the horizontal stress tensor at Yuanping Station (a) The magnitude and direction of the maximum and minimum principal stresses of the original observed raw data compared with the theoretical predictions. The blue arrows are the theoretical predicted values, and the red arrows are the observed values. (b) Comparison of the theoretical values with the corrected observation value. (c) Schematic diagram of observation correction principle. The actual observations come from two contributions: the observed values can be regarded as the sum of the theoretical predicted stress tensor due to the earthquake and a uniform confining pressure of 109 Pa.

为什么以前多次大地震同震应力阶变记录不满足1+3=2+4自检条件?为什么原平地震的阶变却能记录到同震应力变化的应力偏量阶跃,但平面应力部分仍然有109 Pa的差异?这些都是值得深入探讨的问题.

远震的地震波造成的动态应力变化,会对接收点产生多种影响.虽然不像大地震在接近震源地区造成山崩地裂、房倒屋塌的效果,但是也有可能触发地震(Gomberg et al., 1998; Kilb et al., 2002),除了动态库仑应力自身作用外,触发地震的可能机制之一是振动对岩体和断裂带内物质的力学性质产生影响.Taira等(2009)认为远震波动改变了圣安德烈斯断层物质强度,他们利用对圣安德烈斯断层进行的20年(1987—2008)观测资料的研究,提出了可以通过监测断层强度来实现地震现场实际监测的方法.通过对地震散射体的特性随时间变化(可能反映了应力引起的流体运移)和对重复地震序列特征的系统性随时间变化的分析,发现了许多8级以上大地震,包括遥远的2004年苏门答腊—安达曼地震产生的断层强度变化十分重要.地震也可以引起岩石和沉积物的水文力学特性很大变化,因而对很远的地震仍然敏感(Moyer et al., 2018),引起地震后地下水位和水压、水温度、水的化学成分、对固体潮汐反应的变化等(Wang and Manga, 2010, 2014).

本文分析的YRY-4型分量钻孔应变探头直径108 mm,用水泥砂浆固结在约130 mm直径的花岗岩钻孔中44.7 m深处.一般水泥固结后收缩而在探头和岩体间产生空隙,国内外一般需要采用膨胀水泥,但过强膨胀会使探头几个月甚至更长时间仍有明显蠕变.原平台在固结水泥和固结技术上做了精心安排,台站比较快的进入了正常记录观测.

对于以往许多同震阶变不能满足1+3=2+4的物理原因,需要对影响观测记录的三个主要变形传递介质,岩石、胶合水泥和机械探头分别探讨.我国钻孔应变台站都是建立在较完整的岩石上,而且长期变化、潮汐变化,甚至地震波动态记录均能较好满足自检要求,但我国在2008年汶川地震和2011年日本大地震后,和国外一样,阶跃的记录在量值上与理论值相差甚远,且基本都不满足自检条件.首先考虑、也最易于改进的是探头,虽然探头在缓慢变化记录中表现良好,但地震波反复震荡作用下,探头机械部分由于加工精度不够可能发生微小位移而掩盖了真实的同震位移,因此随后更加注意探头的机械加工和安装精度,原平台就是2015年10月投入工作的新建台站之一.预期它应该有更强的抗震动干扰性,但迄今没有大震检验实例,这次震中距很小的ML4.7地震提供了一次检验机会,而且记录的阶跃确实能够满足自检条件.

关于不满足自检条件的其他可能性,董培育等(2015)进行了探讨.他们通过模拟分量式钻孔探头金属外壁和固结水泥环之间产生张裂隙或剪切滑动,发现若钻孔内的探头圆筒和水泥环之间发生剪切滑动、特别是发生张裂会影响测量结果,使得测量结果与1+3=2+4自检条件出现较大差异,严重影响到同震阶变的观测.步玉环等(2011)实验研究显示,膨胀水泥中晶格膨胀剂的加入会使水泥-套管界面胶结强度明显下降,且界面胶结强度抵抗温度和压力变化的能力也变差,他们认为研究一种既能实现膨胀功能、又能改善界面胶结力和水泥石力学性能的膨胀剂具有重要的意义.

为什么国内过去多个大地震多个台站的应变记录阶跃总是不能满足自检条件,暂且不考虑方向,国际上实测阶变数值与理论计算数值为什么仅在大小上就存在数量级的差异?还有一个可能性是地震波会对固结水泥产生影响.应变探头尺度在分米量级,而测量的微小应变仅为10-9~10-10,即位移在10-10~10-11m的量级,大地震面波造成应变记录的幅度变化幅值可达数百纳应变(10-7)量级,在数十纳应变幅度以上的振动持续时间可超过十余分钟甚至几十分钟,相当于应力振动最大幅度为1~0.1 kPa.在这样长时间的反复震荡下,固结探头用的膨胀水泥与探头胶结即使略有瑕疵,只要发生纳米量级的不连续位移,就会造成不满足自检条件.而原平地震震中距只有19 km,面波尚未发育,记录到的S波最大振幅虽然也可达到近百纳应变,但大振幅波动只持续了几秒钟,水泥胶结性能没有受到显著影响.而原平台站花岗岩体完整、探头技术也不断优化,这可能是2015年新建的原平台应变波记录能够满足1+3=2+4自检条件、记录的应力偏量与理论预测一致的主要原因.

为什么原平台站的应变波动记录与理论值比较,仍然会多出109 Pa的围压应力呢?1+3=2+4的自检条件,仅是检测测量区域应力场的应力偏量是否出现了误差,对于均匀围压部分并不能起到此辨别作用.探头测量的面应变结果,可能是反映了岩体内的应力,但也可能只是固结水泥环内水泥膨胀施加于探头的作用,仅从探头读数上是无法辨识的.水泥是一种独特的工程材料,微观研究显示了其复杂的内部结构(Diamond,2004徐超等,2009),如图 6所示.水泥浆中的水泥颗粒,颗粒外层与水接触发生水化反应与水结合形成水化硅酸钙(C-S-H)壳,反应逐步由水泥颗粒外部向内部推进,需要很长时间.平原地震波振动可能加速了残余孔隙水的渗流及与未反应的水泥颗粒的接触和反应,只要水泥环厚度增大了0.1 nm, 就足以产生观测到的均匀围压增大.

图 6 室温下放置100天后的固化的水泥微观结构.水泥颗粒与水发生水合反应使水泥得以凝固,水泥颗粒一般在数十微米大小,较大颗粒的核心部分往往还未发生水合反应(据Diamond,2004) Fig. 6 Microstructure of cured cement after 100 days at room temperature. The hydration reaction between cement particles and water makes the cement solidify. Cement particles are generally tens of microns in size, and the core of larger particles may not yet have been hydrated (According to Diamond, 2004)

因此,今后应该专门研究如何增强探头套管与固结水泥的胶结力,除了水泥本身的改进增强与套管表面的化学胶着力外,探头表面喷砂打毛工艺增强水泥与探头套管接触面上的机械咬合力也值得尝试.水泥的颗粒大小、配方性质、灌浆工艺等都需要进一步细致的探索研究.地应力探头用水泥固结在钻孔中要工作几十年,无法提起更换.因此,这些研究应该在实验室条件下提前部署和安排,把科研成果落实到新一代的地应力台站建设中.

围压增加109 Pa也还可能有其他解释.一种可能机制是,振动引起地下水位变化.压力半周期造成孔隙水被挤出和地下水位的升高,但是张力半周期由于固体颗粒的迟滞反应,固体骨架没有完全复原,水也没有能够全部恢复到原来的孔隙中.多次地震波动的反复作用下,或者是钻孔附近区域地下水位或者是钻孔中的水位升高了几个厘米,都有可能引起探头测到围压变化.Lu和Wen(2018)认为加州钻孔应变仪记录到的一些长周期变化和孔隙弹性有关,对阶变研究也许有借鉴作用.

总之,我国和世界钻孔应变记录的阶变与理论计算的同震应力变化普遍有明显差异(Barbour et al., 2015),是一个世界性的难题.比较有共识的是,动态的地震波对岩石、固结水泥或探头造成了一些局部影响,造成了钻孔应变记录阶变的复杂变化,但认识未能细化.本文研究认为,我国的探头对长期变化、潮汐变化和地震波记录都满足自检条件,这次的震例又显示对阶变也可以做到满足自检条件,说明了我国的探头是可靠的.我国的地应力台站钻孔一般都精心选择在岩体比较完整的场所,例如原平台站在花岗岩上,虽然不排除岩石在振动下发生非弹性变形的可能性,但可能性较小,最大的可能是在较远岩体薄弱的水泥固化环状层中发生了非弹性变形.虽然目前的资料尚不足以得到具体肯定的答案,但应该是今后考察的重点方向.

钻孔应变仪记录的同震阶变经常可以远大于理论估计值一两个数量级甚至更高,使得许多理论预测不会记到阶变的台站经常也会记到阶变(Barbour et al., 2015), 这对应变记录虽然是一种干扰,但反过来,对我们了解远震动态触发却提供了某种启示.如果地震波的持续振动可以对膨胀水泥层造成非弹性变形,地震波对薄弱的活动断裂带也有可能造成同样的影响.这种影响比理论估计值可以大两个数量级以上,这有助于我们进一步探索大地震能在遥远距离触发地震活动(Gomberg and Johnson, 2005)的物理原因.甚至人们可以转变概念,不仅以测量地应力为目标,在完整的岩体里钻孔安装应变探头;也可以今后特地在断裂带内安装应变探头,其目标并不在于准确记录地应力,而是为了了解在远方大地震应力波传来时,断裂带里的记录与完整岩体里的记录有什么差别,探索应力波动下断层物质性质的变化和应力状态的变化有什么特点.

5 结论

钻孔应变仪在地震时记录的同震阶变普遍与理论预测值不吻合,甚至相差一两个数量级,成为世界钻孔应变研究中未解的难题.本文利用山西原平钻孔应力台站得到的其北东19 km发生的2018年4月7日原平ML4.7地震同震阶变的优质记录,计算了原平台实际观测到的同震应力变化,并与基于震源机制的地震错动模型和半无限空间Okada公式,计算了理论预期的同震应力变化.观测结果与理论预期比较发现,原平台的实际观测值可以视为理论预测值加上109 Pa的均匀围压产生的.理论预测的台站处最大主应力方位角为接近北偏东6.7°,实测结果为北偏西16.2°;差应力(σ1-σ2)理论预测值为115 Pa,实测值亦为115 Pa.记录到的主应力方位和应力偏量大小,与理论值大体吻合.可以认为,原平台站成功的记录到了原平地震产生的应力偏量变化.

实测值与理论值相差109 Pa的均匀围压.我们认为这并不反映岩体中真实的应力变化,很可能是把元件固结在钻孔中使用的固结水泥,在地震波动振动下,固结水泥颗粒内部未反应的部分与孔隙水接触并发生了水化反应,水泥膨胀而造成向内挤压元件外壁,发生了不到0.2 nm的微小位移,产生了这相当于109 Pa压应力的微小变化.固结探头时膨胀水泥可以造成数百千帕到1 MPa的压应力,这里振动后发生0.109 kPa变化并不算大,但对于精密的同震应力变化测量仍然是一个突出的干扰因素.另外,对未来远处大震动态应力高达近千帕并持续十余分钟的条件下,仪器是否还能满足自检的准确记录同震应力阶变,尚有待于实践的考验.

地应力测量方法虽然是从国外引进的,但是以地震预报为目标的长期地应力观测是李四光先生在1966年邢台地震后首先投入实践的,是我国有特色的科学创新.这种长期观测与工程部门的观测在灵敏度和稳定性等方面的要求都有很大的不同.早期的地应力观测曾经暴露出许多问题,但近二十多年观测方法日趋成熟,在长期观测、固体潮观测和地震波动观测中都能达到自检要求.唯独在地震同震阶变的观测上仍然存在问题.这也是世界上钻孔应变阶变记录共同遇到的成因尚不完全明朗的问题.本文通过选择山西原平近震的应变波记录分析,表明钻孔应变仪记录到了和理论预测值一致的同震应力偏量变化.这一方面首次通过直接观测验证了地震位错产生同震应力变化计算的可靠性,为今后通过库仑应力计算估计地震危险性提供了观测基础;另一方面也证明了现在钻孔应变测量探头的可靠性,为今后有效实测同震应力阶跃变化积累了经验.由于近几年最新的国产探头是可靠的,阶变观测的异常可以确定是固结水泥或孔外岩石在动态地震波作用下产生的效应,其中可能性最大是由于比岩石脆弱得多的固结水泥中产生了非弹性变形.

本文工作也揭示了,水泥固结探头的技术方法虽然取得了总体的成功,保证了应力长期变化、潮汐变化和地震波动的记录,但是在大地震长达十几分钟到几十分钟的面波重复振动作用下,可能会发生微小的纳米量级的非弹性变形,这种变形虽然不会影响前面谈到的观测,但足以影响同震阶变数据的准确获取.今后在水泥的性质和使用上,在钻孔应变探头的加工工艺上,需要作进一步的细致工作和改进,使我国在地应力预测地震的研究上,继续保持走在世界的前列.而且如果能研究在真空环境中和严酷的温度条件下在数米浅孔内固结探头的方法,也许可以率先在月球上开展月球固体潮、月应力和月震动态及静态应力测量.

致谢  我们感谢三位匿名审稿人的中肯意见.感谢山西省原平钻孔应变台站在提供数据和资料方面的大力支持.
References
Barbour A J, Agnew D C, Wyatt F K. 2015. Coseismic strains on plate boundary observatory borehole strainmeters in Southern California. Bulletin of the Seismological Society of America, 105(1): 431-444. DOI:10.1785/0120140199
Bu Y H, Liu H J, Song W Y. 2011. Experimental study on effects of lattice expansive agent on cement-casing interfacial bonding quality. Acta Petrolei Sinica (in Chinese), 32(6): 1067-1071.
Chen P S, Bai T X. 1991. Quantification relations between the source parameters. Acta Seismologica Sinica (in Chinese), 13(4): 401-411. DOI:10.1007/BF02650539
Chen Y T, Lin B H, Lin Z Y, et al. 1975. The focal mechanism of the 1966 Hsingtai (邢台) earthquake as inferred from the ground deformation observations. Chinese Journal of Sinica (in Chinese), 18(3): 164-182.
Chi S L, Wu H L, Luo M J. 2007. Discussion on strain tidal factor separation and anisotropy-Analysis of first data of borehole component strain-meter of China's digital seismological observational networks. Progress in Geophysics (in Chinese), 22(6): 1746-1753.
Chi S L, Gao R, Wang L, et al. 2008. Data comparison of extensometer and 4-component borehole strainmeter at Yingkou seismological station. Recent Developments in World Seismology (in Chinese), (3): 34-37.
Chinnery M A. 1961. The deformation of the ground around surface faults. Bulletin of the Seismological Society of America, 51(3): 355-372.
Chinnery M A. 1963. The stress changes that accompany strike-slip faulting. Bulletin of the Seismological Society of America, 53(5): 921-932.
Cocco M, Rice J R. 2002. Pore pressure and poroelasticity effects in Coulomb stress analysis of earthquake interactions. Journal of Geophysical Research: Solid Earth, 107(B2): ESE 2-1-ESE 2-17.
Diamond S. 2004. The microstructure of cement paste and concrete-a visual primer. Cement and Concrete Composites, 26(8): 919-933. DOI:10.1016/j.cemconcomp.2004.02.028
Dong P Y, Ren T X, Yang S H, et al. 2015. A problem and explanation for borehole strain meter records of co-seismic strain steps. Journal of Geomechanics (in Chinese), 21(3): 359-370.
Gladwin M T. 1984. High-precision multicomponent borehole deformation monitoring. Review of Scientific Instruments, 55(12): 2011-2016. DOI:10.1063/1.1137704
Gomberg J, Beeler N M, Blanpied M L, et al. 1998. Earthquake triggering by transient and static deformations. Journal of Geophysical Research: Solid Earth, 103(B10): 24411-24426. DOI:10.1029/98JB01125
Gomberg J, Johnson P. 2005. Dynamic triggering of earthquakes. Nature, 437(7060): 830. DOI:10.1038/437830a
Gong Z, Jing Y, Li H B, et al. 2019. Static-dynamic strain response to the 2016 M6.2 Hutubi earthquake (eastern Tien Shan, NW China) recorded in a borehole strainmeter network. Journal of Asian Earth Sciences, 183: 103958. DOI:10.1016/j.jseaes.2019.103958
Gwyther R L, Gladwin M T, Hart R H G. 1992. A shear-strain anomaly following the Loma Prieta earthquake. Nature, 356(6365): 142-144. DOI:10.1038/356142a0
Huang F M, Wang T Y. 1980. The stress field of a dislocating inclined fault. Acta Seismologica Sinica (in Chinese), 2(1): 1-20.
Kanamori H, Anderson D L. 1975. Theoretical basis of some empirical relations in seismology. Bulletin of the Seismological Society of America, 65(5): 1073-1095.
Kilb D, Gomberg J, Bodin P. 2002. Aftershock triggering by complete Coulomb stress changes. Journal of Geophysical Research: Solid Earth, 107(B4): ESE 2-1-ESE 2-14. DOI:10.1029/2001JB000202
King G C P, Stein R S, Lin J. 1994. Static stress changes and the triggering of earthquakes. Bulletin of the Seismological Society of America, 84(3): 935-953.
Kong X Y, Su K Z, Yukio F, et al. 2018. A detection method of earthquake precursory anomalies using the four-component borehole strainmeter. Open Journal of Earthquake Research, 7(2): 124-140. DOI:10.4236/ojer.2018.72008
Li B, Li Z H, Song M Q, et al. 2019. The moment tensor inversion of Yuanping ML4.7 earthquake in Shanxi province. Journal of Taiyuan University of Technology (in Chinese), 50(5): 598-605.
Li Y Z, Yin N, Li X H. 2014. Review of the conversional relationship for different magnitude scales. China Earthquake Engineering Journal (in Chinese), 36(1): 80-87.
Long F, Wen Z X, Xu X W. 2006. Empirical relationships between magnitude and rupture length, and rupture area, for seismogenic active faults in north China. Seismology and Geology (in Chinese), 28(4): 511-535.
Lu Z, Wen L X. 2018. Strong hydro-related localized long-period crustal deformation observed in the Plate Boundary Observatory borehole strainmeters. Geophysical Research Letters, 45(12): 12856-12865. DOI:10.1029/2018GL080856
Ma K F, Chan C H, Stein R S. 2005. Response of seismicity to Coulomb stress triggers and shadows of the 1999 MW=7.6 Chi-Chi, Taiwan, earthquake. Journal of Geophysical Research: Solid Earth, 110(B5): B05S19. DOI:10.1029/2004JB003389
Moyer P A, Boettcher M S, Mcguire J J, et al. 2018. Spatial and temporal variations in earthquake stress drop on Gofar transform fault, East Pacific rise: implications for fault strength. Journal of Geophysical Research: Solid Earth, 123(9): 7722-7740. DOI:10.1029/2018JB015942
Okada Y. 1980. Theoretical strain seismogram and its applications. Bulletin of the Seismological Society of America, 55: 101-168.
Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 75(4): 1135-1154.
Okada Y, Yamamoto E. 1991. Dyke intrusion model for the 1989 seismovolcanic activity off Ito, central Japan. Journal of Geophysical Research: Solid Earth, 96(B6): 10361-10376. DOI:10.1029/91JB00427
Okada Y. 1992. Internal deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 92(2): 1018-1040.
Ouyang Z X, Zhang Z R. 1988. Studies of method for coupling straingauge with borehole wall. //Crustal Tectonic and Crustal Stress (2) (in Chinese). Beijing: Seismological Press. http://www.oalib.com/references/18990729.
Qiu Z H, Chi S L. 2013. Shear strains of p wave observed with YRY-4 borehole strainmeter. Earthquake (in Chinese), 33(4): 64-70.
Qiu Z H. 2017. Theory and Application of Borehole Strain Observation (in Chinese). Beijing: Seismological Press.
Shi Y L, Fan T Y. 2000. Borehole in situ calibration of stress sensors and calculation of variation of stress field during long term observation. Earthquake (in Chinese), 20(2): 101-106.
Simons M, Fialko Y, Rivera L. 2002. Coseismic deformation from the 1999 MW7.1 Hector mine, California, earthquake as inferred from InSAR and GPS observations. Bulletin of the Seismological Society of America, 92(4): 1390-1402. DOI:10.1785/0120000933
Song M Q, Zheng Y, Ge C, et al. 2012. Relocation of small to moderate earthquakes in Shanxi Province and its relation to the seismogenic structures. Chinese Journal of Geophysics (in Chinese), 55(2): 513-525. DOI:10.6038/j.issn.0001-5733.2012.02.014
Taira T, Silver P G, Niu F L, et al. 2009. Remote triggering of fault-strength changes on the San Andreas fault at Parkfield. Nature, 461(7264): 636-639. DOI:10.1038/nature08395
Tang L, Qiu Z H. 2011. Coseismic strain steps of the Wenchuan earthquake observed by borehole four-component strain gauge. Collected Works of Crustal Tectonics and Crustal Stress (in Chinese): 114-124.
Tsuji H, Hatanaka Y, Sagiya T, et al. 1995. Coseismic crustal deformation from the 1994 Hokkaido-Toho-Oki Earthquake Monitored by a nationwide continuous GPS array in Japan. Geophysical Research Letters, 22(13): 1669-1672. DOI:10.1029/95GL01659
Wang C Y, Manga M. 2010. Hydrologic responses to earthquakes and a general metric. Geofluids, 10(1-2): 206-216. DOI:10.1111/j.1468-8123.2009.00270.x
Wang C Y, Manga M. 2014. Earthquakes and Water. //Meyers R ed. Encyclopedia of Complexity and Systems Science. New York, NY: Springer. doi: 10.1007/978-3-642-27737-5_606-1.
Wang C Z, Wu J P, Yang T, et al. 2018. Crustal structure beneath the Taiyuan Basin and adjacent areas revealed by double-difference tomography. Chinese Journal of Geophysics (in Chinese), 61(3): 963-974. DOI:10.6038/cjg2018L0114
Wang X W, Zhao H M, Jin Y K. 2018. Deformation anomalies of Wutaishan fault in Shanxi. Recent Developments in World Seismology (in Chinese), (4): 74-79.
Wu H Y, Li L, Wang X. 2018. Discussion on focal mechanism solutions and seismogenic structure of the M4.2 Yuanping earthquake on 7, April, 2016 in Shanxi province. Bulletin of Science and Technology (in Chinese), 34(3): 40-44.
Wu L H, Li T, Chen Z. 2012. Reliability analysis of deep borehole strainmeter data in Beijing Baishan. Journal of Geodesy and Geodynamics (in Chinese), 32(3): 41-44.
Wu Y, Ding Z F, Wang X C, et al. 2018. Crustal structure and geodynamics of the North China Craton derived from a receiver function analysis of seismic wave data. Chinese Journal of Geophysics (in Chinese), 61(7): 2705-2718. DOI:10.6038/cjg2018L0244
Xu C, Feng Y Y, Huang L. 2009. Microstructure features of voids in cement-bentonite slurries. Hydrogeology and Engineering Geology (in Chinese), 36(4): 90-94.
步玉环, 柳华杰, 宋文宇. 2011. 晶格膨胀剂对水泥-套管界面胶结性能的影响实验. 石油学报, 32(6): 1067-1071.
陈培善, 白彤霞. 1991. 震源参数之间的定量关系. 地震学报, 13(4): 401-411.
陈运泰, 林邦慧, 林中洋, 等. 1975. 根据地面形变的观测研究1966年邢台地震的震源过程. 地球物理学报, 18(3): 164-182.
池顺良, 武红岭, 骆鸣津. 2007. 钻孔应变观测中潮汐因子离散性与各向异性原因探讨——"十五"数字地震观测网络分量钻孔应变仪首批观测资料分析解释. 地球物理学进展, 22(6): 1746-1753. DOI:10.3969/j.issn.1004-2903.2007.06.010
池顺良, 高睿, 王琳, 等. 2008. 营口地震台伸缩仪和分量钻孔应变仪记录资料比对. 国际地震动态, (3): 34-37. DOI:10.3969/j.issn.0253-4975.2008.03.004
董培育, 任天翔, 杨少华, 等. 2015. 钻孔应变观测同震阶变时的一个问题及对策分析. 地质力学学报, 21(3): 359-370. DOI:10.3969/j.issn.1006-6616.2015.03.006
黄福明, 王廷韫. 1980. 倾斜断层错动产生的应力场. 地震学报, 2(1): 1-20.
李斌, 李自红, 宋美琴, 等. 2019. 山西原平ML4.7地震矩张量反演. 太原理工大学学报, 50(5): 598-605.
李莹甄, 殷娜, 李小晗. 2014. 不同震级标度转换关系研究概述. 地震工程学报, 36(1): 80-87. DOI:10.3969/j.issn.1000-0844.2014.01.0080
龙锋, 闻学泽, 徐锡伟. 2006. 华北地区地震活断层的震级-破裂长度、破裂面积的经验关系. 地震地质, 28(4): 511-535. DOI:10.3969/j.issn.0253-4967.2006.04.001
欧阳祖熙, 张宗润. 1988. 钻孔应变仪与井壁耦合方法的研究. //地壳构造与地应力文集(2). 北京: 地震出版社, 1-10.
潘立宙. 1980. 地应力与地应变的定向测量. 力学与实践, (1): 20-26, 33.
邱泽华, 石耀霖. 2003. 地震造成远距离应力阶变的观测实例. 中国科学(D辑), 33(S1): 60-64.
邱泽华, 池顺良. 2013. YRY-4型钻孔应变仪观测的P波剪应变. 地震, 33(4): 64-70.
邱泽华. 2017. 钻孔应变观测理论和应用. 北京: 地震出版社.
石耀霖, 范桃园. 2000. 地应力观测井中元件标定及应力场计算方法. 地震, 20(2): 101-106.
宋美琴, 郑勇, 葛粲, 等. 2012. 山西地震带中小震精确位置及其显示的山西地震构造特征. 地球物理学报, 55(2): 513-525. DOI:10.6038/j.issn.0001-5733.2012.02.014
唐磊, 邱泽华. 2011. 钻孔四分量应变仪观测的汶川地震的同震应变阶. 地壳构造与地壳应力文集: 114-124.
王长在, 吴建平, 杨婷, 等. 2018. 太原盆地及周边地区双差层析成像. 地球物理学报, 61(3): 963-974. DOI:10.6038/cjg2018L0114
王秀文, 赵虎明, 靳玉科. 2018. 山西五台山断裂形变异常的分析研究. 国际地震动态, (4): 74-79.
吴昊昱, 李丽, 王霞. 2018. 2016年4月7日山西原平M4.2地震震源机制解与发震构造探讨. 科技通报, 34(3): 40-44.
吴立恒, 李涛, 陈征. 2012. 北京百善深井钻孔应变仪观测资料可靠性分析. 大地测量与地球动力学, 32(3): 41-44.
武岩, 丁志峰, 王兴臣, 等. 2018. 华北克拉通地壳结构及动力学机制分析. 地球物理学报, 61(7): 2705-2718. DOI:10.6038/cjg2018L0244
徐超, 冯颖彦, 黄亮. 2009. 水泥-膨润土泥浆固结体的微观孔隙结构特征. 水文地质工程地质, 36(4): 90-94.