利用球形地球位错理论研究日本MW9.0地震引起的稳态变形与库仑应力变化

陈伟，刘泰，唐河，佘雅文，付广裕

1 中国地震局地震预测重点实验室（地震预测研究所），北京 100036 2 中国科学院大学计算地球动力学重点实验室，北京 100049 3 中国地质大学（北京）地球物理与信息技术学院，北京 100083

0 引言

2011年3月11日，日本东北部海域发生了灾难性的九级大地震，引发了海啸、核泄漏等次生灾害，造成了大量的人员伤亡和财产损失.该地震发生在太平洋板块与北美板块和欧亚板块之间的俯冲带上，断层走向为201°，倾向NWW，倾角约为10°，深度为 28 km.太平洋板块相对于北美板块以83 mm·a-1的速度向西运动（DeMets et al.,2010），沿着日本海沟俯冲到北美板块的鄂霍次克海板块和欧亚板块的阿穆尔板块之下.该地震产生的破裂面约为500 km×200 km（Ammon et al.,2011;Fujii et al.,2011;Hayes,2011；Zhou et al.,2018），引起了数十米的同震滑动，在日本海沟附近甚至达到50 m以上（Yue and Lay,2011;Simons et al.,2011;Shao et al.,2011;Ito et al.,2011;Fujiwara et al.,2011）.与此同时，在东亚大陆包括中国大陆、韩国和俄罗斯等远场地区的GPS（Global Positioning System，全球定位系统）观测中也可以监测到该地震引起的同震形变（王敏等,2011）.

日本MW9.0地震在造成巨大的同震破裂的同时，也引起了显著的震后变形（刘泰等,2017;Zhao et al.,2018）.震后形变主要包括三种机制：震后余滑、黏弹性松弛效应和孔隙弹性回弹（Miyazaki et al.,2004;Ozawa et al.,2011;Pollitz et al.,2006;Jónsson et al.,2003），其中黏弹性松弛效应是由同震破裂导致的地幔物质应力扰动引起的，在震后长时间存在（刘泰等,2017）.

稳态地震变形指的是地震引起的无穷长时间尺度内震后黏弹性松弛效应之和.考虑到黏弹性松弛效应在震后很长一段时间内都将存在，稳态地震变形一般在震后相当长的时间之后才能趋于收敛.由于地球内部地幔的黏弹性构造，大地震引起的总体影响包括同震瞬时发生的弹性破裂和震后长期存在的地幔黏弹性松弛效应.因此，要研究地震在地球内部引起的总体影响时，就可以利用稳态地震变形的方法来解决，计算大地震导致的地球总体震后形变场，同时还可以分析周边区域对该地震的极限响应.Tang和Sun（2019）借助一致性原理研究了黏弹性地球的稳态地震变形问题.Liu等（2020）基于弹性球形地球模型，考虑地球自重、曲率、可压缩性和无限分层结构，通过数值方法研究了地震引起的地球内部同震变形问题，可以计算有限断层地震位错引起的地球内部任意位置同震库仑应力变化.接着，Liu等（2021）基于黏弹性层状球形地球模型，将互易定理计算的地表位错Love数作为初始值，通过从地球表面到地球内部的数值积分得到内部位错Love数，避免了浅地表变形计算时的精度损失问题，可用来计算地球内部任意位置的震后库仑应力变化.上述工作为研究日本MW9.0地震引起的稳态变形与稳态库仑应力变化奠定了理论基础.

本文首先借鉴Tang和Sun（2019）和Liu等（2020,2021）的研究结果，计算了日本MW9.0地震在日本列岛和中国东部地区产生的稳态水平位移，探讨了该强震对周边形变场的综合效应，然后根据地质资料和前人相关研究成果，建立了中国东北和华北地区主要断裂带模型，最后利用球形地球位错理论（Liu et al.,2020,2021）计算了日本MW9.0地震在上述断裂带上引起的同震与稳态库仑应力变化，力争提升相关计算结果的计算精度，并在此基础上分析该强震对中国东北和华北地区地震危险性的影响.

1 日本MW9.0地震断层滑动模型

日本MW9.0地震发生后，不少学者根据大地测量、地震波、海啸等观测数据，建立了一系列断层破裂模型，精度不断提高（Ozawa et al.,2011;刘泰等,2017;等等）.陈伟等（2021）考虑震后形变中包含的有关同震滑动的信息，联合同震和震后GPS数据，引入黏弹性格林函数反演了日本MW9.0地震的断层滑动模型，反演得到的断层滑动模型如图1所示.该模型结果显示，2011年日本MW9.0地震同震破裂的最大值达到了62.72 m，同震滑动的总地震矩为4.48×1022N·m，相应的矩震级为MW9.03.通过对近场和远场形变的模拟，验证了断层滑动模型的可靠性.

图1 日本MW9.0地震断层滑动模型联合同震和震后数据反演得到的结果，黑色五角星表示震中.Fig.1 Slip model of the 2011 Tohoku-Oki MW9.0 earthquakeSlip model inversed by both co-seismic and post-seismic data.The black star indicates the epicenter.

利用上述断层滑动模型，我们分别计算了日本MW9.0地震在日本列岛（近场）引起的同震水平位移和垂直位移，分别如图2、3所示.图中红色箭头表示GPS观测得到的同震水平位移，蓝色箭头表示利用断层滑动模型模拟的同震水平位移，黑色箭头表示GPS观测和模拟结果之间的残差.从图2、3中可以看出，日本本岛北部出现了显著的同震水平位移，最大位移达到了5.24 m，在日本东部沿海地区出现了较为明显的垂直位移，且大部分为沉降，最大位移为1.13 m.通过残差图可以看到模型正演得到的模拟结果与实际观测值符合得很好，我们计算了同震位移GPS观测值与模拟值之间的均方根误差，发现近场同震水平位移和垂直位移的均方根误差分别为1.76 cm和2.55 cm.为了验证断层滑动模型在远场影响研究的适用性，我们还计算了日本地震在中国东北和华北地区引起的同震位移.结果显示，在该区域内也出现了较为明显的东南向同震水平位移（图4），最大值位于中国东北，达到了35 mm，远场同震水平位移的均方根误差值为0.69 mm.

通过对近场和远场的同震位移的计算和验证，可以看出，图1所示的断层滑动模型在研究日本MW9.0地震引起的近场和远场形变机制时都具有较好的精度.因此，利用该模型研究该地震在中国大陆东部地区产生的影响也有很好的适用性.

除了上述原因之外，本文选用陈伟等（2021）给出的断层滑动模型展开研究，还因为该断层破裂模型在构建过程中考虑了日本MW9.0地震引起的震后形变信息，有助于提升该地震引起的震后库仑应力变化计算结果的可靠性.

2 日本MW9.0地震引起的稳态变形

2.1 稳态变形的基本原理

图2 日本MW9.0地震在近场引起的同震水平位移（a） GPS观测到的同震位移；（b） 地震位错理论模拟的同震位移；（c） 观测值与模拟值之间的残差；黑色五角星表示震中.Fig.2 Coseismic horizontal displacements induced by the 2011 Tohoku-Oki MW9.0 earthquake in near field（a） Coseismic displacements observed by GPS;（b） Coseismic displacements simulated by spherical dislocation theory;（c） Residuals between observed and simulated coseismic displacements;The black star indicates the epicenter.

图3 日本MW9.0地震在近场引起的同震垂直位移（a） GPS观测到的同震位移；（b） 断层滑动模型模拟的同震位移；（c） 观测值与模拟值之间的残差.Fig.3 Coseismic vertical displacements induced by the 2011 Tohoku-Oki MW9.0 earthquake in near field（a）,（b） and （c） show the same as with Fig.2.

图4 日本MW9.0地震在远场（中国东北和华北）引起的同震位移（a） GPS观测到的同震位移；（b） 断层滑动模型模拟的同震位移；（c） 观测值与模拟值之间的残差.Fig.4 Coseismic displacements induced by the 2011 Tohoku-Oki MW9.0 earthquake in far field （Northeast and North China）（a）,（b） and （c） show the same as with Fig.2.

当上述作用于地球的内外力是地震时，相应的稳态变形问题就变成了稳态地震变形问题.稳态地震变形指的是地震引起的无穷长时间尺度内震后黏弹性松弛效应之和.考虑到黏弹性松弛效应在震后很长一段时间内都将存在（刘泰等,2017），稳态地震变形在震后相当长的时间之后才能趋于收敛.

为了验证上述稳态变形问题的正确性，我们以PREM（Preliminary Reference Earth Model，初始参考地球模型）为基础，设置了一个黏弹性地球模型（固态内核-液态外核-固态地幔和弹性地壳的地球模型）来计算震后不同时间尺度的位移格林函数，又设置了一个弹性地球模型（固态内核-液态外核和液态地幔-弹性地壳的地球模型）来计算同震和稳态变形的位移格林函数.两个地球模型均设置弹性层厚度为40 km，震源深度为32 km，为了方便后续计算库仑应力变化，我们计算了四种独立位错源在地下10 km深度处引起的同震、震后以及稳态条件下的位移格林函数（Sun et al.,2009）.y12、y32、y22、y33表示四个独立震源，依次为水平走滑位错、垂直倾滑位错、水平引张位错和垂直引张位错.图5给出了四个独立位错源对应的同震、震后以及稳态条件下的位移格林函数在球坐标系下θ坐标轴方向上的分量uθ，也就是南北向分量（Sun et al.,2009）.上述位移格林函数为无量纲的量，实际应用时应在该结果的基础上乘以地球半径的平方.从图中可以看出，震后初期很短时间内的位移格林函数与同震对应的位移格林函数几乎一致，而当时间尺度变大，震后格林函数逐步向稳态的位移格林函数值演化，当时间达到10000年时，震后位移格林函数与稳态格林函数几乎完全重合.由此可见，随着时间的推移，震后位移的格林函数逐渐向稳态结果逼近.图5所示的结果验证了上述稳态地震变形问题的合理性.

2.2 日本MW9.0地震引起的稳态变形

以上述研究思路为基础，在PREM地球模型的基础上，我们构建了一个黏弹性地球模型（固态内核-液态外核-固态地幔和弹性地壳的地球模型）和一个弹性地球模型（固态内核-液态外核和液态地幔-弹性地壳的地球模型），然后结合已有的断层滑动模型，利用球形地球地震位错理论（Liu et al.,2020,2021），分别计算日本MW9.0地震在近场和远场产生的同震和震后以及稳态变形.

为了验证断层滑动模型在计算震后和稳态变形的可靠性，我们首先对日本地震发生后10年内的震后位移进行了模拟验证.为了消除孔隙弹性回弹和震后余滑的影响，只考虑震后黏弹性松弛效应，我们模拟计算了该地震震后5～10年的年平均位移.根据Diao等（2014）、Yamagiwa等（2015）和Zhao等（2018）的研究成果，设定日本周边地区的弹性层厚度为50 km，地幔黏滞性系数为1.0×1019Pa·s，计算了日本地震在近场引起的震后5～10年的年平均位移，如图6所示.从图6中可以看出，利用模型计算得到的震后5～10年的年平均位移与GPS实测值基本符合，震后位移总体观测值与模拟值之间的均方根误差为8.95 mm.因此，利用该模型继续计算震后100年和1000年甚至更久的震后位移也具有一定的可靠性.

由于时间尺度较大，在计算震后位移时，我们只考虑由黏弹性松弛效应引起的震后形变，忽略了震后余滑和孔隙弹性回弹的影响.图7和图8分别表示近场和远场的模拟计算结果，包括同震位移、震后10年、100年和1000年的震后位移以及稳态水平位移.从图7中可以看出，在日本本岛等近场区域，由于同震位移量级较大，近场形变在震后初期比同震形变的量级小很多；当震后时间达到1000年时，由黏弹性松弛效应引起的累计震后位移才接近同震位移量；当地震变形达到稳态时，累积的震后位移可以达到同震位移的两倍，最大可达到8 m.而在中国大陆等远场区域，震后十年累积的位移就已经接近达到同震位移的两倍，随着时间的增加，达到稳态时震后位移大小则与同震位移相差了几十倍；日本MW9.0地震在中国东北地区产生的稳态水平位移最大值超过0.8 m，在华北地区产生的静态水平位移最大值也将达到0.5 m左右.由此可见，日本MW9.0地震引起的震后形变场，在规模上和范围上，均远远超出我们从前的估计和认识.

图5 四个独立位错源对应的格林函数（水平位移分量uθ）（a） 水平走滑位错源 （y12）；（b） 垂直倾滑位错源（y32）；（c） 水平引张位错源（y22）；（d） 垂直引张位错源（y33）.Fig.5 Green′s functions of four independent dislocation sources （Horizontal displacement component （uθ） ）（a） Strike-slip source （y12）;（b） Dip-slip source （y32）;（c） Horizontal tensile source （y22）;（d） Vertical tensile source （y33）.

图6 日本MW9.0地震在近场引起的震后5～10年的年平均位移（a） GPS观测到的震后位移；（b） 震后黏弹性松弛效应模拟的震后位移；（c） 观测值与模拟值之间的残差.Fig.6 Average postseismic displacements per annum of 5～10 years after the 2011 Tohoku-Oki MW9.0 earthquake in near field（a） Postseismic displacements observed by GPS;（b） Postseismic displacements simulated by viscoelastic relaxation;（c） Residuals between observed and simulated postseismic displacements.

为了评估不同断层滑动模型对稳态地震变形的影响，我们分别利用三种不同的断层滑动模型计算了日本地震在近场引起的稳态变形，结果如图9所示.图9a、b、c分别表示利用本研究采用的断层滑动模型、Iinuma等（2012）断层滑动模型和Wei等（2012）断层滑动模型的稳态变形计算结果.同时计算了不同断层滑动模型稳态变形计算结果之间的差异，图9d、e分别表示稳态变形结果b、c与结果a之间的残差.结果表明，不同的断层模型对近场稳态变形影响不大，原因在于这些断层滑动模型都有非常密集的观测数据的约束.

图7 日本MW9.0地震引起的近场同震、震后和稳态变形（a） 同震形变；（b） 震后形变（10年）；（c） 震后形变（100年）；（d） 震后形变（1000年）；（e） 稳态形变（t=∞）.Fig.7 Coseismic,postseismic,and steady-state deformations in near field induced by the 2011 Tohoku-Oki MW9.0 earthquake（a） Coseismic deformation;（b）,（c）,（d） Postseismic deformations of 10,100,and 1000 years,respectively;（e） Steady-state deformation （t=∞）.

图8 日本MW9.0地震引起的远场同震、震后和稳态变形（a） 同震形变；（b） 震后形变（10年）；（c） 震后形变（100年）；（d） 震后形变（1000年）；（e） 稳态形变（t=∞）.Fig.8 Coseismic,postseismic and steady-state deformations in far field （Northeast and North China）（a）,（b）,（c）,（d）,and （e） show the same as with Fig.7.

图9 不同断层滑动模型计算的近场稳态变形及其差异（a） 本研究断层滑动模型计算结果；（b） Iinuma 等（2012）断层滑动模型计算结果；（c） Wei等（2012）断层滑动模型计算结果；（d） 结果b与结果a的残差；（e） 结果c与结果a的残差.Fig.9 Steady-state deformations of different slip models and their residuals（a） Result calculated by the slip model of this paper;（b） Result calculated by the slip model of Iinuma et al.（2012）;（c） Result calculated by the slip model of Wei et al.（2012）;（d） Residual of results b and a;（e） Residual of results c and a.

稳态地震变形考虑无穷时间尺度的震后黏弹性松弛效应的累加，需要经历无限长时间才能达到稳态条件.因此，在实际地震变形过程中，稳态地震变形的条件很难达到.在达到稳态的过程中，地震的震后效应引起的影响会减弱，逐渐湮没在其他构造因素引起的影响之中.所以，在实际地震事件中，稳态地震变形的研究是一种模拟地震引起的总体震后形变场的有效手段，有利于分析周边地区对地震的极限响应，然而很难用实际的观测资料进行解释.

3 中国东北和华北地区主要断裂带模型

对于计算日本MW9.0地震在东北和华北断裂带上产生的库仑应力变化而言，影响计算精度的最重要因素除了所使用的地震位错理论的完备性以外，目标断层的几何与滑动模型的正确性也至关重要.鉴于此，本节依据地质资料以及前人研究成果，详细构建了东北和华北地区主要断裂带的几何与滑动构造模型.

从大地构造上来说，中国东北和华北地区都位于西太平洋俯冲带的西侧.中国东北地区发育了一系列北东向、北西向和近东西向的活动断裂带，包括郯庐断裂带北段（依兰—伊通断裂带和敦化—密山断裂带）、扶余—肇东断裂、嫩江断裂、第二松花江断裂带、勃利—北安断裂和赤峰—开源断裂等主要活动断裂，这些断裂带控制着中国东北地区主要地震的发生（魏光兴等,1993;徐杰等,1998,2000;章振铨等,1999;邓起东等,2002;施炜等,2003;张培震等,2003;杨清福等,2010;刘晓霞等,2012;殷娜等,2019）.中国华北地区位于大别山和秦岭以北，燕山以南，从鄂尔多斯边缘向东延伸至沿海一带，区域内发育了一系列北东向、北西向、近东西向和近南北向的断裂，主要活动断裂带有张家口—渤海构造带、祁吕构造带、洛阳—石家庄—通州构造带、太行山山前构造带、沧东—唐山断裂带、聊城—兰考断裂带、嘉祥断裂带和郯庐断裂带中南段，华北地区主要地震事件几乎都与这些活动断裂带有关（徐杰等,1998,2000;向宏发等,2000;高战武等，2000；邓起东等，2002;龚伟等,2010;丰成君等,2013;张群伟和朱守彪,2019）.

基于地质资料以及前人的研究，根据东北和华北地区大地构造特征及主要断裂带分布情况，我们简化了区域内的地质构造，考虑了其中地震活动丰富和规模较大的断裂带，建立了中国东北和华北地区主要断裂带模型（图10）.各断裂带模型参数及其依据的参考文献如表1所示.

表1 中国东北和华北地区主要断裂带参数Table 1 Parameters of main faults in Northeast and North China

图10 中国东北和华北地区主要断裂带分布图F1—依兰—伊通断裂带；F2—敦化—密山断裂带；F3—扶余—肇东断裂带；F4—嫩江断裂带；F5—第二松花江断裂带；F6—勃利—北安断裂带；F7—赤峰—开源断裂带；F8—张家口—渤海断裂带；F9—祁吕构造带；F10—洛阳—石家庄—通州构造带；F11—太行山山前断裂带；F12—沧东—唐山断裂带；F13—聊城—兰考断裂带；F14—嘉祥断裂带；F15—郯庐断裂带中南段.Fig.10 Distribution of main faults in Northeast and North ChinaF1—Yilan-Yitong fault;F2—Dunhua-Mishan fault;F3—Fuyu-Zhaodong fault;F4—Nengjiang fault;F5—The second Songhuajiang fault;F6—Boli-Bei′an fault;F7—Chifeng-Kaiyuan fault;F8—Zhangjiakou-Bohai fault;F9—Qilü fault;F10—Luoyang-Shijiazhuang-Tongzhou fault;F11—Taihang Mountain fault;F12—Cangdong-Tangshan fault;F13—Liaocheng-Lankao fault;F14—Jiaxiang fault;F15—The middle and south part of Tan-Lu fault.

4 日本MW9.0地震引起的同震与稳态库仑应力变化

日本MW9.0地震的发生释放了巨大的能量，必然会造成周边地区构造应力环境的变化，对地震危险性产生影响.该地震发生后，许多学者都对该地震在中国大陆产生的影响进行了研究.王敏等（2011）利用插值方法从同震位移场中计算出中国东北和华北地区的水平应变场，发现该地震在郯庐断裂带不同段产生了不同程度的张性应变，其中，在郯庐断裂带北段即依兰—伊通断裂带和敦化—密山断裂带上的张性应变最大，而在该段加载的静态同震库仑应力最大不超过2 kPa，认为该地震对郯庐断裂带北段的地震活动产生了一定的影响，但对郯庐断裂带其他段和中国东部地区其他断裂带影响很小.丰成君等（2013）研究了该地震对中国东北和华北大陆主要活动断裂带的影响，发现在区域内主要活动构造带上的库仑破裂应力变化大部分为正值，但均不超过1 kPa，认为日本大地震对中国东北和华北地区地震活动起一定的加速触发作用，但程度有限.Cheng等（2014）计算了该地震在中国东部地区四条主要活动断裂带上引起的库仑应力变化，发现同震和震后库仑应力变化值均不超过2 kPa，小于有效触发地震的阈值10 kPa（Stein,1999;Yu et al.,2006）.Shao等（2016）研究了日本MW9.0地震对中国东北和华北地区主要断裂带上的影响，发现在大部分北东向的断裂带上产生的库仑应力变化为正值，在北西向、东西向和少部分北东向断裂带上的库仑应力变化为负值.Yu等（2016）研究了中国东北地区主要断裂带上的同震库仑应力变化，发现都小于触发地震的阈值，同时Yu等（2016）统计了中国东北地区2008—2013年之间发生的1级以上的地震，发现依兰—伊通断裂带南段和北段的地震活动性增加，而中段地震活动性下降.Chen等（2020）根据球形地球位错理论（Liu et al.,2020）计算了郯庐断裂带北段的同震库仑应力变化，发现日本地震在依兰—伊通断裂带上产生的影响要大于在敦化—密山断裂带上产生的影响，但库仑应力变化的最大值均不超过3 kPa，证明该地震并不会改变郯庐断裂带北段的应力状态和地震活动性.除Chen等（2020）以外，上述研究均基于平面半空间地震位错理论（Okada,1985;Wang,2006）展开研究，考虑到中国东部地区与日本MW9.0地震之间的距离超过了1500 km，地球曲率的影响不能被忽略，上述计算结果均有可能包含着较为显著的计算误差，有待利用球形地球位错理论展开研究，提升计算结果的准确性.

鉴于此，本文基于中国东北和华北地区主要断裂带模型（表1）和日本MW9.0地震断层破裂模型（图1），利用球形地球位错理论（Liu et al.,2020,2021）计算了该地震在中国东北和华北地区主要断裂带上引起的库仑应力变化.本文利用公式（1）计算库仑应力变化：

ΔCFS=Δτ-μ′Δσn，

（1）

其中，ΔCFS为库仑应力变化值，Δτ和Δσn为断层上剪应力和正应力的变化值，μ′为有效摩擦系数.有效摩擦系数μ′在不同的研究中有差异，通常在0.2～0.8之间取值.Cheng等（2014）和Yu等（2016）在研究日本MW9.0地震对中国东部和东北引起的影响时，都将0.4作为有效摩擦系数μ′的取值，该值也是断层库仑应力变化研究中最常用的有效摩擦系数（King et al.,1994）.因此，本研究将有效摩擦系数μ′取为0.4.

同样的，在PREM地球模型的基础上，我们构建了一个黏弹性地球模型（固态内核-液态外核-固态地幔和弹性地壳的地球模型）和一个弹性地球模型（固态内核-液态外核和液态地幔-弹性地壳的地球模型），然后结合表1给出的断裂带的几何与滑动构造模型，利用球形地球地震位错理论（Liu et al.,2020,2021）计算日本MW9.0地震在中国东北和华北地区主要断裂带上引起的同震与稳态库仑应力变化.计算过程中，断层面上的变形深度设为10 km.计算结果如图11所示.

图11 日本MW9.0地震在中国东北和华北地区主要断裂带上引起的库仑应力变化（a）和（b）分别为同震和稳态库仑应力变化计算结果.Fig.11 Coulomb failure stress changes caused by the 2011 Tohoku-Oki MW9.0 earthquake on the main faults in Northeast and North ChinaSubfigure （a） and （b） show the coseismic and steady-state Coulomb failure stress changes,respectively.

日本MW9.0地震在中国东北和华北地区大部分北东向分布的断裂带上引起的同震库仑应力变化为正值，而在近东西向分布的断裂带上的库仑应力变化多为负值.由图11a可知，同震库仑应力变化在依兰—伊通断裂带北段为正值，最大值约为0.6 kPa，而在南段则为负值；敦化—密山断裂带上的同震库仑应力变化值均为正值，最大值达到了0.8 kPa；在嫩江断裂带上引起的同震库仑应力变化值也是正值，在0.5～0.6 kPa之间；在第二松花江断裂带上的同震库仑应力变化值差异较大，在0.1～0.6 kPa之间变化；在扶余—肇东断裂带和勃利—北安断裂带上的同震库仑应力变化为负值；在赤峰—开源断裂带上的同震库仑应力变化值较小，不超过0.1 kPa；在华北地区大部分北东向的断裂带包括祁吕构造带、洛阳—石家庄—通州构造带、太行山山前断裂带、沧东—唐山断裂带、以及聊城—兰考断裂带南段上均为正值，但数值较小，都在0.3 kPa以内；在张家口—渤海断裂带和郯庐断裂带中南段为负值.总体上，日本MW9.0地震在中国东北地区主要断裂带上引起的同震库仑应力变化大于在华北地区产生的同震库仑应力变化，故该地震对中国东北地区的影响要大于华北地区.日本MW9.0地震在中国东北和华北地区引起的同震库仑应力变化最大值接近1 kPa，小于有效触发地震的阈值10 kPa，这与前人的研究结果基本一致，表明日本MW9.0地震在中国大陆东北和华北地区引起同震影响有限，不会大幅度改变区内主要断裂带的应力状态和地震活动性.

与此同时，我们计算了日本MW9.0地震震后黏弹性松弛效应在中国东北和华北地区引起的库仑应力变化，获取了稳态库仑应力变化.从图11b中可以看出，在依兰—伊通断裂带、敦化—密山断裂带和嫩江断裂带上的库仑应力变化值较大，最大值分别为11.8 kPa、8.0 kPa和10.5 kPa.在华北地区最大的库仑应力变化位于太行山山前断裂带北段，达到7.7 kPa.经比较可知，稳态库仑应力变化可以达到同震库仑应力变化的10倍以上.但是，即便达到稳态，也只在依兰—伊通断裂带北段和嫩江断裂带南段会出现超过10 kPa的库仑应力变化，在其他断裂带上均未出现超过10 kPa这个阈值的库仑应力变化.因此我们得出结论，日本MW9.0地震对我国东北和华北地区地震危险性的总体影响非常有限.

5 讨论

在计算日本地震在近场引起的震后位移（图6）时，我们采用的是均一的地幔黏滞性系数，忽略了地幔黏滞性结构的横向差异.结果发现相比于同震位移，震后位移的实测值与模拟值之间的符合效果稍差.众多前人的研究表明，地幔黏滞性结构的横向差异对震后黏弹性松弛效应有着不可忽略的影响，这也很可能是造成震后位移观测值与模拟值之间差异的主要原因.Hu等（2016）利用三维球形黏弹性有限元模型研究了日本MW9.0地震的震后变形，并对日本地震震源及周边地区的地幔黏滞性结构进行了分析.使用的三维模型中包含弹性的上板块和俯冲板以及黏弹性的地幔楔、海洋地幔和沿板块边界的剪切带，在一个2 km厚的弱剪切带中模拟应力驱动的震后余滑，用瞬态的Burgers流变体来模拟上地幔的黏弹性松弛和剪切带的变形.利用重复地震的余滑估计约束确定50 km以内深度的剪切带的最佳黏滞性系数为1017Pa·s，深部剪切带、地幔楔和海洋上地幔的黏滞性系数分别为5×1017Pa·s、3×1019Pa·s和5×1019Pa·s.梁明等（2018）也对日本MW9.0地震的震后形变进行了研究，利用了GPS和GRACE（Gravity Recovery and Climate Experiment，重力恢复与气候实验）重力卫星观测数据，并对日本地震震源周边地区的地幔黏滞性结构进行了反演.结果表明，震源区的地幔黏滞性结构并不是均一的，在纵向和横向上都存在差异.日本地震断层破裂两侧区域的流变性质也存在差异，位于大陆一侧的地幔黏滞性系数在1.0×1019Pa·s量级，而位于海洋一侧的介质物质的黏滞性系数比大陆一侧略小，约为6.0×1018Pa·s.由此可见，日本地震震源及周边地区地幔黏滞性结构并不是均一的，而日本地震震源及周边地区的地幔黏滞性结构对震后形变有着很大的影响.因此，在后续的研究中需要考虑地幔黏滞性结构的不均匀性，使得模型对震后形变的拟合效果更好.

在计算稳态地震变形时，我们考虑的实际上是无穷时间尺度内震后黏弹性松弛效应的累加.在这里，我们只考虑了2011年日本MW9.0地震的影响，利用反演得到的断层滑动模型计算了该地震在近场和远场引起的稳态地震变形.然而事实上，日本俯冲带是太平洋板块向欧亚板块和北美板块俯冲的区域，是世界上地震活动最丰富的地区之一，该地区强震的孕震周期约为几百年，这远小于震后黏弹性松弛达到稳态所需要的时间.所以在进行稳态地震形变研究时，在达到稳态的过程中该地区就已经发生了很多地震，这就需要考虑这些地震的综合影响.因此，在后续稳态地震变形的研究中，我们需要考虑多个地震的叠加效应，计算这些地震产生的总体地震形变场，进而进行更加深入的研究，包括分析其对区域速度场的影响，考虑其对板块漂移的影响等等.

6 结论

稳态地震变形是指由地震引起的、超长时间尺度内震后黏弹性松弛效应的累加，一般在震后相当长的时间之后才能趋于收敛.通过稳态地震变形，可以计算大地震导致的地球总体震后形变场，有利于分析周边地区对该地震的极限响应.本文首先借鉴Tang和Sun（2019）和Liu等（2020,2021）的研究结果，利用球形地球地震位错理论，高精度计算了日本MW9.0地震在日本列岛和中国东部地区产生的稳态水平位移.结果表明，日本MW9.0地震在日本本岛引起的稳态水平位移最大可达到8 m，大体上是同震水平位移的两倍；在中国东北地区产生的稳态水平位移最大值超过0.8 m，在华北地区产生的静态水平位移最大为0.5 m左右，是同震水平位移的30～40倍.如此规模的日本MW9.0地震震后形变场，远超出我们从前的估计，刷新了我们对该地震影响的认知.

接着，本文根据地质资料和前人研究成果，建立了中国东北和华北地区15条断裂带的几何与滑动构造模型，最后利用最新的球形地球位错理论（Liu et al.,2020,2021）计算了日本MW9.0地震在上述断裂带上引起的同震与稳态库仑应力变化.结果表明，该强震在中国东北和华北地区主要断裂带上引起的同震库仑应力变化值均不超过1 kPa.其中，在依兰—伊通断裂带北段、敦化—密山断裂带和嫩江断裂带上引起的同震库仑应力变化最大，大体在0.5～0.8 kPa之间.日本MW9.0地震引起的稳态库仑应力变化在依兰—伊通断裂带北段和嫩江断裂带南段的一部分出现超过10 kPa的库仑应力变化，在其他断裂带上均未超过10 kPa.稳态库仑应力变化显著超过同震结果，但均难以达到10 kPa这个地震触发阈值（Stein,1999）.因此，日本MW9.0地震对中国东北和华北地区地震危险性影响有限，不会显著改变区内主要断裂带的应力状态和地震活动性.

致谢本文采用的GPS数据来自日本国土地理院和美国内华达大学大地测量实验室.感谢两位匿名审稿人提出了宝贵的修改意见.本文使用GMT软件绘图，在此一并表示感谢.