从同震和震后形变分析1668年M8.5级郯城地震对周边地震活动性的影响

徐 昊，孙玉军，吴中海

(中国地质科学院地质力学研究所，北京 100081)



从同震和震后形变分析1668年M8.5级郯城地震对周边地震活动性的影响

徐 昊，孙玉军，吴中海

(中国地质科学院地质力学研究所，北京 100081)

以山东郯城1668年大地震为例，以前人地表地质调查结果为约束，利用弹性位错理论初步获取了该地震的同震破裂模型；在此基础上，基于粘弹性分层模型分析了该地震的同震和震后形变，同时以主震断层为接收断层计算了库仑应力分布，进一步讨论了地幔不同粘滞性系数对地表形变和库仑应力变化的影响。计算结果显示，该地震是一个右旋走滑为主兼有一定逆冲性质的地震，其同震位移巨大，能量释放较彻底；同震破裂造成震中郯城县西北、东北和南部部分断层库仑应力增加，而震后形变使得这些断层库仑应力进一步增加，在单县、宿迁和日照等地，地震后350 a库仑应力变化量达到+1bar—+1MPa量级；地幔粘滞性系数不同，形变量和库仑应力变化达到稳定的时间不同，但最终趋于稳定的数值基本一致。

郯城大地震；同震和震后形变；库仑应力变化；粘滞性

0 引言

郯庐断裂带是中国东部一条重要的活动断裂带，历史上沿该断裂带发生过多次大地震。1668年7月25日晚，在山东南部发生了一次旷古未有的特大地震，根据历史资料得到该地震的震级约为Ms8.5级，极震区位于山东省郯城、临沭、临沂交界处，震中位置(34.8°N，118.5°E)，极震区烈度达Ⅻ级[1]。由于极震区大部分位于郯城县境内，故称为郯城地震。这是我国大陆东部板块内部最强烈的地震，造成巨大财产损失和人员伤亡。

关于1668郯城地震，我国学者已经对其进行了相关研究，分析了这次地震发震断层的构造特征和破裂机制[2～3]，采用不同手段得到了此次地震震源参数[4～6]，估算了郯城8.5级地震的复发周期[7～8]。但关于此次地震造成的应力场变化研究很少。近年来，越来越多的研究表明，地震发生会随之释放断层上积聚的应力，但应力并不会凭空消失，部分应力转移到其他地区，导致应力积聚，进而诱发后续地震，此即地震的应力触发理论[9]。根据应力触发理论，库仑应力增强相当于断层额外负荷的增加，容易诱发后续地震；反之，库仑应力减弱的区域——应力影区，断层负荷部分卸载，地震发生的概率降低。1668年郯城8.5级特大地震的同震及震后应力特征如何？其造成的应力场变化对周围地区的地震活动性有何影响？本文围绕这些问题，利用弹性位错理论和分层岩石圈模型，计算了郯城1668年地震同震和震后形变以及库仑应力的变化，探讨本次地震对郯庐断裂带及附近区域地震活动性的影响，以期为进一步研究该区域的地震活动性提供参考。

1 区域构造背景

地震的发生与活动断层联系密切，断层活动是地震产生的原因。郯庐断裂带是中国大陆东部一条北北东向的深大断裂，也是中国东部有巨大影响的强震活动带，史料记载曾发生过公元前70年安丘7级地震、1597年渤海7级地震、1668年郯城8.5级地震、1888年渤海7.5级地震、1969年渤海7.4级地震、1975年海城7.3级地震等[10]。其中1668年郯城8.5级地震的地震断层为一高角度右旋走滑逆断层，呈北北东向展布(见图1)，断层面向东倾[6]。

图1 研究区构造背景图(椭圆内加粗断裂为此次地震主震断层)Fig.1 Tectonic background map of Tanlu fault and its surrounding regions

2 计算方法

2.1 分层模型

用于地震研究的弹性位错模型主要可分为3类，第一类主要基于均匀弹性半空间介质模型[11]，第二类主要基于分层介质模型[12]，第三类基于成层球对称地球模型的球形位错模型[13]。本文基于地壳的分层模型采用格林函数解求解这次地震引起的同震位移以及引起的应力场变化。震后很长时间，有许多影响震后应力变化的因素，如震后余滑、孔隙效应和粘性松弛等，其中震后粘性松弛效应对震后形变和应力起到非常重要的作用[14]。Wang[15]提出了利用正交归一法计算地震应力场Green函数方法，并在此基础上建立了粘弹松弛分层模型下的地震同震及震后形变模型，编写了数值计算程序——PSGRN/PSCMP软件，用该程序基于分层、重力作用下的弹性模型来研究由于断层错动引起的同震和震后形变。假设郯城及周边地区震后弹性松弛的过程主要发生在30 km上地壳以下。震后粘弹松弛应力变化是指一定深度下，受高温高压环境的影响，岩石层的力学性质逐渐由脆性向粘塑性转化。由于地震发生时破裂的速度很快，在短时间内粘弹性的中下地壳以及上地幔力学性质表现为弹性体，然后震后随时间推移，中下地壳和上地幔所积累的应力和应变逐渐释放，并向上传输到弹性的上地壳中，导致上地壳中应力状态发生改变[16]。对于地壳和上地幔来说，粘滞性系数不同，对震后形变和库仑应力变化影响不同。粘滞性系数越低，松弛时间越短，震后短期内形变和库仑应力增加的都比较快，很快能达到稳定状态；相反，粘滞性系数越高，震后短期形变和库仑应力增加的都比较慢，但最终的稳定值相对较高[17]。模型假定地壳厚度为30 km，视为完全弹性体，粘滞性系数取无穷大；研究区缺乏长期的形变观测资料和相关的震后形变研究成果，模型复杂化需要更多的参数来约束，因此采用简单的Maxwell体分析地幔粘弹性变形。由于未获取研究区地幔准确粘滞性系数，本文模型中计算所需的粘滞性系数取1.0×1018Pa·s，最后会分析不同地幔粘滞性系数条件下，计算得到的形变量和库仑应力的差异，并结合李勇红等[5]给出的郯城8.5级地震震源区附近地壳上地幔速度结构，给出最终的分层模型及参数(见表1)。

表1 分层模型及参数

2.2 库仑应力变化的计算

随着断层上应力不断累积，最终剪应力能够克服断裂上的静摩擦力，锁闭的断层会发生破裂，断层的破裂符合库伦破裂准则[18]。利用弹性位错理论，可以计算地震同震破裂所引起的应力场变化[19]。所研究的接收断层库仑应力变化为：

(1)

式中：Δσf为计算断层面上库仑应力的变化，Pa；Δτs是断层面上剪应力的变化，Pa；μ′为有效摩擦系数；Δσn是计算断层面上正应力的变化，Pa。关于μ′的取值，不同学者有所差别。King等[20]认为有效摩擦系数对库仑应力变化的结果将产生中等大小的影响，本文在计算过程中取中间值0.4。

通过主震断层模型计算静态库仑应力不但能够有效反映余震的发展态势，同时也能够获取对周围活动断层的影响[17]。由于此次地震发生时间距今已有350 a，本文仅选取一条主震断层计算库仑应力。采用之前学者研究所获得的断层相关参数，走向21.57°、倾角89.51°、滑移角148.9°；起始破裂坐标为34.5°N、118.35°E，断层面起始深度及其实际宽度取0 km和50 km[21]，断层破裂长度为150 km[6]。

3 计算结果与分析

3.1 同震和震后形变

1668年郯城8.5级主震断层为一高角度右旋走滑兼有逆冲性质的断层，呈北北东向展布，断层面东倾，震后观察到水平最大形变量和垂直最大形变量分别为10 m和3 m[6]。尽管计算过程是在半无限空间假设模型下得到的，未考虑相关地形的影响，但通过图像反映出的形变，依然可以判断主震断层为一右旋走滑逆断层。计算得到同震时刻地表水平变形量和350 a后地表水平变形量(见图2a，2c)，同震时刻地表垂直变形量和350 a后地表垂直变形量(见图2b，2d)，其中350 a以后水平方向和垂直方向变形量分别为8.33 m和2.73 m，与王华林等[6]观测数据存在一定误差，由于在模型构建中缺少断层面上精准的滑移数据，因此认为断层模型依然能够在一定程度上反映当时的地震破裂变形。

a—同震时刻地表水平方向形变量；b—同震时刻地表垂直方向形变量(向上为正)；c—350 a后地表水平方向形变量；d—350 a后地表垂直方向形变量(向上为正)图2 同震及震后0 km处地表形变量Fig.2 The coseismic and postseismic deformation of the surface

对比地表处水平方向的形变可以发现，同震时刻和350 a以后断层右旋的性质没有改变，但350 a后水平形变量稍大于同震时刻水平形变量，距离主震断层越远增量越大。对比地表处垂直方向的形变可以发现，同震时刻和350 a以后断层逆冲的性质也未改变，但局部影区形变方向发生了倒转，以右下角区域为例，同震时刻位移方向垂直向上，350 a后则向下变形。推测以上两种变化原因与地壳岩石粘滞性有关：同震时刻断层两盘发生水平方向错动，从而使距离破裂带较近的区域瞬间获得较大的形变量，而震后距破裂带较远的区域水平形变在岩石粘滞性作用下缓慢增加；已知主震断层滑移角为149.5°，同震时刻断层上下盘沿着走向错动，必然会引起所在盘沿滑移方向一侧的抬升，因此同震时刻右下和左上形变均为向上方向；随着时间的推移，由于岩石粘滞性影响，变形逐渐恢复甚至引起反向形变，从而使图中原先向上方向的变形转变为向下变形。关于这种变化的量化分析，会在后来的部分进一步论述。

3.2 库仑应力变化

每一次大地震发生都会引起周围区域的应力增加或降低。一般来说，如果一个断层的库仑应力增加，则该断层发生地震的危险性增强，反之则降低。同震破裂造成震中郯城县西北、东北和南部部分地区库仑应力增加，而震后形变使得这些地区库仑应力进一步增加，在单县、宿迁和日照等地，地震350 a后库仑应力变化量达到+1 bar—+1 MPa量级(见图3)。观察研究区内震后350 a间Ms1—Ms7地震事件发现，地震主要发生在这些库仑应力增加的区域内。虽然地震的发生与构造运动直接相关，但这样的特大型地震发生后对周围断层的应力积累影响十分明显，而且计算结果显示震后350 a该影响依然存在，只是目前已经趋于稳定。

同震库仑应力增加最明显区域主要位于断层破裂带的两端及同震破裂不完全的区域。断层破裂带应力释放明显，而破裂带的两端则为应力集中的区域。本次地震的震源深度为20 km左右，因此，20 km深度同震应力释放较地表同震应力释放更为明显(见图3a，3b)，沿着破裂带主要为应力影区；而在地表有部分滑动位移小，破裂不完全，应力增加的区域(见图3a，3c)。对比同震时刻和350 a的同震库仑应力分布，可以看出震后350 a在地表和20 km处应力增加的区域在面积和数值上都明显增大。将本次地震后该区域的所有地震投影在平面后发现，该区域的大部分地震事件也主要位于应力增加的区域。

为了更清楚地观察各应力影区的形变和库仑应力变化趋势，笔者选取南京、宿迁、单县和日照4座城市计算并绘制其在深度20 km处350 a间位移和库仑应力变化曲线。南京和日照在地震发生后库仑应力都减少，不同的是南京在地震发生约20 a后库仑应力逐渐增加，至今增加到-25 kPa左右(见图4a)；而日照在地震反生后库仑应力稍微减少，随后一直增加，到后来变化值稳定在+1000 kPa左右(见图4b)。宿迁和单县则在地震发生后库仑应力一直增加，最终宿迁库仑应力增加稳定在+1300 kPa左右(见图4c)，而单县仅增加了+150 kPa左右(见图4d)。库仑应力增量表明宿迁和日照等地与郯庐断裂带主震断层相同性质的断层地震危险性增强。

考虑同震和震后粘弹松弛效应对应力场变化的影响时，地幔的粘滞性系数不同取值可能会改变研究区的应力场分布。由于U(x)反映南北方向形变量，断层走向为21.57°且为右旋断层，U(x)的值能清楚反映断层的运动性质。本文选取不同地幔粘滞性系数计算了郯城地震发生后350 a内，南京市深度20 km处的U(x)和库仑应力变化值(见图5)，结果表明，震后形变量和库仑应力变化量初始值相同，受粘滞性系数影响较小。随着时间的推移，不管粘滞性大小，最终形变量稳定值都近似相同，均为-1 m，但粘滞性系数越小，U(x)越快达到稳定值，当粘滞性系数为1×1017Pa·s时，10 a左右就可达到稳定值，当粘滞性系数为7×1018Pa·s时，则需要300 a左右才能达到稳定值(见图5a)。对于库仑应力变化来说，震后初值变化也都相同，均为-17.5 kPa，曲线变化模式类似，都是库仑应力先减少后增加最后趋于稳定；当粘滞性系数为1×1017Pa·s时，50 a左右即可达到稳定值，约为-24 kPa；当粘滞性系数为1×1018Pa·s时，150 a左右达到稳定，约为-22.5 kPa；当粘滞性系数为7×1018Pa·s时，350 a时仍未达到稳定状态，但可以估计其达到稳定状态的值基本相同(见图5b)。

a—地表处同震库仑应力变化分布；b—20 km深处同震库仑应力变化分布；c—地表处震后350 a库仑应力变化分布；d—20 km深处震后350 a库仑应力变化分布图3 同震库仑应力变化分布图(灰色圆圈为1668年郯城地震后至今研究区内发生的M1—M7地震事件)Fig.3 The distribution map of coseismic Coulomb stress change

a—南京市；b—日照市；c—宿迁市；d—单县市△CFS—库仑应力变化量；U(x)—南北向位移(北方向为正)；U(y)—东西向位移(东方向为正)；U(z)—垂向位移(向上为正)图4 研究区内不同城市震后350 a内形变量和库仑应力变化Fig.4 Deformation and Coulomb stress change in different cities in study regions 350 years after the earthquake

图5 不同地幔粘滞性系数对应的震后形变和库仑应力变化(南京市)Fig.5 The postseismic deformation and Coulomb stress change of different mantle viscosities(Nanjing)

4 结论

本文基于粘弹松弛分层模型下的地震同震及震后形变模型，采用数值计算程序——PSGRN/PSCMP软件，利用前人地表观测和反演的地震数据计算了主震断层同震以及震后的地表形变和库仑应力变化量。

通过对比同震和350 a后震中郯城县东南部分位移形变量，从垂直向上的位移量转变为向下位移量，可知震后最终位移方向在受地壳粘弹性效应的影响下会发生倒转。

以主震断层为接收断层计算了同震和震后库仑应力。同震库仑应力增加的区域主要位于断层破裂带的两端及同震破裂不完全的区域，具体分布于震中郯城县西北、东北和东南的部分区域，震后形变使得这些地区库仑应力进一步增加，在单县、宿迁和日照等地，地震350 a后库仑应力变化量达到+1 bar—+1 MPa量级。其他地区库仑应力降低，南京地震350 a后库仑应力变化量为-25 kPa。震后该区域发生的地震主要分布于库仑应力增加的区域。

通过对比地幔不同粘滞性系数对震后形变和库仑应力的影响，发现对于同一区域，粘滞性系数越小，震后形变量和库仑应力变化达到稳定的时间越短，反之则时间越长；不同的地幔粘滞性系数对最终稳定的形变量和库仑应力影响不大。

在讨论此次地震与区域内周边地震活动性关系时，主要依据库仑应力变化量，虽然地震的发生与构造运动直接相关，但这样的特大型地震发生后对周围断层的应力积累影响十分明显，而且计算结果显示震后350 a后该影响依然存在，但目前已经趋于稳定，实际地震的发生还要综合考虑构造活动等其他因素。此外由于1668—1970年之间区域地震事件缺少深度数据，所以在考虑后续地震在0 km和20 km处分布与库仑应力变化的关系时误差较大。

[1] 马玉香，钟普裕．1668年山东郯城8.5级地震综述[J]．国际地震动态，2009，(2)：9～18.

MA Yu-xiang,ZHONG Pu-yu.Summarization of the Tancheng Earthquake with 8.5 in 1668[J]. Recent developments in World Seismology, 2009, (2): 9～18.

[2] 高维明，郑朗荪，李家灵，等．1668年郯城8.5级地震的发震构造[J]．中国地震，1988，4(3)：9～15.

GAO Wei-ming,ZHENG Lang-sun,LI Jia-ling.The triggering structure of theMs8.5 earthquake of Tan Cheng in 1668[J].Earthquake Research in China,1988, 4(3): 9～15.

[3] 鞠林雪．1668年郯城8.5级地震发震断层构造特征与应力场分析[D]．合肥：合肥工业大学，2012.

JU Lin-xue.The research on the structure of triggring faults and stress field in 1668: The Tancheng earthquake[D]. Hefei: Hefei University of Technology, 2012.

[4] 李清河，张元生，鲍海英，等．据地壳速度结构推测1668年山东郯城8.5级大地震震中[J]．中国地震，2014，30(1)：30～42.

LI Qing-he, ZHANG Yuan-sheng, BAO Hai-ying, et al. Crust veolcity structure of the epicenter ofM8.5 earthquake at Tancheng, Shandong, China, in 1668[J]. Earthquake Research in China, 2014, 30 (1): 30～42.

[5] 李永红，胡新亮，许萍，等．以双差定位方法对郯城8.5级地震震中附近现代小震重新定位[J]．西北地震学报，2011，33(1)：71～75.

LI Yong-hong, HU Xin-liang, XU Ping, et al. Relocation of modern small earthquakes near the epicenter of Tancheng 8.5 earthquake using the Double Difference Earthquaake Location Algorithm[J]. Northwestern Seismological Journal, 2011, 33(1): 71～75.

[6] 王华林，耿杰．1668年郯城8.5级地震震源参数及其讨论[J]．地震学刊，1996，(4)：27～33.

WANG Hua-lin, GENG Jie. Discussion about focus parameters of Tancheng earthquake ofMs8.5 in 1668[J]. Journal of Seismology, 1996, (4): 27～33.

[7] 王华林．1668年郯城8.5级地震断裂的全新世滑动速率，古地震和强震复发周期[J]．西北地震学报，1995，17(4)：1～12，22.

WANG Hua-lin. Holocene displacement rate,paleoearthquakes and recurrence intervals of strong earthquakes along the 1668 Tancheng Earthquake (Ms=8.5) Fault[J]. Northwestern Seismological Journal, 1995, 17(4): 1～12.

[8] 王启鸣．三维有限元方法估算郯城8.5级大震重复周期[J]．中国地震，1988，4(3)：54～58.

WANG Qi-ming. The calculation of recurrence intervals of theMs8.5 earthquake in Tancheng by the three dimension finite element method[J]. Earthquake Research in China, 1988, 4(3): 54～58.

[9] 万永革，吴忠良，周公威，等．地震应力触发研究[J]．地震学报，2002，24(5)：533～551.

WAN Yong-ge, WU Zhong-liang, Zhou Gong-wei, et al. Research on seismic stress triggering[J]. Acta Seismologica Sinica, 2002, 24(5): 533～551.

[10] 魏光兴．郯庐带地震活动性研究[M]．北京：地震出版社，1993.

WEI Guang-xing. Research on the seismic activity of Tanlu earthquake belt[M]. Beijing: Seismological Press, 1993.

[11] Okada Y. Surface deformation to shear and tensile faults in a halfspace[J]. Bulletin of the Seismological Society of America, 1985, 75 (4): 1135～1154.

[12] Pollitz F F. Postseismic relaxation theory on a spherical Earth[J]. Bulletin of the Seismological Society of America, 1992, 82(1): 422～453.

[13] 付广裕，孙文科．球体位错理论计算程序的总体设计与具体实现[J]．地震，2012，32(2)：73～87.

FU Guang-yu, SUN Wen-ke. Overall design and special strctures of the computing codes for coseismic deformations on a layered sperical earth[J]. Earthquake, 2012, 32(2): 73～87.

[14] 石耀霖．关于应力触发和应力影概念在地震预报中应用的一些思考[J]．地震，2001, 21(3)：1～7.

SHI Yao-lin. Sress triggers and stress shadows How to apply these concepts to earthquake prediction[J]. Earthquake, 2001, 21(3): 1～7.

[15] Wang R, Lorenzo-Martín F, Roth F. PSGRN/PSCMP: A new code for calculating co-and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory[J]. Computers & Geosciences, 2006, 32: 527～541.

[16] Freed A M. Earthquake triggering by static, dynamic, and postseismic stress transfer[J]. Annual Review of Earth & Planetary Sciences, 2005, 33: 335～367.

[17] 孙玉军，董树文，范桃园，等．从同震和震后形变分析日本东北Mw9.0级大地震对近场地震活动性的影响[J]．地球物理学进展，2013，28(3)：1131～1139.

SUN Yu-jun, DONG Shu-wen, FAN Tao-yuan, et al. The effect of TohotuMw9.0 earthquake on the near field seismic activity from coseismic and postseismic deformation[J]. Progress in Geophys, 2013, 28(3): 1131～1139.

[18] 江国焰．利用库仑破裂准则评估活动断层地震危险性[D]．武汉：武汉大学，2013.

JIANG Guo-yan. Assessing the seismic hazards of active faults using Coulomb Failure Criterion[D]. Wuhan: Wuhan University, 2013.

[19] Okada Y. Internal deformation due to shear and tensile fault in a half space[J]. Bulletin of the Seismological Society of America, 1992, 92: 1018～1040.

[20] King G C P, Stein R S, Lin J. Static stress changes and the triggering of earthquakes[J]. Bull.seismol.soc.amer, 1994, 84: 935～953.

[21] 周翠英，刁桂苓，耿杰，等．1668年郯城大地震震源断层三维特征反演[J]．地球物理学进展，2013，28(6)：2814～2824.

ZHOU Cui-ying, DIAO Gui-ling, GENG jie, et al. 3-D characteristics inversion of hypocenter fault-plane of the 1668 Tancheng great earthquake[J]. Progress in Geophys, 2013, 28(6): 2814～2824.

THE EFFECT OF 1668 TANCHENG M8.5 EARTHQUAKE ON THE SEISMIC ACTIVITY OF THE VICINITY FROM COSEISMIC AND POSTSEISMIC DEFORMATION

XU Hao, SUN Yu-jun, WU Zhong-hai

(Institute of Geomechanics, Chinese Academy of Geological Sciences, Beijing 100081, China)

The Tanlu fault zone is an important active fault zone in east China. Many great earthquakes have occurred along it. The 1668M8.5 Tancheng earthquake is the largest recorded event. With the constraint of surface rupture, we get a coseismic rupture model based on the elastic dislocation theory. In addition, we calculate the coseismic, postseismic deformation and Coulomb stress change of this earthquake with a viscoelastic multilayered model. The results show that the rupture fault of the earthquake is adextral strike-slip fault, which has slight thrust. The coseismic displacement of this event is very large. The accumulated energy is released thoroughly. The coseismic rupture of the 1668 Tancheng earthquake increased the Coulomb stress on the northwest, northeast and southeast faults of Tanlu fault zone. The Coulomb stress increases further due to the postseismic deformation. After 350 years of the Tancheng earthquake, the Coulomb stress changes in Shanxian, Suqian and Rizhaoget to +1bar-+1MPa. When the mantle is imposed different viscosities, the duration time which the deformation and Coulomb changes become stable will be different. However, the stable values of them are generally same.

1668 TanchengM8.5 earthquake; coseismic and postseismic deformation; Coulomb stress change; viscosity

1006-6616(2016)03-0568-09

2016-04-24

中国地质调查局地质调查项目“长江经济带活动构造与区域地壳稳定性调查”(DD20160268)

徐昊(1993-)，男，硕士研究生，研究方向地球动力学。E-mail：xhdky2011@163.com

孙玉军(1983-)，男，博士，副研究员，主要从事地球动力学方面的研究。E-mail：sunyujunabc@163.com

P315.5

A