侯振隆^1,2, 王恩德^1,2
1. 东北大学 资源与土木工程学院, 辽宁 沈阳 110819;
2. 东北大学 深部金属矿山安全开采教育部重点实验室, 辽宁 沈阳 110819
收稿日期：2017-09-18
基金项目：国博士后科学基金资助项目(2017M621151)；中央高校基本科研业务专项资金资助项目(N160103003)；东北大学博士后基金资助项目(20180313)。
作者简介：侯振隆(1988-)，男，辽宁沈阳人，东北大学师资博士后;
王恩德(1957-)，男，辽宁盖州人，东北大学教授，博士生导师。

摘要：为了提高重力勘探中数据成像的效率, 对地球物理勘探中的重力异常数据进行快速三维成像.在相关成像的理论基础上针对立方体元提出利用泰勒级数的算法, 重新定义了异常值的计算模型, 从而得到新的几何函数.该成像方法通过级数展开、积分等方式, 减少了计算量和迭代次数.在理论模型试验中, 证明了提出的方法具有良好的成像能力和抗噪性, 并通过效率分析说明了该算法能够大幅缩短几何函数矩阵的计算时间, 提高成像效率.对文顿盐丘地区的实测重力异常数据进行快速成像, 地质体的位置能够被较好地显示出来, 验证了算法的可行性.
关键词：重力异常数据几何函数快速成像泰勒级数计算效率
Fast Probability Tomography of Gravity Anomaly Data Based on Taylor Series
HOU Zhen-long^1,2, WANG En-de^1,2
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
Corresponding author: HOU Zhen-long, E-mail: houzhenlong@mail.neu.edu.cn
Abstract: In order to improve the data tomography efficiency in the gravity exploration, gravity anomaly data from the geophysical exploration were processed by fast three-dimensional tomography. The algorithm of Taylor series was used for the prism unit based on the related probability tomography theory, in which the calculation model of anomaly data were redefined to obtain new geometric functions. The tomography method reduces the calculation amount and the number of iterations by series expansion and integral. In the theoretical model test, it is proved that the proposed method has good tomography ability and anti-noise characteristic. The efficiency analysis shows that the algorithm can greatly shorten the calculation time of the geometric function matrix and improve the imaging efficiency. Fast tomography was applied for the real gravity anomaly data of Vinton Dome, the location of geological body could be shown well, which verifies the algorithm feasible.
Key words: gravity anomaly datageometrical functionfast tomographyTaylor seriescomputational efficiency
相关成像是由Patella提出的一种位场数据成像方法^[1-2], 该方法对相关系数进行成像, 即利用观测数据及其对应几何函数的相关程度来反映地质体分布情况, 能够取得良好的解释效果; 还可应用于不同类型的观测数据, 例如重力、磁异常及其梯度数据^[3-6]、电阻率异常数据^[7]、电磁异常数据^[8]等; 相关系数还被一些学者用于物性参数的反演^[9-10].在相关系数的计算中需要利用几何函数矩阵, 其规模等于观测点数与地下划分模型数的乘积; 当观测数据维度较大时, 计算该矩阵将占用大量时间.目前的研究主要通过并行计算的手段来缩短几何函数矩阵的运算时间, 例如采取基于单GPU(graphics processing unity)或多GPU的CUDA(computed unified device architecture)方法^[11-12]、基于集群机多节点的MPI(message passing interface)方法^[13]或混合并行方法^[14].
本文针对成像的效率问题, 从重力异常三维相关成像的基本原理出发, 提出基于泰勒级数的相关成像方法(probability tomography based on Taylor series, PTTS).采用立方体模型代替原点元模型, 通过泰勒级数重新定义几何函数; 相比原方法减少了计算量, 提高了成像速度; 提出的立方体模型在一定程度上保障了成像精度.
1 基于泰勒级数的相关成像方法原理在笛卡尔坐标系中, 设坐标系原点为o, xoy平面与z轴正交, 且z轴正方向为垂直向下; 空间中的观测面上分布着若干测点, 其中第i个观测点坐标为(x_i, y_i, z_i).在传统原理中, 测区地下空间被划分成一系列的点元, 则第j个点质量在观测点上引起的重力异常Δg_j为

(1)

式中:γ为万有引力常数; ρ_j, v_j和B_j分别是第j个点元的密度、体积和几何函数.
B_j的定义式为^[3]

(2)

因此, 重力异常的相关系数为^[3]

(3)

由于地下空间是按点元划分的, 若点元数量较少, 则近似认为地下空间是利用“球”模拟的, 而相邻球体间存在缝隙, 不能准确地模拟出地下构造分布情况.在数据解释中, 观测面通常为高度是常数的平面, 在反演中也常将地下介质分割成小立方体; 为了便于与进一步的反演结果作对照, 本文将地下空间也按照立方体进行等分.下面以单个立方体为例推导PTTS方法.
已知当前坐标系中, 任意立方体产生的重力异常Δg为

(4)

式中:(ξ₁, ξ₂), (η₁, η₂)和(ζ₁, ζ₂)分别表示立方体在x, y和z方向上的分布区间; (x, y, z)和(ξ, η, ζ)分别为观测点坐标和立方体内任意一点的坐标; ρ是立方体的密度.
利用式(4)不能直接进行重力异常的计算, 由文献[15-16]可知:

(5)

(6)

(7)

式中:(ξ_i, η_j, ζ_k)为立方体的对顶点坐标, 与(ξ₁, ξ₂)，(η₁, η₂)和(ζ₁, ζ₂)相对应.式(5)中的解常用于重力异常正演, 但其形式较为复杂.
式(4)中的被积函数为

(8)

对式(8)在立方体中心点(ξ₀, η₀, ζ₀)处进行泰勒级数展开, 舍去偏导阶数大于2的项; 将展开式代入式(4)中, 化简得

(9)

式中:Δx, Δy和Δz分别为立方体在x, y和z方向上的棱长; ΔV为立方体的体积; f_ξξ, f_ηη和f_ζζ为展开式中f在x, y和z方向上的二阶偏导数:

(10)

(11)

(12)

利用式(9)可计算立方体引起的重力异常, 类比式(1)和式(2), 可得式(9)的中括号部分:

(13)

由式(2)和式(13)可知:两种几何函数具有相同的形式, 即只包含三方向坐标差值平方的计算, 但式(13)相当于在式(2)基础上增加了含二次偏导数的项.
同理, 式(5)中累和的部分也可作为几何函数, 即

(14)

式中, B′为该算法中的几何函数, 称为PTF(probability tomography based on forwarding).显然, 在基于立方体模型的相关成像方法中, 使用式(13)的计算量要小于式(14)且前者的计算形式更简单.将式(13)或式(14)代入式(3)中可得到对应的相关系数.
多个立方体在各测点上的B或B′构成几何函数矩阵, 矩阵的行数等于测点个数, 列数等于立方体个数, 则第j个立方体在第i个测点上的几何函数位于矩阵的第i行第j列.
2 模型建立与验证2.1 三维地质建模在笛卡尔坐标系下建立理论模型来验证PTTS方法的效果:选定x和y方向上范围在0~2 000 m, 以及深度方向上范围在0~1 000 m的测区, 观测z=0平面, 测区内均匀分布20×20个测点, 地下空间被等分为20层.设地下存在密度为1 g/cm³, 大小为400 m×400 m×400 m, 200 m×1 200 m×400 m的两个三维地质模型, 分别模拟块状地质体和岩墙, 这两个模型为地球物理反演中常用的经典模型, 选择它们可更好地印证PTTS方法的效果, 选取y=1 000 m和z=350 m两个剖面表示模型位置, 见图 1.
图 1(Fig. 1)

图 1 理论模型位置Fig.1 Location of theoretical models (a)—y=1 000 m；(b)—z=350 m.

为验证方法的抗噪能力, 向模型数据中加入5 %的高斯噪声, 产生的重力异常见图 2.
图 2(Fig. 2)

图 2 重力异常的理论值Fig.2 Theoretical value of gravity anomaly (a)—不含噪声；(b)—含5 %高斯噪声.

2.2 数据试验与成像效果分析分别使用PTF和PTTS方法对模型数据进行成像并对比, 使用y=1 000 m, y=500 m和z=350 m三个剖面进行结果分析.在不含噪声的试验结果中(图 3), 分别对比图 3a和3d, 3b和3e, 3c和3f, 可见在不同位置的剖面上, 两种方法均能够对地下的地质体分布进行准确成像.其中, 由于y=1 000 m剖面位于两地质体中央, 故该剖面上能够同时体现出块体和岩墙模型, 由于岩墙能够引起更大的异常, 它所在位置的相关系数也具有较大的值; 而y=500 m剖面只通过岩墙且位于该模型的一端, 因此该剖面上仅能表示岩墙的存在, 结果中呈现的是高值区左“低”右“高”的形态; 在z=350 m剖面上可以清晰地反映出模型的水平位置, 同样地, 岩墙模型所在位置具有更大的相关系数值.
图 3(Fig. 3)

图 3 无噪声PTF和PTTS成像结果Fig.3 Tomography results of PTF and PTTS without noise (a), (d)—y=1 000 m；(b), (e)—y=500 m；(c), (f)—z=350 m.

为进一步分析该方法的成像能力, 含噪声情况下的成像结果见图 4.通过分别对比图 4a和4d, 4b和4e, 4c和4f判断PTF和PTTS的效果.
图 4(Fig. 4)

图 4 有噪声PTF和PTTS成像结果Fig.4 Tomography results of PTF and PTTS with noise (a), (d)—y=1 000 m；(b), (e)—y=500 m；(c), (f)—z=350 m.

由图 4可知PTF和PTTS均具有良好的抗噪性, 可较为清晰地表示块体和岩墙的位置.由于噪声对理论模型具有一定的影响, 图 3中结果的聚焦程度、清晰度等方面要略优于图 4, 但效果相差不大.需要注意的是, 上述结果不能准确地反映模型的底部, 这符合成像特点, 不属于方法的误差.
2.3 计算效率分析数据试验表明基于立方体计算的PTF与PTTS方法具有相近的成像能力, 但通过式(13)和式(14)能够证明二者计算效率不同.因此, 在第2.1节模型基础上另设测区内均匀分布40×40, 80×80个测点, 利用不同规模的观测数据分析效率.本文采用Matlab R2014a软件完成算法实现和计时, 在处理器为i7-6700HQ 2.6-3.5 GHz, 内存为16 GB DDR4, 操作系统为Windows 10, 型号为MSI-GE72-6QF-073XCN的笔记本上电脑运行程序, 时间记录见表 1.
表 1(Table 1)

表 1 PTF与PTTS运算时间Table 1 Calculation time of PTF and PTTS s 方法 20×20 40×40 80×80 PTF 3.712 54.567 881.306 PTTS 1.679 26.774 435.971

表 1 PTF与PTTS运算时间 Table 1 Calculation time of PTF and PTTS s

在表 1中, PTTS与PTF的运行时间之比分别为0.45, 0.49, 0.49, 说明PTTS的计算效率至少比PTF快一倍以上; 随着观测数据量增加, 运行时间比值也增加, 但趋于一稳定值.这是由于两种方法的区别主要在于几何函数矩阵的计算, 该部分是影响两方法运行时间的主要因素.PTF的几何函数利用三个循环来实现且计算中包含对数、三角函数等运算; 而PTTS方法仅包括四则运算, 计算量和难度都比PTF小.当观测数据量较少时, PTTS通过减少循环所体现出的计算效率较为明显; 当数据量增加时, 效率降低并趋于不变.因此, 通过新的几何函数能够在保证结果精度的基础上提升计算效率.
综上, 在解释效果相近的前提下, 采取PTTS方法能够取得更高的成像效率.
3 相关成像在文顿盐丘地区的应用将PTTS方法应用于文顿盐丘地区的重力异常数据.许多学者曾对该地区的实测航空重力全张量梯度数据进行过大量的研究^{[12-14, 17-20]}, 观测异常来源于地下盐丘中的盖岩^[17-19], 其剩余密度约为0.55 g/cm³.重力异常数据是由Bell Geospace公司在地面测量的, 位于WGS84坐标系中, 共20×20个测点均匀分布于地表; x, y轴正方向分别表示东、北, z轴正方向垂直向下.数据在x, y方向范围分别为441.55~443.55 km, 3 332.85~3 334.85 km, z方向最大深度为1 km并等分为20层.成像前通过空间波长为5 000 m的高通滤波器移除异常的背景场.
选取y=3 333.85 km, x=442.55 km和z=350 m三个剖面对PTF与PTTS的成像结果进行解释, 见图 5.
图 5(Fig. 5)

图 5 PTF和PTTS实测结果Fig.5 Testing results of PTF and PTTS (a)，(d)—y=3 333.85 km；(b)，(e)—x=442.55 km；(c)，(f)—z=350 m.

PTF和PTTS两种方法均能够反映盖岩的位置, 结果中相关系数较高的部分即盖岩位置, 其在x, y和z方向的分布范围分别为442.05~442.95 km, 3 333.65~3 334.75 km和150~650 m, 这与其他学者反演的结果相符^[19-20], 证明了PTTS方法能够应用于实测数据解释工作中.
4 结论本文针对重力异常数据提出了基于泰勒级数的相关成像法(PTTS), 该方法在原理论基础上采用三维建模, 通过泰勒级数展开式与积分等运算重新定义了几何函数的计算方法.在理论模型试验中验证了该方法对地下地质体分布具有较好的反映和良好的抗噪性, 在保障结果精度的同时能够有效提升计算效率约一倍以上; 将PTTS应用于文顿盐丘地区实测数据的计算处理, 证明其具有良好的效果.
致谢 感谢Bell Geospace Inc.提供的实测重力异常数据.
参考文献

[1]	Patella D. Introduction to ground surface self-potential tomography[J].Geophysical Prospecting, 1997, 45(4): 653–681.DOI:10.1046/j.1365-2478.1997.430277.x
[2]	Patella D. Self-potential global tomography including topographic effects[J].Geophysical Prospecting, 1997, 45(5): 843–863.DOI:10.1046/j.1365-2478.1997.570296.x
[3]	郭良辉, 孟小红, 石磊, 等. 重力和重力梯度数据三维相关成像[J].地球物理学报, 2009, 52(4): 1098–1106. ( Guo Liang-hui, Meng Xiao-hong, Shi Lei, et al. 3D correlation imaging for gravity and gravity gradiometry data[J].Chinese Journal of Geophysics, 2009, 52(4): 1098–1106.DOI:10.3969/j.issn.0001-5733.2009.04.027)
[4]	郭良辉, 孟小红, 石磊. 磁异常ΔT三维相关成像[J].地球物理学报, 2010, 53(2): 435–441. ( Guo Liang-hui, Meng Xiao-hong, Shi Lei. 3D correlation imaging for magnetic anomaly ΔT data[J].Chinese Journal of Geophysics, 2010, 53(2): 435–441.DOI:10.3969/j.issn.0001-5733.2010.02.022)
[5]	Guo L H, Meng X H, 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]	Guo L H, Shi L, Meng X H. 3D correlation imaging of magnetic total field anomaly and its vertical gradient[J].Journal of Geophysics and Engineering, 2011, 8(2): 287–293.DOI:10.1088/1742-2132/8/2/013
[7]	Mauriello P, Patella D. Resistivity anomaly imaging by probability tomography[J].Geophysical Prospecting, 1999, 47(3): 411–429.DOI:10.1046/j.1365-2478.1999.00137.x
[8]	Mauriello P, Patella D. Principles of probability tomography for natural-source electromagnetic induction fields[J].Geophysics, 1999, 64(5): 1403–1417.DOI:10.1190/1.1444645
[9]	孟小红, 刘国峰, 陈召曦, 等. 基于剩余异常相关成像的重磁物性反演方法[J].地球物理学报, 2012, 55(1): 304–309. ( Meng Xiao-hong, Liu Guo-feng, Chen Zhao-xi, et al. 3D gravity and magnetic inversion for physical properties based on residual anomaly correlation[J].Chinese Journal of Geophysics, 2012, 55(1): 304–309.DOI:10.6038/j.issn.0001-5733.2012.01.030)
[10]	Liu G F, Yan H, Meng X H, et al. An extension of gravity probability tomography imaging[J].Journal of Applied Geophysics, 2014, 102: 62–67.DOI:10.1016/j.jappgeo.2013.12.012
[11]	Chen Z X, Meng X H, Guo L H, et al. GICUDA:a parallel program for 3D correlation imaging of large scale gravity and gravity gradiometry data on graphics processing units with CUDA[J].Computers & Geosciences, 2012, 46: 119–128.
[12]	侯振隆.重力全张量梯度数据的并行反演算法研究及应用[D].长春: 吉林大学, 2016. ( Hou Zhen-long.Research and application of the parallel inversion algorithms based on the full tensor gravity gradiometry data [D].Changchun: Jilin University, 2016.)
[13]	Hou Z L, Huang D N, Wei X H. Fast inversion of probability tomography with gravity gradiometry data based on hybrid parallel programming[J].Journal of Applied Geophysics, 2016, 124: 27–28.DOI:10.1016/j.jappgeo.2015.11.009
[14]	Hou Z L, 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
[15]	Okabe M. Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies[J].Geophysics, 1979, 44(4): 730–741.DOI:10.1190/1.1440973
[16]	Steiner F, Zilahi-Sebess L. Interpretation of filtered gravity maps[M]. Budapest: Akademiai Kiado, 1988: 344-350.
[17]	Coker M O, Bhattacharya J P, Marfurt K J. Fracture patterns within mudstones on the flanks of a salt dome:syneresis or slumping?[J].Gulf Coast Association of Geological Societies, 2007, 57: 125–137.
[18]	Ennen C, Hall S.Structural mapping of the Vinton salt dome, Louisiana, using gravity gradiometry data[C]// SEG Annual Meeting.San Antonio, 2011: 830-835.
[19]	Oliveira V C Jr, Barbosa V C F. 3D radial gravity gradient inversion[J].Geophysical Journal International, 2013, 195(2): 883–902.DOI:10.1093/gji/ggt307
[20]	秦朋波, 黄大年. 重力和重力梯度数据联合聚焦反演方法[J].地球物理学报, 2016, 59(6): 2203–2224. ( Qin Peng-bo, Huang Da-nian. Integrated gravity and gravity gradient data focusing inversion[J].Chinese Journal of Geophysics, 2016, 59(6): 2203–2224.)