融合復雜地質信息的地球物理模型建模技術

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

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

0 引言

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

近十年來，在三維地學建模技術和網格剖分方法逐漸成熟的基礎上，面向地球物理數值計算的復雜四面體網格建模也變得更加高效.隨著三維地學建模技術從可視化到定量化的發展，地質建模軟件可以依據多源數據快速建立實際勘探工區的地形、地下構造及物性特征等多種模型，將其有效整合構建成為實際工區的三維體積模型(Caumon et al., 2009, 2013).三維地學建模技術逐漸運用于成礦預測、城市地下空間等領域(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)能以地質信息為約束建立高質量復雜三維體積模型.軟件內置的Delaunay三角網格剖分算法以及離散光滑插值算法(Discrete Smooth Interpolation，DSI)(Mallet, 1989)能夠在地質信息的約束下生成逼近真實情況的高質量曲面三角網格，從而可以描繪任意精細的地下地質結構幾何形態特征.

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

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

1 建模方法

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

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

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

1.1 構建任意起伏形態地表模型

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

圖2 利用SKUA-Gocad建模平臺生成大規模實際地形表面(a) 基于Delaunay三角剖分建立平面網格(白色框內為研究區域，黃色點為測點)； (b) 研究區域局部加密； (c) 網格形狀調整及測點處局部加密； (d) 基于DEM高程數據插值生成起伏地形.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 構建任意地下地質體模型

地下地質體模型構建包括地層單元、斷層、巖體及礦體的建立，其建模數據源主要是二維矢量剖面圖、地質構造圖、鉆孔數據以及物探解譯數據等.第一步，綜合所有已知數據，從地質剖面圖中推斷出屬于同一地質體的邊界線，將其構建成為線串作為約束地質體外表面范圍和形態的主要控制線.若不同類型數據對于地質體邊界圈定出現差異，則以鉆孔等確定性較強的數據為主要參考確定地質體邊界線.第二步，選取屬于同一地質體的成面線串作為約束，利用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 調整界面網格

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

圖4 相交界面網格調整方法示意圖(a) 斷層面與地層面交叉； (b) 提取公共線、點； (c) 界面網格約束重采樣； (d) 巖體侵入地層面； (e) 提取公共線、點； (f) 界面網格約束重采樣.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 物性參數賦值及重磁解析解公式

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

最后，為了保證精度，我們利用前期推導的基于多面體單元(本文簡化為四面體單元)的重磁解析解公式(圖5)，進行數值計算(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 應用實例

2.1 已有地質信息

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

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

2.2 建立測區地球物理模型

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

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

2.3 正演計算

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

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

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

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

圖8 基于解析解公式的泥河鐵礦區復雜模型重力異常場值計算結果紅色虛線為磁鐵礦分布范圍，黃色虛線為硫鐵礦分布范圍.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 基于解析解公式的泥河鐵礦區復雜模型重力異常梯度值計算結果紅色虛線為磁鐵礦分布范圍，黃色虛線為硫鐵礦分布范圍.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 基于解析解的泥河鐵礦區復雜模型磁感應強度計算結果紅色虛線為磁鐵礦分布范圍，黃色虛線為硫鐵礦分布范圍.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 基于解析解的泥河鐵礦區復雜模型磁異常梯度值計算結果紅色虛線為磁鐵礦分布范圍，黃色虛線為硫鐵礦分布范圍.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 結論

(1) 本文提出了一種基于地質信息的復雜四面體網格建模技術，能夠生成任意起伏地表、任意地層和任意地質體的實際地球物理模型.該技術通過使用SKUA-Gocad內置插件調節界面網格，提高了建模效率和質量.最后基于該技術構建了泥河鐵礦區大規模復雜實際模型，借助四面體網格解析解公式計算了該模型的高精度三維重磁場響應.

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

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