李 楠 宋相龍 楚文楷

(1. 中國地質(zhì)科學(xué)院礦產(chǎn)資源研究所，北京 100037； 2. 中國地質(zhì)大學(xué)(北京)研究生院，北京 100083； 3. 中國地質(zhì)科學(xué)院，北京 100037)

0 引 言

20 世紀(jì)80 年代以來，尋找隱伏礦和深部礦已成為許多國家的主要找礦方向，地質(zhì)找礦難度日漸增大，找礦效果日趨降低(陳建平等，2002，2014a)。隨著礦產(chǎn)資源勘查深度的逐漸加深，傳統(tǒng)的基于二維平面開展的礦產(chǎn)資源預(yù)測評價已不適應(yīng)要求，如何更有效地找礦成為焦點問題，為此各國一直致力于研究新的成礦預(yù)測方法。隨著計算機圖形處理、三維可視化及三維插值技術(shù)的迅速發(fā)展，利用計算機三維地質(zhì)建模進行三維成礦定量預(yù)測的研究逐漸展開，成為礦產(chǎn)資源預(yù)測的一大亮點(陳建平等，2014b；劉靜等，2017) 。

三維地質(zhì)建模是對地學(xué)多元數(shù)據(jù)集進行綜合解譯與可視化表達的過程(肖克炎等，2000)。僅憑有限的鉆孔信息得到的地下三維地質(zhì)情況往往只是“一孔之見”。根據(jù)對研究區(qū)地質(zhì)情況的了解，形成地質(zhì)概念模型，然后在不同剖面上進行地質(zhì)體的圈定，不同的專家可能有不同的地質(zhì)體建模方法。

三維地質(zhì)建模是覆蓋區(qū)與深部找礦預(yù)測的前提和基礎(chǔ)。通過廣泛收集、整理、轉(zhuǎn)換研究區(qū)數(shù)據(jù)，最終形成各類地學(xué)信息在三維空間上的分布模型(肖克炎等，2012)。筆者研究了基于區(qū)域大比例尺且缺少鉆孔信息情況下的一種隱式建模流程，有大量的地表實測數(shù)據(jù)但缺乏足夠的深部數(shù)據(jù)(鉆孔數(shù)據(jù))支持?；诘乇頂?shù)據(jù)，結(jié)合地球物理數(shù)據(jù)，通過三維反演、地質(zhì)統(tǒng)計學(xué)、幾何算法與三維可視化等技術(shù)方法，直觀地刻畫出地下一定深度范圍內(nèi)的構(gòu)造、巖體、地層等的空間展布以及成因和演化關(guān)系，為下一步找礦方向提供指導(dǎo)。

1 區(qū)域大比例尺隱式建模技術(shù)

1.1 三維地質(zhì)建模的尺度

區(qū)域三維地質(zhì)建模是一個多元地學(xué)信息綜合解釋的過程。一個能夠真實反映深部地質(zhì)情況的三維地質(zhì)模型應(yīng)由遵循地質(zhì)規(guī)律的地質(zhì)體及其空間關(guān)系(幾何拓撲)構(gòu)成(Xiao et al., 2015)，在此方面已經(jīng)開展了大量研究工作。根據(jù)研究區(qū)范圍(比例尺與深度)、原始建模數(shù)據(jù)以及使用的建模方法等的不同，將三維地質(zhì)建模分為礦區(qū)、區(qū)域和巖石圈3個不同尺度類型(表1)。從建模數(shù)據(jù)分析，礦區(qū)尺度三維地質(zhì)建模擁有大量地表、鉆孔等數(shù)據(jù)，區(qū)域大比例尺擁有大量地表實測數(shù)據(jù)但缺乏足夠的深部數(shù)據(jù)，巖石圈建模完全依賴于深部地球物理探測技術(shù)手段獲取的數(shù)據(jù)，因此針對上述3個尺度的建模方法存在明顯差異。從對地質(zhì)模型的質(zhì)量要求看，由于受到地球深部探測技術(shù)手段、經(jīng)濟效益和對于深部地質(zhì)演化過程理解的限制，在符合現(xiàn)有地質(zhì)理論的前提下，隨著建模區(qū)域范圍的擴大和深部已知數(shù)據(jù)的減少，對模型精度的要求一般呈降低的趨勢。根據(jù)地質(zhì)找礦理論，不同尺度所面臨的研究任務(wù)也不盡相同。

表1 三維地質(zhì)建模研究的3個尺度

國際上在上述3個尺度均開展了大量的研究工作。影響最大的是澳大利亞聯(lián)合科學(xué)與工業(yè)研究組織(CSIRO)提出的“玻璃地球”計劃，目的是利用三維可視化和地質(zhì)建模等技術(shù)，使從地表到地下1 km的大陸表面“像玻璃一樣透明”，同時通過三維可視化和廣泛的網(wǎng)絡(luò)服務(wù)接口傳播信息填制澳大利亞信息豐富的四維“地圖”，該項工作的主要目標(biāo)之一是服務(wù)于澳大利亞下一代巨型礦床的發(fā)現(xiàn)(陳應(yīng)軍等，2014)。歐美各發(fā)達國家率先響應(yīng)并制定了實施規(guī)劃，開展了一系列研究與實踐：美國計劃在20世紀(jì)90年代繪制的三維地質(zhì)圖的基礎(chǔ)上開展全境范圍的5 km以淺的地質(zhì)填圖工作，主要技術(shù)手段為高分辨率三維地震成像等；英國計劃構(gòu)建全境多比例尺(1∶100萬—1∶1萬)的三維地層模型；俄羅斯、法國、荷蘭、德國等也有類似的三維地質(zhì)建模工作計劃。

1.2 三維地質(zhì)建模的算法

為實現(xiàn)上述三維地質(zhì)建模工作，國外研究人員提出了大量的計算機算法，用于構(gòu)建三維地質(zhì)模型。根據(jù)輸入數(shù)據(jù)與建模過程的不同，算法主要分為顯式建模與隱式建模2類。

1.2.1 顯式建模算法 指在已知點集的基礎(chǔ)上, 通過數(shù)學(xué)函數(shù)或人機交互的方式重構(gòu)地質(zhì)體，即數(shù)據(jù)驅(qū)動型算法(李楠等，2015)，其中已知點集為地質(zhì)體曲面(重構(gòu)曲面)上的全部控制點。顯式建模的算法主要包括3大類。

(1) 應(yīng)用空間Delaunay三角剖分算法重構(gòu)地質(zhì)對象。Cheng等(2000)首先提出了一種基于離散點集凸殼的方法, 其基本思想是基于點集構(gòu)建四面體凸殼，稱作α-shape，并應(yīng)用Delaunay剖分在凸殼的約束下重構(gòu)點集所在的曲面。之后又提出了一些改進算法(Bernardini et al.,1999; Gerritsen et al.,2002)。

(2) 樣條函數(shù)擬合法。Hoppe等(1992)和Grimm等(1995)提出基于B樣條的曲面重構(gòu)算法，Park等(1995)提出基于Bézier曲面的近似曲面重構(gòu)算法。該方法的最大缺點是無法強制要求結(jié)果曲面通過已知地質(zhì)控制點，難以構(gòu)建不連續(xù)或者發(fā)生合理突變等的復(fù)雜地質(zhì)曲面, 例如侵入巖體、斷層網(wǎng)絡(luò)或復(fù)雜的接觸界面等地質(zhì)情況。

(3) 結(jié)合人機交互方式構(gòu)建更加精確的地質(zhì)曲面。最具有代表性的是輪廓線拼接法和離散光滑插值法(Discrete Smooth Interpolation，DSI)：輪廓線拼接法(Tipper,1976)是最為成熟也是應(yīng)用最為廣泛的三維曲面重構(gòu)算法之一，應(yīng)用實例如Surpac、Micromine等礦山勘探軟件，核心思想是將相鄰地質(zhì)剖面上相同屬性的點對通過人機交互的方式連接起來，最終構(gòu)建地質(zhì)體表面模型；Mallet(1992,1997)提出了基于人機交互、離散點空間Delaunay剖分與摩擦函數(shù)(Roughness)的曲面重構(gòu)算法(MLPSR)。通過摩擦函數(shù)，DSI算法使擬合曲面S′通過合理的數(shù)學(xué)模型無限接近于原始曲面S，解決了結(jié)果曲面S′生成時遇到的控制點集(不僅僅是近似)、曲面光滑和不連續(xù)曲面的重構(gòu)等復(fù)雜問題，該算法在區(qū)域地質(zhì)建模領(lǐng)域得到了廣泛的應(yīng)用并推出了商業(yè)化軟件GOCAD。

1.2.2 隱式建模算法 與顯式建模算法相比，隱式建模算法的最大不同是使用了空間插值和等值面提取算法，使該類算法僅需知道部分控制點便可推斷出其他的控制點直至整個地質(zhì)曲面，即數(shù)據(jù)驅(qū)動+知識驅(qū)動算法。隨著帶符號的距離場、Marching Cubes(Lorensen et al.,1987)以及MT(Zhou et al.,1995)等算法的相繼出現(xiàn), 隱式建模算法越來越受到重視，該類方法研究的核心是遵循地質(zhì)規(guī)律的帶符號距離場的空間插值問題。Savchenko等(1995)首次提出隱式曲面重構(gòu)的數(shù)學(xué)表達式，即存在體函數(shù)φ(x,y,z)，當(dāng)φ(x,y,z)=0時表示體的表面。此類方法主要應(yīng)用于醫(yī)學(xué)與地學(xué)領(lǐng)域。

隱式建模廣泛采用的插值算法主要有地質(zhì)統(tǒng)計學(xué)與徑向基函數(shù)(Radial Basis Function, RBF)法等。Lajaunie等(1997)提出應(yīng)用多變量地質(zhì)統(tǒng)計學(xué)算法解決地層曲面的模擬問題；在此基礎(chǔ)上，Calcagno 等(2008)提出了基于地質(zhì)知識的隱式地質(zhì)建模方法模型并討論了單一斷裂等構(gòu)造曲面的提取問題。同時，服務(wù)于隱式建模的軟件GeoModeller應(yīng)運而生。

2 區(qū)域大比例尺隱式建模技術(shù)應(yīng)用

筆者所用的隱式建模方法是基于對偶協(xié)同克里金插值的位勢場方法，該法利用地質(zhì)界面點等式約束和方向數(shù)據(jù)的梯度方向，采取對偶協(xié)同克里金插值法得到標(biāo)量場，其中的一系列等值面視為地質(zhì)界面。為隱式建模技術(shù)的算法基礎(chǔ)，可以運用的地質(zhì)約束歸結(jié)為地質(zhì)界面點的等式約束、方向數(shù)據(jù)的梯度方向約束、地質(zhì)界面點之間的不等式約束等幾類。

根據(jù)隱式建模技術(shù)的特點，三維建模主要使用的數(shù)據(jù)有礦區(qū)地表地質(zhì)信息、音頻大地電磁測深(CSAMT)剖面解譯數(shù)據(jù)、航磁解譯結(jié)果及剖面數(shù)據(jù)、鉆孔數(shù)據(jù)等。建模流程(圖1)由礦區(qū)原始數(shù)據(jù)采集、數(shù)據(jù)處理及數(shù)據(jù)輸入、創(chuàng)建地質(zhì)界面、生成曲面模型和生成實體模型5步組成(吳曉貴等，2020)。

圖1 GeoModeller軟件三維建模流程圖Fig. 1 GeoModeller 3D modeling flowchart

2.1 原始數(shù)據(jù)收集

研究區(qū)礦床為一中低溫?zé)嵋盒徒鸬V床，礦體主要產(chǎn)于向斜的南翼地層1(炭質(zhì)板巖)中，少量賦存于地層2中(圖2)，受兩側(cè)侵入巖體擠壓，在向斜南翼東段形成弧形構(gòu)造，地層1及部分地層2的炭質(zhì)板巖提供地球化學(xué)障，形成容礦空間。需要構(gòu)建的賦礦地層為地層1和地層2。

成礦作用方式為含金流體在Ⅲ號脆韌性剪切構(gòu)造帶(圖2)的控礦構(gòu)造中充填交代，關(guān)鍵控礦因素是構(gòu)造變形，根據(jù)音頻大地電磁測深(CSAMT)L1剖面、航磁解譯T圖及航磁解譯剖面S1—S5(圖3)，向斜南翼核部段可能存在第二個成礦富集段，應(yīng)作為深部找礦的主要目標(biāo)。在地層1和地層2的建?；A(chǔ)上，向斜核部段的Ⅲ號脆韌性剪切帶與地層1和地層2疊加區(qū)域應(yīng)作為主要建模范圍。

圖2 研究區(qū)地形地質(zhì)圖1-地層1；2-地層2；3-地層3；4-地層4；5-地層5；6-地層6；7-物探剖面；8-脆性斷裂帶；9-脆韌性斷裂帶Fig. 2 Topographic and geological map of the study area

2.2 地質(zhì)模型的建立

所用軟件為GeoModeller，支持廣泛的數(shù)據(jù)格式。建模所用數(shù)據(jù)有礦區(qū)1∶10 000地形地質(zhì)圖、1∶25 000航磁解譯圖、1幅CSAMT剖面圖、1套航磁原始數(shù)據(jù)、5幅航磁剖面圖、1套鉆孔測斜數(shù)據(jù)。

2.2.1 工程創(chuàng)立 首先定義建模所需的工程范圍，然后給出建模準(zhǔn)確度及導(dǎo)入該地區(qū)地形數(shù)據(jù)。主要圍繞礦區(qū)控礦構(gòu)造和賦礦地層進行建模(圖4)，面積為44.85 km2，深度為2.8 km。

圖4 工程定義示意圖Fig. 4 Engineering definition interface

2.2.2 地層構(gòu)建 在建模前需先構(gòu)建地質(zhì)模塊，如地層、巖基、構(gòu)造等，并依據(jù)從老到新的準(zhǔn)則對整個礦區(qū)構(gòu)建一個序列。建模的地層從老到新為地層6—地層1(圖5)，各地層之間的位置關(guān)系按地層的屬性onlap(沉積)設(shè)置。

圖5 地層構(gòu)建示意圖Fig. 5 Formation construction interface

2.2.3 地表數(shù)據(jù)導(dǎo)入 地表數(shù)據(jù)主要包括地表控制線和產(chǎn)狀信息。根據(jù)礦區(qū)地形地質(zhì)圖添加/設(shè)置各類型地質(zhì)界線(圖6)。在添加/設(shè)置各類型地質(zhì)界線時對每個地質(zhì)體依次描繪，逐步完善地質(zhì)體界線及其產(chǎn)狀信息。在描繪過程中要穿插進行計算調(diào)試工作，查看地質(zhì)界線的準(zhǔn)確性，進行校準(zhǔn)。在地表控制線和產(chǎn)狀信息全部輸入完成后，可先計算一個初步模型，由于目前控制數(shù)據(jù)過少，所以模型計算的結(jié)果往往偏差較大，需進一步添加控制數(shù)據(jù)。

圖6 地表數(shù)據(jù)導(dǎo)入Fig. 6 Surface data import

圖3 地球物理綜合解譯圖(a)CSAMT-L1剖面；(b)航磁T圖；(c)S1—S5航磁剖面Fig. 3 Comprehensive geophysical interpretation(a) CSAMT-L1 profile; (b) Aeromagnetic T-diagram; (c) S1-S5 aeromagnetic profiles

圖7 剖面數(shù)據(jù)添加(a)剖面添加；(b)剖面1添加；(c)L1剖面添加Fig. 7 Profile data addition(a) Section addition; (b) Section 1 addition;(c) L1 section addition

2.2.4 剖面數(shù)據(jù)增加 為提高模型準(zhǔn)確度，還需增加剖面數(shù)據(jù)控制以降低模型的不確定性。首先將所有需要的剖面數(shù)據(jù)的控制點坐標(biāo)投射到地表上(圖7a)，并對這些剖面依次進行命名，使其得以在二維窗口中單獨顯示；然后結(jié)合模型形態(tài)特征，在每條剖面的二維界面中添加剖面控制線和產(chǎn)狀信息等。以剖面1(圖7b)為例，由于該礦區(qū)深部還存在隱伏花崗閃長巖體，在剖面二維視圖中添加音頻大地電磁測深(CSAMT)L1剖面(圖7c)，以推斷隱伏巖體的控制線和產(chǎn)狀信息等，并根據(jù)CSAMT反演的賦礦地層及控礦構(gòu)造情況輔助修改模型，提高模型準(zhǔn)確度。

2.2.5 三維模型計算與顯示 將上述設(shè)置好的數(shù)據(jù)導(dǎo)入三維模型進行計算，生成地質(zhì)單元的三維地質(zhì)模型。按照礦區(qū)實際情況，還可對剖面上推斷的數(shù)據(jù)進行調(diào)整，以逼近真實地質(zhì)場景。同樣，在此部分還可設(shè)置計算網(wǎng)格的大小，也能手動輸入需要計算的區(qū)域(圖8a)，對于同樣的模型，計算精度與消耗的時間成反比。

最終建立的賦礦地層——地層1、地層2及隱伏花崗閃長巖基底的三維地質(zhì)模型見圖8b。在X、Y、Z方向上均可進行不同距離的切割，也可在任意方向進行旋轉(zhuǎn)，可進行單個層位和多個層位的顯示。分別沿X、Y、Z方向切割相應(yīng)的距離，可見模型內(nèi)部任意位置的組合形態(tài)。模型可以導(dǎo)出其他形式的網(wǎng)格文件，為模型后續(xù)改進提供可能。

2.2.6 脆韌性剪切帶建模 將圖3c的S1—S5航磁剖面導(dǎo)入地表，并在二維顯示窗口中結(jié)合模型形態(tài)特征，在每條剖面的二維界面中添加Ⅱ號脆性斷裂帶和Ⅲ號脆韌性剪切構(gòu)造帶的剖面控制線和產(chǎn)狀信息等，將剖面信息放入三維模型中計算，建立Ⅱ號脆性斷裂帶和Ⅲ號脆韌性剪切構(gòu)造帶的三維模型結(jié)構(gòu)。圖9所示模型即為褶皺構(gòu)造南北兩翼地層賦礦地層(地層1、地層2)和控礦構(gòu)造(Ⅲ號脆韌性剪切構(gòu)造帶)復(fù)合顯示圖。

圖8 三維模型計算與顯示Fig. 8 3D model calculation and display(a) Calculation adjustment interface; (b) 3D model of ore-controlling strata

圖9 賦礦地層-控礦構(gòu)造三維模型Fig. 9 3D model of ore-hosting strata and ore-controlling structure

3 結(jié) 論

針對一些礦區(qū)在區(qū)域大比例尺范圍鉆孔數(shù)據(jù)過少的問題，提出了一種三維地質(zhì)建模方法，得出下列結(jié)論。

(1) 以礦區(qū)數(shù)據(jù)為例，基于地表地質(zhì)測量數(shù)據(jù)，結(jié)合高精度地球物理數(shù)據(jù)，運用GeoModeller軟件，通過三維反演、三維可視化等技術(shù)，可直觀刻畫出該區(qū)域一定深度范圍內(nèi)的構(gòu)造、巖體、地層等的空間展布以及成因和演化關(guān)系。

(2) 為在區(qū)域上對控礦構(gòu)造(脆韌性剪切帶)和賦礦地層(地層1、地層2)的結(jié)合提供了一個新的成礦預(yù)測思路。