基于多點統(tǒng)計學(xué)的地質(zhì)體三維重構(gòu)及其在地鐵工程中的應(yīng)用*

劉成軍，劉恒光，侯衛(wèi)生，3，4，王典，陳曉丹，張華，陳勇華

1. 廣州地鐵設(shè)計研究院股份有限公司，廣東 廣州 510010

2. 中山大學(xué)地球科學(xué)與工程學(xué)院，廣東 珠海 519082

3. 南方海洋科學(xué)與工程廣東省實驗室（珠海），廣東 珠海 519082

4. 廣東省地球動力學(xué)作用與地質(zhì)災(zāi)害重點實驗室，廣東 珠海 519082

經(jīng)濟發(fā)展促進了城市尤其是大城市、城市群的快速擴張。大力發(fā)展地下軌道交通是解決城市高速發(fā)展中所面臨的交通擁擠、實現(xiàn)土地集約利用等問題的一個良好途徑。然而，漫長地質(zhì)過程不僅形成了形態(tài)復(fù)雜多變的地質(zhì)體，也改變了巖石的結(jié)構(gòu)和力學(xué)性質(zhì)。建立合理的、高精度、大比例尺的三維地質(zhì)結(jié)構(gòu)模型，明確地下軌道交通沿線的地質(zhì)條件，是地鐵工程的設(shè)計和施工中一個重要的數(shù)據(jù)框架和基礎(chǔ)。

三維地質(zhì)模型是基于現(xiàn)代空間信息理論和技術(shù)建立的一種具有地質(zhì)意義的數(shù)學(xué)表達，給出了易于理解的地質(zhì)體空間結(jié)構(gòu)三維透視圖，為基礎(chǔ)地質(zhì)調(diào)查、工程設(shè)計、施工和決策等提供基礎(chǔ)數(shù)據(jù)和決策支持［1-6］。假設(shè)數(shù)據(jù)分布服從二階平穩(wěn)條件，基于變差圖的空間分布函數(shù)，兩點地質(zhì)統(tǒng)計學(xué)方法可充分利用軟數(shù)據(jù)和硬數(shù)據(jù)約束，實現(xiàn)復(fù)雜地質(zhì)空間結(jié)構(gòu)的三維重構(gòu)［7-9］。然而，基于變差圖的空間關(guān)系函數(shù)缺乏對數(shù)據(jù)整體空間相關(guān)性的度量，使得兩點地質(zhì)統(tǒng)計學(xué)方法在模擬非線性地質(zhì)結(jié)構(gòu)上存在顯著偏差，難以準確揭示地質(zhì)結(jié)構(gòu)的空間變化［7，10-13］。

與兩點統(tǒng)計學(xué)相對應(yīng)的是，多點統(tǒng)計學(xué)（MPS，multiple-point statistics）方法從訓(xùn)練圖像（TI，training image）中提取多點間不同尺度模式結(jié)構(gòu)，通過隨機模擬過程重組這些空間特征。MPS 著重關(guān)注空間中多點之間相關(guān)性的表達，可在模擬過程中引入地質(zhì)知識，使模擬結(jié)果不僅忠于建模數(shù)據(jù)，還能有效描述地質(zhì)體的復(fù)雜空間幾何形態(tài)，更好地服務(wù)于各種科學(xué)研究和應(yīng)用領(lǐng)域，如油氣藏勘探［14-17］、煤田地質(zhì)［18］、水文地質(zhì)［19-21］、遙感制圖［22，23］、礦產(chǎn)資源評價［24］、微觀結(jié)構(gòu)特征［25-28］等。MPS 模擬方法有兩個關(guān)鍵步驟：一是當前數(shù)據(jù)事件搜索，即從TI 中搜索出符合當前數(shù)據(jù)事件條件的空間模式；二是數(shù)據(jù)事件更新，即對所搜索出的模式進行更新和粘貼。從粘貼像素的數(shù)目來看，已有MPS 方法大致可以分成基于點模式和基于模板的方法兩大類。前者在每次模擬過程中僅粘貼一個像素點，常用方法主要有SNESIM［29］、DS［7］、HOSIM［30-31］和IMPALA［32］等?；谀０宓姆椒ㄔ诿看文M中粘貼一個模板大小的數(shù)據(jù)，如GOSIM［33］、SIMPAT［34］、FILTERSIM［35］、DISPAT［36］和MS-CCSIM［37］等。除GOSIM 外，其余方法均以序貫?zāi)M為基礎(chǔ)。這種基于序貫?zāi)M的方法在數(shù)據(jù)事件更新過程中所產(chǎn)生的誤差難以得到修正，并在模擬過程中不斷累積以致模擬結(jié)果偏差較大。引入具有更低敏感度的柵格路徑方法可以減小誤差，但誤差累積仍蘊含其內(nèi)。結(jié)合圖像重構(gòu)中的多尺度迭代計算方法，GOSIM 克服了傳統(tǒng)迭代MPS 方法難以收斂的問題，并在模擬過程中不斷修正誤差，獲得了更符合實際的模擬結(jié)果［33］。

本質(zhì)上，三維地質(zhì)建模的核心在于尋找并明確已知地質(zhì)數(shù)據(jù)之間的空間關(guān)系，推測未采樣位置的地質(zhì)屬性值。與區(qū)域地質(zhì)調(diào)查不同的是，受城市建設(shè)的影響，地鐵工程地質(zhì)調(diào)查的橫向?qū)挾认鄬^窄，鉆孔數(shù)據(jù)卻相對密集，形成了具有數(shù)據(jù)量少、但局部密度高的數(shù)據(jù)特點。常規(guī)的三維地質(zhì)建模方法難以從已知數(shù)據(jù)中提取出地質(zhì)屬性的空間分布結(jié)構(gòu)。因此，數(shù)據(jù)量少、局部密度高數(shù)據(jù)條件下的地質(zhì)結(jié)構(gòu)三維重建成為三維地質(zhì)模型在地鐵工程中深入應(yīng)用的一個瓶頸。本文以Extended GOSIM 算法［38］為核心，結(jié)合地鐵工程數(shù)據(jù)特點，以二維地質(zhì)剖面作為建模的基礎(chǔ)數(shù)據(jù)集，從分級融合的角度探討面向地鐵工程應(yīng)用的地質(zhì)結(jié)構(gòu)三維重建方法。

1 基于MPS 的地質(zhì)結(jié)構(gòu)三維重建方法

1.1 MPS的基本概念

MPS 使用數(shù)據(jù)模板掃描TI 獲取地質(zhì)數(shù)據(jù)的空間模式以及地質(zhì)變量的統(tǒng)計量，并利用這些統(tǒng)計量隨機模擬出與TI 相似的特征。TI 是對現(xiàn)實世界的一個隨機過程的抽象描述，是描述地下介質(zhì)結(jié)構(gòu)的二維或三維圖像?？梢姡琓I 的來源也可以多種多樣，如露頭照片的解釋，地質(zhì)剖面圖，或者是基于對象或基于過程的模擬實現(xiàn)。本研究采用了二維鉆孔地質(zhì)剖面作為TI進行三維地質(zhì)模擬。

數(shù)據(jù)模板T是由空間中n個向量組成的一種幾何結(jié)構(gòu)。假定模板中心位置為u，那么模板就可以表示為：T=｛u，u1，u2，…，un-1｝=｛u+h0，u+h1，u+h2，…，u+hm｝（如圖1所示），其中h0，h1，h2，…，hm為滯后距離。

圖1 一個二維的5×5的數(shù)據(jù)模板示意圖Fig.1 Schematic diagram of two dimensional data template with size of 5×5

MPS 在獲取TI 的空間模式時，從TI 的某個位置開始，以數(shù)據(jù)模板大小為窗口，獲取窗口內(nèi)的數(shù)據(jù)屬性分布狀態(tài)；然后，滑動窗口，再獲取新窗口范圍內(nèi)的屬性分布，直至TI內(nèi)所有的結(jié)構(gòu)均有效獲取，其過程如圖2（a）所示。掃描結(jié)束后，MPS即可獲得如圖2（b）所示的空間模式集。通常而言，對TI 進行搜索的時間相對較長。為了加速搜索模式，本研究在獲取空間模式集后，構(gòu)建了模式數(shù)據(jù)庫，避免了模擬過程中反復(fù)對TI進行全局搜索。

圖2 MPS的空間模式獲取過程示意圖Fig.2 Schematic diagram of extracting spatial patterns from TI in the MPS

考慮到一個MPS 的實現(xiàn)（模擬結(jié)果）R與訓(xùn)練圖像集合UTI在空間結(jié)構(gòu)上具有一定的相似性，那么R中就應(yīng)該有盡可能多的模式來自UTI。令PRi表示R中的一個模式，PTI為UTI中的模式，R和UTI之間的差異就可以定義為［39］

其中D（·）為廣義距離函數(shù)。對于任一個PiR∈R，找到與它最接近（相似）的PTI∈UTI，R與UTI相對應(yīng)的模式距離之和就代表了R和UTI之間的整體不相似程度。可見，d（R，UTI）是MPS 模擬中需要最小化的目標函數(shù)。式（1）給出了基于MPS的地質(zhì)結(jié)構(gòu)三維重建過程中的全局核函數(shù)，相應(yīng)的目標函數(shù)就是實現(xiàn)式（1）的最小化，即

其中Rfinal為最終的實現(xiàn)結(jié)果。令P1TI，...，PnTI表示UTI中的所有模式。對于PR∈R和PTI∈UTI，P_TI是PR的“最近鄰”，如果它們滿足

式（3）表明P_TI是與PR最相似的模式。如果P_TI與PR的距離大于PR與它的最近鄰的距離值，則稱P_TI是PR的“近似最近鄰”。這樣，我們就可以把問題從實現(xiàn)的整體相似性轉(zhuǎn)換到空間模式相似性的求解上。

1.2 基于MPS的地質(zhì)結(jié)構(gòu)三維重建方法

Rfinal值的求取，需要同時知道R中的模式以及UTI中與它們最相似的模式。然而，式（2）中的R和Rfinal都是未知的，簡單方法難以實現(xiàn)該問題的求解。前人研究表明，期望最大化（EM，expectationmaximization） 算法非常適合解決這類優(yōu)化問題［40-41］。EM 算法采用迭代計算的思路來實現(xiàn)最小化目標函數(shù)過程，包含了兩個主要步驟：E步驟和M 步驟。E 步驟為每個PR∈R找到最相似的PTI∈UTI，即模擬網(wǎng)格中每個數(shù)據(jù)事件在UTI中的近似最相似模式；M 步驟則根據(jù)找到的包含該點的近似最近鄰的對應(yīng)點計算加權(quán)平均值，更新該點的屬性值。EM 算法實現(xiàn)基本過程如圖3 所示。由此，EM 迭代可以實現(xiàn)最小化目標函數(shù)。

圖3 EM迭代算法實現(xiàn)流程Fig.3 The implementation process of the iterative EM algorithm

地鐵工程大多采用鉆孔勘察方式來直接獲取地質(zhì)體的空間分布。這些鉆孔一般沿地鐵線路呈條帶狀分布，在線路的橫截面上分布較少。僅通過鉆孔數(shù)據(jù)難以獲得足夠的地質(zhì)結(jié)構(gòu)的空間延展信息。依托已有地質(zhì)調(diào)查成果，根據(jù)鉆孔巖芯數(shù)據(jù)所繪制的二維地質(zhì)剖面不僅再現(xiàn)了單個鉆孔內(nèi)的地質(zhì)結(jié)構(gòu)信息，還融合了地質(zhì)專家對地質(zhì)對象空間展布特征的認知。因此，本文以二維鉆孔地質(zhì)剖面作為建模數(shù)據(jù)集進行地質(zhì)結(jié)構(gòu)三維重建。

然而，二維剖面并不能直接提供地質(zhì)對象的三維空間展布模式。Extended GOSIM算法可以有效實現(xiàn)以二維鉆孔地質(zhì)剖面為UTI的三維地質(zhì)結(jié)構(gòu)重構(gòu)［38］?？紤]到破碎帶中存在不同類型的巖性，以及地層可能存在的沉積旋回，本研究以Extended GOSIM算法為基礎(chǔ)，提出了一個融合分級建模思想的多尺度MPS地質(zhì)結(jié)構(gòu)三維重建算法（見圖4）。即，首先將具有分級特征的地質(zhì)對象作為一個整體進行特征結(jié)構(gòu)的重建，獲得該地質(zhì)對象分布大致的范圍，然后，以其分布范圍為約束，再進行次級結(jié)構(gòu)的細化模擬。

圖4 融合分級建模思想的多尺度MPS地質(zhì)結(jié)構(gòu)三維重建算法實現(xiàn)過程Fig.4 The multiple-scale MPS implementation process of 3D geological structure reconstruction with hierarchical modeling idea

在本研究中，有效獲取地質(zhì)對象的三維模式是實現(xiàn)地質(zhì)結(jié)構(gòu)合理重建的基礎(chǔ)和前提。利用Extended GOSIM 算法所提出的二維剖面擴展步驟［38］，我們以隨機的模式將二維剖面擴展到一個模板大小厚度的三維剖面，然后再從擴展剖面中提取出地質(zhì)結(jié)構(gòu)的三維模式，并建立三維空間模式數(shù)據(jù)庫。地質(zhì)剖面從二維擴展到三維的過程是逐層進行的。首先將二維剖面按照實際位置放置于三維模擬網(wǎng)格中。此時，除剖面所占據(jù)的網(wǎng)格外，其余網(wǎng)格節(jié)點的屬性均未知。當對二維剖面進行擴展時，待擴展點的屬性值就從外層中與當前待擴展點相鄰的已賦值或剖面點的屬性值中隨機挑選一個。如此依次進行，在填充完一層后再填充下一層，按照設(shè)定的層數(shù)，從外層訓(xùn)練圖像向內(nèi)逐層擴展，最終得到的三維訓(xùn)練圖像是具有一個模板大小厚度的地質(zhì)剖面，如圖5 所示?？梢?，在擴展的三維地質(zhì)剖面時，并不是簡單地將二維剖面屬性復(fù)制到模擬網(wǎng)格節(jié)點中，而是以二維剖面為基礎(chǔ)的隨機賦屬性值過程。圖6給出了剖面⑤在西、中和東側(cè)3個位置上的切片，其中黑色橢圓的虛線范圍標識了3個切片的屬性值差異位置。在剖面的不同位置上，屬性分布差異主要出現(xiàn)在巖性變化的交界處。

圖5 擴展后的三維地質(zhì)剖面圖Fig.5 Three-dimensional(3D)geological cross-sections expended from 2D cross-sections

圖6 剖面⑤經(jīng)過拓展后在不同位置的切片圖（虛線橢圓標注了不同切片的差異）Fig.6 Slices in different positions that cut from different position of the extended geological cross-sections(Elliptic curves in dash line show differences between slices)

為了加快從模式庫中搜索模式，在完成二維剖面擴展后，利用移動立方體掃描三維地質(zhì)剖面，建立三維模式數(shù)據(jù)庫。所采用的移動立方體大小為建模過程中所采用的模板一致，其過程與1.1中所闡述的基本一致。所不同的是，移動立方體所得到的是三維空間結(jié)構(gòu)。在建立模式數(shù)據(jù)庫時，先利用漢明距離［42］來計算三維模式之間的相似度，通過聚類的方法，將具有較小漢明距離的模式劃分至相同的子模式類中。當所有模式都被聚類并建立索引后，就可獲得建模所需的三維模式數(shù)據(jù)庫。

2 地鐵站點地質(zhì)數(shù)據(jù)分析

本研究以廣州市軌道交通十一號線江泰路站址為研究實例（見圖7）。研究區(qū)位于曉港-赤崗碎屑巖臺地與沖-海積平原地貌交界地帶［43］。鉆孔揭示的地層主要為白堊系和第四系。第四系主要有雜填土、第四紀殘積土（Q4el）、全新統(tǒng)（Q4）和上更新統(tǒng)（Q3）河流三角洲相和海陸交互相淤泥質(zhì)土、砂層、粉質(zhì)黏土，其中全新統(tǒng)一般幾米到十幾米，第四系總厚度可達20 余米。研究區(qū)揭露的白堊系地層主要為上白堊統(tǒng)三水組康樂段（K2s1）和下白堊統(tǒng)白鶴洞組廣鋼段（K1b2）。康樂段（K2s1）屬內(nèi)陸湖泊相為主的粗砂-細砂碎屑巖建造，巖性主要為砂巖、含礫粗砂巖、礫巖，以及鈣質(zhì)粉砂巖和鈣質(zhì)泥巖。在顆粒組成上表現(xiàn)為下粗上細，泥質(zhì)膠結(jié)為主，粉細粒結(jié)構(gòu)，中厚層狀構(gòu)造。受巖體膠結(jié)物及顆粒組成不同影響，風化巖中常出現(xiàn)不同風化程度的夾層，導(dǎo)致部分地段不同風化程度的風化巖交替出現(xiàn)。廣鋼段（K1b2）為內(nèi)陸湖泊相為主的粗砂-細砂碎屑巖夾不純碳酸鹽巖，巖性主要為粉砂質(zhì)泥巖、泥質(zhì)粉砂巖、粉砂巖、細砂巖、含礫砂巖、粉砂質(zhì)泥灰?guī)r和泥灰?guī)r互層。中厚層狀構(gòu)造，泥質(zhì)、鈣質(zhì)膠結(jié)，并含微-薄層及星點狀-網(wǎng)狀石膏。該段頂部被三水組康樂段平行不整合覆蓋。

圖7 研究區(qū)基巖地質(zhì)圖（改自文獻［44］）Fig.7 Bedrock geological map of the study area(modified from ref. [44])

鉆孔巖芯還揭露出英安斑巖、安山巖等火成巖。火成巖主要從上白堊統(tǒng)三水組康樂段與下白堊統(tǒng)白鶴洞組廣鋼段之間侵入，呈柳葉刀形，后受廣三斷裂切割，呈兩段分離的扇葉型。其中，英安斑巖（ηγ53-3）巖石呈斑狀結(jié)構(gòu)，基質(zhì)為霏細結(jié)構(gòu)，斑晶25%～40%，其中石英為15%，其余為鉀長石、斜長石、角閃石、黑云母等。安山巖青灰色夾暗紫紅色，巖體破碎，巖芯呈碎塊狀夾塊狀，斑狀結(jié)構(gòu)，塊狀構(gòu)造，見氣泡，巖質(zhì)較堅硬。

前人研究表明，廣三斷裂在本站與線路相交，走向基本上為近東西向或略有偏轉(zhuǎn)，傾向南或略有偏轉(zhuǎn)，傾角50o～85o［45］。在該站址的多處鉆孔巖芯中發(fā)現(xiàn)有斷層泥、斷層角礫和碎裂巖等斷層現(xiàn)象。構(gòu)造巖以碎裂巖為主，該斷裂南東端明顯錯開白堊系地層，呈逆時針扭動方向，屬斜沖正斷裂。已有研究表明，廣三斷裂在第四紀早期曾有過多次活動，但晚更新世以來未見地表或近地表的活動跡象。

本研究在江泰路站址13 288 m2范圍內(nèi)收集了38 個鉆孔的巖芯數(shù)據(jù)，并對鉆孔數(shù)據(jù)進行分層解釋，結(jié)合淺層地震和井間電阻率CT 解釋剖面，繪制了6 條交叉剖面（見圖8）。其中，剖面⑤為鉆孔CT 測量位置，剖面⑥為聯(lián)絡(luò)剖面以約束斷裂形態(tài)和走向。

圖8 研究區(qū)地質(zhì)數(shù)據(jù)分布示意圖Fig.8 Geological data collected in the study area

3 模擬結(jié)果分析

以圖8中的鉆孔地質(zhì)剖面作為建模數(shù)據(jù)，本研究利用前述方法重建了建模區(qū)域內(nèi)的三維地質(zhì)結(jié)構(gòu)。模擬過程采用了2個尺度，每個尺度的網(wǎng)格大小是后一個尺度的8倍。在最精細尺度下，模擬網(wǎng)格數(shù)為288×352×604，各方向上網(wǎng)格精度為0.25 m。

除火成巖和破碎帶外，其余地質(zhì)體呈近水平分布，地層之間沒有出現(xiàn)倒轉(zhuǎn)和穿插現(xiàn)象（圖9）。其中，雜填土、粉土和粉質(zhì)黏土呈層狀水平連續(xù)分布于研究區(qū)域中，其平均厚度分別為4.09、3.94和3.23 m（見圖9：d～f）。中、粗砂以及淤泥質(zhì)土少量、零星分布在模型中（見圖9f）。破碎帶主要分布在建模區(qū)的東北部，主體呈東西方向延展，在中部向下尖滅于模型底部。破碎帶頂部起伏較大，高程變化范圍為-27.86～-0.58 m。破碎帶中依巖性特征，可以進一步細分為斷層角礫、斷層泥和碎裂巖。其中，斷層角礫是破碎帶主要巖性，占總體體積的81.23%；斷層泥和碎裂巖并沒有大面積、連續(xù)性分布，僅零星分布于破碎帶中（見圖10：b～d）。可見，以已建立的破碎帶整體三維結(jié)構(gòu)為約束，再利用隨機模擬方法可以很好地重構(gòu)破碎帶次級巖性分布。

圖9 研究區(qū)三維地質(zhì)模型（單位：m）Fig.9 Three-dimensional（3D)geological model of the study area(unit:m)

為驗證模型的有效性，本研究給出了地質(zhì)剖面和模擬結(jié)果中各地質(zhì)體的占比情況（見表1）。除了破碎帶與或火成巖外，各地質(zhì)體在模型中整體占比與在剖面中的占比基本接近。破碎帶所包括的斷層角礫、斷層泥和斷層巖整體分布在建模區(qū)的北部地區(qū)，分布并不均勻，模型占比均小于其在剖面整體中的占比。此外，在剖面中只出現(xiàn)在破碎帶周圍的火成巖，在建模結(jié)果中也表現(xiàn)出相似特性，其占比也遠小于剖面中整體中占比。

表1 研究區(qū)三維模型與地質(zhì)剖面中各地質(zhì)對象占比Table 1 The proportions of different geological bodies in cross-sections and the final 3D model %

從鉆孔所揭露的巖性看，工程力學(xué)性質(zhì)相對較差的主要為破碎帶和風化火成巖。站址與破碎帶模型的空間位置關(guān)系（圖10a）表明：站址范圍內(nèi)涵蓋了部分破碎帶，破碎帶體積約占站址范圍的3.2%。站址的中東部底板與破碎帶相交。為更好分析站址底部巖性分布情況，圖11 給出了站址底部以下20 m 范圍內(nèi)的巖性分布情況。在100 m 以東區(qū)域，站址底板分布有大量斷層角礫、斷層泥和火成巖。在70 m 以東區(qū)域，站址下方具有厚度不等的破碎帶分布，以西區(qū)域以砂巖為主?？梢?，在進行站址底板設(shè)計過程中，應(yīng)結(jié)合站址底部巖性變化特征，分段設(shè)計底板構(gòu)筑物結(jié)構(gòu)。

圖10 破碎帶模型及地鐵站址關(guān)系（單位：m）Fig.10 The spatial relationship of the metro station and fault zone model(unit:m)

圖11 站址底部以下20 m內(nèi)巖性分布情況（單位：m）Fig.11 Lithology distribution within 20 m beneath the bottom of the metro station(unit:m)

4 結(jié) 論

針對地鐵工程應(yīng)用中數(shù)據(jù)集中、但難以獲取地質(zhì)對象空間分布模式的難題，本論文以二維鉆孔地質(zhì)剖面為建模數(shù)據(jù)集，依托Extended GOSIM算法，提出了一種融合分級建模思想的多尺度MPS 地質(zhì)結(jié)構(gòu)三維重建方法。廣州某地鐵站址三維地質(zhì)模型的構(gòu)建實例表明：分級建模的思想可以很好地實現(xiàn)斷裂中不同類型巖性的空間分布，也可以將其推廣到沉積旋回地層的三維重建中；所提出的三維建模方法能充分利用有限數(shù)據(jù)內(nèi)所蘊含的地質(zhì)對象間空間關(guān)系，有效識別、提取和重建地質(zhì)對象的三維空間模式?？梢?，基于MPS的隱函數(shù)建模方法能構(gòu)建出高精度的三維地質(zhì)結(jié)構(gòu)，可以為地鐵工程設(shè)計、施工提供精細數(shù)據(jù)基礎(chǔ)。以隨機方式將二維地質(zhì)剖面擴展成三維剖面在一定程度上可以解決三維地質(zhì)結(jié)構(gòu)難以獲取的問題，由于擴展范圍只有一個模板大小，簡單隨機外推方法所獲得的三維地質(zhì)剖面與二維地質(zhì)剖面間的差異過小，也難以獲得地質(zhì)對象的全局結(jié)構(gòu)特征。