李 楠，雷玲玲，肖克炎，

(1．中國地質(zhì)科學(xué)院礦產(chǎn)資源研究所，北京 100037;2．長(zhǎng)江大學(xué)地球科學(xué)學(xué)院，湖北 荊州434023)

一種基于反距離加權(quán)方法的層狀三維地質(zhì)界面擬合算法

李 楠1，雷玲玲2，肖克炎1，2

(1．中國地質(zhì)科學(xué)院礦產(chǎn)資源研究所，北京 100037;2．長(zhǎng)江大學(xué)地球科學(xué)學(xué)院，湖北 荊州434023)

三維空間中地質(zhì)曲面模擬是礦產(chǎn)勘查、大比例尺成礦預(yù)測(cè)等地學(xué)領(lǐng)域的重要研究?jī)?nèi)容之一。與地球化學(xué)、地球物理數(shù)據(jù)的曲面擬合相比，地質(zhì)界面擬合更加復(fù)雜且具有一般性。從功能需求角度出發(fā)，地質(zhì)界面擬合需要處理層狀地層、不整合地層以及存在構(gòu)造的地層等多種情況;從算法實(shí)現(xiàn)角度講，該擬合過程不同于物化探數(shù)據(jù)通過某一屬性值進(jìn)行曲面擬合的情況，因?yàn)樵诘刭|(zhì)界面模擬過程中，不存在這樣的特征值進(jìn)行曲面擬合，因此不能直接使用物化探數(shù)據(jù)等值面提取的思想模擬地質(zhì)接觸界面。針對(duì)上述問題，提出了一種基于反距離加權(quán)法移動(dòng)立方格算法的層狀地層曲面擬合算法。該方法模型適合絕大多數(shù)層狀地學(xué)三維曲面模擬。

地質(zhì)界面;空間插值;移動(dòng)立方格算法;等值面;隱式地質(zhì)建模

0 引言

在地質(zhì)勘探過程中，勘探人員通過鉆孔、槽探、坑道、填圖和編錄等技術(shù)手段可以獲取地質(zhì)體和地質(zhì)界面的三維控制點(diǎn)信息。在計(jì)算機(jī)領(lǐng)域，這些控制點(diǎn)被稱為空間散亂點(diǎn)或離散點(diǎn)數(shù)據(jù)。由離散點(diǎn)集出發(fā)，重構(gòu)該點(diǎn)集所屬的曲面，稱為曲面重建。筆者所討論的地質(zhì)界面重建屬于該類問題。曲面重建，即構(gòu)造出對(duì)象幾何模型，是當(dāng)今三維數(shù)據(jù)場(chǎng)建模的核心問題之一。曲面重建問題抽象為：給定三維空間曲面S上的一組采樣點(diǎn)P，通過數(shù)學(xué)模型找到一個(gè)曲面 S'，使得 S'合理地逼近 S(Hoppe et al，1992)。其主要包括兩大類算法：一類是參數(shù)曲面模型法，另一類是曲面提取法。

參數(shù)曲面模型法主要包括逼近樣條曲面和插值樣條曲面兩種類型。其中，插值樣條通過離散點(diǎn)，而逼近樣條不過離散點(diǎn)。該方法的優(yōu)點(diǎn)是曲面的獲得不是由用戶給定某一特征值得到的，而是由某一數(shù)學(xué)模型擬合得到的。但無論上述哪種方法，生成的地質(zhì)界面均為單值面，不能很好地滿足地質(zhì)實(shí)際需要。文獻(xiàn)對(duì)上述方法做了詳細(xì)分析(毛先成，2006)。

曲面提取方法是地學(xué)曲面重構(gòu)的主要研究?jī)?nèi)容。其核心思想是按照給定的閾值從體數(shù)據(jù)中抽取空間曲面，然后再按照傳統(tǒng)方法進(jìn)行繪制。這類方法主要包括輪廓線法(Keppel，1975)，等值面法(Lorensen et al，1987)，幾何變換(Boissonnat，1984)等方法。提取等值面法由于其操作簡(jiǎn)單、精確度較高等原因一直是曲面重建研究的重要內(nèi)容。這種方法的核心思想是用戶給定一組值，算法根據(jù)這組值提取對(duì)應(yīng)的等值面。

地質(zhì)體接觸界面提取與其他地學(xué)數(shù)據(jù)等值面提取的主要區(qū)別是接觸界面提取的對(duì)象往往是一組地質(zhì)認(rèn)識(shí)信息而不是某一特征值，如，地球化學(xué)、地球物理等數(shù)據(jù)的分析值等。

綜上所述，在上述兩種曲面擬合算法思想的基礎(chǔ)上，提出基于反距離加權(quán)方法和移動(dòng)立方格(Marching Cubes，MC)算法的層狀三維地質(zhì)界面擬合算法。該算法的核心思想：首先由離散點(diǎn)進(jìn)行空間插值得到三維柵格模型;其次建立增量函數(shù)f(Δ z)=zIDW－z，求出最優(yōu)擬合曲面值 f(Δz)=ziso;最后，基于Marching Cubes算法和ziso追蹤得到等值面S，即為層狀地質(zhì)界面。

1 算法思想

1．1 算法基本概念

介紹一些相關(guān)的概念及術(shù)語定義以便后文論述使用。

(1)曲面擬合(Edward，2003)。所謂曲面擬合，就是根據(jù)實(shí)際試驗(yàn)測(cè)試數(shù)據(jù)，求取函數(shù)f(x，y，z)與變量x，y，z之間的解析式，使其所確定的曲面通過或近似通過所有的實(shí)驗(yàn)測(cè)試點(diǎn)。也就是說，是所有實(shí)驗(yàn)測(cè)試點(diǎn)擬合的曲面M'能夠合理逼近控制點(diǎn)所在曲面M。

(2)地理學(xué)第一定律(Tobler，1970)。任何事物之間均相關(guān)，而距離較近的事物總是比距離較遠(yuǎn)的事物影響更大。這里的距離不僅包含空間上的，也包含時(shí)間上的距離。

(3)移動(dòng)立方格算法(Lorensen et al，1987)。移動(dòng)立方格算法本質(zhì)是從一個(gè)三維數(shù)據(jù)場(chǎng)中提取出一個(gè)等值面，所以也稱為“等值面”算法。其算法思想是在規(guī)模大小為Nx×Ny×Nz的體數(shù)據(jù)場(chǎng)中，由每8個(gè)頂點(diǎn)組成1個(gè)體元，每個(gè)頂點(diǎn)都有1個(gè)數(shù)值，體元每條棱上的數(shù)值都是線性變化的。等值面的生成就是根據(jù)兩相鄰頂點(diǎn)的增量函數(shù)值生成等值點(diǎn)，然后將每個(gè)體元中的所有等值點(diǎn)按組合方式連接起來，生成指定值的等值面的過程。在標(biāo)量場(chǎng)曲面的增量函數(shù)建立完成后，可以使用上述方法追蹤地質(zhì)體曲面。

(4)K-近鄰查詢。給定查詢點(diǎn)和正整數(shù)K，從數(shù)據(jù)集中距離查詢點(diǎn)最近的K個(gè)數(shù)據(jù)，當(dāng)K=1時(shí)，為最近領(lǐng)域查詢。

1．2 算法主要思想

筆者提出基于反距離加權(quán)方法和移動(dòng)立方格算法的層狀三維地質(zhì)界面擬合算法。反距離加權(quán)方法假設(shè)未知點(diǎn)的值受較近控制點(diǎn)的影響比較遠(yuǎn)控制點(diǎn)影響更大。影響程度用點(diǎn)之間距離乘方的倒數(shù)表示。這種特性使得該方法可以應(yīng)用于模擬連續(xù)曲面，最經(jīng)典的例子是構(gòu)建DTM模型。在三維空間中，盡管空間插值的結(jié)果是規(guī)則網(wǎng)格，但其數(shù)學(xué)本質(zhì)不會(huì)發(fā)生變化，即較近控制點(diǎn)的影響比較遠(yuǎn)控制點(diǎn)影響更大。因此，如果待插值為高程值z(mì)，則空間插值的結(jié)果表示由采樣點(diǎn)擬合的曲面M。此時(shí)就可以使用移動(dòng)立方格算法提取曲面M'，使之能夠合理逼近控制點(diǎn)所在曲面M。

1．3 算法流程

如1．2節(jié)所述，總結(jié)筆者提出的算法總體流程如下。

步驟1：基于反距離加權(quán)(IDW)方法對(duì)地層采樣點(diǎn)的高程值z(mì)進(jìn)行空間插值。

步驟2：建立插值結(jié)果與待插點(diǎn)實(shí)際高程值z(mì)'的增量函數(shù)f(Δz)。

步驟3：選擇增量函數(shù)f(Δz)的某個(gè)值，作為曲面提取值。

步驟4：應(yīng)用移動(dòng)立方格(MC)法，提取f(Δz)的某個(gè)值處的等值面。

步驟5：計(jì)算曲面光照、設(shè)置顏色、重建曲面拓?fù)浣Y(jié)構(gòu)等。

綜上所述，根據(jù)算法的主要思想和總體流程，總結(jié)算法流程如圖1所示。

圖1 算法總體流程圖

2 算法實(shí)現(xiàn)

本算法主要包含兩個(gè)主要內(nèi)容：(1)快速三維空間插值算法;(2)移動(dòng)立方格算法。

在三維空間中，隨著點(diǎn)數(shù)量的增加，K-近鄰查詢非常耗時(shí)。例如，在三維空間中對(duì)12 000個(gè)點(diǎn)進(jìn)行反距離加權(quán)空間插值，設(shè)待插點(diǎn)達(dá)到1 000萬，如果使用線性方法進(jìn)行方法其時(shí)間復(fù)雜度將達(dá)到0(12 000×10 000 000)。可以看出，上述方法的時(shí)間復(fù)雜度是用戶無法接受的。

筆者提出使用八叉樹和R-Tree樹兩種空間索引結(jié)構(gòu)，完成三維散亂點(diǎn)云的K-近鄰查詢，提高三維空間插值算法的效率。

2．1 算法1：三維離散點(diǎn)快速空間插值算法

輸入：離散點(diǎn)集合 P{(x，y，z)|x，y，z為任意實(shí)數(shù)}

輸出：柵格體S

步驟1：建立三維點(diǎn)云P的八叉樹空間索引。步驟2：建立八叉樹空間索引的三維R-Tree樹索引。

步驟3：通過矩陣計(jì)算獲得待插點(diǎn)的搜索范圍。

步驟4：應(yīng)用反距離加權(quán)方法對(duì)地層分界點(diǎn)的高程值z(mì)進(jìn)行空間插值，得到S。

使用移動(dòng)立方格算法，提取地質(zhì)曲面。

在移動(dòng)立方格方法中，由于每個(gè)體元存在8個(gè)頂點(diǎn)，每個(gè)頂點(diǎn)可能有正和負(fù)兩種狀態(tài)(等值面經(jīng)過0狀態(tài))，因此，每個(gè)體元按照其頂點(diǎn)的正、負(fù)分布，共有28=256種不同的組合狀態(tài)。盡管判斷等值面將與哪些體元相交在理論上很容易理解，但是，根據(jù)這256種狀態(tài)求出每個(gè)體元的等值面是個(gè)相當(dāng)繁瑣且容易出錯(cuò)的過程。實(shí)際工作中，根據(jù)互補(bǔ)對(duì)稱性和旋轉(zhuǎn)對(duì)稱性可以將這256種狀態(tài)簡(jiǎn)化為15種情況(蔣先剛，2003)(圖2)。但是，簡(jiǎn)化的模型是一個(gè)存在二義性的模型，在文中采用漸近判別法(Nielson et al，1991)來消除二義性。

圖2 模型簡(jiǎn)化的15種情況

2．2 算法2：提取地質(zhì)分界曲面

輸入：柵格體S，待提取等值面的值z(mì)iso

輸出：三角網(wǎng)Surface

步驟1：查找包含值z(mì)iso的立方格。

步驟2：按照?qǐng)D2的方式生成三角面片。

步驟3：計(jì)算每一個(gè)三角形的法向量。

步驟4：結(jié)束。

3 實(shí)驗(yàn)與分析

使用西藏某礦區(qū)三維地層采樣數(shù)據(jù)作為實(shí)驗(yàn)數(shù)據(jù)，該礦區(qū)由鉆孔獲得的地層分界界面采樣數(shù)據(jù)是416個(gè)。實(shí)驗(yàn)環(huán)境、硬件開發(fā)環(huán)境：CPU主頻為2.16G Hz和內(nèi)存為1G的PC機(jī);軟件開發(fā)環(huán)境：MS Windows XP操作系統(tǒng)和Visual C++7．0。算法效果如圖3、圖4、圖 5、圖 6、圖7 所示。

4 結(jié)論

在已有曲面擬合算法思想的基礎(chǔ)上，筆者提出了基于反距離加權(quán)方法和移動(dòng)立方格算法的層狀三維地質(zhì)曲面擬合算法。本算法的優(yōu)點(diǎn)在于以數(shù)學(xué)原理與地理學(xué)第一定律為基礎(chǔ)，實(shí)現(xiàn)層狀地層的快速、精確擬合。該方法模型能夠適合絕大多數(shù)層狀地學(xué)三維曲面模擬。不足之處是不能解決不整合、構(gòu)造等復(fù)雜情況，這也是本算法下一步研究的方向。

圖3 地層接觸界面采樣點(diǎn)

圖4 柵格化結(jié)果

圖5 地層分界界面擬合結(jié)果

圖6 地層分界界面擬合與剖切結(jié)果對(duì)比

圖7 地質(zhì)界面與化探異常以及礦體復(fù)合顯示

毛先成．2006．三維數(shù)字礦床與隱伏礦體立體定量預(yù)測(cè)研究［D］．長(zhǎng)沙：中南大學(xué)．

蔣先剛．2003．三維數(shù)據(jù)場(chǎng)重構(gòu)與顯示工程軟件設(shè)計(jì)［M］．北京：萬水出版社．

BOISSONNAT J D．1984．Geometric structures for three-dimensional shape represent ation［J］．ACM Trans Graphics，3(4)：266－286．

EDWARD ANGLE．2003．Interactive Computer Graphics：A Top-Down Approach Using OpenGL［M］．Beijing：Higher Education Press．

HOPPE H，DEROSE T，DUCHAMP T，et al．1992．Surface reconstruction from unorganized points［J］．Computer Graphics，26(2)：71－78

KEPPEL E．1975．Approximating complex surfaces by triangulation of contour lines［J］．IBM Journal of Research and Development，19(1)：2 －11．

LORENSEN W E，CLINE H E．1987．Marching Cubes：A high Resolution 3D Surface Construction Algorithm ［J］．Computer Graphics，21(4)：163 －169．

NIELSON G M，HAMANN B．1991．The asymptotic decider：resolving the ambiguity in marching cubes［C］//Proceedings of Visulization'91，San Diego：CA，83－91．

TOBLER W．1970．A computer movie simulating urban growth in the Detroit region［J］．Economic Geography，46(2)：234－240．

Stratiform 3D geological interface fitting algorithm based on inverse distance weighted method

LI Nan1，LEI Ling-ling2，XIAO Ke-yan1，2

(1．Institute of Mineral Resources，Chinese Academy of Geological Sciences，Beijing 100037，China;2．School of Geosciences，Yangtze University，Jingzhou 434023，Hubei)

Geological curved surface modeling was an important research content in fields of mineral exploration and large-scale mineral assessment．Geological curved surface modeling was more complex，and of generality compared with that of other geological data like geochemical data or geophysical data．Geological curved surface modeling needed to solve many complex situations such as unconformity or faults．In response to these issues，the authors proposed an algorithm of fitting layered geological surface that was based on two algorithms of inverse distance weighted method(IDW)and marching cubes(MC)．The algorithms were characteristic of speedy，accurate fitting for stratified layers on the foundations of mathematical principles and the first theorem of geography，suited most 3D curved surface modeling of the layered earth sciences．

Geological interface;Spatial interpolation;Marching cubes algorithm;Iso-surface;Implicit geological modeling

TP391

A

1674－3636(2012)03－0291－05

10．3969/j．issn．1674-3636．2012．03．291

2012－06－06;編輯：陸李萍

國家高技術(shù)研究發(fā)展計(jì)劃(863)“礦產(chǎn)資源評(píng)價(jià)數(shù)字礦床模型及可視化技術(shù)研究”(2006AA06Z114)、國家科技支撐計(jì)劃“西部?jī)?yōu)勢(shì)礦產(chǎn)資源潛力評(píng)價(jià)技術(shù)及應(yīng)用研究”(2006BAB01A01)、礦產(chǎn)勘查三維預(yù)測(cè)評(píng)價(jià)信息平臺(tái)開發(fā)及應(yīng)用示范研究(1212010012013)、中央國家機(jī)關(guān)基本業(yè)務(wù)費(fèi)基金項(xiàng)目等基金資助

李楠(1980— )，男，副研究員，博士，主要從事地圖制圖與地理信息工程、大比例尺礦產(chǎn)資源評(píng)價(jià)研究，E-mail：superln1980@yahoo．com．cn