三维立体成矿预测中区域大比例尺建模技术研究与应用

2022-08-04宋相龙楚文楷

地质学刊 2022年2期

李楠宋相龙楚文楷

(1. 中国地质科学院矿产资源研究所，北京 100037； 2. 中国地质大学(北京)研究生院，北京 100083； 3. 中国地质科学院，北京 100037)

0 引言

20 世纪80 年代以来，寻找隐伏矿和深部矿已成为许多国家的主要找矿方向，地质找矿难度日渐增大，找矿效果日趋降低(陈建平等，2002，2014a)。随着矿产资源勘查深度的逐渐加深，传统的基于二维平面开展的矿产资源预测评价已不适应要求，如何更有效地找矿成为焦点问题，为此各国一直致力于研究新的成矿预测方法。随着计算机图形处理、三维可视化及三维插值技术的迅速发展，利用计算机三维地质建模进行三维成矿定量预测的研究逐渐展开，成为矿产资源预测的一大亮点(陈建平等，2014b；刘静等，2017) 。

三维地质建模是对地学多元数据集进行综合解译与可视化表达的过程(肖克炎等，2000)。仅凭有限的钻孔信息得到的地下三维地质情况往往只是“一孔之见”。根据对研究区地质情况的了解，形成地质概念模型，然后在不同剖面上进行地质体的圈定，不同的专家可能有不同的地质体建模方法。

三维地质建模是覆盖区与深部找矿预测的前提和基础。通过广泛收集、整理、转换研究区数据，最终形成各类地学信息在三维空间上的分布模型(肖克炎等，2012)。笔者研究了基于区域大比例尺且缺少钻孔信息情况下的一种隐式建模流程，有大量的地表实测数据但缺乏足够的深部数据(钻孔数据)支持。基于地表数据，结合地球物理数据，通过三维反演、地质统计学、几何算法与三维可视化等技术方法，直观地刻画出地下一定深度范围内的构造、岩体、地层等的空间展布以及成因和演化关系，为下一步找矿方向提供指导。

1 区域大比例尺隐式建模技术

1.1 三维地质建模的尺度

区域三维地质建模是一个多元地学信息综合解释的过程。一个能够真实反映深部地质情况的三维地质模型应由遵循地质规律的地质体及其空间关系(几何拓扑)构成(Xiao et al., 2015)，在此方面已经开展了大量研究工作。根据研究区范围(比例尺与深度)、原始建模数据以及使用的建模方法等的不同，将三维地质建模分为矿区、区域和岩石圈3个不同尺度类型(表1)。从建模数据分析，矿区尺度三维地质建模拥有大量地表、钻孔等数据，区域大比例尺拥有大量地表实测数据但缺乏足够的深部数据，岩石圈建模完全依赖于深部地球物理探测技术手段获取的数据，因此针对上述3个尺度的建模方法存在明显差异。从对地质模型的质量要求看，由于受到地球深部探测技术手段、经济效益和对于深部地质演化过程理解的限制，在符合现有地质理论的前提下，随着建模区域范围的扩大和深部已知数据的减少，对模型精度的要求一般呈降低的趋势。根据地质找矿理论，不同尺度所面临的研究任务也不尽相同。

表1 三维地质建模研究的3个尺度

国际上在上述3个尺度均开展了大量的研究工作。影响最大的是澳大利亚联合科学与工业研究组织(CSIRO)提出的“玻璃地球”计划，目的是利用三维可视化和地质建模等技术，使从地表到地下1 km的大陆表面“像玻璃一样透明”，同时通过三维可视化和广泛的网络服务接口传播信息填制澳大利亚信息丰富的四维“地图”，该项工作的主要目标之一是服务于澳大利亚下一代巨型矿床的发现(陈应军等，2014)。欧美各发达国家率先响应并制定了实施规划，开展了一系列研究与实践：美国计划在20世纪90年代绘制的三维地质图的基础上开展全境范围的5 km以浅的地质填图工作，主要技术手段为高分辨率三维地震成像等；英国计划构建全境多比例尺(1∶100万—1∶1万)的三维地层模型；俄罗斯、法国、荷兰、德国等也有类似的三维地质建模工作计划。

1.2 三维地质建模的算法

为实现上述三维地质建模工作，国外研究人员提出了大量的计算机算法，用于构建三维地质模型。根据输入数据与建模过程的不同，算法主要分为显式建模与隐式建模2类。

1.2.1 显式建模算法指在已知点集的基础上, 通过数学函数或人机交互的方式重构地质体，即数据驱动型算法(李楠等，2015)，其中已知点集为地质体曲面(重构曲面)上的全部控制点。显式建模的算法主要包括3大类。

(1) 应用空间Delaunay三角剖分算法重构地质对象。Cheng等(2000)首先提出了一种基于离散点集凸壳的方法, 其基本思想是基于点集构建四面体凸壳，称作α-shape，并应用Delaunay剖分在凸壳的约束下重构点集所在的曲面。之后又提出了一些改进算法(Bernardini et al.,1999; Gerritsen et al.,2002)。

(2) 样条函数拟合法。Hoppe等(1992)和Grimm等(1995)提出基于B样条的曲面重构算法，Park等(1995)提出基于Bézier曲面的近似曲面重构算法。该方法的最大缺点是无法强制要求结果曲面通过已知地质控制点，难以构建不连续或者发生合理突变等的复杂地质曲面, 例如侵入岩体、断层网络或复杂的接触界面等地质情况。

(3) 结合人机交互方式构建更加精确的地质曲面。最具有代表性的是轮廓线拼接法和离散光滑插值法(Discrete Smooth Interpolation，DSI)：轮廓线拼接法(Tipper,1976)是最为成熟也是应用最为广泛的三维曲面重构算法之一，应用实例如Surpac、Micromine等矿山勘探软件，核心思想是将相邻地质剖面上相同属性的点对通过人机交互的方式连接起来，最终构建地质体表面模型；Mallet(1992,1997)提出了基于人机交互、离散点空间Delaunay剖分与摩擦函数(Roughness)的曲面重构算法(MLPSR)。通过摩擦函数，DSI算法使拟合曲面S′通过合理的数学模型无限接近于原始曲面S，解决了结果曲面S′生成时遇到的控制点集(不仅仅是近似)、曲面光滑和不连续曲面的重构等复杂问题，该算法在区域地质建模领域得到了广泛的应用并推出了商业化软件GOCAD。

1.2.2 隐式建模算法与显式建模算法相比，隐式建模算法的最大不同是使用了空间插值和等值面提取算法，使该类算法仅需知道部分控制点便可推断出其他的控制点直至整个地质曲面，即数据驱动+知识驱动算法。随着带符号的距离场、Marching Cubes(Lorensen et al.,1987)以及MT(Zhou et al.,1995)等算法的相继出现, 隐式建模算法越来越受到重视，该类方法研究的核心是遵循地质规律的带符号距离场的空间插值问题。Savchenko等(1995)首次提出隐式曲面重构的数学表达式，即存在体函数φ(x,y,z)，当φ(x,y,z)=0时表示体的表面。此类方法主要应用于医学与地学领域。

隐式建模广泛采用的插值算法主要有地质统计学与径向基函数(Radial Basis Function, RBF)法等。Lajaunie等(1997)提出应用多变量地质统计学算法解决地层曲面的模拟问题；在此基础上，Calcagno 等(2008)提出了基于地质知识的隐式地质建模方法模型并讨论了单一断裂等构造曲面的提取问题。同时，服务于隐式建模的软件GeoModeller应运而生。

2 区域大比例尺隐式建模技术应用

笔者所用的隐式建模方法是基于对偶协同克里金插值的位势场方法，该法利用地质界面点等式约束和方向数据的梯度方向，采取对偶协同克里金插值法得到标量场，其中的一系列等值面视为地质界面。为隐式建模技术的算法基础，可以运用的地质约束归结为地质界面点的等式约束、方向数据的梯度方向约束、地质界面点之间的不等式约束等几类。

根据隐式建模技术的特点，三维建模主要使用的数据有矿区地表地质信息、音频大地电磁测深(CSAMT)剖面解译数据、航磁解译结果及剖面数据、钻孔数据等。建模流程(图1)由矿区原始数据采集、数据处理及数据输入、创建地质界面、生成曲面模型和生成实体模型5步组成(吴晓贵等，2020)。

图1 GeoModeller软件三维建模流程图Fig. 1 GeoModeller 3D modeling flowchart

2.1 原始数据收集

研究区矿床为一中低温热液型金矿床，矿体主要产于向斜的南翼地层1(炭质板岩)中，少量赋存于地层2中(图2)，受两侧侵入岩体挤压，在向斜南翼东段形成弧形构造，地层1及部分地层2的炭质板岩提供地球化学障，形成容矿空间。需要构建的赋矿地层为地层1和地层2。

成矿作用方式为含金流体在Ⅲ号脆韧性剪切构造带(图2)的控矿构造中充填交代，关键控矿因素是构造变形，根据音频大地电磁测深(CSAMT)L1剖面、航磁解译T图及航磁解译剖面S1—S5(图3)，向斜南翼核部段可能存在第二个成矿富集段，应作为深部找矿的主要目标。在地层1和地层2的建模基础上，向斜核部段的Ⅲ号脆韧性剪切带与地层1和地层2叠加区域应作为主要建模范围。

图2 研究区地形地质图1-地层1；2-地层2；3-地层3；4-地层4；5-地层5；6-地层6；7-物探剖面；8-脆性断裂带；9-脆韧性断裂带Fig. 2 Topographic and geological map of the study area

2.2 地质模型的建立

所用软件为GeoModeller，支持广泛的数据格式。建模所用数据有矿区1∶10 000地形地质图、1∶25 000航磁解译图、1幅CSAMT剖面图、1套航磁原始数据、5幅航磁剖面图、1套钻孔测斜数据。

2.2.1 工程创立首先定义建模所需的工程范围，然后给出建模准确度及导入该地区地形数据。主要围绕矿区控矿构造和赋矿地层进行建模(图4)，面积为44.85 km2，深度为2.8 km。

图4 工程定义示意图Fig. 4 Engineering definition interface

2.2.2 地层构建在建模前需先构建地质模块，如地层、岩基、构造等，并依据从老到新的准则对整个矿区构建一个序列。建模的地层从老到新为地层6—地层1(图5)，各地层之间的位置关系按地层的属性onlap(沉积)设置。

图5 地层构建示意图Fig. 5 Formation construction interface

2.2.3 地表数据导入地表数据主要包括地表控制线和产状信息。根据矿区地形地质图添加/设置各类型地质界线(图6)。在添加/设置各类型地质界线时对每个地质体依次描绘，逐步完善地质体界线及其产状信息。在描绘过程中要穿插进行计算调试工作，查看地质界线的准确性，进行校准。在地表控制线和产状信息全部输入完成后，可先计算一个初步模型，由于目前控制数据过少，所以模型计算的结果往往偏差较大，需进一步添加控制数据。

图6 地表数据导入Fig. 6 Surface data import

图3 地球物理综合解译图(a)CSAMT-L1剖面；(b)航磁T图；(c)S1—S5航磁剖面Fig. 3 Comprehensive geophysical interpretation(a) CSAMT-L1 profile; (b) Aeromagnetic T-diagram; (c) S1-S5 aeromagnetic profiles

图7 剖面数据添加(a)剖面添加；(b)剖面1添加；(c)L1剖面添加Fig. 7 Profile data addition(a) Section addition; (b) Section 1 addition;(c) L1 section addition

2.2.4 剖面数据增加为提高模型准确度，还需增加剖面数据控制以降低模型的不确定性。首先将所有需要的剖面数据的控制点坐标投射到地表上(图7a)，并对这些剖面依次进行命名，使其得以在二维窗口中单独显示；然后结合模型形态特征，在每条剖面的二维界面中添加剖面控制线和产状信息等。以剖面1(图7b)为例，由于该矿区深部还存在隐伏花岗闪长岩体，在剖面二维视图中添加音频大地电磁测深(CSAMT)L1剖面(图7c)，以推断隐伏岩体的控制线和产状信息等，并根据CSAMT反演的赋矿地层及控矿构造情况辅助修改模型，提高模型准确度。

2.2.5 三维模型计算与显示将上述设置好的数据导入三维模型进行计算，生成地质单元的三维地质模型。按照矿区实际情况，还可对剖面上推断的数据进行调整，以逼近真实地质场景。同样，在此部分还可设置计算网格的大小，也能手动输入需要计算的区域(图8a)，对于同样的模型，计算精度与消耗的时间成反比。

最终建立的赋矿地层——地层1、地层2及隐伏花岗闪长岩基底的三维地质模型见图8b。在X、Y、Z方向上均可进行不同距离的切割，也可在任意方向进行旋转，可进行单个层位和多个层位的显示。分别沿X、Y、Z方向切割相应的距离，可见模型内部任意位置的组合形态。模型可以导出其他形式的网格文件，为模型后续改进提供可能。

2.2.6 脆韧性剪切带建模将图3c的S1—S5航磁剖面导入地表，并在二维显示窗口中结合模型形态特征，在每条剖面的二维界面中添加Ⅱ号脆性断裂带和Ⅲ号脆韧性剪切构造带的剖面控制线和产状信息等，将剖面信息放入三维模型中计算，建立Ⅱ号脆性断裂带和Ⅲ号脆韧性剪切构造带的三维模型结构。图9所示模型即为褶皱构造南北两翼地层赋矿地层(地层1、地层2)和控矿构造(Ⅲ号脆韧性剪切构造带)复合显示图。