重力数据在笛卡尔坐标系x、y、z方向上的分量为：

式中，G为万有引力常数；r为测量点（x，y，z）到地质体（ξ，η，ζ）的距离；ρ为地质体密度。传统的边界识别方法是利用了重力数据3个方向导数的组合，如总水平导数（total horizontal derivative，THD）、解析信号振幅（analytic signal amplitude，ASM）、斜导数和θ导数。总水平导数为水平导数的平方和（${\mathit{T}}_{\mathrm{T}\mathrm{H}\mathrm{D}}$），为：

解析信号振幅为x、y、z方向导数的平方和（${\mathit{T}}_{\mathrm{A}\mathrm{S}\mathrm{M}}$），为：

斜导数（${\mathit{T}}_{\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}}$）为：

θ导数（${\mathit{T}}_{\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{s}}$）为：

式中，g_x、g_y、g_z分别为重力异常g的水平x方向、y方向及垂向z方向的异常导数。

从式（4）和式（5）可以看出，斜导数、θ导数方法是以重力g的x、y、z方向而建立的角度边界识别。但是这种单方向的角度边界识别法都存在着同样的问题，即在较深的地质体边界识别中误差大，且当地质体相邻很近时，其分辨能力不足。为了更加准确地探测目标体，可利用全张量梯度的边界识别方法。重力全张量梯度数据为重力向量g=（g_x，g_y，g_z）在3个正交方向（x，y，z）空间的导数，重力梯度全张量成分如图 1所示。重力全张量梯度可以表示为：

图  1  重力梯度全张量成分

Figure 1.  Components of Full Tensor Gravity Gradient

式中，${\mathit{g}}_{\alpha \beta }=\frac{\partial {\mathit{g}}_{\alpha }}{\partial \beta }$，$\alpha \mathrm{、}\beta \in \left\{x,y,z\right\}$。其中，

在x、y、z方向中，共有9个重力梯度张量。g_xx、g_yy、g_zz满足拉普拉斯方程：${\mathit{g}}_{xx}+{\mathit{g}}_{yy}+{\mathit{g}}_{zz}=\mathrm{0}$，且有${\mathit{g}}_{xy}={\mathit{g}}_{yx}$，${\mathit{g}}_{xz}={\mathit{g}}_{zx}$，${\mathit{g}}_{yz}={\mathit{g}}_{zy}$。

类似于式（4）和式（5），令水平面内所有水平分量g_yy、g_xy、g_yx、g_xx为水平导数，g_xz、g_yz为垂直导数，可建立公式如下：

本文所用的地形数据来自Sandwell和Smith海深数据，如图 2所示，数据分辨率为1$\mathrm{\text{'}}$×1$\mathrm{\text{'}}$。空间重力异常资料同样来自于Sandwell和Smith全球卫星重力异常，数据分辨率为1$\mathrm{\text{'}}$×1$\mathrm{\text{'}}$，在海域数据精度可达到3.03×10^-5 m/s²。

图  2  南海海深数据图

Figure 2.  Bathymetric Data of the South China Sea

本文对南海自由空气重力异常采用洋陆分离的方法进行布格校正。在海洋区域，采用传统Parker-Oldenburg正演将地层分为海水层、沉积层以及基底3层来进行计算；在陆地，则是采用传统地面地形校正^[25]来进行计算。地层密度取2 670 kg/m³，海水密度取1 030 kg/m³。将陆地和海域两个结果进行叠加得出了南海及其周边地区的布格重力异常，如图 3所示。

图  3  南海布格重力异常图

Figure 3.  Bouguer Gravity Anomalies of the South China Sea

由图 3可见，由洋盆区域向边界海、大陆架以及陆地区域方向，布格重力异常值在慢慢地降低。一些大的盆地构造，如莺歌海盆地、中建南盆地、南沙海槽盆地、曾母盆地等，布格重力异常表现出高值，这直观地反映了南海地区布格重力异常与基底起伏呈负相关。从布格重力异常走向上来看，南海北部陆架由于受到华南大陆和海区构造的影响，布格异常散乱，整体呈NE向分布。北部陆坡向中央海盆异常值增加，存在布格重力异常急剧增大的梯度带。南部陆坡布格重力异常存在不连续的重力梯度带，走向以NE、NW向为主。西部陆坡布格异常变化比较平缓，总体呈SN走向，部分异常呈NW走向。东部陆坡布格异常变化梯度大。西南海盆区内布格重力异常走向明显，呈NE向展布。东部海盆异常区异常走向不明显，大部分区域为正异常。

基于前人研究成果，比较多种常用的重力梯度边界识别方法，包括直接利用重力梯度三分量法和全张量梯度组合法，分析它们的优缺点，并将其用于南海异常地质体边界和断裂位置的圈定。图 4给出了重力梯度全张量异常。

图  4  全张量重力梯度

Figure 4.  Full Tensor Gravity Gradient Data

由图 4可见，g_xx可以有效地反映地质构造南北向的特征，g_yy可以更好地反映地质构造东西向的特征，g_xy可以更多地反映地质构造4个方位角的特征，g_zx可以识别地质体南北边界，g_zy可以识别地质体的东西边界。由于直接由全张量梯度分量识别的地质体边界较模糊，分辨能力较差，因此研究了不同梯度分量组合的边界识别法。传统方法计算的斜导数和θ导数如图 5所示。

图  5  传统三分量导数组合的南海边界识别结果

Figure 5.  Edge Detection Results of the South China Sea Based on the Classical Three-Component Gradient Combination

由图 5可以看到，传统的边界识别方法产生额外的错误边界信息，获得的边界不够收敛，会对实际地质解释造成干扰。针对全张量重力梯度数据信息量大、信号频率高，能更好地描述小的异常特征等特点，本文采用全张量梯度数据改进的T_tan值以及改进的T_cos值进行南海地质体边界识别。图 6为全张量梯度组合的边界识别结果。

图  6  全张量梯度组合的边界识别结果

Figure 6.  Edge Detection Results Based on the Combination of Full Tensor Gravity Gradient Data

由图 6可见，改进后的两种算法能有效地识别出地质体线性边界信息，而且两种方法识别的边界一致。与传统方法比较，基于全张量梯度的边界探测器能很好地识别地质体的边界，不仅可以有效地避免产生错误的边界信息，而且获得的边界结果更加准确、收敛。基于T_tan值和T_cos值特征法获取的结果可以清楚地展现出南海NE向的展布特征以及南海海盆区域、中西沙区域、南沙群岛区域、华南陆架陆坡区域、华南大陆区域、苏禄海区域、西部南北弧形区域、吕宋岛弧等主要构造单元的边界。

Canny边界提取算法具有信噪比好、定位精度高、最简响应等优点，可以检测到不易识别的弱边界^[26-27]。基于Canny算法提取断裂平面位置的流程如下。

1）使用自适应高斯函数平滑布格重力异常：

式中，（x，y）为点位坐标；σ为标准差；f（x，y）为原始数据；I（x，y）为平滑后数据；*为卷积运算。

2）计算全张量梯度的幅值和方向。传统Canny算法采用2×2大小的模板来计算梯度幅值和方向。这种方法对噪声敏感，容易丢失真正边界信息。本文在传统方法的基础上利用重力全张量梯度的组合计算梯度的幅值和方向，并且给出了具体的公式。

3）非极大值抑制。即将求解构造边界问题转变为求解幅值局部最大值的问题。若中心处的幅值比沿梯度线方向上的相邻点处的幅值小，则赋零，认为不是边界点；反之，确定为边界点。由此过程生成的边界仍包含许多由噪声而产生的虚假边界。

4）通过高低阈值结合连接边界。高阈值一般为低阈值的2倍。首先，利用高阈值寻找构造边界，由此获取的边界中虚假边界很少，但有间断点；然后，利用低阈值将边界断点处闭合，实现完整边界的连接。结合南海区域地质构造给出南海区域断裂识别结果。

图 7（a）给出了Canny算法提取得到的精确的断裂平面位置，黄色线条为改进的Canny算法提取的断裂位置，底图为${\mathit{T}}_{\mathbf{c}\mathbf{o}\mathbf{s}}$值边界识别结果。

图  7  南海及邻区断裂平面位置分布

Figure 7.  Plane Location of Faults in the South China Sea and Its Adjacent Areas

从图 7可以看出，改进后的Canny边缘检测法可直接获取断裂的平面位置，定位精度高，提取的边界基本由构造边界中心点连接而成。对于每条边界只标识一次，不响应虚假边界，最大限度地抑制了噪声边缘。比较图 7（a）和图 7（b），可以看出，得到的断裂走向基本一致。主要断裂有：F1-1和F1-4为海陆交接边界，F1-2为台西南盆地、珠江口盆地和琼东南盆地的边界，F1-3为穿过莺歌海延伸到万安盆地的断裂，F1-5和F1-10为中央海盆的边界，F1-6和F1-7为海盆大洋中脊，F1-8为中央海盆与周围地块的边界，F1-9为马尼拉海沟俯冲带，F1-10为巴拉望盆地、礼乐盆地、文莱沙巴盆地和南沙海槽盆地的边界，F1-11为北康盆地和曾母盆地的边界。

本文推断南海断裂以NE向走向和NW向走向为主，以NEE向、NW向、EW向和近SN向走向为辅。南海海盆区域呈菱形分布，其东西两侧的断裂构造走向近SN向，其他线性构造走向均为NE-NEE向。南海西部为越东陆架，SN向的线性异常与海岸线大致平行，有一系列平直的阶梯状正断层，具有剪切-拉张特征，由西向东断裂，多为走滑型断裂；南海东部断裂主要呈近SN向分布，与吕宋岛弧、马尼拉海沟展布方向相一致，呈向西凸出的弧形，南海东部构造呈现挤压型边缘特征；南海北部陆缘区线性异常方向大致平行于NE向走向的海岸线，其内有一系列阶梯状断层，为拉张型边缘，反映了该区中生代俯冲增生带构造特征；南海南部构造特征主要呈现NE、NW及向南突出的弧形走向，以海沟俯冲型断裂为主。NE向走向以南沙海槽为主，NW向走向以廷贾断裂以及乌鲁根断裂为主，弧形褶断褶带以武吉米辛等缝合带和主卢帕尔为主。

南海中部呈现三角形的线性异常带大致展示出中央海盆的边界。海盆内部的断裂构造为锯齿状张裂，多为NE向和近EW向走向。NE向张裂为规模较大的岩石圈断裂和地壳断造，构成南海海盆的基本构造轮廓。NE向断裂主要分布在隆起的岛块及其边缘，延伸方向与华南板块的主要断裂平行，对南海海盆形成明显的控制作用。近EW向断裂包括NEE向走向和EW向走向两种。NEE向主要分布在陆坡区、北部陆架，以基底断裂为主。EW向断裂以岩石圈断裂为主，中央海盆区的EW向断裂通常构成海底扩张轴。

近年来，位场数据的特征值被广泛应用于地质体的边界识别，其中，重力梯度全张量的特征值是一个重要方法。本文对重力张量数据的边界识别法进行了深入研究，分析了重力梯度三分量组合法和全张量梯度组合法的优缺点，并将其用于南海异常地质体边界和断裂位置的圈定。

1）全张量重力梯度异常可以获得更多的重力细节特征，为地质解释提供了更多的依据。由于全张量重力梯度信息量大、含有更高频的信号成分，联合不同梯度张量组合能更加准确地进行位场解释。通过比较多种边界识别方法得到的结果，可以看出，利用全张量重力梯度分量组合的T_cos值和T_tan值识别的边界一致，定位精准且收敛更好，提高了描述断裂构造或者地质体边界的水平位置的准确性。

2）基于改进的Canny边界提取算法精确获取南海断裂的平面位置，得到的南海断裂边界与T_cos值和T_tan值结果一致，精确给出了南海构造单元边界。南海海盆边缘特征各不相同，西部为走滑型断裂，东部为海沟俯冲型断裂，南部为挤压型断裂，北部为拉张型断裂。南海断裂以NE向走向和NW向走向为主，以NEE向、NW向、EW向和近SN向走向为辅。