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

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

0 引 言

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

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

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

1 GemPy概述

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

圖1 GemPy軟件框架

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

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

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

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

2 建模理論

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

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

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

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

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

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

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

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

基于上述理論,GemPy建模需要有已知的地質(zhì)界面采樣點(surfaces_point)和方位(orientation)兩類數(shù)據(jù)。界面點坐標確定地質(zhì)界面的位置,方位數(shù)據(jù)用于約束地質(zhì)體的形狀和趨勢。界面點坐標可以通過鉆孔直接獲取,而方位數(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)三點坐標來計算產(chǎn)狀;周文輝等[24]將產(chǎn)狀數(shù)據(jù)的求取轉(zhuǎn)換為求取剖面線所在近似平面和剖面線上任意處的切向量,然后借助過渡向量,將切向量轉(zhuǎn)換為法向量。受已有地質(zhì)資料所限,本文根據(jù)原始地質(zhì)剖面圖提取所選巖性段處的方位角和傾角,其中方位角為剖面線方向,傾角為巖性段處地層分界線切線方向與水平面的夾角?；陂_源軟件webplotDigitizer在二維剖面下提取離散點數(shù)據(jù),在Python環(huán)境下運用SciPy庫將提取的離散點數(shù)據(jù)進行三次樣條曲線擬合,再求取其一階導(dǎo)數(shù)作為切向量。計算切向量與Y軸的夾角后,將其轉(zhuǎn)換為法向與Y軸的角度。圖3為在Matplotlib中二維剖面切向量提取結(jié)果。

3.2 建模流程

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

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

表1 點坐標數(shù)據(jù)表結(jié)構(gòu)

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

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

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

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

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

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

(7) 拓撲分析及重力正演。GemPy內(nèi)置模型拓撲關(guān)系分析功能,使用函數(shù)gp.topology_compute( )可進行拓撲分析。通過創(chuàng)建centered_grid,使用GravityPreprocessing( )函數(shù)可進行重力的正演計算。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ù)準備。本文所提出的產(chǎn)狀數(shù)據(jù)提取方法需要有已知的剖面圖且產(chǎn)狀數(shù)據(jù)獲取較為繁瑣,無法獲得未被剖面穿過區(qū)域的數(shù)據(jù),為驗證GemPy在復(fù)雜地質(zhì)條件下的建模能力,示例采用的點數(shù)據(jù)及方位數(shù)據(jù)是在其他軟件平臺創(chuàng)建三維模型后提取。模型創(chuàng)建選取深度較深并能較完整反映場地地層情況的鉆孔,使用了點坐標數(shù)據(jù)235條,方位數(shù)據(jù)139條。

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

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

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

表3 地層層序配置

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

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

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

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

圖6 三維地形面模型

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

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

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

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

5 結(jié) 語

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

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

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

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

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