地下目标体重力张量数据的边界识别与增强算法的研究

张 超，王庆宾，黄佳喜，冯进凯

(信息工程大学 地理空间信息学院，河南 郑州 450052)

地下目标体重力张量数据的边界识别与增强算法的研究

张 超，王庆宾，黄佳喜，冯进凯

(信息工程大学 地理空间信息学院，河南 郑州 450052)

边界增强与识别在重力数据处理中占据重要地位，与传统重力异常数据相比,重力张量及其高阶分量对于直接反映异常体的边界具有更高的精度。当异常数据中的所有网格点的值均较低时,通过Sigmoid变换,可以实现高异常值网格数据的拉升,同时压缩低灰度级像素,从而凸显地质体边界,提高边界增强后图像的识别效果。文中利用张量及其分量构建常用的边界识别算法,通过组合体模型进行多种边界识别算法的试算，以比较分析各自的效果,并对结果进行Sigmoid变换。结果表明:对于张量高阶分量组合形式,水平梯度模、解析信号能基本反映浅异常体的边界,gzz水平梯度模能较好反映浅异常体边界,但三者均不能识别深异常体边界;Tilt梯度、Theta和ITA3效果不理想;ITA2能在有效均衡不同强度异常信号的同时,清晰地识别不同深度异常体的边界;采用Sigmoid变换,可以明显提高边界识别的识别效果。

重力张量;边界识别;水平梯度;解析信号;Sigmoid变换

重力勘探方法属于应用地球物理地质勘探与矿藏寻找技术的一个重要分支。由于矿产资源或地下地质构造的密度以及磁性一般同附近的岩体有着显著的不同，因此，重力勘探在寻找矿产资源与地质构造研究中具有重要的应用价值。重力勘探具有快速、经济、范围大等优点，是进行构造划分、异常圈定和发掘矿产资源等快速而有效的手段[1]。地质体边界是地质体的断裂带、相异地质体之间的分界线。在这些位置处，由于重力能在一定程度上表现出数值的不连续性，呈现出重力异常的较大变化。利用这一特点，可以通过重力异常的较大变化率寻找对应的地质体边界线，从而为地下地质结构的分析提供支持。

边界增强与识别在重力数据处理中占据重要地位。原始实测重力异常数据对地下异常体信息反映并不突出，直接根据此数据并不能准确划分出地质体的边界。因此需要对原始异常数据进行相应滤波处理来增强对地质体边界的识别。重力异常水平梯度的极大值和垂直梯度的零值基本反映地质体边界，边界提取算法也是基于此提出的[2]。常见的位场数据边界提取算法主要包括Pedersen等提出的两个张量不变量组合、水平梯度模法(The total horizontal derivatives,THDR)[3]、解析信号法(The analytic signal amplitude,ASA)[4]、Tilt梯度法[5]、Tilt梯度的水平导数法[6]、Theta图法[7]、gzz的水平梯度模法等等，然后利用图中的极值点、零点或其他特征点来识别边界。本文利用张量及其高阶分量构建了常用的边界识别算法，通过组合体模型进行多种边界识别算法的试算，比较分析各自的效果，并对结果进行Sigmoid变换。

1 重力异常及梯度张量

建立一个局部坐标系，以地面某一点O为坐标原点，坐标轴ζ轴沿重力方向垂直向下，ξ，η轴在水准面上，以下所有的计算都在局部坐标系下进行。设异常体剩余密度为σ，则引力位

(1)

重力异常就是引力位沿ζ轴方向的导数，即

(2)

(3)

由于在地球外部，位场满足Laplace方程Uxx+Uyy+Uzz=0 ，以及重力场的无旋性，即Uxy=Uyx，Uxz=Uzx，Uyz=Uzy，因此重力梯度张量中只有5个独立分量。

2 边界增强与识别算法

2.1 边界识别算法

常用的异常数据处理滤波器主要都是基于水平梯度和垂直导数的各种组合来构建，以增强异常地质体的边界特征，从而实现边界识别。与传统重力异常数据相比，重力张量数据对于直接反映异常体的边界具有更高的精度。通过梯度张量及其高阶分量的不同组合，把多个分量的信息集中增强，突出边界信息，从而辅助数据处理和解释。

对于各张量组合形式的表达式见表1和表2。

表1 各张量组合形式的表达式

表2 常见张量高阶组合形式的表达式

2.2 重力张量数据对比度增强

采用对比度增强方法对重力张量网格数据进行对比度拉伸，增强数据的对比度，扩展其中的强异常与弱异常区域的对比效果，可以增强边界识别和检测的效果。这里主要运用Sigmoid变换对网格值进行处理。

Sigmoid变换方程式原型为

(4)

Sigmoid曲线呈现为S状，通过将其作用于异常数据，可以实现高异常值网格数据的拉伸，同时压缩低灰度级像素[8]。这在一定意义上符合异常体边界初步增强后，进一步提高边界对比度的要求，从而凸显地质体的边界，提高边界增强后图像的视觉效果。通过引入合适的参数，可以在一定程度上改变响应曲线的形状，从而使变换算法更加实用。引入参数后的模型为

(5)

式中：fmax为异常数据最大值；I(i,j) 为原位于网格(i,j)的异常值；f(i,j)为网格(i,j)处理结果异常值；α，β为调整变换模型的参数。其中α为用于控制形状的参数，β则用于控制相位。在实际应用时，可以根据当前处理的异常数据的不同特性选择适当的α，β值，从而达到有针对性的改善异常数据分布质量，实现对比度拉升。参数中的α值决定变换函数的陡峭程度，而β值则确定变换函数在坐标水平轴上的平移分量。分析可得：当β取值很小时，可以实现压缩低异常值部分，拉伸高异常值部分。当β取值较大时，则全图的所有网格点的值均会得以拉升，其中高异常值的像素的拉升度要强于低灰阶的拉升度。由该函数性质特点可知，当异常数据中的所有网格点的值均较低时，采用这种方式得到的增强效果比较明显。

3 模型试算

为充分验证表2边界识别算法对地下目标体的识别效果，这里通过建立规则的地下空洞模型进行计算分析。本文所选地面观测数据区域为401×401点距(点距单位为m)的平面网格面积，并通过规则模型正演计算得到重力异常数据。本文布设一个由球体和长方体构成的组合体模型：其中球体半径为30 m，中心点坐标为(-110,-110),上顶面埋深为30 m,与围岩密度差为2.67 g/cm2；长方体的长宽高为120 m×30 m×30 m，中心点坐标为(50,50)，上顶面埋深为15 m，与围岩密度差为2.67 g/cm2。

3.1 常见张量高阶组合形式的边界识别效果分析

分别利用水平梯度模法、解析信号法、Tilt梯度法、Theta图法、gzz的水平梯度模法、改进二阶倾斜角等方法ITA2对重力异常进行处理，结果见图1。

图1 模型重力张量组合形式边界识别效果

从以上试算结果中可以看出：水平梯度模和解析信号能基本反映出浅异常体平面位置和边界，但基本不能反映深异常体的边界；虽然Tilt梯度和Theta图能有效均衡反映不同深度的异常体信号，但Tilt梯度不能反映异常体的边界，而Theta图反映的边界模糊且发散，并且两者均出现边界效应；gzz的水平梯度模能够较好地反映出浅异常体的边界，但不能反映深异常体的边界；改进二阶倾斜角ITA2能有效均衡反映不同深度的异常体信号，识别出的边界略大于真实边界，对异常体边界反映效果最好。

3.2 改进二阶倾斜角ITA2和改进三阶倾斜角ITA3对比

由3.1中试算结果可知，ITA2对不同深度的异常体边界识别效果最好。理论上，重力张量高阶导数对异常体的边界识别具有更高的分辨率，能更加突出异常体边界。这里构建ITA3，并比较ITA2和ITA3对不同深度的异常体模型的边界识别效果。

(6)

试验模型与3.1中模型一样，结果如图2所示。

图2 ITA2和ITA3比较结果

从图2可以看出：与ITA2相比，ITA3所识别的浅异常体的边界形状清晰准确，识别的深异常体边界清晰，但却明显小于真实异常体边界，而且产生较大的噪声和边界效应。总的来说，识别效果较差，明显不如ITA2识别效果，可能是由于高阶导数的计算增大了噪声的干扰。

3.3 Sigmoid变换分析

分别对前文中水平梯度模、解析信号、Tilt梯度、Theta图法、gzz的水平梯度模、改进二阶倾斜角ITA2和ITA3进行Sigmoid变换，结果如图3所示。

图3 Sigmoid变换结果

从图3对比可以看出：水平梯度模、解析信号、Tilt梯度、Theta图法、gzz的水平梯度模、改进二阶倾斜角ITA2和ITA3等异常数据的对比度得到明显的提高，凸显异常体的边界，提高边界识别图像的识别效果；但是从(d)和(f)可以看出，Sigmoid变换在增强有用异常信号的同时也使得数据中的噪声被增大。

4 结 论

本文通过对组合体模型试算分析，全面比较多种张量高阶组合形式的边界识别方法的应用效果，并对各种边界识别算法得到的张量异常数据进行Sigmoid变换以增强数据对比度。分析结果表明：水平梯度模和解析信号能基本反映出浅异常体平面位置和边界，但不能反映出深异常体的边界；Tilt梯度和Theta图能有效均衡反映不同深度的异常体信号，但Tilt梯度不能反映异常体的边界，而Theta图反映的边界模糊且发散，并且两者均出现边界效应；gzz的水平梯度模能够较好地反映出浅异常体的边界，但不能反映深异常体的边界；ITA3所识别的浅异常体的边界形状清晰准确，但由于高阶导数的计算增大噪声的干扰，从而产生较大的噪声和边界效应；改进二阶倾斜角ITA2能有效均衡反映不同深度的异常体信号，识别出的边界略大于真实边界，对异常体边界反映效果最好；当异常数据中的所有网格点的值均较低时，采用Sigmoid变换可以凸显异常体的边界，提高边界图像的识别效果。ITA2和Sigmoid变换能够在异常体的边界识别工作中发挥更大的作用，从而促进重力反演解释工作的顺利进行。

[1] 马国庆．位场(重&磁)及其梯度异常自动解释方法研究[D]．长春：吉林大学，2013．

[2] 吴亮，刘长弘，王庆宾，等．重力异常对地下异常体边界的识别[J]．测绘科学，2015，40(12)：257-262．

[3] 马国庆，李丽丽，杜晓娟．磁张量数据的边界识别和解释方法[J]．石油地球物理勘探，2012，47(5)：815-821．

[4] NABIGHIAN M N.The analytic signal of two dimensional magnetic bodies with polygonal cross section：its properties and use for automated anomaly interpretation [J]．Geophysics，2012，37(3)：507．

[5] 王明，郭志宏，骆遥，等．Tilt-Euler方法在位场数据处理及解释中的应用[J]．物探与化探，2012(2)：126-132．

[6] VERDUZCO B，FAIRHEAD J D，GREEN C M，et al．The meter reader-new insights into magnetic derivatives for structural mapping [J]，Leading Edge，2012 (2)：116．

[7] WIJNS C，PEREZ C，KOWALCZYK P．Theta map：edge detection in magnetic data[J]．Geophysics，2005，70(4)：39-43．

[8] 高云军，刘召芹，万文辉，等.DEM网格间距对重力近区地形改正精度影响验证[J].测绘工程，2016，25(8)：14-18.

[9] 贾鲁，王世忠，王庆宾，等.基于Android平台的两种重力垂直梯度解算方法研究[J].测绘工程，2016，25(7)：57-61.

[责任编辑：张德福]

Research on underground anomaly boundary recognition and enhancement algorithms based on gravitational tensor

ZHANG Chao, WANG Qingbin, HUANG Jiaxi, FENG Jinkai

(School of Surveying and Mapping, Information Engineering University, Zhengzhou 450052,China)

Boundary enhancement and recognition occupies the important position in magnetic data processing. Compared with the traditional gravity anomaly data, gravitational tensor data has a higher accuracy in directly reflecting the boundary of the anomaly. High abnormal values can be pulled up with low gray-scale pixels compressed at the same time by Sigmoid transformation, when all the values of the grid points are lower, thus highlighting geological boundary and enhancing image recognition effect. This paper calculates many kinds of gravity anomaly processing algorithms based on tensor components, including Sigmoid transformation, and then analyzes their effects through combination model. The result shows that: the horizontal gradient，analytic signal and the horizontal gradient of,gzzcan reflect the boundary of the shallow anomaly with the horizontal gradient of,gzza better effect, but none of them can recognize the boundary of deep anomaly. The result of Tilt gradient and Theta diagram method is not ideal, so it is with the ITA3. By contrast, the ITA2can effectively balance the different intensity anomaly and at the same time clearly identify the different depth of anomaly boundary. Sigmoid transformation is proposed, which can distinctly improve the recognition effect of boundary identification.

gravitational tensor;boundary recognition;horizontal gradient;analytical signal;sigmoid transformation

引用著录：张 超，王庆宾，黄佳喜，等.地下目标体重力张量数据的边界识别与增强算法的研究[J].测绘工程，2017，26(4)：16-21.

10.19349/j.cnki.issn1006-7949.2017.04.004

2016-06-09

973国家高科技计划资助项目(6132220202)

张 超(1990-)，男，硕士研究生.

P223

A

1006-7949(2017)04-0016-06