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

（1.同濟(jì)大學(xué) 土木工程學(xué)院，上海 200092；2.同濟(jì)大學(xué) 巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200092；3.同濟(jì)大學(xué) 測繪與地理信息學(xué)院，上海 200092）

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

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

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

1.1 貝葉斯克里金基本原理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 煤層表面估計(jì)與結(jié)果分析

2.1 工程地質(zhì)條件

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

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

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

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

為了合理地利用BK方法，在估計(jì)之前，地質(zhì)剖面和地震等高線所代表的軟數(shù)據(jù)需要離散化，考慮計(jì)算效率，其離散間距尺寸為200m×200m，而輸出網(wǎng)格尺寸為20m×20m.以煤層C13底板高程數(shù)據(jù)為例，剔除部分C13煤層缺失的鉆孔資料，其中可以利用的硬數(shù)據(jù)總數(shù)為26個(gè)，軟數(shù)據(jù)總數(shù)為238個(gè)，共有264個(gè)軟硬數(shù)據(jù)參與估計(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)計(jì)指標(biāo).硬數(shù)據(jù)所代表的高程均值為-809.54m，而軟數(shù)據(jù)所代表的高程均值為-931.09m.對(duì)比圖1可看出鉆孔主要分布在該圖的中間偏上部位，這也與該處煤層深度較淺、方便鉆井施工有關(guān).通過高程的標(biāo)準(zhǔn)差和變化范圍可以看出軟數(shù)據(jù)具有較大的空間變異性和不確定性.

表1 高程的軟硬數(shù)據(jù)統(tǒng)計(jì)指標(biāo)參數(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ù)所代表的隨機(jī)函數(shù)的數(shù)學(xué)期望，因此，首先利用一個(gè)傾向東北、傾角為25°的二階多項(xiàng)式曲面方程來擬合軟數(shù)據(jù)一階矩μM（x）分布，其代表了軟數(shù)據(jù)隨位置不同而變化的期望函數(shù)，如圖2所示.具體的軟數(shù)據(jù)所代表的隨機(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）（一般通過主觀假設(shè)和專家經(jīng)驗(yàn)獲得），才能完整地求解γM（xi，xj），γZ｜M（xi-xj）以及BK 方程組，而對(duì)于該煤層地質(zhì)模型，由式（6）和式（7）可以看出，剖面及地震數(shù)據(jù)的不確定性是通過在γM（xi，xj）中的σM（x）來體現(xiàn)的，σM（x）表達(dá)了不同位置處的剖面及地震數(shù)據(jù)的離散性，即不確定性.一般來說，由于地質(zhì)剖面大多數(shù)依賴地質(zhì)工程師的專家經(jīng)驗(yàn)和鉆孔之間的插值，地震數(shù)據(jù)解譯過程中的速度遷移模型、時(shí)深轉(zhuǎn)換以及分辨率等因素都會(huì)帶來解譯結(jié)果的不確定性，且深度越大不確定性越大，因此假定σM（x）隨著深度線性增加，如表2所示.

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

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

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

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

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

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

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

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

3 估計(jì)方法的比較與驗(yàn)證

3.1 估計(jì)方法的比較

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

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

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

表4 2種方法的整體定量比較標(biāo)準(zhǔn)Tab.4 Quantitative comparison criteria of these two methods in the overall area

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

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

3.2 估計(jì)方法的驗(yàn)證

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

圖6 C13煤層底板高程不同估計(jì)方法的局部驗(yàn)證結(jié)果Fig.6 The validation of the estimation results about the coal seam C13

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

表5 局部估計(jì)偏差的定量比較標(biāo)準(zhǔn)Tab.5 Quantitative comparison criteria of the discrepancies between the observations and the estimation in the local area

4 煤層三維模型的構(gòu)建

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

5 結(jié)論

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

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

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

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

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

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

［1］陳學(xué)習(xí)，吳立新，車德福，等.基于鉆孔數(shù)據(jù)的含斷層地質(zhì)體三維建模方法［J］.煤田地質(zhì)與勘探，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］姜在炳.煤層動(dòng)態(tài)建模技術(shù)及應(yīng)用［J］.煤炭學(xué)報(bào)，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ì)及三維煤層建模［J］.煤炭學(xué)報(bào)，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］鐘登華，李明超，楊建敏.復(fù)雜工程巖體結(jié)構(gòu)三維可視化構(gòu)造及其應(yīng)用［J］.巖石力學(xué)與工程學(xué)報(bào)，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］高美娟，朱慶忠，張淑華.利用貝葉斯-克里金估計(jì)技術(shù)進(jìn)行儲(chǔ)層參數(shù)預(yù)測［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］張仁鐸.空間變異理論及應(yīng)用［M］.北京：科學(xué)出版社，2005.ZHANG Renduo.The theory and application of spatial variability［M］.Beijing：Science Press，2005.