随着数字化地震台站布设密度的增加，地震波到时的测量精度得到提升。天然地震体波层析成像技术已成为研究壳幔结构的重要工具，但利用双差层析成像方法对沂沭断裂带及邻区进行成像分析的研究较少。本文在对沂沭断裂带及邻区进行地震重新精确定位的基础上，利用地震波走时数据对34°~37°N、117°~120°E范围进行精细双差层析成像分析，以期获得更可靠的地下介质速度结构信息，进而分析研究区速度分布及不同构造单元速度结构，从而揭示沂沭断裂带及邻区的地下介质性质。此外，为明确研究区地震分布与低速层的关系，从地壳三维速度分布出发，重点探讨沂沭断裂带等典型构造对地震活动的控制作用与孕震环境。

1 研究方法原理

双差层析成像(tomoDD)是一种体波走时层析成像技术，由Zhang等^[1-2]将双差定位方法^[3]引入地震层析成像理论而提出的。根据射线理论，如果2个相邻地震到同一台站的路径相似，那么震源i到地震台站k的观测到时可表示为：

$ T_{k}^{i}=\tau^{i}+\int_{i}^{k} \delta u \mathrm{~d} s $

(1)

式中，τⁱ为地震i的发震时刻，δu为慢度矢量，ds为径积分元。

若地震i被台站k所记录，则：

$ r_{k}^{i}=\sum _{l=1}^{3} \frac{\partial T_{k}^{i}}{\partial x_{l}^{i}} \Delta x_{l}^{i}+\Delta \tau^{i}+\int_{i}^{k} \delta u \mathrm{~d} s $

(2)

式中，r_kⁱ为理论到时与观测到时的时间差。震源参数(x₁ⁱ, x₂ⁱ, x₃ⁱ)、发震时刻、慢度场、射线路径均为未知参量。若地震j也被台站k所记录，则：

$ r_{k}^{j}=\sum _{l=1}^{3} \frac{\partial T_{k}^{j}}{\partial x_{l}^{j}} \Delta x_{l}^{j}+\Delta \tau^{j}+\int_{j}^{k} \delta u \mathrm{~d} s $

(3)

地震事件i、地震事件j在同一台站k与理论走时差的残差，即双差为：

$ \begin{gathered} r_{k}^{i}-r_{k}^{j}=\sum\nolimits_{l=1}^{3} \frac{\partial T_{k}^{i}}{\partial x_{l}^{i}} \Delta x_{l}^{i}+\Delta \tau^{i}+\int_{i}^{k} \delta u \mathrm{~d} s- \\ \sum\nolimits_{l=1}^{3} \frac{\partial T_{k}^{j}}{\partial x_{l}^{j}} \Delta x_{l}^{j}-\Delta \tau^{j}-\int_{j}^{k} \delta u \mathrm{~d} s \end{gathered} $

(4)

与传统地震层析成像方法相比，双差层析成像可改善地震波传播路径的差异和介质速度结构的空间不均匀性，能更好地约束地震的空间位置。双差数据的引入可有效提高震源区的成像分辨率，因此双差层析成像相比传统的走时成像具有明显的优势。在反演计算过程中，将相对走时残差和绝对走时残差的2阶范数作为约束函数，可迭代得到P波及S波三维速度结构和震源位置的最优解。

2 数据与方法 2.1 数据

本文搜集整理山东及周边省份地震台网提供的2008-10~2017-12震相观测资料和182个地震台站的资料信息，包含1 872个天然地震事件，最大震级为M4.2。图 1为研究区地震M-T图及地震深度分布情况，统计结果表明，M2.0以上地震1 371个，震源深度集中分布在5~25 km。图 2(a)为原始数据的走时-震中距关系曲线(简称时距曲线)，由于存在人工拾取等产生的误差，根据时距曲线剔除原始数据中离散程度较大的数据，结果如图 2(b)所示。经过筛选，满足条件的地震事件共1 190个，采用tomoDD程序中ph2dt工具将走时数据转化为绝对走时数据和走时差数据，最终获得P波绝对走时数据14 652条，S波绝对走时数据14 650条，相对走时数据89 722条。

2.2 初始模型选择

图 3为研究区二维射线路径分布情况，根据二维射线的密集程度，本文将经度、纬度、深度3个方向的网格个数分别设置为20、18、8。以0.2°×0.2°网格划分经纬度，以0.5°间隔对边缘进行划分，在深度方向上以5 km为间隔，分别在5 km、10 km、15 km、20 km、25 km、30 km、35 km深度开展研究。计算过程中，初始速度模型参考文献[4]，具体参数信息见表 1。

2.3 正则化参数选择

双差层析成像方法采用阻尼LSQR算法，为提高反演结果的稳定性，在反演过程中对引入的光滑因子和阻尼参数进行权衡分析^[5-6]，其中阻尼因子可在迭代过程中保证稳定迭代和快速收敛，而光滑因子则能解决迭代过程中模型的突变情况。本文分别选择10~800和0.1~2 000作为阻尼因子和光滑因子权重的测试范围，通过L曲线法选择最佳平滑因子20和阻尼因子200作为反演过程中的控制参数^[7]，结果如图 4所示。

3 反演结果可靠性评价

为更好地检验研究区速度结构反演空间特征的真实性，本文采用Lévěque等^[8]提出的棋盘检测板进行测试。检测板分辨率测试主要是通过正负方向的异常扰动来反演结果的真实性。将定义的速度模型节点设置为空间相间的棋盘分布，在计算理论走时过程中加入±5%的扰动速度来检测模型的稳定性，最后根据理论和实际反演的模型对正负速度异常的相对变化进行评估。在前期重定位研究中发现，研究区深度相对较浅，因此本文在检测板测试中仅给出20 km深度范围内的检测结果，以5 km作为1个距离区间，10 km层深的横向分辨率可达0.2°(图 5)。

4 地壳速度结构成像及分析 4.1 不同深度地壳介质水平向P波速度结构

对研究区内参与计算的地震事件进行重新定位，精定位前后的走时残差直方分布如图 6所示。图 6(a)和6(b)分别为初始定位和重新定位后的走时残差分布，对比发现，精定位后走时残差得到明显改善，定位精度显著提高。成像结果的可信度评估目前主要基于2个参数，分别为检测板测试恢复度RES(resolution, 最高恢复度为1)和反演网格节点处的射线空间分布DWS(distribution of derivative weight sum)。研究表明^[9]，当DWS≥100时，反演结果的可信度较高。经过综合考虑，本文选择DWS≥100且RES≥0.8的区域用于最终成像分析。

整体来看，研究区主要涉及3个大地构造单元(图 7(a))，主要包括鲁西地块、胶南地块和胶辽地块。鲁西地块和胶南地块分别处于沂沭断裂带西侧和东侧，其中胶辽地块为一个NNE向的狭长台隆陆块。图 7(b)~7(e)为不同深度的P波速度成像结果，其中胶南地块为苏鲁大别山造山带的东延区域，受构造作用的影响，该地区地壳变形强烈，上地壳内普遍存在低速结构；而鲁西块体速度分布极不均匀，这可能与其发育典型的幔源岩浆活动有关，尤其在泰莱凹陷区域出现不同规模的低速结构，而在泰安上地壳存在明显的高速结构。受限于反演分辨率，胶辽地块西南部分反演结果的可靠性较差，总体速度高于胶南地块，在此不作深入分析。对于沂沭断裂带而言，在不同深度均呈现南北两侧速度低、中间速度高的特征，体现出沂沭断裂带的分段特性。

图 7(b)~7(e)为横切剖面上下2 km地震的投影分布，从图中可以看出，高速-低速异常过渡带地区更容易发生地震。根据Artyushkov^[10]的观点，高速-低速过渡区多发中强地震的主要原因为流变边界处由于下地壳在地幔热流驱动作用下发生挤压，导致上地壳脆性岩石产生形变，当应力达到破裂阈值时，浅源地震随即发生。本文根据李清河等^[11]的相关研究，在图 8(e)中用红色五角星标识出1668年郯城8.5级大地震(34.8°N，118.6°E, 震源深度约为23 km)震中，发现该地震基本处于高速结构与低速结构之间的过渡带中。

4.2 垂向剖面成像结果

图 8分别为沿鲁西地块、沂沭断裂带、胶南地块和胶辽地块纵向和横向展布的4条P波速度剖面，分别表示为A-A′、B-B′、C-C′和D-D′。图 8中震源位置为离剖面垂直距离30~40 km范围内的地震投影。

剖面A-A′经过沂沭断裂带内部，从P波速度剖面可以看出，断裂带中上地壳速度结构分段特性明显，表现为南北段地壳速度低、中段速度高(郯城至莒县)。申金超等^[12]研究发现，郯城-莒县段内存在明显的地壳减薄带和高泊松比区域。剖面B-B′沿沂沭断裂带走向穿过鲁西地块，考虑到成像分辨率的可靠性，主要分析上地壳速度分布情况。从剖面可以看出，鲁西地块南部区域P波速度普遍高于北部。从地震活动性分析，多数地震发生在高速结构边缘及其向上延展的过渡带内。人工地震测深结果结合重磁异常表明，鲁西地块内具有深部岩浆向上运移的良好通道，考虑可能存在地幔热流上涌，使得地球深部应力传导至地壳，造成地壳脆性岩石破裂，继而引发地震。

剖面C-C′与D-D′近垂直于沂沭断裂带，就C-C′剖面而言，鲁西地块中上地壳P波速度普遍高于胶南地块，其中郯城大地震发生于高速结构边缘区域。郯城大地震震源区下方存在代表岩石圈拆沉的高波速异常，也存在代表热物质上涌的低波速异常，表明郯城大地震的发生与地幔热物质上涌有关，该结论与Lei等^[13]的最新研究结果一致，这种由“大地幔楔”中岩石圈拆沉引起热物质上涌至地壳浅部的结论也与Ma等^[14]关于高密集地震台阵的背景噪声成像研究及Tian等^[15]关于远震地壳结构研究的结果吻合。

D-D′剖面起点为泰山地区，该区域中上地壳表现为高速结构。孟亚锋^[16]研究认为，泰安区域地壳存在明显的隆生现象，其地壳高速结构对应的区域与隆升位置接近，主要原因为鲁西地块地幔热流的上涌直接影响至地壳。此外，泰山区域地震活动多沿高速-低速过渡带发生，原因可能为中地壳属于脆韧性转化的拆离面，该拆离面易于积累应力并引发地震。

5 结语

本文对沂沭断裂带及邻区的双差层析成像结果进行研究，主要得到以下结论：

1) 运用双差层析成像方法，基于中国地震台网中心提供的研究区内地震事件的观测报告和地震台站信息，获得研究区可靠的中上地壳三维速度结构，且大部分地区横向分辨率可达0.2°。

2) 沂沭断裂带及邻区速度结构整体具不均匀性，断裂带内速度结构具有明显的分段特性，郯城大地震发生的位置处于高速结构边缘区域，存在代表热物质上涌的低波速异常。郯城以南和莒县以北呈低速特征，而中间段具有高速结构，其中郯城和莒县中地壳存在明显的低速结构。沂沭断裂带两侧的胶南地块和鲁西地块速度结构存在差异，胶南地块中上地壳内P波速度普遍低于鲁西地块。鲁西地块内速度分布情况较为复杂，上地壳存在不同规模且连续性不强的低速结构，而泰山地区中上地壳表现为高速结构，其原因可能为鲁西地块地幔热流上涌直接影响至地壳，而泰山地区地壳高速结构对应的区域与隆升位置接近。

3) 本文揭示了三维P波速度结构与地震活动性的关系。多数地震均发生在高速结构边缘或高速-低速过渡带内，其中穿过1668年郯城8.5级地震震中附近的A-A′和C-C′剖面显示，震中基本处于高速结构边缘的过渡带上，发震原因主要为地幔热流上涌使地球深部应力传导至地壳，造成地壳脆性岩石破裂，继而引发地震。

致谢: 文章在撰写过程中应用到由中国科学技术大学张海江教授提供的tomoDD计算程序，在此表示感谢。