基于貝葉斯克里金的地下空間多源數(shù)據(jù)建模

李曉軍，李培楠，朱合華，劉 俊

（1.同濟大學 土木工程學院，上海 200092；2.同濟大學 巖土及地下工程教育部重點實驗室，上海 200092；3.同濟大學 測繪與地理信息學院，上海 200092）

隨著地下工程信息化水平的不斷提高，三維地質建模技術在地下工程中得到越來越多的應用［1］.借助成熟的三維地質建模技術和軟件，建立描述地質體的形態(tài)特征和構造要素空間關系的地下工程三維地質模型，不僅能全面地描述研究區(qū)域地層的空間幾何形態(tài)及其相互關系［2］，還能為地質體邊界限定、資源評估、地下構建筑物的規(guī)劃、隧道或巷道開挖以及后期的力學分析提供幫助［3］.為了減少地質模型的不確定性，應集成和融合多種地質勘探數(shù)據(jù)，如鉆孔、地質剖面以及礦業(yè)生產中的各種等高線信息來建立更加準確的地質模型［4］.傳統(tǒng)的多源地質建模方法主要是基于鉆孔、剖面和地形地質圖等數(shù)據(jù)采用地質統(tǒng)計學來建立地質模型［5］.例如Martelet等［6］利用鉆孔、剖面以及斷層邊界范圍數(shù)據(jù)來推斷地質界面并建立模型；Wu等［7］針對地質數(shù)據(jù)大多比較稀疏和低采樣率的特征，提出逐步細分的多源數(shù)據(jù)集成建模方法，對初始地質模型不斷修正；鐘登華等［8］利用多源地質數(shù)據(jù)建立了水利水電工程壩區(qū)的三維地質模型.上述研究將各種不同來源的信息離散化，在插值和建模時進行硬合成，實現(xiàn)多源信息的集成，以提高模型的精確度.然而上述方法的問題在于未考慮各種地質數(shù)據(jù)的不確定性差異，因而無法準確評價估計結果的不確定性.

為了更有效地利用不同類型和質量的地質數(shù)據(jù)，并保證其在地質建模過程中的一致性，本文以礦業(yè)生產為背景，考慮地質數(shù)據(jù)的不確定性程度，把勘探數(shù)據(jù)分為硬數(shù)據(jù)（如鉆孔、露頭和開挖信息）和軟數(shù)據(jù)（如剖面、地震、地質解譯圖件和專家經(jīng)驗等）［9］，利用區(qū)域變量理論研究這2類信息的空間結構變化特征，并采用貝葉斯克里金（BK）方法集成2種不同類型的數(shù)據(jù)來估計煤層表面高程，據(jù)此建立煤層三維地質模型，并通過與普通克里金（OK）方法的比較對該方法的正確性和有效性進行了驗證.

1 基于貝葉斯克里金方法的煤層建模

1.1 貝葉斯克里金基本原理

用z（x）表示區(qū)域D上的一個區(qū)域性變量，相應的隨機函數(shù)記為｛Z（x），x∈D｝，在這里Z（x）代表觀測數(shù)據(jù)所對應的隨機函數(shù)，簡稱硬數(shù)據(jù)；用M（x）表示區(qū)域D上的另一個區(qū)域性變量，相應的隨機函數(shù)記為｛M（x），x∈D｝，這里的M（x）代表猜測數(shù)據(jù)，簡稱軟數(shù)據(jù).

假定隨機函數(shù)Z（x），x∈D 的一組觀測集［Z（xi），i＝1，2，…，N］.由此定義一個隨機函數(shù)為

式中：μM（x）為 M（x）的數(shù)學期望值.對任意一組觀測值｛ZT（xi）＝Z（xi）-μM（xi），i＝1，2，…，N｝，BK估計希望得到的估計值Z＊（x0）具有如下形式：

式中：x0為區(qū)域D 中一點；λi（i＝1，2，…，N）為待定的加權系數(shù).

依據(jù)估計的無偏性和估計方差的最小性，利用拉格朗日方法得到關于權系數(shù)λi的BK方程組.

通過求解BK方程組，可以得到λi和β值.求得λi后，還可以進一步求得最小的估計誤差方差，即BK方差.

式中：Nh為數(shù)據(jù)對的個數(shù)，h為不同數(shù)據(jù)點之間的空間分離距離.γM（xi，xj）不僅依賴于xi和xj之間的距離，還依賴于它們本身的位置.因此，它應該從更大的一類函數(shù)族中去選取.Omre［10］證明了當取有如下形式時，它是條件正定的.

式中：σM（x）是定義于D 上的非零函數(shù)；γS（·）為標準的變異函數(shù)形式，且［1-γS（xi-xj）］是正定函數(shù).上式還可以改寫成

式中：CS（xi-xj）＝1-γS（xi-xj）.CS（·）可以看成是一個空間相關函數(shù)，由于通常假設｛M（x）；x∈D｝的一、二階矩是已知的，因此（x）就是軟數(shù)據(jù)的方差（不確定性的定量表示）.如果軟數(shù)據(jù)的（x）在整個區(qū)域D內是常量，那么這種情況就可以簡化成基于隨機殘余函數(shù)｛ZT（x），x∈D｝的普通克里金估計，也對應于減去預定義漂移的泛克里金方法.特別地，如果軟數(shù)據(jù)中不含任何不確定量，即方差Var［M（x）］＝0，x∈D，則γZ（xi，xj）＝γZ｜M（xixj），這就退化到了簡單克里金系統(tǒng)［11］.

1.2 數(shù)據(jù)網(wǎng)格化過程中的斷層處理

當研究區(qū)存在許多斷層時，構造非常復雜，在數(shù)據(jù)網(wǎng)格化的時候就要考慮到斷層的影響.目前常用的方法主要有：分塊法、斷層恢復法、斷面法、最小繞射距離權法和斷層軌跡法.這些方法中，斷層恢復法的效果較好且簡便，其基本思想是將地質體恢復到?jīng)]有被斷層切割前的狀態(tài)，對其直接插值，再疊加受斷層影響所產生的局部位移，最終得到真實的地層或煤層面網(wǎng)格數(shù)據(jù)，過程如下.

（1）把零位移線和斷層軌跡圍成的區(qū)域看作是地層面的移位部分，其他部分作為非移位部分.根據(jù)零位移線和斷層軌跡上的落差可以內插出在移位區(qū)域網(wǎng)格點上的落差值.非移位部分的網(wǎng)格點上的落差假定為零，由此得到垂直位移網(wǎng)格數(shù)據(jù)體.

（2）利用垂直位移網(wǎng)格數(shù)據(jù)把受斷層影響區(qū)域附近的層位復原到斷裂發(fā)生前的位置.用網(wǎng)格化的落差數(shù)據(jù)在位于地層面的移位部分的每個已知數(shù)據(jù)點上內插一個值，這個值可以認為是地層面在這一點上的落差.不在移位部分的每個數(shù)據(jù)點上的內插值為零.把所有的內插值從原始數(shù)據(jù)中減去就消除了由斷裂引起的落差，這相當于得到地層斷裂前的觀測數(shù)據(jù).顯然，只有受斷層影響的數(shù)據(jù)才被替換.

（3）建立復原層面的網(wǎng)格數(shù)據(jù)體.如果位移網(wǎng)格數(shù)據(jù)是正確的，即數(shù)據(jù)點處的層位已復原到斷裂前的位置上，可用不考慮斷層影響的常規(guī)網(wǎng)格化技術求得復原地層的網(wǎng)格數(shù)據(jù)體.即便某些斷塊上只有很少的觀測數(shù)據(jù)或根本沒有觀測數(shù)據(jù)，也可以進行網(wǎng)格化，因為這時已沒有斷層的影響了.

（4）把落差網(wǎng)格和復原地層網(wǎng)格相加，即可得到實際地層的網(wǎng)格數(shù)據(jù)體.

2 煤層表面估計與結果分析

2.1 工程地質條件

以安徽淮南望峰崗煤礦地下采掘工程為背景，主要利用鉆孔和地震數(shù)據(jù)建立煤層地質模型.該礦共含煤37層，平均總厚40.55m.其中，可采煤層共17層，主要可采煤層12層.井田含有38條斷層，其中逆斷層5條，正斷層33條；落差大于等于100m的7條，100～50m的3條，50～30m的5條，落差在30m以下的23條.井田范圍內總共有37層煤層，其平均總厚度為40.55m.考慮到煤層厚度、深度以及可持續(xù)開采程度的限制，可采煤層有11層，分別為：A1，A3，B4b，B6，B7，B9b，B10，B11b，C13和C15.井田勘測數(shù)據(jù)主要包括99個原始勘探鉆孔，形成勘探剖面23個.每條勘探剖面上有相應的煤層分布信息，以煤層界面信息為主.

本文所建模型在勘探線VI以北地區(qū)，面積為3.6km×1.8km，其中共有30個鉆孔數(shù)據(jù)、11條勘探剖面圖以及11層煤層底板等高線圖.其中煤層底板高程等值線是由地震勘探的解譯信息推斷所繪，地質剖面圖是基于地震解譯生成，并由鉆孔數(shù)據(jù)修正.其中包含的信息主要有煤層與地層的交界面和斷層線，在三維地質建模中將作為主要的約束信息參與地層模型的插值.其中，C13煤層以其良好的開采條件作為該煤礦的首采煤層，其中煤厚分布相對較平穩(wěn)，在7～9m附近波動，且埋藏深度較淺.因此，以C13煤層底板高程來驗證BK算法的可靠性.

2.2 基于貝葉斯克里金方法的煤層底板高程估計

硬數(shù)據(jù)是指在確切位置上的精確地質信息，如鉆孔數(shù)據(jù).軟數(shù)據(jù)是指那些可能的或模糊的地質采樣信息，如剖面圖和底板等高線.這些原始的輸入數(shù)據(jù)可以用如下3種類型的數(shù)學形式來表達：硬數(shù)據(jù)，Z（xi）＝Zi，其中Z（xi）為空間任意一點xi處的屬性值；間隔型軟數(shù)據(jù)，ai≤M（xi）≤bi；分布型軟數(shù)據(jù)，F(xiàn)xi（m）＝P ［M （xi）≤m ］，其中M（xi）為空間屬性值所代表的隨機變量，［ai，bi］為間隔型軟數(shù)據(jù)的變化區(qū)間，F(xiàn)xi（m）為分布型軟數(shù)據(jù)的累積頻率，P為概率分布.因此，不論是間隔型軟數(shù)據(jù)，還是分布型軟數(shù)據(jù)，都可以利用σM（x）表示地質信息的離散程度并對其不確定性進行編碼.

為了合理地利用BK方法，在估計之前，地質剖面和地震等高線所代表的軟數(shù)據(jù)需要離散化，考慮計算效率，其離散間距尺寸為200m×200m，而輸出網(wǎng)格尺寸為20m×20m.以煤層C13底板高程數(shù)據(jù)為例，剔除部分C13煤層缺失的鉆孔資料，其中可以利用的硬數(shù)據(jù)總數(shù)為26個，軟數(shù)據(jù)總數(shù)為238個，共有264個軟硬數(shù)據(jù)參與估計，軟硬數(shù)據(jù)的空間分布情況如圖1所示，其中，NW-SE表示西北-東南方向，SW-NE 表示西南-東北方向.

圖1 投影到水平面上的數(shù)據(jù)空間位置Fig.1 Map of data locations projected onto a horizontal plane

表1給出了不同數(shù)據(jù)源下的C13煤層底板高程數(shù)據(jù)的整體概率統(tǒng)計指標.硬數(shù)據(jù)所代表的高程均值為-809.54m，而軟數(shù)據(jù)所代表的高程均值為-931.09m.對比圖1可看出鉆孔主要分布在該圖的中間偏上部位，這也與該處煤層深度較淺、方便鉆井施工有關.通過高程的標準差和變化范圍可以看出軟數(shù)據(jù)具有較大的空間變異性和不確定性.

表1 高程的軟硬數(shù)據(jù)統(tǒng)計指標參數(shù)Tab.1 Statistics for hard and soft data of coal floor elevation

為了求解BK線性方程組，需要建立2種變異函數(shù)模型，即利用軟數(shù)據(jù)得到的變異函數(shù)模型和綜合軟硬數(shù)據(jù)的條件變異函數(shù)模型.由于變異函數(shù)的確定是求解BK方程的首要問題，而條件變異函數(shù)的獲取首先需要得到軟數(shù)據(jù)所代表的隨機函數(shù)的數(shù)學期望，因此，首先利用一個傾向東北、傾角為25°的二階多項式曲面方程來擬合軟數(shù)據(jù)一階矩μM（x）分布，其代表了軟數(shù)據(jù)隨位置不同而變化的期望函數(shù)，如圖2所示.具體的軟數(shù)據(jù)所代表的隨機函數(shù)的期望函數(shù)μM（x）的二次曲面函數(shù)擬合公式為

圖2 基于軟數(shù)據(jù)的C13煤層底板高程云圖等高線Fig.2 Map of C13coal seam floor elevation contours used as soft data for the model

求出了μM（x），還需要得到σM（x）（一般通過主觀假設和專家經(jīng)驗獲得），才能完整地求解γM（xi，xj），γZ｜M（xi-xj）以及BK 方程組，而對于該煤層地質模型，由式（6）和式（7）可以看出，剖面及地震數(shù)據(jù)的不確定性是通過在γM（xi，xj）中的σM（x）來體現(xiàn)的，σM（x）表達了不同位置處的剖面及地震數(shù)據(jù)的離散性，即不確定性.一般來說，由于地質剖面大多數(shù)依賴地質工程師的專家經(jīng)驗和鉆孔之間的插值，地震數(shù)據(jù)解譯過程中的速度遷移模型、時深轉換以及分辨率等因素都會帶來解譯結果的不確定性，且深度越大不確定性越大，因此假定σM（x）隨著深度線性增加，如表2所示.

表2 分配給特定深度處軟數(shù)據(jù)的標準差Tab.2 Standard deviations assigned to the soft data

考慮到煤層在走向和傾向上的區(qū)域化變量的差異比較明顯，因此建立各向異性的變異函數(shù)模型來對區(qū)域化變量的結構特征進行描述.這2個方向上的變異函數(shù)理論模型的參數(shù)是采用加權多項式回歸方法來擬合采樣變異函數(shù)所得到，最終通過交叉驗證來選取最合適的理論變異函數(shù)模型及其參數(shù)［12］.以26個鉆孔資料為驗證數(shù)據(jù)源，以平均誤差（ME）、均方差（MSE）、平均克里金方差（MKV）和標準化克里金方差（SKV）為統(tǒng)計學指標，其判斷標準為：ME盡可能接近于零；MSE和MKV盡可能小，且兩者盡量接近；SKV盡可能接近于1［12］.表3給出了變異函數(shù)各擬合模型指標參數(shù)的比較.

表3 理論變異函數(shù)模型的精確度評價Tab.3 The assessment on precision of theoretical model

通過計算，軟數(shù)據(jù)的各向異性變異函數(shù)的套合結構模型（圖3a）最終用球狀模型描述，其中變異函數(shù)的塊金值為542.5m2，基臺值為1710.2m2，NW-SE方向變程為1611.988m，SW-NE方向變程為829.565m，各向異性比為1.95.軟數(shù)據(jù)的各向異性變異函數(shù)結構套合模型為

集成軟硬數(shù)據(jù)的條件各向異性變異函數(shù)模型（圖3b）的塊金值為5985m2，基臺值為23488m2，NW-SE方向變程為1542.154m，SW-NE方向變程為873.043m，各向異性比為1.76.集成軟硬數(shù)據(jù)的條件各向異性變異函數(shù)結構套合模型為

各向異性分離距離h的套合公式如下所示，其中hNW-SE為西北-東南的分離距離，hSW-NE為西南-東北的分離距離，K為各向異性比.

圖3 各向異性變異函數(shù)模型Fig.3 Sample and modeled variogram along main anisotropy directions of normal

利用BK方法進行估計的C13煤層底板高程和估計標準差分布情況如圖4所示.在加上了地震數(shù)據(jù)之后估計誤差標準差在鉆孔附近偏小，遠離鉆孔的地方標準差較大，不確定性在遠離鉆孔的研究區(qū)域底部和右上角達到最大.BK方法并不是把軟數(shù)據(jù)簡單地離散作為硬數(shù)據(jù)的補充.同時，在鉆孔數(shù)據(jù)稀少的右下角區(qū)域，該等值線主要反映了地震數(shù)據(jù)的特征，由此認為，BK方法克服了鉆孔數(shù)據(jù)過少帶來的估計失真和外推精度較差的問題.

3 估計方法的比較與驗證

3.1 估計方法的比較

為了探討B(tài)K方法在集成多源數(shù)據(jù)上的優(yōu)勢，將BK估計結果與軟數(shù)據(jù)直接離散化后加入硬數(shù)據(jù)中的OK估計結果進行比較，以ME，MSE和皮爾遜相關系數(shù)r作為評價指標，評價結果如表4所示.

同樣，基于26個鉆孔數(shù)據(jù)，此處的交叉驗證法所得到的BK方法估計的結果與表3中球狀模型的ME和MSE指標相同，這是因為本文通過球狀變異函數(shù)模型來進行空間區(qū)域性變量的估計.然后，再利用這26個測量值與估計值之間的誤差作頻率直方圖，其擬合曲線如圖5所示.

圖4 C13煤層底板高程BK估計結果Fig.4 The estimation map and standard deviation map in the coal seam C13obtained by BK procedure

表4 2種方法的整體定量比較標準Tab.4 Quantitative comparison criteria of these two methods in the overall area

圖5 BK和OK的估計誤差分布直方圖Fig.5 Histograms of estimation errors for BK and OK

由表4可見，無論是平均誤差還是均方誤差，OK方法都比BK方法大，因此可以認為，相對于BK方法，OK方法估計值的不確定性更大.這是因為OK方法雖然簡單地集成了軟硬數(shù)據(jù)，但沒有考慮軟硬數(shù)據(jù)的不確定性程度區(qū)別，在軟數(shù)據(jù)存在較大的測量、解譯以及推斷誤差的情況下，即在深部地震解譯信息的處理中，如果沒有給出合適的誤差波動范圍，估計結果會和鉆孔或開采巷道處的真實信息存在較大的偏差.由交叉驗證生成的整體估計誤差分布如圖5所示.從圖5中2條誤差頻率直方圖擬合曲線可以看出，利用OK方法所產生的誤差分布主要集中于［-11.25，13.75］，利用BK方法所產生的誤差分布主要集中于［-7.00，9.00］，后者誤差分布范圍更集中.由此說明，OK估計方法會產生較大的誤差離散程度，導致其不確定性更大，對于其估計結果的精確度更難控制.而BK方法能在較少硬數(shù)據(jù)的條件下，通過合理利用軟數(shù)據(jù)作為補充，給出更為精確的估計結果.

3.2 估計方法的驗證

為了更進一步驗證BK方法的有效性，基于煤層開采中所獲得的局部實測信息，即煤巷中的見煤點和地質素描數(shù)據(jù)，利用OK和BK方法分別進行對比驗證.煤層估計值與實測值比較曲線見圖6a，估計值與實測值之間的偏差見圖6b.其中實測曲線由首采區(qū)工作面附近煤巷中的見煤點數(shù)據(jù)擬合而成，隆起部位是由于鄰近區(qū)域內有鉆孔.

圖6 C13煤層底板高程不同估計方法的局部驗證結果Fig.6 The validation of the estimation results about the coal seam C13

對比2條估計曲線（圖6a），BK預測相對于OK預測更加接近于實際情況.OK方法把軟數(shù)據(jù)當成硬數(shù)據(jù)，嚴格地通過每一個軟數(shù)據(jù)點，由于在局部采區(qū)該深度范圍內，相對于鉆孔數(shù)據(jù)來說地震數(shù)據(jù)平均有2～3m的誤差，因此把軟數(shù)據(jù)直接作為硬數(shù)據(jù)的補充有一定的局限性，所產生的偏差比BK方法更大，其中OK方法平均偏差為0.72m，BK方法平均偏差為0.28m（圖6b）.相反，BK方法對地震數(shù)據(jù)賦予一定的不確定性，并在估計中加以考慮，所得結果與實際值更加吻合.2種方法局部估計偏差的評價指標如表5所示.OK方法的平均誤差和均方誤差的絕對值都比BK方法大，因此同樣可以認為，在局部尺度下BK方法相對于OK方法能夠提供更加準確的估計結果.在隧道的兩端，BK方法估計結果與真實值的相關程度更大，這表明在鉆孔較稀少的區(qū)域，BK方法能結合地震數(shù)據(jù)得到更合理的外推結果.

表5 局部估計偏差的定量比較標準Tab.5 Quantitative comparison criteria of the discrepancies between the observations and the estimation in the local area

4 煤層三維模型的構建

由于C13煤層被4條大型斷層（落差為15～30 m）、2條中型斷層以及3條小型斷層所切割，因此，在數(shù)據(jù)網(wǎng)格化過程中需要利用相應的斷層處理技術來處理斷層對空間變量的影響.即：①建立地層面的垂直位移網(wǎng)格；②用已得到的垂直位移數(shù)據(jù)對數(shù)據(jù)點作內插，使數(shù)據(jù)點處的層面復原到斷裂發(fā)生前的位置；③用復原后的數(shù)據(jù)建立一個網(wǎng)格數(shù)據(jù)體；④把2個網(wǎng)格相加，即生成了考慮斷層影響的網(wǎng)格數(shù)據(jù)體.本文采用BK方法完成步驟③之后，再加上由斷層引起的垂直位移，最終通過約束Delaunay算法生成一系列不規(guī)則三角網(wǎng)，如圖7所示.圖8為最后利用煤層上下表面三角網(wǎng)的拓撲對應關系生成的煤礦區(qū)域三維地質模型.

5 結論

圖7 C13煤層底板的不規(guī)則三角網(wǎng)Fig.7 A TIN mesh of C13coal floor surface

圖8 3D煤層地質模型Fig.8 The 3Dgeological model of the coal seams

通過BK方法使鉆孔數(shù)據(jù)和地震數(shù)據(jù)有機結合，而不是單純從地震剖面上獲取地層在平面各點處的深度，以此彌補了鉆孔采樣數(shù)據(jù)的不足.其估計結果在硬數(shù)據(jù)比較稀少的地區(qū)主要反映了軟數(shù)據(jù)的特征，而在硬數(shù)據(jù)比較稠密的地區(qū)主要取決于硬數(shù)據(jù).該方法不僅有效地集成了鉆孔和地震數(shù)據(jù)，還克服了稀少鉆孔數(shù)據(jù)的外推缺陷.通過與把軟數(shù)據(jù)直接硬化后加入待估數(shù)據(jù)集的OK方法進行對比，BK方法在合理利用軟硬數(shù)據(jù)方面具有的優(yōu)越性和合理性表現(xiàn)如下：

（1）地震深度數(shù)據(jù)是描述地質體空間分布特征的信息，它在反映地層界面在一個地區(qū)的整體變化趨勢方面具有比較好的效果，但就局部范圍而言，其精度又較差，因此把這2類數(shù)據(jù)不加區(qū)別地混合在一起使用，必然會造成較大的誤差.

（2）與OK方法相比，BK方法能夠得到更小的平均誤差和均方誤差，所以在集成大范圍精度較差的地震數(shù)據(jù)工作中BK方法能夠給出更加令人滿意的結果.

（3）BK方法的優(yōu)勢在于能夠集成直接觀測信息和含有不確定性的全局性先驗信息，利用σM（x）合理地分配軟數(shù)據(jù)的全局不確定性，并基于γM（xi，xj）在插值過程中進行計算.在觀測點密集的區(qū)域，直接觀測數(shù)據(jù)對其計算結果的影響較大，而在遠離觀測集的區(qū)域內，軟數(shù)據(jù)及其不確定性σM（x）對最終計算結果也有顯著的影響.因此BK方法給出的估計結果更加符合地質數(shù)據(jù)的實際情況，且在鉆孔數(shù)據(jù)稀少且存在非精確地震數(shù)據(jù)作為補充的情況下具有較為廣闊的應用前景.

［1］陳學習，吳立新，車德福，等.基于鉆孔數(shù)據(jù)的含斷層地質體三維建模方法［J］.煤田地質與勘探，2005，33（5）：5.CHEN Xuexi，WU Lixin，CHE Defu，et al.3D modeling method of geological bodies including faults based on borehole data［J］.Coal Geology ＆Exploration，2005，33（5）：5.

［2］姜在炳.煤層動態(tài)建模技術及應用［J］.煤炭學報，2006，31（1）：40.JIANG Zaibing.Dynamic modelling technology of coal seam and its application［J］.Journal of China Coal Society，2006，31（1）：40.

［3］Caumon G，Collon-Drouaillet P，le Carlier de Veslund，et al.Surface-based 3D modeling of geological structures ［J］.Mathematical Geosciences，2009，41（8）：927.

［4］Kaufmann O， Martin T.3D geological modelling from boreholes，cross-sections and geological maps，application over former natural gas storages in coal mines［J］.Computers ＆Geosciences，2008，34（3）：278.

［5］李曉軍，胡金虎，朱合華，等.基于Kriging方法的煤層厚度估計及三維煤層建模［J］.煤炭學報，2008，33（7）：765.LI Xiaojun，HU Jinhu，ZHU Hehua，et al.The estimation of coal thickness based on Kriging technique and 3D coal seam modeling［J］.Journal of China Coal Society，2008，33（7）：765.

［6］Martelet G，Calcagno P，Gumiaux C，et al.Integrated 3D geophysical and geological modelling of the Hercynian Suture Zone in the Champtoceaux area（south Brittany，F(xiàn)rance）［J］.Tectonophysics，2004，352：117.

［7］Wu Q，Xu H，Zou X.An effective method for 3D geological modeling with multi-source data integration［J］.Computers ＆Geosciences，2005，31（1）：35.

［8］鐘登華，李明超，楊建敏.復雜工程巖體結構三維可視化構造及其應用［J］.巖石力學與工程學報，2005，24（4）：575.ZHONG Denghua，LI Mingchao，YANG Jianmin.3D visual construction of complex engineering rock mass structure and its application［J］.Chinese Journal of Rock Mechanics and Engineering，2005，24（4）：575.

［9］高美娟，朱慶忠，張淑華.利用貝葉斯-克里金估計技術進行儲層參數(shù)預測［J］.石油地球物理勘探，1999，34（4）：390.GAO Meijuan，ZHU Qingzhong，ZHANG Shuhua.Reservoir parameter prediction using Bayesian kriging estimation technique［J］.Oil Geophysical Prospecting，1999，34（4）：390.

［10］Omre H.Bayesian kriging merging observations and qualified guesses in kriging［J］.Mathematical Geology，1987，19：25.

［11］Nobre M M，Sykes J F.Application of Bayesian kriging to subsurface characterization［J］.Canadian Geotechnical Journal，1992，29（4）：589.

［12］張仁鐸.空間變異理論及應用［M］.北京：科學出版社，2005.ZHANG Renduo.The theory and application of spatial variability［M］.Beijing：Science Press，2005.