基于GemPy的隱式三維地質(zhì)建模方法

趙 勇,許 國,盧 鵬,文 詩 寶,黃 梅 婷

(1.南寧市勘測設(shè)計(jì)院集團(tuán)有限公司,廣西 南寧 530022; 2.南寧市淺表地質(zhì)大數(shù)據(jù)工程技術(shù)研究中心,廣西 南寧 530022)

0 引 言

三維地質(zhì)模型不僅可以直觀展現(xiàn)地質(zhì)情況,還能輔助技術(shù)決策,在實(shí)際工作中發(fā)揮重要的作用[1]。創(chuàng)建模型的方法可分為顯式建模和隱式建模兩類。受傳統(tǒng)地礦勘察設(shè)計(jì)思路、方法及工作習(xí)慣的約束,目前三維地質(zhì)建模研究與應(yīng)用多集中于顯式建模[2];然而,傳統(tǒng)顯式建模工作量大、過程繁瑣、效率較低、模型更新重構(gòu)困難并存在拓?fù)溴e(cuò)誤、表面粗糙、棱角尖利等問題[3-4]。隱式建?；诤瘮?shù)插值空間離散點(diǎn)生成三維地質(zhì)模型,具有建模速度快、實(shí)時(shí)動(dòng)態(tài)更新、人工干預(yù)少、適應(yīng)數(shù)據(jù)能力強(qiáng)、結(jié)果精度高等特點(diǎn)[5],但隱式建模軟件費(fèi)用昂貴且閉源,限制了隱式建模的推廣和研究。國內(nèi)外學(xué)者提出了多種用于隱式建模的空間插值方法,如集成在GOCAD、Leapfrog Geo、GeoModeller等商業(yè)軟件中的離散光滑插值(DSI)[6]、快速徑向基插值(FastRBF)[7]、勢場理論插值[8]以及處于研究階段的HRBF插值[9]、GRBF插值[10]、楊赤中推估法[11]等,但提供這些插值方法的地質(zhì)建模軟件均未對公眾免費(fèi)開放。為減少對商業(yè)軟件的依賴,部分基于開源軟件Blender、FreeCAD的研究能實(shí)現(xiàn)礦山三維自動(dòng)建模及數(shù)字化應(yīng)用[12]、水閘工程三維參數(shù)化建模和有限元結(jié)構(gòu)計(jì)算[13],但軟件未融入地質(zhì)概念,用于地質(zhì)建模容易出現(xiàn)違背地質(zhì)規(guī)律的問題。少數(shù)學(xué)者利用Python研發(fā)地質(zhì)建模軟件,但因數(shù)據(jù)處理能力、算法優(yōu)化、建模方式等局限不能完成復(fù)雜地質(zhì)模型的創(chuàng)建[14]。

GemPy作為一款免費(fèi)開源的隱式三維地質(zhì)建模軟件,能創(chuàng)建褶皺、斷層、不整合等復(fù)雜地質(zhì)模型,其插值算法可與商業(yè)軟件媲美,并引入機(jī)器學(xué)習(xí)庫Theano使大規(guī)模矩陣運(yùn)算得以優(yōu)化,且考慮建模的不確定性,利用貝葉斯推理框架實(shí)現(xiàn)隨機(jī)地質(zhì)建模和貝葉斯反演[15]。國外學(xué)者利用GemPy強(qiáng)大的隱式建模能力輔助科學(xué)研究,如Fandel等[16]將GemPy與SKS(隨機(jī)巖溶模擬)、SWMM(暴雨洪水管理模型)集成,用于探索巖溶系統(tǒng)中管道網(wǎng)絡(luò)結(jié)構(gòu)與水力參數(shù)對地下水流動(dòng)特性的影響;Thomas等[17]考慮到隱式建模能快速重構(gòu)模型的特點(diǎn),用Petrel生成二維曲面并提取曲面點(diǎn)數(shù)據(jù)導(dǎo)入GemPy中創(chuàng)建不同分辨率的三維地質(zhì)模型;Jessell等[18]將GemPy作為三維地質(zhì)建模引擎,研發(fā)從數(shù)據(jù)提取到模型構(gòu)建全過程自動(dòng)化的三維地質(zhì)建模方法。國內(nèi)對GemPy的應(yīng)用與研究相對較少,孫丙陽[19]曾用GemPy創(chuàng)建三維地質(zhì)模型,研究了基于MOOSE的地質(zhì)模型體素化處理方法以及巖體工程開挖對連續(xù)-非連續(xù)變形的影響規(guī)律,但其采用的建模數(shù)據(jù)極少,且生成的模型較為簡單和理想化;王博等[20]僅將GemPy用于地層表面結(jié)構(gòu)數(shù)據(jù)計(jì)算。

雖然GemPy能創(chuàng)建復(fù)雜的三維地質(zhì)模型,且模型能滿足可視化及后續(xù)研究需求,有較強(qiáng)的實(shí)用性和拓展性,但目前仍多用于地層層序簡單的三維地質(zhì)模型創(chuàng)建;建模數(shù)據(jù)提取方法與數(shù)據(jù)結(jié)構(gòu)化是使用GemPy的關(guān)鍵,但前人研究論文中鮮有提及,且未能梳理出清晰的建模流程。為此,本文以廣西巖灘水電站為例,提出適用于GemPy的建模數(shù)據(jù)獲取方法和詳細(xì)可行的建模流程,創(chuàng)建受斷層和風(fēng)化剝蝕等地質(zhì)作用影響的復(fù)雜三維地質(zhì)模型,并分析其建模功能、原理、優(yōu)點(diǎn)及局限性,旨在為隱式三維地質(zhì)建模方法應(yīng)用與研究提供一些思考。

1 GemPy概述

GemPy是基于Python語言編寫的一款開源隱式三維建模軟件,它融合了地層層序、斷層、褶皺、不整合接觸等地質(zhì)要素,基于勢場理論,采用克里金插值方法,能構(gòu)造復(fù)雜的三維地質(zhì)模型。GemPy通過調(diào)用第三方庫來實(shí)現(xiàn)三維地質(zhì)模型的創(chuàng)建和展示。軟件框架如圖1所示,可概括為以下4個(gè)模塊。

圖1 GemPy軟件框架

(1) 數(shù)據(jù)讀取與導(dǎo)入模塊。GemPy將原始數(shù)據(jù)存儲(chǔ)和轉(zhuǎn)化為pandas數(shù)據(jù)結(jié)構(gòu)用于后續(xù)計(jì)算,將原始數(shù)據(jù)的創(chuàng)建、導(dǎo)出等操作封裝在類(Class)DateManagement中來初始化數(shù)據(jù),以減少計(jì)算量。

(2) 地層層序及構(gòu)造關(guān)系定義模塊。該模塊包含Series、Fault、Surface。Series用于定義地層層序、Fault用于指定斷層與地層的切割關(guān)系、Surface用于定義地層面。

(3) 插值運(yùn)算模塊。使用Theano優(yōu)化并高效計(jì)算出克里金參數(shù),使用scikit-image庫提供的移動(dòng)立方體算法進(jìn)行等值面提取和網(wǎng)格模型構(gòu)建。

(4) 成果展示和分析模塊。調(diào)用PyVista實(shí)現(xiàn)三維模型展示、調(diào)用Matplotlib實(shí)現(xiàn)二維成果展示,提供模型拓?fù)潢P(guān)系分析及根據(jù)地質(zhì)模型與指定密度進(jìn)行重力正演,使用PyMC3進(jìn)行模型貝葉斯反演。

2 建模理論

GemPy三維地質(zhì)建模應(yīng)用的基本原理是勢場理論[21-23]。該理論已成熟應(yīng)用于商業(yè)建模軟件GeoModeller。運(yùn)用此理論進(jìn)行地質(zhì)建模,有以下3種假設(shè):① 一個(gè)地質(zhì)界面劃分兩個(gè)地層;② 地層內(nèi)的方位數(shù)據(jù)與地質(zhì)界面有關(guān);③ 地質(zhì)界面是同方向無限組平行的曲面。

該理論的原理是在三維空間中創(chuàng)建任意點(diǎn)p=(x,y,z)的標(biāo)量函數(shù)T(p),使得所有T(p)=tk(tk為勢場的某個(gè)未知值)的點(diǎn)形成一個(gè)等勢面,將這個(gè)等勢面作為一個(gè)地質(zhì)分界面,則標(biāo)量場中的每一個(gè)等勢面將代表一個(gè)沉積界面,兩等勢面間形成一個(gè)地質(zhì)建造(地層),圖2為勢場原理示意。

圖2中,紅藍(lán)點(diǎn)為地層分界面上的采樣點(diǎn),箭頭為方位數(shù)據(jù)。采樣點(diǎn)通過的等勢線(紅藍(lán)線)為地層分界線,地層分界線之間的綠色等勢線可視為同一地層在不同時(shí)期的沉積界面(可認(rèn)為是地質(zhì)體趨勢或地層內(nèi)部層理軌跡),由紅藍(lán)地層分界線分隔出3個(gè)地層單元。從圖可觀察到方位數(shù)據(jù)可以是采樣點(diǎn)上的,也可以是其他位置的。通過插值采樣點(diǎn)和方位數(shù)據(jù)獲得地層分界面。

根據(jù)勢場理論,一個(gè)勢場被定義為一個(gè)獨(dú)立的地層序列(Series),在同一勢場內(nèi)等勢面之間相互平行,各獨(dú)立單元中的地質(zhì)模型只由原始數(shù)據(jù)及插值算法決定。然而現(xiàn)實(shí)地質(zhì)情況非常復(fù)雜,因沉積、侵蝕、侵入、中斷等地質(zhì)作用,地層之間相互切割、相互重疊,需要將勢場理論與地質(zhì)規(guī)律相結(jié)合來構(gòu)建復(fù)雜的三維地質(zhì)模型,可以通過定義多個(gè)勢場的方式解決。

基于協(xié)同克里金實(shí)現(xiàn)勢場理論插值算法,其函數(shù)表達(dá)方式如下:

式中:μα和νβ是p和p0的函數(shù),由協(xié)同克里金確定的權(quán)重值。pα和p′α是位于同一地質(zhì)界面上的2個(gè)采樣點(diǎn),M是增量T(pα)-T(p′α)的總個(gè)數(shù);?T(pβ)/?uβ是N個(gè)方位數(shù)據(jù)點(diǎn)pβ在任意方向的偏導(dǎo),代表標(biāo)量場的變化率。勢場理論插值算法實(shí)現(xiàn)模型重構(gòu)的過程可概述為將位于同一等值面上的采樣點(diǎn)及方位數(shù)據(jù)代入上式計(jì)算出權(quán)重系數(shù)后,反算出空間中未知點(diǎn)p與固定點(diǎn)po的距離增量。增量T*(p)-T*(p0)>0表示p點(diǎn)位于等值面外側(cè),T*(p)-T*(p0)<0表示p點(diǎn)位于等值面內(nèi)側(cè),T*(p)-T*(p0)=0表示位于等值面上。為了計(jì)算出增量為零的點(diǎn),采用Marching Cube移動(dòng)立方體算法計(jì)算出立方體8個(gè)頂點(diǎn)的數(shù)值從而提取等值面,最終實(shí)現(xiàn)模型重構(gòu)。

3 數(shù)據(jù)提取方法及建模流程

3.1 數(shù)據(jù)提取方法研究

基于上述理論,GemPy建模需要有已知的地質(zhì)界面采樣點(diǎn)(surfaces_point)和方位(orientation)兩類數(shù)據(jù)。界面點(diǎn)坐標(biāo)確定地質(zhì)界面的位置,方位數(shù)據(jù)用于約束地質(zhì)體的形狀和趨勢。界面點(diǎn)坐標(biāo)可以通過鉆孔直接獲取,而方位數(shù)據(jù)(包括傾向、傾角和極性)獲取比較困難,是使用GemPy建模的關(guān)鍵。方位數(shù)據(jù)最好是在有代表性、連續(xù)露頭的地質(zhì)界面處測量的傾向和傾角數(shù)據(jù),極性用±1表示,沉積關(guān)系的地層+1表示指向年代較晚的沉積方向,侵入體+1表示由外指向內(nèi)。

在方位數(shù)據(jù)的獲取方法方面,王博等[20]利用最小二乘法擬合出能反映全局產(chǎn)狀的地層平面并任取該平面內(nèi)三點(diǎn)坐標(biāo)來計(jì)算產(chǎn)狀;周文輝等[24]將產(chǎn)狀數(shù)據(jù)的求取轉(zhuǎn)換為求取剖面線所在近似平面和剖面線上任意處的切向量,然后借助過渡向量,將切向量轉(zhuǎn)換為法向量。受已有地質(zhì)資料所限,本文根據(jù)原始地質(zhì)剖面圖提取所選巖性段處的方位角和傾角,其中方位角為剖面線方向,傾角為巖性段處地層分界線切線方向與水平面的夾角?；陂_源軟件webplotDigitizer在二維剖面下提取離散點(diǎn)數(shù)據(jù),在Python環(huán)境下運(yùn)用SciPy庫將提取的離散點(diǎn)數(shù)據(jù)進(jìn)行三次樣條曲線擬合,再求取其一階導(dǎo)數(shù)作為切向量。計(jì)算切向量與Y軸的夾角后,將其轉(zhuǎn)換為法向與Y軸的角度。圖3為在Matplotlib中二維剖面切向量提取結(jié)果。

3.2 建模流程

由于GemPy缺少交互式操作界面,需要編寫Python代碼來創(chuàng)建三維地質(zhì)模型。運(yùn)行GemPy首先應(yīng)配置軟件所需的Python環(huán)境,然后借助PyCharm或Jupyter軟件來實(shí)現(xiàn)。建模的具體步驟如下。

(1) 數(shù)據(jù)準(zhǔn)備及初始化。GemPy按照預(yù)定的數(shù)據(jù)結(jié)構(gòu)進(jìn)行計(jì)算,因軟件無數(shù)據(jù)庫管理系統(tǒng),需要在外部以.csv文件按表1及表2結(jié)構(gòu)存儲(chǔ)原始數(shù)據(jù)。設(shè)置模型范圍及分辨率,由于軟件的限制,每個(gè)軸的分辨率設(shè)置值最大不超過100,即網(wǎng)格精度為模型各邊長度范圍值除以100。

表1 點(diǎn)坐標(biāo)數(shù)據(jù)表結(jié)構(gòu)

表2 方位數(shù)據(jù)表結(jié)構(gòu)

(2) 創(chuàng)建Grid網(wǎng)格。Grid網(wǎng)格為模型插值運(yùn)算提供空間坐標(biāo)數(shù)據(jù),是進(jìn)行空間分析、數(shù)值模擬等地質(zhì)應(yīng)用的基礎(chǔ),是屬性模型的載體[25]。在初始化階段選擇網(wǎng)格類型,GemPy提供5種網(wǎng)格類型(regular、topography、custom、section、centered)。regular為默認(rèn)網(wǎng)格,用于常規(guī)網(wǎng)格的建立;custom能設(shè)置任意形狀的網(wǎng)格;section用于二維剖面;centered用于計(jì)算物理屬性等。其中topography能利用基于GDAL(開源柵格空間數(shù)據(jù)轉(zhuǎn)換庫)的數(shù)據(jù)進(jìn)行地形網(wǎng)格創(chuàng)建。

(3) 地層序列的創(chuàng)建。GemPy在創(chuàng)建地層序列時(shí)遵循地層層序關(guān)系和斷層時(shí)序關(guān)系,使用函數(shù)gp.map_stack_to_surfaces( )定義地層的順序,地層新老關(guān)系按函數(shù)定義的先后順序確定。斷層需與地層同時(shí)設(shè)置且必須設(shè)置在地層序列的前面。對于地質(zhì)情況復(fù)雜的模型,可以定義多個(gè)序列(Series),由Series中不同巖層組合關(guān)系來處理剝蝕、侵入等地質(zhì)現(xiàn)象。

(4) 斷層面的創(chuàng)建。在上述地層層序中定義斷層面后,使用函數(shù)model.set_is_fault( )將斷層激活,激活后可設(shè)置斷層的接觸關(guān)系。GemPy對斷層的處理實(shí)質(zhì)是在克里金函數(shù)中添加一個(gè)偏移項(xiàng),斷層面創(chuàng)建所需要的數(shù)據(jù)與地層一樣,以斷層軌跡離散點(diǎn)作為點(diǎn)數(shù)據(jù)、方位數(shù)據(jù)作為梯度進(jìn)行插值。模型生成后,受斷層影響的地層面會(huì)貼合斷層面并產(chǎn)生急劇的彎曲來保證地層面的連續(xù)性,所以并未實(shí)質(zhì)性將地層切斷。

(5) 插值運(yùn)算。三維建模需要使用插值方法將離散的數(shù)據(jù)轉(zhuǎn)換為空間連續(xù)的數(shù)據(jù),基于克里金函數(shù)的插值是GemPy的核心,使用全局插值使同一勢場內(nèi)的兩個(gè)面不產(chǎn)生交叉。

(6) 二維、三維成果的展示。插值運(yùn)算完成后,調(diào)用Matplotlib查看二維剖面、調(diào)用PyVista查看三維模型。

(7) 拓?fù)浞治黾爸亓φ荨emPy內(nèi)置模型拓?fù)潢P(guān)系分析功能,使用函數(shù)gp.topology_compute( )可進(jìn)行拓?fù)浞治?。通過創(chuàng)建centered_grid,使用GravityPreprocessing( )函數(shù)可進(jìn)行重力的正演計(jì)算。GemPy建模流程如圖4所示。

圖4 GemPy建模流程

4 示 例

4.1 研究區(qū)地形地貌及巖土層情況

廣西巖灘水電站壩址工程所在地區(qū)為低山峽谷地形,兩岸高程在500～750 m之間,底部河流受區(qū)域構(gòu)造控制,河床狹窄,水流湍急。周圍地層巖性可概化為石炭系灰?guī)r、二疊系輝綠巖,其中二疊系輝綠巖按形成的年代由老到新分為5層(用1_1～1_5表示),場地地層被一條逆掩斷層切斷。

4.2 模型建立

(1) 數(shù)據(jù)準(zhǔn)備。本文所提出的產(chǎn)狀數(shù)據(jù)提取方法需要有已知的剖面圖且產(chǎn)狀數(shù)據(jù)獲取較為繁瑣,無法獲得未被剖面穿過區(qū)域的數(shù)據(jù),為驗(yàn)證GemPy在復(fù)雜地質(zhì)條件下的建模能力,示例采用的點(diǎn)數(shù)據(jù)及方位數(shù)據(jù)是在其他軟件平臺(tái)創(chuàng)建三維模型后提取。模型創(chuàng)建選取深度較深并能較完整反映場地地層情況的鉆孔,使用了點(diǎn)坐標(biāo)數(shù)據(jù)235條,方位數(shù)據(jù)139條。

(2) 初始化數(shù)據(jù)。按上述數(shù)據(jù)結(jié)構(gòu)存儲(chǔ)建模數(shù)據(jù),根據(jù)場區(qū)的范圍設(shè)置x,y,z3個(gè)方向的起止值,并將分辨率設(shè)置為軟件建議的最大值。采用默認(rèn)的regular網(wǎng)格來創(chuàng)建地質(zhì)模型,選擇topography網(wǎng)格并通過gdal讀取ASCII Grid(.asc)文件創(chuàng)建地形面,并預(yù)設(shè)2個(gè)剖面網(wǎng)格(section)。初始化數(shù)據(jù)后,可以通過PyVista在三維坐標(biāo)下查看數(shù)據(jù)分布情況(見圖5)。

圖5 數(shù)據(jù)空間分布

(3) 創(chuàng)建地層層序。根據(jù)場地地層及構(gòu)造情況將模型的地層分為3個(gè)序列,Series1為斷層上盤(5層輝綠巖)、Series2為斷層下盤(5層輝綠巖)、Basement為軟件默認(rèn)的基巖面(石炭系灰?guī)r),并將斷層上下盤的巖性更改為相同的顏色,如表3所列。

表3 地層層序配置

(4) 地層接觸關(guān)系處理。由研究區(qū)地質(zhì)資料可知,該地區(qū)地層被逆掩斷層錯(cuò)斷后,上盤地層又經(jīng)受了風(fēng)化剝蝕,所以底部巖層在地表出露。若采用GemPy提供的斷層功能,無法重構(gòu)巖層被剝蝕出露的狀態(tài)。受限于軟件的功能,從建模效果考慮并未采取讓斷層切割地層的方式,而是單獨(dú)將地層分為上盤和下盤兩個(gè)地層序列(Series)進(jìn)行建模。

(5) 插值運(yùn)算。上述步驟完成后,對數(shù)據(jù)進(jìn)行克里金插值和模型計(jì)算,最終運(yùn)行結(jié)果參數(shù)如表4所列。

表4 模型計(jì)算參數(shù)

(6) 模型展示。GemPy將PyVista內(nèi)置在程序內(nèi)部,通過函數(shù)gp.plot_3d( )可以查看創(chuàng)建的三維模型。PyVista提供簡單的三維操作功能,用戶可以進(jìn)行模型的旋轉(zhuǎn)、地層面隱藏與顯示等。圖6～8分別是通過GemPy生成的三維地形面模型、三維地層模型與三維地質(zhì)模型。

圖6 三維地形面模型

調(diào)用第三方庫ipywidgets,按單元來展示地質(zhì)模型在剖面上的地層分布,圖9～10分別為沿Y軸與沿X軸方向單元格(cell)為40,50,60位置處的地層剖切。

圖9 沿Y軸不同位置處的剖切面

圖10 沿X軸不同位置處的剖切面

基于勢場理論,GemPy利用點(diǎn)坐標(biāo)與方位數(shù)據(jù)通過克里金函數(shù)自動(dòng)插值生成地形面與地層面,進(jìn)而形成三維地質(zhì)模型。建模過程無需過多的人工交互,建模效率和自動(dòng)化程度較高。從實(shí)際建模效果來看,創(chuàng)建的地層曲面光滑,可視化效果好;創(chuàng)建的地質(zhì)模型與實(shí)際地質(zhì)情況相似,并能反映出逆掩斷層、地層剝蝕出露等復(fù)雜的地質(zhì)現(xiàn)象;創(chuàng)建的剖切面可以觀察不同方位巖層的空間分布情況。

5 結(jié) 語

通過對GemPy的功能、建模原理進(jìn)行分析并結(jié)合實(shí)際工程進(jìn)行三維地質(zhì)建模,總結(jié)分析出GemPy的優(yōu)缺點(diǎn)如下。

(1) 能實(shí)現(xiàn)復(fù)雜地質(zhì)情況的三維地質(zhì)模型及任意方向的二維剖切面創(chuàng)建與可視化。因軟件采用機(jī)器學(xué)習(xí)及貝葉斯推理,使創(chuàng)建的模型可信度較高,能反映實(shí)際地質(zhì)情況。

(2) 建模數(shù)據(jù)獲取困難,但軟件并沒有提供解決方案。通過在外部提取已有地質(zhì)資料的方法,效率很低,精度也不夠,如何獲取建模數(shù)據(jù)成為使用GemPy建模的關(guān)鍵。模型分辨率低,網(wǎng)格分辨率受限,軟件不能設(shè)置精細(xì)網(wǎng)格。

(3) 無交互式界面,需要用Python代碼實(shí)現(xiàn)交互功能,導(dǎo)致建模操作不夠便捷。沒有獨(dú)立的數(shù)據(jù)庫管理系統(tǒng),數(shù)據(jù)的存儲(chǔ)與調(diào)用不方便。

(4) 軟件算法需進(jìn)一步優(yōu)化,由于三維地質(zhì)建模處理的數(shù)據(jù)量非常大,GemPy在進(jìn)行克里金插值運(yùn)算時(shí)需占用較大內(nèi)存,無法制作分辨率較精細(xì)的模型。同時(shí),利用Python語言編寫的軟件進(jìn)行插值運(yùn)算非常耗時(shí),運(yùn)算效率較低。