利用断层应力积累探讨大地震发生的地点——以日本2011年东北9级大地震为例

谢周敏，蔡永恩，吉岡祥一，阿部大毅

1 国家自然灾害防治院，北京 100085 2 北京大学地球与空间科学学院理论与应用地球物理研究所，北京 100871 3 日本神戸大学都市安全研究中心，神戸 657-8501 4 日本神户大学研究生院行星科学系，神户 657-8501

0 引言

地震的孕育是震源区一个复杂的物理、化学过程，它的发生是地下岩石突然破裂的结果.即使科学技术高度发达的今天，预报地震的地点、时间和大小(震级)仍然是一个公认的难题.其问题主要在于地球内部的“不可入性”、大地震的“非频发性”和地球物理过程的复杂性(陈运泰，2009).

王仁(1997)认为：岩层断裂的矛盾双方,是促使岩层破裂的推动力同岩层抵抗破裂的阻力的矛盾.前者是地壳在各种外力作用下的应力分布,后者是岩层的强度.地震的孕育和发生就是这对矛盾对立统一的发展过程.因此，从力学的角度看，解决地震预报问题，需要研究震源岩石的强度和应力随时间的变化过程.前者，一般通过岩石力学实验可以得到不同深度、不同温压条件下的一些认识；后者由于受技术条件的限制，除了地表浅层外，很难进行震源深度处的直接测量，只能通过应力变化引起的各种前兆现象，例如地震活动、地表变形、地震波速异常、大地电磁场、地下水位和温度场等，反演震源的运动学特征和动力学环境，进行间接预测.

地震发生地点的预测是地震预报首先要解决的问题.一个地震活跃的区域，其历史地震活动图像反映了这个区域应变能积累、释放和转移的过程.根据板块学说，在地质时间尺度内，板块运动的速率是稳定和连续的.板块间的边界是最明显的地震活动带，在适当长的时间间隔内，地震带内的大地震的破裂区往往会连接、无重叠地展布于地震带.如果地震带内某个地区长期没有发生大地震，那么，这个地区就被视为未来最有可能发生大地震的地区，通常把这个地区称为地震空区.茂木(Mogi，1968)认为，空区是地震应变能的积累区，大地震是当震源区积累的应变能达到极限时，导致岩石突然破裂而发生的.

地震空区用来作为预测地震发生的位置，最初是由Fedotov(1965)提出来的.他把发生在日本—库页岛—勘察加地震带的浅源大地震的震源区画出来后,发现在这些破裂区之间存在多年没有破裂的区域，并且推断这些区域就是未来大地震可能发生的位置.几年后，在他所预报的一些地区先后发生了3次7.8级以上的大地震.

如何确定地震空区是预测板间大地震位置的关键.Kelleher等(1973)首次提出的一个准则是：空区是走滑或逆冲断层地震带的一部分，并且这部分至少30年没有破裂.但有些板块边界地段的历史地震记录并没有逆冲大地震，对于这样的地区，历史地震的平静要么代表一个长期的平稳过程，要么是一种暂时现象，其中大地震的重复时间超过了历史记录的长度.因此，需要区分不同类型的地震空区.后来又对这一标准进行了补充，即把板块边界某些部位有一次或多次大地震的历史记录作为未来大地震在这个部位发生的证据.

预测板块边界地震带大地震位置的地震空区方法，也被用于板块内部的地震(陈章立等，1981；魏光兴等，1981；陆远忠等，1985；梁春涛等，2018).

利用地震空区预测大地震发生位置，方法虽然简单，但是茂木(Mogi，1968)发现，当使用80 年的地震记录时，日本太平洋沿岸对于大于 7 级的事件几乎没有空区.这与使用 30 年得到的结果存在显著差别.因此，判断利用多长时间的地震记录得到的地震空区才是可靠的，至关重要.此外，由于板块边界存在无震滑动区，所以，空区不一定是未来大地震的发生区(Mogi，1968).

为了克服地震空区预测大地震位置的缺欠，早期地表三角测量方法被用来观测大地震前的地表变形，并将计算得到的应变的显著异常区作为可能的大地震位置.1923 年日本关东大地震前几十年就发现东京湾西部的日本佐贺三原地区的构造应变比通常的大好几倍，而在其他地区未发现任何显著的应变(Fujii and Nakane,1979).这个异常是由于亚洲板块和菲律宾海板块间的深部界面的无震逆冲滑动与浅部闭锁段加速下行，在地表产生的挤压导致的.这个方法也被用于研究阿拉斯加Shumagin地震空区的应变积累(Savage and Lisowski,1986).

三角测量方法的优点是通过直接观测地表变形，定量地得到地表应变，但是工作量大，测点稀疏、精度低，受地形影响较大.20世纪70年代全球定位系统(GPS)问世，以其全天候、高精度、自动化、高效益等显著特点，给地表变形监测带来了革命性的变化，使得人们能利用GPS数据直接研究地震的变形场和断层的破裂过程.无论是真实的地震空区还是震前地表应变高异常区，都是由断层闭锁区引起的，因此利用GPS观测数据直接反演地表变形和地震断层的闭锁区，就可能直接预测大地震的发生位置.

断层闭锁程度通常用断层闭锁区 “后滑(backslip)”模型(Savage,1983)反演地表变形得到.这个模型基于断层弹性位错理论.假设俯冲板块沿其长度以汇聚速率均匀地向软流圈滑动，通过在预期的闭锁区施加与俯冲汇聚速率相反的滑动速率，同时在预期的断层闭锁区外的断层指定零位错，就可以实现断层的闭锁.由于稳定状态位移解，理论上在地表不产生应变或地形变化，因此，一般模拟俯冲带断层闭锁，只要在预期的地震断层滑动部位上施加与同震位错反方向的位错(正断层位错)就可以了.这个反位错也称为滑动亏损，其值越大，意味着应变积累越大.这个方法需要预先人为指定反演地震断层闭锁区和板块的运动方向，并且需要对位错方向进行约束.有研究者利用这个方法和GPS数据反演了2011年日本东北大地震前断层的闭锁区(Ito et al.,2000;Nishimura et al.,2004;Suwa et al.,2006;Loveless and Meade,2010),得到滑动亏损深达40～100 km;使用层状黏弹性模型反演得到的滑动亏欠深度在10～40 km，峰值在35 km.这些模型得到的滑动亏损区和震间蠕滑区(靠近海沟)与这次地震的破裂区明显不符.

反演断层闭锁区方法的优点是由于正演使用了地震位错理论模型，位移场有理论解，计算格林函数方便、简单，但是不能很好地反映介质的不均匀性和复杂的断层几何形态.

地震应力模型是一种新的直接用来反演断层应力变化的动力学模型(Xie and Cai,2018)，这个模型由断层和弹性块体组成，考虑到地震的破裂区厚度远远小于断层的长度和宽度，断层被视为两个弹性块体之间的接触面.与断层仅有剪切位错运动的地震位错模型不同，地震应力模型的断层面上不但有剪应力变化，而且还有正应力变化.应力矢量在断层接触面上大小相等，方向相反.断层面上的应力增加和减少分别代表应变能的积累和释放；零剪应力线是断层破裂区的边界.用地震应力模型反演地表变形和断层应力变化,不需要像位错模型反演断层滑动或闭锁那样需要预先人为地指定断层的滑移区或闭锁区，也不需要约束应力变化方向.这个方法已经成功地用来反演2011年日本东北9级大地震的同震断层应力降(图1a)及其震前7年地表变形和断层应力的积累区(Xie et al.,2019).反演得到的3个震前应力积累区，最大的就是主震的破裂区，其余两个分别是主震后30 min内发生在主震破裂区南面和北面的MW7.8和MW7.4大余震区(图1c).

图1 反演得到的断层应力变化(Xie and Cai,2018;Xie et al.,2019)(a)和(b)分别是断层同震剪应力和正应力变化.主震破裂区南部是主震后30 min内发生的M7.8大余震的破裂区,计算得到的主震和M7.8大余震的矩震级与USGS的一致;白色曲线为零剪应力变化等值线.(c)和(d)分别是震前7年(2004—2010年)内断层剪应力和正应力变化.显著的剪应力异常区与同震破裂区一致，紫色线为图1c中同震破裂区边界，主震异常区的南部和北部分别是M7.8和M7.4两个大余震震前的剪应力积累区.Fig.1 Fault stress changes obtained by inversion (Xie and Cai,2018;Xie et al.,2019)(a)and (b)are coseismic shear and normal stress changes on the fault,respectively.The southern part of the main shock rupture zone is the rupture zone of the M7.8 aftershock that occurred within 30 minutes after the main shock.The calculated moment magnitudes of the main shock and the M7.8 aftershock are consistent with those obtained by the USGS.The white curve line is the contour of zero-shear stress change.(c)and (d)are the fault shear and normal stress changes within 7 years (2004—2010)before the main shock,respectively.The significant shear stress anomaly zone is consistent with the coseismic rupture zone.The purple line is the boundary of the coseismic rupture zone in Fig.1c.The shear stress anomaly areas on the south and north sides of the main shock are the shear stress accumulation areas before the two large aftershocks of M7.8 and M7.4,respectively.

本文是2011年日本9级大地震前7年断层应力变化反演工作的继续.为了进一步探讨这次大地震应力积累过程和检验震前7年反演得到的孕震区的稳定性，我们在原来工作的基础上(2004—2010年)，利用日本GNSS数据，反演震前14年到震前8年(1997—2003年)的断层应力变化，同时结合反演得到的前7年断层的应力变化，探讨这次大地震断层应力积累过程和孕震区随时间演化.

1 反演方法

利用日本1997—2003年的GNSS数据,反演断层应力变化所使用的方法与我们提出的方法(Xie and Cai,2018)相同，略述如下.

有限元方法用于正演问题的计算.整个模型的断层面分为有限个子断层，这些子断层都是待反演的未知应力变化区.子断层的数值格林函数解由下面定解问题确定：

Δσji,i=0,

(1)

Δσji=λΔεkkδij+2μΔεji,

(2)

Δεji=(Δuj,i+Δui,j)/2,

(3)

njΔσji|Γ1=Ti,

(4)

Δσji|Γ2=0,

(5)

Δuini|Γ3=0,

(6)

Δσji|Γ3=0,(i≠j).

(7)

上述方程(1)至(3)分别为平衡方程、本构方程和几何方程，其中Δσji,Δεji和Δui分别是地震应力模型内部的应力、应变和位移变化，λ和μ是Lame常数，δij是Kronecker 符号，下角标相同的项表示求和；方程(4)至(7)为模型的边界条件；Ti是法线方向为ni(i=1，2，3)的子断层面上的单位应力分量，Г1,Г2和Г3分别代表模型的断层面、上表面、侧表面和底面边界.

利用有限元方法可以得到上述边值问题的每一个子断层上单位应力分量载荷的位移场，由此可以合成数值格林函数G，进而构造如下的满足GNSS观测数据d(dx，dy，dz)的线性方程组：

(8)

β(Wh-Wf)T=0,

(9a)

λLT=0.

(9b)

(8)和(9)式中，T(Tx，Ty，Tz)是待求的由所有子断层面上应力矢量构成的地震应力模型的断层面上的总应力变化矢量,下标h和f分别代表断层的上盘和下盘；(9)式中的W代表断层面上的法向位移矢量，其中(9a)式和(9b)式分别代表断层法向位移约束条件和解的平滑约束条件(Masterlack,2003)，其中β是约束参数，λ是平滑系数，L是拉普拉斯微分算子.

断层面上的总应力变化矢量T，可以利用带有约束的最小二乘方法，通过构造如下的目标函数(Xie and Cai,2018)：

Sobj(T)=‖dh-GhT‖+‖df+GfT‖

+β2‖(Wh-Wf)T‖+λ2‖LT‖

(10)

求得.

2 结果与讨论

由于本文的研究是以前工作的继续，所以研究区域、地震应力模型和有限元模型与文献(Xie and Cai,2018)相同，如图2所示，模型力学参数和反演中使用的平滑因子、断层法向位移约束参数与我们同震和震前7年反演工作所使用的参数(Xie and Cai,2018;Xie et al.,2019)相同.使用的GNSS数据来自于http:∥www.gsi.go.jp/ENGLISH/page_e30233.html.

图2 (a)2011年日本东北大地震研究区域;(b)建立在弹性理论和有限元方法上的地震应力模型.图中红线是断层线;(c)断层的上盘与下盘,Tr,Ts和Tn表示应力分量变化;(d)断层面子断层划分，字母“B”表示较大的子断层，远离震源Fig.2 (a)Study area of the 2011 great Tohoku-Oki earthquake in Japan;(b)Earthquake stress model based on elastic theory and finite element method.The red line in the picture denotes the fault trace;(c)Hanging wall and foot wall of the fault,Tr,Ts and Tn represent changes of stress components,respectively;(d)Sub-fault division of the fault plane,the letter “B”indicates a larger sub-fault,far away from the source

为方便起见，我们把1997年到2003年反演得到的新结果称为后期结果；把2004年到2010年反演得到的结果称为前期结果.

图3给出了14年间反演得到的后期和前期上盘断层面剪应力的每年变化结果和每期的平均结果.从图中可以看出，震前应力的变化逐年增加，但不是均匀的，其增加区和减小区与2004—2010年应力变化类似.每期平均剪应力增加最显著的地区(图3h和图3p)与图1a中的同震破裂区基本一致.平均剪应力在主震初始破裂点位置(五角星处)1997—2003年为0.03 MPa，2004—2010年为0.04 MPa.平均剪应力结果表明，剪应力积累速率不是均匀的，有变快的趋势.2001年后,在主震的M7.8大余震区的剪应力明显增加,但是在2008年不增反降，与其他年份不同.后来查明，这个地方在2008年发生了一次M6地震.由于此地震的变形场没有从GNSS数据中剔除，所以上述剪应力降低区正是M6地震的应变能释放导致的.这个结果无意中证明了反演方法的可靠性.2007年大余震区震前的剪应力为何也有些下降，是否是2008年M6地震的前震影响，目前原因尚不清楚，需要进一步研究.从图3中剪应力的年平均值可以看出，剪应力明显下降区大致出现在北纬36°到40°，经度在140°到142°，最浅深度在零剪应力变化线处，离海平面约29 km.剪应力降低意味着应变能释放，这个剪应力降低区可以用前震、小的重复地震或者20到60 km深处的加速滑动解释(Wang and Bilek,2014;Mavrommatis et al.,2014;Yokota and Koketsu,2015).

图3 反演得到的断层剪应力变化(a)—(g)1997—2003年每年剪应力变化；(h)7年剪应力变化的平均值；(i)—(o)2004—2010年每年剪应力变化；(p)7年剪应力变化的平均值.Fig.3 Shear stress changes obtained on fault by inversion(a)—(g)Yearly change of shear stress from 1997 to 2003;(h)Average of shear stress change in 7 years;(i)—(o)Yearly change of shear stress from 2004 to 2010;(p)Average of shear stress change in 7 years.

图4给出了反演得到的后期和前期14年间上盘断层面正应力的每年变化结果.从图中可以看出，正应力的主要增加区也基本是稳定的，从1997年开始，在主震破裂区总的趋势是逐渐增加，意味着在这些区域剪切强度增大.主震初始破裂点位置的平均正应力在1997—2003年为0.001 MPa，2004—2010年为0.006.这个结果表明，正应力积累速率也不是均匀的，有变快的趋势，其增加值比剪应力小近一个量级，并且正应力增加区与剪应力的不一致.近南北的零正应力等值线深度(20 km)要浅于图3中零剪应力等值线的深度(29 km).这意味着剪应力增加区包含一部分正应力降低区(零正应力变化等值线下方至零剪切应力变化等值线上方的区域).由于剪切强度与正应力成正比，可以推断在剪应力增加区内存在剪切强度降低区，而在剪应力增加区内正应力增加的区域剪切强度增加.由此可见，零正应力变化等值线的位置就是断层强度降低和增加的分界线.如果孕震区在整个孕震期间是基本稳定的，其有效摩擦系数基本相同，由于断层初始破裂点下方区域的正应力随着时间逐年降低，意味着那里的剪切强度也随时间降低，这就解释了为何这次地震在破裂开始的时候先向断层下方传播的原因(Ide et al.,2011).值得注意的是零正应力变化等值线不但在主震的初始破裂点(地震波定的震源位置)附近穿过，而且还在主震的两个大余震的初始破裂点附近穿过，这一点目前还没有看到有关文献报道.这似乎不能用巧合来解释.这可能是由于这些破裂点基本都处在剪应力增加最显著的地方，在其上方由于正应力增加显著，导致剪切强度明显增加，不利于在那里开始破裂；而在其下方虽然正应力降低，但是幅度较小，因此剪切强度降低不明显，虽然该区处在剪应力增加区，但是其增加不明显，因此在该区开始破裂的条件不如在零应力变化等值线附近.

图4 反演得到的断层正应力变化(a)—(g)1997—2003年每年正应力变化；(h)7年正应力变化的平均值；(i)—(o)2004—2010年每年正应力变化；(p)7年正应力变化的平均值.Fig.4 Normal stress changes obtained on fault by inversion(a)—(g)Yearly change of normal stress from 1997 to 2003;(h)Average of normal stress change in 7 years;(i)—(o)Yearly change of normal stress from 2004 to 2010;(p)Average of normal stress change in 7 years.

图5给出了使用最佳光滑因子和最佳约束系数都为0.05的地表观测位移和反演得到的位移的匹配情况.从图中可以看到，它们符合得很好.这表明反演得到的断层应力变化是可靠的，能够满足观测数据.

图5 符合观测位移的最佳地震应力模型(1997—2003)(a)—(g)水平位移矢量；(h)—(n)垂直位移矢量.红色和黑色箭头分别代表反演和观测得到的位移.Fig.5 Displacements best fitted to the earthquake stress model of the Tohoku-Oki earthquake (1997—2003)(a)—(g)Horizontal displacement vectors;(h)—(n)Vertical displacement vectors.The red and black arrows represent the inverted and observed displacements,respectively.

图6解释了俯冲带逆冲大地震的孕震体与地震空区、地表应变积累区和断层面应力的关系.图中红色A区代表孕震体,应力或应变能增加区，白色B区代表应力减小区或应变能释放区，黄色C区代表地震空区和由三角网得到的应变积累区，A和B间的粗黑线是孕震体的边界，即零剪应力变化等值线，白色区域的黑箭头表示B区的运动方向，蓝色区的黑箭头表示俯冲板块的运动方向；在零剪应力上方的绿线是零正应力变化等值线，其上的圆圈是震源体破裂点在海底的投影(海底震中).图中白色箭头便是剪切位移(与板块尺寸不成比例)，白箭头联线与过断层的垂直白线之间的夹角是剪应变大小，其与该点剪应力成正比.

图6 板块俯冲孕育大地震的示意图Fig.6 Schematic drawing of a large earthquake bred by plate subduction

3 结论

1997—2010年先后两期共14年反演得到的每年断层剪应力显著增加区与断层同震的破裂区结果基本一致，表明利用地震应力模型反演得到的剪应力显著增加区可以用来预测俯冲带大地震的发生位置.在观测数据足够多的地区，预测的可靠性更有保障.反演得到的14年间的主震的孕震体虽然是稳定的，但是孕震体的大小在整个孕震期间，或主震的复发间隔内，是否也是稳定的，本文结果不能回答这个问题，尚需进一步研究.大地震一般孕震时间很长，譬如这次大地震，复发周期长达千年(宍仓正展,2012).如果震源体积是随时间变化的，并且观测资料在时间上不足够长，在空间上不足够密，就不能保证所预测的震源体大小的可靠性，尽管如此，预测的地震大概位置还是有参考价值的.

地震产生的变形不但与断层剪应力有关，而且还与正应力有关.对于倾滑或逆冲型地震，正应力的影响尤为重要.利用地震位错模型得到的闭锁区(滑动亏损区)明显与主震破裂区不符的原因，可能是地震位错模型中没有考虑断层面上法向应力变化对闭锁区深度和大小的影响(Xie et al.,2019)，使得所预测的地震破裂区位置过深.

我们发现零剪应力和零正应力变化的等值线不重合，前者在断层面上的深度大于后者，这意味着在剪应力增加的区域存在强度降低区.震源初始破裂点容易发生在零正应力变化等值线附近.本文结果表明，正应力变化对俯冲带大地震的破裂有影响.

与地震空区方法、地表三角测量寻找应变(二维)异常区方法和反演断层滑动亏损方法预测地震发生地点相比，断层应力反演方法是最直接的动力学方法，抓住了预测大地震位置的关键因素或主要矛盾.随着空间技术的日新月异，这种方法可望在地震动力学机制研究和地震位置预报中发挥重要的作用.