重力异常反映的郯庐带中南段构造边界与应力场

朱益民 余 腾 王 鑫 张 鹏

1 中国地质大学(北京)土地科学技术学院，北京市学院路29号，100083 2 宿迁学院建筑工程学院，江苏省宿迁市黄河南路399号，223800 3 宿迁学院地质环境与工程健康监测重点实验室, 江苏省宿迁市黄河南路399号，223800 4 应急管理部国家自然灾害防治研究院，北京市安宁庄路1号，100085 5 南京工业大学交通运输工程学院，南京市浦珠南路30号，211816

通过分析重力异常的梯度变化能够确定异常体的边缘与规模，进而研究地壳内密度异常变化及其与地震活动性的联系[1]。部分学者通过垂线偏差反演构造应力场，探讨应力状态与地质体运动的基本特征[2-3]。国内外学者利用重力梯度数据识别地质体边界、划分构造单元[4-5]，借助重力总水平梯度和垂线偏差研究区域构造应力场[6-7]，基于重力场模型计算不同深度下的地壳密度异常与构造应力场[8]。

针对郯庐带的构造特征研究多基于GPS位移分析[9]、地震反射剖面[10]、地应力测量[11]、地壳速度结构及泊松比[12]、地震重定位与震源机制解以及地震活动性[13-15]等数据，针对郯庐带及邻区的重力异常数据处理多采用小波多尺度分解方法[16-17]以及Tilt导数法[18]，鲜有基于重力全张量梯度数据进行的郯庐带中南段构造特征研究。基于重力场一阶导数定义的传统方法仅适用于浅层构造和矿产勘探目标界定[19]，且会因为同时存在正负异常而产生错误边界；而重力全张量梯度数据已在南海断裂圈定[20]、四川盆地地质构造水平位置界定[21]等研究与应用中展现出较好的适用性。因此，基于重力全张量研究郯庐带中南段地质体交接关系与构造单元划分具有一定的实用意义。

本文通过重力全张量数据重新定义水平导数THDR和解析信号振幅ASM，并提出改进的倾斜导数NTDM。经过模型验证，该方法可均衡识别不同埋深下的地质体边界，具有较好的稳定性与收敛性。同时，利用垂线偏差数据反演研究区构造应力场与壳幔密度差发现，二者在空间上具有较强的相关性。研究区最大主应力方向总体呈EW向，与前人研究结果一致。

1 研究区概况

1.1 郯庐带中南段地形及构造特征

郯庐带是我国东部规模最大的NNE走向活动断裂，其中新沂以北至渤海为郯庐带中段，新沂以南为郯庐带南段。郯庐带中南段贯穿本文研究区(32.7°～35.2°N、117.0°～119.5°E)，研究区内以浅源地震为主，震群较为集中。区域内多条斜列的NWW、EW向断裂与郯庐带中南段交接，组成复杂的断陷带(图1)。第四纪以来，在西太平洋板块弧后扩张与印度板块向北碰撞的共同作用下，郯庐带表现为逆冲兼右旋走滑的活动性质。

图1 研究区地形与主要构造特征[22]Fig.1 Topography and main structural features of the study area[22]

1.2 郯庐带中南段布格重力场、莫霍面埋深及地壳厚度特征

本文收集由BGI官方网站发布的2′×2′格网大小的WGM2012模型计算布格重力异常(图2(a))，此外还收集了Sandwell V31的1′×1′垂线偏差数据，用于反演研究区构造应力场、壳幔密度差和重力总水平梯度。地壳厚度与莫霍面埋深数据来自Crust1.0(图2(b))。

郯庐带是一条重要的地球物理分隔带。由图2(a)可见，NNE走向的重力异常分隔带呈现西低东高的特征，反映出华北板块和扬子板块基底性质的不同，深部构造复杂；西侧的华北板块表现为数值偏低的布格异常区，变化趋势为EW向。F10(淮阴-响水断裂)和F11(洪泽-沟墩断裂)共同形成弧形串珠状低重力异常区，分隔下扬子板块。

郯庐带也是一条NNE走向的莫霍面突变带。由图2(b)可见，总体上东侧莫霍面埋深较浅，西侧莫霍面埋深较深。莫霍面埋深由北至南逐渐增大，局部还存在较为明显的分段性。

2 数据与方法

2.1 重力场方向导数法与全张量梯度

对布格重力异常网格数据进行x、y、z三个坐标轴方向的求导，得到重力矢量g在笛卡尔坐标系中沿坐标轴方向的分量gx、gy、gz，可反映浅层地质体的交接关系，初步揭示区域构造格局。重力全张量梯度为重力向量沿坐标轴方向分量的一阶导数，其在三维直角坐标系中可表示为：

(1)

总水平导数THDR表达式为[1]：

(2)

解析信号振幅ASM表达式为[23]：

(3)

图2 研究区布格重力异常和地壳厚度与莫霍面埋深Fig.2 Bouguer gravity anomaly and crustal thickness and Moho face buried depth of the study area

由于式(2)、(3)是基于一阶导数定义的边界识别方法，在圈定较深地质构造的边界时存在较大误差，因此为使浅层与深部地质体边界圈定达到均衡状态，本文利用重力全张量数据重新定义式(2)、(3)：

(4)

(5)

对THDRM的总水平导数与ASMM的垂向导数进行融合处理，定义基于重力全张量的倾斜角法NTDM：

(6)

THDRM和ASMM的局部极值和最值处可以界定地质体边界，倾斜导数NTDM在±π/2的振幅范围内随倾斜角发生显著变化。零交叉点接近异常体结构边缘，局部极值能有效圈定其边界，同时也能判定地质体在平面位置的延伸方向。经过3种不同实验环境下的模型验证后认为，该方法在均衡探测深浅断裂的同时，对边界识别也具有较好的稳定性与收敛性，能清晰区分异常体区域与边界，便于边界解译。图3为式(2)～(5)的数值模拟对比图，图中黑色虚线为地质体水平投影所对应的边界。

2.2 构造应力场与壳幔密度差

构造应力是在正压力(静岩应力)上附加的一种应力，文献[2]在Turcotte公式的基础上，基于垂线偏差推导其与构造应力场的转换公式。基于此，本文推导出如下公式：

(7)

(8)

(9)

由垂线偏差定义的最大主应力方向为：

(10)

图3 数值模拟对比Fig.3 Comparison of numerical simulation

3 地质体交接关系识别与分析

3.1 重力异常方向导数

对布格重力数据进行方向导数处理，可突出其局部梯度变化信息，位场解析延拓也可使目标场源特征更加明显。

图4(a)～(c)为经方向导数处理后的重力异常与研究区断裂位置的对应情况。由图4(a)可见，重力场x方向导数gx表现出近NS向构造边界，横向上具有显著差异，西侧为正梯度带，东侧为负梯度带，变化幅度为0.4×10-3mGal/km。由图4(b)可见，重力场y方向导数gy表现出近EW向地质体边界，纵向上的正负异常梯度交替出现。以F6(无锡-宿迁断裂)为界，研究区中部以低异常梯度为主，分隔南北正异常梯度带。由图4(c)可见，重力场z方向导数gz表现出垂向变化特征，与区域布格重力异常的相似性较高，西北侧的华北板块转变为平缓的低异常梯度，无明显变化。浅层地壳的刚体性质较为明显，郯庐带中南段处在北宽南窄的正异常梯度带上。

(a)、(b)中的断裂分别为筛选后的NS向、EW向断裂图4 郯庐带中南段布格重力场方向导数与向上延拓情况Fig.4 Directional dervatives and upward continuation of Bouguer gravity anomaly in the middle-southern segment of Tan-Lu fault zone

3.2 基于全张量梯度的NTDM边界识别方法

3.2.1 模型验证

为验证NTDM方法的可行性，建立4个不同体积与底部埋深的棱柱体模型(图5)，对THDRM、ASMM、NTDM以及重力全张量gxx、gxy、gxz、gyy、gyz、gzz进行模拟计算与验证。 模型基本参数如表1所示。

表1 棱柱体模型参数

图6为模型的重力梯度异常(图中红、橙、黑、绿实线位置表示模型M1、M2、M3、M4在x轴上的投影，图7、8同理)，其中gxx、gyy、gzz可识别特定方向的边界。埋深较深的模型M2、M4的异常振幅与周围存在明显差异，能界定异常体的规模。gxx对边界位置识别较为发散，gxy可反映异常体角点信息,gxz、gyz可对应边界信息。与THDRM相比，ASMM对不同埋深的异常振幅大小的均衡探测效果较好，二者均对埋深较深的模型异常体规模具有较好的揭示效果。对于埋深最浅的模型M1，THDRM未能对其进行有效界定。剖面图能更直观地展示异常振幅的差异和不同方法的分辨率特征。如图6(g)中黑色虚线和红色箭头所示，THDRM的局部极小值可界定异常体边界；而ASMM主要以最大值圈定边界，均衡圈定的能力较好，且ASMM边界识别的收敛性强于THDRM。由图6(i)可见，NTDM具有边界增强作用，边界识别的分辨率高于图6(a)～(h)，NTDM剖面图的极大值可对应异常体边界。

图5 模型剖面图、三维图、重力异常Fig.5 Section, 3D graph, gravity anomaly of model

图6 模型重力异常的全张量及经全张量组合后的边界识别结果Fig.6 Full tensors of gravity anomaly and the boundary identification results after the combination of the full tensors

本文在重力梯度分量中加入5%的高斯随机噪声，以验证NTDM的抗噪能力(图7)。从图7(g)～(i)可以看出，NTDM具有较强的抑噪能力，剖面图中虚线位置与模型的水平投影位置对应较好，NTDM在均衡探测异常体边界与抗噪能力上均最优。但由于NTDM的表达式中包含三阶导数，会放大异常体周边的噪声影响，因此剖面线的异常扰动较为明显。

通过向上延拓的方法对数据进行滤波处理，延拓后的数据可抑制较浅的异常体。由图8(g)～(i)所示，模型M1底部埋深最浅，THDRM、ASMM、NTDM对其边界信息的识别总体较弱，其中NTDM的识别效果优于THDRM、ASMM，这是因为THDRM和ASMM的识别易受到周边异常体的干扰，结果较为发散。3个边界识别方法对埋深较大异常体的规模圈定效果较好，NTDM在地质体周边的正负异常对比使其具有更稳定的边界识别能力。

图7 加入5%高斯噪声重力异常的全张量及经全张量组合后的边界识别结果Fig.7 Full tensors of gravity anomalies with 5% Gaussian noise and the boundary identification results after the combination of the full tensors

图8 加入5%高斯噪声并向上延拓1 km后重力异常的全张量及经全张量组合后的边界识别结果Fig.8 Full tensors of gravity anomaly with 5% Gaussian noise and 1 km upward continuation, and the boundary identification results after the combination of the full tensors

经3种实验环境验证，THDRM、ASMM、NTDM的识别效果总体上优于单个重力全张量梯度。THDRM主要以局部极小值确定边界，ASMM、NTDM主要以局部极大值和最大值圈定边界。3种方法均具有较好的抗噪性与收敛性，不会产生错误边界，其中NTDM的边界识别结果更清晰。

3.2.2 实际数据应用

为验证本文方法在实际重力场数据中的应用效果，利用上述方法对郯庐带中南段及邻区布格重力异常数据进行处理，确定地质体的交接关系与构造边界(图9)。由图可见，重力梯度gxx中央的低异常区与郯庐带中南段对应较好，但异常区规模更大、东西向更宽；gyy、gyz能识别EW向断裂，同时斜切的EW和NW向断裂也使得郯庐带中南段具有分段性；受到F6(无锡-宿迁断裂)、F7(铜石-甘霖断裂)的阻隔作用，gzz可反映郯庐带中南段深部构造的不连通性。

相较于ASMM，THDRM边界识别的分辨率更高，南部的串珠状弧形分隔带有利于揭示边界信息，但ASMM在圈定地质体规模上具有一定的优势。新方法NTDM在均衡深部和浅部构造特征上更具优势，反映出的信息与已有成果对应较好，还可以展示周边更丰富的细节信息，不受噪声干扰，可为地质构造单元的划分提供依据。

图9(i)边界识别结果显示，以郯庐带为界，研究区中部的地质体东西两侧受到郯庐带的阻挡作用，延伸角度的弧度多为负值，方向呈NW、EW向，与郯庐带交接；研究区南部地质体延伸方向发生明显扭转，在水平向与垂直向上均发生变化。F11(洪泽-沟墩断裂)的东南侧存在一个由EW向转为NE向的异常地质体边界，可作为构造单元划分的依据。

图9 郯庐带中南段重力异常边界识别结果Fig.9 Identification results of gravity anomaly boundary in the middle-southern Tan-Lu fault zone

综合图9(g)～(i)的识别结果，结合前人研究资料[24-26]，对郯庐带中南段及其周边地区进行划分，共得到14个构造单元(图10)。由图可见，区域Ⅰ郯庐带作为二级构造单元，分隔东西不同的地块，区域Ⅱ可能为胶辽地块南沿部分。F8(海泗断裂)具有北延迹象，其北界与F9(邵店-桑墟断裂)和郯庐带交接位置可作为区域Ⅲ苏北胶南断块的划分依据。 F10(淮阴-响水口断裂)不仅是华北陆块区和秦祁昆造山系的分界线，同时也是鲁西陆块和大别-苏鲁造山带的分界线。根据边界识别结果可知，北延的次级断裂可划分区域Ⅳ连云港-泗洪断块。区域Ⅴ盐阜坳陷、区域Ⅵ建湖隆起和区域Ⅶ东台坳陷可由研究区东南方位的弧形异常边界进行划分，且异常边界与F11(洪泽-沟墩断裂)的对应较好。研究区西侧呈现出弧形的NW向异常带(方位角弧度值为负)，分隔了区域Ⅷ鲁中隆起和区域Ⅸ鲁西断块，以及区域Ⅹ邳县-睢宁断块和区域Ⅺ徐淮坳褶带。区域ⅩⅢ灵璧-泗县凹陷属于区域Ⅻ淮北断陷带，可由NW向陡变为EW向的异常边界进行划分。ⅩⅢ灵璧-泗县凹陷和区域ⅩⅣ蚌埠隆起可由EW向的弧形异常边界进行划分。

图10 郯庐带中南段构造单元划分Fig.10 Division of tectonic units in the middle-southern segment of Tan-Lu fault zone

4 基于垂线偏差的构造应力场反演

根据§2.2理论基础，借助重力垂线偏差数据反演研究区构造应力场、壳幔密度差和重力总水平梯度(图11)。

蓝色实线为沿郯庐带剖面线作垂直于剖面线左右各30 km的剖面(图中等间隔的白色实线)，绘制成上图中蓝色曲线图11 构造应力场、壳幔密度差、重力总水平梯度Fig.11 Tectonic stress field, crust-mantle density difference, total horizontal gradient of gravity

由图11可见，构造应力场、壳幔密度差和重力总水平梯度东西向差异明显，水平方向上很不均匀，由图中曲线可以看出，郯庐带是一条构造应力差异变化分隔带。在构造应力变化明显的区域，壳幔密度和重力总水平梯度变化也较大。断裂带两侧构造应力等值线密集，多呈规则的环形，地震大多发生在2个环形等值线的交接处。同时，环形构造应力等值线或许可以作为局部构造单元的划分依据，如多个环形等值线与区域ⅩⅢ灵璧-泗县凹陷、区域ⅩⅣ蚌埠隆起的对应较好，研究区东南侧的弧形变化区与F11(洪泽-沟墩断裂)的延伸方向也具有相似性。

图12为郯庐带中南段构造应力场矢量图，由图可见，最大主应力方向整体以EW向为主。中段东北角主应力方向由NW转为SW，该区域构造应力值较大、等值线密集，受到明显的挤压作用，最大主应力方向为NEE-EW；南段以F10(淮阴-响水口断裂)为分界线，东南角由NW向一直延伸至蚌埠一带，淮安-洪泽一带由NW和SW向转为EW向并延伸至郯庐带及以西地区，最大主应力变化区也与高构造应力值对应较好。郯庐带中南段应力场整体方向呈EW向，与受到太平洋板块向欧亚大陆俯冲和印度洋板块共同作用形成的中国大陆东部的应力场一致。表2为前人对最大主应力方向的研究成果。

实线为郯庐带中段与南段分界线图12 郯庐带中南段最大主应力方向Fig.12 The direction of the maximum stress in middle-southern segment of Tan-Lu fault zone

表2 最大主应力方向对比

5 结 语

本文根据重力全张量定义的总水平导数与解析信号振幅，提出新的倾斜导数边界识别方法NTDM，并将其应用到郯庐带中南段布格重力异常数据的处理与解释中。借助垂线偏差可反演得到郯庐带中南段的构造应力场特征，得出以下结论：

1)当正负异常同时出现并含有噪声时，NTDM能表现出良好的抑噪性与收敛性，能清晰识别异常体边界、均衡不同深部的异常体振幅。相比于单一的全张量数据，THDRM、ASMM在界定局部边界时具有一定的优势，NTDM能更准确地圈定边界。

2)依据边界识别结果并结合相关资料，以区域Ⅰ郯庐带中南段为界，大致可将该区域划分为区域Ⅱ胶辽地体、区域Ⅲ苏北胶南断块、区域Ⅳ连云港-泗洪断块、区域Ⅴ盐阜坳陷、区域Ⅵ建湖隆起、区域Ⅶ东台坳陷、区域Ⅷ鲁中隆起、区域Ⅸ鲁西断块、区域Ⅹ邳县-睢宁断块、区域Ⅺ徐淮坳褶带、区域Ⅻ淮北断陷带、区域ⅩⅢ灵璧-泗县凹陷、区域ⅩⅣ蚌埠隆起等14个构造单元。

3)构造应力场、壳幔密度差与重力总水平梯度均存在明显的以郯庐带中南段为界的水平差异，三者与区域重力异常的相关性较大，构造应力值范围为-28～-12 Mpa，东西向最大差异约为16 Mpa，壳幔密度差水平向最大差异约为0.2 g/cm3。构造应力等值线密集区的构造运动强烈，地震频发。

4)研究区最大主应力方向总体呈EW向，其中郯庐带中段为NEE-EW，郯庐带南段为EW-NWW，与中国大陆东部应力场保持较高的一致性。

致谢：应急管理部国家自然灾害防治研究院和宿迁市地震服务中心提供相关资料；开源Python库GEOIST(https:∥cea2020.gitee.io/geoistdoc/)对本文模型验证提供支持；文中部分图件由GMT绘制，并参考有关实例，在此一并表示感谢。