边界识别在地球物理位场数据解释中占有重要的位置.目前有许多方法用来进行地质体的边界检测.其中大多数都是使用重力梯度的部分分量来定义边界探测器，如总水平导数THD和解析信号振幅(Cordell, 1979; Cordell and Grauch, 1985).Marson和Klingele(1993)指出总水平导数识别结果的分辨率优于解析信号振幅.然而，当多个地质体位于不同埋藏深度时，这些方法仅能对浅部异常边界识别较为准确，对深部异常识别效果较差.Miller和Singh(1994)使用垂向导数与总水平导数的比来定义倾斜角边界识别法，其零等值线对应地质体边界.它是第一个被提出的均衡滤波器，能均衡不同振幅强度的异常信息.Verduzco等(2004)用倾斜角的总水平导数THDR来进行边界检测.Wijns等(2005)提出了Theta图滤波器进行边界探测，它是用解析信号振幅来对总水平导数进行归一化处理.Cooper和Cowan(2006)总结了多种不同形式的边界识别方法，提出了TDX边界识别滤波器进行异常体边界识别.这些均衡边界探测器虽然都能同时显示不同振幅大小的异常的边界信息，但识别出来的边界分辨率较低.Ma和Li(2012)通过对总水平导数的归一化，能同时获得不同振幅大小的异常边界，且分辨率有所提高.然而，当测量数据中同时包含正负异常体时, 上述均衡边界识别滤波器就可能带来错误的边界信息.结构张量特征值是一种图像处理技术，其呈现的是一个n维空间的局部方向(Weickert, 1999a; 1999b).近年来，该方法被应用于圈定地质体构造边界，其最大特征值描述地质体边界，最小特征值描述地质体的角.Sertcelik和Kafadar(2012)利用该技术进行重力异常边界和角的提取.该方法在结构张量中使用了高斯包络，可以有效降低干扰噪声影响，虽然识别出来的边界结果更加清晰，并且该方法可以消除正负密度出现的假异常，但不能均衡不同大小的振幅异常.

南海地处太平洋板块、欧亚板块和印——澳板块之间(Leloup et al., 1995)，地质构造复杂.许多学者对这一特殊的地质构造单元进行了研究(刘光鼎，1992；李家彪等, 2011a, 2011b；吴招才等，2017；于传海等，2017)，但仍存在许多的问题和争议.其中包括红河断裂在南海后续延伸问题、古南海缝合带位置问题和南海地质构造单元划分问题.

红河断裂是介于华南板块和印支板块之间的一条断裂.起始于我国的青藏高原地区，中间经过我国云南省与越南部分地区，最终延伸到南海的莺歌海处.其走向为ES向，长达1000多公里.但是该断裂在进入南海之后，其走向问题一直是大家争议的焦点.姚伯初等(1994)认为在西沙群岛的西永一井获得的岩石基底与印支半岛的昆嵩隆起处基底都是前寒武纪变质岩基底，指出红河断裂与西沙群岛属于同一地块，从而认为红河断裂向东延伸；刘宝民等(2006)认为红河断裂与南海西缘断裂都是走滑性质断裂，且与Tapponnier等(1982, 1986)提出的碰撞-挤出-逃逸模式相吻合，故而认为红河断裂向南延伸；万玲等(2000)认为南海的中建南断裂与红河断裂在走向均是NW向，而南海断裂的整体呈NE向，故这两条断裂是同一构造性质作用下产生的，并且中建南断裂的莫霍面等值线与红河断裂的莫霍面等值线刚好一致，因此提出红河断裂向ES延伸.

南沙海槽处于加里曼丹的北部，呈NE走向.关于南沙海槽是否为古南海俯冲带位置学术界有不同的观点.第一种观点，Hall(1996)提出南沙海槽是一个大陆边缘，而古南海正是发育在其南部，后来由于礼乐地块从华南陆块脱离往南移动，使得古南海消亡在南沙海槽处，且基底为洋壳；Hinz(1985, 1989)提出南沙海槽并不是古南海的俯冲带，而是巴拉望仰冲楔状体推覆前缘，基底属于陆壳性质.

对于南海地质构造分区，从1980年以来，许多地质学家在不同的角度对其展开了相应的研究.金庆焕(1989)以不同年代地壳基底作为划分依据，将南海以及周围地区分成了7个古老基底的微板块区，以及性质不同的34个块体；刘昭蜀等(2002)将南海周边大陆边缘的地堑系分成了5个二级构造单元和24个三级构造单元；宋海斌等(2002)同样将南海海域划分为了四块，分别为：南海北缘、中西沙、南沙以及南海海盆四个块体；姚伯初(2006)根据岩石圈不同性质将南海看成南海圈一个二级单元块体，然后将其划分了4个三级单元块体；鲁宝亮等(2015)通过南海地壳类型、构造层序列和沉积建造将南海分为了8个二级构造单元和9个三级构造单元.

本文针对上述问题，提出一种新的重力边界识别方法Tilt-Eigen法，利用结构张量最大特征值定义Tilt angle函数实现.该方法不仅能降低噪声影响，均衡不同振幅大小异常的边界信息，而且还能避免产生额外的虚假边界异常.将Tilt-Eigen方法用于南海实际重力数据处理，获得了20条比较清晰的边界结果.根据这些边界结果，本文推测出红河断裂入海后的走向以及延伸情况；分析了古南海俯冲带位置；并对南海地区进行地质构造单元划分.

重力梯度张量矩阵为引力位G在三维空间x，y，z三个方向上的二阶偏导数，张量矩阵形式：

(1)

原始结构张量是由高斯包络和引力位水平梯度张量卷积而成(Sertcelik and Kafadar, 2012)，其表达式为：

(2)

其中，

(3)

δ_x和δ_y是高斯包络在x和y方向上的标准差，.结构张量矩阵M的齐次方程为：

(4)

最大特征值为：

(5)

最小特征值为：

(6)

最大特征值λ₁的最大值描述地质体的边界，最小特征值λ₂描述地质体的棱角.

基于上述方法存在的一些缺陷，本文做出了一些改进，利用结构张量最大特征值定义Tilt angle函数，称为Tilt-Eigen，其表达式为:

(7)

其中，Tilt-Eigen的最大值对应地质体的边界.根据反正切三角函数的特性，其取值范围是-π/2到π/2.由于结构张量矩阵的特征值是非调和函数，其最大特征值的垂向导数通过Florio等(2006)提出垂直导数算子计算实现，Matlab函数为dzder.m.水平导数则直接通过Matlab内置函数gradient.m计算得到.

为了验证本文Tilt-Eigen方法的准确性，我们将其与一些传统的方法进行对比分析.它们是总水平导数THD、TDX(Cooper and Cowan, 2006)、倾斜角总水平导数TAHG(Ferreira et al., 2013)和拉普拉斯算子滤波器(王明等, 2015).

建立一个深度分别为10 m和20 m的两个棱柱体模型，其厚度均为40 m，剩余密度为0.2 g·cm^-3.模型的重力梯度分量G_zx，G_zy，G_zz异常见图 1a到图 1c.图 1d-1i已分别展示的是THD、λ₁(δ_x=δ_y=0.5)、TDX、TAHG、拉普拉斯算子和Tilt-Eigen的边界识别结果，其中黑色虚线即为两个棱柱体的真实水平位置.

图 1

Fig. 1

图 1 理轮模型1的重力梯度分量和不同探测器所得的边界 (a)水平导数Gzx；(b)水平导数Gzy；(c)垂直导数Gzz；(d) THD边界结果；(e) λ1边界结果；(f) TDX边界结果；(g) TAHG边界结果；(h)拉普拉斯算子边界结果；(i) Tilt-Eigen边界结果. Fig. 1 Gravity gradients of Synthetic model 1 and edge results identified by different methods (a) Horizontal gradient Gzx; (b) Horizontal gradient Gzy; (c) Vertical gradient Gzz; (d) THD edge result; (e) λ1 edge result; (f) TDX edge result; (g) TAHG edge result; (h) Laplacian edge result; (i) Tilt-Eigen edge result.

从图 1d和1e可以看出，THD和λ₁均可清晰地识别出浅部异常幅值大的边界，但对于深部地质体引起的小异常振幅，识别出来的边界结果相对模糊.因此，这两种边界识别滤波器不能有效的均衡不同振幅大小的异常.然而，从图 1f-1i可以看出，TDX、TAHG、拉普拉斯算子以及Tilt-Eigen均可很好地均衡深部和浅部的异常，其中TDX、TAHG和Tilt-Eigen通过最大值来圈定异常体边界，而拉普拉斯算子则通过其零等值线圈定.图 1h中的黑色实为拉普拉斯算子识别结果的零等值线.尽管这些方法都能均衡不同振幅大小的异常，但是他们的分辨率却有所不同.从图 1f中可以看出TDX圈定的边界位置超出了真实边界范围，这主要是因为棱柱体的埋深相对其宽度较大.众所周知，随着埋深的增加，重力异常会变得越来越小，振幅也会越来越光滑，这就会使得测量边界大于真实边界.而TAHG、拉普拉斯算子和Tilt-Eigen识别的边界结果与真实边界结果一致性较好，更加准确.

为了更加直观的不同方法的分辨率特征，提取过场源的中心的剖面进行单独分析.剖面位置见图 1a中的黑色实线，结果见图 2.图 2中的垂直黑色虚线表示两个棱柱体的真实边界，而图 2e中的水平实线表示零等值线.可以看出，THD和λ₁的最大值可以准确地圈定边界，但对于深部异常的边界结果振幅较小，在图形显示过程中会被浅部异常边界压制致使结果相对模糊.对于传统的TDX，其最大值的宽度与地质体真实边界差别较大，识别的边界宽度较大，如图 2c.而对于TAHG、拉普拉斯算子和Tilt-Eigen，其识别的边界宽度与真实宽度一致性较好.

图 2

Fig. 2

图 2 图一模型的剖面结果 (a) THD边界结果；(b) λ1边界结果；(c) TDX边界结果；(d) TAHG边界结果；(e)拉普拉斯算子边界结果；(f) Tilt-Eigen边界结果; (g—h)模型横截面. Fig. 2 Profile results of model 1 shown in Fig. 1 (a) THD edge result; (b) λ1 edge result; (c) TDX edge result; (d) TAHG edge result; (e) Laplacian edge result; (f) Tilt-Eigen edge result; (g—h) Cross section of model.

为了更好的验证Tilt-Eigen的探测能力，建立了一个更加复杂的地质体模型，包含三个埋深度不同、大小完全相同的棱柱体组成，埋藏深度分别为4 km、1 km和2.5 km，厚度均为1 km.其中，1号和2号棱柱体的剩余密度为0.2 g·cm^-3，3号棱柱体的剩余密度为-0.2 g·cm^-3.图 3a和图 3b分别为模型2的平面图和3D图.模型2的重力梯度分量G_zx、G_zy、G_zz异常见图 4a-4c.图 4d-4i为不同方法识别的边界结果.图中白色的虚线代表棱柱体的真实边界.可以看出，图 4d和图 4e所表示的THD和λ₁(δ_x=δ_y=0.5)都不能清晰地展现深部异常边界.而图 4f和图 4h中的TDX和拉普拉斯算子虽然可以清晰圈定边界，但结果产生了额外的虚假异常边界.这主要是因为同时存在的正负异常体之间会产生额外的零等值线.图 4g和图 4i中TAHG和Tilt-Eigen既可以清晰地圈定边界，且没有引入虚假异常边界.对比上述结果可以看出，TAHG和Tilt-Eigen可以获得更好的效果，应用范围相对更广.

图 3

Fig. 3

图 3 模型2的平面图和3D图 (a)平面图；(b) 3D图. Fig. 3 Plan view and 3D view of the model 2 (a) Plan view; (b) 3D view.

图 4

Fig. 4

图 4 模型2的重力梯度分量及其边界结果 (a)水平导数Gzx；(b)水平导数Gzy；(c)垂直导数Gzz；(d) THD边界结果；(e)λ1边界结果；(f) TDX边界结果；(g) TAHG边界结果；(h)拉普拉斯算子边界结果；(i) Tilt-Eigen边界结果. Fig. 4 Gravity gradients of Synthetic model 2 and edge results identified by different methods (a) Horizontal gradient Gzx; (b) Horizontal gradient Gzy; (c) Vertical gradient Gzz; (d) THD edge result; (e) λ1 edge result; (f) TDX edge result; (g) TAHG edge result; (h) Laplacian edge result; (i) Tilt-Eigen edge result.

为了进一步验证方法对噪声的抗干扰能力，本文对模型2的每个梯度分量分别加入其最大异常振幅的10%的高斯随机噪声.图 5a-5c为加噪声后的重力梯度分量异常图.图 5d-5i分别为THD、λ₁(δ_x=δ_y=1.5)、TDX、TAHG、拉普拉斯算子和Tilt-Eigen针对噪声数据的边界识别结果.这里针对高斯包络参数δ_x和δ_y的选择主要采用多组参数进行处理，最终选取最优解.通过与图 4对比可以看出，传统方法THD和λ₁受噪声影响较小.而传统的均衡滤波器TDX、TAHG和拉普拉斯算子所识别的边界结果受噪声影响较大，主要因为它们都使用了重力梯度张量的高阶导数，放大了噪声的影响.而本文方法Tilt-Eigen因使用了高斯包络，起到了滤波效果，可以很好的降低噪声影响.

图 5

Fig. 5

图 5 理论模型2乘以10%的高斯噪音后的结果 (a)水平导数Gzx；(b)水平导数Gzy；(c)垂直导数Gzz；(d) THD边界结果；(e) λ1边界结果；(f) TDX边界结果；(g) TAHG边界结果；(h)拉普拉斯算子边界结果；(i) Tilt-Eigen边界结果. Fig. 5 The result of synthetic model 2corrupted by 10% Gaussian noise (a) Horizontal gradient Gzx; (b) Horizontal gradient Gzy; (c) Vertical gradient Gzz; (d) THD edge result; (e) λ1 edge result; (f) TDX edge result; (g) TAHG edge result; (h) Laplacian edge result; (i) Tilt-Eigen edge result.

为了验证方法在实际中的应用效果，本文对中国南海及其周边地区(104°E-122°E、2°N-24°N)的重力异常数据进行处理.重力数据来源于美国加利福尼亚大学斯克里普斯海洋研究所提供的全球高分辨率测高卫星空间重力异常，网格间距为5 km，http://topex.ucsd.edu/cgi-bin/get_data.cgi.在NOAA National Geophysical Data Center(NGDC)获得了对应区域的地形数据和沉积层厚度数据，数据网隔大小为5′×5′.数据投影方式为WGS84中的墨卡托投影.在此基础上，笔者利用geosoft软件对南海地区的陆地和岛屿做地形校正和中间层校正，而对海洋区域进行相应的海水层和沉积层校正，实现对整个区域的布格重力校正，异常图见图 6.

图 6

Fig. 6

图 6 中国南海地区的布格重力异常 Fig. 6 The bouger gravity of south China sea

针对南海地区布格重力异常进行边界识别，其对应的重力梯度G_zx，G_zy，G_zz采用数值计算获取，见图 7a-7c.图 7d-7i分别为THD、λ₁(δ_x=δ_y=0.5)、TDX、TAHG、拉普拉斯算子和Tilt-Eigen圈定的边界结果.可以看出，THD和λ₁不能清晰的反映深部的构造信息；TDX、TAHG和拉普拉斯算子虽然可以均衡深部和浅部构造信息，但其结果相对发散，且受噪声干扰较大.而本文方法Tilt-Eigen不仅可以均衡深部和浅部构造特征，所获得边界结果清晰、准确，收敛较好，受噪声干扰小.

图 7

Fig. 7

图 7 中国南海地区不同方法的边界结果 (a)水平导数Gzx；(b)水平导数Gzy；(c)垂直导数Gzz；(d) THD边界结果；(e) λ1边界结果；(f) TDX边界结果；(g) TAHG边界结果；(h)拉普拉斯算子边界结果；(i) Tilt-Eigen边界结果. Fig. 7 The edge results of south China sea (a) Horizontal gradient Gzx; (b) Horizontal gradient Gzy; (c) Vertical gradient Gzz; (d) THD edge result; (e) λ1 edge result; (f) TDX edge result; (g) TAHG edge result; (h) Laplacian edge result; (i) Tilt-Eigen edge result.

结合南海区域地质图，对本文识别的边界结果进行划分和编号，见图 9.其中1号、5号、17号和18号分别表示华南大陆、印支大陆和婆罗洲与其临近海洋的重力梯度带边界；2号和3号表示珠江口盆地和台西南盆地内两条NE走向的大型断裂，在图 8中可以看出这两条断裂后期被多次切割，说明此处构造活动强烈；4号表示珠江口盆地台西南盆地与临近莺歌海盆地和琼东南盆地等的重力梯度带边界；6号为穿越莺歌海往南延伸到万安盆地的一条大型走滑断裂；7号和11号表示中央海盆与周围地块的边界，具有很强的重力梯度变化；8号表示菲律宾地块往西运动的俯冲带(马尼拉海沟俯冲带)；9号和10号分别表示中央海盆与西南次海盆的大洋中脊；12号表示北康盆地、南沙海槽盆地、礼乐盆地与曾母盆地、文莱沙巴盆地和巴拉望盆地的重力梯度带边界；13号和14号为苏禄海与南沙之间的两条NE向断裂；15和16号表示苏禄海和苏维拉西海的边界，同样具有很强的重力梯度变化；19号表示曾母盆地与其他地块的重力梯度带边界；20号为廷贾断裂，呈NWW走向.

图 8

Fig. 8

图 8 结合南海底图的Tilt-Eigen边界结果 Fig. 8 Tilt-Eigen edge results on map of south sea

从图 9中可以看出，在南沙海槽东南部存在着一条边界线.其起始于文莱沙巴北部的西南端，朝NE向延伸，在文莱沙巴盆地的中间，突然往正北转变到达南沙海槽南部的东北端，继而继续朝NE向延伸至北巴拉望盆地西北侧.并且从图 6的布格重力异常中也可以显示，在这条边界的两侧梯度变化明显，并且大于南沙海槽西南端的梯度变化.此外，该区域重力最大值可达260 mGal，分布在文莱沙巴与南沙海槽之间.这和姚伯初提出的南沙海槽西南部为洋壳，东北部陆壳相吻合.因此本文推断南沙海槽西南侧的并不是古南海缝合带，缝合带应该位于文莱沙巴盆地附近，而南沙海槽的东北部分则为古南海缝合带的位置.

图 9

Fig. 9

图 9 南海区域主要边界结果 1号、5号、17号和18号表示的是陆地与海洋的交接边界；2号和3号表示珠江口盆地和台西南盆地内的两条大的断裂；4号表示珠江口盆地、台西南盆地与临近莺歌海盆地和琼东南盆地等的边界；6号为穿越莺歌海往南延伸到万安盆地的一条大型断裂；7号和11号表示中央海盆与周围地块的边界；8号表示马尼拉海沟俯冲带；9号和10号表示南海海盆大洋中脊；12号表示北康盆地、南沙海槽盆地、礼乐盆地与曾母盆地、文莱沙巴盆地和巴拉望盆地的边界；13号和14号为苏禄海与南沙之间的两条断裂；15和16号表示苏禄海和苏拉威西海的边界；19号表示曾母盆地与巽它地块的边界；20号为廷贾断裂. Fig. 9 The main edge results of the south China sea Line 1、line 5、line 17 and line 18 represents the boundary between land and sea; line 2 and line 3 represent two major faults in the Pearl river estuary basin and the southwestern Taiwan basin; line 4 indicates the boundary between the southwest Taiwan basin, the Pearl river estuary basin and the adjacent Yinggehai basin and southeast Hainan basin; line 6 is a large fracture extending southward from Yinggehai to the wanan basin; line 7 and line 11 represent the boundary between the central basin of the south China sea and the surrounding land; line 8 represents the subduction zone of the Philippine plot to the west; line 9 and line 11 represent the mid-ocean ridge of the central basin; line 12 indicates the boundary between the Beikang basin, the Nansha trough basin, the Liyue basin and the Zengmu basin, the Brunei Sabah basin and the Palawan basin; line 13 and line 14 are the two faults between the Sulu sea and Nansha plot; line 15 and line 16 indicate the boundary between the Sulu sea and Suvila sea; line 19 represents the boundary between the Zengmu basin and other plots.line 20 is Tinja′s fault.

另外从图 8中还可以看出，在南海西北角存在了一条从莺歌海起始朝ES方向延伸的断裂，其在海南岛南部附近突然转为朝正南方向延伸.本文推断为红河断裂延伸到海南岛南部附近以后，开始往南进行延伸，进而与越东万安断裂相连接.结再合图 6南海布格重力异常，发现红河断裂在南海部分的两侧重力异常值明显不同，其中，NE侧的布格重力异常值显示偏高，变化较为急剧；WS侧的布格重力异常则整体偏低，甚至出现很多负异常，且变化也较为缓慢.而在越东万安断裂两侧，同样可以观察到两侧具有不同的布格重力异常，其E侧为高值重力异常，变化急剧，从80 mGal左右突然增加到了280 mGal；而其W侧的布格重力异常则是偏低，也出现大量的负异常，变化较为宽缓.综合本文识别的边界结果和南海布格重力异常分析，认为越东万安断裂才是红河断裂之后的延伸断裂.

根据图 8识别出的边界结果，本文对南海地区进行构造划分，分成8个地块，见图 10.华南地块为中国南部内陆的一个巨大稳定区域，其主要处于陆壳基底之上，位于1号边界的北部.琼台地块则是在华南地块与南海洋盆之间的一个过渡带，包含了北部湾盆地、珠江口盆地和台西南盆地等多个大型沉积盆地，位于在1号和4号边界内.菲律宾地块位于南海与太平洋板块之间的一个地块，由菲律宾海的大洋岩石圈和菲律宾大陆基底共同构成，其与南海交界处有一个板块俯冲带(马尼拉海沟俯冲带)，位于8号边界东侧，在图 8中可以很清晰的显示.中西沙和南沙地块位于南海海盆两侧，由于新南海的张开而分离开来，故而将中西沙和南沙地块划分为一块，其位于4号、6号、7号、12号和20号边界内.中央海盆地块则是由于新南海张开后所形成的菱形海盆区，主要以大洋地壳为基底，位于7号、8号和11号之间.印支地块是由印支半岛的稳定陆块所组成，位于中西沙南沙地块的西部，华南地块的西南部，在5号边界的西侧.苏禄海地块是由两个新生成的海(苏禄海和苏拉威西海)及其周边盆地所组成，地处菲律宾地块南部，中西沙南沙地块东部，位于14号和17号边界的东南侧.巽它地块则是由巽它大陆及其增生带和临近海域组成，主要是以新生代沉积为主，位于中西沙南沙地块的南部.

图 10

Fig. 10

图 10 南海区域构造划分 Ⅰ区为华南地块；Ⅱ区为琼台地块；Ⅲ区为菲律宾地块；Ⅳ区为中央海盆地块；Ⅴ区为中西沙南沙地块；Ⅵ区为印支地块；Ⅶ区为苏禄海地块；Ⅷ区为巽它地块. Fig. 10 Regional tectonic division of the south China sea Area Ⅰ is south China land block; Ⅱ area for Hainan and Taiwan block; Ⅲ area for the Philippines plot; Ⅳ area for the central basin of the south China sea land; Ⅴ area for the Zhongsha, Xisha and Nansha block; Ⅵ area for the Indosinian block; Ⅶ area for the Sulu sea block; Ⅷ area for Xunta block.

本文通过对重力梯度结构张量特征值求取倾斜导数定义了一种新的边界识别滤波器Tilt-Eigen.通过对不同埋藏深度的棱柱体模型、同时包含正负密度异常的不同深度的棱柱体模型、含噪声的不同深度的正负密度棱柱体模型进行试验，并将其应用到南海实际重力数据处理与解释中.本文得出以下结论：

(1) Tilt-Eigen法比传统方法能准确地圈定地质体边界，不会出现边界发散的情况；能均衡不同深度异常体的振幅；在有正负密度地质体同时存在的情况下，该方法不会产生虚假的异常边界；由于在计算特征向量的算法中引用了高斯包络，使其对噪音的影响比传统边界探测器显著降低.

(2) 推断出南沙海槽的东北部分为古南海缝合带的位置，而南沙海槽西南侧的并不是古南海缝合带，且俯冲带应该位于文莱沙巴盆地附近.

(3) 针对红河断裂在南海海域后续延伸的问题，推断出红河断裂往东南延伸到海南岛南部(109°E, 16°N)之后，开始往南延伸，与越东万安断裂相连接.

(4) 将南海地区分为华南地块、琼台地块、菲律宾地块、中央海盆地块、中西沙南沙地块、印支地块、苏禄海地块以及巽它地块这8个块体.

感谢审稿专家和编辑部的大力支持.感谢黄大年老师生前指导和对本文的巨大贡献.