沂沭断裂带是郯庐断裂带在山东境内的统称，将整个山东地区分割成鲁西隆起和胶辽隆起两大块体(图 1)，是我国著名的地震活动带(方仲景等，1986；李家灵等，1994；王小凤等，2000；张鹏等，2015)，也是岩石圈尺度的构造边界带(高维明等，1980；刘昌铨等，1983；马杏垣等，1991；卢造勋和夏怀宽，1992；Yin and Nie, 1993；陈沪生等，1993；万天丰等，1996；朱光等，2002；牛漫兰等，2005；Chen et al., 2006；Li et al., 2014；严乐佳等，2014；刘保金等，2015；Huang et al., 2015)，呈“两堑夹一垒”的构造特点(陈国星等，1988；高维明等，1988；晁洪太等，1998；李洪奎等, 2009, 2017；张鹏等，2010).沂沭断裂带位于鲁西隆起、胶辽隆起与苏鲁褶皱带之间(图 1)，呈NNE-SSW走向，北起莱州湾，向南至新沂一带，主要由昌邑-大店、沂水-汤头、白芬子-浮来山、鄌郚-葛沟以及安丘-莒县五条大型深断裂组成，是郯庐断裂带地震活动最为强烈的部分，1668年郯城8½级强震就发生在这条断裂上.沂沭断裂带作为活动地块的边界带，第四纪以来以右旋走滑兼具逆冲为主，其两侧活动块体间的相互作用及其差异运动对地震活动具有直接控制作用(张培震等，2003；张国民等，2004；王志才等，2005).近几年以来，沂沭断裂带及相邻地区地震活动较为活跃(图 2)，相继在沂沭断裂带北部及东侧地区发生显著的持续时间较长的长岛震群、莱州震群以及乳山震群活动.特别是乳山震群的地震活动从2013年10月1日开始一直持续至今，已发生地震近14000次，其中M_L＜0.0地震3980次、0.0≤M_L≤0.9地震7983次、1.0≤M_L≤1.9地震1662次、2.0≤M_L≤2.9地震262次、3.0≤M_L≤3.9地震28次、4.0≤M_L≤4.9地震3次和M_L≥5.0地震1次，该震群活动频度与持续时间在我国大陆东部地区较为罕见.2011年3月11日日本9.0级巨震发生以来，沂沭断裂带区域被中国地震局确定为重点防御区，因此，弄清该地区深部结构，建立尽可能详细的上地幔顶部速度结构，对于进一步了解该地区地震孕育发生的动力学背景具有重要科学意义.

图 1

Fig. 1

图 1 研究区主要区域地质构造示意图 黑色实线为主要活动断裂.Taihang orogen，太行山造山带；North China basin，华北盆地；Luxi uplift，鲁西隆起；Jiaoliao uplift，胶辽隆起；NYSB，北黄海盆地；Sulu belt，苏鲁褶皱带；NSYSB，南黄海北部盆地；SSYSB，南黄海南部盆地；Yishu Fault，沂沭断裂带；Tanlu Fault，郯庐断裂带.虚线代表块体边界, 圆圈代表 7级以上地震, 红色五角星代表1668年郯城8½级地震. Fig. 1 Sketch map of regional geological tectonics in the study region Black solid lines mark major active faults. The dashed lines represent block boundaries, whereas the circles denote epicenters of earthquakes with magnitudes greater than 7.0. The red star is the 1668 Tancheng M8½ earthquake.

图 2

Fig. 2

图 2 2008-2019年期间地震震中分布(中国地震台网中心编目系统http://10.5.160.18/console/index.action) 蓝色方框为乳山地区. Fig. 2 The epicenter distribution of the earthquakes that occurred during 2008-2019 (China Earthquake Networks Center Cataloging System http://10.5.160.18/console/index.action) The blue box represents Rushan area.

上地幔顶部是壳幔物质交换与能量传递的关键区域(李志伟等，2011)，上地幔顶部Pn波速度变化与上地幔压力、温度、物质构成有关(Hearn，1996；McNamara et al., 1997；Liang et al., 2004；He et al., 2017)，而各向异性反映了地幔物质的运动及形变特征(Hearn，1996；Liang et al., 2004；Bamford，2010).因此，开展沂沭断裂带及周边区域上地幔顶部Pn波速度和各向异性结构研究，对于进一步理解该地区地震孕育的深部地球物理环境、发震机制和地球动力学过程具有重要科学意义和现实意义.

近年来，已有大量学者利用Pn波层析成像方法获得了多个中国大陆地区上地幔速度结构模型(如汪素云等, 2001, 2003；裴顺平等，2004；Hearn et al., 2004；Liang et al., 2004；Pei et al., 2007；Zhou and Lei, 2016).围绕沂沭断裂带也开展了一些相关研究，比如吕子强等(2016)和顾勤平等(2016, 2020)分别反演了环渤海地区及郯庐断裂带中南段Pn波速度结构及各向异性，殷伟伟等(2019)利用2008-2016年2级以上地震的Pn波到时资料，获得郯庐断裂带Pn波成像结果.这些研究尽管获得了有意义的研究结果，但由于资料的局限性其反演结果的分辨率相对低些.在本研究中，我们采用手工拾取研究区域2008-2019年期间2级以上地震波形的高质量Pn波到时资料，并结合2008-2019年期间观测报告中的Pn波到时资料.另外，前人工作中仅使用了山东省内的地震台站接入中国地震台网的40个台站到时资料，而本研究还增加了未接入中国地震台网的88个台站数据资料.这88个台站的加入，使得台站密集地分布在沂沭断裂带两侧，大大提高了研究区射线交叉性覆盖，对于揭示沂沐断裂带细深部特征具有重要意义.本研究所用的26598条Pn波到时数量远多于前人研究中所用到的资料数量，使得整个沂沭断裂带及周边地区Pn射线交叉性覆盖更为密集，为获取更高分辨率结构特征奠定了基础.

为使沂沐断裂带能被更多的地震射线有效覆盖，本研究将资料搜集范围扩大至112°E-124°E、32°N-42°N，采用手工识别方式提取到2008-2019年期间2级以上地震事件的33397条Pn射线，并结合地震观测报告2008-2019年期间2级以上地震事件的37253条Pn射线.在合并两个数据集时，如遇重复射线时以手工拾取的到时为准.为保证数据的质量，我们采用以下标准来选取资料：(1)震中距为1.8°~12°，确保首先到达的是Pn震相；(2)震源深度小于33 km；(3)每个地震事件至少有6个台站记录到Pn震相，每个台站至少要记录到6个事件的Pn震相；(4)走时残差在±3.0 s内.按上述标准作为条件，我们最终挑选出研究区内290个台站记录到1665个地震事件的26598条Pn到时资料.图 3为Pn波射线路径分布图，可以看出整个研究区内的射线交叉覆盖较为密集.

图 3

Fig. 3

图 3 研究区内290个台站(绿色钻石号)记录到1665个地震(红色圆圈)的26598条Pn波射线(蓝色细线)分布 Fig. 3 Distribution of 26598 Pn ray-paths (blue thin lines) from 1665 earthquakes (red circles) recorded at 290 seismic stations (green diamonds)

本研究采用Hearn(1996)的成像方法，反演沂沭断裂带及周边区域上地幔顶部Pn波速度和方位各向异性结构.该方法将上地幔顶部划分成二维网格，每个网格的速度和方位各向异性是待反演项.Pn波走时残差t_ij可表示为

(1)

式中a_i和b_j分别为台站i和地震事件j的静延迟，d_ijk为事件j到台站i经过网格k的射线长度，s_k是第k个网格的慢度扰动，A_k、B_k是第k个网格的各向异性系数，φ为台站相对于事件的方位角(即反方位角).各向异性大小可表示为，各向异性快波方向为1/2arctan(B_k/A_k)+90°(Pei et al., 2007；Lei et al., 2014).反演网格采用0.2°×0.2°，根据CRUST1.0模型(Laske et al., 2013)，研究区上地幔顶部Pn波平均速度为7.99 km·s^-1，壳内P波平均速度为6.17 km·s^-1，莫霍面平均埋深为33 km，使用阻尼最小二乘迭代算法来求解式(1).在求解过程中，选取一组合适的速度阻尼因子与各向异性阻尼因子，以平衡高分辨率和低误差之间的折中问题，两者之间的比值控制着速度和各向异性量值的相对大小.这两个阻尼因子确定过程如下：首先不包含各向异性，将不同速度阻尼因子进行多次反演得到速度均方根与走时残差均方根并绘制L曲线，在该曲线的最大曲率处70被确定为最佳速度阻尼因子(图 4a)；其次，计算速度和各向异性的相对折中以确定各向异性阻尼因子.其中，速度对各向异性的折中定义为反演获得的各向异性扰动的均方根与速度扰动的均方根的比值(RMS_Ani/RMS_Vel)，而各向异性对速度的折中定义为反演获得的速度扰动的均方根与各向异性扰动的均方根的比值(RMS_Vel/RMS_Ani).保持速度阻尼因子70不变，改变各向异性阻尼因子进行反演，当速度对各向异性的折中与各向异性对速度的折中相等时的各向异性阻尼因子110被确定为最佳各向异性阻尼因子(图 4b)，因此，本研究将速度阻尼因子70和各向异性阻尼因子110应用于反演过程中.通过反演，得到了沂沭断裂带及周边区域上地幔顶部Pn波速度结构与各向异性分布图像.图 5为反演前后走时残差分布图，可以看到反演后的走时残差相对于反演前得到了明显收敛.

图 4

Fig. 4

图 4 速度阻尼因子(a)与各向异性阻尼因子(b)的选取.有关细节请参阅正文 Fig. 4 Damping parameters chosen for velocity (a) and anisotropy (b). For details, see the text

图 5

Fig. 5

图 5 资料走时残差分布图 (a) 反演前；(b) 反演后. Fig. 5 Distribution of travel time residuals (a) Before the inversion; (b) After the inversion.

为了检验反演结果的可靠程度，使用检测板实验方法对使用数据的分辨能力进行检验.初始模型Pn波速度设为7.99 km·s^-1，在此基础上在每个网格单元上交替添加一定的速度大小扰动和各向异性方向扰动，采用与实际资料相同的资料数与反演方法进行反演，计算不同异常尺度的分辨率模型.图 6和图 7分别展示出研究区0.5°×0.5°的速度异常恢复结果和1.0°×1.0°各向异性异常恢复结果，可以看出绝大部分区域异常的形态和幅度均得到比较好的恢复，只有东部海域边缘地区恢复的异常存在变形，这与射线交叉性覆盖程度(图 3)相一致，表明研究区域速度异常分辨率达到0.5°×0.5°，而各向异性异常分辨率达到1.0°×1.0°，因此本文仅讨论与此分辨率相当的异常尺度的结构特征.

图 6

Fig. 6

图 6 异常尺度为0.5°×0.5°的Pn波速度检测板分辨率实验结果 Fig. 6 Results of the checkerboard resolution test for Pn velocity with an anomaly size of 0.5°×0.5°

图 7

Fig. 7

图 7 异常尺度为1.0°×1.0°的Pn波各向异性检测板分辨率实验结果 Fig. 7 Results of the checkerboard resolution test for Pn anisotropy with an anomaly size of 1.0°×1.0°

图 8给出了研究区域Pn波速度结构(模型结果数据参见电子附件)，其中蓝色和红色分别表示高速异常和低速异常.由图可以看出，研究区域速度结构存在明显的横向不均匀性，与区域地质构造明显相关，且异常形态与地质构造比较吻合.华北盆地显示出明显的高速异常，异常形态清晰地勾勒出华北盆地构造块体的轮廓，且盆地内部速度结构也存在显著的横向不均匀性.相对于前人的研究结果(汪素云等，2003；Huang and Zhao, 2004；李志伟等，2006；Lei et al., 2008；于湘伟等，2010；Li et al., 2011；吕子强等，2016；殷伟伟等，2019)，本研究可以更清晰地展现出华北盆地一些细结构特征，比如华北盆地内部强高速异常表现出明显的分块现象，这在以往前人研究结果中从未揭示出来.华北盆地南端表现为低速异常，苏鲁褶皱带呈现低速异常但不连续，显示出两个低速异常块体并沿断裂带的走向相间排列，其两侧的速度异常也明显不同，分别表现出高速异常和低速异常，表明断裂带两侧上地幔顶部的介质强度不同，具有不同的构造环境特征.太行山造山带呈现低速异常，且异常形态沿山脉走向分布.沂沭断裂带东西两侧的速度异常差别明显，并表现出明显的分段特征，大致以鲁西隆起和郯城地震为界线，鲁西隆起以北、沂沭断裂带西侧地区表现为高速异常，其东侧表现为低速异常；郯城地震以南地区则相反，其西侧表现为低速异常，东侧则表现为高速异常.这些结构特征表明，沂沭断裂带可能深延达上地幔顶部(叶高峰等，2009；徐义刚等，2009；牛漫兰等，2000；漆家福等，2010；黄耘等，2011；Lei，2012).朱日祥等(2012)研究指出太平洋板块的俯冲是华北克拉通破坏的一级动力学因素，郯庐断裂带在克拉通破坏之前就存在，其岩石圈强度较低，有利于岩浆运移，在克拉通破坏期间，成为强烈的伸展带和岩浆带.最近的华北地区上地幔顶部研究结果(Lü et al., 2020)与本研究呈现很相似的反演结果，且也支持该推论.鲁西隆起表现为低速异常，可能源于太平洋板块西向深俯冲形成“大地幔楔”结构(Lei and Zhao, 2005, 2006; Lei et al., 2020)导致深部地幔热物质沿沂沭断裂带向上涌动形成，其伸展构造形成过程中幔源岩浆活动及热侵蚀起到重要作用(燕守勋等，1996)，这种热物质上涌还得到最近地壳结构成像结果(孟亚锋等，2019；Ma et al., 2020; Tian et al., 2020)的支持.沂沭断裂带东侧的胶辽隆起表现为低速异常，与鲁西隆起的异常程度表现不同，暗示这两个隆起区上地幔顶部具有不同的介质强度.黄海海域盆地几乎全部表现为高速异常，说明盆地内部上地幔介质刚度较高，结构稳定.胶东半岛与南黄海北部为低速异常且相连，说明两区域的地幔介质性质可能相近.

图 8

Fig. 8

图 8 获得的研究区Pn波速度和各向异性结果 红色和蓝色分别表示低速和高速异常，黑色短线方向为快波方向，长度与各向异性的大小成正比，速度与各向异性色标位于图顶.四个白色椭圆分别代表了四个特定区域.其他标示与图 1中的相同. Fig. 8 Pn velocity and anisotropy inferred in the study region The red and blue colors denote low- and high-velocity anomalies. The black short lines show the fast direction of Pn propagation, the length of which is proportional to the amount of anisotropy. The scales of velocity and anisotropy are shown on the top. Four white ellipses denote four specific areas. The other labeling is the same as that in Fig. 1.

图 8展示出研究区7级以上的地壳内强震震中分布，红色五角星代表1668年郯城8½级地震位置.由图可以看出，强震发生的区域与Pn波速度分布密切相关，多数强震发生在高速和低速异常交界区域上方的地壳内，特别是速度变化剧烈的地带.雷建设等(2009)对龙门山断裂带地壳精细结构与汶川地震发震机理研究指出，汶川地震余震主要发生在速度表现出比较强烈的横向不均匀的地区，且汶川主震震源区下方存在有显著的低波速异常体.也有学者采用地震层析成像研究揭示华北大震大多都发生在高速与低速区的过渡带上(王椿镛等，2017).魏斯禹和滕吉文(1985)对我国几个大地震地区地壳与上地幔顶部结构关系的分析指出，上地幔顶部的局部隆起地带是地壳大地震孕育和发生的深部构造背景，造成局部隆起的物理机制是由于上地幔结构和物质组成的不均匀性以及上地幔热物质的上涌，而上地幔顶部Pn波速度的变化与上地幔压力、温度、物质构成密切相关，上述这些研究结果也都印证了本文的观点，说明壳内强震的孕育发生与上地幔顶部速度结构分布有关.低速异常说明介质温度较高(汪素云等，2003)，容易产生韧性形变，导致应变能积累；高低速异常交界区域意味着相应区域介质性质存在明显差别，应力在流变边界集中(Artyushkov，1973；Zhang et al., 2009)，容易诱发地震.地震活动十分活跃的乳山和长岛地区就位于低速异常区域的边缘上(图 8)，暗示乳山震群和长岛震群的发生可能与Pn波速度结构有关.前人研究结果表明，通常发生地震的断裂带两侧的速度异常存在较为明显差别(雷建设等，2009)，沂沭断裂带和苏鲁褶皱带交汇处速度结构横向差异明显，容易导致应力集中，产生应力差，1668年郯城8½级地震震中的位置就位于此处，这可能是郯城地震发生的主要成因之一(图 8)，这一认识在其他地区的研究结果中也得到证实(雷建设等，2009；黎源和雷建设，2012；Lei et al., 2014；Zhou and Lei, 2015).

Pn波各向异性通常被认为是由于上地幔顶部物质发生变形导致橄榄岩晶格优势排列引起的(Hess，1964；Christensen，1984；Nicolas and Christensen, 1987)，反映最近一次构造变形的特征.总体上，在简单剪切形变情况下，各向异性的快波方向表现为与剪切带走向平行，如板块边缘和穿透地壳的剪切带以及岩石圈地幔的剪切拖曳带；在纯剪切形变情况下，快波方向与最大拉伸方向一致，如岩石圈中没有断裂活动的地区(Savage and Silver, 1993；Silver，1996；Hearn, 1996, 1999).图 8为研究区Pn波各向异性(模型结果数据参见电子附件)分布图，反映了研究区上地幔顶部的形变特征，其中短线方向表示快波传播方向，其长短表示各向异性的强弱.可以看出，Pn波各向异性与区域地质构造和速度结构具有一定的关联.尽管个别稳定地块因为保留了较早形成的“化石”各向异性，其各向异性呈现出一定程度的强各向异性，如黄海盆地(图 8中的白色大椭圆1)，但多数呈现出低速异常的太行山造山带、华北盆地南部以及胶辽隆起区各向异性相对较强(图 8)，这些区域上地幔介质韧性大，容易发生形变，这与青藏高原东缘地区上地幔顶部定量统计结果(Lei et al., 2014)相一致.在沂沭断裂带内部Pn波各向异性快波优势方向与该断裂带走向基本一致，说明沂沭断裂带的活动以简单剪切为主，剪切深度达上地幔顶部.华北盆地内部各向异性变化较为明显，西部与太行山造山带交界区域Pn波各向异性快波方向为近E-W向，盆地中部弱低速区各向异性快波方向为近N-S向(图 8中的白色大椭圆2).太行山造山带地区(图 8中的白色大椭圆3)和胶辽隆起区域(图 8中的白色大椭圆4)的Pn波各向异性快波优势方向为近NNE-SSW向，与断裂带的走向趋于一致，与该区域地质构造的伸展方向趋于一致，而SKS快波方向为NW-SE方向(常利军等, 2009, 2012；苗庆杰等，2011；吴萍萍等，2012)，这种显著差异暗示该区域岩石圈在一定深度可能存在局部解耦现象.由于郯庐断裂大规模剪切运动使东侧南黄海扬子块体向北运动嵌入中朝块体，在南黄海地区形成若干NE-SW或NNE-SSW向伸展的隆起和坳陷，这与利用重力异常资料推测得到的黄海地区断裂体系是由N-S向和NE-SW或NNE-SSW向两组断裂构成的研究结果(郝天珧等，2002；王谦身等，2003)相一致.在南黄海北部盆地Pn波各向异性快波方向为近NE-SW向，而南黄海南部盆地Pn波各向异性快波优势方向为近N-S方向(白色椭圆1)，与区域断裂走向基本平行，与构造伸展方向趋于一致，暗示南黄海地区部分断裂的深度已达上地幔顶部成为岩石圈尺度的断裂.

为了解手工提取资料与观测报告中资料对成像结果的差别，本研究将手工提取的33397条Pn到时资料(图 9a)和观测报告中的37253条Pn射线(图 10a)分别开展反演研究，图 9b和图 10b展示出这两种资料集获得的速度与各向异性结果.由图中可以看出，两种资料集成像的结果较为接近，均显示出沂沭断裂带呈现明显的分段特征，其两侧的速度异常差别比较明显，华北盆地展现明显的高速异常，且盆地内部速度结构存在明显的横向不均匀性，太行山造山带、鲁西隆起以及胶辽隆起区域表现为低速异常，苏鲁褶皱带呈现不连续的低速异常，黄海海域盆地几乎全部表现为高速异常.太行山造山带地区和胶辽隆起区域的Pn波各向异性快波优势方向为近NNE-SSW向，南黄海北部盆地Pn波各向异性快波方向为近NE-SW向，南黄海南部盆地Pn波各向异性快波方向为近N-S向.然而，这两种资料的反演结果也存在着一些细节差别.譬如，在渤海海域，利用观测报告资料获得的结果显示出以西地区与断裂带内部为高波速异常、而以东地区为弱低波速异常(图 10)；利用人工拾取资料获得的反演结果显示出郯庐断裂带以西为明显高波速异常、而以东地区为明显低波速异常的分界特征(图 9).这些结果说明，相对于观测报告(图 10a)，利用人工拾取的Pn到时资料具有较好一致性、射线有较好的交叉性覆盖(图 9a)，可获得更为有意义的研究结果.

图 9

Fig. 9

图 9 手工提取Pn波射线分布(a)及相应的Pn波速度与各向异性成像结果(b) 其他标示分别与图 3和图 8中的相同. Fig. 9 Distributions of manually picked Pn ray-paths (a) and their results of velocity and anisotropy (b) The other labeling is the same as that in Figs. 3 and 8, respectively.

图 10

Fig. 10

图 10 观测报告中的Pn波射线分布(a)及相应的速度与各向异性成像结果(b) 其他标示分别与图 3和图 8中的相同. Fig. 10 Distributions of Pn ray-paths from observational bullletins (a) and their results of velocity and anisotropy (b) The other labeling is the same as that in Figs. 3 and 8, respectively.

通过手工拾取大量高质量的Pn波走时数据，利用Pn波层析成像方法，重建了沂沭断裂带及周边区域上地幔顶部高分辨率Pn波速度与各向异性结构模型，该模型明显地展示出与地表地质构造相关性，且揭示出前人所未揭示的结构特征，为深入理解沂沐断裂带深部结构与动力学过程提供了重要地震学依据.

(1) Pn波速度与各向异性结构和区域地质构造具有一定的关联.华北盆地显示明显的高速异常，且盆地内部强高速异常表现出明显的分块现象.苏鲁褶皱带展现两个低速异常块体并沿断裂带的走向相间排列，其两侧的速度异常也不相同，表明断裂带两侧具有不同的构造环境特征.沂沭断裂带东西两侧速度结构差别明显，且表现出明显的分段特征，大致以鲁西隆起和郯城地震为界线.构造活动强烈的太行山造山带、鲁西隆起区、胶辽隆起区的速度结构呈现低速异常，表明沂沭断裂带下方存在地幔热物质上涌.黄海海域盆地几乎全部表现为高速异常，说明其构造活动不受华北地区地幔热扰动的影响.

沂沭断裂带Pn波各向异性快波方向与断裂走向基本一致，说明沂沭断裂带的活动以简单剪切为主，且剪切深度达上地幔顶部.华北盆地内部各向异性差异较大.太行山造山带、胶辽隆起区Pn波各向异性快波方向与地质构造的伸展方向基本一致，与区域SKS快波方向存在明显差异，说明研究区内岩石圈在一定深度可能存在局部解耦现象.南黄海盆地Pn波各向异性快波方向与区域断裂的走向基本一致，暗示扬子块体向北挤压在南黄海地区形成的部分断裂的深度已达上地幔，成为岩石圈尺度的断裂带.

(2) 强震的空间分布与上地幔顶部速度分布密切相关.地壳强震主要发生在上地幔顶部的高速和低速异常的交界区域，表明壳内强震的孕育发生受到上地幔顶部构造环境影响，这可能与地幔热物质上涌、上地幔顶部介质横向不均匀性相关.郯城8½级地震位于沂沭断裂带和苏鲁褶皱带的交汇处，该区域速度结构横向差异大，是构造应力易于集中增强的地区，说明郯城地震的发生明显地受到了深部构造控制作用.