序貫指示模擬法在三維地質(zhì)模型中的應(yīng)用
——以西遼河平原為例

楊明武 龍文華

1.中煤科工集團(tuán)北京華宇工程有限公司 北京 100120

2.廣東省地質(zhì)環(huán)境監(jiān)測(cè)總站 廣東 廣州 510510

地質(zhì)體三維建模這一概念由加拿大的學(xué)者Simon W.Houdling在1993年提出后，近年來(lái)隨著芯片技術(shù)的飛速發(fā)展，計(jì)算機(jī)軟件運(yùn)算能力大幅提高，因而，計(jì)算機(jī)圖形學(xué)理論的日趨漸完善，地學(xué)領(lǐng)域的地質(zhì)專家對(duì)三維地質(zhì)結(jié)構(gòu)模型的研究和計(jì)算機(jī)領(lǐng)域的研發(fā)專家對(duì)三維可視化的研究也取得了極大的突破。目前,國(guó)內(nèi)外專家開(kāi)發(fā)出了諸多的地學(xué)方面商業(yè)化應(yīng)用軟件,如:靈圖的VRMap系統(tǒng)（中國(guó)）、加拿大金康公司開(kāi)發(fā)的GemCom軟件、美國(guó)考格尼塞斯公司的GeoSec3D軟件、T-surf公司的GOCAD軟件（法國(guó)）、以色列帕拉代姆地球物理公司的GeoDepth系列軟件、美國(guó)CTech公司開(kāi)發(fā)的環(huán)境地質(zhì)可視化系統(tǒng)(EVS&MVS)及地下水模擬系統(tǒng)GMS(美國(guó))等。從這些軟件的使用功能來(lái)看，它們大多都能利用區(qū)域地震資料、場(chǎng)區(qū)鉆井資料和地表地形資料生成三維地質(zhì)模型解釋成果,且各有各的特色。但是以上軟件多針對(duì)油氣藏、物探及采礦行業(yè)，僅EVS&MVS和GMS在水文地質(zhì)、環(huán)境地質(zhì)三維建模方面較為成熟，但應(yīng)用上多集中于局部范圍內(nèi)的研究，對(duì)于區(qū)域松散層的建模研究較少。

GOCAD軟件作為新一代地質(zhì)建模軟件的代表,具有通用地質(zhì)模型所設(shè)想的功能——全三維、全拓?fù)洹⑹褂每臻g插值算法、擁有完善的地質(zhì)統(tǒng)計(jì)分析功能等特征，達(dá)到了半智能化地質(zhì)建模的世界最高水平,在國(guó)外受到眾多用戶青睞[1]。

由于在大區(qū)域且存在復(fù)雜相變的松散沉積盆地的三維模型構(gòu)建時(shí)，用非常有限的離散點(diǎn)數(shù)據(jù)進(jìn)行數(shù)學(xué)插值的方法描述空間巖性分布格局時(shí)會(huì)產(chǎn)生較大的模型失真。本文試圖應(yīng)用GOCAD中的序貫指示模擬功能對(duì)西遼河平原地層結(jié)構(gòu)建立三維地質(zhì)模型，探索一種能夠真實(shí)反應(yīng)在大區(qū)域盆地內(nèi)構(gòu)建三維地質(zhì)模型有效的方法，同時(shí)對(duì)西遼河平原的地質(zhì)、水文地質(zhì)結(jié)構(gòu)提供可視化的三維認(rèn)識(shí)，為各含水層組的分布、變化特征的深入研究提供了依據(jù)，為水文地質(zhì)概念模型的建立奠定了基礎(chǔ)。

1 序貫指示模擬原理

序貫?zāi)M是地質(zhì)統(tǒng)計(jì)學(xué)的一個(gè)重要組成部分，是將序貫思路與克立格插值相結(jié)合的條件模擬方法，將沿著隨機(jī)路徑序貫求出各地質(zhì)網(wǎng)格結(jié)點(diǎn)的條件累積分布函數(shù)值作為主要開(kāi)發(fā)思路，并從中取得數(shù)值模擬值。常用的有序貫高斯模擬和序貫指示模擬。序貫高斯模擬主要適用于滿足高斯分布的連續(xù)型數(shù)據(jù)場(chǎng)，而序貫指示模擬則主要適用于離散型數(shù)據(jù)場(chǎng)，也可用于離散化的連續(xù)型數(shù)據(jù)場(chǎng)，其基本原理為[2,3]：

利用線性無(wú)偏估計(jì)克里格法,來(lái)估計(jì)條件概率的公式為：

基于上，由于K閾值已經(jīng)給定，那么對(duì)于任何一個(gè)待估位置，各個(gè)閾值都有一個(gè)方程組和其對(duì)應(yīng)，在變量Z(x)的變化范圍內(nèi)，用K閾值對(duì)該估計(jì)范圍離散化，通過(guò)求解每個(gè)閾值對(duì)應(yīng)的克里格方程組，可以得到每個(gè)閾值下對(duì)應(yīng)的累積分布函數(shù)F{Zk,x|(n)}，然后，采用Monto Carlo法來(lái)求得隨機(jī)函數(shù)Z(x)在該位置的一個(gè)具體實(shí)現(xiàn)值。當(dāng)隨機(jī)變量Z(x)為分類變量(如巖相、沉積微相等)時(shí),只需將指示指標(biāo)變量的定義變?yōu)?

式中，tk為分類變量類型k，對(duì)應(yīng)于Z(x)的k種分類變量類型由閾值的個(gè)數(shù)K表示,而估計(jì)的條件概率公式變?yōu)?

2 序貫指示數(shù)值模擬的步驟[4]

首先將模擬變量劃分為K個(gè)閾值(設(shè)有n個(gè)初始值,L個(gè)網(wǎng)格節(jié)點(diǎn)),同時(shí)，定義一個(gè)通過(guò)所有網(wǎng)格節(jié)點(diǎn)的隨機(jī)路徑，然后按如下步驟進(jìn)行序貫指示模擬的實(shí)現(xiàn)。

(1)定義查找范圍內(nèi)的條件數(shù)據(jù)，條件數(shù)據(jù)包括鄰域內(nèi)的原始指示數(shù)據(jù)和模擬得到的指示值；

(2)對(duì)任意k=1,2,…,首先要由K所對(duì)應(yīng)的閾值指示協(xié)方差數(shù)值模型建立普通克里格模擬方程組,然后求解,最后由指示條件數(shù)據(jù)的線性組合計(jì)算該閾值的分布范圍函數(shù)；

(3)由Monto Carlo法得到每個(gè)節(jié)點(diǎn)上的模擬值；

(4)按上述已定義的K個(gè)閾值,將模擬值轉(zhuǎn)換為K個(gè)指示估算值；

(5)將序貫?zāi)M結(jié)果歸入條件指示數(shù)據(jù)集中；

(6)重復(fù)上述(1)-(5)項(xiàng)，沿著隨機(jī)路徑對(duì)下一個(gè)節(jié)點(diǎn)進(jìn)行序貫?zāi)M，直到所有的結(jié)點(diǎn)都被模擬。

3 應(yīng)用實(shí)例

3.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古自治區(qū)的中東部，包括吉林省和遼寧省的部分地區(qū)，屬西遼河平原。區(qū)內(nèi)的總體地形走勢(shì)為南北高、中間低；最高點(diǎn)高程為670.0m，最低點(diǎn)高程119.0m。

區(qū)內(nèi)地層主要由第四系粘性土、砂和礫石層等構(gòu)成；其地層厚度變化規(guī)律為：中心厚，邊緣薄，做為沉積中心的開(kāi)魯一帶，地層厚度可達(dá)200m。地層巖性總體變化規(guī)律為：水平方向，自山前向平原，從上游向下游顆粒由粗變細(xì)，變化過(guò)程為砂礫石—中粗砂—細(xì)砂、粉砂，且粘土夾層增多，厚度增大。

上部地層由全新統(tǒng)、上更新統(tǒng)顧?quán)l(xiāng)屯組、排頭營(yíng)子組的中細(xì)砂、粉砂、細(xì)砂和泥質(zhì)細(xì)砂組成；中部地層由中更新統(tǒng)大青溝組沖洪積、沖湖積的粗砂、中砂、細(xì)砂和淤泥質(zhì)粉砂組成；下部地層則由中更新統(tǒng)白土山組冰水、冰磧的礫石、砂礫石、泥質(zhì)砂礫石和粗砂、中細(xì)砂等組成。

3.2 三維地質(zhì)模型的建立

3.2.1 巖性概化

本次工作收集到研究區(qū)共641個(gè)鉆孔，通過(guò)分析發(fā)現(xiàn)，鉆孔年代不一，巖性描述復(fù)雜且缺乏統(tǒng)一性，考慮到序貫數(shù)值模擬要與水文地質(zhì)模型相結(jié)合的原則，根據(jù)第四系地層巖性的顆粒粒徑、透水性等水文地質(zhì)特性，將模擬區(qū)內(nèi)的第四系地層巖性概化為粘土、粉土、細(xì)砂、粗砂和礫石四類。

3.2.2 建立三維地質(zhì)體模型

根據(jù)地面高程的DEM數(shù)據(jù)，確定三維地質(zhì)體模型的頂板，同時(shí)利用鉆孔資料及地質(zhì)剖面圖確定其底板，由此建立研究區(qū)的三維地質(zhì)體模型（見(jiàn)圖1），該模型為進(jìn)行巖性隨機(jī)模擬確定出空間范圍，這樣可以從很大程度上避免頂?shù)装褰徊娴膯?wèn)題。

圖1 研究區(qū)三維地質(zhì)體模型

3.2.3 序貫指示模擬

3.2.3.1 頻率及閾值計(jì)算

將第四系巖性作為類型變量進(jìn)行指示變異函數(shù)計(jì)算。根據(jù)巖性概化，分別用模型變量1、2、3、4來(lái)代表粘土類、粉土類、細(xì)砂類、粗砂和礫石類。在模擬進(jìn)行指示變異函數(shù)計(jì)算時(shí)，分別用模擬類型變量1、2、3、4作為閾值，見(jiàn)表1。

表1 各類地層巖性頻率及其閾值

3.2.3.2 變異函數(shù)參數(shù)計(jì)算

選擇平面上的45°方向，135°方向，垂直方向的滯后距0.01，水平方向滯后距0.025、容差角20°、帶寬0.2，分別計(jì)算四個(gè)類型變量的指示變異函數(shù)參數(shù)，計(jì)算結(jié)果見(jiàn)表2。

表2 各類地層巖性變異函數(shù)參數(shù)計(jì)算結(jié)果表

3.2.3.3 序貫指示模擬

將計(jì)算所得的數(shù)值模擬參數(shù)輸入到GOCAD的序貫指示模擬模塊中進(jìn)行序貫指示模擬，即可得到西遼河平原區(qū)松散層三維地質(zhì)模型。圖2可視化呈現(xiàn)了模擬區(qū)第四系地層巖性分布，圖3以柵欄圖的方式呈現(xiàn)了研究區(qū)第四系三維地質(zhì)結(jié)構(gòu)。從圖2、圖3可見(jiàn)研究區(qū)內(nèi)地層巖性顆粒粗細(xì)分布與變化，區(qū)內(nèi)地層巖性以細(xì)砂類、粗砂類為主；向東到中部巖性粒徑由粗變細(xì)，粘粒物質(zhì)的含量逐漸增加，砂層與粘土層交錯(cuò)分布。

圖2 研究區(qū)三維地質(zhì)模型

圖3 研究區(qū)立體柵欄圖

3.2.4 模型驗(yàn)證

為了驗(yàn)證序貫指示模擬地質(zhì)模型的效果與精度，本文選擇了一條實(shí)際勘探地質(zhì)剖面圖與本次數(shù)值模擬模型結(jié)果進(jìn)行比較（見(jiàn)圖4）。從模擬地質(zhì)剖面中可以看出，研究區(qū)西部的第四系地層巖性以細(xì)砂類、粗砂類為主，向東粘粒物質(zhì)的含量逐漸增大，逐漸變成砂土層與粘土層交錯(cuò)沉積，而實(shí)際地質(zhì)剖面的西南部以細(xì)砂、中細(xì)砂及中砂為主，向北東粒度逐漸變細(xì)，巖性變?yōu)榧?xì)砂、中細(xì)砂夾多層粉土、粉質(zhì)粘土，這說(shuō)明模型模擬的結(jié)果與實(shí)際情況基本相符,所建立的三維地質(zhì)模型能夠比較真實(shí)地反映研究區(qū)的第四系地質(zhì)結(jié)構(gòu)。

圖4 序貫指示模擬剖面與實(shí)測(cè)地質(zhì)剖面對(duì)比圖

4 結(jié)論

利用GOCAD的序貫指示模擬模塊建立了西遼河平原第四系三維地質(zhì)模型，真實(shí)、直觀的反映了區(qū)內(nèi)第四系的地質(zhì)結(jié)構(gòu)，解決了以往利用二維剖面和等值線圖不能真實(shí)、直觀的反應(yīng)西遼河平原第四系地層結(jié)構(gòu)的問(wèn)題，同時(shí)為該區(qū)分析水文地質(zhì)條件、建立水文地質(zhì)概念模型和地下水流動(dòng)態(tài)評(píng)價(jià)模型提供了地質(zhì)依據(jù)。