李 楠 宋相龍 楚文楷

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

0 引 言

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

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

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

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

1.1 三維地質建模的尺度

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

表1 三維地質建模研究的3個尺度

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

1.2 三維地質建模的算法

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

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

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

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

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

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

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

2 區(qū)域大比例尺隱式建模技術應用

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

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

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

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

研究區(qū)礦床為一中低溫熱液型金礦床，礦體主要產(chǎn)于向斜的南翼地層1(炭質板巖)中，少量賦存于地層2中(圖2)，受兩側侵入巖體擠壓，在向斜南翼東段形成弧形構造，地層1及部分地層2的炭質板巖提供地球化學障，形成容礦空間。需要構建的賦礦地層為地層1和地層2。

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

圖2 研究區(qū)地形地質圖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 地質模型的建立

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

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

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

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

圖5 地層構建示意圖Fig. 5 Formation construction interface

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

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

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

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

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

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

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

3 結 論

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

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

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