侯振隆^1,2, 王恩德^1,2, 唐水亮³, 袁园⁴
1. 东北大学 资源与土木工程学院, 辽宁 沈阳 110819;
2. 东北大学 深部金属矿山安全开采教育部重点实验室, 辽宁 沈阳 110819;
3. 吉林大学 地球探测科学与技术学院, 吉林 长春 130026;
4. 国家海洋局 第二海洋研究所, 浙江 杭州 310012
收稿日期：2017-10-30
基金项目：中国博士后科学基金资助项目(2017M621151)；中央高校基本科研业务费专项资金资助项目(N160103003)；东北大学博士后基金资助项目(20180313)。
作者简介：侯振隆(1988-)，男，辽宁沈阳人，东北大学师资博士后研究人员;
王恩德(1957-)，男，辽宁盖州人，东北大学教授，博士生导师。

摘要：为了提高重力勘探中边界识别方法的分辨率, 增强对地质体边界的划分和检测能力, 本文在方向总水平导数法的基础上进行了改进.提出了基于重力梯度数据的方向总水平导数与加强方向总水平导数的垂向导数、综合x, y和z方向的方向总水平导数及其垂直导数等6种新型边界探测器.通过理论模型与噪声试验, 证明了改进的方法具有较好的边界划分能力; 与原方法相比, 改进后的方法识别出的边界与模型更好地吻合.将文顿盐丘地区的实测数据引入改进后的方法中, 更清晰地识别出地下盖岩的边界, 验证了方法的可行性.
关键词：重力梯度数据方向总水平导数垂向导数边界识别文顿盐丘
Improvements of the Directional Total Horizontal Derivatives and Their Application in Edge Detection
HOU Zhen-long^1,2, WANG En-de^1,2, TANG Shui-liang³, YUAN Yuan⁴
1. School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China;
2. Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819, China;
3. College of Geo-exploration Science and Technology, Jilin University, Changchun 130026, China;
4. Second Institute of Oceanography, State Oceanic Administration, Hangzhou 310012, China
Corresponding author: HOU Zhen-long, E-mail: houzhenlong@mail.neu.edu.cn
Abstract: An improvement was made on the basis of directional total horizontal derivative to increase the edge detection resolution in the gravity exploration and enhance the ability of dividing and detecting geological body boundary. Six new kinds of edge detectors were presented, including the vertical derivative of directional total horizontal derivative and enhanced directional total horizontal derivative, directional total horizontal derivative and its vertical derivative integrated in x, y and z direction, which are based on gravity gradiometry data. Through the theoretical model and noise test, it is proved that the improved methods have the good ability for boundary division. Compared with the ordinary ones, the detected edges agree well with that from the model. The real data from Vinton Dome are introduced into the proposed methods and the edges of underground cap rock are recognized more clearly, and the feasibility of the proposed methods is verified.
Key words: gravity gradiometry datadirectional total horizontal derivativesvertical derivativeedge detectionVinton Dome
边界识别是重力勘探中重要的解释方法之一, 可用于确定地质体边界、断层线、区分不同物性的岩石, 常用方法有总水平导数法(total horizontal derivative, THD)^[1]、解析信号法(analytic signal, AS)^[2]、倾斜角法(tilt angle, TA)^[3]等.Marson等^[4]通过研究发现利用重力垂向梯度对密度差造成的边界进行定位比重力异常的效果好, 且水平导数结果的分辨率比解析信号、欧拉反褶积的高; 重力梯度数据也被用于成像、反演等工作中, 被证实具有良好的空间分辨能力^[5-6].文献[7-8]在此基础上进一步提出方向总水平导数(directional total horizontal derivative, DTHD)与加强方向总水平导数(enhanced DTHD, EDTHD), 分别结合x, y方向的DTHD和EDTHD可以得到新的边界探测器:总水平导数型边界探测器(edge detector of total horizontal derivative, EDT)和增强型EDT(enhanced EDT, EEDT).这两种方法就是通过重力梯度数据及利用振幅极大值进行边界的检测; 虽然EDT和EEDT能够较准确地识别出地质体边界, 但其分辨率仍待提升.
本文在EDT和EEDT基础上作出改进, 联合x, y和z方向上的DTHD和EDTHD以及垂向导数等方法, 进一步提高边界识别能力; 利用双立方体模型与实测数据试验, 分析得出了6种探测器的识别结果, 以检验其边界检测的水平.
1 改进的方向总水平导数原理1.1 EDT和EEDT的定义式根据DTHD的定义^[7], 在x, y和z方向上的表达式分别为

(1)

(2)

(3)

其中, V表示重力位, 则V_{αβ(α, β=x, y, z)}表示重力位在x, y和z方向上的导数, 即重力梯度数据.式(3)也是THD方法的定义式, 其值等于重力异常在x, y方向上导数平方和的平方根.联合式(1)和式(2)定义边界探测器EDT^[7]:

(4)

对DTHD应用高阶导数可以进一步增强边界识别的分辨率, 得到x, y和z方向上的表达式分别为^[8]

(5)

(6)

(7)

其中, n为大于或等于1的阶数; 当n等于1时, 式(5)~式(7)和式(1)~式(3)相同; 通常情况下n的值为1或2(本文取值为2).和EDT类似, 联合式(5)和式(6)定义边界探测器EEDT^[8]:

(8)

式中, n=2.
式(1)~式(8)均是通过极大值进行边界定位的.
1.2 改进的边界探测器1.1节中EDT和EEDT采取的是联合x, y方向上DTHD和EDTHD的方法, 主要利用了THD_x与ETHD_xⁿ, THD_y与ETHD_yⁿ能够描绘x, y方向上边界的特性; 而z方向上的THD_z及和ETHD_zⁿ可以直接确定边界, 但分辨能力低于EDT和EEDT.因此, 分别联合THD_x, THD_y及THD_z, ETHD_xⁿ, ETHD_yⁿ和ETHD_zⁿ得到改进的EDT(improved EDT, IEDT)和改进的EEDT(improved EEDT, IEEDT):

(9)

(10)

分别对式(4), 式(8)~式(10)计算垂直方向的一次导数(简称垂直导数)以改善识别效果.

(11)

(12)

(13)

(14)

其中, DEDT, DEEDT, DIEDT和DIEEDT分别表示EDT, EEDT, IEDT和IEEDT的导数，式(12)和式(14)是当n=2时的结果; 式(11)~式(14)中使用了重力全张量梯度数据(full tensor gravity gradiometry, FTG)中的V_xz, V_xy, V_yz, V_zz 4个分量.上述公式中隐含了重力位在z方向的三阶导数V_zzz, 为避免计算高阶导数带来的噪声影响, 使用式(15)计算V_zzz:

(15)

综上, 式(9)~式(14)为本文提出的6种边界识别方法, 也是通过极大值定位边界.
2 双立方体模型试验2.1 理论模型为检验方法的有效性, 建立双立方体模型进行试验, 并加入5%的高斯噪声模拟实际情况:测区在x, y方向范围均为0~200 m, 观测点均匀分布于z=0的平面上; 两立方体模型的剩余密度(1 g/cm³)与体积(30 m×30 m×30 m)分别相等, 埋深分别为10 m和15 m; 模型水平位置在图 1, 图 2中用矩形框表示, 其中右上角的立方体埋深较大.
图 1(Fig. 1)

图 1 未加噪声情况下, 不同方法的边界识别结果示意图(矩形框表示立方体边界)Fig.1 In the case without noise, the edge detection results of different methods(the rectangular are the boundaries of prisms) (a)—EDT; (b)—DEDT; (c)—EEDT; (d)—DEEDT; (e)—IEDT; (f)—DIEDT; (g)—IEEDT; (h)—DIEEDT.

图 2(Fig. 2)

图 2 加噪声情况下, 不同方法的边界识别结果示意图(矩形框表示立方体边界)Fig.2 In the case with noise, the edge detection results of different methods(the rectangular are the boundaries of prisms) (a)—EDT; (b)—DEDT; (c)—EEDT; (d)—DEEDT; (e)—IEDT; (f)—DIEDT; (g)—IEEDT; (h)—DIEEDT.

2.2 结果分析对双立方体模型分别使用EDT，EEDT和文中提出的6种方法进行边界识别, 并对比识别的效果(图 1).从结果中发现, 浅部立方体的边界被很好地识别出来, 各方法对边界与角点的识别具有不同的特点, 而对较深的立方体边界仅能识别其范围.分别比较图 1a和1b, 1c和1d, 1e和1f, 1g和1h, 发现计算导数后的方法具有更高的边界识别能力:利用式(11)~式(14)得出的结果能够使边界更好地在模型的真实位置处收敛, 即极大值集中于实际的边界上, 而原方法得出的极大值更多地集中在角点位置; 识别出的极大值与背景值的差值更大, 使边界能够更清晰地在图中呈现.在图 1b, 图 1d~图 1h中, 基于THD的DEDT, IEDT和DIEDT方法的结果图像呈光滑分布; 而基于ETHD的DEEDT和DIEEDT方法在立方体的四周产生了极小值, 即高阶导数对结果产生了一定的影响; 基于ETHD的方法能够更好地识别较深地质体的边界; 在本次试验中, IEEDT的识别效果最好.
对包含噪声的数据利用高斯滤波器进行去噪, 采取相同的方法实行边界识别, 结果见图 2.和未加噪声的试验结果相似, 各方法对浅部模型识别效果最好, 但噪声与去噪处理的引入使背景值出现“斑状噪点”.通过结果对比发现, 提出的6种方法识别效果优于EDT和EEDT, 极大值更集中在立方体边界上; 该次试验中, DEEDT和DIEEDT结果中, 立方体周围没有出现明显的极小值虚假异常, 由于局部极大值更集中在深部立方体边界附近, 认为DIEEDT的效果最好, IEEDT和DIEDT等方法次之.
通过数据试验证明了提出的6种边界识别方法能够较好地识别出立方体边界, 并能同时识别出不同埋深的立方体位置; 和原方法相比, 具有更好的识别效果.
3 文顿盐丘实测数据试验将上述方法应用于美国文顿盐丘地区的实测航空FTG数据.文顿盐丘是由岩盐及其上方的盖岩组成的, 观测的梯度数据主要来自于盖岩^[9-10].许多学者对该区域的地下密度分布进行了反演研究^[10-13], 计算的盖岩剩余密度约为0.55 g/cm³.
数据位于WGS84坐标系, 令x轴正方向为东、y轴正方向为北, 选取x, y方向范围分别为440.50~444.50 km, 3 332.78~3 336.78 km的区域进行研究; 飞行测量时的平均高度为85 m; 数据已经过密度为2.2 g/cm³的地形改正, 在边界识别前仍需利用高斯滤波器进行噪声压制.根据式(9)~式(14), 选取FTG数据中的V_xz, V_xy, V_yz和V_zz进行计算.边界识别的结果见图 3.结合其他学者的研究^[9-11], 认为红色的高值区域表示盖岩的位置.极大值区域的边界更为清晰, 但是蓝色的背景值区域呈现出斑状的噪声现象, 说明在实际工作中, 过高阶数的导数不能取得最理想的效果, 这是由于在地球物理数据处理中求高阶导数相当于对数据滤波, 导数阶数过高会引入其他频率的干扰信息.分别对比图 3a和图 3e, 图 3c和图 3g, 加入z方向数据后改进方法在识别能力上有一定的改善, 计算的结果区间增加, 使极大值与背景值的差值增加, 增强了对边界的显示效果.综上, 认为IEEDT在文顿盐丘的数据试验中能够更准确地识别出盖岩边界.从图 3g中可发现盖岩呈弧状分布, 并主要集中于研究区域的中部, 根据Hou等对密度的反演结果^[12-13]可知, 盖岩的密度高值区位于盖岩的东南部; 初步判断盖岩在x, y方向的分布范围分别为441.80~443.00 km，3 333.80~3 334.80 km.需要注意的是, 已知研究成果显示盖岩顶部形状呈东南高、西北低, 即东南部的埋深较西北部的浅, 这和红色区域中的极大值位置对应, 说明该处表示的是盖岩的浅部.
图 3(Fig. 3)

图 3 实测数据试验结果Fig.3 Results of real data test (a)—EDT; (b)—DEDT; (c)—EEDT; (d)—DEEDT; (e)—IEDT; (f)—DIEDT; (g)—IEEDT; (h)—DIEEDT.

4 结论本文在方向总水平导数的基础上, 为EDT和EEDT方法引入z方向的分量, 并分别计算了各自的垂向导数, 提出了IEDT等6种边界检测器.在双立方体模型和文顿盐丘实测数据试验中证实了以上方法的有效性和可行性, 能够较为准确地识别出模型或盖岩的边界位置, 可用于确定地下构造或地质体边界等.在试验中还验证了过高阶数的导数对结果的影响, 在实际数据解释中应根据情况选择适当的方法.
致谢
感谢Bell Geospace Inc提供的实测FTG数据.
参考文献

[1]	Cordell L, Grauch V J S.Mapping basement magnetization zones from aeromagnetic data in the San Juan Basin, New Mexico[C]// 1982 SEG Technical Program Expanded Abstracts.Tulsa: Society of Exploration Geophysicists, 1982: 246-247.
[2]	Roest W R, Verhoef J, Pilkington M. Magnetic interpretation using the 3-D analytic signal[J].Geophysics, 1992, 57(1): 116–125.DOI:10.1190/1.1443174
[3]	Miller H G, Singh V. Potential field tilt-a new concept for location of potential field sources[J].Journal of Applied Geophysics, 1994, 32(2/3): 213–217.
[4]	Marson I, Klingele E E. Advantages of using the vertical gradient of gravity for 3-D interpretation[J].Geophysics, 1993, 58(11): 1588–1595.DOI:10.1190/1.1443374
[5]	Guo L, Meng X, Shi L. 3D correlation imaging of the vertical gradient of gravity data[J]. Journal of Geophysics and Engineering, 2011, 8(1): 6–12.DOI:10.1088/1742-2132/8/1/002
[6]	Zhdanov M S, Ellis R, Mukherjee S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J].Geophysics, 2004, 69(4): 925–937.DOI:10.1190/1.1778236
[7]	Yuan Y, Geng M.Directional total horizontal derivatives of gravity gradient tensor and their application to delineat the edges[C]// 76th EAGE Conference and Exhibition.Amsterdam, 2014: 1-3.
[8]	袁园.全张量重力梯度数据的综合分析与处理解释[D].长春: 吉林大学, 2015. ( Yuan Yuan.Comprehensive analysis, processing and interpretation of the full tensor gravity gradient data[D].Changchun: Jilin University, 2015.http://cdmd.cnki.com.cn/Article/CDMD-10183-1015587927.htm)
[9]	Ennen C, Hall S.Structural mapping of the Vinton salt dome, Louisiana, using gravity gradiometry data[C]// 2011 SEG Annual Meeting.Tulsa: Society of Exploration Geophysicists, 2011: 830-835.
[10]	Oliveira V C Jr, Barbosa V C F. 3-D radial gravity gradient inversion[J].Geophysical Journal International, 2013, 195(2): 883–902.DOI:10.1093/gji/ggt307
[11]	Geng M, Huang D, Yang Q, et al. 3D inversion of airborne gravity-gradiometry data using cokriging[J].Geophysics, 2014, 79(4): 37–47.DOI:10.1190/geo2013-0393.1
[12]	Hou Z, Huang D, Wei X. Fast inversion of probability tomography with gravity gradiometry data based on hybrid parallel programming[J]. Journal of Applied Geophysics, 2016, 124: 27–38.DOI:10.1016/j.jappgeo.2015.11.009
[13]	Hou Z, Huang D. Multi-GPU parallel algorithm design and analysis for improved inversion of probability tomography with gravity gradiometry data[J].Journal of Applied Geophysics, 2017, 144: 18–27.DOI:10.1016/j.jappgeo.2017.06.009