序貫指示模擬法在三維地質模型中的應用
——以西遼河平原為例

楊明武 龍文華

1.中煤科工集團北京華宇工程有限公司 北京 100120

2.廣東省地質環境監測總站 廣東 廣州 510510

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

GOCAD軟件作為新一代地質建模軟件的代表,具有通用地質模型所設想的功能——全三維、全拓撲、使用空間插值算法、擁有完善的地質統計分析功能等特征，達到了半智能化地質建模的世界最高水平,在國外受到眾多用戶青睞[1]。

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

1 序貫指示模擬原理

序貫模擬是地質統計學的一個重要組成部分，是將序貫思路與克立格插值相結合的條件模擬方法，將沿著隨機路徑序貫求出各地質網格結點的條件累積分布函數值作為主要開發思路，并從中取得數值模擬值。常用的有序貫高斯模擬和序貫指示模擬。序貫高斯模擬主要適用于滿足高斯分布的連續型數據場，而序貫指示模擬則主要適用于離散型數據場，也可用于離散化的連續型數據場，其基本原理為[2,3]：

利用線性無偏估計克里格法,來估計條件概率的公式為：

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

式中，tk為分類變量類型k，對應于Z(x)的k種分類變量類型由閾值的個數K表示,而估計的條件概率公式變為:

2 序貫指示數值模擬的步驟[4]

首先將模擬變量劃分為K個閾值(設有n個初始值,L個網格節點),同時，定義一個通過所有網格節點的隨機路徑，然后按如下步驟進行序貫指示模擬的實現。

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

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

(3)由Monto Carlo法得到每個節點上的模擬值；

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

(5)將序貫模擬結果歸入條件指示數據集中；

(6)重復上述(1)-(5)項，沿著隨機路徑對下一個節點進行序貫模擬，直到所有的結點都被模擬。

3 應用實例

3.1 研究區概況

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

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

上部地層由全新統、上更新統顧鄉屯組、排頭營子組的中細砂、粉砂、細砂和泥質細砂組成；中部地層由中更新統大青溝組沖洪積、沖湖積的粗砂、中砂、細砂和淤泥質粉砂組成；下部地層則由中更新統白土山組冰水、冰磧的礫石、砂礫石、泥質砂礫石和粗砂、中細砂等組成。

3.2 三維地質模型的建立

3.2.1 巖性概化

本次工作收集到研究區共641個鉆孔，通過分析發現，鉆孔年代不一，巖性描述復雜且缺乏統一性，考慮到序貫數值模擬要與水文地質模型相結合的原則，根據第四系地層巖性的顆粒粒徑、透水性等水文地質特性，將模擬區內的第四系地層巖性概化為粘土、粉土、細砂、粗砂和礫石四類。

3.2.2 建立三維地質體模型

根據地面高程的DEM數據，確定三維地質體模型的頂板，同時利用鉆孔資料及地質剖面圖確定其底板，由此建立研究區的三維地質體模型（見圖1），該模型為進行巖性隨機模擬確定出空間范圍，這樣可以從很大程度上避免頂底板交叉的問題。

圖1 研究區三維地質體模型

3.2.3 序貫指示模擬

3.2.3.1 頻率及閾值計算

將第四系巖性作為類型變量進行指示變異函數計算。根據巖性概化，分別用模型變量1、2、3、4來代表粘土類、粉土類、細砂類、粗砂和礫石類。在模擬進行指示變異函數計算時，分別用模擬類型變量1、2、3、4作為閾值，見表1。

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

3.2.3.2 變異函數參數計算

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

表2 各類地層巖性變異函數參數計算結果表

3.2.3.3 序貫指示模擬

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

圖2 研究區三維地質模型

圖3 研究區立體柵欄圖

3.2.4 模型驗證

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

圖4 序貫指示模擬剖面與實測地質剖面對比圖

4 結論

利用GOCAD的序貫指示模擬模塊建立了西遼河平原第四系三維地質模型，真實、直觀的反映了區內第四系的地質結構，解決了以往利用二維剖面和等值線圖不能真實、直觀的反應西遼河平原第四系地層結構的問題，同時為該區分析水文地質條件、建立水文地質概念模型和地下水流動態評價模型提供了地質依據。