地球物理学报  2017, Vol. 60 Issue (11): 4327-4340   PDF    
地震的应变张量观测与应用前景
和泰名, 李世愚     
中国地震局地球物理研究所(地震观测与地球物理成像重点实验室), 北京 100081
摘要:地震发生时的动态应变场,在研究地震触发、地震破裂、地面破坏、水文和岩浆变化等方面都具有重要应用意义.地震的应变张量观测和现有的惯性地震仪观测的物理量不同.前者可以直接记录到地震发生时震源辐射的应变(应力)波,而后者记录到的是位移、速度或加速度.地震频率的应变测量在地震学中的应用前景主要表现在:①测量震源机制解理论预言的辐射4象限分布;②测量库仑应力变化;③换算成动态应力以评估地震烈度;④测量地震波的能量密度;⑤测量地震断层形变加速和形变局部化过程.用惯性地震仪的记录虽然在理论上也可以解算出动态应变值,然而种种原因导致计算结果的误差很大,往往不可接受.应变张量地震仪若能与现有的惯性地震仪配套起来,形成大规模台阵,则有可能推动应变地震学的诞生,在地震观测和地震学科领域引起重大革新.
关键词: 应变张量地震仪      惯性地震仪      地震观测      应力      应变地震学     
The seismological application prospect of strain tensor meters
HE Tai-Ming, LI Shi-Yu     
Key Laboratory of Seismic Observation and Geophysical Imaging, Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
Abstract: The dynamic strains associated with seismic waves may play a significant role in earthquake triggering, hydrological and magmatic changes, earthquake damage, and ground failure. The physical quantity recorded of strain tensor meter is completely differ to that recorded of existing inertia seismometer. The former can directly measure the strains associated with seismic waves, and yet the latter measurement is displacement, velocity and acceleration. The application prospects of dynamic strains meter with seismic wave frequency are mainly (1) measurement of the radiation quadrant distribution of earthquake source mechanism; (2) measurement of Coulomb stress changing; (3) evaluation of earthquake intensity by converting dynamic strains to dynamic stresses; (4) measurement of energy intensity of seismic wave similarly; (5) measurement the process of strain acceleration of seismic fault and deform localization. Although the data recorded of the inertia seismometer can also be calculated to strains, however the error of the results are often too large due to the dependence of theory model and other reasons, and often be unacceptable. If the fiber laser strain tensor meters can be assorted to inertial seismometers and form large-scale seismic array, the strain seismology may be born and cause grave revolution in the area of earthquake observation and seismology.
Key words: Strain tensor seismometer    Inertial seismometer    Seismological observation    Stress    Strain seismology    
1 引言

测量地震发生时的动态应变场,在研究地震触发、地震破裂、地面破坏、水文和岩浆变化等方面都具有重要应用意义.这种测量的完全实现将取决于应变张量地震仪的普遍应用和规模化.应变张量地震仪和现有的惯性地震仪所记录的物理量完全不同.前者可以直接记录到地震发生时震源辐射的应变(应力)波,而后者记录到的是位移、速度或加速度.惯性地震仪和应变仪的记录内容不相同这一事实一直受到地震学家的关注.据Gomberg和Agnew(1996)的记载,对两种地震仪器记录的关系的研究最早可以追溯到Benioff(1935).该文论述了一些基本方程,还设计了一个实用的应变地震仪.然而在当时的技术条件下,要想实现地震的应变张量观测显然是不那么容易的.该文以及后来Benioff和古登堡(1952)合作的论文讨论了不同类型的地震波在线应变与地面位移在方位响应上的不同.他们提出,将两种记录结合起来可以区分不同阶段的地震波的震相.Romney(1964)对此进一步作了发展,他提出应用应变记录仪和惯性地震仪的不同记录来区分噪声和P波.这项工作激发了后来大量的研究工作,Fix(1973)做了综述.Mikumo和Aki(1964)提出新的想法.他们指出,由于应变和速度的比例与相速度有关,因而可以将应变和单台位移的记录联系来估计相速度,而不必应用台阵. Mikumo和Aki采用的数据只限于单方向应变和水平向位移.Sacks等(1976)指出这些技术的应用在甚长周期(>100 s)的范围内遇到的困难是由于水平向惯性地震仪的倾斜造成的.他们主张采用体(膨胀)应变仪和垂直向运动数据,至少对于这些波的类型来说这些数据是非零的.然而,他们并没有处理任何数据.Borcherdt等(1988)处理了一些局部和地方地震的P波和S波的地面运动和应变数据,但这些数据仅局限于体应变的单一分量.Borcherdt(1989)通过理论的研究,认为上述数据的处理结果和理论相符.

2 用惯性地震仪的记录能否精确地换算成动态应变值

Gomberg和Agnew(1996)指出,动态应变值的精确测量是十分重要的.在一些情况下,动态应变场的很小差异就能在过程中导致后果很大的不同.精确评估动态应变的能力在工程方面也具有重要意义.例如,地震辐射的面波在震中距为几十到上百公里的地方就能造成地震破坏.

问题是,用惯性地震仪记录的动态位移能否精确地换算成动态应变值?

Gomberg和Agnew(1996)对此做了尝试,他们比较了美国加州Pinon Fla Observatory(PFO)台的三分量长基线激光应变仪记录和同一台站惯性地震仪的数据,对美国加州的三个地震的记录数据进行了处理.其中PFO激光应变仪的测量基线长度为732 m, 应变的分辨率在地震频段达到10-4 με.系统的应变线性限辐值为±3.2 με, 记录系统的限辐值为±0.2 με,对应的当地加速度响应值为0.1 g.这个值要大于讨论的这几个地震的当地记录.在频率为0.05~1 Hz区间可将信号与地面噪音分离.用干涉条纹方法,应变测量的精准率达到0.1%.记录系统响应的频率域平坦部分为为0~1 MHz, 但由于基线长度的影响,在10 Hz以上部分幅值压低.在记录输出部分采用4极点,极值为2 Hz的低通滤波,这样得到的数据在频率低于0.5 Hz部分的幅值响应是平坦的,相位延迟为-0.3 s.数据记录在10 Hz, 因此不会有明显压低.惯性地震仪采用Streckeisen STS-1系统,由IDA/IRIS组操作.仪器的带宽为0.003~0.5 Hz, 数据记录在20 Hz.

对惯性地震仪的数据的处理采用交互谱(Cross-Spectral)的方法.该方法由Berger等(1979)提出.该方法主要应用于接近平面波的传播, 且频散不大的情况.其原理推导可参见Jaeger和Cook(1979).假定平面波沿x轴方向传播,其位移的一般解为达朗伯公式:

(1)

其中fF为与初始条件有关的函数, 并分别代表沿x的正方向及负方向传播的波;c为波速.在c为接近常数,即频散不大的情况下,质点的速度和位移有如下普遍关系:

(2)

或写成应变

(3)

交互谱方法的具体步骤是:(1)将地震仪的记录数据校正以得到地面运动速度; (2)将水平向分量旋转至径向得到径向线应变εrr, 剪应变ε和环向线应变εφφ,其中φ为台站相对震中的方位角; (3)将数据进行低通滤波(使用巴特沃斯滤波器,为4阶分频斜率,拐角频率为1 Hz)取十分之一为10 Hz(以便和应变数据匹配); (4)将数据做傅里叶变换,并将傅氏谱除以基本模型的相速度,其中径向取瑞利波,剪切向取勒夫波的2倍波速;(5)进行傅里叶反变换来得到应变时间序列的估计值.该方法要求观测点的震中距比较远;振相比较单一;或者其它模式的相速度的没有太大区别,这样仅用单一模型不会产生很大误差.

Gomberg和Agnew(1996)处理的这三个地震的震中距分别为470, 345和206 km.在这些距离上可以将地震波近似视为平面波.他们通过带通滤波处理了低于10 Hz频率的面波数据.分析结果表明,在多数情况下,惯性地震仪的振相到时与应变记录符合得很好,偏差在±0.1~±0.2个周期之间,但是幅值则有系统性偏离,大约超出20%;在有些情况下偏离两倍之多.另外还发现虽然εφφ在理论上用惯性地震仪的计算数值应为零,但应变仪实际记录并非如此.Gomberg和Agnew(1996)认为,二者幅值上的偏差主要来源于相速度存在频散,而εφφ的非零记录与台站附近的构造不均匀性、地形地貌的起伏等有关.

Roeloffs等(2007)在评估报告中对Gomberg和Agnew(1996)的工作结果的评价是:一方面,这个发现意味着高频记录的应变可能会带来地震波的新知识.另一方面,这些内容很难理解, 还可能会对用地震频率应变数据来校正应变仪形成障碍.

3 地震频率应变测量的应用意义

Gomberg和Agnew(1996)指出,地震频率的应变数据之所以重要有两个原因:

(1) 应变与应力联系更密切,而一般认为应力才是决定动态触发地震的相关物理量;

(2) 地面运动的应变值在工程上是对预测地震破坏有用的量.

以下给出一些具体的例子.

3.1 地震断层面解

首先我们来比较同一个纯走滑地震断层辐射的地震波场(库仑应力变化)的零频附加位移(图 1)和应力分量等值线, 这里只显示了最大主应力σ1c和平均主应力(图 2)(李世愚等, 2010, 2015).其中应力降为Δτ=-67.7 MPa;u为与断层走向平行(x轴)的位移值;v为与断层走向垂直(y轴)的位移值;μ为介质的剪切模量.

图 1 一个纯走滑断层的附加位移等值线.应力降Δτ=-67.7 MPa.破裂面两侧的箭头指断层面错动的方向 (a)与断层走向平行的位移值2μ·uc; (b)与断层走向垂直的位移值2μ·vc. Fig. 1 Numerical calculation of the displacement contours for a strike slip fault, where stress drop Δτ=-67.7 MPa and the arrows indicate the sense of dislocation (a) The displacement 2μ·uc parallel to fault strike. (b) The displacement 2μ·vc perpendicular to fault strike.
图 2 一个纯走滑断层的附加应力场(库仑应力变化).相应的应变分量可以用张量应变地震仪测得.应力降Δτ=-67.7 MPa (a)最大主应力σ1c. (b)平均主应力 Fig. 2 Numerical calculation of the σ1 contours for a strike slip fault, where stress drop Δτ=-67.7 MPa Corresponding strain can be directly recorded by strain tensor meters. (a) Maximum main stress σ1c; (b) Average main stress

图 1的内容可以利用惯性地震仪或GPS直接测到;图 2的应变值则首先由张量应变仪测出应变量,然后用广义胡克定律换算得到.

可以直观地看出,图 1图 2的内容完全不同.

图 3a表示一个震源理论模型, 走滑断层FF′两盘水平运动, 震源附近的区域被断层面FF′和与之正交的辅助面AA′分为四个象限.在这些象限里, P波的初动是被压缩(以实心圆圈●表示)或膨胀(以空心圆圈○表示).图 3b为相应的震源机制解的“沙滩球”.在表达震源机制解的参数时,压应力轴(P轴)位于初动是膨胀(-)的象限,张应力轴(T轴)位于初动是压缩(+)的象限,见图 3c(陈运泰,顾浩鼎, 待出版).

图 3 由一个纯走滑断层的运动产生的地震波初动的压缩与膨胀的分布 (a)断层面FF′与辅助面AA′; (b)震源机制示意; (c)地震断层释放的应力的主轴. Fig. 3 Distributions of compression and dilatation first motions for a strike slip fault (a) The fault FF′ and auxiliary AA′ plane; (b) Focal mechanism; (c) The pressure axis and the tension axis.

可以看出,图 3的震源周围压缩区和膨胀区可以应用张量应变仪直接测出体应变θθ=εiic.可以将θ换算成应力球张量σ0=σiic/3=.右上角的符号c表示同震的附加应力或应变.假定观测点位于地面的自由表面上,σ3c=σzz=0, σ0就可以表示为

(4)

其中K为体积模量.p值分布见图 2b.然而迄今为止,所有的震源机制解都不是这样直接测出来的,而是根据惯性地震仪位移(或速度)的垂直向初动方向测量,通过理论模型推算出来的.在断层面运动情况比较复杂的情况下,例如存在膨胀分量、线性补偿偶极子等情况下,惯性地震仪的垂直向记录存在矛盾符号.另外,在井孔台站的深度与震源深度相同或更浅(例如矿山)的情况下,利用垂直向初动方向的判断P波初动是压缩还是膨胀就更加困难.在上述情况下,所有的震源机制解算结果都是依赖于理论模型,为了解决结果的多解性和不确定性, 不得不对理论模型加以某种修补.事实上,地震震源不仅可以用运动源(例如位错源)来表征,也可以用等效的双力偶源来表征(Maruyama, 1963; Burridge and Knopoff, 1964; 陈运泰和顾浩鼎,2017).后面我们会说到,震源辐射的地震波,不仅可以表征为运动方程,更可以表征为应力波方程.因此,单纯用运动源来表述震源机制,用运动方程来表述地震波辐射,在理论上至少是不完整的,在操作上也是不简洁的.震源理论的现状并不是出于必然,而是出于一种无奈,即由于缺少应变张量的直接观测数据造成的“巧妇难为无米之炊”.

3.2 库仑应力变化

地震破裂的库仑应力变化对于判断震后其它地区地震发生的危险性是否加重具有重要意义.然而迄今为止,库仑应力变化的计算仍然是利用惯性地震仪的位移测量数据,通过理论模型推算出来的, 而不是应用张量应变仪直接测得的.这样一来,推算结果的可靠性和精确性就要依赖于理论模型的选择.由于震源机制解的复杂性,反演结果的不唯一性, 加上场地的地质构造和介质的不均匀性、各向异性,可能造成库仑应力变化计算结果的误差,在缺乏应变张量测量的情况下,计算结果的可靠性难以验证.例如,目前大多数计算是应用伏尔泰拉模型计算得到的.该模型假定地震断层两盘的位错量是个常数(即在断层面上处处相等).李世愚等(2010)和泰名等(2011)李世愚等(2015)通过理论计算证明,用伏尔泰拉模型的结果会造成部分断层面的应力降出现奇异性, 和破裂端部前缘应力奇异性构成双重奇异性, 这是一个物理上的重大缺陷.他们主张用滑动弱化模式,或近似用线弹性断裂力学模式及其相应的索米亚那位错模式来取代伏尔泰拉位错模式.但效果如何,有待于应变张量观测数据来验证.

在未来,如果能用合理数量和密度分布的应变张量地震仪测到同震的附加应变场εijc,就可以利用广义胡克定律得到库仑应力变化,假定介质可以看成均匀、连续、各向同性的,就可以计算为

(5)

其中μ为剪切模量; λ为拉梅常数; δij为克罗内克符号,.

另外,在得到动态的库仑应力变化σijc后,还需要在事先知道当地的背景应力(或称为本底应力) σijb, 再将二者叠加,得到

(6)

才能用σij值(绝对应力)来进行破坏的可能性估计.由此可见,仅仅测量库仑应力变化是不够的,还需要测量当地的背景应力,但这同样需要足够密度分布的应力张量测量台网.

3.3 地震时地面烈度估计

世界上一些早期的烈度表曾经把等价的地面加速度作为烈度表的一部分.但近期的研究表明,烈度与地面运动的加速度等参数不存在简单的换算关系,由地面加速度换算成烈度还没有找到科学依据.用加速度谱进行烈度估计有很大的离散性,上下限的数值可以相差几倍到几十倍.有证据表明, 地面加速度不是影响烈度的非常重要的参数.因此,新的《欧洲地震烈度表》中就不再把地面运动参数(如加速度)包括到烈度表中来(鄢家全和郝玉芹,2005).

从物理本质上,地震时介质的破坏是由动态应力决定的.造成破坏的根本原因并不是加速度本身,而是介质加速度的突发性、不均匀性和不协调性造成的局部附加应力.还要加上背景应力当地的绝对应力, 见式(6).这种不协调性,主要由场地局部的地质构造、地形地貌、建筑物的结构、几何形状造成的对地震波的响应所决定.例如,火箭在加速时并不能产生箭体的破坏,就是因为火箭的严密结构使得各部位在加速时处于协调状态.由于我们现在只有运动参数(位移、速度或加速度)的记录,因此用运动参数来预测震害烈度实际上是一种不得已的办法.应变张量地震仪未来在敏感部位的普遍铺设,有希望测出该部位在接近破裂极限条件下的应力状态,为烈度估计提供最需要的动力参数.在进行烈度估计时,我们需要在几种判据中做出选择:

(1) 对于场地及建筑结构,可以利用计算出叠加后σij的最大主张应力σ1,以及σ1的峰值σ1top的持续时间t(σ1top).以及在地下出现σ1>0的情况.在进行烈度估计时还需要测量σij峰值的持续时间t(σijtop).

(2) 对于相对完整的岩体,可以利用库仑-莫尔准则;

(3) 对于已有的断层或薄弱带,可以利用拜尔利定律.

众所周知, 这些准则的内容主要都是以应力观测值来表述的.因此应变张量测量无疑是必要的.

在一些特殊情况下,应变张量的测量对于烈度评估具有不可替代的优势, 例如在一些慢振动中,地面的运动位移、速度或加速度都很小,但某些场地的应力-应变却可能已经超出介质的承受极限而出现破坏,这时的惯性地震仪并没有明显记录, 我们可以用拔河比赛双方僵持时的绳子状态来比喻.在1975年海城地震和1976年唐山地震后在一些地区报告的“有感无震”现象可能就当属此类.

3.4 地震波能量

根据克拉贝龙定理,地震波携带的应变能密度为

(7)

对于均匀、各向同性介质,应变能密度增量还可以写成

(8)

可以看出,其中所有的应变分量的增量都可以由张量应变仪直接测出.单台记录的地震波总能量由整个地震事件的测到的ES

(9)

其中

式中t0为首波到时,tm为最后震相结束时刻.c=ckek是波传播的速度; u是观测点的位移,dS曲面S上的面积元.Ek称为动能密度(即单位体积中的动能),Ep称为位能密度.可以证明(Jaeger and Cook, 1979; 郭自强, 1982), Ek=Ep.即波在传播过程中动能密度与位能密度总是等分的.而地震波的总能量

(10)

从式(9)可以看出,Ek可以用惯性地震仪测到,而Ep可以由应变地震仪直接测出.因此,应变地震仪和惯性地震仪在测定地震波能量方面可以同样起作用.在这个问题上,还不清楚哪个方法更优越.但两者结合起来,至少可以起到互相校验的作用, 使得观测体系更加可靠、完整.另外,由于震源辐射波的各向异性,Es在曲面S上的分布为倾角(δ)和方位角(φ)的函数,因此应变地震仪台网中台站需要有足够的数量和合理的分布.

3.5 地震前物理场

大地震前物理场的测量关系到地震前兆和地震预测.根据岩石断裂力学, 地震断层破裂的起始判据,主要是断层面上的应力是否超出了断层介质的最大承受能力.而其它物理表现都是派生现象.因此,断层及其附近介质的应力-应变状态应该是震前物理场测量的核心内容.近年来,顾国华等(2009, 2015)等利用GPS资料,池顺良等(2009), 邱泽华等(2010)利用应变测量都在这方面进行过一些有益的尝试,并提出过类似的看法.不过由于国内相关数据的稀少,还需要更加有说服力的证据.

由于缺少动态应变张量观测,其它派生物理表现的观测就失去了核心性参照.例如地震的活动性.现有的数据表明,这些派生表现(如大地电磁、水化学)和震前中小地震的活动不存在时间上的同步性.另外, 从中小地震活动的各种图像和参数中,对于如何识别前震的特征标志至今没有确切的判据.1975年海城MS7.3地震有直接前震,1976年唐山MS7.8地震以及其它大多数大地震前都没有公认的前震,表现为震前平静.最近的看法就是前震不是大地震前的必然表现.即使是少数有前震的大地震,识别小震群是否属于前震始终没有决定性的判据.不要说在震前,就是在震后也拿不出来识别的依据.例如2008年汶川MS8.0之前, 发生在汶川附近都江堰一带的3.9级小震群是否前震, 有没有区别于一般震群的特征?至今没有拿出令人信服的论证,就是一个例子.

地震前的形变加速或局部化现象一定伴随小地震的密集,这个误判来自对一些小样品实验的误读.图 4测量的是完整岩样的表面应变场,显示了在临近破坏前的形变局部化(许向红, 2005).图 5显示的是在小样品岩石破裂过程中不同阶段获得的声发射(AE)定位, 在临近破坏前出现AE的高度集结(成核) (Lockner and Madden, 1991).这两个典型实验曾成为地震孕育的成核模式的重要依据, 但对它们的不适当的宣传也到了误导作用.应当指出,这两个实验的观测系统和观测内容是完全不同的.前者为应变测量系统,和大尺度应变测量相似,后者为声发射测量系统,和惯性地震仪的观测相似.对于公里尺度情况来说,两种效应之间有时候是等效的,但只是有时候,实际上不一定都是如此.长期以来,我们总是把这两个实验结果等同起来,并且无条件地外推到一切类型的地震,从而得出大地震必然有前震的推论.这种惯性思维往往容易得出一些错误推断.1976年唐山MS7.8地震颠覆了这个传统理念, 使我们反思这个理念在公里尺度条件下被颠覆的原因.

图 4 用白光数字散斑方法测量的岩石表面应变场ε22* (许向红, 2005) 材料:大理岩.尺寸: 20×20×40 mm3. Fig. 4 Rock surface strain ε22* field measured by Digital Speckle Correlation Method(Xu Xianghong, 2005) Material: marble.Size: 20×20×40 mm3.
图 5 (A) Westerly花岗岩破裂过程中不同阶段获得的AE定位.上图从断层走向方向显示AE事件,下图从正对断层面显示AE事件; (B)相应位移-应力曲线, 变形过程中(a)—(f)阶段中记录到的声发射源的位置见上图(Lockner and Madden, 1991; 陈颙等,2009) Fig. 5 (A) Acoustic emissions (AE) hypocenter locations during fault formation of initially intact Westerly granite.Time processes from left to right.Upper plots show events viewed along-strike of eventual fault plane. Lower show same events when fault plane is viewed on face-on.(B) Accompanying displacement-stress curve indicates segments of the experiment from which acoustic emission plots are made (Lockner and Madden, 1991; Chen Yong et al., 2009)

(1) 上述实验是用脆性样品进行的,而根据现有大地震的震源深度,大地震的发生前,孕震区的深度往往在脆-塑转换带.因此对这两个实验的等效性外推到公里尺度不满足相似性第三定理中的物理状态相似条件.

(2) 实验岩石破裂实验证明,岩石内部约99%的破裂都不伴随声发射,尤其是塑性破裂和蠕变.真正伴随声发射的只有不到1%的脆性破裂(Lockner and Madden, 1991; 蒋海昆和张流, 1998; 李世愚和和雪松,2003).

(3) 完整岩样的破坏实验和已有断层的摩擦过程不完全等同,例如声发射(AE)的b值,完整岩样在破坏前b值以下降为主,而摩擦过程(耿乃光, 1986)在滑动之前b值以升高为主.由于凹凸体闭锁等原因,摩擦过程的粘滑阶段完全有可能缺少中小地震而表现为震前平静.实际的地震过程既有岩石破裂,又有断层面的摩擦,两者的前兆有相同的表现,也有相互为负面的效应.两种机制的交叉造成了前兆的复杂性和多样性.

以上论述表明,地震活动性图像并不能全部反映地下的形变情况,甚至不能全部反映地下断层的破裂情况.由此证明,离开应变应力测量,单纯地把捕捉地震前兆的希望寄托在直接前震上的思路早已走入死胡同.大地测量的反演结果表明,唐山MS7.8地震前存在大规模地面形变加速现象(陈运泰等, 1979尹祥础和郑天愉,1982), 其后的一些震例表明这个现象有可能是大地震前的普遍现象, 而这个过程不一定伴随中小地震的震群.

遗憾的是,在缺乏前震的情况下,现有的惯性地震仪对于慢事件反映不敏感,即使记录到了也很难被认可.惯性地震仪组成的台网对于识别这种形变加速的现象无能为力;而GPS虽然对于慢速位移的记录比较擅长,但由于GPS台网比较稀疏,有可能捕捉不到形变局部化形成的破裂成核带.这样,即使前震出现,也很难做出判断(Brodsky and Lay, 2014).由此可以看出, 研究地震前物理场变化的战略重点应该将地震活动性与地下介质的动态应变或应力状态图像变化的跟踪观测结合起来.

在断层破坏的滑动-弱化模式中,我们不仅需要测量断层的应力降,也要测量断层的分段滑动情况,但这需要将应力(应变)测量和滑动位移(或滑动速率)测量巧妙地结合起来.因此,如何在活断层上布置测量台阵,是一个需要深入研究和反复试验的重大课题.我们不妨利用滑坡场地作为中尺度试验场,或借鉴现有滑坡问题的研究成果.

4 地震的应变张量测量的若干问题 4.1 地震应变张量测量的理论基础

我们从最基本的应力形式的波动方程开始.根据牛顿第二定律,有

(11)

其中fi为体力,ρ为介质密度.按照小变形的假设,可将应力按泰勒展开,并略去二阶及二阶以上的高阶小量.我们还假定物体未变形时内部没有应力,得到应力应变的关系:

(12)

在介质为均匀、各项同性和线弹性的假设前提下,应力应变的关系才能写成

(13)

其中为体积应变,δij为Kroneker符号.

在小变形的假设前提下,应变与位移的关系才能写成

(14)

将广义胡克定律式(13)和位移-应变关系式(14)代入波动方程(11), 并写成相应的位移向量形式为

(15)

其中ρ为介质密度,f为体力,为体积应变.利用恒等式Δ2u=Δ(Δ·u)-Δ×Δ×u,式(15)又可以写成

(16)

假定体力为保守力,或是常力, 即fΩ, 其中Ω为势函数, 满足Δ2Ω=0.对式(16)两边求散度, 利用恒等式Δ(Δ×u)=div rot u=0,得

(17)

上式就是用应变表示的涨缩波、纵波或P波.

对式(16)两边求旋度, 并利用恒等式Δ×(ΔΩ)=0, 得

(18)

其中2ω表示为2ω=rot(u)=Δ×u.式(18)叫做旋转波、横波或S波.为了导出式(17)和(18)的应变表示,我们考虑平面波的传播.不失一般性,假定平面波沿x1方向传播,那么位移分量u1, u2, u3就只是x1t的函数, 因而式(17)和(18)可以写成

(19)

(20a)

(20b)

由于位移分量u1, u2, u3只是x1t的函数, , 因此, 波动方程(19)、(20a)和(20b)就可以写成应变形式

(21)

(22)

其中分别为纵波和横波的波速.

由于实际的应变测量只能测到线应变,因此我们还需要给出式(22)的线应变换算式.这里举一种测量方案为例.

假定直角坐标Ox1x2x3的三个坐标轴,不妨假定x1轴为地理坐标的东西向,x2轴为北南向,x3轴为垂直向向下.OP, OQ, OR分别位于Ox1x2, Ox2x3Ox3x1平面内,与Ox1, Ox2Ox3轴的夹角分别为δ, φ, ψ, 见图 6.分别测量ε11, ε22, ε33, 以及εP, εQεR.由坐标变换,

图 6 应变张量测量 Fig. 6 Strain tensor measuring method

(23)

因而

(24)

同理,

(25)

(26)

特别地,当δ=φ=ψ=45°,

(27)

对于沿x1轴传播的平面波,ε22=ε33=εQ=0,

(28)

由于应变张量为对称张量,ε12=ε21, ε23=ε32, ε31=ε13, 因此其中只有6个独立分量.这样,Ox1x2x3坐标系中εij的9个应变分量就全部得到.

进一步计算矩阵εij的特征值ε和特征向量,即可算出应变张量的三个主应变εi和主方向Ni(i=1, 2, 3).

在实际布设中, 往往还加上一些“额外”分量,例如在Ox1x2平面中增加εS,其中OSOP, 见图 6.这种结构主要用于“自洽”分析,其原理是ε11+ε22+ε33=εP+εS+ε33=II为应变张量第一不变量,因而只要测量系统附近介质保持连续,就恒有ε11+ε22=εP+εS.当测量系统基线附近介质出现损伤时,这种“自洽”关系便被改变而出现异常.

从以上讨论可以看出:

(1) 三维的应变张量有6个独立分量.仅设置三个水平方向的测量是不够的,还需要增设3个其它非水平向,包括一个垂直向和两个倾斜向.

(2) 对于地表布设测量,OQ, OR两个分量可以省略,但考虑到地表往往不是理想水平的,或地表形状不太平整,或台基往往布设在山洞里,并不是真的在地表,因此在地面布设中仍然需要保留OQ, OR两个分量.

(3) 应变张量测量系统的“自洽”性质出现异常与测量系统附近的活动构造破裂“成核”的过程是否关联?这个问题需要慎重考察.从原理上看,只有与应变张量的测量系统耦合的介质在基线范围内遭到损伤时,才会导致测量系统的“自洽”性质出现异常.实际的测量系统附近的构造很多,无论大小、远近,微破裂的演化总是存在的,需要距离多远,损伤程度多小才不关联?目前缺少准确的说法.

4.2 应变张量地震观测设备的研制和使用现状

如前所述,用现有的惯性地震仪得到动态应变结果不够精确,人们自然将希望寄托在应变张量的观测上.地形变测量主要包括GPS和应变测量(体应变、应变张量等)、InSAR技术以及地倾斜等.据了解,顾国华等(2009, 2015)用GPS资料对于汶川8.0地震、东日本MW9.0大地震前、同震及震后地壳水平运动进行了有益的分析.GPS观测在近年来也已经采用地震频率采样.这里主要针对应变张量的观测设备的研制与使用.

邱泽华和张宝红(2002)介绍了我国钻孔应变(应力)观测现状; 邱泽华和石耀霖(2004)对国外钻孔应变观测的发展现状进行了综述.在国际上,1968年9月,美国卡奈基研究所的Sacks和德克萨斯大学的Evertson研制出世界上第一台钻孔(体)应变仪(Sacks and Evertson, 1971).中国首先研制出RZB-1型分量应变仪.之后Gladwin研制出张量应变仪.30多年来,钻孔应变仪先后在中国、日本和澳大利亚被研制出来,并得到迅速发展和日益完善.欧阳祖熙(2011)着重介绍了美国PBO(Plate Boundary Observatory)计划建立的钻孔应变仪台网.

Roeloffs等(2007)对美国地质调查局1985—2004年所收集的钻孔应变仪数据进行了评估.该报告的主要篇幅是针对与潮汐和气象变化的周期相近的每小时一次的记录.至于这些高频记录数据在地震学中的意义,在Roeloffs等(2007)的报告中没有作为重点.这里需要说明GTSM制造的Gladwin应变仪在配置上存在的几个问题:

(1) 加利福尼亚的GTSM没有高频数据记录.据Roeloffs等(2007)报道,Gladwin本人认为这些仪器不是为记录高频数据研制的.Roeloffs等(2007)则举出不少例子来说明,在“高频”范围内,所有其它的应变仪都没有出现任何问题,即使记录长期变化不太准确的Motocross和Chabot体应变仪也有着非常好的远距离地震记录。因此,安装的好坏对钻孔应变仪在这个频段内获取数据影响不大.因此Roeloffs等(2007)建议GTSM增加高频记录.

(2) 没有垂直向应变记录.本文作者在地震局地球所2014年交流的过程中曾经向Gladwin本人咨询过这个问题.首先是有没有因为技术上的困难.Gladwin的回答是不存在技术上的障碍.

(3) 与上一个问题同时存在的就是GTSM应变仪在应变换算中实际上默认观测点的应变是在一个自由表面的问题,即有两个主应变都在水平向,而垂直向的应变εzz为主应变,这只适用理想水平地面的情况,在100~200 m深的钻孔内一般情况下不普遍成立.这样,GTSM应变仪的应变换算在理论上就存在严重失误.

近年来,光纤激光技术的引入,由于光纤激光地震仪灵敏度高, 频带宽, 动态范围大, 使得地震仪的性能有了飞速发展.2011年9月,中国地震局在普洱完成数字化中深井综合观测系统(刘文义等,2012).不过,其中的井下光纤光栅综合地震仪只是z分量的加速度检波计, 而不是应变计.

据了解,应变测量除地壳形变外,在工程抗震、滑坡监测等都已经普遍应用,体应变仪在部分观测台站中已经得到应用,但多设置为分钟或更低频率记录.

从以上情况可以看出,应变张量地震仪在地震观测中的发展遇到的障碍主要还不是在技术上,而是在于其它几种因素:

(1) 行业之间沟通不足.

应变张量地震仪的发展应该是参与研制仪器的科研人员、力学专业、工程专业、地震学专业研究人员和管理部门相互沟通和共同协作的结果.目前最大的问题是地震学科研人员对此关注不够,参与不足.这种情况的不良后果就是观测系统的研发缺少需求导向和理论导向,造成仪器的结构和性能与地震应变观测的要求不符.

(2) 管理体制和数据交流共享上的问题.

迄今为止,张量应变仪在只是作为钻孔应变仪系列的一个组分,作为GPS和InSAR的补充来观测地壳变形和地应力场,并不是现有的地震台网的配套部分.在没有形成规模化的台网的情况下,少量的、孤立的数据起不了太大作用.目前最大的问题是应变测量没有实现与地震台网的数据交换和数据共享,更没有实现和惯性地震仪的共台基、共时间轴的同步观测,应变数据对外是封闭的.有鉴于此,必须尽早改革管理体制,并制订未来应变张量测量数据的交换与共享制度.

(3) 应变仪真的需要普遍布设在井下吗?

已有的应变仪几乎都放置在钻孔内,深度达200 m之多,而惯性地震仪多布设在地面.两者的施工成本几乎相差一个数量级.据有关文献(Roeloffs et al., 2007)给出的理由是钻孔内的噪声比较小.但这是针对潮汐频率说的,对于地震频率真的如此吗?我们举一个矿震记录波形为例.图 7显示的是2007年1月1日兖州鲍店煤矿一次矿震垂直向记录.其中蔡家厂(LD1T)和廿里铺(LD2T)是地面流动台,其余都是井下地震计.可以看出, 北风井(BFJS)和南风井(NFJS)的井下记录波形背景噪声较大, 初动清晰度都很差,其它井下记录(GGQS工广台, XLKS兴隆矿)也不太理想,而蔡家厂(LD1T)和廿里铺(LD2T)的地面流动台虽然离震中的距离相对比较远,但记录波形质量明显好于其它台.这些井下地震计在安装最初半年的状况曾经是比较好的,但半年后逐渐转差,和矿区地下环境不稳定,耦合变差有关.

图 7 2007年1月1日记录到的兖州鲍店煤矿一次矿震记录垂直向波形,ML=0,其中LD1T-蔡家厂,LD2T-廿里铺 Fig. 7 Vertical waveforms of mining induced seismic event of January 1, 2007 with ML=0 at Yanzhou Baodian coal mine, where LD1T represent Caijiachang and LD2T represent Ershilipu

可以看出, 台网质量并不完全取决于是否设置井孔地震仪.井孔地震仪的耦合条件一旦变差,维护成本高.地面布设只要适当选点,信噪比未必差,另外维护成本低.

5 地震应变张量观测的创新发展空间

近年来,由于管理部门对地震预测研究的重视,在中国地震局地震观测技术研究院(E-研究院)的主导下,应变张量地震仪的研发已获得进展.应变张量地震观测技术创新正面临重要发展契机.

前面说过, 用应变张量观测地震目前还存在一些问题,而这些问题的改进和解决恰恰是它的发展空间.

(1) 增设非水平向记录

对于山区、地表倾斜或有局部构造地区,非水平分量的观测就更不可忽略.

(2) 设计观测点的位置

应变地震仪记录的多半是台站的局部变化(欧阳祖熙, 2011),因此,必须事前明确应变仪布设的位置与观测目的的关系.如果是为了知道大区域的地壳变形,就必须将应变仪布设在远离断层和局部构造的基岩之内,而如果是为了了解局部构造、断层端部附近、大型建筑的敏感点或地形地貌的敏感点对地震波的响应,就应该反其道而行之,把观测点布设在这些敏感点(其实噪声和气象本身也是需要考察的影响因素).

(3) 以地面与惯性地震仪共点布设为主

地面选点原则和惯性地震仪应在同一观测点,以便相互比对和联合分析.与惯性地震仪共点布设有利于形成应变张量地震仪台网.少量的、孤立的数据不能起到供地震学研究的作用.只有形成规模化张量应变地震仪台网,才能用来反演震源机制,预测强地面运动的烈度,评估地震能量和震级等参数.

(4) 推动应变地震学理论的创新与发展

根据地震学原理,柱面波可由平面波叠加而成,而球面驻波也可以看成空间中有所有方向传播的平面波叠加而成(傅淑芳和朱仁益,1997).既然平面波可以用应变张量来表示,那么柱面波和球面驻波在原则上也可以由应变张量来表示.但从Gomberg和Agnew(1996)的论文可以看出,对于近场问题,这个问题在理论上有关的可行性研究还没有完全解决.

从字面上,根据位移表示定理, 位错点源的辐射波场既然可以用位移来表示,就可以直接表示成应变张量场,但蔡永恩教授认为,这个问题看似容易,其实在理论上并没有很好解决.

6 讨论

有关应变地震学的理论还会遇到许多新的问题.随着应变张量地震仪的普遍使用和观测数据的获得,这些问题的陆续出现解决势必倒逼地震学的理论工作者作为迫切问题进行深入的研究.

可以预言,应变张量的直接观测,将能使震源机制反演、地震参数的计算、地面运动的评估可以得到直接验证,避免间接推算带来的误差,而且操作上都可以得到简化.因此,实现应变张量观测的规模化,并与惯性地震仪共点同步观测,势必引起地震观测的革命,推动应变地震学的诞生, 从而带来地震学的根本性的创新.

在开展有关仪器研发和解决场地观测技术的问题的同时,需要开展更加深入的理论方面的预研究.这些预研究,对于明确应变张量的研究需求,确定观测技术的研发目标无疑能起到理论的引导作用.

致谢

本文在完成过程中与陈运泰院士、陈慧忠研究员、许忠淮研究员、李小军研究员、蔡永恩教授、雷军教授、李丽研究员、刘文义研究员等进行过有益的讨论,李丽研究员提供了重要的参考资料,在此一并表示感谢.

参考文献
Berger J, Agnew D C, Parker R L, et al. 1979. Seismic system calibration: 2. Cross-spectral calibration using random binary signals. Bull. Seismol. Soc. Am., 69(1): 271-288.
Brodsky E E, Lay T. 2014. Recognizing foreshocks from the 1 April 2014 Chile earthquake. Science, 344(6185): 700-702. DOI:10.1126/science.1255202
Burridge R, Knopoff L. 1964. Body force equivalents for seismic dislocations. Bull. Seismol. Soc. Amer., 54(6A): 1875-1888.
Chen Y, Huang T F, Liu E R. 2009. Rock Physics (in Chinese). Hefei: University of Science and Technology Press.
Chen Y T, Huang L R, Lin B H, et al. 1979. A dislocation model of the Tangshan earthquake of 1976 from the inversion of geodetic data. Acta Geophysica Sinica (in Chinese), 22(3): 201-217.
Chen Y T, Gu H D. 2017. Source Theories (In Press).
Chi S L, Chi Y, Deng T, et al. 2009. The necessity of building national strain-observation network from the strain abnormality before Wenchuan earthquake. Recent Developments in World Seismology (in Chinese)(1): 1-13.
Fu S F, Zhu R Y. 1997. Higher Seismology (in Chinese). Beijing: Seimological Press: 50-51.
Geng N G. 1986. The development of b-value simulated experiments and the begining of b-value simulated experiments in China. Acta Seismologica Sinica (in Chinese), 8(3): 330-333.
Gladwin M T. 1984. High precision multi-component borehole deformation monitoring. Rev. Sci. Instru., 55(12): 2011-2016. DOI:10.1063/1.1137704
Gomberg J S, Agnew D C. 1996. The accuracy of seismic estimates of dynamic strains: An evaluation using strainmeter and seismometer data from Piñon Flat Observatory, California. Bull. Seismol. Soc. Am., 86(1A): 212-220.
Gu G H, Wang W X, Xu Y R, et al. 2009. Horizontal crustal movement before the 2008 Wenchuan MS8.0 earthquake obtained from GPS observations in regional network. Acta Seismologica Sinica (in Chinese), 31(6): 597-605.
Gu G H, Wang W X, Zan W, et al. 2015. Preseismic, coseismic and postseismic horizontal crustal movements of the Mw9.0 Tohoku earthquake in Japan, 2011. Geomatics and Information Science of Wuhan University (in Chinese), 40(12): 1669-1676.
Guo Z Q. 1982. Wave in Solids (in Chinese). Beijing: Seimological Press: 100.
He T M, Li S Y, Zhang H K, et al. 2011. Coulomb failure stress change in slip-weakening model and remote triggering of earthquakes. Acta Seismologica Sinica (in Chinese), 32(2): 165-186.
Jaeger J C, Cook N G W. 1979. Fundamentals of Rock Mechanics. 3rd ed. London: Chapman and Hall.
Jiang H K, Zhang L. 1998. Recent studies of acoustic emission for micro-fracture of source in rock. Translated World Seismology (in Chinese)(3): 1-7.
Kuksenko V, Tomilin N, Damaskinskaya E, et al. 1996. A two stage model of fracture of rocks. Pure and Applied Geophysics, 146(2): 253-263. DOI:10.1007/BF00876492
Li S Y, He X S. 2003. Accurately appraising the information amount of seismicity patterns. Recent Developments in World Seismology (in Chinese)(4): 29-31.
Li S Y, He T M, Yin X C. 2010. Introduction of Rock Fracture Mechanics (in Chinese). Hefei: University of Science and Technology Press.
Li S Y, He T M, Yin X C. 2015. Rock Fracture Mechanics (in Chinese). Beijing: Science Press.
Liu W Y, Zhang W T, Li L, et al. 2012. Optical fiber sensory technology: future direction for earthquake precursor monitoring. Earthquake (in Chinese), 32(4): 92-102.
Lockner D A, Madden T R. 1991. A multiple crack model of brittle fracture: 2. Time-dependent simulations. J. Geophys. Res., 96(B12): 19643-19654.
Maruyama R. 1963. On the force equivalents of dynamical elastic dislocations with reference to the earthquake mechanism. Bull Earthq. Res. Inst., Tokyo Univ., 41: 467-486.
Ouyang Z X. 2011. U.S. PBO Project: Borehole strainmeter networks face a challenge. Recent Developments in World Seismology (in Chinese)(10): 19-28.
Qiu Z H, Zhang B H. 2002. Current situation of the seismological networks for borehole stress-strain precursor surveilance in China. Recent Developments in World Seismology (in Chinese)(6): 5-9.
Qiu Z H, Shi Y L. 2004. Developments of borehole strain observation outside China. Acta Seismologica Sinica (in Chinese), 26(S1): 162-168.
Qiu Z H, Zhang B H, Tang L, et al. 2011. Abnormal strain changes observed at Guza before the Wenchuan earthquake. Sci. China Earth Sci., 54(2): 233-240. DOI:10.1007/s11430-010-4057-1
Xu X H. 2005. Effects of mesoscopic heterogeneity on macroscopic deformation, damage and catastrophic failure of rock [Ph.D.](in Chinese). Beijing: Graduate School of Chinese Academy of Sciences.
Yan J Q, Hao Y Q. 2005. Intensity and acceleration applied in engineering fortification against earthquake. Cities and Disaster Reduction (in Chinese)(6): 17-18.
Yin X C, Zheng T Y. 1983. A rheological model for the process of preparation of an earthquake. Science China Chemistry, 26(3): 285-296.
Yin X C. 1985. Solid Mechanics (in Chinese). Beijing: Earthquake Press.
Yin X C. 2011. Solid Mechanics (in Chinese). Beijing: Earthquake Press.
陈颙, 黄庭芳, 刘恩儒. 2009. 岩石物理学. 合肥: 中国科学技术大学出版社.
陈运泰, 黄立人, 林邦慧, 等. 1979. 用大地测量资料反演的1976年唐山地震的位错模式. 地球物理学报, 22(3): 201–217.
陈运泰, 顾浩鼎. 2017. 震源理论(待发表).
池顺良, 池毅, 邓涛, 等. 2009. 从5.12汶川地震前后分量应变仪观测到的应变异常看建设密集应变观测网络的必要性. 国际地震动态(1): 1–13.
傅淑芳, 朱仁益. 1997. 高等地震学. 北京: 地震出版社: 50-51.
耿乃光. 1986. b值模拟实验的进展和我国b值模拟实验的开端. 地震学报, 8(3): 330–333.
顾国华, 王武星, 徐岳仁, 等. 2009. 区域网GPS观测得到的2008年汶川MS8.0地震前的地壳水平运动. 地震学报, 31(6): 597–605.
顾国华, 王武星, 占伟, 等. 2015. 东日本Mw9. 0大地震前、同震及震后地壳水平运动. 武汉大学学报(信息科学版), 40(12): 1669–1676.
郭自强. 1982. 固体中的波. 北京: 地震出版社: 100.
和泰名, 李世愚, 张洪魁, 等. 2011. 滑动弱化模型下的库仑应力变化与远程触发问题. 地震学报, 33(2): 165–186.
蒋海昆, 张流. 1998. 岩石微破裂震源过程的声发射研究进展. 世界地震译丛(3): 1–7.
李世愚, 和雪松. 2003. 应正确评估地震活动图像的信息量. 国际地震动态(4): 29–31.
李世愚, 和泰名, 尹祥础. 2010. 岩石断裂力学导论. 合肥: 中国科学技术大学出版社: 1-468.
李世愚, 和泰名, 尹祥础. 2015. 岩石断裂力学. 北京: 科学出版社.
刘文义, 张文涛, 李丽, 等. 2012. 光纤传感技术—未来地震监测的发展方向. 地震, 32(4): 92–102.
欧阳祖熙. 2011. 美国PBO计划:钻孔应变仪台网遭遇挑战. 国际地震动态(10): 19–28. DOI:10.3969/j.issn.0235-4975.2011.10.005
邱泽华, 张宝红. 2002. 我国钻孔应力-应变地震前兆监测台网的现状. 国际地震动态(6): 5–9.
邱泽华, 石耀霖. 2004. 国外钻孔应变观测的发展现状. 地震学报, 26(S1): 162–168.
邱泽华, 张宝红, 唐磊, 等. 2010. 汶川地震前姑咱台观测的异常应变变化. 中国科学:地球科学, 40(8): 1031–1039.
Roeloffs E, Hodgkinson K, Bryan C, et al. 2007. 对美国地质调查局1985—2004年所收集的钻孔应变仪数据的评估. 地壳构造与地壳应力(2): 1-66.
许向红. 2005. 岩石的细观非均匀性对其宏观变形、损伤和灾变的影响[博士论文]. 北京: 中国科学院研究生院.
鄢家全, 郝玉芹. 2005. 工程抗震设防中的烈度与加速度. 城市与减灾(6): 17–18.
尹祥础, 郑天愉. 1982. 地震孕育过程的流变模式. 中国科学, B辑(10): 922–930.
尹祥础. 1985. 固体力学. 北京: 地震出版社.
尹祥础. 2011. 固体力学. 北京: 地震出版社.