基于开源GemPy的城市地下空间三维隐式势场建模方法研究

Research on the 3D implicit potential field modeling method for urban underground space based on the open-source GemPy

摘要

三维地质建模是城市地下空间规划和管理的重要基础工作。在地质科学大数据时代,当前三维建模逐渐转向三维隐式建模。常见的三维地质建模软件如GeoModeller、Surpac等都提供了隐式建模技术,但是上述商业软件均封装完好、并不开源。本文基于开源GemPy建模平台,重点实现了三维隐式势场建模技术与针对地层尖灭的建模方案,将算法方案加入GemPy底层建模架构,并以某城市地下空间为例,对三维地质精细化建模工程流程进行了验证。论文构建了基于研究区钻孔数据的建模数据集,采用势场法和泛克立格插值法建立城市三维地质模型,并针对城市地质中常见的地层尖灭进行了单独的布尔运算处理。为细致直观展示研究区三维地质模型地层界面的拟合情况,对共计18个工程地质层分别进行可视化。最后,对模型的精度采用分层K折交叉验证方法进行验证,并用皮尔逊相关系数(CC)等4种指标对建模效果进行衡量。对研究区的实验结果表明,基于开源GemPy的三维隐式势场地质建模方法针对城市地质构造具有一定通用性,并在精度上具有较好效果,可以为城市地下空间开发利用提供决策依据。

关键词

三维地质建模 / 隐式势场法 / GemPy / 地层尖灭 / 城市地下空间

Key words

three-dimensional geological modeling / implicit potential field method / GemPy / stratigraphic pinch-outs / urban subterranean space

0 引言

城市化进程的不断加快^[1]导致城市规划用地日益紧俏,启用、开发城市地下空间成为当前城市规划的最优解之一^[2]。城市地下空间规划利用和实际开发极具复杂因素,在进行实际开发前需要做好详实的规划^[3]。三维地质模型增强了城市地下空间的真实感、可触感^[4]和立体感。对模型的多重分析、仿真模拟、综合设计和决策基础^[5],能够支撑城市地下空间的资源评估和空间评价、规划、开发和利用,为智慧城市建设服务^[6]。

地质建模方法可分为表面重建和体重建^[7],以表面重建为主。常见的三维地质建模方法有两种:显式建模和隐式建模^[8⇓⇓-11]。较之显式建模会受到数据密度、地质特征复杂性、模型分辨率和建模剖面方向的强烈影响^[12],三维隐式地质建模具备以下优势^[13]:可根据应用需求生成不同尺度的高质量三维模型,通过调整约束规则获得建模结果,对多种不同类型的地质数据同时进行建模,融合多域模型且保证模型之间较好的接触关系,融合多种不同方法获得建模的结果,如将显式建模的结果融合到隐式建模模型中^[14]。

三维隐式建模的核心技术是隐函数插值方法^[15]。离散平滑插值方法(DSI)^[16]适用于从地震数据或激光扫描仪获得的分散的、无组织的点数据,隐式重建具有复杂几何形状的断层表面。该方法主要应用于油气藏领域。相较于城市地质建模应用,基于径向基函数(RBF)的方法^[17-18]更适合被应用于点云重建隐式曲线和曲面。基于Hermite径向基函数(HRBF)的方法^[19]允许通过界面从建模区域自动提取三维地质模型统一基于RBF的曲面重构方法。然而,非监督地下建模的合理性和准确性将极大地受限于地质专家的解译和修改。而带部分约束的HRBF方法^[20]正是为地质专家绘制剖面约束线提供灵活有力的手动交互式工具,可将剖面线的几何特征抽象为坐标和法向量,用来共同计算HRBF方程的隐含地质表面函数参数。但是由于此处的法向量是使用该点及其附近点的近似平面计算的,对城市地质中具有较大变化的剖面线的情况难以处理。移动最小二乘法(MLS)^[21]在曲面拟合阶段虽可自动忽略异常值,但二阶多项式缺乏先验知识将大幅降低重建的可信度。快速多极法(FMM)^[22]在拟合过程中尽管带来显著的压缩和进一步的计算优势,但并未解决RBF算法本身的固有缺陷。快速径向基函数法(FastRBF)^[23]从一系列轮廓切片构造曲面,但对模糊条件的统一化处理在针对断层等不整合构造时,效果并不理想。提出于20世纪90年代的等值面提取法^[24]认为,可通过定义地质特征之间的接触关系引入矿化趋势和构造趋势。

当前大量的商业化软件,如GeoModeller、Surpac、DataMine集成了上述隐式插值方法,具备了隐式建模能力,但底层代码并未公开。部分基于开源软件FreeCAD、Blender的工作实现了水闸工程三维参数化自动建模^[25]、智慧矿山三维自动建模应用场景^[26],但受限于并未融入地质学背景,容易出现违背地质规律的问题^[27]。GemPy是一套开源的三维地质建模工具包,与技术成熟、功能强大的商业建模软件GeoModeller相仿^[28],融入了构造地质学概念,同时代码和底层结构完全开源,有较强的实用性与拓展性,相较于封装完整的商业建模软件,建模核心算法有较大改善和补充空间,可较好地灵活应对城市地下空间的复杂地质情形。GemPy的三维隐式势场建模方法可有效构建简单层状地质模型,自动化处理断层等构造,但受限于隐式势场法的全局插值属性,所得到的连续趋势面会忽略地层尖灭这类普遍存在于城市地质环境中的特殊构造。

为了进一步研究隐式建模技术与开源GemPy在城市地下空间中的应用,本文介绍了基于GemPy的三维隐式势场地质建模方法。此项工作从研究区的城市地下空间钻孔勘探原始数据的输入开始,对钻孔勘测得到的原始数据首先进行初始处理,针对常见的浅孔深部缺乏控制等情形进行数据预处理,得到建模数据集;采用势场法和泛克立格插值法构建标量场模型,并使用移动立方体算法构建几何块体模型;针对城市地质第四系地层中广泛存在的地层尖灭问题提出布尔切割运算的建模方案;对建模数据集进行二维可视化、三维地质建模,最后采用分层K折交叉验证方法对模型精度进行验证。

1 基于GemPy的三维建模原理

1.1 建模流程

GemPy具备创建复杂地质模型的能力,可以实现层序特征和结构特征的地质模型构建任务,例如,背斜、向斜类的褶皱构造,考虑相互接触关系的复杂断层和各种不整合的复杂模型,实现多个地质体之间的拓扑关系建模。

整体工作流程如图1,获取研究区的工程钻孔勘探原始数据,对源数据进行建模原始处理,并考虑由于钻探实际工作等多方面现实因素导致的、数据集中可能出现的浅孔深部缺少控制问题进行预处理,构建建模数据集。以布尔切割运算为建模方案处理地层尖灭,以隐式势场法、移动立方体法和泛克立格插值法作为建模核心方法进行三维地质建模。在模型分析过程中,对地质模型进行可视化;对不同体元大小的格网数条件下模型的平滑程度进行实验,并对三维地质模型各地层及部分地层的深度图进行逐层展示。最后,对模型精度进行量化衡量。

1.2 构建建模数据集

1.2.1 钻孔源数据的预处理

此次工作以中国某城市作为研究区,收集的钻孔数据共240个,以较为分散均匀的空间分布陈列在研究区范围内。钻孔位置分布如图2所示。表1列出了此次工作中对研究区地质地层的划分,并描述了所处的地层单位、岩性等。共划分为18个工程地质层,为便于后续建模输入地层层序关系,分别命名为SA1_0至SA10(SA:Study Area)。

数据类型划分上,研究区钻孔点数据包括3个工作表:(1)钻孔地理位置工作表。主要表征和记录的是研究区钻孔所在的地理位置,包括Hole ID(钻孔编号)、East(东坐标)、North(北坐标)、Elevation(高程)四项数据,东坐标、北坐标为高斯克吕格坐标。(2)钻孔详细信息工作表。主要表征和记录的是研究区钻孔的深度及相对应的地层序列号,包括Hole ID(钻孔编号)、From(钻孔起始深度)、To(钻孔抵达深度)、Geology(所属地层序号)四项数据,第二、三项数据差值即为钻孔的深度。(3)钻孔附加信息工作表。主要表征和记录的是研究区钻孔的倾角、方位角、长度信息。包括Hole ID(钻孔编号)、Dip(倾角)、Azimuth(方向角)、Depth(深度)四项数据。

为尽可能精确描述地层形态,选用钻孔底界深度以及地层产状两类数据决定地层的走向姿态。建模数据集的几何数据类型包括两类:(1)地层交接点(surface points);(2)地层产状测量值(orientation measurement)。钻孔与每一层地层底面的交接点作为地层交接点,钻孔与上、下两地层底面的交接点深度差值即下地层的厚度,具体对应规则如下:(1)交接点的xy坐标=所在钻孔口的xy坐标;(2)交接点的高程坐标值=钻孔口的绝对高程-钻孔向下的深度;(3)交接点所属地层=钻深抵达的地层。

1.2.2 建模数据集的构建

由于地形复杂程度、钻探的施工目的、矿区实际情况等等因素的影响,钻孔数据深度不一,部分钻孔缺乏在深部地层的数据控制,如图3所示,由于钻孔Zk3未能深入至地层D,无地层控制点,因此建模软件直接输出黑色实线的建模效果:地层C在右侧下倾,与原本地层D所在位置重合。对照实际地质解释情况,参考手工绘制的钻孔剖面图可知,红色虚线所示才是推断的地层走向形态。

结合研究区数据,如下图4(a)选取了模型中较为典型的一处举例,由于钻孔N92_ZK18未能深入至地层SA8_1(黄色),因此建模结果显示,地层SA8_1在钻孔N101_ZK18控制点处呈现突然下陷的走向,连锁反应导致地层SA10(绿色)出现异常凸起,如图4(a)中的红圈所示。结合研究区实际地质背景,地层SA8_1的地层单位属于第四系中更新统下组,成因是冲湖积,岩性杂色黏性土层,可硬塑,局部含砾;地层SA10的地层单位属于前第四系,岩性是基岩,岩性由不同时代不同类型的岩石决定,根据表1,地层SA8_1、SA10应呈现严谨的上下沉积层序关系。因此,需要对数据做浅孔深部控制的预处理,增加虚拟控制点。

DeepFit算法^[29]结合了神经网络来学习加权最小二乘多项式曲面拟合得到的逐点权重,学习的权重充当表面点的邻域。该算法基准法线和曲率估计数据集上取得了较好的拟合结果,展示了其在异常值消除、三维模型曲面拟合等方面的应用。采取DeepFit算法预处理原始数据集中存在的浅孔深部控制问题后,得到的效果如图4(b)所示。

为了清晰地看到模型在深部地层处的处理效果,对模型进行了z轴方向上的翻转,红色点表示追加的虚拟控制点。预处理后的地层起伏形态较为平滑,与图3剖面图中所示的红色曲线展示的效果相似。

1.3 三维地质自动化建模原理

建模工作原理如图5所示。建模核心原理为:输入地层分界面控制点数据、地层产状数据,以及克立格参数,整合成系列地层数据子集后,建模平台的概率建模功能^[30]将数据集概率分布采样得到的扰动输入数据参与建模。用泛克立格法进行插值得到标量场模型,通过设置格网参数分割标量场模型得到离散块体模型,并加入地层系列参数构建最终的标量场模型。定义地层尖灭点并使用布尔切割运算处理地层尖灭,最后利用移动立方体算法得到三维地质构造模型。

1.3.1 构建标量场模型

1.3.1.1 势场法^[31]

理想状态下,城市地质沉积环境是由一系列水平地层序列组成的,地层序列按照一般规律从老到新依次向上累叠^[32]。然而,真实的地层常常出现断层、侵入、褶皱等特殊地质构造。因此,生成三维地质模型也是为了表征具有显著物理性质变化的相邻岩性界面。常见的三维建模是针对每一个地层单独进行插值,往往出现地层穿层等情况,而势场法利用全局插值器优势在于:(1)在同一沉积环境中,地层之间的位置会相互影响,这使得处于同一势场的两层不可能出现地层穿插;(2)全局插值器能够利用目标地层之间的数据,拓展可能的测量范围,便于在后续插值中使用。全局插值器的使用前提是在同一势场下,两个相邻地层面不会交叉沉积,这牵涉到势场法的概念。

设有构造插值函数G(x),其中x表示连续三维空间(x,y,z)∈R³中的任何点,将域D描述为标量场。标量场中的每个等值面代表该地层对应的其他同步沉积地层面。梯度将垂直于等值面,并且可以位于空间中的任何位置。标量场的梯度为∂G/∂u,其中t是任意单位向量。为表征插值中的标量场,用到两种类型的参数:地层界面控制点x_α描述了对应的相邻两个地层之间的分界面;标量场的梯度x_β垂直于地层面,并可以位于三维空间(x,y,z)∈R³中的任意位置。标量场Z及其梯度之间的关系如下:

(1) \frac{\partial G}{\partial t}(x)= \underset{\rho \to 0}{\mathrm{l}\mathrm{i}\mathrm{m}} \frac{G(x+t)-G\left(x\right)}{\rho }

其中,标量场的属性是无量纲的,唯一的数学约束是该值必须在梯度方向上增加,即描述了地层沉积过程。因此,存在两个需要保留的约束项:(1)属于确定控制地层分界面 x_{{\alpha }_i}^q的所有控制点都必须具有相同的标量场值,即存在一个标量场等值面连接了同属地层分界面的所有控制点,如公式(2),其中是地层分界面的参考点;(2)标量场将垂直于三维空间中任何位置的梯度x_β。

(2)G( x_{{\alpha }_i}^q)-G( x_{{\alpha }_0}^q)=0

图6展示了标量场及其梯度在地层面中的关系。输入数据由分布在3个地层的8个控制点以及3个标量场梯度构成。一个等值面连接几个控制点,标量场垂直于梯度。

1.3.1.2 泛克立格插值法

克里格法是地统计学中一种经典的插值方法,适用条件是区域化变量存在空间相关性,而根据空间场是否存在偏移(drift)可将克里格法分为普通克里格法和泛克立格法。泛克立格法主要针对非平稳插值问题^[33],其主要优点是能够利用地质空间中不同地层位置采样的方向进行估计。在大多数地质环境中数据呈线性趋势,即一层的平均厚度在整个地质区呈线性变化。这种趋势可以用多项式漂移函数在普通克立格插值的方程中描述,即泛克立格插值方法。上文提到了构造插值函数G(x),若该函数变量在研究区内的数学期望E[G(x₀)]是变化的,则G(x)是非平稳的:

(3)

E[G(x)]=w(x)

非平稳随机变量包括偏移和波动。偏移是指随机变量G(x)在点x处的数学期望,即公式(4)中的w(x),波动R(x)是指点x处随机变量G(x)与偏移w(x)的差。两者的关系为公式(4)^[34]:

(4)

G(x)=w(x)+R(x)

偏移w(x)可以理解为随机变量G(x)在研究区某种明确的变化趋势,而波动R(x)则是随机变量G(x)在w(x)附近的随机摆动,波动R(x)的数学期望值为0。

(5)

E[R(x)]=E[G(x)]-E[w(x)]=0

偏移w(x)一般使用一次或二次多项式表示:

(6) \mathrm{w} \mathrm{(x)} = \stackrel{k}{\sum _{l=0}} {\mu }_l f_l \mathrm{(x)}

式中f_l(x)为已知函数,μ_l为未知系数。

考虑到地质建模的本质是对地层分界面进行曲面插值,漂移通常采用公式(5)表示:

(7)

w(x,y)=μ₀+u₁x+u₂y+u₃xy+u₄x²+u₅y²

随机变量G(x)在预测点x₀处的真实值G(x₀)未知,因此需要用x₀领域内的已知测量值的加权平均值来估计,λ_i为点x_i处的权重。

(8) G^{\mathrm{*}} \left( x_0 \right) = \stackrel{n}{\sum _{i=1}} {\lambda }_i G \left( x_i \right)

而为追求插值的精确性,预测点x₀出的真实值G(x₀)与估计值G^*(x₀)之差应满足无偏条件:

(9)

E[G^*(x₀)-G(x₀)]≡0

结合公式(8)、公式(9)可得

(10) \begin{array}{l}\mathrm{E}\left[G^{\mathrm{*}}\right(x_0\left)\right]=\mathrm{E}\left[\stackrel{n}{\sum _{i=1}}\right]{\lambda }_{i}G\left(x_i\right)=\stackrel{n}{\sum _{i=1}}{\lambda }_{i}w\left(x_i\right)\\ =\stackrel{k}{\sum _{l=0}}{\alpha }_l\left\{\stackrel{n}{\sum _{i=1}}{\lambda }_i{f}_l\left(x_i\right)\right\}=\stackrel{k}{\sum _{l=0}}{\alpha }_l{f}_l\left(x_0\right)=w\left(x_0\right)\end{array}

若要使得公式(10)对任何一组系数均成立,则必须满足

(11) \stackrel{n}{\sum _{i=1}} {\lambda }_i f_l \left( x_i \right) = f_l \left( x_0 \right) , l = \mathrm{0,1} , \dots , k

公式(11)即为用k+1个方程表示的无偏条件。

用变异函数表示估计方差表达式如下, {\sigma }_E^2为估计方差。

(12) {\sigma }_E^2 = - \stackrel{n}{\sum _{i=1}} \stackrel{n}{\sum _{j=-1}} {\lambda }_i {\lambda }_j \gamma ( x_i , x_j ) + 2 \stackrel{n}{\sum _{i=1}} {\lambda }_i \gamma ( x_i , x_j )

采用拉格朗日乘数法,将无偏条件公式(12)代入公式(13)可得目标函数式:

(13) \mathrm{F} = {\sigma }_E^2 - 2 \stackrel{k}{\sum _{l=0}} {\mu }_l \left[ \stackrel{n}{\sum _{i=1}} {\lambda }_i f_l \right( x_i ) - f_l ( x_0 \left) \right]

将公式(13)对λ_i及μ_l求一阶偏导,令其为0,便可得到泛克立格方程组。

(14) \left\{\begin{array}{l}\stackrel{n}{\sum _{j=1}}{\lambda }_j\gamma (x_i,x_j)+\stackrel{k}{\sum _{l=1}}{\mu }_l{f}_l\left(x_i\right)+{\mu }_0=\gamma (x_i,x_0)\\ \stackrel{n}{\sum _{i=1}}{\lambda }_i{f}_l\left(x_i\right)=f_l\left(x_0\right),l=\mathrm{1,2},\dots ,n\\ \stackrel{n}{\sum _{i=1}}{\lambda }_i=1\end{array}\right.

1.3.2 移动立方体算法构建地质模型

三维地质建模的目标是运用泛克立格插值法得到随机函数与其梯度,通过插值得到空间中任一点x的标量场的值,构建标量场模型,分割得到地质构造的表面模型,用以定义地质结构的空间分布。为了分割标量场模型,采用Lorensen等^[35]在Scikit-Image图像处理库中提供的移动立方体算法(Marching cube)等值面追踪算法,获得地质构造的几何表面模型。为在三维离散数据场中通过线性插值逼近目标等值面应做如下工作:(1)首先,体素化标量场模型,使其分散在三维体素小立方体中,逐个计算比对要提取的地质等值面的值是否落在了单个体素立方体的表面、体内或体外;(2)对落在体内的体素立方体的顶点进行插值,计算体素立方体顶点与目标等值面之间交点;(3)分析这些交点,并比较8个顶点值与目标等值面值的关系,根据移动立方体算法的索引表进行加权体素计算;(4)生成立方体的三角面和法线,合并上述三角面,输出地质界面、断层等构造面,最后输出几何块段模型。

1.3.3 地层尖灭复杂构造的建模方案

GemPy建模平台可以针对断层等复杂地质情况给出相应的建模方案,实现层序特征和结构特征的地质模型构建任务^[31]。上述隐式势场建模方法采取全局插值器,优势是自动插值形成连续趋势场,避免所建模型出现穿层交叉。但为得到连续的势场面,除非单独处理,否则隐式势场法将自动补齐地层中由于尖灭导致的缺失部分,无法体现沉积地层的真实地况。与此同时,包括此研究区在内的城市地质环境大多处于第四纪地质年代,地质条件复杂,地层尖灭现象较为普遍。要使基于GemPy的三维隐式地质建模方法应用于城市地下空间,必须单独设计算法解决对地层尖灭的建模问题。

此次工作利用GemPy代码开源、算法易补充完善的优势,在核心建模架构中提出针对地层尖灭特殊构造的建模方案。地层尖灭定义了一个地层块体向一侧结束延伸、在地质学上闭合并消失在另一个地层表面上的几何形状^[36]。地层的不连续性会导致上、下边界在地质剖面图上突然消失并生成尖灭,因此,可以将对尖灭的研究转化为两对相邻钻孔之间地层边界的研究。布尔切割运算在处理地层边界问题时十分有效。该方法的核心思想是将出现地层尖灭的区域按照地层层序分割成上下两个地层,并在地质模型中分别建立这两个地层。

具体算法思路如下:(1)尖灭定义。在获取建模数据集后,借鉴GemPy的三维隐式势场建模方法处理断层问题的思路,首先需要明确定义建模数据集中出现地层尖灭的区域。这包括定义地层尖灭处的几何形状、地层尖灭点的属性和参数。通过给相应的地层控制点赋予“pinch-out”属性,并以地层尖灭点的深度值作为下一步中切割地层的阈值。(2)逻辑判断与属性分配。对于位于地层尖灭区域的模型控制点,通过比较其深度属性,确定其位于地层尖灭点上方还是下方,并相应地分配地层属性“formation”。在尖灭点上方,该点的地层属性取自上层地层模型;在尖灭点下方,该点的地层属性则取自下层地层模型。(3)布尔切割运算。在建模完成后,通过布尔取交集的运算(AND运算用于确定每个数据点是否同时满足尖灭点深度、尖灭区域xy阈值等多个条件;WHERE运算将尖灭点上下两个地层进行组合),实现在尖灭点处对上下两个地层模型进行精确地切割。上层模型中大于阈值的部分被移除(remove),将移除部分对应的属性值设为特定值NaN(Not a Number),而下层模型则覆盖了被移除部分。在尖灭点上方,保留了上层地层模型;在下方,保留了布尔切割后的下层地层。

布尔切割运算实际上是一种对地层属性进行动态调整的方法,通过逻辑关系在尖灭点上下两侧分别保留或移除地层属性值。这一步骤模拟了地层尖灭点的位置,使得在地质模型中能够准确地捕捉到地层的变化和尖灭现象,通过逻辑运算来决定最终模型中该点的属性。整个模型在地层尖灭点上、下两侧形成了清晰的分界线,准确地反映了地层尖灭导致部分地层缺失的现象。

此外,基于GemPy三维隐式势场建模方法引入了随机场模型,允许在地层模型中引入一定程度的不确定性和变异性,以更好地模拟断层的复杂性。在处理地层尖灭的算法中,同样引入随机场模型,在使用布尔切割运算得到地层尖灭点的形态后,其领域内相对连续的地层场模型由随机场模型进行控制,以保证模型的整体性。

图7展示了对研究区进行三维地质建模时的局部模型剖面,地层SA6_3_4与地层SA7_1的分界面处存在地层尖灭点,使用上述算法思路进行处理,得到地层尖灭点与其分割的上、下地层。同时,对地层尖灭的处理结果在图11每一层地层的建模结果中也得以呈现,地层尖灭导致地层面的部分缺失。

2 实验结果

2.1 二三维建模可视化结果

2.1.1 建模数据集二维剖面可视化

根据研究区建模数据集进行二维剖面可视化,在x、y、z方向上根据自定义切片轨迹或单元格数绘制模型截面,输入的控制点数据和控制地层产状的数据可以包括在内,控制点的箭头即代表该点处的地层产状,如图8。图8(a)、(b)分别展示了从x,y轴方向进行纵切得到的研究区建模数据集的二维剖面,图中不同颜色的线表示了各地层的勘探线,部分锯齿状源于模型网格精度限制。图9展示了研究区建模数据集在三维空间的分布情况。

2.1.2 模型三维可视化

由于钻孔深度取值区间为[0.3,142.2],相较于x,y轴数值而言分辨率较小,为了使模型在z轴方向的可视化效果更清晰、更细致地展示模型内部结构,将三维地质模型在z轴方向上的垂向深度值等比例夸张显示(等比例放大系数A=100,夸张显示系数e=2),即该模型显示的垂向深度值为实际地层厚度的200倍。如图10,分别展示了模型无透视三维视角下,x、y、z轴顶部、底部的视角下的可视化结果,建立的模型较好地反映了研究区的地质构造特征。

同时,图11(a)提供了研究区三维地质模型的相应地层面。为了便于理解模型的构建原理和地层的空间分布,本文对实体模型提供了多种显示模式,如线框模式(图11(b))。为清晰显示每一层地层面的拟合建模情况,将其分别单独可视化,如图12,其中,地层SA2、SA3、SA4_2、SA5_3和SA6_3_4的部分缺失是由于地层尖灭导致的,与GeoModeller所建模型进行对比效果基本一致。

此外,模型的平滑度取决于基本建模元素的精度(即体元大小)和模型的三角形面片数量。在本研究中通过自定义模型体元大小(即设置模型格网的分辨率)来控制网格模型的精度。图13展示了不同的平面精度值(取决于DEM网格点的个数)生成的网格模型。显然,(a)显示的模型平滑程度和精度均优于(b)。不同平面精度的结果表明,提高体元的平面精度,可以提高模型的可视化质量。然而,这种改进需要更高的内存消耗。因此,综合模型平滑度与内存消耗两个标准作为三维地质建模精度约束,一般在城市工程实际应用中控制分割体元立方体数的范围落在[62 000,100 000]区间较为适宜。

深度图是地质学家可视化地层空间分布的重要工具。为细致直观地展现研究区三维地质模型每一层地层底面的曲面拟合情况,图14展示了模型SA0至SA2共4个地层面的深度图,z轴代表的深度值以图例颜色表示。

2.1.3 断裂构造与地震效应相关分析

考虑到三维地质建模的目的是服务于城市地下空间的实际开发利用工程应用,而城市地质环境中的断裂构造与地震效应关联安全^[37],因此有必要针对本模型结果开展相关分析。

断裂构造是岩层受地应力作用发生破裂的特殊地质构造,对于城市地下空间的开发具有明显的不利影响。对于计划穿过断裂破碎带的暗挖的区间隧道或车站,将被列为设计和施工的重点区段,开挖和支护需要采取必要的措施防止塌方。断裂破碎带中,尤其是断层上盘常富含流通性好的地下水,若实际开发施工时仅以排水为主,极有可能导致周围地层泱水积攒而产生严重下沉或局部地陷。因此,需要依据断裂构造所处位置“因地施策”处理所带来的富集地下水:对于出现在深基坑侧下方的部分,需要防范承压水冲破基坑冒出的风险;对于出现在深基坑底部的部分,应规避采取抽排水措施,因为将极易引起基坑旁的地面不均匀下沉;对于出现在竖井中的部分,需要对竖井壁建设好防坍护壁并着重考虑地下水。其次,在断裂破碎带区域对基坑进行开发开挖也极易受到地质断裂带中沿岩石裂隙面滑动的滑动力的不利影响,从而带来风险。断裂活动可能对研究区内第四系覆盖区的全新统可液化砂层和可能发生震陷的淤泥层产生重要影响,甚至造成地基失效。此外,断裂构造发生地带多伴有隐伏溶沟、溶槽、漏斗等地质“空洞”,会改变地质应力分布状态,使土体经开挖后处于松散状态而发生坍塌。

在实际评估工作时,将距断裂和地裂缝带的距离作为重要指标,并结合三维距离场分析方法将断裂带作为中心区作三维缓冲区,对待评估区域进行开发适宜性强度的划分。活动断裂与地震活动紧密相关,一般采用地震液化(严重液化、中等液化、轻微液化、不液化)、活动断裂(强烈全新活动断裂、微弱全新活动断裂、中等全新活动断裂、非全新活动断裂、无活动断裂)、地震烈度(>9度,9度,7~8度,6度)等指标进行评估^[38]。

对研究区三维地质模型(图10、图11)的内部形态及地质模型每一层地层的具体形态(图12)进行观察可知,区内模型整体走向较连贯,无明显断裂带存在;结合研究区历史地质资料和基础地理数据综合分析,该区属冲海相沉积平原区,无深大断裂通过,也无重大和频繁的地质灾害发生,地下空间开发条件总体良好。

2.2 模型精度的交叉验证

为评估上述所建立的研究区三维地质模型的精度,本研究对比深度图中提供的三维模型中各地层控制点的深度值与原始建模数据集对应数值,选用分层K折交叉验证方法(Stratified K-Fold),将研究区建模数据集随机划分12个子集(即12折),每个子集包含20个钻孔数据,每折包含每个目标样本数据近似相同的百分比。在其中任意11个子集上迭代训练算法,进行建模,同时使用剩余1个子集作为测试集。最后,利用误差精度量化指标对比训练集和测试集数据。通过这种方式,可以在未参与建模的数据集上测试模型。为尽力避免偶然性,每一层地层中所有控制点数据均参与验证,具体形式是将建模得到的具体每个地层中的深度值与原数据集中的测试数据进行对比,并使用CC、RMSE、bias和MAE^[39]4种指标进行量化衡量,如表2。表2显示,实验得到的研究区三维地质模型具有实际应用精度水平,且接近地表的年轻地层的建模精度明显高于较老地层,这是由于钻孔勘测实际工作中,受制于地质条件和钻工深度,在较深地层位置的控制点较少,导致经过插值等建模过程后的精度较差。后续工作中,为进一步控制三维地质模型的精度和质量,可考虑从数据入手,控制缺乏钻孔控制点的深部地层和断层构造等形态走向,比如结合现场调查数据、地球物理数据、虚拟钻孔控制数据与多源数据整合等。

3 讨论

3.1 建模方法的优势

随着隐式插值函数的相关算法不断改进,三维地质建模在助力城市地质精细化建模方面的贡献也更加深入。

建模方法需具有一定的通用性。但基于移动四面体算法的隐式地质面建模方法^[40]对钻孔数据的预处理并未充分考虑,可能会忽视钻孔勘测数据本身对模型精度的负面影响以及对实际地质情况的刻画不佳。而此次工作的整体建模流程采用势场法和泛克立格插值法作为隐式建模核心算法,兼具考虑了由于钻孔勘探技术有限、研究区复杂地质条件影响、钻探的施工目的等因素可能导致的钻孔采集数据和标准底层表分层不匹配的情况,通过钻孔源数据的初始处理、浅孔深部控制预处理等过程将工程所获的钻孔勘测数据转换为建模数据集,并进行建模数据集二维剖面可视化、三维地质建模。

基于GemPy开源平台加入了针对地层尖灭处理的建模方案。GemPy三维隐式地质建模的优势之一是自动化处理断层等地质构造,但在处理地层尖灭时存在问题。城市地下空间的三维地质建模主要服务于地下工程建设,工程地质层一般为第四系地层,其中各岩土层的尖灭现象较为普遍,需要单独处理,此次工作在GemPy核心建模算法中加入针对地层尖灭的建模方案,使得整套建模工程流程更适用于城市地质领域。

建模方法所得到的地质模型具有较好的精度。根据2.2小节中交叉验证的结果,表2显示通过该建模方法得到的研究区三维地质模型精度较好,CC、RMSE、bias和MAE 4种指标均显示对研究区各地层的建模达到了一定精度,且由于结构地质建模中采用了全局插值器可规避各地层之间的相互关系对建模效果的负面影响。总体而言,较深地层的钻孔控制点较少会导致相应建模精度,而拥有较多控制点的上层地层显示,研究区三维地质模型的精度较好。

概率地质建模可在一定程度上解决三维地质建模的不确定性。由于钻孔数据在数量及钻深尺度等参数上的局限性,三维地质建模的不确定性普遍存在,传统上无法量化。GemPy迭代至2.1版本之后,增添的概率地质建模功能采用概率分布采样得到的扰动输入数据参与构建地质模型集合,以此模拟不确定性。相较于Borghi等提出的随机建模方法^[41]对数据的严格要求(例如需要密集的钻孔或剖面),具有一定优势。

3.2 可深入开展的研究

基于GemPy的三维隐式势场插值地质建模方法具有完整的建模流程,可以有效生成完整的城市地下空间三维地质模型,但是仍存在值得深入研究的问题。例如,图13表明,提高体元的平面精度可以提高模型的可视化质量,然而这种改进计算模型耗时较长,如何在提高模型平滑程度和可视化质量的同时尽可能节省内存消耗是未来需要持续关注的关键突破口。另外,构建建模数据集时,对钻孔源数据中可能出现的其他特殊情况未做充分的预处理工作,需要在日后的工作中深入探讨对这些情形的预处理方案。

4 结论

当前地质大数据技术与城市地下空间开发的结合对三维隐式地质建模的需求日益增加,相较于显式建模而言,三维隐式地质建模基于隐函数插值空间中的离散点生成模型,具有模型实时动态更新、人工操作少、数据适用度强等优势。同时,GemPy作为开源平台集成了三维隐式势场建模方法,相较封装完整的商业软件而言,算法改进难度更小,应用灵活性更大。本文对钻孔勘测得到的原始数据首先进行初始处理,并采用结构地质建模中的泛克里格作为核心建模算法,充分利用开源GemPy平台的优势,针对城市地质中常见的地层尖灭构造进行单独的建模方案设计并加入核心建模架构中;对研究区进行城市地下空间的三维地质建模,并对建模数据集进行二维可视化且对模型进行网格模式下的三维可视化等。最后,对得到的三维地质模型进行精度验证,采用分层K折交叉验证方法将建模数据集分为验证集和测试集,采用CC、RMSE、bias和MAE 4种指标来衡量模型精度。结果显示模型精度符合实际应用精度水平,对地层尖灭等特殊构造的单独处理使模型内部细节进一步精细化,尤其对拥有较充足控制点的沉积地层而言,建模精度更好。后续加入虚拟控制钻孔、多源数据融合、地质规则约束等方案可再次提升模型精度质量。

本文介绍的建模方法基于建模工作流程和钻孔数据的提取,体现了GemPy三维隐式势场建模方法在算法效率、整体平滑程度等相较显式建模的优势,而且充分利用GemPy开源平台的优势,提出针对城市地质中普遍存在的地层尖灭构造的建模方案,能够有效地自动模拟城市地下空间地层模型,对城市地下空间开发利用工作可提供有力数据支撑与决策依据。此外,对地质模型的各种空间分析功能可以辅助城市土地利用设计和规划,提高经济效益和社会效益。

参考文献

[1] ELMQVIST T, ANDERSSON E, MCPHEARSON T, et al. Urbanization in and for the Anthropocene[J]. Urban Sustainability, 2021, 1(1): 6.

[2] 郝卉. 城市地下空间发展现状与前景[J]. 砖瓦, 2021(11): 85-86.

[3] XIE H, ZHANG Y, CHEN Y, et al. A case study of development and utilization of urban underground space in Shenzhen and the Guangdong-Hong Kong-Macao Greater Bay Area[J]. Tunnelling and Underground Space Technology, 2021, 107: 103651.

[4] OSWALD C, RINNER C, ROBINSON A. Applications of 3D printing in physical geography education and urban visualization[J]. Cartographica: The International Journal for Geographic Information and Geovisualization, 2019, 54(4): 278-287.

[5] OLIEROOK H K H, SCALZO R, KOHN D, et al. Bayesian geological and geophysical data fusion for the construction and uncertainty quantification of 3D geological models[J]. Geoscience Frontiers, 2021, 12(1): 479-493.

[6] 田玉培, 王俊, 刘文文. 城市地下空间三维地质建模的研究及实现[J]. 现代矿业, 2021, 37(8): 49-52, 55.

[7] JING H, HANHAN H, GUISEN Z, et al. 3D geological modelling of superficial deposits in Beijing City[J]. Geology in China, 2019, 46(2): 244-254.

[8] BIRCH C. New systems for geological modelling-black box or best practice?[J]. Journal of the Southern African Institute of Mining and Metallurgy, 2014, 114(12): 993-1000.

[9] COWAN E J, BEATSON R K, FRIGHT W R, et al. Rapid geological modelling[C]//Applied structural geology for mineral exploration and mining, international symposium, Kalgoorlie, Australia: Kalgoorlie. 2002, 9: 23-25.

[10] COWAN E J, BEATSON R K, ROSS H J, et al. Practical implicit geological modelling[C]//Proceedings of the 5th international mining geology conference. Bendigoe: Australasian Institute of Mining and Metallurgy Melbourne, 2003, 8: 89-99.

[11] VOLLGGER S A, CRUDEN A R, AILLERES L, et al. Regional dome evolution and its control on ore-grade distribution: insights from 3D implicit modelling of the Navachab gold deposit, Namibia[J]. Ore Geology Reviews, 2015, 69: 268-284.

[12] VON HARTEN J, DE LA VARGA M, HILLIER M, et al. Informed local smoothing in 3D implicit geological modeling[J]. Minerals, 2021, 11(11): 1281.

[13] 王权, 邹艳红. 基于轮廓线层间形态插值的三维地质隐式曲面重构[J]. 地质科技通报, 2023, 42(5): 293-300.

[14] 王权. 三维地质隐式建模形态不确定性分析与定量评价研究[D]. 长沙: 中南大学, 2022.

[15] ZHONG D, WANG L, LIN B I, et al. Implicit modeling of complex orebody with constraints of geological rules[J]. Transactions of Nonferrous Metals Society of China, 2019, 29(11): 2392-2399.

[16] FRANK T, TERTOIS A L, MALLET J L. 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[J]. Computers & Geosciences, 2007, 33(7): 932-943.

[17] CUOMO S, GALLETTI A, GIUNTA G, et al. Reconstruction of implicit curves and surfaces via RBF interpolation[J]. Applied Numerical Mathematics, 2017, 116: 157-171.

[18] SKALA V. RBF interpolation with CSRBF of large data sets[J]. Procedia Computer Science, 2017, 108: 2433-2437.

[19] MACÊDO I, GOIS J P, VELHO L. Hermite radial basis function simplicits[C]//Proceedings of computer graphics forum. Oxford: Blackwell Publishing Ltd, 2011, 30(1): 27-42.

[20] GUO J, WU L, ZHOU W, et al. Section-constrained local geological interface dynamic updating method based on the HRBF surface[J]. Journal of Structural Geology, 2018, 107: 64-72.

[21] FLEISHMAN S, COHEN-OR D, SILVA C T. Robust moving least-squares fitting with sharp features[J]. ACM Transactions on Graphics (TOG), 2005, 24(3): 544-552.

[22] CARR J C, BEATSON R K, CHERRIE J B, et al. Reconstruction and representation of 3D objects with radial basis functions[C]//Proceedings of the 28th annual conference on computer graphics and interactive techniques (SIGGRAPH’01). New York: Association for Computing Machinery, 2001: 67-76.

[23] STOCH B, ANTHONISSEN C J, MCCALL M J, et al. 3D implicit modeling of the Sishen Mine: new resolution of the geometry and origin of Fe mineralization[J]. Mineralium Deposita, 2018, 53(6): 835-853.

[24] JONES M W, CHEN M. A new approach to the construction of surfaces from contour data[C]//Proceedings of computer graphics forum. Edinburgh: Blackwell Science Ltd, 1994, 13(3): 75-84.

[25] 曹勇, 苏晓慧, 孙梦梦, 开源建模软件在水利工程BIM领域中的应用[J]. 水利规划与设计, 2021(8): 96-101, 126.

[26] 张志武. 三维矢量数据矿山模型算法研究及在Blender中的应用[D]. 武汉: 武汉科技大学, 2016.

[27] 赵勇, 许国, 卢鹏, 基于GemPy的隐式三维地质建模方法[J]. 人民长江, 2023, 54(10): 98-104.

[28] WU J, SUN B. Discontinuous mechanical analysis of manifold element strain of rock slope based on open source Gempy[C]//E3S Web of conferences. Hangzhou: EDP Sciences, 2021, 248: 03084.

[29] BEN-SHABAT Y, GOULD S. DeepFit: 3D surface fitting via neural network weighted least squares[C]//Computer Vision-ECCV 2020: 16th European conference. Glasgow: Springer International Publishing, 2020: 20-34.

[30] SCHAAF A, VARGA M, WELLMANN F, et al. Constraining stochastic 3-D structural geological models with topology information using approximate Bayesian computation in GemPy 2.1[J]. Geoscientific Model Development, 2021, 14(6): 3899-3913.

[31] DE LA VARGA M, SCHAAF A, WELLMANN F. GemPy 1.0: open-source stochastic geological modeling and inversion[J]. Geoscience Model Development, 2019, 12 (32): 1-32.

[32] CHEN G, ZHU J, QIANG M, et al. Three-dimensional site characterization with borehole data: a case study of Suzhou area[J]. Engineering Geology, 2018, 234: 65-82.

[33] FISCHER M, PROPPE C. Enhanced universal kriging for transformed input parameter spaces[J]. Probabilistic Engineering Mechanics, 2023, 74: 103486.

[34] 熊前进, 柴小婷, 万翔, 对数泛克立格法在东天山土屋—延东地区Cu、Au异常信息提取中的应用[J]. 资源环境与工程, 2021, 35(4): 536-540.

[35] LORENSEN W E, CLINE H E. Marching cubes: a high resolution 3D surface construction algorithm[J]. ACM SIGGRAPH Computer Graphics, 1987, 21(4): 163-169.

[36] SANTOLARIA P, VENDEVILLE B C, GRAVELEAU F, et al. Double evaporitic décollements: influence of pinch-out overlapping in experimental thrust wedges[J]. Journal of Structural Geology, 2015, 76: 35-51.

[37] 苏栋, 黄茂隆, 韩文龙, 考虑城市地质环境影响的深圳市地下空间开发适宜性评价[J]. 地学前缘, 2023, 30(4): 514-524.

[38] 蒋旭, 王婷婷, 穆静. 地下空间开发利用适宜性与资源量的应用研究[J]. 地下空间与工程学报, 2018, 14(5): 1145-1153.

[39] RATA M, DOUAOUI A, LARID M, et al. Comparison of geostatistical interpolation methods to map annual rainfall in the Chéliff watershed, Algeria[J]. Theoretical and Applied Climatology, 2020, 141: 1009-1024.

[40] GUO J, WANG X, WANG J, et al. Three-dimensional geological modeling and spatial analysis from geotechnical borehole data using an implicit surface and marching tetrahedra algorithm[J]. Engineering Geology, 2021, 284: 106047.

[41] BORGHI A, RENARD P, FOURNIER L, et al. Stochastic fracture generation accounting for the stratification orientation in a folded environment based on an implicit geological model[J]. Engineering Geology, 2015, 187: 135-142.

基金资助

国家重点研发计划项目“煤矿灾害融合监控与决策数字化关键技术装备及示范应用(2022YFC3000083)”