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

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

0 引言

網(wǎng)格單元是組成地球物理計(jì)算模型的基礎(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; 葉益信和李予國(guó), 2016; Dunham et al., 2018; Liu et al., 2018; 底青云等, 2019).然而，由于基于四面體網(wǎng)格的數(shù)值計(jì)算方法對(duì)于模型網(wǎng)格質(zhì)量的過(guò)高要求，建立復(fù)雜地質(zhì)模型需要花費(fèi)很大的時(shí)間成本，嚴(yán)重制約了該方法在野外實(shí)際勘探中的進(jìn)一步應(yīng)用.

近十年來(lái)，在三維地學(xué)建模技術(shù)和網(wǎng)格剖分方法逐漸成熟的基礎(chǔ)上，面向地球物理數(shù)值計(jì)算的復(fù)雜四面體網(wǎng)格建模也變得更加高效.隨著三維地學(xué)建模技術(shù)從可視化到定量化的發(fā)展，地質(zhì)建模軟件可以依據(jù)多源數(shù)據(jù)快速建立實(shí)際勘探工區(qū)的地形、地下構(gòu)造及物性特征等多種模型，將其有效整合構(gòu)建成為實(shí)際工區(qū)的三維體積模型(Caumon et al., 2009, 2013).三維地學(xué)建模技術(shù)逐漸運(yùn)用于成礦預(yù)測(cè)、城市地下空間等領(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ì)信息的約束下生成逼近真實(shí)情況的高質(zhì)量曲面三角網(wǎng)格，從而可以描繪任意精細(xì)的地下地質(zhì)結(jié)構(gòu)幾何形態(tài)特征.

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

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

1 建模方法

現(xiàn)有建模方法往往利用Gmsh和RingMesh等外部工具對(duì)交叉界面的三角網(wǎng)格進(jìn)行重采樣來(lái)處理斷層以及巖體侵入等界面交叉問(wèn)題.為了進(jìn)一步提高建模效率，本文基于SKUA-Gocad三維建模平臺(tái)的內(nèi)置程序Surface Resample和Beautify Triangles實(shí)現(xiàn)了界面網(wǎng)格的約束重采樣，并以此提出了一種新的建模方法(圖1).該方法無(wú)需借助外部工具，不僅能夠快速解決界面交叉問(wèn)題，且可以實(shí)現(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)格剖分兩個(gè)主體部分.首先，利用地理信息系統(tǒng)(GIS)軟件將地質(zhì)剖面圖等信息數(shù)字化，使其抽象為線(xiàn)串和點(diǎn).再利用SKUA-Gocad和Facetmodeller建模平臺(tái)，以這些線(xiàn)串和點(diǎn)作為約束建立地形表面、地下地層單元、斷層以及礦體等要素模型.接著通過(guò)界面網(wǎng)格調(diào)整將這些地質(zhì)要素模型連接為一個(gè)綜合整體，從而構(gòu)成三維體積模型以進(jìn)行四面體網(wǎng)格剖分.期間，為生成能夠用于地球物理數(shù)值計(jì)算的高質(zhì)量復(fù)雜四面體網(wǎng)格，需要不斷改善畸形網(wǎng)格和實(shí)施局部網(wǎng)格加密.最后利用Tetgen進(jìn)行四面體網(wǎng)格剖分以及物性參數(shù)賦值.

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

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

圖2 利用SKUA-Gocad建模平臺(tái)生成大規(guī)模實(shí)際地形表面(a) 基于Delaunay三角剖分建立平面網(wǎng)格(白色框內(nèi)為研究區(qū)域，黃色點(diǎn)為測(cè)點(diǎn))； (b) 研究區(qū)域局部加密； (c) 網(wǎng)格形狀調(diào)整及測(cè)點(diǎn)處局部加密； (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ì)體的邊界線(xiàn)，將其構(gòu)建成為線(xiàn)串作為約束地質(zhì)體外表面范圍和形態(tài)的主要控制線(xiàn).若不同類(lèi)型數(shù)據(jù)對(duì)于地質(zhì)體邊界圈定出現(xiàn)差異，則以鉆孔等確定性較強(qiáng)的數(shù)據(jù)為主要參考確定地質(zhì)體邊界線(xiàn).第二步，選取屬于同一地質(zhì)體的成面線(xiàn)串作為約束，利用DSI插值算法生成地層分界面(圖3a)，同時(shí)利用Facetmodeller軟件建立巖、礦體輪廓(圖3b).最后，利用軟件接口將Facetmodeller中的巖、礦體模型鑲嵌進(jìn)SKUA-Gocad建立的地層單元當(dāng)中，組成初步的三維體積模型.

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

圖4 相交界面網(wǎng)格調(diào)整方法示意圖(a) 斷層面與地層面交叉； (b) 提取公共線(xiàn)、點(diǎn)； (c) 界面網(wǎng)格約束重采樣； (d) 巖體侵入地層面； (e) 提取公共線(xiàn)、點(diǎn)； (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)格并建立空間拓?fù)潢P(guān)系連續(xù)的封閉三維體積模型后，利用Tetgen剖分程序獲得三維體積模型的四面體網(wǎng)格以及密度-磁化率等物性參數(shù)賦值.為盡可能與實(shí)際保持一致，我們采用基于野外樣本測(cè)量統(tǒng)計(jì)的平均值作為各個(gè)地質(zhì)單元的密度-磁化率參數(shù)值.

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

2.1 已有地質(zhì)信息

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

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

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

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

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

2.3 正演計(jì)算

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

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

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

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

圖8 基于解析解公式的泥河鐵礦區(qū)復(fù)雜模型重力異常場(chǎng)值計(jì)算結(jié)果紅色虛線(xiàn)為磁鐵礦分布范圍，黃色虛線(xiàn)為硫鐵礦分布范圍.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ù)雜模型重力異常梯度值計(jì)算結(jié)果紅色虛線(xiàn)為磁鐵礦分布范圍，黃色虛線(xiàn)為硫鐵礦分布范圍.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)強(qiáng)度計(jì)算結(jié)果紅色虛線(xiàn)為磁鐵礦分布范圍，黃色虛線(xiàn)為硫鐵礦分布范圍.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ù)雜模型磁異常梯度值計(jì)算結(jié)果紅色虛線(xiàn)為磁鐵礦分布范圍，黃色虛線(xiàn)為硫鐵礦分布范圍.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í)際地球物理模型.該技術(shù)通過(guò)使用SKUA-Gocad內(nèi)置插件調(diào)節(jié)界面網(wǎng)格，提高了建模效率和質(zhì)量.最后基于該技術(shù)構(gòu)建了泥河鐵礦區(qū)大規(guī)模復(fù)雜實(shí)際模型，借助四面體網(wǎng)格解析解公式計(jì)算了該模型的高精度三維重磁場(chǎng)響應(yīng).

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

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