Ｋｒｉｇｉｎｇ算法在含水量三維屬性模型構建中的應用研究

摘 要：Kriging算法通過構造區(qū)域化變量的變異模型，并據(jù)此求取未知數(shù)據(jù)的最優(yōu)線性無偏估計量。結合構建工程地質(zhì)體含水量參數(shù)三維模型的應用，強調(diào)在三維環(huán)境下分析樣本的空間變異結構特征，以K-D樹建立樣本數(shù)據(jù)的空間索引快速搜索插值鄰域。針對插值過程中存在的負權值問題，采用線性規(guī)劃的方法加以解決。研究表明，該算法綜合考慮了各向異性以及不同變異尺度對算法的影響，比基于分層的二維Kriging方法更為合理。

關鍵詞：Kriging算法;各向異性;變異尺度;負權值; K-D樹;線性規(guī)劃;含水量

中圖分類號：TP391 文獻標志碼：A 文章編號：1001-3695(2008)08-2554-03

Research on application of Kriging in construction of

three-dimensional property model about water content

YAO Ling-qing1，PAN Mao1，QU Hong-gang1，LAN Xiang-rong1，CONG Wei-qing1，LIU Xue-qing2，YU Chun-lin2

（1.Key Laboratory of Orogenic Belts Crustal Evolution of Education Ministry，Peking University， Beijing 100871 China; 2. Beijing Bureau of Geology Mineral Prospecting， Beijing100050， China)

Abstract:Kriging algorithm builds varation model of regional variables and through which seeks for best linear unbiased estimator (BLUE) of unknown property. This paper combined with application of constructing 3D property model of water content within a certain engineering geological body， emphasized analyzing variation of regional variable under 3D circumstance， and adopted K-D tree to acquire neighborhoods quickly.Resolved negative weights during interpolation by linear program. The results show that 3D Kriging behaves better than 2D Kriging based on dividing layer for consideration of anisotropy and different scales on variation model.

Key words: Kriging algorithm;anisotropy;variation scale;negative weight; K-D tree;linear program;water content

0 引言

地質(zhì)統(tǒng)計學在礦產(chǎn)勘查、石油勘探、環(huán)境監(jiān)測、圖像處理等科學領域都有著廣泛的應用。Kriging插值方法作為地質(zhì)統(tǒng)計學中的核心算法，是一種綜合考慮變量的結構性和隨機性，對未采樣點屬性進行最優(yōu)無偏估計的估值方法。

迄今為止，Kriging方法在二維場合的應用非常普遍，相比之下三維環(huán)境下的實際應用較為少見。理論上地質(zhì)統(tǒng)計學可以在任意的完備度量空間進行，當然也就包括了三維歐氏空間。而之所以出現(xiàn)Kriging方法在三維空間中應用的不足，究其原因，筆者認為可以歸納為三點：a）二維形態(tài)的研究對象涵蓋了大量的應用場合，如研究煤層在分布區(qū)域上的厚度變化；b）對不同層位的研究區(qū)進行二維Kriging插值，形成三維空間的一系列截面，可間接地表達對象屬性分布的三維形態(tài)；c）計算機硬件條件的制約使三維地質(zhì)信息系統(tǒng)發(fā)展相對滯后，間接影響了地質(zhì)統(tǒng)計學在三維空間的應用。

利用Kriging方法構建研究對象的屬性場，相關的文獻根據(jù)應用領域的不同，研究的側重點也有所差異^[1~3]。文獻[1]以油層溫度場為研究對象，側重分析回歸模型和相關模型對Kriging算法的影響，變異結構以二維分析為主，插值結果通過曲面表示；文獻[2]通過Kriging方法構造二維的鉆孔溫度等深面，通過等深面疊加形成規(guī)則體元數(shù)據(jù)并進行體視化；文獻[3]利用Kriging方法建立水文地質(zhì)層三維模型，思路與文獻[2]類似，側重于方法的實現(xiàn)過程以及體視化的表達。

上述文獻在實現(xiàn)思路上仍然以二維空間的Kriging算法為主。本文以建立工程地質(zhì)體含水量參數(shù)三維模型的應用為契機，重點研究Kriging方法在三維空間中的拓展，深入探討隨之產(chǎn)生的變異結構各向異性特征、變異尺度對算法的影響、三維鄰域搜索算法以及插值運算中產(chǎn)生的負權值問題。

1 Kriging算法思想

Kriging算法的前提條件是構造區(qū)域化變量的變異模型，變異模型通過變異函數(shù)或協(xié)方差函數(shù)描述。算法的核心致力于滿足無偏和最小估計方差的要求，求取局部鄰域內(nèi)樣品的權值系數(shù)，數(shù)學上等價于一個求解條件極值的問題，通過附加一個拉格朗日乘數(shù)可列出求解權值的線性方程組，稱為Kriging方程組。方程組的系數(shù)矩陣借助協(xié)方差函數(shù)確立?？梢宰C明，在二階平穩(wěn)條件下，協(xié)方差函數(shù)和變異函數(shù)存在簡單的線性關系^[4]，因此Kriging方程組可以方便轉(zhuǎn)換為以變異函數(shù)值來表達系數(shù)矩陣。在協(xié)方差函數(shù)不存在的條件下，根據(jù)變異函數(shù)建立的Kriging方程組依然成立。基于變異函數(shù)在Kriging算法中的通用性，本文選擇變異函數(shù)為工具對區(qū)域化變量進行變異結構分析。

2 Kriging算法在含水量參數(shù)三維模型中的應用

實例的原始數(shù)據(jù)來源于工程地質(zhì)體鉆孔資料，該地質(zhì)體主要成分為粉粘土。地質(zhì)體分布為薄層狀結構，水平方向展布約16×12 km2，垂直方向的厚度僅有5~10 m，地勢走向基本上為西高東低。展布形態(tài)如圖1所示。

根據(jù)鉆孔資料取樣得到含水量參數(shù)的采樣數(shù)據(jù)（計量單位為百分數(shù)），利用Kriging插值結果可以構建出含水量參數(shù)在整個工程地質(zhì)體分布的三維模型。通過分析地質(zhì)體形態(tài)以及數(shù)據(jù)分布的特點，將地質(zhì)體離散化為分辨率128×128×128的規(guī)則體元數(shù)據(jù)。

2.1 采樣數(shù)據(jù)特點

采樣數(shù)據(jù)來源于工程鉆孔資料，共計48口鉆孔，取樣數(shù)共計345個。鉆孔分布在水平方向上分布稀疏，間隔在800~1 200 m，垂向上采樣較為稠密，間距在0.2~1 m。其分布形態(tài)如圖 2所示。參數(shù)分布近似為正態(tài)分布，符合Kriging插值方法應用的條件。

2.2 應用流程與關鍵問題

本節(jié)重點論述Kriging方法構建工程地質(zhì)體屬性三維模型的流程以及存在的關鍵問題。在此基礎上探討其與常見二維場合中Kriging方法的差異，并通過實驗進行對比與分析。

2.2.1 變異結構分析

綜合考慮方向和尺度因素對區(qū)域化變量變異性的影響，獲取一個相對合理的變異函數(shù)，稱為變異結構分析。

1) 各向異性分析 經(jīng)過對不同方向的采樣數(shù)據(jù)計算形成實驗變異函數(shù)（式（1）），發(fā)現(xiàn)實驗數(shù)據(jù)在水平方向上的變異表現(xiàn)為各向同性，在垂向上的變異趨勢與水平方向不同。

γ(h)=1 / (2N(h)) N(h)i=1[Z(xi)-Z(xi+h)]2（1）

2) 變異尺度分析 在進行變異分析的過程中務必要考慮變量在不同尺度上發(fā)生的變異，而不能籠統(tǒng)地采用一種尺度對變異結構進行衡量。這一點可以借助采樣理論^[5]加以解釋，低頻的采樣會丟失小尺度的變化。

本文將水平方向上的實驗變異函數(shù)滯后距設為800 m，距離容許范圍為200 m；垂直方向滯后距設為0.5 m，距離容許范圍為0.2 m；不同方向上的角度容差范圍都設為15°角。對兩個方向上的實驗變異函數(shù)采用球形理論模型擬合，擬合結果如圖 3、 4所示。擬合后的理論模型在水平方向上為

r(h)=13.51+12.01×sph(3.0×103 m)

垂直方向的理論模型為

r(h)=26.2×sph(4.48 m)

2.2.2 基于K-D樹索引搜索三維鄰域

K-D樹索引最初由Bentley提出^[6]，用于解決數(shù)據(jù)庫檢索的問題。K-D樹本質(zhì)上是二叉樹。具體而言，K維空間數(shù)據(jù)每個維度分別對應一個關鍵字作為空間劃分的標準，小于該關鍵字的歸類到左子樹，否則歸類到右子樹。依次按維度循環(huán)進行上述流程，直到葉子節(jié)點包含的數(shù)據(jù)小于指定的個數(shù)則剖分完成。圖 5為一簡單的二維空間K-D樹數(shù)據(jù)組織示意圖，圖中的矩形單元代表其中的節(jié)點。

文獻[7]是改進的K-D樹，選擇節(jié)點中對應劃分維度的關鍵字的中位數(shù)作為分界點，可以優(yōu)化索引的搜索效率。本文采用上述改進的K-D樹作為采樣數(shù)據(jù)鄰域搜索的索引。構建K-D樹的算法可參考文獻[7]。

插值的鄰域設置為以不同軸向的變程為半徑的長方體區(qū)域，窗口搜索的算法描述如下（C++語言）：

struct KDNode_Base

{enum {NODE， LEAF} tp;}//節(jié)點類型

struct KDNode : public KDNode_Base

{

size_t dividing_dim;//當前劃分維度

float dividing_val;//當前劃分界點

struct KDNode lchild;//左子節(jié)點

struct KDNode rchild;//右子節(jié)點

};

struct KDLeaf : public KDNode_Base

{Point* pnts;}//樣本數(shù)據(jù)

//鄰域搜索函數(shù)主體

Range_Query(KDNode_Base kdn， float lo[3]， float hi[3])

{

if(kdn->tp==KDNode_Base::LEAF){

//將葉子節(jié)點中滿足要求的數(shù)據(jù)提取出來

GetNeighbors(KDBLeaf(kdn)， float lo[3]， float hi[3]);

}else{KDNode tmp = (KDNode)kdn;

size_t dim = tmp->dividing_dim;

float dval = tmp->dividing_val;//遍歷左子節(jié)點

if(hi[dim]< dval‖lo[dim] dval)

Range_Query(tmp->lchild， lo， hi);//遍歷右子節(jié)點

if(hi[dim] >=dval || lo[dim] >=dval)

Range_Query(tmp->rchild， lo， hi);}

}

2.2.3 負權值問題

關于負權值的問題，已有相關文獻進行過研究，但在實際插值過程中卻常被忽略。Kriging算法應用的大量對象不容許出現(xiàn)負屬性值，如本文實例所研究的含水量參數(shù)。Kriging插值的負權值除了可能產(chǎn)生無效的負屬性值之外，還可能導致屬性值的畸高。

對于負權值的消除，簡單的方案可以將負權值系數(shù)置為零，然后重新計算權值^[4]，但缺乏合理的理論依據(jù)。Herzfeld^[8]提出基于Kuhn-Tucher條件利用二次規(guī)劃解決負權值問題的方案，但計算較為復雜。胡小榮^[9]提出利用線性規(guī)劃解決權值非負約束問題的方法，簡化了運算過程，但相對消元法求解計算復雜度依然偏高。

本文采用綜合的方法：僅在屬性估計值或期望方差為負時采用線性規(guī)劃的方案重新求解，而其他情況則利用消元法。由于上述情形并不常見，可以節(jié)省大量的計算時間。其合理性在2.3.3節(jié)詳細分析。

2.3 實驗

2.3.1 二維分層Kriging方法與三維Kriging方法的比較

為了對比二維分層Kriging方法與三維Kriging方法，本文利用原始數(shù)據(jù)交叉驗證對兩者插值的結果進行了比較。采用二維分層的Kriging方法誤差為[-20.77，13.49]，誤差絕對值均值為3.99； 三維Kriging方法誤差為[-16.65， 15.50]，誤差絕對值均值為3.42。顯然三維Kriging方法從誤差分布到誤差平均絕對值的大小都優(yōu)于二維分層的方法。進一步的比較結果如圖 6所示，可見兩者的估計誤差主要集中于(-5，5)，但是三維Kriging方法在較低誤差處有更高的百分比。

分析其原因，對于本實例數(shù)據(jù)，源于采樣密度的緣故，參數(shù)的變異在垂直方向上有更顯著的空間相關性，而在水平方向上卻有明顯的躍遷，存在塊金效應^[4]，即隨機性較強。按照垂直方向分層的方式對數(shù)據(jù)進行Kriging插值，由于未能考慮到相關性較強的垂向變化，不能充分體現(xiàn)Kriging算法根據(jù)待估計量與樣本相關程度分布權值的優(yōu)點。

少量數(shù)據(jù)點估計存在較大的偏差，達到10以上，經(jīng)過與原始數(shù)據(jù)的對比觀察發(fā)現(xiàn)，偏差較大的數(shù)據(jù)均位于邊界或采樣數(shù)據(jù)稀疏的地方，表明插值數(shù)據(jù)的構形對估計的結果有較大的影響，在采樣數(shù)據(jù)充分的前提下三維Kriging方法可以取得理想的結果。

2.3.2 K-D樹搜索鄰域與變異結構的關系

K-D樹在低維空間中有很好的查詢性能^[10]，因此在二維和三維空間中都可以作為離散數(shù)據(jù)索引組織的方式。插值的結果顯示，K-D樹能很好地應用于Kriging算法的鄰域搜索中。

值得注意的是，搜索鄰域范圍指定應參考變異函數(shù)的變程來確定。位于變程之外的數(shù)據(jù)由于缺少相關性，對估計點的貢獻不大，可排除在搜索范圍之外。這樣可減少Kriging方程組的元數(shù)，降低計算復雜度以及提高求解權值系數(shù)的穩(wěn)定性。

2.3.3 負權值處理方法分析

在利用消元法求解Kriging方程組的過程中，負權值出現(xiàn)的情況較為常見。通過采用簡單將負權值置零的方法或利用線性規(guī)劃的方法均可消除其影響。

經(jīng)進一步研究發(fā)現(xiàn)，負權值導致出現(xiàn)負的屬性估計值或期望方差的情況較為少見。負權值在大多數(shù)情況下不影響插值結果的合理性，相反在這種情況下采用其他方法消除負權值反而帶來更高的估計方差。而在消元法求取的插值結果無實際意義（估計值或方差小于0）的情況下，通過求取負權值處（消元法）估計值的期望方差，發(fā)現(xiàn)線性規(guī)劃求解法的方差普遍低于簡單方法，說明后者在負權值處理上更為合理。對比結果如表 1所示。

表 1 負權值消除方法對應方差

負權值處理方差最小值最大值平均值簡單置零（插值結果無實際意義時消除）53.63365.25560.413線性規(guī)劃（插值結果無實際意義時消除）22.51930.50428.205消元法不作處理（插值結果有意義）21.08754.54728.053線性規(guī)劃（插值結果有意義）21.21454.86828.732綜上所述，在常規(guī)方法求解估計值無意義的前提下，利用線性規(guī)劃消除負權值的影響對插值結果有明顯的改進，比簡單將負權值置零的方法表現(xiàn)出更低的方差。在權值都為正數(shù)的情況下，線性規(guī)劃和常規(guī)方法取得相同的權值系數(shù)^[9]。在其他情況下，線性規(guī)劃消除負權值反而會帶來略為更高的方差（插值結果相近），簡單置零的方法尤為顯著。

由于出現(xiàn)無意義插值結果的情況僅為少數(shù)，而線性規(guī)劃運算量高于消元法。將兩者綜合運用，僅在消元法求解不能滿足實際要求時采用線性規(guī)劃方法，可以節(jié)省大量計算時間。測試硬件配置為Mobile DualCore Intel Merom，1866 MHz，2 GB內(nèi)存，顯卡NVIDIA GeForce 7900 GS(512 MB)。本實例對比結果有:消元法為21 s；線性規(guī)劃法為101 s；綜合法為22 s.可見綜合法耗時與消元方法大致相同，但是比單純利用線性規(guī)劃法求解要節(jié)省約4倍的時間。

2.3.4 三維模型結果

綜合利用前文提到的方法及分析結果，通過Kriging算法插值建立了工程地質(zhì)體含水量的最終三維模型如圖 7所示。

為了探明含水量是否與地勢存在一定關系，在模型體視化時有意加大了縱向的拉伸比例（150：1）。從整體趨勢來看，含水量的變化與地勢的走向有一定關系。地勢較高的地方含水量較低，地勢低洼處含水量較高，呈現(xiàn)出西低東高的走勢。在水平方向上的變化呈區(qū)域性分布，如圖 8所示，代表了含水層不同物源遷移聚集的結果。

3 結束語

本文的研究表明，在三維環(huán)境下應用Kriging算法，可以充分發(fā)掘樣本數(shù)據(jù)的空間變異特征，比傳統(tǒng)二維分層的Kriging方法更具優(yōu)勢，在應用的過程中需要考慮各向異性以及觀測尺度對算法的影響。本文選擇K-D樹索引用于離散樣本數(shù)據(jù)的近鄰搜索，搜索的范圍參考不同方向的變程確定，對于常規(guī)的二維和三維數(shù)據(jù)均可以采用。Kriging方程組求解中存在的負權值問題常被忽略，本文進行了詳細的分析——利用線性規(guī)劃可以消除負權值的影響，在負權值足以引起無實際意義的屬性估計值或方差時，可以取得較好的結果；反之其合理性有待商榷，計算時間也大大增加。關于負權值的數(shù)學意義還有待進一步的研究。最終的插值結果顯示，三維Kriging方法用于三維工程參數(shù)屬性模型的構建可以取得較好的效果。

參考文獻:

[1]杜宇健， 蕭德云. Kriging 算法在溫度場計算中的應用分析[J]. 計算機輔助設計與圖形學學報， 2004，16(8):1153-1158.

[2]尚慶生，郭建文，李新，等.基于Kriging 插值的鉆孔地溫數(shù)據(jù)體視化[J].遙感技術與應用， 2006，21(4): 302-306.

[3]顏輝武，祝國瑞，徐智勇，等.基于Kriging 水文地質(zhì)層的三維建模與體視化[J].武漢大學學報:信息科學版，2004，29(7):611-614.

[4]孫洪泉.地質(zhì)統(tǒng)計學及其應用[M]. 徐州：中國礦業(yè)大學出版社， 1990.

[5]BENEDETTO J， PAULO J S， FERREIRA G. Modern sampling theory: mathematics and applications [M]. Boston: Birkhauser Boston， 2001.

[6]BENTLEY J L. Multidimensional binary search trees used for associative searching [J]. Communications of the ACM，1975，18(9):509-517.

[7]FRIEDMAN J H， BENTLEY J L， FINKEL R A. An algorithm for finding best matches in logarithmic expected time[J]. ACM Trans of Mathematical Software，1977，3(3): 209-226.

[8]HERZFELD U C. A note on programs performing Kriging with nonnegative weights[J]. Mathmatical Geology，1989，21(3):391-392.

[9]胡小榮. 一種考慮權值非負約束的克立格算法[J].地質(zhì)與勘探，1999，35(4): 29-32.

[10]CHEN Y S， HUNG Y P， YEN T F， et al. Fast and versatile al-gorithm for nearest neighbor search based on a lower bound tree [J]. Pattern Recognition，2007，40(2):360-375.

注：本文中所涉及到的圖表、注解、公式等內(nèi)容請以PDF格式閱讀原文