基于SKUA的淺層三維地質建模方法與實例研究

杜 睿，楊歧焱

DU Rui, YANG Qi-yan

河北地質大學 京津冀城市群地下空間智能探測與裝備自然資源部重點實驗室 河北省戰略性關鍵礦產資源重點實驗室，河北 石家莊 050031

Hebei GEO University, Shijiazhuang 050031, China

0 引言

傳統的勘察技術及數據處理手段單一，得到的二維矢量圖（平面圖、剖面圖）不能直觀、準確、清晰地表達地下地質情況。三維地質建模是在三維環境中，利用計算機軟硬件系統，對空間中的地質數據進行采集、處理、解釋并與地學統計、空間分析、地質地層預測、三維可視化等相關工具結合起來，并用于地質方面分析的一門技術[1]。20世紀80年代，已有學者開始對三維地質建模進行研究，Haldorson提出了隨機模擬儲層建模[2]方法和Calson E地下結構有關的三維概念模型[3]，A.B.Ekoule解決了三維物體非凸輪廓線的重建問題[4]；Bak、Smith等提出了一些三維模型和地學信息的三維表示[5]。直到90年代初“三維地質建模”的概念才被加拿大學者Simon W.Houlding正式提出[6]；之后，針對地質數據本身具有的不確定性和不連續性的特點，Mallet提出了“離散光滑插值”算法[7,8]，使得構造建模技術的地質曲面技術得到了突破，美國T-surf公司設計基于離散光滑插值算法研發了GOCAD（Geological Object Computer Aided Design）三維地質建模軟件；隨后，Carl Youngman（1990），Molenaar Marien（1991），Fritsch D，AaperJ.F（1992），Molenaar M，TurnerA.K（1992），Thomas R.Fisher（1993），Rongxingli（1994），Lattuada R（1995），DekempE.A（1999；2003），HerbertM J（2001），Saini-Eidukat B（2002）等人又進行了大量的研究，主要包括三維數據的可視化、數據的空間模型與空間結構、矢量化地圖的三維數據結構等方面[9-11]；2006年GOCAD被加拿大Paradigm公司收購；2007年Paradigm公司推出SKUA-GOCAD，也稱為SKUA（Subsurface Knowledge Unified Approach，SKUA）軟件；2011年Paradigm公司進行二次開發，將SKUA集成到Epos數據庫平臺上（Paradigm公司內部資料），推出了可以建立任何復雜地質模型（地層、尖滅、斷裂等）的SKUA軟件，為三維地質建模理論的發展做出了貢獻。

1 空間插值擬合

由于工程鉆孔價格昂貴，能獲得的有用鉆孔較少。因此通常需要對原始采樣數據進行插值擬合。所謂插值就是指滿足某種變異函數的前提下，尋找有限區域已知離散數據附近的擬合點，從而在整個研究區范圍內生成足夠多的連續、分布均勻的數據點[12]。由于不同的插值方法有不同的適應條件，本文考慮到樣本點和待估點之間的位置關系，還考慮到待估點之間具有的空間相關性，擬采用普通克里金插值和離散光滑插值。

1.1 克里金插值

克里金（Kriging）計算權值的方式是通過引入一個以距離為自變量的半變差函數來實現，確定權值后即可確定待測值點的位置，從而確定待測點處的最佳線性估計值[13]。克里金插值考慮了所有觀測點的空間分布結構特點，與傳統插值方法相比，其插值結果更接近實際情況[14,15]。

其算法如下：設xi（i=1，2，…，n）為空間數據點的坐標，zi（i=1，2，…，n）表示其對應的觀測值。對于某點x0處的值記為Z（x0），則Z*（x0）可用克里金插值公式

權重λi由方差最小和無偏性估計來唯一確定一組λi值：

本文利用確定出的權值進行克里金插值計算，多次試驗后發現普通克里金插值能夠很好地擬合研究區地質信息的變化規律，并對準確估計未知點變量。

1.2 離散光滑插值

離散光滑插值（Discrete Smooth Interpolation，簡稱DSI）是Mallet教授針對CAD平面制圖軟件在地質構造數據內插上的不足提出的，其基本內容是：首先建立一個相互聯系的自然體網格，然后利用網格結點的鄰接關系，完成從已知結點函數值到未知結點函數值的估算[16]。離散光滑插值作為SKUA的特色插值技術，在SKUA軟件的構造建模和速度建模中得到了成熟的應用。

現將DSI算法的理論基礎簡述如下：

假設G=G（S）為曲線多邊形，Ω=Ω（S）為曲線多邊形上的頂點集合，N=|Ω|為在Ω上的點，φ（k）為定義在k∈Ω上的函數，只在Ω的子集上φ（k）為已知的，φ為一個n維列矩陣φ=[φ（1），（2），…，φ（n）]T。

為了方便討論，將φ矩陣分為矩陣φI和φL（其中φI為已知矩陣，φL為未知矩陣）。

首先，要找出R*(?) =R(?) +ρ(?) = min的φ（k）函數（其中，R（φ）為全局粗糙度，ρ（φ）為線性約束），并利用φ（k）函數對曲面進行估計。

a.粗糙度

R（φ|k）是局部粗糙度，R（φ）為全局粗糙度

式中，{να(k)}為給定的權系數，其值為整數；

N（k）為節點k的鄰域；

μ（k）在Ω上給定的非負權函數。

從上述式中可以得到R（φ）是根據{να(k)}和μ（k）確定的。

b.最小平方的線性約束ρ（φ）

對于n維矩陣Ai和一個常數bi，定義有：

c.問題解決

由φ（k）函數（k∈Ω）插值所組成的值域中確定出使R*(?) =R(?) +ρ(?)值最小的一個數，然后將R*(?)展開有：

式中，W*為對稱半正定矩陣，Q為列矩陣，β2為非負系數。

其中，

依據上述分析把φ矩陣劃分為φI和φL兩個矩陣，能夠得到矩陣W*和矩陣Q兩個矩陣。

綜合分析上面的數學求解過程，可以知道離散光滑插值（DSI）的基本思想是找出一個φ（k）函數使R*（φ）達到最小值，其值完全依賴于權函數{να(k)}和μ（k）的取值，其中II方陣秩的唯一性決定了DSI插值法解的唯一性。

2 建模思路

三維地質建模主要分為兩類：一類是構造和地層（Structure & Stratigraphy modeling）建模，適用于鉆孔數據較為豐富的區域，構造和地層模型能夠表達地層面和斷層面的位置和相互關系，能夠直接觀測到地表的跌宕起伏；另一類是三維儲層柵格結構[17]（3D Reservoir Grid Building）建模，適用于高光譜遙感圖像、數字高程圖像等方法的區域性建模。儲層建模能反映地質體的物性參數，如巖石的滲透率、孔隙度等。其中構造建模是儲層建模的基礎，本文采用SKUA軟件和構造建模為主闡述快速三維地質建模基本流程。

2.1 建模原始資料

本文在滄州市高新區某工程區10.2 km2的范圍內，根據工程的實際需要，收集的數據情況如下（表1）。理論上講，鉆孔的間距越小數量越多、工程地質剖面圖越詳細、對區域地質的認識越深入，越有利于地質建模工作的開展[18]。

表1 滄州市某工區建模數據一覽表Table 1 List of modeling data for a work area in Cangzhou city

2.2 鉆孔數據庫的建立與導入

本文以60張鉆孔柱狀圖，10張工程地質剖面圖為原始數據，依據地震反射同相軸分層的原則將地層分為4層。首先將圖、表等數據導入EXCEL整理后導入ArcGIS，進行坐標轉換，轉換成功后導出坐標系并保存為SKUA-GOCAD能夠直接導入的TXT格式。SKUA-GOCAD軟件對鉆孔數據的錄入格式如下：井軌跡包含鉆孔編號（WellName）、測深（MD）及坐標（X，Y，Z）等5列數據；鉆孔測斜表包括鉆孔編號（WellName）、傾角（Dip）、方位角（Azimuth）、測深（MD）、分層名稱（Markername）等5列數據[19]。數據導入之前最好對數據進行全面排查，尤其是小數位數，正負號等，避免建模過程產生各種意想不到的誤差。

鉆孔數據導入SKUA過程如下：（1）導入井軌跡數據：Import→Well Data→Paths→column-Based File→選擇Wellpath數據→Next→選中X-Y-TVDSSMD→選中Meter→選中Use a column→Next→為每列匹配相應的參數；（2）導入鉆孔測斜表：Import→Well Data→Markers→column-Based File→選擇Markers數據→Next→選中By column→Next→選擇有無MD和Dip and Azimuth→Next→選中Overwrite it→Next→為每列匹配對應的參數。

2.3 建立點模型

原始點模型的建立可以直接用井分層文件來建立。也可以用各個地層面的點數據來建立，過程如下：PointsSet→New→From Well Markers→Locations→選中所有的井文件→選擇相應的Marker’s feature→Name。由于鉆孔信息提供的原始數據點較少，為了快速準確的進行建模需要進行空間插值，本文分別對每一層的原始數據點進行克里金插值，以A3層為例，插值前后效果如下圖1所示

圖1 （a） 原始數據點； （b）克里金插值后數據點Fig.1 （a）Original data points；（b）Data points after kriging interpolation

2.4 建立面模型

面模型是由點模型生成的，面模型不僅可以表現層面的高低起伏，還可以展現地質體結構和地質體形態，是構造建模中的核心部分[20]。過程如下：Surface→New→From PointsSets→PointsSet→選擇相應的點文件，SKUA會自動給Surface命名與點文件同名的名字，點擊apply生成面模型（圖2（a））。

從圖3中的面模型可以看出，部分插值點與地層面存在不一致的情況，盡管克里金插值優于其他算法，但它們都有一個共同的缺點，插值數據生成的曲面與原始控制點并不完全一致[21]。為此，SKUA軟件在DSI理論的基礎上提供了一種幾何適應（fit geometry）的功能，它可以將插值曲面與不在曲面上的控制點進行重新匹配，從而不斷修正幾何畸變，得到更加符合實際的曲面。操作過程如下：Surface→Constraints→Control Points→Set Control Points→選擇相應的地層面和控制點→點擊Fit Geometry，多次重復，直到控制點完全落在層面上（圖2（b））。

圖2 （a）由插值數據生成的面模型； （b）幾何適應后的面模型Fig.2 （a）Surface model generatedby the interpolation data；（b）Surface model after fit geometry.

2.5 建立構造模型

點擊Structure & Stratigraphy→選擇深度域→選中所有的地層→然后建立地層柱→選中所有的well→Check Data Consistency→Next→定義感興趣的區域（感興趣的應大于最頂層和最底層的高度）→Next→選中所有的地層→Preview。建立地層之前先預覽各層面，預覽地層面時可以進行如下操作：1）逐層建立層面，忽略不一致的數據點；2）檢測是否存在交叉層面；3）根據剖面圖對層面局部區域進行微調。預覽后再Build All建立各地層→Compute Paleo-coordinates→Next→Building the Geologic Grid→Compute Unit Areal Extension→Compute Unit Thickness，最小單元不同和厚度不同均會影響計算效率，本文根據推薦采用的Cell size為39.822 m，Set Number of Cells為10。Build Geologic Grid建立地質格網如圖3。

圖3 研究區三維地質構造模型Fig.3 3D geological structure model of the study area

2.6 地質剖面與解釋剖面做對比

地質剖面是地質圖件的一種。它是沿地球表面某一方向的假想垂直面與地形相切得到的剖面圖，稱為地質剖面圖[22]。它能直觀地反映出剖面方向、地形及地層的厚度、巖性、斷層和傾角等信息。其操作如下：General→Cross Section→選擇構造模型→New→From Digitized Polyline→從左往右任意畫一條線，生成的剖面如圖4（a）。

圖4 （a）地質構造模型剖面圖； （b）某測線地震反射時間剖面Fig.4 （a）Geological structure model profile view ；（b）Seismic reflection timeprofile of a seismic line

隨機從該工程區三條地震測線中取出一條，其測線剖面如圖4（b），從地震勘探地質解釋剖面上可以看出，該測線有多組較強的反射層位，比較明顯的反射層在90 m、150 m、260～300 m、400～420 m的深度，推測其分別為T1、T2、T3、T4。將圖4（a）與圖4（b）對比可知，構造模型剖面與任意一條測線的剖面大致相吻合，構造模型較好。

3 總結

本文提出了在地層情況較為簡單，但鉆孔數據不足的情況下，結合空間插值擬合技術進行快速建模的方法。針對滄州某工區插值后的控制點進行幾何適應，去除離群值，幾乎都能落在層面上，插值效果較好，建模方法簡便快捷，可以廣泛應用于工程地質、油氣開發、地質災害治理、地球物理勘探等領域，與二維平面圖相比能提供一種更直觀的表現形式。但是，三維地質建模相關理論和軟件在國內地質領域應用還較少，有不少復雜問題還需進一步研究。