王麗芳，吳湘濱，張寶一,2，李小麗，楊 莉

1.中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室/地球科學(xué)與信息物理學(xué)院,長(zhǎng)沙 410083 2.中國(guó)地質(zhì)大學(xué)地質(zhì)過程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430074

?

一種可存儲(chǔ)路徑的三維非均質(zhì)空間最短距離場(chǎng)生成算法

王麗芳1，吳湘濱1，張寶一1,2，李小麗1，楊 莉1

1.中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室/地球科學(xué)與信息物理學(xué)院,長(zhǎng)沙 410083 2.中國(guó)地質(zhì)大學(xué)地質(zhì)過程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430074

針對(duì)成礦定量預(yù)測(cè)強(qiáng)調(diào)三維空間信息綜合研究的特殊需求和地質(zhì)要素定量表達(dá)難、深層次信息挖掘難等關(guān)鍵問題，筆者提出了一種基于三維非均質(zhì)地質(zhì)空間的最短距離場(chǎng)算法，并且可回溯出目標(biāo)體元到源體元的最短路徑。應(yīng)用該方法實(shí)現(xiàn)了對(duì)某銅礦區(qū)三維地質(zhì)實(shí)體模型深層次和組合信息的挖掘，在反映斷層對(duì)地質(zhì)空間分劃作用的基礎(chǔ)上，建立了侵入巖體的最短距離場(chǎng)，可用于對(duì)侵入巖體熱力場(chǎng)和圍巖蝕變等問題的分析。

距離場(chǎng)；三維；非均質(zhì)；路徑；成礦預(yù)測(cè)

0 引言

隨著現(xiàn)代成礦理論的發(fā)展、強(qiáng)穿透物化探技術(shù)的革新和地質(zhì)勘查程度的深入，二維空間地質(zhì)信息表達(dá)方式的局限性愈來愈明顯，需要從真三維空間角度來分析和解決礦產(chǎn)資源定量預(yù)測(cè)問題，發(fā)展綜合地質(zhì)數(shù)據(jù)的深層次三維空間分析手段，從而提升對(duì)礦化空間分布規(guī)律性的認(rèn)識(shí)來指導(dǎo)找礦，解決隱伏找礦工作中直接指示信息少、干擾大、三維展布復(fù)雜、數(shù)據(jù)挖掘難等問題。近年來，三維地質(zhì)建模已成為三維地理信息系統(tǒng)(3D GIS)、地球科學(xué)、三維計(jì)算機(jī)圖形學(xué)等各相關(guān)領(lǐng)域的研究熱點(diǎn)和前緣，一大批三維地質(zhì)建模軟件在地質(zhì)、石油、采礦等行業(yè)得到推廣使用，如GoCAD、Lynx、Surpac、Landmark、Mincom、EarthVision、Vulcan、Micromine、Datamine等。所謂三維地質(zhì)建模，即采用適當(dāng)?shù)臄?shù)據(jù)結(jié)構(gòu)在計(jì)算機(jī)中建立起能反映地質(zhì)構(gòu)造形態(tài)和各要素之間關(guān)系，以及地質(zhì)體物理、化學(xué)屬性空間分布等地質(zhì)特征的數(shù)學(xué)模型。Houlding[1]最早闡述了地下三維可視化技術(shù)的核心理論和基本方法，Lynx軟件系統(tǒng)就是基于該理論和方法體系實(shí)現(xiàn)的；Mallet[2]也提出了一種新的地質(zhì)目標(biāo)計(jì)算機(jī)輔助設(shè)計(jì)方法以適應(yīng)地質(zhì)、地球物理和油藏工程的需要，其理論和方法體系已經(jīng)成為GoCAD 軟件系統(tǒng)的主要技術(shù)支撐；武強(qiáng)等[3]提出了一種“逐步優(yōu)化”(stepwise refinement)的多源數(shù)據(jù)整合三維地質(zhì)建模方法，來漸進(jìn)、有效地提高三維模型的精度；吳立新等[4]將三維地質(zhì)建模的應(yīng)用分為三維可視化、三維查詢、三維分析、工程應(yīng)用4個(gè)目標(biāo)層次；潘懋等[5]概括了三維地質(zhì)建?！岸嘣磾?shù)據(jù)-多方法集成-多層次干預(yù)”的“三多”建模思路。三維地質(zhì)建模是將不同專業(yè)、期次、來源、類型的數(shù)據(jù)在統(tǒng)一環(huán)境下進(jìn)行集成，建設(shè)多源地學(xué)數(shù)據(jù)倉(cāng)庫(kù)的有效途徑，為后續(xù)三維可視化的數(shù)值分析提供了數(shù)據(jù)和技術(shù)支撐[6-7]，可以在工程勘察[8]、礦產(chǎn)勘查[9-10]和地下水調(diào)查[11]等領(lǐng)域開展更深入的研究。

建立起囊括地層、構(gòu)造、物化探等要素的區(qū)域三維地質(zhì)空間數(shù)據(jù)模型是開展三維地學(xué)模擬的前提條件[11]。三維地質(zhì)空間數(shù)據(jù)模型可以分為基于實(shí)體的模型和基于場(chǎng)的模型兩類[12-13]：前者將三維地質(zhì)空間表示為離散的點(diǎn)、線、面、體四類幾何實(shí)體的組合，側(cè)重于表達(dá)地質(zhì)構(gòu)造形態(tài)和各要素之間的關(guān)系，也稱為三維幾何模型；后者則將三維地學(xué)空間剖分為規(guī)則的或不規(guī)則的體元(voxel)構(gòu)成的連續(xù)場(chǎng)，側(cè)重于表達(dá)三維地質(zhì)空間內(nèi)物理、化學(xué)屬性的空間分布特征，也稱為三維屬性模型。三維地質(zhì)建模的位勢(shì)場(chǎng)方法(potential-field method)就采用了協(xié)同克里格方法將地質(zhì)接觸面的位置數(shù)據(jù)和構(gòu)造場(chǎng)的產(chǎn)狀數(shù)據(jù)插值成一個(gè)連續(xù)的三維位勢(shì)場(chǎng)標(biāo)量函數(shù)來表達(dá)地質(zhì)實(shí)體的幾何形狀，將地質(zhì)接觸位置集合作為等值面的參考點(diǎn)，而產(chǎn)狀數(shù)據(jù)為標(biāo)量函數(shù)的梯度，將斷層表示為位勢(shì)場(chǎng)中的不連續(xù)[14-16]。在進(jìn)行三維空間分析時(shí)，由于場(chǎng)數(shù)據(jù)便于空間疊加、空間相關(guān)和空間模擬等空間分析，因而通常使用三維場(chǎng)數(shù)據(jù)進(jìn)行三維空間分析和工程應(yīng)用。

三維柵格(3D-raster)模型是一種有代表性的規(guī)則數(shù)據(jù)場(chǎng)表示方法，它是一種無采樣約束的面向場(chǎng)物質(zhì)(如重力場(chǎng)、磁場(chǎng)、礦化品位)的連續(xù)空間標(biāo)準(zhǔn)分割方法，一般用八叉樹(octree)模型對(duì)其進(jìn)行數(shù)據(jù)壓縮及空間索引[17]，廣泛應(yīng)用于資源、水體、污染和環(huán)境等問題的建模。三維柵格地質(zhì)模型可以看作為一幅含有空間目標(biāo)(如地質(zhì)體、構(gòu)造、礦體等)的三維圖像，這幅三維圖像中不僅有空間目標(biāo)的分布范圍，還可利用GIS三維空間分析方法對(duì)其進(jìn)行深加工，從而衍生新的數(shù)據(jù)場(chǎng)來開展統(tǒng)計(jì)分析。

在地質(zhì)研究中，經(jīng)常需要對(duì)地表或地下的某些地質(zhì)體(如斷層、巖體等)進(jìn)行距離分析，因而距離場(chǎng)的提取及其形態(tài)特征的表達(dá)對(duì)于眾多地質(zhì)問題的解決具有特殊的意義。例如，在研究地下成礦元素運(yùn)移時(shí)，一個(gè)重要的問題就是分析成礦元素運(yùn)移到達(dá)的空間區(qū)域，這就需要考慮最短距離和運(yùn)移路徑。Payne和Toga[18]提出了一個(gè)從完整的封閉表面來計(jì)算距離場(chǎng)的算法，并檢驗(yàn)了距離場(chǎng)在表面移動(dòng)、內(nèi)插和縫合中的應(yīng)用；李芳玉[19]在Dijkstra算法與距離變換的基礎(chǔ)上，提出了一種基于三維柵格的最短距離算法，并將其應(yīng)用于非均質(zhì)的三維緩沖體的生成；李瑩等[20]以立方體品位模型為基礎(chǔ)，提出了一種按同級(jí)品位提取三維礦體輪廓的算法，從三維的角度分析礦體品位分布的不均勻性；陳建平等[21]以Micromine軟件為平臺(tái)建立了新疆阿勒泰地區(qū)可可托海3號(hào)偉晶巖脈9個(gè)共生結(jié)構(gòu)帶的三維實(shí)體模型和鉆孔三維品位模型，并利用鉆孔數(shù)據(jù)開展隱伏礦體的定位和儲(chǔ)量估算；安聰榮等[22]提出了一種可應(yīng)用于非均質(zhì)三維地質(zhì)空間中、允許存在多個(gè)源點(diǎn)的柵格模型最短距離等值面算法；毛先成等[23]以安徽銅陵鳳凰山礦田為例，采用8個(gè)模板進(jìn)行三維歐氏距離變換生成了柵格歐氏距離場(chǎng)，并利用開閉濾波或閉開濾波的數(shù)學(xué)形態(tài)學(xué)分析提取了地質(zhì)體表面的形態(tài)趨勢(shì)。

三維距離分析常用于構(gòu)建三維地質(zhì)變量，是礦體與控礦要素間關(guān)聯(lián)關(guān)系提取及成礦三維定量預(yù)測(cè)模型建立的基礎(chǔ)。毛先成等[24-25]分別以廣西大廠錫礦和安徽銅陵鳳凰山礦田為例，按照“地質(zhì)數(shù)據(jù)集成--成礦信息定量提取--立體定量預(yù)測(cè)”的流程，從場(chǎng)分析的觀點(diǎn)建立控礦地質(zhì)因素的場(chǎng)模型以及礦化分布場(chǎng)模型，制作統(tǒng)計(jì)數(shù)據(jù)散點(diǎn)圖得到控礦地質(zhì)因素與礦化指標(biāo)的關(guān)聯(lián)關(guān)系，建立了反映這種關(guān)聯(lián)關(guān)系的礦體立體定量預(yù)測(cè)模型；張寶一等[10]以遼寧紅透山銅礦為例，建立了成礦信息場(chǎng)、地層形態(tài)場(chǎng)、斷裂影響場(chǎng)、地球物理場(chǎng)等多源地學(xué)空間數(shù)據(jù)場(chǎng)，并提取了三維定量預(yù)測(cè)的地質(zhì)控礦指標(biāo)集和建立了預(yù)測(cè)模型。

為了滿足成礦定量預(yù)測(cè)強(qiáng)調(diào)三維空間信息綜合研究的特殊需求，解決地質(zhì)要素定量表達(dá)難、深層次信息挖掘難等關(guān)鍵問題，筆者提出了一種基于三維非均質(zhì)地質(zhì)空間的最短距離場(chǎng)算法，其可回溯出目標(biāo)體元到源體元的最短路徑；應(yīng)用該方法建立了某銅礦區(qū)侵入巖體的最短距離場(chǎng)，可用于對(duì)侵入巖體的熱力場(chǎng)和圍巖蝕變等問題的分析，場(chǎng)的非均質(zhì)性可反映出斷層對(duì)熱力場(chǎng)的分劃作用。

1 三維地質(zhì)空間場(chǎng)的建立

以面向成礦預(yù)測(cè)應(yīng)用為例，利用礦區(qū)構(gòu)造或礦體的總體趨勢(shì)來確定體元的三軸方向u、v、w(默認(rèn)分別與坐標(biāo)軸x、y、z平行)，每個(gè)體元i用其在3個(gè)軸向上的序號(hào)組合(ui,vi,wi)來唯一標(biāo)識(shí)，三維地質(zhì)空間場(chǎng)定義為按規(guī)則格網(wǎng)陣列記錄的三維數(shù)組。三維地質(zhì)空間場(chǎng)的結(jié)構(gòu)特征(如格網(wǎng)起始點(diǎn)、分辨率及方向等)直接影響著地質(zhì)空間數(shù)據(jù)場(chǎng)的構(gòu)建和空間相關(guān)場(chǎng)的計(jì)算精度。按照一定體元大小(格網(wǎng)分辨率)對(duì)所描述的地質(zhì)對(duì)象進(jìn)行重建時(shí)，格網(wǎng)分辨率直接確定了數(shù)據(jù)量和計(jì)算精度，一般通過計(jì)算礦化品位各方向的變差函數(shù)，根據(jù)變程值、礦體形態(tài)及誤差范圍綜合確定格網(wǎng)分辨率。

筆者使用SURPAC軟件建立了某銅礦區(qū)地質(zhì)體三維實(shí)體模型(圖1a)，該模型是僅以不規(guī)則三角網(wǎng)(TIN)面表示的空間閉合實(shí)體。在進(jìn)行地質(zhì)控礦指標(biāo)、礦化分布指標(biāo)提取時(shí)，按照單元尺寸為10 m×10 m×10 m的體元將建模范圍內(nèi)的三維地質(zhì)空間剖分為規(guī)則場(chǎng)，并結(jié)合礦區(qū)的數(shù)字高程模型(DEM)剔除地表以上的體元(見附錄A1)，如圖1b所示。三維地質(zhì)空間場(chǎng)涵蓋了從礦山收集到的地質(zhì)、物化探、勘探資料的三維空間范圍，統(tǒng)一的三維空間場(chǎng)結(jié)構(gòu)便于地、物、化、遙等各專業(yè)信息的綜合分析。

2 三維地質(zhì)空間場(chǎng)的賦值

三維地質(zhì)空間場(chǎng)的結(jié)構(gòu)確定以后，首要任務(wù)就是結(jié)合三維地質(zhì)實(shí)體模型，對(duì)三維場(chǎng)體元進(jìn)行多源地學(xué)信息屬性賦值，主要包括：將三維地質(zhì)實(shí)體模型轉(zhuǎn)換為體元模型，這一過程可看作是對(duì)地質(zhì)體的三維柵格化，圖2a和圖2b分別為將三維侵入巖體和斷層體元實(shí)體模型作為約束條件進(jìn)行柵格化得到的體元模型；將三維地質(zhì)實(shí)體模型轉(zhuǎn)換為三維位勢(shì)場(chǎng)的標(biāo)量場(chǎng)，整合及假整合地層界面表示為等值面，斷層及其他不整合面作為不連續(xù)面；利用勘探工程(槽、井、坑、鉆)采樣數(shù)據(jù)進(jìn)行重采樣和空間插值(如地質(zhì)統(tǒng)計(jì)學(xué)方法、距離反比加權(quán)法等)，對(duì)三維地質(zhì)空間場(chǎng)開展主要成礦元素的屬性估值，對(duì)三維物化、鉆井化探等數(shù)據(jù)進(jìn)行場(chǎng)屬性賦值，其他二維物、化探及遙感數(shù)據(jù)只進(jìn)行地表體元賦值。接下來需設(shè)計(jì)三維最短距離場(chǎng)算法，計(jì)算非均質(zhì)地質(zhì)空間中任意體元相對(duì)于多個(gè)“地質(zhì)源實(shí)體”的最短距離，將地質(zhì)源實(shí)體(如斷裂、侵入巖體等)的影響也表示為標(biāo)量場(chǎng)。

圖1 某銅礦區(qū)三維地質(zhì)實(shí)體模型(a)與三維地質(zhì)空間場(chǎng)(b)Fig.1 3D geological entity model (a) and 3D geological spatial field (b) of a copper deposit area

a.侵入巖體；b.斷層。圖2 某銅礦區(qū)三維地質(zhì)體元模型Fig. 2 3D geological voxel model of a copper deposit area

3 加權(quán)最短距離場(chǎng)的生成算法

3.1 最短距離場(chǎng)的定義

最短距離場(chǎng)可表示為每個(gè)體元k到含有m個(gè)“源實(shí)體”(可以是地質(zhì)界面或感興趣區(qū)域)的集合S中所有實(shí)體的距離中的最小值的集合：

(1)

式中：Fd為最短距離場(chǎng)的場(chǎng)值；d(,)為距離計(jì)算函數(shù)。

定義出發(fā)點(diǎn)所在的體元為源體元，目標(biāo)點(diǎn)所在的體元為目標(biāo)體元，它們處在同一格網(wǎng)中。距離量算就是要計(jì)算目標(biāo)體元與源體元的距離。設(shè)體元的大小為1單元，任意兩個(gè)體元之間的距離可以有多種定義方式，最常用的一種是明考夫斯基(Minkowsky)距離：

(2)

式中：(ui,vi,wi)和(uj,vj,wj)表示任意兩個(gè)體元在3個(gè)軸向上的序號(hào)；當(dāng)q=2時(shí)，d2(i,j)表示的就是歐氏距離；當(dāng)q=1時(shí)，d1(i,j)表示的是絕對(duì)值距離；當(dāng)q→∞時(shí)，d∞(i,j)表示的是切比雪夫距離。

在實(shí)際應(yīng)用中，為了精確表達(dá)地質(zhì)體及地質(zhì)現(xiàn)象、滿足大比例尺定量預(yù)測(cè)需求，格網(wǎng)的分辨率一般要非常精細(xì)，一個(gè)礦區(qū)的體元數(shù)量甚至可以達(dá)到上百萬個(gè)。如果對(duì)每個(gè)體元求取其中心點(diǎn)到一個(gè)或多個(gè)三維表面的距離，計(jì)算量很大，因此通常采用近似的場(chǎng)模型的距離量算來替代。本文設(shè)計(jì)了一種“種子生長(zhǎng)算法”，可同時(shí)實(shí)現(xiàn)最短距離場(chǎng)的生成及路徑搜索，并允許在三維地質(zhì)空間中存在多個(gè)“源實(shí)體”，以及體元可加載“費(fèi)用”權(quán)值(默認(rèn)等于1)來表達(dá)距離場(chǎng)的非均質(zhì)性。

3.2 鄰接體元的編碼

三維場(chǎng)模型的距離是以各體元的中心點(diǎn)為節(jié)點(diǎn)、以體元的鄰接關(guān)系作為連接進(jìn)行計(jì)算的。三維場(chǎng)中對(duì)每個(gè)體元有一個(gè)3×3×3的以該體元為中心的鄰域，用0--26對(duì)其鄰域內(nèi)的體元進(jìn)行編碼，該編碼也用來表示從源體元指向鄰接體元的方向，方向編碼(Dir)為0則指向自身。每個(gè)體元最多與26個(gè)體元相鄰，體元之間的鄰接形式有3種：面鄰接6個(gè)(編碼1--6)，與源體元共用一個(gè)面，距離源體元1單元；邊鄰接12個(gè)(編碼7--18)，與源體元共用一條邊，距離源體元2(1/q)單元；點(diǎn)鄰接8個(gè)(編碼19--26)，與源體元共用一個(gè)頂點(diǎn)，距離源體元3(1/q)單元，如圖3所示。

圖3 26個(gè)鄰接體元的編碼Fig.3 Coding of 26 adjacent voxels

3.3 加權(quán)最短距離場(chǎng)的生成

設(shè)源體元k和鄰接體元k′加載的“費(fèi)用”權(quán)值分別為Ck和Ck′，則源體元k到其目標(biāo)體元k′的加權(quán)距離為

(3)

應(yīng)用“種子生長(zhǎng)算法”同時(shí)生成最短距離場(chǎng)及最短路徑搜索的具體步驟(圖4)如下：

第一步，將離散的體元分成2個(gè)集合，特征體元集合S和目標(biāo)體元集合S′。集合S為“源實(shí)體”占據(jù)的體元，其距離場(chǎng)值為0(近似值)，其向量的Dir為0；集合S′為三維場(chǎng)中其他體元集合，其距離場(chǎng)值默認(rèn)為∞，其向量的Dir為無效值-1。

第二步，以集合S中的每個(gè)體元o為“種子”，向其26個(gè)相鄰體元“生長(zhǎng)”，按公式(3)計(jì)算相鄰體元o′的加權(quán)距離場(chǎng)值：

(4)

如果FdC(o′)小于o′原距離場(chǎng)值，則用FdC(o′)替換原距離場(chǎng)值(反之，沿用原距離場(chǎng)值)，并將向量o′o(從體元o′中心節(jié)點(diǎn)指向體元o中心節(jié)點(diǎn)的向量)的Dir存儲(chǔ)在體元o′中。

第三步，從集合S′中選取距離場(chǎng)值最小的體元k作為“種子”，向體元k的26個(gè)相鄰體元“生長(zhǎng)”，按公式(3)計(jì)算相鄰體元k′的加權(quán)距離場(chǎng)值：

(5)

如果FdC(k′)小于k′原距離場(chǎng)值，則替換之(反之，沿用原距離場(chǎng)值)，并將向量k′k的Dir存儲(chǔ)在體元k′中，然后將體元k從集合S′移動(dòng)到集合S中。

第四步，循環(huán)執(zhí)行第三步，直至集合S′為空時(shí)停止，成功生成最短距離場(chǎng)。

圖4 “種子生長(zhǎng)法”加權(quán)最短距離場(chǎng)生成算法流程圖Fig. 4 Flowchart of ‘seed-growing’ weighted shortest distance field generation

第五步，對(duì)任意體元t′執(zhí)行路徑搜索時(shí)，只需讀取其向量方向編碼，獲取指向的體元t；再讀取體元t的向量方向編碼，獲取其指向的體元。循環(huán)調(diào)用該過程，直至最終體元的Dir為0(源實(shí)體所占據(jù)的體元)，就得到了體元t′到達(dá)“源實(shí)體”的一條最短路徑。

a.全局圖；b.柵狀圖。圖5 某銅礦區(qū)侵入巖體的最短距離場(chǎng)Fig.5 Intrusive rocks’ shortest distance field of a copper deposit area

a.方向編碼；b.路徑回溯。圖6 根據(jù)方向編碼進(jìn)行路徑回溯Fig.6 Tracing back a path according to direction code

4 計(jì)算結(jié)果

為了分析研究區(qū)內(nèi)與成礦作用有關(guān)的侵入巖體的熱力場(chǎng)分布和圍巖蝕變情況，采用文中的“種子生長(zhǎng)算法”對(duì)三維地質(zhì)模型中的侵入巖體計(jì)算了最短距離場(chǎng)，將定性數(shù)據(jù)轉(zhuǎn)變?yōu)槎繑?shù)據(jù)，實(shí)現(xiàn)對(duì)空間數(shù)據(jù)的深加工，見附錄A2和附錄A3。最短距離場(chǎng)可克服三維GIS緩沖區(qū)分析只能識(shí)別出距離侵入巖體某一范圍內(nèi)的鄰域(且該鄰域是均質(zhì)的)、且無法反映出侵入巖體熱力場(chǎng)依距離漸變性的缺陷(圖5a)；設(shè)置斷層對(duì)應(yīng)體元的“費(fèi)用”權(quán)值為10(其他體元費(fèi)用權(quán)值默認(rèn)為1)，使斷層所在單元在距離場(chǎng)計(jì)算中起到了屏蔽作用，以模擬斷層對(duì)熱力場(chǎng)的空間分劃作用(圖5b)。

每個(gè)體元除存儲(chǔ)其到最近的源體元的距離值外，還存儲(chǔ)有用于追蹤距離變換的上一級(jí)體元的方向向量編碼，如圖6a所示；通過逐級(jí)追蹤直至源為止，可以得到一條體元到源體元的最短路徑(見附錄A4)，如圖6b所示。

5 結(jié)論

本文提出對(duì)與成礦作用有關(guān)的侵入巖體計(jì)算最短距離場(chǎng)，能夠表達(dá)斷層對(duì)地質(zhì)空間分劃作用的非均質(zhì)性；在生成距離場(chǎng)的同時(shí)通過記錄上級(jí)體元方向編碼，使每個(gè)單元均可回溯出一條到最近源體元的最短路徑。進(jìn)一步研究將包括提取地質(zhì)空間場(chǎng)內(nèi)巖體封閉邊界通過的體元，從而支持同時(shí)計(jì)算巖體內(nèi)部和外部的距離場(chǎng)值(內(nèi)部場(chǎng)值設(shè)置為負(fù)值，外部場(chǎng)值設(shè)置為正值，邊界處場(chǎng)值為零值)，為更精細(xì)地分析侵入巖體的熱力場(chǎng)分布提供支撐。利用場(chǎng)分析技術(shù)對(duì)三維地質(zhì)實(shí)體模型進(jìn)行深層次和組合信息挖掘，既能借助三維可視化仿真環(huán)境，將地、物、化、遙等多源地學(xué)空間數(shù)據(jù)集成在統(tǒng)一的三維空間數(shù)據(jù)場(chǎng)構(gòu)架下，又為成礦定量預(yù)測(cè)模型中三維地質(zhì)變量的構(gòu)造和選擇提供了支持。

紅透山礦業(yè)有限公司石長(zhǎng)巖副總經(jīng)理、遼寧有色地質(zhì)局101隊(duì)地勘院楊鐵軍副院長(zhǎng)在資料收集、坑道調(diào)查工作中給予了細(xì)致周密的安排和熱情協(xié)助，國(guó)家地理信息系統(tǒng)工程技術(shù)研究中心與中南大學(xué)共建“MAPGIS實(shí)驗(yàn)室”為本研究工作提供了MAPGIS軟件，在此一并表示感謝！

附錄：Microsoft SQL Server 2012源代碼

A1 建立初始三維空間數(shù)據(jù)場(chǎng)

--建立地質(zhì)空間數(shù)據(jù)場(chǎng)表格HTS101010F

CREATE TABLE HTS101010F(

ID INT NOT NULL UNIQUE,--唯一的ID

IX INT NOT NULL,--X方向序號(hào)

IY INT NOT NULL,--Y方向序號(hào)

IZ INT NOT NULL,--Z方向序號(hào)

X INT NOT NULL,--X坐標(biāo)，m

Y INT NOT NULL,--Y坐標(biāo)，m

Z INT NOT NULL,--Z坐標(biāo)，m

TO_SOURCE CHAR(1),--是否是源體元，1-源體元/0-其他體元

WEIGHTNESS SMALLINT,--權(quán)重值，默認(rèn)值為1，斷層體元為10

PATH_TO SMALLINT,--到上一級(jí)體元的方向編碼

DISTANCE DECIMAL(9,3),--最短距離場(chǎng)值，m

CUR_STATE CHAR(1)--是否計(jì)算距離場(chǎng)，1-已計(jì)算/0-未計(jì)算

);

--對(duì)唯一的ID建立索引

CREATE UNIQUE INDEX HTS101010FID ON HTS101010F(ID);

--對(duì)Z方向序號(hào)建立索引

CREATE INDEX HTS10101FIZ ON HTS101010F(IZ);

--對(duì)Y方向序號(hào)建立索引

CREATE INDEX HTS10101FIY ON HTS101010F(IY);

--對(duì)X方向序號(hào)建立索引

CREATE INDEX HTS10101FIX ON HTS101010F(IX);

--向地質(zhì)空間數(shù)據(jù)場(chǎng)表格中插入記錄的存儲(chǔ)過程

CREATE PROCEDURE [dbo].[InsertAllReords]

@minx int, --X坐標(biāo)最小值

@miny int, --Y坐標(biāo)最小值

@minz int, --Z坐標(biāo)最小值

@maxx int, --X坐標(biāo)最大值

@maxy int, --Y坐標(biāo)最大值

@maxz int, --Z坐標(biāo)最大值

@minix int,--X方向序號(hào)最小值

@miniy int,--Y方向序號(hào)最小值

@miniz int,--Z方向序號(hào)最小值

@maxix int,--X方向序號(hào)最大值

@maxiy int,--Y方向序號(hào)最大值

@maxiz int,--Z方向序號(hào)最大值

@nx int,--X方向單元數(shù)

@ny int,--Y方向單元數(shù)

@step int--單元尺寸，步長(zhǎng)

AS

BEGIN

Declare @ix int = @minix, @iy int = @miniy, @iz int = @miniz

Declare @id int, @x int, @y int, @z int

Declare @dem float

-- 刪除三維地質(zhì)空間數(shù)據(jù)場(chǎng)表格中的記錄

DELETE FROM HTS101010F

set @iy = @miniy

While (@iy <= @maxiy)

Begin --開始Y方向循環(huán)

set @y = @miny + (@iy-@miniy)*@step

set @ix = @minix

While (@ix <= @maxix)

Begin --開始X方向循環(huán)

set @x = @minx + (@ix-@minix)*@step

--根據(jù)X、Y坐標(biāo)從DEM數(shù)據(jù)表DZT5KL01E中獲取地表高程值

select @dem = (Select dem from DZT5KL01E WHERE X=@x and Y=@y)

set @iz= @miniz

While (@iz <= @maxiz)

Begin --開始Z方向循環(huán)

set @id = (@iz-1)*@nx*@ny + (@iy-1)*@nx + @ix

set @z = @minz + (@iz-@miniz)*@step

if(@z > @dem) break --如果高于地表高程，不插入表記錄

--向三維地質(zhì)空間數(shù)據(jù)場(chǎng)表格HTS101010F中插入記錄

Insert into HTS101010F

values(@id, @ix, @iy, @iz, @x, @y, @z, 0,1,null,null,0)

set @iz = @iz+1

END --結(jié)束Z方向循環(huán)

set @ix = @ix+1

END --結(jié)束X方向循環(huán)

set @iy = @iy+1

END --結(jié)束Y方向循環(huán)

END

A2 設(shè)置源體元和距離場(chǎng)權(quán)值

-- 將侵入巖體通過的體元設(shè)置為源體元

update HTS101010f

set to_source = '1', PATH_TO = 0, Distance = 0, CUR_STATE = 0

where hts101010f.id in

(select hts101010.id from HTS101010 where HTS101010.xp > 0);

-- 設(shè)置其他體元的的距離場(chǎng)值為初始值99999

update HTS101010f

set DISTANCE = 99999

where DISTANCE IS NULL

-- 將斷層通過的體元的權(quán)值設(shè)置為10

update HTS101010f

set weightness = 10

where hts101010f.id in

(select hts101010.id from HTS101010 where HTS101010.f > 0);

A3 建立非均質(zhì)的最短距離場(chǎng)

-- 生成最短距離場(chǎng)的存儲(chǔ)過程

CREATE PROCEDURE[dbo].[DistanceFields]

AS

BEGIN

DECLARE @MINDIST REAL = 0

WHILE @MINDIST!=-999 --返回距離值無效時(shí)結(jié)束循環(huán)

BEGIN

SET @MINDIST = -999

--對(duì)當(dāng)前最小距離的體元運(yùn)行生長(zhǎng)算法

EXEC MinDistanceGrowth @MINDIST OUTPUT --返回當(dāng)前最小距離

END

END

-- 當(dāng)前最小距離值體元生長(zhǎng)的存儲(chǔ)過程

CREATE PROCEDURE [dbo].[MinDistanceGrowth]

@mindist real = -9999 out --返回當(dāng)前最小距離

AS

BEGIN

-- 獲取距離場(chǎng)計(jì)算狀態(tài)未知的最小距離值的體元

set @mindist = (select min(distance) from HTS101010F where distance is not null and CUR_STATE=0) --獲得當(dāng)前的最小距離值

-- 更新當(dāng)前最小距離值的體元的距離場(chǎng)計(jì)算狀態(tài)為完成

UPDATE HTS101010F SET CUR_STATE =1 WHERE DISTANCE = @mindist

-- 定義光標(biāo)獲取當(dāng)前最小距離值的所有體元的ID和權(quán)值

DECLARE pIDs CURSOR FOR select ID,WEIGHTNESS from HTS101010F where distance= @mindist order by ID

OPEN pIDs --打開光標(biāo)

DECLARE @ID INT= 0, @WEIGHT SMALLINT = 0

fetch next from pIDs into @ID, @WEIGHT--讀取第一個(gè)體元的ID和權(quán)值

while @@fetch_status = 0 --如果當(dāng)前光標(biāo)無效，結(jié)束循環(huán)

BEGIN

-- 更新該體元鄰近26個(gè)體元的距離場(chǎng)值

EXEC UpdateNeighborDistance @ID,@mindist,@WEIGHT

fetch next from pIDs into @ID, @WEIGHT --讀取下一個(gè)體元的ID和權(quán)值

END

CLOSE pIDs --關(guān)閉光標(biāo)

Deallocate pIDs --銷毀光標(biāo)

END

-- 根據(jù)當(dāng)前體元距離值更新鄰近體元的值

CREATE PROCEDURE [dbo].[UpdateNeighborDistance]

@ID INT = 0,--當(dāng)前體元ID

@DIST real = -999,--當(dāng)前體元距離場(chǎng)值

@WEIGHT SMALLINT = -1--當(dāng)前體元費(fèi)用權(quán)值

AS

BEGIN

DECLARE @IX INT=0, @IY INT=0, @IZ INT=0

-- 根據(jù)ID獲取X、Y、Z方向上的序號(hào)

EXEC GetXYZfromID @ID,@IX OUTPUT,@IY OUTPUT,@IZ OUTPUT

-- 定義光標(biāo)選取當(dāng)前體元鄰近的26個(gè)體元

DECLARE pNeighbors CURSOR FOR SELECT

ID,IX,IY,IZ,CUR_STATE,DISTANCE,WEIGHTNESS,PATH_TO

FROM HTS101010F

WHERE IX >=@IX-1 AND IX <=@IX+1

AND IY>=@IY-1 AND IY <=@IY+1

AND IZ>=@IZ-1 AND IZ <=@IZ+1

AND CUR_STATE = 0

AND DISTANCE > @Dist -- 大于當(dāng)前體元距離

ORDER BY ID

OPEN pNeighbors --打開光標(biāo)

DECLARE @NID INT = 0, @NIX INT=0, @NIY INT=0, @NIZ INT=0,

@NSTATE SMALLINT=0, @NDIST REAL = 0,@NWEIGHT SMALLINT =0, @NPATH SMALLINT =0

fetch next from pNeighbors into --讀取第一個(gè)鄰近體元

@NID,@NIX,@NIY,@NIZ,@NSTATE,@NDIST,@NWEIGHT,@NPATH

while @@fetch_status = 0 --如果當(dāng)前光標(biāo)無效結(jié)束循環(huán)

BEGIN

DECLARE @TEMPDIST REAL = 99999

IF @NSTATE = 0 --如果該鄰近體元未完成距離場(chǎng)計(jì)算

BEGIN

-- 根據(jù)鄰近體元的位置設(shè)置其距離場(chǎng)值@TEMPDIST(略)

-- 根據(jù)鄰近體元的位置設(shè)置其方向編碼@NPATH(略)

IF @TEMPDIST < @NDIST --小于原距離值

BEGIN

-- 更新鄰近體元的距離場(chǎng)值和方向編碼

UPDATE HTS101010F SET

PATH_TO = @NPATH,

DISTANCE = @TEMPDIST

WHERE ID = @NID

END

END

--讀取下一個(gè)鄰近體元

fetch next from pNeighbors into

@NID, @NIX,@NIY,@NIZ,@NSTATE,@NDIST,@NWEIGHT,@NPATH

END

CLOSE pNeighbors --關(guān)閉光標(biāo)

Deallocate pNeighbors --銷毀光標(biāo)

END

A4 到源體元的最短路徑回溯

-- 獲得當(dāng)前體元到源體元的最短路徑

CREATE PROCEDURE GetPathByID

@id int

AS

BEGIN

DECLARE @code int = 0

--獲取當(dāng)前體元的方向編碼

SET @code = (select PATH_TO from HTS101010F where ID = @id)

--建立臨時(shí)表用于存儲(chǔ)路徑所經(jīng)體元的ID和方向編碼

DECLARE @PathTable TABLE(ID int,SrcCode int)

--將當(dāng)前體元添加到臨時(shí)表

INSERT INTO @PathTable VALUES (@id,@code)

WHILE @CODE !=0 --如果當(dāng)前體元是源體元，結(jié)束循環(huán)

BEGIN

declare @nid int = 0

-- 獲得上一級(jí)體元的ID

exec GetNeighborByCode @id,@code,@nid output

-- 獲取上一級(jí)體元的的方向代碼

SET @code = (select PATH_TO from HTS101010F where ID = @nid)

-- 將上一級(jí)體元添加到臨時(shí)表

INSERT INTO @PathTable VALUES (@nid,@code)

set @id = @nid -- 設(shè)置上一級(jí)體元為當(dāng)前體元

END

-- 將路徑經(jīng)過的體元從三維數(shù)據(jù)場(chǎng)表格中選擇出來

select * from HTS101010F where ID in ( select ID from @PathTable )

ORDER BY DISTANCE DESC

END

[1] Houlding S W. 3D Geoscience Modeling: Computer Techniques for Geological Characterization[M]. Berlin: Springer-Verlag, 1994.

[2] Mallet J L. Geomodeling[M]. New York: Oxford University Press, 2002.

[3] Wu Q, Xu H, Zou X K. An Effective Method for 3D Geological Modeling with Multi-Source Data Integration[J]. Computers and Geosciences, 2005, 31(1): 35-43.

[4] 吳立新, 車德福,郭甲騰. 面向地上下無縫集成建模的新一代三維地理信息系統(tǒng)[J]. 測(cè)繪工程, 2006, 15(5): 1-6. Wu Lixin, Che Defu, Guo Jiateng. The New 3DGIS for Seamless Integration of Terrain Overground and Underground Entities[J]. Engineering of Surveying and Mapping, 2006, 15(5): 1-6.

[5] 潘懋, 方裕,屈紅剛. 三維地質(zhì)建模若干基本問題探討[J]. 地理與地理信息科學(xué), 2007, 23(3): 1-5. Pan Mao, Fang Yu, Qu Honggang. Discussion on Several Foundational Issues in Three-Dimensional Geological Modeling[J]. Geography and Geo-Information Science, 2007, 23(3): 1-5.

[6] Zhang B Y, Wu X B, Wang L F, et al. The Preliminary Research of Feature-Based 3D Geological Modeling[C]//2nd Conference on Environmental Science and Information Application Technology (ESIAT’2010). Wuhan: IEEE Computer Society, 2010: 321-325.

[7] Zu X F, Hou W S, Zhang B Y, et al. Overview of Three-Dimensional Geological Modeling Technology[C]//Lee G. International Conference on Future Computer Supported Education (FCSE). Seoul: IERI Procedia, 2012: 921-927.

[8] Zu X F, Zhang B Y, Luo J, et al. 3D Modeling of Tunnel Engineering Based on Geological Body[C]//Lee G. International Conference on Future Computer Supported Education (FCSE). Seoul: IERI Procedia, 2012: 928-934.

[9] 張寶一, 尚建嘎, 吳鴻敏, 等. 三維地質(zhì)建模及可視化技術(shù)在固體礦產(chǎn)儲(chǔ)量估算中的應(yīng)用[J]. 地質(zhì)與勘探, 2007, 43(2): 76-81. Zhang Baoyi, Shang Jianga, Wu Hongmin, et al. Application of 3D Geological Modeling and Visualization in Solid Mineral Resource Estimation[J]. Geology and Prospecting, 2007, 43(2): 76-81.

[10] 張寶一, 吳湘濱, 王麗芳, 等. 紅透山銅礦外圍隱伏礦體三維定量預(yù)測(cè)[J]. 中國(guó)有色金屬學(xué)報(bào), 2012, 22(3): 863-871. Zhang Baoyi, Wu Xiangbin, Wang Lifang, et al. 3D Quantitative Prediction of Concealed Ore-Body in Surrounding Areas of Hongtoushan Copper Deposit[J]. The Chinese Journal of Nonferrous Metals, 2012, 22(3): 863-871.

[11] 張寶一, 吳湘濱, 王麗芳, 等. 三維地質(zhì)建模及應(yīng)用實(shí)例[J]. 地質(zhì)找礦論叢, 2013, 28(3): 328-336. Zhang Baoyi, Wu Xiangbin, Wang Lifang, et al. Three-Dimensional Geological Modeling and the Application Cases[J]. Contributions to Geology and Mineral Resources Research, 2013, 28(3): 328-336.

[12] 程朋根, 劉少華, 王偉, 等. 三維地質(zhì)模型構(gòu)建方法的研究及應(yīng)用[J]. 吉林大學(xué)學(xué)報(bào):地球科學(xué)版, 2004, 34(2): 309-313. Cheng Penggen, Liu Shaohua, Wang Wei, et al. Study and Application of a New 3D Geological Model Construction Method[J]. Journal of Jilin University: Earth Science Edition, 2004, 34(2): 309-313.

[13] 程朋根, 王承瑞, 甘衛(wèi)軍, 等. 基于多層DEM與QTPV的混合數(shù)據(jù)模型及其在地質(zhì)建模中的應(yīng)用[J]. 吉林大學(xué)學(xué)報(bào):地球科學(xué)版, 2005, 35(6): 806-811. Cheng Penggen, Wang Chengrui, Gan Weijun, et al. A Hybrid 3D Data Model Based on Multi-DEMs and QTPVs and Its Application in Geological Modeling[J]. Journal of Jilin University:Earth Science Edition, 2005, 35(6): 806-811.

[14] Calcagno P, Chilès J P, Courrioux G, et al. Geological Modelling from Field Data and Geological Knowledge: Part I: Modelling Method Coupling 3D Potential-Field Interpolation and Geological Rules[J]. Physics of the Earth and Planetary Interiors, 2008, 171(1/2/3/4): 147-157.

[15] De Kemp E A. 3-D Visualization of Structural Field Data: Examples from the Archean Caopatina Formation, Abitibi Greenstone Belt, Quebec, Canada[J]. Computers and Geosciences, 2000, 26(5): 509-530.

[16] Jessell M. Three-Dimensional Geological Modelling of Potential-Field Data[J]. Computers and Geosciences, 2001, 27(4): 455-465.

[17] Shi W. Development of A Hybrid Model for 3D GIS[J]. Geo-Spatial Information Science, 2000, 3(2): 6-12.

[18] Payne B A,Toga A W. Distance Field Manipulation of Surface Models[J]. Computer Graphics and Applications, IEEE, 1992, 12(1): 65-71.

[19] 李芳玉. 基于三維柵格的GIS距離分析算法研究與應(yīng)用[J]. 計(jì)算機(jī)工程與應(yīng)用, 2008, 44(15): 246-248. Li Fangyu. Research and Application of Distance Analysis Based on 3D Raster GIS[J]. Computer Engineering and Applications, 2008, 44(15): 246-248.

[20] 李瑩, 肖克炎,陳建平. 基于立方體模型的三維礦體模擬與資源評(píng)估[J]. 地質(zhì)通報(bào), 2010, 29(10): 1547-1553. Li Ying, Xiao Keyan, Chen Jianping. Ore Body Simulation and Resource Assessment Based on Three-Dimensional Cube Model[J]. Geological Bulletin of China, 2010, 29(10): 1547-1553.

[21] 陳建平, 陳勇, 朱鵬飛, 等. 數(shù)字礦床模型及其應(yīng)用:以新疆阿勒泰地區(qū)可可托海3號(hào)偉晶巖脈稀有金屬隱伏礦預(yù)測(cè)為例[J]. 地質(zhì)通報(bào), 2011, 30(6): 630-641. Chen Jianping, Chen Yong, Zhu Pengfei, et al. Digital Ore Deposit Model and Its Application: A Case Study of the Prognosis of the Koktokay No.3 Pegmatite Dike Concealed Rare Metal Deposit in Altay Area of Xinjiang[J]. Geological Bulletin of China, 2011, 30(6): 630-641.

[22] 安聰榮, 劉展, 黃榮剛, 等. 基于三維柵格模型的最短距離等值面提取[J]. 計(jì)算機(jī)工程, 2011, 37(8): 7-9. An Congrong, Liu Zhan, Huang Ronggang, et al. Shortest Distance Isosurface Extraction Based on Three-Dimensional Raster Model[J]. Computer Engineering, 2011, 37(8): 7-9.

[23] 毛先成, 唐艷華,鄧浩. 地質(zhì)體的三維形態(tài)分析方法與應(yīng)用[J]. 中南大學(xué)學(xué)報(bào):自然科學(xué)版, 2012, 43(2): 588-595. Mao Xiancheng, Tang Yanhua, Deng Hao. Three-Dimensional Morphological Analysis Method for Geologic Bodies and Its Application[J]. Journal of Central South University:Science and Technology, 2012, 43(2): 588-595.

[24] 毛先成, 戴塔根, 吳湘濱, 等. 危機(jī)礦山深邊部隱伏礦體立體定量預(yù)測(cè)研究:以廣西大廠錫多金屬礦床為例[J]. 中國(guó)地質(zhì), 2009, 36(2): 424-435. Mao Xiancheng, Dai Tagen, Wu Xiangbin, et al. The Stereoscopic Quantitative Prediction of Concealed Ore Bodies in the Deep and Marginal Parts of Crisis Mines: A Case Study of the Dachang Tin Polymetallic Ore Deposit in Guangxi[J]. Geology in China, 2009, 36(2): 424-435.

[25] 毛先成, 鄒艷紅, 陳進(jìn), 等. 危機(jī)礦山深部、邊部隱伏礦體的三維可視化預(yù)測(cè):以安徽銅陵鳳凰山礦田為例[J]. 地質(zhì)通報(bào), 2010, 29(2/3): 401-413. Mao Xiancheng, Zou Yanhong, Chen Jin, et al. Three-Dimensional Visual Prediction of Concealed Ore Bodies in the Deep and Marginal Parts of Crisis Mines: A Case Study of the Fenghuangshan Ore Field in Tongling, Anhui, China[J]. Geological Bulletin of China, 2010, 29(2/3): 401-413.

A Path-Savable Shortest-Distance Field Generating Algorithm in Three-Dimensional Heterogenetic Space

Wang Lifang1, Wu Xiangbin1, Zhang Baoyi1,2, Li Xiaoli1, Yang Li1

1.MOEKeyLaboratoryofMetallogenicPredictionofNonferrousMetalsandGeologicalEnvironmentMonitoring/SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,Changsha410083,China2.StateKeyLaboratoryofGeologicalProcessesandMineralResources,ChinaUniversityofGeosciences,Wuhan430074,China

A comprehensive study on three-dimensional (3D) spatial information is highly demanded for a quantitive prediction of mineralization. The difficulties are often confronted in geological feature quantitive expression and deep-seated information extraction. A shortest-distance field generating algorithm in 3D heterogenetic space is developed, which can track back to the shortest path from an object voxel to the source voxel. The algorithm has been used to a copper mine district for its deeper and comprehensive information based on the 3D geological entity model. A shortest-distance field to the intrusive rocks is calculated to analyze the geological problems, such as thermodynamic field and the alteration of surrounding rock in consideration of the roles played by the faults in dividing 3D geological space.

distance field; three-dimensional; heterogenetic; path；metallogenic prediction

10.13278/j.cnki.jjuese.201504305.

2014-10-29

國(guó)家自然科學(xué)基金項(xiàng)目(41302259)；“十一五”國(guó)家科技支撐計(jì)劃項(xiàng)目(2008BAB34B02)

王麗芳(1979--)，女，博士研究生，工程師，主要從事三維地質(zhì)建模研究，E-mail:csuwlf@139.com

張寶一(1979--)，男，在站博士后，副教授，主要從事地理信息科學(xué)研究，E-mail:zhangbaoyi@csu.edu.cn。

10.13278/j.cnki.jjuese.201504305

P628.3

A

王麗芳，吳湘濱，張寶一，等. 一種可存儲(chǔ)路徑的三維非均質(zhì)空間最短距離場(chǎng)生成算法.吉林大學(xué)學(xué)報(bào):地球科學(xué)版,2015,45(4):1257-1268.

Wang Lifang, Wu Xiangbin, Zhang Baoyi, et al. A Path-Savable Shortest-Distance Field Generating Algorithm in Three-Dimensional Heterogenetic Space.Journal of Jilin University:Earth Science Edition,2015,45(4):1257-1268.doi:10.13278/j.cnki.jjuese.201504305.