基于Kriging方法的三維地質(zhì)可視化系統(tǒng)構(gòu)建

陳 靚， 蔣春曦， 楊 欣， 于書媛， 馬 犇， 劉世靖

（1. 安徽省地震局， 合肥 230031； 2. 合肥工業(yè)大學(xué)， 合肥 230009）

0 引言

近年來(lái)， 安徽省地震局積極推進(jìn) “一場(chǎng)一帶一站” 建設(shè)， 并將相關(guān)研究工作納入科技發(fā)展戰(zhàn)略。 各地市地震小區(qū)劃和重大工程地震安全性評(píng)價(jià)是省局長(zhǎng)期開展的重點(diǎn)項(xiàng)目， 此類項(xiàng)目都收集了豐富的資料和鉆孔數(shù)據(jù)， 利用已有數(shù)據(jù)構(gòu)建研究區(qū)域內(nèi)三維地質(zhì)體能補(bǔ)充和完善該項(xiàng)目建設(shè)成果， 不僅對(duì)地下層巖的形態(tài)和空間分布有較為直觀生動(dòng)的展示， 更可以為正在開展的合肥活斷層探測(cè)項(xiàng)目等今后 “一場(chǎng)一帶一站” 建設(shè)中開展的其他地下空間研究提供高效、 便捷的技術(shù)手段，為開發(fā)集成地形、 遙感等多源數(shù)據(jù)集成的三維地質(zhì)可視化系統(tǒng)奠定基礎(chǔ)。

三 維 地 質(zhì) 建 模（3D Geological Modeling，3DGM）是一門運(yùn)用現(xiàn)代空間信息理論來(lái)研究地層及其環(huán)境的信息處理、 數(shù)據(jù)組織、 空間建模與數(shù)字表達(dá)， 并運(yùn)用科學(xué)可視化（SciV）來(lái)對(duì)地層及其環(huán)境進(jìn)行真三維再現(xiàn)和可視化交互的科學(xué)與技術(shù)[1]。目前， 國(guó)內(nèi)外許多學(xué)者開展了與三維地質(zhì)模擬有關(guān) 的 研 究， 出 現(xiàn) 了 諸 如Micromine、 SURPAC、VULCAN、 DATAMINE 等一些具有三維地質(zhì)模擬功能的商品化地質(zhì)采礦軟件[2]， 這些系統(tǒng)一般是針對(duì)開采礦山和估算儲(chǔ)量而專業(yè)開發(fā)研制的， 對(duì)于需要通過(guò)建模提高可視化的用戶群體來(lái)說(shuō)這些軟件掌握難度大、 受眾群體小。 此外， 利用當(dāng)前研究比較多的八叉樹模型、 線框模型、 四面體格網(wǎng)模型、 實(shí)體模型、 廣義三棱柱模型和面向?qū)ο蟮捏w元拓?fù)淠Ｐ偷炔煌目臻g數(shù)據(jù)結(jié)構(gòu)來(lái)開發(fā)三維地質(zhì)體建模系統(tǒng)， 也可以實(shí)現(xiàn)3D 地質(zhì)模型可視化， 但該方法實(shí)現(xiàn)周期長(zhǎng)、 難度大、 可視化效果不如成熟的建模軟件。

鑒于上述建模方法不足， 本文選取某地震工程項(xiàng)目數(shù)據(jù)為研究對(duì)象， 利用Kriging 方法對(duì)三維地質(zhì)模型構(gòu)建進(jìn)行研究， 在ArcGIS 平臺(tái)上實(shí)現(xiàn)了三維地質(zhì)體的快速、 準(zhǔn)確建立， 并以ArcEngine 為二次開發(fā)組件， 完成地質(zhì)體三維展示系統(tǒng)的開發(fā)，實(shí)現(xiàn)三維展示、 專題圖出圖和三維分析等功能。

2 研究方法及技術(shù)方案

2.1 地統(tǒng)計(jì)插值方法

地統(tǒng)計(jì)（Geostatistics）又稱地質(zhì)統(tǒng)計(jì)， 是以區(qū)域化變量為基礎(chǔ)， 借助變異函數(shù)， 研究既具有隨機(jī)性又具有結(jié)構(gòu)性， 或具有空間相關(guān)性和依賴性的自然現(xiàn)象的一門科學(xué)。 凡是與空間數(shù)據(jù)的結(jié)構(gòu)性和隨機(jī)性， 或空間相關(guān)性和依賴性， 或空間格局與變異有關(guān)的研究， 并對(duì)這些數(shù)據(jù)進(jìn)行最優(yōu)無(wú)偏內(nèi)插估計(jì)， 或模擬這些數(shù)據(jù)的離散性、 波動(dòng)性時(shí)， 皆可應(yīng)用地統(tǒng)計(jì)學(xué)的理論與方法[3]。 克里格插值（Kriging）又稱空間局部插值法， 是以變異函數(shù)理論和結(jié)構(gòu)分析為基礎(chǔ)， 在有限區(qū)域內(nèi)對(duì)區(qū)域化變量進(jìn)行無(wú)偏最優(yōu)估計(jì)的一種方法， 是地統(tǒng)計(jì)學(xué)的主要內(nèi)容之一。 地統(tǒng)計(jì)學(xué)的理論核心為變異函數(shù)（variogram）， 具體計(jì)算公式為：

圖1 實(shí)驗(yàn)半變異函數(shù)（點(diǎn)）與擬合半變異函數(shù)模型[4]Fig.1 The model of experimental semivariogram (points) and fitted semivariogram

在地震工程項(xiàng)目應(yīng)用中， 地統(tǒng)計(jì)插值的數(shù)據(jù)源主要是已知具有三維空間坐標(biāo)和屬性的點(diǎn)， 來(lái)源于勘探點(diǎn)鉆孔數(shù)據(jù)， 建模首要任務(wù)是收集勘探點(diǎn)鉆孔數(shù)據(jù)。 由于Kriging 方法要求樣本點(diǎn)屬性值服從正態(tài)分布， 需要對(duì)收集整理后的鉆孔數(shù)據(jù)進(jìn)行探索分析， 主要包括原始數(shù)據(jù)的統(tǒng)計(jì)分析和相關(guān)的數(shù)據(jù)變換， 統(tǒng)計(jì)分析用以判斷數(shù)據(jù)是否服從正態(tài)分布。 空間數(shù)據(jù)滿足正態(tài)分布的前提下， 便可基于Kriging方法基本原理及變異函數(shù)擬合方法， 對(duì)數(shù)據(jù)進(jìn)行空間變異結(jié)構(gòu)分析， 繪制實(shí)驗(yàn)變異函數(shù)曲線圖， 并根據(jù)變異函數(shù)的理論模型擬合實(shí)驗(yàn)變異函數(shù)。 地統(tǒng)計(jì)插值的主要步驟見(jiàn)圖2 所示。

圖2 地統(tǒng)計(jì)插值的主要步驟Fig.2 The main steps of geostatistical interpolation

2.2 地層模型構(gòu)建方法

地層模型的構(gòu)建主要有不規(guī)則三角網(wǎng)（TIN）和規(guī)則格網(wǎng)法（Grid）兩種方法。 TIN 是一種由許多相鄰但又不互相重疊的三角形組成的對(duì)地形表面的連續(xù)鋪蓋。 按照空間分布的幾何結(jié)構(gòu)， 可以有一般三角網(wǎng)和Delaunay 之分。 與TIN 相對(duì)應(yīng)， Grid 是由空間中若干連續(xù)的、 形狀相同的正方形構(gòu)成的組合面， Grid 的數(shù)據(jù)來(lái)構(gòu)造連續(xù)、 平順、 光滑的曲面是可能的， 且數(shù)學(xué)模型的構(gòu)造也較簡(jiǎn)單[1]。 本文采用的地層模型構(gòu)建方法是TIN 的生成與拉伸。

2.3 技術(shù)方案

本文所采用的數(shù)據(jù)包含多種格式文件， 在構(gòu)建三維模型之前需要先進(jìn)行數(shù)據(jù)格式的轉(zhuǎn)換和幾何校正。 對(duì)鉆孔數(shù)據(jù)進(jìn)行讀圖整理、 屬性表設(shè)計(jì)、巖層數(shù)據(jù)的檢查與補(bǔ)充， 構(gòu)建鉆孔空間數(shù)據(jù)庫(kù)。之后， 利用地統(tǒng)計(jì)學(xué)的Kriging 插值、 Tin 的生成與拉伸等方法實(shí)現(xiàn)三維地質(zhì)體的構(gòu)建。 并在此基礎(chǔ)上， 進(jìn)行巖層的三維可視化表達(dá)。 最后， 本文以ArcEngine 為二次開發(fā)組件， 完成地質(zhì)體三維展示系統(tǒng)的開發(fā)， 總體技術(shù)流程見(jiàn)圖3 所示。

圖3 技術(shù)流程Fig.3 Technique flowchart

3 三維地質(zhì)體的構(gòu)建

3.1 數(shù)據(jù)處理與地統(tǒng)計(jì)探索

3.1.1 數(shù)據(jù)預(yù)處理

數(shù)據(jù)預(yù)處理階段主要工作是將研究區(qū)域內(nèi)的勘探點(diǎn)數(shù)據(jù)進(jìn)行坐標(biāo)整理、 通過(guò)讀取鉆孔柱狀圖來(lái)實(shí)現(xiàn)數(shù)據(jù)的分層及屬性信息的完善， 共完成了研究區(qū)98個(gè)鉆孔點(diǎn)的整理和入庫(kù)。 數(shù)據(jù)屬性表結(jié)構(gòu)見(jiàn)表1 所示。

3.1.2 定義投影

在對(duì)本身有投影信息的數(shù)據(jù)采集完成時(shí)， 為了保證數(shù)據(jù)的完整性和易交換性， 要對(duì)數(shù)據(jù)定義投影， 針對(duì)鉆孔的分層數(shù)據(jù)在ArcToolbox—Data Management Tools—Projections and Transformations中將已獲取的數(shù)據(jù)統(tǒng)一定義坐標(biāo)系統(tǒng)， 并進(jìn)行坐標(biāo)系統(tǒng)的統(tǒng)一轉(zhuǎn)換。

表1 勘探點(diǎn)數(shù)據(jù)屬性表結(jié)構(gòu)Table 1 Attribute table structure of borehole data

3.1.3 地統(tǒng)計(jì)探索

地統(tǒng)計(jì)探索是利用一系列的圖形工具和適用于數(shù)據(jù)的插值方法， 確定統(tǒng)計(jì)數(shù)據(jù)屬性、 探測(cè)數(shù)據(jù)分布、 全局和局部異常值（過(guò)大值或過(guò)小值）、尋求全局的變化趨勢(shì)、 研究空間自相關(guān)和理解多種數(shù)據(jù)集之間相關(guān)性。 利用直方圖、 QQPlot 分布圖等工具來(lái)進(jìn)行數(shù)據(jù)的空間探索， 如檢驗(yàn)出原始數(shù)據(jù)沒(méi)有顯示出正態(tài)分布， 則可根據(jù)需要選擇適當(dāng)?shù)淖儞Q方式將數(shù)據(jù)進(jìn)行轉(zhuǎn)換， 使之服從正態(tài)分布。

QQPlot 圖提供了一種度量數(shù)據(jù)正態(tài)分布的方法， 利用QQPlot 圖， 可以將現(xiàn)有數(shù)據(jù)的分布與標(biāo)準(zhǔn)正態(tài)分布進(jìn)行對(duì)比， 如果數(shù)據(jù)越接近一條直線，則它越接近于服從正態(tài)分布。 本文對(duì)地表層勘探點(diǎn)數(shù)據(jù)進(jìn)行地統(tǒng)計(jì)探索時(shí)， 發(fā)現(xiàn)其不太符合正態(tài)分布， 經(jīng)過(guò)log 變換后， 數(shù)據(jù)近似符合正態(tài)分布，見(jiàn)圖4 所示。 同理， 對(duì)其他地層勘探點(diǎn)數(shù)據(jù)也可以用上述方法進(jìn)行地統(tǒng)計(jì)探索。

圖4 經(jīng)log 變換后的地表層正態(tài)QQPlot 圖Fig.4 Normal QQPlot graph of surface layer after log transformation

3.2 三維模型的構(gòu)建

3.2.1 等高線生成地統(tǒng)計(jì)探索之后， 對(duì)地表層的鉆孔數(shù)據(jù)選擇合適的變換方法， 進(jìn)行Kriging 插值， 得到等高線， 并生成預(yù)測(cè)圖， 見(jiàn)圖5 所示。 重復(fù)kriging 插值步驟， 可以獲得各地層等高線預(yù)測(cè)圖。

3.2.2 Tin 的生成及拉伸

在ArcScene 下的3D Analyst 的“創(chuàng)建TIN” 工具， 加載巖層等高線作為輸入要素， 選擇研究區(qū)范圍， 應(yīng)用3D 分析模塊進(jìn)行兩個(gè)巖層間的模型拉伸， 得到了一層的地質(zhì)體， 對(duì)區(qū)域內(nèi)所有地層進(jìn)行上述重復(fù)操作， 即得到研究區(qū)范圍內(nèi)的三維模型， 如圖6 所示。 重復(fù)以上步驟， 即可得到該工程場(chǎng)地所有地層三維模型。

圖5 地表層等高線預(yù)測(cè)圖Fig.5 Prediction map of ground surface contour

圖6 粉質(zhì)粘土層地質(zhì)體Fig.6 Geologic body of silty clay layer

3.2.3 地質(zhì)體的三維可視化

數(shù)字高程模型（DEM）與遙感影像是建立地表形態(tài)的數(shù)據(jù)來(lái)源， 遙感影像所含的數(shù)據(jù)信息量遠(yuǎn)高于普通地形圖所含信息量， 再加上數(shù)字高程模型表現(xiàn)出的地形起伏特征， 三維可視化的效果要好于傳統(tǒng)方法對(duì)地形的表現(xiàn)[5]。 常用的DEM 有兩種形式： 一種稱為高程矩陣， 另一種DEM 的高程數(shù)據(jù)布滿覆蓋整個(gè)區(qū)域的三角網(wǎng)網(wǎng)點(diǎn)， 實(shí)際上就是TIN 數(shù)據(jù)[6]。 將遙感數(shù)據(jù)導(dǎo)入ArcScene 中， 將TIN 設(shè)置為遙感影像的Base Height， 實(shí)現(xiàn)地形的三維可視化（圖7）。

圖7 含地表紋理三維地質(zhì)體Fig.7 The 3D geological body with surface texture

3.2.4 模型的裁切

如果想查看地質(zhì)體內(nèi)部結(jié)構(gòu)， 可以在已實(shí)現(xiàn)的三維地質(zhì)體中， 任意做剖面， 采用相交運(yùn)算后，得到裁切后的剖面體， 見(jiàn)圖8 所示。

4 三維地質(zhì)可視化系統(tǒng)的開發(fā)

圖8 剖面體Fig.8 Section body

在實(shí)現(xiàn)地質(zhì)體的三維構(gòu)建基礎(chǔ)上， 采用C# 作為開發(fā)語(yǔ)言， 以Visual Studio 2010 為開發(fā)環(huán)境，ArcEngine 為二次開發(fā)組件， 設(shè)計(jì)開發(fā)了三維可視化系統(tǒng)， 實(shí)現(xiàn)三維展示、 專題圖出圖和三維分析等功能（圖9）。

5 結(jié)論

本文對(duì)利用ArcGIS 軟件進(jìn)行地質(zhì)體三維建模方法進(jìn)行了研究， 在建立巖層三維模型時(shí)， 充分采用了Kriging 插值地統(tǒng)計(jì)方法對(duì)數(shù)據(jù)進(jìn)行了正態(tài)分布探索， 使插值后的三維曲面呈現(xiàn)出較為光滑真實(shí)的效果。 基于ArcEngine 開發(fā)組件進(jìn)行三維地質(zhì)可視化系統(tǒng)的開發(fā)， 具有較高的開發(fā)環(huán)境優(yōu)勢(shì)，界面友好， 執(zhí)行效率較穩(wěn)定。 本文的研究為地震工程場(chǎng)地地形三維可視化和分析提供了科學(xué)手段，有良好的可視化效果， 具有一定的實(shí)用價(jià)值。 本系統(tǒng)相關(guān)技術(shù)升級(jí)后， 其數(shù)據(jù)管理和可視化分析功能， 可以為地震地質(zhì)調(diào)查、 地震安評(píng)等提供有力的技術(shù)支持。 由于數(shù)據(jù)源等客觀條件限制， 目前此方法在對(duì)大范圍地質(zhì)體的構(gòu)建和可視化程度上效率不高， 三維開發(fā)分析功能還有待進(jìn)一步探索和研究。

圖9 系統(tǒng)主界面Fig.9 Main interface of system

[1] 張凱選， 武文波， 白云峰， 等. 基于ArcGIS 的地質(zhì)三維可視化[J]. 遼寧工程技術(shù)大學(xué)學(xué)報(bào)， 2007， 26（3）：345-347.

[2] 趙州. 基于剖面的三維地質(zhì)建模研究[D]． 西安： 西安科技大學(xué)， 2004.

[3] 湯國(guó)安. ARCGIS 地理信息系統(tǒng)空間分析實(shí)驗(yàn)教程[M].北京： 科學(xué)出版社， 2006.

[4] 喬金海， 潘 懋， 金 毅， 等.基于Kriging 方法的天然地基承載力三維模擬及分析[J]. 北京大學(xué)學(xué)報(bào)（自然科學(xué)版）， 2011， 47（5）： 812-818.

[5] 陳靚， 路海濤， 謝慶勝， 等. 基于3D GIS 和RS 的安徽省地震監(jiān)測(cè)臺(tái)網(wǎng)電子沙盤制作的實(shí)現(xiàn)[J]. 西北地震學(xué)報(bào)， 2010， 32（4）： 388-392.

[6] 劉南， 劉仁義. 地理信息系統(tǒng)[M]. 北京： 高等教育出版社， 2002.

[7] 肖海紅. 基于ArcScene 的三維地形可視化及其應(yīng)用[J].工程地質(zhì)計(jì)算機(jī)應(yīng)用， 2007（3）： 12-16.

[8] 宋仁波， 江南. 基于GIS 的三維地質(zhì)可視化系統(tǒng)設(shè)計(jì)與實(shí)現(xiàn)[J]. 工程地質(zhì)計(jì)算機(jī)應(yīng)用， 2011， 64（4）： 22-27.

[9] 張敏， 李自紅， 劉鴻福， 等. 基于ArcGIS 的太原斷陷盆地第四系三維地質(zhì)建模[J]. 震災(zāi)防御技術(shù)， 2007， 2（3）： 243-248.

[10] 葛偉鵬， 袁道陽(yáng)， 郭華. 對(duì)城市活斷層探測(cè)項(xiàng)目中地震地質(zhì)數(shù)據(jù)建模的探討[J].西北地震學(xué)報(bào)， 2006， 28（2）： 134-139.

[11] 李自紅， 劉鴻福， 張敏， 等. 地震與活動(dòng)斷裂空間關(guān)系的三維可視化建模[J]. 地震地質(zhì)， 2013， 35（3）：565-575.

[12] 劉玉娟. 地震活斷層及其所在地質(zhì)空間的三維建模與可視化[D]． 北京： 中國(guó)地震局地殼應(yīng)力研究所， 2007.

[13] 黃忠. 地質(zhì)三維建模與可視化技術(shù)的研究與應(yīng)用[D]． 合肥： 合肥工業(yè)大學(xué)， 2008.

[14] 戴晟. 地質(zhì)體三維可視化研究與系統(tǒng)實(shí)現(xiàn)[D]． 上海： 華東師范大學(xué)， 2008.

[15] 程慧. 基于地震資料的地質(zhì)建模及效果評(píng)價(jià)[D]． 青島：中國(guó)石油大學(xué)（華東）， 2009.

[16] 呂亞楠. 地質(zhì)體三維可視化應(yīng)用研究[D]． 西安： 長(zhǎng)安大學(xué)， 2010.

[17] 張洋洋， 周萬(wàn)蓬， 吳志春， 等. 三維地質(zhì)建模技術(shù)發(fā)展現(xiàn)狀及建模實(shí)例[J].華東理工大學(xué)學(xué)報(bào)（社會(huì)科學(xué)版）， 2013， 32（3）： 403-408.