0 引言
重力勘探是大面积快速探测中最常用的地球物理探测技术之一，其主要目的是推断引起重力异常的地质结构的几何形态与密度异常。然而，和其他位场方法一样，重力场异常的地质解释也是非唯一的，因为多种重力密度与几何分布情况可以引起相同的观测异常[1]。虽然使用球体、圆柱或者棱镜等简单的几何形状剖分进行重力数值模拟的方法通常是合理的[2]，但是高分辨率的重力与磁法勘探的数值模拟也非常重要[3]。
随着地球物理勘探成为矿产普查与勘探的重要工具，新的技术被引入到二维观测中，并使用重磁异常曲线约束二维地质结构[4]。目前为止，出现了很多可以用位场数据建立二维或者三维地质模型的工具[5-10]。Plouff[11]建立了有限水平板重力异常的数值解，对地形的校正也有非常好的效果；Pignatelli等[12]用一组具有恒定磁化率和密度的棱柱形单元模拟地下密度与磁化强度的分布；Tontini等[13]基于三维快速傅里叶变换(FFT)推导了重力与磁场的三维正演模拟公式。为了利用重磁场数据研究地下结构，已经开发出各种软件。Blakely[14]将它们分为三类：正演数值模拟、反演计算和数据增强与显示方法。重力数值模拟方法基于地质和地球物理的判断，为源体构造初始模型，然后计算模型的重磁响应，以便与观测到的异常进行比较[15]。因此，地球物理正演模拟在重力异常的直接解释中起着重要作用。断层、褶皱等地质构造往往复杂，密度、磁化强度等地球物理参数也往往不均匀。张岭等[16]在2006年针对任意形状的二维模型进行了Delaunay的网格剖分，将模型剖分成了若干个非结构的三角形，并计算了模型的重力异常值，但是仅仅适用于二维的重力数值模拟。考虑到二维或2.5维建模具有一定的局限性，Barnett和Okabe[17-19]将均匀多面体引入了重力与磁法的数值模拟中。使用多边形单元构成真实的任意形状地质体是非常容易的，逼近真实程度取决于多边形面的数量和顶点的选择，对地形的校正也有非常好的效果[20]。
一个非均匀的地质体可以被划分成几个更小的同质地质体组合。因此，由均质多面体引起的重力异常的计算比起其他类型的方法，特别是对于高分辨率的重力测量，更具有实用性和重要性[21]。随着计算设施、软件技术和数值方法的快速发展，应用数学、地球物理数值模拟的能力得到大大增强[22]。由于地下介质具有复杂的几何结构，任何实际反演问题的求解都需要求解类似地下结构的正演问题。因此，很多数值方法被引入到地球物理领域，例如有限差分、有限元、积分式、体积积分、边界积分和混合方法等[23-24]。李振海等[25]在2012年采用三维Delaunay将三维目标地质体分解为若干四面体，建立了剩余密度值与重力异常值的线性式组，进行了数值模拟，但该计算方法在程序上实现较复杂，且模型相对简单，没有局部加密、刻画边界等特点。
将三维地质体分解成有限个四面体后，重力场等于所有四面体对地表产生的重力响应的总和。这种方法在复杂不规则区域建模中的优点包括：可以使用非均匀的网格计算，能精细刻画地质体，计算结果准确；可以根据模型需要，对模型的边界进行刻画、局部区域进行网格加密，而其他区域的网格密度可以根据其特点进行疏松处理，进而节省数值模拟的计算时间。
在实际应用中，往往要求生成的Delaunay三角网满足一些约束(限定)条件：在二维情况下有约束点、约束线，在三维情况下还包括约束面。在上述约束条件下生成Delaunay三角网格的过程称作约束Delaunay三角剖分，也称为限定Delaunay三角剖分。
本文利用约束Delaunay三角网格剖分的办法将三维复杂地质体进行四面体剖分，该方法实现了使用Delaunay三角网格对复杂地质体建模的目标；同时，提出并实现了使用重心体积的方法进行三维带地形的重力数值模拟；最后，用数值实例和野外实例进行了验证。
1 基于Delaunay网格的重力正演方法将地质体模型进行三维Delaunay三角网格剖分之前，将地质体模型离散成赋有地质体信息的空间分布的点的集合，可根据需要调整空间分布的点的密度。根据空间点集的位置信息，生成Delaunay三角网格组成的地质体模型[26]，然后进行三维重力正演计算(图 1)。
|
| P (x, y, z)是计算异常的观测点；E (ξ, η, ζ)是四面体的重心；A、B、C、O是四面体的4个顶点；Δh为目标体重心与地表测量点的垂直距离。 图 1 重力正演原理图 Fig. 1 Gravity forward diagram |
|
|
首先计算每一个四面体对地表的重力响应，并假设单元四面体内密度值是均匀的。将所有四面体计算的数值求和，就能得到该模型在地表任意点的重力响应。
单个四面体对地表任意一点重力响应的表达式为
(1) 式中：Δgq代表地下的地质体q在地表某一点重力异常值的垂直分量；G为有引力常数；Δρq为目标地质体q的密度；Vq为目标地质体体积；(x, y, z)为地表观测点坐标；(ξ, η, ζ)为四面体的重心点坐标。观测点可以位于地表起伏的界面上，公式(1)可以实现起伏地表的重力异常计算。
因此Q个地质体的地表重力响应可以表示为
(2) 根据式(2)可得：
(3) 式中，Δh=(z—ζ)。因此，求取式(3)每个单元的体积，就可以计算地质体地表重力响应。
四面体越接近等边四面体，Delaunay三角网格质量越高，数值模拟的精度也越高[26]。主要原因为当异常体与观测点距离较小时，若单元质量较差，会出现观测点与重心之间距离太近的情况，对Δg的计算结果会造成很大影响。当异常体与观测点距离较远时，四面体可以近似看作质点对观测点的重力异常响应，网格质量较高的四面体的性质更接近质点。
任意四面体体积y表示为
(4) 令O点坐标为(xO, yO, zO)，A点坐标为(xA, yA, zA)，B点坐标为(xB, yB, zB)，C点坐标为(xC, yC, zC)，可得
(5) 考虑一个简单的“十”字模型，如图 2a所示。该模型剩余密度为400 kg/m3，上表面高度为-600 m，下表面高度为-1 000 m，短边为800 m，长边为1 600 m。观测点位于高度-100~100 m的连续随机起伏界面上，且观测点为每个Delaunay三角网格的顶点。使用本文提出的重心体积方法计算重力响应之后，需要对计算结果进行检验，否则不能保证该方法获得重力场的可靠性。由于该模型为规则形体，用长方体网格剖分计算的数值解可以等效认为是解析解，如图 2b所示。同时对该模型分别采用不同密度的点集进行刻画，然后进行网格剖分，如图 2c、d所示。3种网格剖分的结果如表 1所示。如图 3所示，使用四面体网格剖分模型的正演结果与长方体计算结果进行对比可以看出，细化网格计算的重力场更接近解析解，粗糙网格的误差相对较大，但仍然满足数值模拟的要求。因此，在计算实际数据时，我们可以采用局部点集加密的方法提高重力场的准确性。这不仅可以满足对局部地质体的刻画，同时还能大量节省重力数值计算的时间。
|
| a.异常体属性；b.长方体网格；c. Delaunay粗糙网格；d. Delaunay细化网格。 图 2 异常体模型不同网格剖分结果 Fig. 2 Results of different mesh generation of abnormal body models |
|
|
| 剖分方法 | 网格数 | 内角最小值/(°) | 内角最大值/(°) | 均方差 | CPU计算时间/s |
| 粗糙网格 | 721 | 72.649 542 | 132.863 07 | 0.501 506 | 0.703 125 |
| 细化网格 | 8 736 | 72.558 988 | 138.892 73 | 0.127 238 | 4.781 250 |
| 长方体 | 93 | 0.012 5 |
|
| a.长方体解析解；b.Delaunay粗糙网格解析解；c.Delaunay细化网格解析解；d. a与b的差；e. a与c的差。 图 3 异常体模型不同网格剖分后的计算结果 Fig. 3 Calculation results of abnormal body by different mesh generation |
|
|
通过图 4与表 1可以看出：细化网格的四面体数相比粗糙网格多出10余倍，数值模拟时间高出几倍。模型被粗糙网格离散后，其内角范围在72°~133°之间，被细化网格离散后，内角范围在72°~139°之间。其中，细化网格离散模型后内角最大值较大，原因是细化网格刻画模型后会产生更多的四面体，产生质量较差的四面体的机会较大。然而，从直方图图 4可以看出，细化网格模型相比粗糙网格模型具有更大比例的小角度内角四面体。对于本文提出的重力数值模拟的方法，四面体内角角度越小，四面体越接近等边四面体，数值计算结果也越准确。所以，细化网格离散的模型虽然存在一些内角角度较大的四面体，但是相比粗糙网格离散的模型仍然具有更好的刻画精度。因此，适当增加四面体网格的密度，对重力数值模拟具有重要意义。
|
| a. Delaunay细化网格；b. Delaunay粗糙网格。 图 4 四面体最大内角角度分布直方图 Fig. 4 Histogram of maximum inner angle of tetrahedron |
|
|
金川含矿超基性侵入岩体位于西南边缘龙首山隆起带的东北端，矿区位于纬度38°N与39°N、经度102°E与103°E之间。龙首山隆起带北部的深断层(F1)与阿拉善地块的潮汐断陷盆地相邻，而南部的深断层则与北祁连早古生代褶皱带的走廊过渡带相邻。金川矿区自西向东被称为Ⅲ、Ⅰ、Ⅱ、Ⅳ矿区(图 5a), 各矿区岩体的规模、形态、产状及矿石存储量都不同，本文选择金川Ⅱ号矿区开展研究[27-28]。
|
| a.金川矿区平面图；b.金川矿区34号断面图；c. Delaunay网格剖分断面；d.金川Ⅱ号矿区1号矿体模型；e. Delaunay网格剖分模型。 图 5 研究区域图及地质模型 Fig. 5 Study area map and geologic model |
|
|
Ⅱ矿区侵入体大约长3 000 m，介于断层F16西南方向，东部被第四系覆盖，其余地段均出露地表[27]。由东向西逐渐变宽，在F17附近宽度最大，约530 m，再向西则又变窄。该侵入体走向约N50°W，倾角50°~60°，东部倾角相对变缓。西端受断层影响有明显的偏转，延深较大，最大延深超过1 000 m，其中赋存规模巨大的富矿体，即为1号矿体。东段岩体较浅, 出露较宽, 最大宽度可达530 m，即为2号矿体。1和2号两个主矿体含镍量约占本矿区的99%以上。图 5a中紫线为图 5b断面位置，并被离散为图 5c。
首先对金川铜镍矿的Ⅱ号局部矿区进行重力建模，如图 5d所示。根据断面图所示，矿化带及贫富矿体被埋深在围岩中，该断面由钻井控制，其准确性可以保证。同时，该地区存在其他断面及钻井信息，虽然其他大部分断面矿体及矿化带延伸到钻井控制范围以外，缺少部分数据信息，但仍然可以作为建模参考标准。因此采用该断面建模具有一定的科学性。建模后将该模型使用Delaunay网格进行网格剖分。观测点位于起伏的地表，每个观测点为Delaunay四面体网格位于起伏地表上的顶点，如图 5b与5c所示。为了更好地刻画岩矿体边界及岩矿体内部精细构造，增加相应位置顶点数，能够得到网格密度不均匀的重力模型。如图 5e所示，根据已知资料建立金川Ⅱ矿区重力模型：绿色区域密度约为2 970 kg/m3，主要成分为基性橄榄岩；红色区域是高品质矿体，主要成分为海绵陨铁状硫化物纯榄岩，密度为3 030 kg/m3左右；黄色区域为低品质矿体，主要成分为超基性岩型硫化物表外矿石，密度为3 030 kg/m3；灰色区域是围岩，主要包括大理岩、混合岩等，密度为2 670 kg/m3。将2 670 kg/m3视为背景场，这些密度参数主要由野外采集试验获得。
矿体模型中四面体的剩余密度为金川矿区地质体密度与围岩密度的差。因此，高品质与低品质矿体的剩余密度为360 kg/m3，橄榄岩的剩余密度为300 kg/m3。最后使用本文提出的算法进行正演模拟，可以得到剩余密度对地表产生的重力异常响应(图 6)。
同时我们对实际数据与数值模拟结果进行对比(图 7)。可以看出，使用本文算法得到的计算结果与真实数据比较接近，效果良好，真实数据与模拟结果最大异常数值均接近2 000×10-8 m/s2。由于金川Ⅱ号矿区仅为金川矿区一部分，因此实际数据曲线与数值模拟有些偏差，或是该位置深部仍然存在未发现的矿体。
|
| 图 7 实际数据与计算结果 Fig. 7 Measured data and calculation results |
|
|
本文在重力数值模拟的过程中引入了Delaunay网格生成技术, 利用四面体的重心与体积数据计算获得了地表重力异常场。通过模型试验及计算结果分析得出了以下结论：
1) Delaunay网格是离散目标地质体的一种高效的网格生成算法。这种算法不仅考虑到目标体的复杂构造，还能够进行地质模型的网格分辨率控制。
2) 在Delaunay网格的生成过程中可以根据需要对网格进行优化控制。合成数据结果显示，细化网格刻画的模型有较大比例的高质量四面体，为数值模拟的精度提供了保证。
3) 在实际应用中，使用Delaunay网格剖分方法能够有效地剖分复杂模型，同时，数值模拟计算出的异常曲线与实际异常曲线有较好的一致性。
| [1] |
Juan G A. 3D Forward and Inverse Modeling of Total-Field Magnetic Anomalies Caused by a Uniformly Magnetized Layer Defined by a Linear Combination of 2D GaussianFunctions[J]. Geophysics, 2007, 73(1): 11. |
| [2] |
Bhattacharyya B K. A Generalized Multibody Model for Inversion of Magnetic Anomalies[J]. Geophysics, 1980, 45(2): 255-270. DOI:10.1190/1.1441081 |
| [3] |
Gunn P J. Quantitative Methods for Interpreting Aeromagnetic Data:A Subjective Review[J]. AGSO Journal of Australian Geology & Geophysics, 1997, 17(2): 105-113. |
| [4] |
Nettleton L L. Gravity and Magnetic Calculations[J]. Geophysics, 1942, 7(3): 293. DOI:10.1190/1.1445015 |
| [5] |
Caratori T F, Cocchi L, Carmisciano C. Rapid 3-D Forward Model of Potential Fields with Application to the Palinuro Seamount Magnetic Anomaly (Southern Tyrrhenian Sea, Italy)[J]. Journal of Geophysical Research, 2009, 114(B2): B02103. |
| [6] |
Hamilton D E, Jones T A. Computer Modeling of Geologic Surfaces and Volumes[C]//AAPG Computer Applications in Geology. Tulsa: [s.n.], 1992: 297.
|
| [7] |
Christian J T. 3D Geoscience Modeling:Computer Techniques for Geological Characterization[M]. Berlin: Berlin Springer-Verlag, 1996.
|
| [8] |
Jessell M. Three-Dimensional Geological Modelling of Potential-Field Data[J]. Computers & Geosciences, 2001, 27(4): 455-465. |
| [9] |
Reinhard P, John W H. Computer Graphics in Geology[J]. Lecture Notes in Earth Sciences, 1992, 41(4): 613-614. |
| [10] |
Talwani M. Computation with the Help of a Digital Computer of Magnetic Anomalies Caused by Bodies of Arbitrary Shape[J]. Geophysics, 1965, 30(5): 797. DOI:10.1190/1.1439654 |
| [11] |
Plouff D. Gravity and Magnetic Fields of Polygonal Prisms and Applications to Magnetic Terrain Corrections[J]. Geophysics, 1976, 41: 727-741. DOI:10.1190/1.1440645 |
| [12] |
Pignatelli A, Nicolosi I, Carluccio R, et al. Graphical Interactive Generation of Gravity and Magnetic Fields[J]. Computers & Geosciences, 2011, 37(4): 567-572. |
| [13] |
Tontini C. Rapid Interactive Modeling of 3D Magnetic Anomalies[J]. Computers & Geosciences, 2012, 48: 308-315. |
| [14] |
Blakely R J. Potential Theory in Gravity and Magnetic Applications[M]. London: Cambridge University Press, 1995.
|
| [15] |
Shin Y H, Choi K S, Xu H. Three-Dimensional Forward and Inverse Models for Gravity Fields Based on the Fast Fourier Transform[J]. Computers & Geosciences, 2006, 32(6): 727-738. |
| [16] |
张岭, 郝天珧. 基于Delaunay剖分的二维非规则重力建模及重力计算[J]. 地球物理学报, 2006, 49(3): 877-884. Zhang Ling, Hao Tianyao. 2-D Irregular Gravity Modeling and Computation of Gravity Based on Delaunay Triangulation[J]. Chinese Journal of Geophysics, 2006, 49(3): 877-884. DOI:10.3321/j.issn:0001-5733.2006.03.033 |
| [17] |
Barnett C T. Theoretical Modeling of the Magnetic and Gravitational Fields of an Arbitrary Shaped Three-Dimensional Body[J]. Geophysics, 1976, 41: 1353-1364. DOI:10.1190/1.1440685 |
| [18] |
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 |
| [19] |
Pohanka V. Optimum Expression for Computation of the Gravity Field of Homogeneous Polyhedral Body[J]. Geophysical Prospecting, 1988, 36(7): 733-751. DOI:10.1111/j.1365-2478.1988.tb02190.x |
| [20] |
Liu S, Hu X, Xi Y, et al. 2D Inverse Modeling for Potential Fields on Rugged Observation Surface Using Constrained Delaunay Triangulation[J]. Computers & Geosciences, 2015, 76: 18-30. |
| [21] |
Luo Yao, Yao Changli. Forward Modeling of Gravity, Gravity Gradients, and Magnetic Anomalies Due to Complex Bodies[J]. Journal of China University of Geosciences, 2007, 18(3): 280-286. DOI:10.1016/S1002-0705(08)60008-4 |
| [22] |
Roy K K. Potential Theory in Applied Geophysics[M]. Calcutta: Springer Science & Business Media, 2007.
|
| [23] |
刘海飞, 柳建新, 郭荣文, 等. 起伏地形三维激电连续介质模型快速反演[J]. 吉林大学学报(地球科学版), 2011, 41(4): 1212-1218. Liu Haifei, Liu Jianxin, Guo Rongwen, et al. Efficient Inversion of 3D IP Data for Continuous Model with Complex Geometry[J]. Journal of Jilin University (Earth Science Edition), 2011, 41(4): 1212-1218. |
| [24] |
林家勇, 汤井田, 丁茂斌, 等. 复杂地形条件下激发极化有限单元法三维数值模拟[J]. 吉林大学学报(地球科学版), 2010, 40(5): 1183-1187. Lin Jiayong, Tang Jingtian, Ding Maobin, et al. Three-Dimension Numerical Simulation of Induced Polarization & Finite Element Method Under Complicated Terrain[J]. Journal of Jilin University (Earth Science Edition), 2010, 40(5): 1183-1187. |
| [25] |
李振海, 罗志才, 钟波. 基于3D Delaunay剖分算法的重力建模与分析[J]. 地球物理学报, 2012, 55(7): 2259-2267. Li Zhenhai, Luo Zhicai, Zhong Bo. Gravity Modeling and Analyzing Based on 3D Delaunay Triangulation Algorithm[J]. Chinese Journal of Geophysics, 2012, 55(7): 2259-2267. |
| [26] |
郑耀, 陈建军. 非结构网格生成:理论、算法和应用[M]. 北京: 科学出版社, 2016. Zheng Yao, Chen Jianjun. Unstructured Mesh Generation:Theories, Algorithms and Applications[M]. Beijing: Science Press, 2016. |
| [27] |
Tang Z. Main Genetic Types of Ni Ore Deposits in China and Their Relations to Paleo-Plate Tectonics[J]. Geochemistry, 1984, 3(2): 102-114. |
| [28] |
段俊, 钱壮志, 焦建刚, 等. 甘肃龙首山岩带西井镁铁质岩体成因及其构造意义[J]. 吉林大学学报(地球科学版), 2015, 45(3): 832-846. Duan Jun, Qian Zhuangzhi, Jiao Jiangang, et al. Genesis of Xijing Intrusion from Longshoushan Terrane and the Tectonic Significance[J]. Journal of Jilin University (Earth Science Edition), 2015, 45(3): 832-846. |