2. 武汉大学地球空间环境与大地测量教育部重点实验室, 武汉 430079;
3. 武汉大学测绘遥感信息工程国家重点实验室, 武汉 430079
2. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China;
3. State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, Wuhan 430079, China
地球重力场信息在地球构造研究、国土资源调查、矿产勘探以及国民经济等领域中发挥着重大作用.重力的正问题与反问题是解释重力资料的核心研究内容,其中前者主要研究不同形状、产状和场源密度等场源体或地质体所引起的重力异常特征及分布等,其计算过程为正演计算;后者主要根据重力异常的分布来计算地质体或场源体的密度分布,其计算过程为反演计算[1].目标地质体的重力建模建立了地质体的剩余密度值与地面重力异常值的关系,既可利用重力正演公式计算地质体所引起的重力异常,也可以应用于求解重力反问题.因此,重力建模是重力正反演问题的关键环节,建模方法的优劣则直接影响到重力正反演的计算精度.
随着对重磁资料解释的深入研究、重力测量精度的提高及计算机的快速发展,重力正反演研究由过去的二维逐步发展为三维.对于三维连续密度变化的目标体,目前最常用的建模方法是块体方法,即将目标地质体剖分为规则的构造单元如长方体,计算所有单元产生的异常并叠加便得到目标地质体所产生的异常值.楼海等[2]给出了基于矩形网格模型的三维连续密度分布目标地质体的重力解析公式,并在实际应用中得到较好的结果,但当地质形状比较复杂时,其计算结果并不理想,如果加密构造单元,则又会增加计算量及存储空间.盛国平[3]采用水平截面法对非规则几何形状的目标体进行分割,并假设每一个质面层的密度值为常数,给出了较严密的数学解析表达式.张岭等[4]针对截面为任意形状的二度体问题,利用二维Delaunay剖分方法,将截面分割成若干三角形,二度体被分解成若干三棱柱的组合来计算二度体的重力异常值,这种方法表达了二维的连续密度变化,但仅限制于二度体的重力计算问题.
对于任意形状和变密度的三度体,现有重力建模方法均存在不足,为此,本文提出基于3D Delaunay 剖分算法的重力建模方法.采用3D Delaunay剖分算法将三维目标地质体分解为若干变密度四面体体元,推导基于四面体体元的重力正演公式;比较分析常规块体方法和本文方法应用于重力正演的计算效果,并采用共轭梯度法加密度约束条件对非规则形状变密度的地质体进行反演计算,验证本文重力建模方法的正确性和有效性.
2 3D Delaunay剖分的重力建模 2.1 地质体的3D建模地质体的3D 建模步骤为:首先根据地质体的三维坐标信息及物性信息等相关资料将目标地质体离散化为空间离散点;其次,采用合理的建模算法组织空间离散点,比如传统重力建模时采用的块体方法,地学领域2D、2.5D 及3D Delaunay 剖分方法等;最后将所有体元结构(例如四面体、三棱柱体、六面体和长方体等)组合起来,根据体元顶点属性特征采用颜色或纹理对模型进行可视化.对于规则形状的地质体建模,采用传统块体方法即可对地质体进行快速建模,并真实逼近模型结构及地质体的物性特征,但是对于非规则形状变密度地质体的建模,如褶皱、断层等,传统建模方法则会改变地质体的结构、密度分布和其他物性特征,所以本文引进3D Delaunay剖分方法来逼近真实的地质体.
3D Delaunay剖分算法由二维Delaunay三角剖分算法演化而来,1908年G.Voronoi为了限定二维平面离散点的有效作用范围,首先定义了二维平面上的Voronoi图,1934 年由B.Delaunay将Voronoi图演化出了更易于分析应用的Delaunay 三角网. Voronoi图和Delaunay三角网为目前普遍接受和广泛采用的分析研究区域离散数据的有力工具,已应用到石油勘探、地质、矿业、城市规划和环境监测等领域[5].Lewis[6]实现了3D Delaunay网生成算法,Chen[7] 对边界恢复算法进行改进,陈晓勇[8]、李清泉[9]等研究了四面体格网结构(TEN)模型的生成算法.
3D Delaunay剖分算法的剖分结果为四面体体元,四个顶点为空间离散点,包含坐标信息和物性特征.四面体格网数据模型是二维三角形网(Triangulated Irregular Network,TIN)数据结构在三维空间上的扩展,与之相应的数据结构为四面体格网结构(TEN).TEN 模型以3D Delaunay三角剖分为基础,将目标空间用紧密排列但不重叠的非规则四面体的组合来表示.3D Delaunay 四面体应具有如下特点[10-11]:①生成的四面体之间无重叠部分; ②所形成的四面体组合可以覆盖整个目标体的三维空间;③四面体的外接球不包含空间离散点集中除该四面体四个顶点外的任一点.这些性质保证了组成四面体的三角形近似等边或等角,及四面体体元的组合更逼近真实的目标地质体.
为了比较常规块体方法及3D Delaunay剖分算法对重力建模的影响,本文设计了两个目标地质体. 目标体Ⅰ为连续密度变化的矩形体,上顶面埋深为1km,范围为10km×10km×1km(长×宽×高),八个角点的剩余密度值依次为100、300、300、500、 0、200、200和400kg/m3,如图 1a所示;目标体Ⅱ为倾斜台阶组合模型体,八个角点的剩余密度值与目标体Ⅰ相同,具体参数参见图 1b.
对于目标体Ⅰ 的建模,首先将目标体离散化为若干空间点位,再分别采用常规块体剖分方法与四面体剖分方法对这些离散点位建模.在x、y和z方向取分块数为3×3×3,空间点数为4×4×4,对目标体Ⅰ的剖分结果见图 2.比较图 2 和图 1a可知: 两种方法的建模结果均能准确表达地质体,理论上得到的重力正演计算结果应该相差不大.
采用传统块体方法对目标体Ⅱ的建模结果见图 3a.当采用3D Delaunay剖分算法对目标体Ⅱ 建模时,如果直接利用空间离散点进行3D Delaunay四面体构网,则会出现“跨越"现象,即左右倾斜台阶会出现错位TEN 结构,所以本文首先采用2D Delaunay算法建立左右倾斜台阶接触面的TIN 数据结构,作为组合模型的特征面,然后进行整体3D Delaunay四面体构网,其建模结果如图 3b所示.从图 3b可以看出,具有约束的3D Delaunay算法可以避免“跨越"现象,从而得到合理的三维Delaunay模型.比较图 3和图 1b可知:块体方法的建模结果与真实地质体有一定的差异,而3D Delaunay剖分算法的建模结果能更好地逼近真实地质体,因此理论上两种方法的重力正演计算结果应该存在较大差异;块体方法的建模结果与分块数密切相关,而3D Delaunay剖分算法受分块数的影响较小,随着图 3a 分块数的增加,块体方法的建模结果将逐渐逼近目标体的形状及密度分布,但计算量也将大幅度增加.
四面体体元为3D Delaunay剖分算法的基本结构单元,非规则形状变密度地质体对地面点产生的重力异常为所有四面体体元产生的重力异常之和. 任意一个四面体体元(见图 4)对地面点P产生的重力异常为
(1) |
式中,ξ =x-xP,η =y-yP,ζ =z-zP,ρ表示在四面体内任意点Q(x,y,z)的剩余密度,由四个顶点的剩余密度值线性插值得到,满足关系式:
(2) |
式中,λ1、λ2、λ3 和λ4 表示四面体体元四个顶点的剩余密度值.
对于任意形状的四面体,利用式(1)计算地面点的重力异常值,不易得到x、y和z三个方向的积分域,因此也不利于重力异常的计算.本文借鉴文献[3]提出的水平截面法,首先将四面体切分为若干个水平面层,计算每一层所产生的重力异常,然后叠加得到一个四面体所引起的重力异常值.与文献[3]的不同之处在于:四面体四个顶点的密度是不同的,其每一个面层也是变密度的.图 5 为变密度的水平截面被划分为若干个三角形的示意图.
如图 5所示,P′ 为地面计算点P在目标地质体的一个面层上的投影,则面层ABC被分为n个三角形,以ΔP′BC为例,计算一个三角形对P点产生的重力异常为
(3) |
G为万有引力常数,Z=Z0 -ZP,Z0 为截面埋深,σ 为面密度,ξ 和η 与式(1)相同.
将式(2)代入式(3),得到:
(4) |
将式(4)展开并整理得到:
(5) |
式中[I]为系数阵;[C]由四个顶点坐标组成;R=
对于I1 与I4 的计算,只需将直角坐标转换为极坐标,推导并整理为如下形式:
(6) |
(7) |
式中,W=+1 (mi<0),W=-1 (mi>0),S=+1 (pi<0),S=-1 (pi>0),
对于I2 与I3 的计算,可以利用格林公式将二重积分化为三角形边界的线积分,计算公式如下:
(8) |
(9) |
式中,
在计算系数[I]时,需要注意以下三点:
(1) 当线段平行于x轴(η =A0)时,I2 与I3 的组成项变为
当线段平行于y轴(ξ = B0)时,I2 与I3 的组成项变为
(2) P′、B和C三点在一条直线上,即mi=0,则ΔgP′BC=0;
(3) 根据格林公式求解I2 与I3 时,定义逆时针为正,要求I1、I4 与I2、I3 一致,即mi>0或pi>0 时表示顺时针,系数W=-1,S=-1.
采用公式(6)、(7)、(8)和(9)计算出系数阵[I] 后,即可计算一个三角形质面对地面P点产生的重力异常,按逆时针方向计算组成一个截面的n个三角形对P点产生的异常,将其叠加即为水平截面对地面P点产生的重力异常.假设四面体被分为了m层水平截面,vi、vi+1 和vi+2 分别为相应埋深Zi、Zi+1 和Zi+2 的水平截面所产生的重力异常值,则一个四面体对地面P点产生的重力异常为
(10) |
式中,
综上所述,最终可建立目标体M个剩余密度值[m]与N个地面重力异常[Δg]的关系式,即
(11) |
分别采用常规块体剖分方法和3D Delaunay剖分方法对目标体Ⅰ建模,并计算了位于目标体Ⅰ 正上方范围为10km×10km×3km 的空间格网点的重力异常,常规块体方法的正演计算公式可以参见文献[2],两种方法的计算结果如图 6所示.从图 6c 可以看出:两种方法的计算结果相差较小,最大差值为4.12×10-7 m/s2.
对目标体Ⅰ设计分块数为2-6时,分别采用常规块体剖分方法和3D Delaunay 剖分方法对其建模,正演计算了位于目标体Ⅰ 正上方1km 处100 个格网点的重力异常,其差值的比较结果如表 1 所示.从表 1可以看出:随着分块数的增加,两种方法的计算结果并无明显差异,最大差值在0.04×10-5m/s2 以内,均方误差在0.029×10-5m/s2 以内.因此,对于规则形状且连续密度变化的目标体,采用常规块体剖分方法和3D Delaunay剖分方法均能获得较好的重力正演结果.
分别采用常规块体剖分方法和3D Delaunay剖分方法对目标体Ⅱ建模,计算了与3.1 节相同区域的重力异常,计算结果如图 7所示.从图 7c可以看出:对于非规则形状变密度的目标体Ⅱ,两种方法的计算结果差异较大,其差值在(0.12~9.36)×10-6m/s2 范围之内,主要差异体现在双倾斜台阶组合模型体的公共斜平面上,其原因是常规块体方法改变了目标体的几何结构和物性特征,对地质体的建模不够准确.
类似于目标体Ⅰ,对目标体Ⅱ 设计分块数为2-6,分别采用常规块体方法和3D Delaunay剖分方法对目标体建模并计算其产生的重力异常,两种方法的比较结果见表 2.结果表明:随着分块数的增加,常规块体剖分方法对目标体Ⅱ 的建模逐渐逼近真实地质体,使得两者计算结果的差值逐渐减小.
从表 1和表 2 可以看出:平均值与均方误差量级相当,说明常规块体方法和3D Delaunay剖分方法之间存在明显的系统误差,并且系统误差随着分块数的增加而减小.这主要是目标体的建模误差和式(10)中划分层数m的取值所导致.根据文献[1] 和[3],m的取值应根据研究区域及其地质条件来确定,本文m的取值为1000.
3.3 目标体Ⅱ的反演计算3D Delaunay剖分算法不仅可以应用于地质体的重力正演计算,也可以应用于地质体的三维物性反演.针对地球物理的反演问题,Buckus和Gilbert[12-13] 最早对连续线性问题做了深入研究;Li和Oldenburg [14-15] 引入深度加权函数,克服了三维反演结果的“趋附" 现象;Braile等[16]将光滑度引入到正则化反演;在重磁反演求解大型线性方程组时,为了不降低计算效率及消耗过多计算机资源,Pilkington[17-18]采用共轭梯度法进行三维反演.本文基于3D Delaunay剖分算法,采用共轭梯度法[19-23]加密度约束条件对目标体Ⅱ进行三维物性反演.
地质体剩余密度值与地面重力异常值之间的线性关系式(11)往往是欠定的,不易得到真实物性反演结果,为此本文增加了密度约束条件Φm.定义目标函数的表达式为
(12) |
式中Φm为密度约束条件;Φd对应于式(11),其具体表达式如下:
(13) |
这里d为观测值,m为实例模型结点剩余密度值,G为该线性方程组的系数阵,d、m与G分别对应于式(11)的Δg、m与I′,ma和mb分别表示密度约束的下界和上界,在本算例的数值分别为100kg/m3、 500kg/m3.
采用最小二乘法求解式(13)时,其矩阵求逆时间较长,而且有可能出现病态方程,因此本文利用共轭梯度法求解.假定该区域的地质条件已知,背景密度为0kg/m3,以3.2 节基于本文方法正演计算的重力异常作为地面观测值,总共86个地面重力异常值,其等值线图见图 8a;采用86个地面重力异常值反演目标体Ⅱ的128 个剩余密度值,其反演结果见图 8d;利用密度反演结果计算得到的重力异常见图 8b;图 8c为图 8a和图 8b的重力异常之差,用于评价三维密度反演效果.
图 8d的结果表明,采用共轭梯度法加密度约束条件得到的反演结果非常接近于真实目标地质体Ⅱ 的密度分布;从图 8c可以看出:利用反演密度值计算的重力异常(见图 8b)与模拟重力异常(见图 8a) 相差很小,其量级为10-13m/s2.以上结果初步验证了本文的重力建模方法和三维密度反演方法的正确性和有效性.
4 结论与展望针对任意形状和变密度的三度体,本文提出了基于3D Delaunay剖分算法的重力建模方法.采用3D Delaunay剖分算法将三维目标地质体分解为若干变密度四面体体元,推导了基于四面体体元的重力正演公式;以规则形状且连续密度变化的地质体(如长方体)和非规则形状变密度的地质体(如倾斜台阶组合模型体)为例,比较分析了常规块体方法和3D Delaunay剖分算法应用于重力正演的计算效果,并采用共轭梯度法加密度约束条件对非规则形状变密度的地质体进行了反演计算.模拟结果表明: 对于规则形状且连续密度变化的地质体,基于常规块体方法和3D Delaunay剖分算法的重力建模均能获得满意的计算结果;对于非规则形状变密度的地质体,基于3D Delaunay剖分算法的重力建模能获得更加合理的计算结果;随着对地质体分块数的增加,尽管常规块体方法的建模结果将逐渐逼近真实地质体,但其消耗的计算机资源也随之大幅度增加. 计算结果初步验证了基于3D Delaunay剖分算法的重力建模更适合于存在褶皱、断层、裂缝等复杂地质体的重力正反演计算.
本文仅以长方体和倾斜台阶组合模型体为例验证了基于3D Delaunay剖分算法的重力建模方法的正确性和有效性,后续工作需要利用实际观测资料对该方法进行深入研究.实际应用本文方法时需要结合如地质物性数据、地震数据和地磁数据等先验信息对非规则变密度地质体进行重力建模,同时也需要进一步研究约束条件与重力正反演计算的最优化结合问题、数据离散化问题及三维物性反演的全局优化反演技术方法等.
[1] | 曾华霖. 重力场与重力勘探. 北京: 地质出版社, 2005 . Zeng H L. Gravity Field and Gravitational Prospecting (in Chinese). Beijing: Geological Publishing House, 2005 . |
[2] | 楼海, 王椿镛. 三维连续密度分布的重力计算及应用. 地震学报 , 1999, 21(3): 297–304. Lou H, Wang C Y. Gravity effect calculation of three-dimensional linear density distribution and its application. Acta Seismologica Sinica (in Chinese) (in Chinese) , 1999, 21(3): 297-304. |
[3] | 盛国平. 利用水平截面法对任意形状三度体的重力异常的快速正演. 青海地质 , 1993(1): 40–46. Sheng G P. Rapid direct method of gravity anomaly of arbitrary three-dimensional body by horizontal section. Geology of Qinghai (in Chinese) (in Chinese) , 1993(1): 40-46. |
[4] | 张岭, 郝天珧. 基于Delaunay剖分的二维非规则重力建模及重力计算. 地球物理学报 , 2006, 49(3): 877–884. Zhang L, Hao T Y. 2D irregular gravity modeling and computation of gravity based on Delaunay triangulation. Chinese J. Geophys. (in Chinese) , 2006, 49(3): 877-884. |
[5] | 武晓波, 王世新, 肖春生. Delaunay三角网的生成算法研究. 测绘学报 , 1999, 28(1): 28–34. Wu X B, Wang S X, Xiao C S. A new study of Delaunay triangulation creation. Acta Geodaetica Et Cartographica Sinica (in Chinese) (in Chinese) , 1999, 28(1): 28-34. |
[6] | Lewis R W, Zheng Y, Gethin D T. Three-dimensional unstructured mesh generation: part 3. volume meshes. Computer Methods in Applied Mechanics and Engineering , 1996, 134(3-4): 285-310. DOI:10.1016/0045-7825(95)00918-3 |
[7] | Chen J J, Zheng Y. Redesign of a conformal boundary recovery algorithm for 3D Delaunay triangulation. Journal of Zhejiang University-Science A , 2006, 7(12): 2031-2042. DOI:10.1631/jzus.2006.A2031 |
[8] | Chen X, Murai S. Analysis and visualization of 3D geospatial data by using Delaunay tetrahedral tessellation. Proceedings of Asian Conference on Remote Sensing (ACRS). 1999 . |
[9] | Li Q Q. Algorithms for tetrahedral network (TEN) generation. Geo-spatial Information Science , 2000, 3(1): 11-16. DOI:10.1007/BF02826801 |
[10] | 苏幸, 黄临平. 基于不规则四面体的三维离散数据地质建模算法. 物探与化探 , 2008, 32(2): 192–195. Su X, Huang L P. A three-dimensional dispersed data model based on the irregular tetrahedron and its generation algorithm. Geophysical and Geochemical Exploration (in Chinese) (in Chinese) , 2008, 32(2): 192-195. |
[11] | 吴江斌, 朱合华. 基于Delaunay构网的地层3D TEN模型及建模. 岩石力学与工程学报 , 2005, 24(24): 4581–4587. Wu J B, Zhu H H. 3D TEN model of strata and its realization based on Delaunay triangulation. Chinese Journal of Rock Mechanics and Engineering (in Chinese) (in Chinese) , 2005, 24(24): 4581-4587. |
[12] | Backus G E, Gilbert F. Numerical applications of a formalism for geophysical inverse problems. Geophys. J. R. Astr. Soc. , 1967, 13(1-3): 247-276. DOI:10.1111/gji.1967.13.issue-1-3 |
[13] | Backus G E, Gilbert F. The resolving power of gross earth data. Geophys. J. R. Astr. Soc. , 1968, 16(2): 169-205. DOI:10.1111/gji.1968.16.issue-2 |
[14] | Li Y, Oldenburg D W. 3-D inversion of magnetic data. Geophysics , 1996, 61(2): 394-408. DOI:10.1190/1.1443968 |
[15] | Li Y, Oldenburg D W. 3-D inversion of gravity data. Geophysics , 1998, 63(1): 109-119. DOI:10.1190/1.1444302 |
[16] | Braile L W, Keller G R, Peeples W J. Inversion of gravity data for two-dimensional density distributions. J. Geophys. Res. , 1974, 79(14): 2017-2021. DOI:10.1029/JB079i014p02017 |
[17] | Pilkington M. 3-D magnetic imaging using conjugate gradients. Geophysics , 1997, 62(4): 1132-1142. DOI:10.1190/1.1444214 |
[18] | Pilkington M. 3D magnetic data-space inversion with sparseness constraints. Geophysics , 2009, 74(1): L7-L15. DOI:10.1190/1.3026538 |
[19] | 陈召曦. 三维盐丘模型重力异常正演与约束反演研究[硕士论文]. 北京: 中国地质大学, 2009 . Chen Z X. The research on the forward calculation and constrained inversion of the 3D salt dome model’s gravity anomaly (in Chinese). Beijing: China University of Geosciences, 2009 . |
[20] | 李军. 基于Marquardt算法重磁复杂形体三维全局反演研究. 成都: 成都理工大学, 2005 . Li J. Study on 3-D full-scale inversion of complex gravity and magnetic model based on Marquardt algorithm (in Chinese). Chengdu: Chengdu University of Technology, 2005 . |
[21] | 姚长利, 郝天珧, 管志宁. 重磁反演约束条件及三维物性反演技术策略. 物探与化探 , 2002, 26(4): 253–257. Yao C L, Hao T Y, Guan Z N. Restrictions in gravity and magnetic inversions and technical strategy of 3D properties inversion. Geophysical and Geochemical Exploration (in Chinese) (in Chinese) , 2002, 26(4): 253-257. |
[22] | 吴小平, 徐果明. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报 , 2000, 43(3): 420–427. Wu X P, Xu G M. Study on 3-D resistively inversion using conjugate gradient method. Chinese J. Geophys. (in Chinese) , 2000, 43(3): 420-427. |
[23] | Tartaras E. LC1DINV: laterally-constrained 1-D inversion for processing airborne EM data. British Geological Survey, IR/05/118, 2005: 19. |