赵 勇,许 国,卢 鹏,文 诗 宝,黄 梅 婷

(1.南宁市勘测设计院集团有限公司,广西 南宁 530022; 2.南宁市浅表地质大数据工程技术研究中心,广西 南宁 530022)

0 引 言

三维地质模型不仅可以直观展现地质情况,还能辅助技术决策,在实际工作中发挥重要的作用[1]。创建模型的方法可分为显式建模和隐式建模两类。受传统地矿勘察设计思路、方法及工作习惯的约束,目前三维地质建模研究与应用多集中于显式建模[2];然而,传统显式建模工作量大、过程繁琐、效率较低、模型更新重构困难并存在拓扑错误、表面粗糙、棱角尖利等问题[3-4]。隐式建模基于函数插值空间离散点生成三维地质模型,具有建模速度快、实时动态更新、人工干预少、适应数据能力强、结果精度高等特点[5],但隐式建模软件费用昂贵且闭源,限制了隐式建模的推广和研究。国内外学者提出了多种用于隐式建模的空间插值方法,如集成在GOCAD、Leapfrog Geo、GeoModeller等商业软件中的离散光滑插值(DSI)[6]、快速径向基插值(FastRBF)[7]、势场理论插值[8]以及处于研究阶段的HRBF插值[9]、GRBF插值[10]、杨赤中推估法[11]等,但提供这些插值方法的地质建模软件均未对公众免费开放。为减少对商业软件的依赖,部分基于开源软件Blender、FreeCAD的研究能实现矿山三维自动建模及数字化应用[12]、水闸工程三维参数化建模和有限元结构计算[13],但软件未融入地质概念,用于地质建模容易出现违背地质规律的问题。少数学者利用Python研发地质建模软件,但因数据处理能力、算法优化、建模方式等局限不能完成复杂地质模型的创建[14]。

GemPy作为一款免费开源的隐式三维地质建模软件,能创建褶皱、断层、不整合等复杂地质模型,其插值算法可与商业软件媲美,并引入机器学习库Theano使大规模矩阵运算得以优化,且考虑建模的不确定性,利用贝叶斯推理框架实现随机地质建模和贝叶斯反演[15]。国外学者利用GemPy强大的隐式建模能力辅助科学研究,如Fandel等[16]将GemPy与SKS(随机岩溶模拟)、SWMM(暴雨洪水管理模型)集成,用于探索岩溶系统中管道网络结构与水力参数对地下水流动特性的影响;Thomas等[17]考虑到隐式建模能快速重构模型的特点,用Petrel生成二维曲面并提取曲面点数据导入GemPy中创建不同分辨率的三维地质模型;Jessell等[18]将GemPy作为三维地质建模引擎,研发从数据提取到模型构建全过程自动化的三维地质建模方法。国内对GemPy的应用与研究相对较少,孙丙阳[19]曾用GemPy创建三维地质模型,研究了基于MOOSE的地质模型体素化处理方法以及岩体工程开挖对连续-非连续变形的影响规律,但其采用的建模数据极少,且生成的模型较为简单和理想化;王博等[20]仅将GemPy用于地层表面结构数据计算。

虽然GemPy能创建复杂的三维地质模型,且模型能满足可视化及后续研究需求,有较强的实用性和拓展性,但目前仍多用于地层层序简单的三维地质模型创建;建模数据提取方法与数据结构化是使用GemPy的关键,但前人研究论文中鲜有提及,且未能梳理出清晰的建模流程。为此,本文以广西岩滩水电站为例,提出适用于GemPy的建模数据获取方法和详细可行的建模流程,创建受断层和风化剥蚀等地质作用影响的复杂三维地质模型,并分析其建模功能、原理、优点及局限性,旨在为隐式三维地质建模方法应用与研究提供一些思考。

1 GemPy概述

GemPy是基于Python语言编写的一款开源隐式三维建模软件,它融合了地层层序、断层、褶皱、不整合接触等地质要素,基于势场理论,采用克里金插值方法,能构造复杂的三维地质模型。GemPy通过调用第三方库来实现三维地质模型的创建和展示。软件框架如图1所示,可概括为以下4个模块。

图1 GemPy软件框架

(1) 数据读取与导入模块。GemPy将原始数据存储和转化为pandas数据结构用于后续计算,将原始数据的创建、导出等操作封装在类(Class)DateManagement中来初始化数据,以减少计算量。

(2) 地层层序及构造关系定义模块。该模块包含Series、Fault、Surface。Series用于定义地层层序、Fault用于指定断层与地层的切割关系、Surface用于定义地层面。

(3) 插值运算模块。使用Theano优化并高效计算出克里金参数,使用scikit-image库提供的移动立方体算法进行等值面提取和网格模型构建。

(4) 成果展示和分析模块。调用PyVista实现三维模型展示、调用Matplotlib实现二维成果展示,提供模型拓扑关系分析及根据地质模型与指定密度进行重力正演,使用PyMC3进行模型贝叶斯反演。

2 建模理论

GemPy三维地质建模应用的基本原理是势场理论[21-23]。该理论已成熟应用于商业建模软件GeoModeller。运用此理论进行地质建模,有以下3种假设:① 一个地质界面划分两个地层;② 地层内的方位数据与地质界面有关;③ 地质界面是同方向无限组平行的曲面。

该理论的原理是在三维空间中创建任意点p=(x,y,z)的标量函数T(p),使得所有T(p)=tk(tk为势场的某个未知值)的点形成一个等势面,将这个等势面作为一个地质分界面,则标量场中的每一个等势面将代表一个沉积界面,两等势面间形成一个地质建造(地层),图2为势场原理示意。

图2中,红蓝点为地层分界面上的采样点,箭头为方位数据。采样点通过的等势线(红蓝线)为地层分界线,地层分界线之间的绿色等势线可视为同一地层在不同时期的沉积界面(可认为是地质体趋势或地层内部层理轨迹),由红蓝地层分界线分隔出3个地层单元。从图可观察到方位数据可以是采样点上的,也可以是其他位置的。通过插值采样点和方位数据获得地层分界面。

根据势场理论,一个势场被定义为一个独立的地层序列(Series),在同一势场内等势面之间相互平行,各独立单元中的地质模型只由原始数据及插值算法决定。然而现实地质情况非常复杂,因沉积、侵蚀、侵入、中断等地质作用,地层之间相互切割、相互重叠,需要将势场理论与地质规律相结合来构建复杂的三维地质模型,可以通过定义多个势场的方式解决。

基于协同克里金实现势场理论插值算法,其函数表达方式如下:

式中:μα和νβ是p和p0的函数,由协同克里金确定的权重值。pα和p′α是位于同一地质界面上的2个采样点,M是增量T(pα)-T(p′α)的总个数;∂T(pβ)/∂uβ是N个方位数据点pβ在任意方向的偏导,代表标量场的变化率。势场理论插值算法实现模型重构的过程可概述为将位于同一等值面上的采样点及方位数据代入上式计算出权重系数后,反算出空间中未知点p与固定点po的距离增量。增量T*(p)-T*(p0)>0表示p点位于等值面外侧,T*(p)-T*(p0)<0表示p点位于等值面内侧,T*(p)-T*(p0)=0表示位于等值面上。为了计算出增量为零的点,采用Marching Cube移动立方体算法计算出立方体8个顶点的数值从而提取等值面,最终实现模型重构。

3 数据提取方法及建模流程

3.1 数据提取方法研究

基于上述理论,GemPy建模需要有已知的地质界面采样点(surfaces_point)和方位(orientation)两类数据。界面点坐标确定地质界面的位置,方位数据用于约束地质体的形状和趋势。界面点坐标可以通过钻孔直接获取,而方位数据(包括倾向、倾角和极性)获取比较困难,是使用GemPy建模的关键。方位数据最好是在有代表性、连续露头的地质界面处测量的倾向和倾角数据,极性用±1表示,沉积关系的地层+1表示指向年代较晚的沉积方向,侵入体+1表示由外指向内。

在方位数据的获取方法方面,王博等[20]利用最小二乘法拟合出能反映全局产状的地层平面并任取该平面内三点坐标来计算产状;周文辉等[24]将产状数据的求取转换为求取剖面线所在近似平面和剖面线上任意处的切向量,然后借助过渡向量,将切向量转换为法向量。受已有地质资料所限,本文根据原始地质剖面图提取所选岩性段处的方位角和倾角,其中方位角为剖面线方向,倾角为岩性段处地层分界线切线方向与水平面的夹角。基于开源软件webplotDigitizer在二维剖面下提取离散点数据,在Python环境下运用SciPy库将提取的离散点数据进行三次样条曲线拟合,再求取其一阶导数作为切向量。计算切向量与Y轴的夹角后,将其转换为法向与Y轴的角度。图3为在Matplotlib中二维剖面切向量提取结果。

3.2 建模流程

由于GemPy缺少交互式操作界面,需要编写Python代码来创建三维地质模型。运行GemPy首先应配置软件所需的Python环境,然后借助PyCharm或Jupyter软件来实现。建模的具体步骤如下。

(1) 数据准备及初始化。GemPy按照预定的数据结构进行计算,因软件无数据库管理系统,需要在外部以.csv文件按表1及表2结构存储原始数据。设置模型范围及分辨率,由于软件的限制,每个轴的分辨率设置值最大不超过100,即网格精度为模型各边长度范围值除以100。

表1 点坐标数据表结构

表2 方位数据表结构

(2) 创建Grid网格。Grid网格为模型插值运算提供空间坐标数据,是进行空间分析、数值模拟等地质应用的基础,是属性模型的载体[25]。在初始化阶段选择网格类型,GemPy提供5种网格类型(regular、topography、custom、section、centered)。regular为默认网格,用于常规网格的建立;custom能设置任意形状的网格;section用于二维剖面;centered用于计算物理属性等。其中topography能利用基于GDAL(开源栅格空间数据转换库)的数据进行地形网格创建。

(3) 地层序列的创建。GemPy在创建地层序列时遵循地层层序关系和断层时序关系,使用函数gp.map_stack_to_surfaces( )定义地层的顺序,地层新老关系按函数定义的先后顺序确定。断层需与地层同时设置且必须设置在地层序列的前面。对于地质情况复杂的模型,可以定义多个序列(Series),由Series中不同岩层组合关系来处理剥蚀、侵入等地质现象。

(4) 断层面的创建。在上述地层层序中定义断层面后,使用函数model.set_is_fault( )将断层激活,激活后可设置断层的接触关系。GemPy对断层的处理实质是在克里金函数中添加一个偏移项,断层面创建所需要的数据与地层一样,以断层轨迹离散点作为点数据、方位数据作为梯度进行插值。模型生成后,受断层影响的地层面会贴合断层面并产生急剧的弯曲来保证地层面的连续性,所以并未实质性将地层切断。

(5) 插值运算。三维建模需要使用插值方法将离散的数据转换为空间连续的数据,基于克里金函数的插值是GemPy的核心,使用全局插值使同一势场内的两个面不产生交叉。

(6) 二维、三维成果的展示。插值运算完成后,调用Matplotlib查看二维剖面、调用PyVista查看三维模型。

(7) 拓扑分析及重力正演。GemPy内置模型拓扑关系分析功能,使用函数gp.topology_compute( )可进行拓扑分析。通过创建centered_grid,使用GravityPreprocessing( )函数可进行重力的正演计算。GemPy建模流程如图4所示。

图4 GemPy建模流程

4 示 例

4.1 研究区地形地貌及岩土层情况

广西岩滩水电站坝址工程所在地区为低山峡谷地形,两岸高程在500～750 m之间,底部河流受区域构造控制,河床狭窄,水流湍急。周围地层岩性可概化为石炭系灰岩、二叠系辉绿岩,其中二叠系辉绿岩按形成的年代由老到新分为5层(用1_1～1_5表示),场地地层被一条逆掩断层切断。

4.2 模型建立

(1) 数据准备。本文所提出的产状数据提取方法需要有已知的剖面图且产状数据获取较为繁琐,无法获得未被剖面穿过区域的数据,为验证GemPy在复杂地质条件下的建模能力,示例采用的点数据及方位数据是在其他软件平台创建三维模型后提取。模型创建选取深度较深并能较完整反映场地地层情况的钻孔,使用了点坐标数据235条,方位数据139条。

(2) 初始化数据。按上述数据结构存储建模数据,根据场区的范围设置x,y,z3个方向的起止值,并将分辨率设置为软件建议的最大值。采用默认的regular网格来创建地质模型,选择topography网格并通过gdal读取ASCII Grid(.asc)文件创建地形面,并预设2个剖面网格(section)。初始化数据后,可以通过PyVista在三维坐标下查看数据分布情况(见图5)。

图5 数据空间分布

(3) 创建地层层序。根据场地地层及构造情况将模型的地层分为3个序列,Series1为断层上盘(5层辉绿岩)、Series2为断层下盘(5层辉绿岩)、Basement为软件默认的基岩面(石炭系灰岩),并将断层上下盘的岩性更改为相同的颜色,如表3所列。

表3 地层层序配置

(4) 地层接触关系处理。由研究区地质资料可知,该地区地层被逆掩断层错断后,上盘地层又经受了风化剥蚀,所以底部岩层在地表出露。若采用GemPy提供的断层功能,无法重构岩层被剥蚀出露的状态。受限于软件的功能,从建模效果考虑并未采取让断层切割地层的方式,而是单独将地层分为上盘和下盘两个地层序列(Series)进行建模。

(5) 插值运算。上述步骤完成后,对数据进行克里金插值和模型计算,最终运行结果参数如表4所列。

表4 模型计算参数

(6) 模型展示。GemPy将PyVista内置在程序内部,通过函数gp.plot_3d( )可以查看创建的三维模型。PyVista提供简单的三维操作功能,用户可以进行模型的旋转、地层面隐藏与显示等。图6～8分别是通过GemPy生成的三维地形面模型、三维地层模型与三维地质模型。

图6 三维地形面模型

调用第三方库ipywidgets,按单元来展示地质模型在剖面上的地层分布,图9～10分别为沿Y轴与沿X轴方向单元格(cell)为40,50,60位置处的地层剖切。

图9 沿Y轴不同位置处的剖切面

图10 沿X轴不同位置处的剖切面

基于势场理论,GemPy利用点坐标与方位数据通过克里金函数自动插值生成地形面与地层面,进而形成三维地质模型。建模过程无需过多的人工交互,建模效率和自动化程度较高。从实际建模效果来看,创建的地层曲面光滑,可视化效果好;创建的地质模型与实际地质情况相似,并能反映出逆掩断层、地层剥蚀出露等复杂的地质现象;创建的剖切面可以观察不同方位岩层的空间分布情况。

5 结 语

通过对GemPy的功能、建模原理进行分析并结合实际工程进行三维地质建模,总结分析出GemPy的优缺点如下。

(1) 能实现复杂地质情况的三维地质模型及任意方向的二维剖切面创建与可视化。因软件采用机器学习及贝叶斯推理,使创建的模型可信度较高,能反映实际地质情况。

(2) 建模数据获取困难,但软件并没有提供解决方案。通过在外部提取已有地质资料的方法,效率很低,精度也不够,如何获取建模数据成为使用GemPy建模的关键。模型分辨率低,网格分辨率受限,软件不能设置精细网格。

(3) 无交互式界面,需要用Python代码实现交互功能,导致建模操作不够便捷。没有独立的数据库管理系统,数据的存储与调用不方便。

(4) 软件算法需进一步优化,由于三维地质建模处理的数据量非常大,GemPy在进行克里金插值运算时需占用较大内存,无法制作分辨率较精细的模型。同时,利用Python语言编写的软件进行插值运算非常耗时,运算效率较低。