融合復(fù)雜地質(zhì)信息的地球物理模型建模技術(shù)

任政勇， 王祥， 湯井田， 吳啟紅

1 中南大學有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室， 長沙 410083 2 中南大學有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室， 長沙 410083 3 中南大學地球科學與信息物理學院， 長沙 410083 4 成都大學建筑與土木工程學院， 成都 610106

0 引言

網(wǎng)格單元是組成地球物理計算模型的基礎(chǔ).非結(jié)構(gòu)化四面體網(wǎng)格相比結(jié)構(gòu)化網(wǎng)格能夠更加精確地模擬地形起伏和復(fù)雜地質(zhì)體輪廓，因此在地球物理正反演方法中得到了廣泛的應(yīng)用(Lelièvre et al., 2012; Ren et al., 2013, 2017; Lu et al., 2016; 葉益信和李予國, 2016; Dunham et al., 2018; Liu et al., 2018; 底青云等, 2019).然而，由于基于四面體網(wǎng)格的數(shù)值計算方法對于模型網(wǎng)格質(zhì)量的過高要求，建立復(fù)雜地質(zhì)模型需要花費很大的時間成本，嚴重制約了該方法在野外實際勘探中的進一步應(yīng)用.

近十年來，在三維地學建模技術(shù)和網(wǎng)格剖分方法逐漸成熟的基礎(chǔ)上，面向地球物理數(shù)值計算的復(fù)雜四面體網(wǎng)格建模也變得更加高效.隨著三維地學建模技術(shù)從可視化到定量化的發(fā)展，地質(zhì)建模軟件可以依據(jù)多源數(shù)據(jù)快速建立實際勘探工區(qū)的地形、地下構(gòu)造及物性特征等多種模型，將其有效整合構(gòu)建成為實際工區(qū)的三維體積模型(Caumon et al., 2009, 2013).三維地學建模技術(shù)逐漸運用于成礦預(yù)測、城市地下空間等領(lǐng)域(Liu et al., 2012; Deng et al., 2016; Liu et al., 2016; Qin and Liu, 2018; Di et al., 2020; Guo et al., 2020a, 2020b; Jin et al., 2020).目前，三維建模軟件SKUA-Gocad?(Paradgim, 2014)能以地質(zhì)信息為約束建立高質(zhì)量復(fù)雜三維體積模型.軟件內(nèi)置的Delaunay三角網(wǎng)格剖分算法以及離散光滑插值算法(Discrete Smooth Interpolation，DSI)(Mallet, 1989)能夠在地質(zhì)信息的約束下生成逼近真實情況的高質(zhì)量曲面三角網(wǎng)格，從而可以描繪任意精細的地下地質(zhì)結(jié)構(gòu)幾何形態(tài)特征.

然而面對越來越精細化的三維地質(zhì)模型，如何將其轉(zhuǎn)化為高質(zhì)量地球物理四面體網(wǎng)格模型且最大程度保留地質(zhì)信息成為了當前的熱點問題.面向地球物理數(shù)值計算的三維模型在進行四面體網(wǎng)格剖分時一般要滿足三個要求：1)界面三角網(wǎng)格不能畸形，需保證足夠的質(zhì)量以使計算過程穩(wěn)定；2)測點、源點等處需局部網(wǎng)格加密；3) 保證計算區(qū)域中各個地質(zhì)體單元三角網(wǎng)格封閉且空間拓撲關(guān)系連續(xù)，以分配四面體單元物理屬性.為達到上述要求，需要對三維體積模型的界面網(wǎng)格不斷進行調(diào)整，才能生成用于數(shù)值模擬計算的高質(zhì)量四面體網(wǎng)格模型.Zehner(2011)最先提出了面向有限元計算的地下幾何模型半自動化構(gòu)建流程，模擬計算了地下水的流動以及儲層中的CO2儲存情況.Zehner等(2015)基于Gocad、Gmsh(Geuzaine and Remacle, 2009)、Tetgen(Si, 2015)等軟件給定了一套四面體網(wǎng)格邊界上線和點的控制調(diào)整方法，可以通過表面三角網(wǎng)重采樣大幅度提高網(wǎng)格剖分的質(zhì)量，并基于此框架實現(xiàn)了復(fù)雜模型的三維瞬變電磁正演模擬.緊接著Ansari等(2020)對其進行了進一步的修改與完善，利用RINGMesh(Pellerin et al., 2017)調(diào)整網(wǎng)格，利用建模工具Facetmodeller(Lelièvre et al., 2018)來添加有限元計算邊界，建立了加拿大Lalor富集硫礦的真實復(fù)雜三維地球物理模型并計算了該模型的三維大地電磁響應(yīng).但是上述建模方法中的界面網(wǎng)格調(diào)整的策略都需要像Gmsh和RingMesh這樣的外部工具來進行重采樣，網(wǎng)格需要在兩個軟件中來回轉(zhuǎn)換，建模效率較低，需要進一步改進.

本文在前人研究基礎(chǔ)上，提出了一種新的復(fù)雜地球物理模型建模技術(shù)，能夠基于地質(zhì)剖面圖等先驗信息，構(gòu)建包含大區(qū)域任意起伏地形，任意起伏地層和任意形態(tài)異常體的高質(zhì)量地球物理模型.利用SKUA-Gocad建模軟件完成三維體積模型的建立以及界面網(wǎng)格的調(diào)整，一定程度上提高了建模效率和模型質(zhì)量.而后利用本方法建立了安徽泥河鐵礦區(qū)復(fù)雜地球物理四面體模型，并基于前期推導(dǎo)的重磁解析解公式計算其重磁場響應(yīng)，計算結(jié)果與地質(zhì)背景基本吻合，說明本方法建立的地下復(fù)雜模型實用性和準確性較高，為解決大規(guī)模復(fù)雜實際地質(zhì)模型的反演問題提供了一定的技術(shù)基礎(chǔ).

1 建模方法

現(xiàn)有建模方法往往利用Gmsh和RingMesh等外部工具對交叉界面的三角網(wǎng)格進行重采樣來處理斷層以及巖體侵入等界面交叉問題.為了進一步提高建模效率，本文基于SKUA-Gocad三維建模平臺的內(nèi)置程序Surface Resample和Beautify Triangles實現(xiàn)了界面網(wǎng)格的約束重采樣，并以此提出了一種新的建模方法(圖1).該方法無需借助外部工具，不僅能夠快速解決界面交叉問題，且可以實現(xiàn)局部網(wǎng)格加密以及改善地形畸形網(wǎng)格等功能，提高了模型質(zhì)量.

圖1 基于地質(zhì)信息的復(fù)雜地球物理模型建模流程Fig.1 Modeling flow of complex geophysical model based on geological information

本文建模方法主要包含三維地質(zhì)建模以及非結(jié)構(gòu)四面體網(wǎng)格剖分兩個主體部分.首先，利用地理信息系統(tǒng)(GIS)軟件將地質(zhì)剖面圖等信息數(shù)字化，使其抽象為線串和點.再利用SKUA-Gocad和Facetmodeller建模平臺，以這些線串和點作為約束建立地形表面、地下地層單元、斷層以及礦體等要素模型.接著通過界面網(wǎng)格調(diào)整將這些地質(zhì)要素模型連接為一個綜合整體，從而構(gòu)成三維體積模型以進行四面體網(wǎng)格剖分.期間，為生成能夠用于地球物理數(shù)值計算的高質(zhì)量復(fù)雜四面體網(wǎng)格，需要不斷改善畸形網(wǎng)格和實施局部網(wǎng)格加密.最后利用Tetgen進行四面體網(wǎng)格剖分以及物性參數(shù)賦值.

1.1 構(gòu)建任意起伏形態(tài)地表模型

基于高分辨率數(shù)字高程數(shù)據(jù)(DEM)，利用SKUA-Gocad三維建模平臺實現(xiàn)了大區(qū)域的實際地形面構(gòu)建(圖2).首先，根據(jù)研究要求確定研究區(qū)域的網(wǎng)格單元大小，采用Delaunay三角剖分技術(shù)建立研究范圍平面網(wǎng)格，并對測點、源點等處進行局部網(wǎng)格加密.再通過SKUA-Gocad平臺中的Beautify Triangles工具改善畸形網(wǎng)格，提高網(wǎng)格質(zhì)量.接下來將研究區(qū)域DEM高程數(shù)據(jù)導(dǎo)入建模平臺，利用該平臺內(nèi)置的DSI插值算法，以DEM高程數(shù)據(jù)為約束對平面網(wǎng)格進行插值賦值，從而生成任意起伏形態(tài)的地形三角網(wǎng)格曲面(圖2d).

圖2 利用SKUA-Gocad建模平臺生成大規(guī)模實際地形表面(a) 基于Delaunay三角剖分建立平面網(wǎng)格(白色框內(nèi)為研究區(qū)域，黃色點為測點)； (b) 研究區(qū)域局部加密； (c) 網(wǎng)格形狀調(diào)整及測點處局部加密； (d) 基于DEM高程數(shù)據(jù)插值生成起伏地形.Fig.2 Generation of large-scale realistic topography surface based on SKUA-Gocad software(a) The surface grids generated by Delaunay triangulation (the white rectangle is the area of interesting and the yellow points are measure points); (b) Local refinement in the area of interesting; (c) Adjustment of mesh shape and local refinement at measure points; (d) Generation of undulating topography by interpolation of DEM data.

1.2 構(gòu)建任意地下地質(zhì)體模型

地下地質(zhì)體模型構(gòu)建包括地層單元、斷層、巖體及礦體的建立，其建模數(shù)據(jù)源主要是二維矢量剖面圖、地質(zhì)構(gòu)造圖、鉆孔數(shù)據(jù)以及物探解譯數(shù)據(jù)等.第一步，綜合所有已知數(shù)據(jù)，從地質(zhì)剖面圖中推斷出屬于同一地質(zhì)體的邊界線，將其構(gòu)建成為線串作為約束地質(zhì)體外表面范圍和形態(tài)的主要控制線.若不同類型數(shù)據(jù)對于地質(zhì)體邊界圈定出現(xiàn)差異，則以鉆孔等確定性較強的數(shù)據(jù)為主要參考確定地質(zhì)體邊界線.第二步，選取屬于同一地質(zhì)體的成面線串作為約束，利用DSI插值算法生成地層分界面(圖3a)，同時利用Facetmodeller軟件建立巖、礦體輪廓(圖3b).最后，利用軟件接口將Facetmodeller中的巖、礦體模型鑲嵌進SKUA-Gocad建立的地層單元當中，組成初步的三維體積模型.

圖3 建立地層分界面和礦體輪廓面(a) 基于控制線約束插值生成地層界面； (b) 基于控制點、線建立礦體輪廓.Fig.3 Generation of stratigraphic surface and deposit outline(a) Generation ofstratigraphic surface based on interpolation constrained by control strings; (b) Generation of deposit outline based on interpolation by control strings and points.

1.3 調(diào)整界面網(wǎng)格

為利用Tetgen剖分程序生成高質(zhì)量復(fù)雜四面體網(wǎng)格，要對初步三維體積模型的界面網(wǎng)格進行調(diào)整，包含改善畸形三角網(wǎng)格、進行局部網(wǎng)格加密、相交界面網(wǎng)格修正及添加有限元計算邊界等操作.修正相交界面的網(wǎng)格是為了保證三維體積模型中各個地質(zhì)體單元三角網(wǎng)格封閉且空間拓撲關(guān)系連續(xù)，從而能夠利用Tetgen進行四面體剖分同時給不同的地質(zhì)體單元賦屬性參數(shù).當存在斷層或者巖體侵入的復(fù)雜情況時(圖4a、d)，地層面、斷層面和巖體外表面需要分別獨立構(gòu)建，故導(dǎo)致二者相交處的三角網(wǎng)格空間拓撲關(guān)系是不連續(xù)的.因此，本文將基于SKUA-Gocad平臺中的內(nèi)置插件Surface Resample和Beautify Triangles，以相交界面的邊界線和公共交線為約束進行網(wǎng)格重采樣，從而使二者空間拓撲關(guān)系連續(xù)(圖4b、c、e、f).該界面網(wǎng)格調(diào)整方法不需要外部工具，適用性較廣，大幅度地提高了建模效率.

圖4 相交界面網(wǎng)格調(diào)整方法示意圖(a) 斷層面與地層面交叉； (b) 提取公共線、點； (c) 界面網(wǎng)格約束重采樣； (d) 巖體侵入地層面； (e) 提取公共線、點； (f) 界面網(wǎng)格約束重采樣.Fig.4 Schematic diagram of grid adjustment method for intersecting interface(a) The intersection of fault and Stratigraphic plane; (b) Extraction of common lines and points; (c) Resampling with interface mesh constraints; (d) The intrusion of rock mass; (e) Extraction of common lines and points; (f) Resampling with interface mesh constraints.

1.4 物性參數(shù)賦值及重磁解析解公式

經(jīng)反復(fù)調(diào)整界面網(wǎng)格并建立空間拓撲關(guān)系連續(xù)的封閉三維體積模型后，利用Tetgen剖分程序獲得三維體積模型的四面體網(wǎng)格以及密度-磁化率等物性參數(shù)賦值.為盡可能與實際保持一致，我們采用基于野外樣本測量統(tǒng)計的平均值作為各個地質(zhì)單元的密度-磁化率參數(shù)值.

最后，為了保證精度，我們利用前期推導(dǎo)的基于多面體單元(本文簡化為四面體單元)的重磁解析解公式(圖5)，進行數(shù)值計算(Ren et al.,2017a,2018)：

(1)

圖5 基于四面體單元的重磁解析解公式示意圖Fig.5 Schematic diagram of calculation of gravity and magnetic analytical solution formula based on tetrahedral element

(2)

(3)

(4)

2 應(yīng)用實例

2.1 已有地質(zhì)信息

泥河鐵礦床位于安徽省廬樅火山巖盆地西北邊緣，整體處在羅河—缺口斷裂東南側(cè)，北東向羅河—缺口成礦帶上，南西向距羅河鐵礦約4 km，于2006年被發(fā)現(xiàn).礦體隱伏于地下600 m的深處，共探明鐵礦石資源量18379.72萬t，TFe平均品位為29.69%.

礦區(qū)內(nèi)出露的地層由老到新為下白堊統(tǒng)磚橋組(K1z)和雙廟組(K1s)火山巖，楊灣組(K1y)砂巖及第四系(Q).賦礦地層為磚橋組(J3z)輝石粗安巖，閃長玢巖為侵入巖，主要礦石成分為磁鐵礦、黃鐵礦、硫鐵礦和硬石膏(圖6a)，礦體形態(tài)為環(huán)狀與透鏡狀，深度在-600 m到-1200 m之間(董樹文等, 2010; 范裕等, 2012; Lü et al., 2013).

2.2 建立測區(qū)地球物理模型

首先，確定研究范圍為泥河鐵礦區(qū)東西向X=531000 m-535000 m，南北向Y=3434800 m～3438000 m.緊接著建立研究區(qū)域地質(zhì)信息數(shù)據(jù)庫，包括測區(qū)30 m分辨率DEM數(shù)字高程數(shù)據(jù)，1 幅高分辨率平面地質(zhì)構(gòu)造圖，12 幅NW-SE向平行地質(zhì)剖面圖，1 幅NE-SW向垂直交叉地質(zhì)剖面圖，15個鉆孔點物性數(shù)據(jù)等.然后利用GIS軟件對上述數(shù)據(jù)進行矢量化、數(shù)字化，使其抽象成為一系列的數(shù)據(jù)點和線串，并作為先驗地質(zhì)信息實現(xiàn)對模型的約束(圖6a、b).

泥河鐵礦區(qū)地球物理模型的具體建立過程如圖7所示，按照建模流程，首先基于SKUA-Gocad平臺建立了高質(zhì)量的真實地形及地層界面幾何網(wǎng)格(圖7a、b、c)，再基于Facetmodeller建立磁鐵礦、硫鐵礦礦體輪廓幾何模型(圖7d).接著將其嵌入背景幾何模型，對邊界以及接觸部分進行控制與調(diào)整，整合成為封閉的泥河鐵礦三維體積模型.接下來采用開源軟件Tetgen對體積模型進行四面體網(wǎng)格剖分，賦予四面體單元密度、磁化率及電阻率等物性參數(shù)，最終生成12.39 km×10 km×10 km泥河鐵礦區(qū)復(fù)雜綜合地球物理模型(圖7e、f).

2.3 正演計算

為測試泥河鐵礦區(qū)大規(guī)模帶地形復(fù)雜地球物理模型的適用性，本文采用基于非結(jié)構(gòu)四面體網(wǎng)格的解析解公式計算該模型的重磁場響應(yīng).泥河鐵礦區(qū)計算模型大小為12.39 km×10 km×5 km(去除空氣層)，含四面體單元146422個，研究區(qū)域內(nèi)地表高程在10～53 m之間，磁鐵礦和硫鐵礦賦存在研究區(qū)域內(nèi)-600～-1200 m的深度，選取觀測平面位于礦體正上方，觀測點均勻分布在x=[531000 m, 535000 m]、y=[3434800 m,3438000 m]、z=100 m的范圍內(nèi)，線距100 m，點距100 m，共1353個測點(圖7b).

依據(jù)鉆孔巖芯以及野外采樣巖石的物性參數(shù)統(tǒng)計(楊波等, 2019; 周文月等, 2021)，將每個區(qū)域單元的平均密度值和平均磁化率值作為該區(qū)域內(nèi)的所有四面體單元的密度值和磁化率值，具體數(shù)值如表1所示.

表1 泥河鐵礦模型地層單元及礦體物性Table 1 Stratigraphic units and physical properties of Nihe Iron deposit

圖8、9和圖10、11為泥河鐵礦區(qū)復(fù)雜模型的高精度重磁場響應(yīng)計算結(jié)果，包括場值和梯度值響應(yīng).從圖8、圖10場值計算結(jié)果可以看出，垂向上的重力場gz和磁場Bz都呈現(xiàn)出明顯的局部異常.且由于磁鐵礦和硫鐵礦的磁化率差異相對密度差異更加顯著，使磁異常中心更靠近磁鐵礦.圖9、圖11中的梯度值計算結(jié)果顯示，垂向重磁場梯度都能較好的顯示出礦體的水平邊界，與地質(zhì)背景基本一致.為了探測該礦深部潛在的密度和磁性礦體，必須精確計算和減去復(fù)雜圍巖產(chǎn)生的重磁響應(yīng).總的來說，泥河鐵礦區(qū)復(fù)雜模型的高精度重磁場計算結(jié)果反映出本文所提出的復(fù)雜四面體網(wǎng)格建模技術(shù)具有良好的適用性和可行性.

圖8 基于解析解公式的泥河鐵礦區(qū)復(fù)雜模型重力異常場值計算結(jié)果紅色虛線為磁鐵礦分布范圍，黃色虛線為硫鐵礦分布范圍.Fig.8 Calculation result of gravity anomaly field of Nihe iron deposit based on analytical solution formulaThe red dash line is the distribution range of magnetite, and the yellow dash line is pyrite.

圖9 基于解析解公式的泥河鐵礦區(qū)復(fù)雜模型重力異常梯度值計算結(jié)果紅色虛線為磁鐵礦分布范圍，黃色虛線為硫鐵礦分布范圍.Fig.9 Calculation result of gravity anomaly gradient of Nihe iron deposit based on analytical solution formulaThe red dash line is the distribution range of magnetite, and the yellow dash line is pyrite.

圖10 基于解析解的泥河鐵礦區(qū)復(fù)雜模型磁感應(yīng)強度計算結(jié)果紅色虛線為磁鐵礦分布范圍，黃色虛線為硫鐵礦分布范圍.Fig.10 Calculation result of magnetic anomaly field of Nihe iron deposit based on analytical solution formulaThe red dash line is the distribution range of magnetite, and the yellow dash line is pyrite.

圖11 基于解析解的泥河鐵礦區(qū)復(fù)雜模型磁異常梯度值計算結(jié)果紅色虛線為磁鐵礦分布范圍，黃色虛線為硫鐵礦分布范圍.Fig.11 Calculation result of magnetic anomaly gradient of Nihe iron deposit based on analytical solution formulaThe red dash line is the distribution range of magnetite, and the yellow dash line is pyrite.

3 結(jié)論

(1) 本文提出了一種基于地質(zhì)信息的復(fù)雜四面體網(wǎng)格建模技術(shù)，能夠生成任意起伏地表、任意地層和任意地質(zhì)體的實際地球物理模型.該技術(shù)通過使用SKUA-Gocad內(nèi)置插件調(diào)節(jié)界面網(wǎng)格，提高了建模效率和質(zhì)量.最后基于該技術(shù)構(gòu)建了泥河鐵礦區(qū)大規(guī)模復(fù)雜實際模型，借助四面體網(wǎng)格解析解公式計算了該模型的高精度三維重磁場響應(yīng).

(2) 泥河鐵礦實際模型垂向方向的重力場及其梯度、磁場及其梯度計算結(jié)果能清晰地分辨出磁鐵礦和硫鐵礦的水平邊界，表明本文建模技術(shù)可靠性較高.這意味著利用已知信息可建立實際地球物理模型并精確計算其圍巖的重磁場響應(yīng)，接著將其從總體重磁響應(yīng)中去除，便可以探測該礦深部潛在的密度和磁性礦體.該研究成果進一步驗證了利用已知地質(zhì)信息建立高精度復(fù)雜地球物理模型的必要性.

(3) 本文的研究成果為復(fù)雜實際模型的反演提供了一定的技術(shù)基礎(chǔ).且文中所提出的建模流程同樣適用于直流電、電磁法以及地震方法中的非結(jié)構(gòu)四面體網(wǎng)格復(fù)雜模型構(gòu)建，具有廣泛的應(yīng)用前景.