0.   引言

郯庐断裂带是中国东部一条重要的活动断裂带，历史上沿该断裂带发生过多次大地震。1668年7月25日晚，在山东南部发生了一次旷古未有的特大地震，根据历史资料得到该地震的震级约为M_s 8.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  研究区构造背景图(椭圆内加粗断裂为此次地震主震断层)

Figure  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×10¹⁸ Pa·s，最后会分析不同地幔粘滞性系数条件下，计算得到的形变量和库仑应力的差异，并结合李勇红等^[5]给出的郯城8.5级地震震源区附近地壳上地幔速度结构，给出最终的分层模型及参数(见表 1)。

表  1  分层模型及参数

Table  1.  Parameters of the multilayer crustal model

序号	深度范围/km	P波速度/(km·s-1)	S波速度/(km·s-1)	密度/(kg·m-3)	粘滞系数/(Pa·s)
1	0~10	5.00	2.89	2100	∞
2	10~15	6.16	3.56	2500	∞
3	15~20	6.20	3.58	2800	∞
4	20~30	6.35	3.67	2850	∞
5	30~35	6.32	3.66	2950	1.0×1018
6	35~40	6.50	3.76	3000	1.0×1018
7	＞40	8.18	4.73	3100	1.0×1018

下载: 导出CSV 

| 显示表格

2.2   库仑应力变化的计算

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

式中：Δσ_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]观测数据存在一定误差，由于在模型构建中缺少断层面上精准的滑移数据，因此认为断层模型依然能够在一定程度上反映当时的地震破裂变形。

图  2  同震及震后0 km处地表形变量

a—同震时刻地表水平方向形变量；b—同震时刻地表垂直方向形变量(向上为正)；c—350 a后地表水平方向形变量；d—350 a后地表垂直方向形变量(向上为正)

Figure  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间M_s1—M_s7地震事件发现，地震主要发生在这些库仑应力增加的区域内。虽然地震的发生与构造运动直接相关，但这样的特大型地震发生后对周围断层的应力积累影响十分明显，而且计算结果显示震后350 a该影响依然存在，只是目前已经趋于稳定。

图  3  同震库仑应力变化分布图(灰色圆圈为1668年郯城地震后至今研究区内发生的M1—M7地震事件)

a—地表处同震库仑应力变化分布；b—20 km深处同震库仑应力变化分布；c—地表处震后350 a库仑应力变化分布；d—20 km深处震后350 a库仑应力变化分布

Figure  3.  The distribution map of coseismic Coulomb stress change

下载: 全尺寸图片 幻灯片

同震库仑应力增加最明显区域主要位于断层破裂带的两端及同震破裂不完全的区域。断层破裂带应力释放明显，而破裂带的两端则为应力集中的区域。本次地震的震源深度为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)。库仑应力增量表明宿迁和日照等地与郯庐断裂带主震断层相同性质的断层地震危险性增强。

图  4  研究区内不同城市震后350 a内形变量和库仑应力变化

a—南京市；b—日照市；c—宿迁市；d—单县市△CFS—库仑应力变化量；U(x)—南北向位移(北方向为正)；U(y)—东西向位移(东方向为正)；U(z)—垂向位移(向上为正)

Figure  4.  Deformation and Coulomb stress change in different cities in study regions 350 years after the earthquake

下载: 全尺寸图片 幻灯片

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

图  5  不同地幔粘滞性系数对应的震后形变和库仑应力变化(南京市)

Figure  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处分布与库仑应力变化的关系时误差较大。