2020, Vol. 40
2. 黑龙江省地震局，哈尔滨市鸿翔路24号，150000
茅山断裂带位于苏南隆起的北部，属于江苏省西南部的下扬子断块。大尺度范围内受到的构造应力主要来自于太平洋板块向欧亚板块俯冲形成的挤压应力及从贝加尔经大华北直到琉球海沟大范围的引张应力[1]。目前普遍认为，茅山地区是由一系列外来岩片组成的推覆构造带，并伴有二次正断层作用，控制上白垩统和下第三系的沉积[2]。茅山断裂带及邻区(以下简称研究区)内构造复杂，发育数条切割深度达下地壳的大断裂，其中NNE向的茅山断裂带最大，对该地区构造活动有控制作用。如图 1(图中圆圈表示2008~2018年期间M≥1.0地震分布)所示，由NNE向的茅山断裂带(F1)、近EW向的幕府山-焦山断裂(F2)和NW向的南京-湖熟断裂(F3)相交形成倒三角形区域(以下简称“三角区”)，其中“三角区”的顶点区域是历史地震和现代破坏性地震的多发区[3]。
|
图 1 溧阳震区主要断裂分布 Fig. 1 Distribution of main fractures in Liyang area |
1974-04-22、1979-07-09江苏溧阳相继发生MS5.5和MS6.0地震，其中1979年的MS6.0地震是现代发生在江苏陆域内震级最高的地震，两次极震区基本重合，可视为同源地震。在5 a内同一地区重复发生两次破坏性地震在东部板内地震中非常少见[4]，具有十分重要的研究意义。就发震断裂而言，叶洪等[5]、Chung等[6]认为，MS5.5地震发震断裂应为NWW向断裂。刘万琴等[7]认为，MS6.0地震是NW向断层向SE方向发生的左旋走滑破裂。Chung等[6]利用长周期P波、SH波研究发现，NNE向的茅东断裂为发震断裂。侯康明等[3]研究认为，这两次地震的发震主断裂为金坛-南渡断裂。储飞等[4]通过计算同震库仑应力变化认为，两次地震的发震断层均为NW走向。
本文尝试采用双差层析成像技术，对研究区小尺度范围进行精细三维速度结构成像研究，更好地剖析地质体或者断裂构造的速度结构特征，探讨溧阳两次破坏性地震以及中小地震活动与速度分布的关系。
1 理论方法与数据 1.1 双差层析成像理论方法双差层析成像(tomoDD)的理论原理是：假设两个地震事件i、j在空间距离上相近，地震发生后由同一地震台站k记录到，而且它们到同一个台站的路径相似，两者可以组成一个地震对，则两次地震事件在同一台站记录的观测到时与理论到时之差的残差称之为双差，用
| $ \begin{array}{c} {\rm d}r_k^{ij} = \sum\limits_{m = 1}^3 {\frac{{\partial T_k^i}}{{\partial x_m^i}}} \cdot x_m^i + {\tau ^i} - \\ \sum\limits_{m = 1}^3 {\frac{{\partial T_k^j}}{{\partial x_m^j}}} \cdot x_m^j - {\tau ^j} \end{array} $ |
式中，Tki和Tkj为地震对的观测到时，xmi和xmj为地震对的震源参数，τi和τj为地震对的发震时刻。详细理论公式推导见文献[8-9]。在计算过程中，在地震对彼此相近的情况下，在震源区内利用到时数据进行精细速度结构的反演。在震源区外，其速度结构通过绝对走时数据反演进行计算。所以，在双差层析成像方法中，通过引入绝对走时克服在双差定位中台站到地震事件之间速度恒定的假设，避免了常规地震定位中的弥散现象[10]，减小了系统偏差，因此地震定位结果更准确，同时能够获得几百米尺度的高分辨率成像结果。
双差层析成像通过节点法对模型参数化，采用球状介质中的伪弯曲法进行射线追踪[8, 11]，反演过程中采用阻尼最小二乘算法，在3个方向采取相同的光滑权重对模型进行光滑约束，以相对走时残差和绝对走时残差的二阶范数作为目标函数，多次迭代直至得到稳定的解。
1.2 数据本文收集整理由中国地震台网中心提供的2008~2018年30°~34°N、117°~121°E范围内的近震初至Pg、Sg波震相走时数据，剔除爆破、塌陷等非天然地震事件，得到共计天然地震事件516个。地震活动性如图 2所示，灰色柱状图表示地震发生的月频次，蓝色实心点表示每个地震的震源深度，重定位前地震震源深度主要集中分布在5~20 km。
|
图 2 研究区地震随时间分布和M1.0以上地震深度分布 Fig. 2 Temporal distribution of earthquakes and depths of M≥1.0 earthquakes |
采用hypoDD方法[10]提取相对走时与绝对走时数据，具体选择参数如下：事件对到台站最大距离maxdist=400 km；地震对之间的最大距离maxsep=10 km，最小距离minsep=0.1 km；每个事件的最大邻居数maxngh=10；每个地震事件对所需要的震相数的最小值minobs=8。最大值maxobs=50；一个事件如果能成为另一个事件的邻居，那么该事件对需要的相同震相的个数的最小值minlink=8。图 3(a)给出了原始数据中走时-震中距关系图(简称时距曲线)，剔除初始数据中离散度较大(图 3(a)中绿色曲线以外)的走时信息后，满足条件的地震事件共计460个，参与计算的台站82个，时距曲线如图 3(b)所示。最终提取P波绝对走时数据6 267条，S波绝对走时数据6 267条，相对走时数据共28 306条。
|
图 3 P波和S波时距曲线 Fig. 3 Time-distance curves of P and S wave |
根据地震射线的密集程度选择网格大小，如图 4所示。分别沿经度、纬度和深度3个方向设置反演网格节点，研究区域内以0.2°×0.2°网格进行划分，边缘以0.5°进行划分，即沿经度方向设置为117.50°、118.00°、118.20°、118.40°、118.60°、118.80°、119.00°、119.20°、119.40°、119.60°、119.80°、120.00°、120.20°、120.40°、120.60°、120.80°、121.00°；沿纬度方向设置为30.50°、31.00°、31.20°、31.40°、31.60°、31.80°、32.00°、32.20°、32.40°、32.60°、32.80°、33.00°、33.50°；沿深度方向设置为1.00 km、5.00 km、10.00 km、15.00 km、20.00 km、25.00 km、30.00 km、35.00 km。
|
图 4 研究区二维射线路径分布 Fig. 4 Distribution of P-wave ray paths in the studied area |
李细兵等[12]基于多震相联合反演得到江苏地区的一维速度模型，结果表明，中上地壳20 km以内的速度变化较小，整个地壳的平均速度为6.0 km/s，上地幔顶部速度为8.06 km/s。如表 1所示，本文采用该模型作为初始P波速度模型，初始波速比设置为：σ=VP/VS=1.731，得到初始S波速度，然后将一维速度模型转换为三维网格点速度点阵进行后续反演计算。
|
表 1 双差层析成像采用的初始速度模型 Tab. 1 Initial velocity model used in tomoDD |
双差层析成像反演计算过程中，共进行4次迭代计算，每次速度结构与地震定位联合反演之后增加1次定位反演，以减小联合反演时速度收敛快于定位计算的影响。双差层析成像采用带阻尼的最小二乘QR分解算法，该算法中阻尼系数damp和平滑系数smooth的选择对反演结果的稳定性有较大影响。其中，平滑权重因子在反演过程中用来约束慢度的变化量，而阻尼因子同时约束震中位置和慢度的变化量。因此在反演之前，必须对两者进行权衡分析[13]。本文通过L-curve方法[14]计算最优的参数值。如图 5所示，根据归一化模型与走时残差关系曲线，选择最佳平滑因子20和阻尼因子120作为反演过程中的控制参数。
|
图 5 平滑因子和阻尼因子均衡曲线 Fig. 5 Trade-off curves of damping and smoothing weight parameters |
地震分布、台站分布、地震射线密集程度和射线交叉等情况均会影响速度结构反演的分辨率[15]。为了检验在实际数据和网格模型下反演得到的速度结构的可靠性及空间分辨能力，需要预先进行检测板测试[16]。主要步骤如下：1)在反演网格点上正负交错添加5%的速度扰动；2)以此作为理论模型通过正演计算理论走时；3)根据该理论走时和实际反演的初始速度模型对未添加扰动的网格点速度进行反演，根据扰动值的恢复情况对成像结果进行评价。检测板测试控制参数与真实数据时参数一致，平滑因子和阻尼因子采用上文的20和120，结果如图 6所示。
|
图 6 不同深度剖面检测板测试结果 Fig. 6 Checkboard resolution test results of each depth |
由图 6可见，研究区内不同深度地层切片上均能恢复正负相间的速度异常变化，在5 km、10 km、15 km、20 km多数区域的P波速度横向分辨率可达0.2°。随着深度增加，分辨率有所降低，主要原因是天然地震和地震台站分布不均匀，导致射线覆盖性有所差异。在实际计算中需要选择恢复度高的区域，以确保成像结果的可靠性。
3 成像结果与讨论 3.1 不同深度地壳介质水平向P波速度结构双差层析成像过程中对参与计算的地震事件进行双差精定位。图 7给出重定位前后地震分布与走时残差分布的变化。由图 7看出，重定位后的地震发生深度分布在5~20 km；重新定位前走时残差主要分布在0.2~0.5 s，平均走时残差为0.47 s，重新定位后绝对大多数地震走时残差分布在-0.2~0.2 s，平均走时残差为0.14 s，地震定位得到明显的改善。结合上节检测板测试结果，本文主要分析检测板测试恢复度RES(resolution)高于0.8(最高恢复度为1)的上地壳三维速度结构成像结果。另外需要说明的是，速度成像的精度一般利用反演网格节点处的射线空间分布进行评估。Thurber等[17]定义参量DWS(distribution of derivative weight sum)表征节点周围平均的相对射线密度，当DWS≥100时，反演结果精度可信度较高。本文通过RES和DWS共同约束成像结果，以保证最终结果的可靠性。
|
图 7 重定位前后地震分布及走时残差变化 Fig. 7 Epicenter distribution and mean square of travel-time residual before and after relocation |
图 8为不同深度下的P波速度成像结果，图中的震中分布表示的是该切片层上下2 km精定位地震的投影情况。根据ISC国际地震中心给出的结论[5]，1974年溧阳MS5.5地震发震深度为16 km，将其投影至图 8(c)；1979年溧阳MS6.0地震发震深度为12 km，将其投影到图 8(b)。从发震位置的速度结构分析，MS5.5地震明显处于高速-低速体过渡带中，偏向低速体的边缘，而MS6.0发震位置的速度结构变化性不大，相对也处于高速-低速体过渡带中，接近高速体边缘。另外，整体来看，多数中小地震也发生在低速体内或高速-低速过渡带之内。不同层深的速度结构显示研究区存在明显的横向不均匀性，以下从不同层深进行速度结构的分析。5 km深度初始速度为5.88 km/s，反演之后的速度变化率为-11.0%~11.9%(图 8(a))，低速区处于茅山断裂带中段和“三角区”以内。NNE向的茅山断裂带与近EW向的幕府山-焦山断裂交汇处以及与NW向的南京-湖熟断裂交汇处均显示高速结构，并且在陈家堡-小海断裂附近也出现低速体。综合来看，沿着茅山断裂带走向自北向南基本呈现出低速区和高速区相间分布的情况。10 km层深的初始速度为5.88 km/s，反演后的速度变化率为-10%~12.9%(图 8(b))。幕府山-焦山断裂北部地区和陈家堡-小海断裂处于低速区；南京-湖熟断裂和溧阳-南渡断裂两侧速度存在差异，南侧高于北侧，“三角区”顶点区域处在高低速过渡区。上地壳基底15 km深度的初始速度为6.05 km/s，反演后的速度变化率为-8.4%~11.7%(图 8(c))。茅山断裂带与溧阳-南渡断裂交汇处出现明显的低速体结构，而与南京-湖熟断裂交汇处出现高速体，两者由茅山断裂带分开。其他区域基本与5 km和10 km层深的速度分布具有一定的一致性。20 km层的初始速度为6.25 km/s，反演得到的速度变化率为-8.9%~11.5% (图 8(d))。本层由于射线密度和恢复度的限制，研究区内的西南区域速度结构不可靠，在此不作分析。反演结果表明，该层速度结构较上地壳而言，陈家堡-小海断裂邻区由低速体转换为高速体。
|
图 8 研究区不同深度P波速度成像结果(RES≥0.8且DWS≥100) Fig. 8 Map views of the P-wave velocity perturbations of different depth(RES≥0.8 and DWS≥100) |
通过反演计算得到不同深度下网格点的速度值，最终利用三维数据线性插值算法得到研究区整体的P波速度结构，如图 9所示。研究区三维P波速度结构反映出上地壳基底与中地壳具有比较明显的速度界面，上地壳基底速度约为6.1 km/s，而中地壳速度约为6.5 km/s，并且界面起伏较大。如图 9所示，AA′剖面基本沿茅山断裂带走向穿过“三角区”，BB′剖面近平行于AA′ 剖面穿过溧阳震源区，CC′剖面近垂直于AA′剖面穿过震源区。图 10给出3条剖面的垂向速度成像结果，其中地震表示为剖面两侧20 km范围内的地震投影。
|
图 9 研究区三维立体P波速度结构 Fig. 9 3D P wave velocity structures |
|
图 10 深度剖面的速度结构 Fig. 10 Vertical cross-sections of velocity structure profile |
溧阳震源区大致划分为3层：盖层构造层、基底构造层及下地壳[18]，如图 10中虚线所示。其中，基底构造层(8~20 km)包括浅变质岩系和深变质岩系，两者分界深度约为12 km，在深变质岩系内18 km处存在低速体层界面；上地壳与下地壳分界深度约20 km。根据垂向剖面结果与震源区的地壳划分，可以发现一些基本特征：1)从图 10中发现，整体震源区的速度分布与胡连英等[18]划分的地壳有一定的相似性，但是也存在一些差异，就本文的成像结果而言，深变质岩系的深度应该在22~25 km左右；2)从图 10(b)发现，震源区内存在拆离的高速体；3)震源区基底构造层是由浅变质岩系和深变质岩系组成，前者属于低速、低阻、易破碎的上-中远古界浅变质岩系，而后者是高速、高阻、不易破碎的下远古-上太古深变质岩系，反演结果显示，不同岩性的分界面不平整，起伏较大。
图 10(a)AA′剖面显示了沿着茅山断裂带走向的“三角区”内地壳速度变化，可以看出，茅山断裂带东北段速度低于西南段。在江苏省上沛镇下方15 km处出现较大规模的高速体，而溧阳两次破坏性地震震源深度处于高速体的边缘部分。从图 10(b)BB′剖面可以更明显地看出两次地震的异同，MS5.5地震发生在深变质岩系内低速层上，而MS6.0地震发生在深变质岩系与浅变质岩系边界上，这是由于地壳内高速体和低速体之间介质参数差异大，更容易发生较大地震。图 10(b)显示，震源区下方的高速体上升至盖层基底附近；图 10(c)也显示，在震源深度附近出现低速体向上扩展的形态。依据速度结构出现的特征分析，溧阳两次地震的孕震环境主要是受上地幔物质向上热涌受阻，逐步积累弹性应变能形成震源体，在脆弱-韧性转换低强度部位出现震源破裂点。另外，据地表考察、精密磁测、显微和超显微构造等研究表明[18]，茅山断裂带内存在玄武岩喷发带和茅东盆地中轴含幔源包体玄武岩喷发带，历经多次构造作用，含幔源包体玄武岩已形成脆-韧性的复合部位，从速度结构成像的角度也说明了这一现象的存在。
4 结语1) 采用中国地震台网中心提供的516个近震事件，结合周边82个测震台站信息，通过双差层析成像方法获得研究区地壳的三维P波速度结构和地震重新定位结果，地震定位前后走时残差由0.47 s降低至0.14 s，明显改善了地震定位精度。反演计算之前，通过L曲线法获得反演过程中最优正则化参数：平滑因子20和阻尼因子120；另外，结合检测板测试，选择恢复度RES≥0.8和射线密度DWS≥100指标约束三维速度结构反演结果，最终发现，大部分地区的横向分辨率可以达到0.2°。
2) 反演后的速度与初始速度模型之间的变化率基本在±10%左右。从不同深度P波速度成像结果来看，P波的速度结构横向差异变化比较明显。5 km深度的低速区处于茅山断裂带中段西侧的“三角区”内，其两端均显示明显高速体结构；10 km深度下，“三角区”顶点区域处在高低速过渡区；在15 km深度，溧阳震源区在茅山断裂带东西两侧分别出现明显的低速体和高速体。P波速度剖面结构与前人划分的溧阳震源区地壳3层构造具有相似性，但是各层界面不平整，并且在震源深度附近出现低速体向上扩展的形态。此外，深变质岩系的深度(也就是中地壳深度)应该在22~25 km范围。
3) 从地震投影的情况分析，多数地震发生在低速体和高速体-低速体的过渡带内。本文成像结果表明，1974年溧阳MS5.5地震发生在深变质岩系内的低速区，而1979年溧阳MS6.0地震发生在深变质岩系与浅变质岩系边界上，两次地震均存在一定的孕震环境，本文从速度结构上说明了发震深度的可靠性。而发震原因可能在于上地幔物质向上热涌受阻，逐步积累弹性应变能形成震源体。16 km深度处于低速区内，岩石具有塑性流变，具备发生较大地震的孕震环境，而溧阳地区12 km深度处于速度梯度较大的转换复合位置，介质参数差异大，故有利于大震的发生。
致谢: 中国科学技术大学张海江教授提供tomoDD计算程序，中国地震台网中心提供地震观测报告，在此表示衷心感谢。
| [1] |
徐纪人, 赵志新. 苏鲁-大别造山带岩石圈应力场、构造运动特征以及超高压变质带折返机制的研究[J]. 岩石学报, 2007, 23(12): 3317-3324 (Xu Jiren, Zhao Zhixin. Study on Characteristics of Lithospheric Stress Field, Tectonic Motions and Its Exhumation Mechanism of the Ultrahigh-Pressure Metamorphic Belt in and around the Sulu-Dabie Orogenic Belt[J]. Acta Geologica Sinica, 2007, 23(12): 3317-3324 DOI:10.3969/j.issn.1000-0569.2007.12.026)
(  0) |
| [2] |
滕龙, 邸兵叶, 张俊. 江苏茅山断裂带中段重力异常多尺度分析及其构造识别[J]. 资源调查与环境, 2015, 36(1): 49-56 (Teng Long, Di Bingye, Zhang Jun. Multi-Scale Analysis of Gravity Anomalies and Its Discernment of Structures in the Middle Segment of the Maoshan Faulted Zone[J]. Resources Survey and Environment, 2015, 36(1): 49-56 DOI:10.3969/j.issn.1671-4814.2015.01.007)
( 0) |
| [3] |
侯康明, 熊振, 李丽梅. 对江苏省溧阳2次破坏性地震发震构造的新认识[J]. 地震地质, 2012, 34(2): 303-312 (Hou Kangming, Xiong Zhen, Li Limei. New Insights into the Seismogenic Structures of the Two Destructive Earthquakes in Liyang, Jiangsu Province[J]. Seismology and Geology, 2012, 34(2): 303-312 DOI:10.3969/j.issn.0253-4967.2012.02.009)
( 0) |
| [4] |
储飞, 王韶稳, 张毅, 等. 1974年和1979年溧阳两次地震同震库仑应力变化及其发震构造研究[J]. 地震研究, 2016, 39(3): 397-403 (Chu Fei, Wang Shaowen, Zhang Yi, et al. Research on Coseismic Coulomb Stress Changes Caused by Two Liyang Earthquakes in 1974 and 1979 and Their Seismogenic Tectonic[J]. Journal of Seismological Research, 2016, 39(3): 397-403 DOI:10.3969/j.issn.1000-0666.2016.03.006)
( 0) |
| [5] |
叶洪, 张文郁, 于之水, 等. 1979年溧阳6级地震震源构造的研究[J]. 地震地质, 1980, 2(4): 27-38 (Ye Hong, Zhang Wenyu, Yu Zhishui, et al. On the Source Tectonics of 1979 Liyang Earthquake of Magnitude 6[J]. Seismology and Geology, 1980, 2(4): 27-38)
( 0) |
| [6] |
Chung W Y, Wei B Z, Brantley B J. Faulting Mechanisms of the Liyang, China, Earthquakes of 1974 and 1979 from Regional and Teleseismic Waveforms——Evidence of Tectonic Inversion under a Fault-Bounded Basin[J]. Bulletin of the Seismological Society of America, 1995, 85(2): 560-570 DOI:10.1016/0009-2541(94)00129-V
( 0) |
| [7] |
刘万琴, 黄家正. 用瑞雷波广义方向性函数研究5~6级地震的破裂过程[J]. 地震学报, 1982(1): 27-34 (Liu Wanqin, Huang Jiazheng. A Study of Fracture Process of Earthquakes with Magnitudes MS=5-6 by the Generalized Directivity Function of Rayleigh Wave[J]. Acta Seismologica Sinica, 1982(1): 27-34)
( 0) |
| [8] |
Zhang H J, Thurber C H. Double-Difference Tomography: The Method and Its Application to the Hayward Fault, California[J]. Bulletin of the Seismological Society of America, 2003, 93(5): 1875-1889 DOI:10.1785/0120020190
( 0) |
| [9] |
Zhang H J, Thurber C. Development and Applications of Double-Difference Seismic Tomography[J]. Pure and Applied Geophysics, 2006, 163(2-3): 373-403 DOI:10.1007/s00024-005-0021-y
( 0) |
| [10] |
Waldhauser F, Ellsworth W L. A Double-Difference Earthquake Location Algorithm: Method and Application to the Northern Hayward Fault, California[J]. Bulletin of the Seismological Society of America, 2000, 90(6): 1353-1368 DOI:10.1785/0120000006
( 0) |
| [11] |
郭雨帆, 罗佳宏, 宫伟, 等. 日本-千岛海沟转折处地震层析成像研究[J]. 大地测量与地球动力学, 2018, 38(6): 639-645 (Guo Yufan, Luo Jiahong, Gong Wei, et al. P-Wave Seismic Tomography of the Japan-Kuril Transition Zone[J]. Journal of Geodesy and Geodynamics, 2018, 38(6): 639-645)
( 0) |
| [12] |
李细兵, 赵启光, 朱峰, 等. 基于多震相联合反演江苏地区一维速度模型研究[J]. 中国地震, 2018, 34(2): 312-321 (Li Xibing, Zhao Qiguang, Zhu Feng, et al. The Study of 1D Velocity Model Based on Multi-Phase Joint Inversion in Jiangsu Province[J]. Earthquake Research in China, 2018, 34(2): 312-321)
( 0) |
| [13] |
Ma X, Westman E C, Fahrman B P, et al. Imaging of Temporal Stress Redistribution Due to Triggered Seismicity at a Deep Nickel Mine[J]. Geomechanics for Energy and the Environment, 2016(5): 55-64 DOI:10.1016/j.gete.2016.01.001
( 0) |
| [14] |
Hansen P C, O'Leary D P. The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems[J]. SIAM Journal on Scientific Computing, 1993, 14(6): 1487-1503 DOI:10.1137/0914086
( 0) |
| [15] |
刘建明, 李志海, 冯雪玲, 等. 2012年新源和静6.6级地震前后P波速度结构的演化[J]. 大地测量与地球动力学, 2017, 37(5): 492-496 (Liu Jianming, Li Zhihai, Feng Xueling, et al. Evolution of P Wave Velocity Structure before and after the Xinyuan-Hejing MS6.6 Earthquake in 2012[J]. Journal of Geodesy and Geodynamics, 2017, 37(5): 492-496)
( 0) |
| [16] |
Lévěque J, Rivera L, Wittlinger G. On the Use of the Checker-Board Test to Assess the Resolution of Tomographic Inversions[J]. Geophysical Journal International, 1993, 115(1): 313-318 DOI:10.1111/j.1365-246X.1993.tb05605.x
( 0) |
| [17] |
Thurber C, Eberhart-Phillips D. Local Earthquake Tomography with Flexible Gridding[J]. Computers and Geosciences, 1999, 25(7): 809-818 DOI:10.1016/S0098-3004(99)00007-2
( 0) |
| [18] |
胡连英. 江苏溧阳地震震源构造[J]. 江苏地质, 1997(3): 43-50 (Hu Lianying. Hypocentral Structure of Liyang Earthquakes, Jiangsu Provience[J]. Jiangsu Geology, 1997(3): 43-50)
( 0) |
2. Heilongjiang Earthquake Agency, 24 Hongxiang Road, Harbin 150000, China