基于二維馬爾可夫鏈的武漢長江公鐵隧道地層識別

程利力， 陳 健， 陳 睿， 魏林春

(1. 華中科技大學(xué) 土木與水利工程學(xué)院， 湖北 武漢 430074； 2. 中國科學(xué)院武漢巖土力學(xué)研究所， 湖北 武漢 430071； 3. 上海隧道工程有限公司， 上海 200032)

巖土體在長期形成過程中受到地質(zhì)成分、沉積環(huán)境等因素的作用，導(dǎo)致地層分布存在一定的不確定性，即地層變異性[1]；同時地下工程不同于結(jié)構(gòu)工程，這種不確定性很難進行準確分析，傳統(tǒng)方法往往根據(jù)已有鉆孔數(shù)據(jù)進行簡單線性插值從而獲取未知區(qū)域的地層分布，這種方法簡單方便但計算效果有待提高[2]。

目前，一些學(xué)者提出了通過概率預(yù)測方法對地層變異性進行模擬。譚卓英等[3]采用一種地層識別方法，對鉆孔指標進行有效檢測。趙磊等[4]基于支持向量機的方法，根據(jù)測井?dāng)?shù)據(jù)建立相關(guān)地層識別模型。Cao等[5]利用貝葉斯方法對靜力觸探試驗數(shù)據(jù)進行識別，從而統(tǒng)計了土層的厚度和邊界。劉毅等[6]對地質(zhì)數(shù)據(jù)采用主成分分析法進行地層識別，從而獲取有效的地質(zhì)信息。Ching等[7]基于靜力觸探試驗數(shù)據(jù)，采用小波變換識別最大值，從而對土層進行劃分。這些研究主要針對單個鉆孔進行地層識別，然而對于多鉆孔的地層識別研究較少。Qi等[8]提出了一種基于有限鉆孔數(shù)據(jù)的多變量自適應(yīng)回歸樣條法，從而可以預(yù)測較為合理的地質(zhì)剖面。Li等[9]基于二維馬爾可夫鏈模型，為研究鉆孔之間的地層識別提供了一種較為簡便的方法；其中，地層各方向的轉(zhuǎn)移概率是十分重要的模型參數(shù)，縱向轉(zhuǎn)移概率矩陣可以由鉆孔資料獲得，而橫向轉(zhuǎn)移概率矩陣受鉆孔數(shù)據(jù)的限制，往往不足以直接估計[10]。

基于上述討論，為推進地層識別在巖土工程中的應(yīng)用，本文采用二維馬爾可夫鏈理論模型，間接獲取了橫向轉(zhuǎn)移概率矩陣，給出了地層識別的算法流程圖；最后，以武漢長江公鐵隧道武昌岸右線部分鉆孔資料為例，進行了有效的地層識別并驗證了該方法的可靠性，為以后實際地鐵工程施工提供有效的地質(zhì)參考資料，從而進行相應(yīng)的工程分析。

1 工程基本概況

1.1 工程概況

武漢長江公鐵隧道是國內(nèi)第一條橫貫長江的公鐵兩用隧道，隧道包括左右兩線，單線隧道長2590 m，由武昌工作井始發(fā)，經(jīng)過武昌堤防后穿越長江，到達漢口工作井。圖1為武漢長江公鐵隧道地理位置圖。

圖1 武漢長江公鐵隧道地理位置

1.2 水文地質(zhì)條件

武漢長江公鐵隧道位于長江河床及兩岸一級階地區(qū)。長江江底環(huán)境復(fù)雜，其中水、泥、沙較多。長江兩岸建筑物密集，路網(wǎng)發(fā)達，局部有暗塘、暗溝分布，地下管線、排水管道密布，其中盾構(gòu)段穿越粉質(zhì)黏土、粉細砂土、泥巖等復(fù)合土層，地質(zhì)環(huán)境復(fù)雜，屬于中等復(fù)雜程度，工程重要性程度為一級，因此，地質(zhì)資料的準確與否對隧道工程施工安全十分重要。圖2為武漢長江公鐵隧道水文地質(zhì)圖。

圖2 武漢長江公鐵隧道水文地質(zhì)

1.3 地質(zhì)鉆孔信息

由于漢口岸一帶建筑物密集，地下管線眾多，對場地環(huán)境擾動過大，因此，本文采用武漢長江公鐵隧道武昌岸右線部分地質(zhì)鉆孔資料；按城市軌道交通巖土工程勘察規(guī)范[11]，該工程為中等復(fù)雜場地，地下工程勘探點間距宜取30～50 m，因此選取區(qū)域大小為 186 m×60 m的工程場地；同時考慮到盾構(gòu)掘進一環(huán)約為2 m，故選取大小為2 m×1 m的單元網(wǎng)格進行劃分。該工程地質(zhì)鉆孔的相對位置如圖3所示，為簡便起見，將鉆孔編號為Z1～Z6，將各地質(zhì)鉆孔進行投影，得到各鉆孔的地層信息剖面如圖4所示，該地層主要有6種土體類

圖3 鉆孔相對位置

圖4 各鉆孔反映的地層信息

型，其中A為雜填土，B為3-1粉質(zhì)粘土，C為3-2粉質(zhì)粘土，D為粉細砂土，E為中粗砂土，F(xiàn)為泥巖。

2 二維馬爾可夫鏈模型

2.1 二維馬爾可夫鏈的基本原理2.1.1 二維馬爾可夫鏈的條件概率

地質(zhì)沉積過程可以看作馬爾可夫過程[12]，利用馬爾可夫鏈的前提條件是土體轉(zhuǎn)移狀態(tài)符合一階馬爾可夫性；巖土體在橫向和縱向上的沉積不同，假設(shè)分別由兩個不同的一維馬爾可夫鏈組成，然后將其按Elfeki等[13]的方法進行耦合，其在二維空間的分布如圖5所示，左右兩側(cè)單元可以由鉆孔資料獲取，按行從左到右依次進行模擬得到未知單元的狀態(tài)，其中具體步驟見文獻[14]。在已知單元(i-1，j)，(i，j-1)，(Nx，j)的狀態(tài)分別為Zp，Zq，Zw的情況下，單元Zn的條件概率為：

圖5 二維馬爾可夫鏈狀態(tài)序列

p(Si,j=Zn|Si-1,j=Zp,Si,j-1=Sq,SNx,j=Zw)

(1)

2.1.2 縱向轉(zhuǎn)移概率矩陣

縱向轉(zhuǎn)移計數(shù)矩陣可由鉆孔資料統(tǒng)計獲得，記為Cv：

(2)

則縱向轉(zhuǎn)移概率矩陣為Rv為：

(3)

(4)

2.1.3 橫向轉(zhuǎn)移概率矩陣

Walther相序定律[15]認為同一剖面地層在各方向的沉積順序是相同的，但橫向的沉積規(guī)模一般大于縱向的沉積規(guī)模，不能直接采用縱向轉(zhuǎn)移概率代替橫向轉(zhuǎn)移概率。通常，地層沉積往往由走向、傾向、傾角三要素表示，而其中地層傾角的作用可量化為橫向與縱向延伸長度之比，記為K[16]，由于地層延伸長度在橫向遠大于縱向，于是在橫向上地層向自身發(fā)生轉(zhuǎn)移的概率比較大，則橫向轉(zhuǎn)移計數(shù)矩陣可表示Ch為：

(5)

橫向轉(zhuǎn)移概率矩陣Rh為：

(6)

(7)

2.1.4K值估計

地層在橫向與縱向延伸長度之比K可以利用地質(zhì)鉆孔進行估計獲得，其基本思路類似于極大似然估計，不同K值得到的土體類型不同，將各鉆孔區(qū)間的轉(zhuǎn)移似然度相乘即為總的土體轉(zhuǎn)移似然度Q；當(dāng)Q取得最大值時，K值最優(yōu)。

Q=Q1…QN-1

(8)

Qz=R(Zi,1=Zm,Zj,1=Zo)R(Zi,2=Zm+1,
Zj,2=Zo+1)…R(Zi,Ny=Zn,Zj,Ny=Zp)

(9)

式中：Zm,…,Zn為第z個地質(zhì)鉆孔剖面區(qū)間內(nèi)對應(yīng)的土體類型；Zo,…,Zp為第z+1個地質(zhì)鉆孔剖面區(qū)間內(nèi)對應(yīng)的土體類型；Qz為鉆孔z與鉆孔z-1兩側(cè)鉆孔之間的土體類型轉(zhuǎn)移發(fā)生概率的似然度；Q為鉆孔1到鉆孔N所有相鄰兩鉆孔之間轉(zhuǎn)移發(fā)生概率的似然度相乘之積。圖6為各地質(zhì)鉆孔剖面。

圖6 地質(zhì)鉆孔剖面

2.2 二維馬爾可夫鏈的計算流程

本文運用MATLAB軟件，結(jié)合理論公式推導(dǎo)，通過二維馬爾可夫鏈進行模擬，相應(yīng)的流程如圖7所示。

(1)確定建模區(qū)域，讀取鉆孔數(shù)據(jù)，確定地質(zhì)剖面土體類型；

(2)利用適當(dāng)?shù)牟蓸娱g隔對橫向和縱向區(qū)域進行離散化，得到馬爾可夫鏈網(wǎng)格；

(3)將地層信息映射到相應(yīng)的地質(zhì)單元，計算橫向和縱向轉(zhuǎn)移概率矩陣；

(4)估計地質(zhì)剖面內(nèi)橫向和縱向延伸長度之比K；

(5)對地質(zhì)剖面進行二維馬爾可夫鏈模擬得到各地質(zhì)單元的土體類型；

(6)輸出模擬得到的地質(zhì)剖面結(jié)果。

3 結(jié)果分析與討論

3.1 轉(zhuǎn)移概率矩陣估計3.1.1 縱向轉(zhuǎn)移概率矩陣

本文針對該工程區(qū)域的6個鉆孔地質(zhì)信息，以鉆孔Z1，Z3，Z4，Z6作為約束鉆孔，鉆孔Z2，Z5作為參照鉆孔，由式(2)可得縱向轉(zhuǎn)移計數(shù)矩陣，如表1所示；將表1的計算結(jié)果代入式(3)可得相應(yīng)的縱向轉(zhuǎn)移概率矩陣，如表2所示。

表1 縱向轉(zhuǎn)移計數(shù)矩陣

表2 縱向轉(zhuǎn)移概率矩陣

3.1.2 橫向轉(zhuǎn)移概率矩陣

在求解橫向轉(zhuǎn)移概率矩陣之前需先確定K值，分別取1，2，…，20等20個不同的K值，由式(5)可得橫向轉(zhuǎn)移計數(shù)矩陣，如表3所示；其中，K取12時對應(yīng)的橫向轉(zhuǎn)移概率矩陣如表4所示。

表3 橫向轉(zhuǎn)移計數(shù)矩陣

表4 橫向轉(zhuǎn)移概率矩陣(K=12)

利用式(8)(9)依次得到不同K值相應(yīng)的轉(zhuǎn)移似然度Q，并繪制Q隨K的變化曲線，如圖8所示。

圖8 lg(Q)隨K的變化曲線

由圖8可知，土體轉(zhuǎn)移概率似然度Q隨K值的增大先增大，然后逐漸趨于穩(wěn)定，其中， 當(dāng)K為12時，似然度Q最大。

3.2 地層識別結(jié)果

采用上述Z1，Z3，Z4，Z6作為約束鉆孔進行二維馬爾可夫鏈模擬，以鉆孔Z2，Z5作為參照，用于對以下地質(zhì)剖面實現(xiàn)的準確性評估，運用蒙特卡羅思想[17]，對不同K值進行500次實現(xiàn)；由于K=4與K=20之間似然度比較接近，因此圖9給出了K取1，4，8，12，16，20時的一次典型地層模擬結(jié)果。

圖9 不同K值時的一次典型地層模擬

為描述不同K值局部鉆孔Z2，Z5所對應(yīng)的土體類型效果，采用模擬地層與實際地層分布的誤差率來進行定量評價：

(10)

式中：Ni為i號鉆孔單元總數(shù)；Zi為i號鉆孔模擬單元與實際單元的誤差個數(shù)。

圖10顯示了預(yù)測鉆孔不同K值的誤差率。由圖9，10可知，當(dāng)K=1，4，8時，在橫向上同一土體類型轉(zhuǎn)移概率較小，地層鉆孔之間缺乏較好的橫向延伸性，因此預(yù)測鉆孔Z2，Z5的模擬效果不佳，誤差率較大；當(dāng)K=16，20時，地層在橫向上同一土體類型發(fā)生轉(zhuǎn)移的概率會增大，因此不同土體類型在橫向上很難進行轉(zhuǎn)移，這對預(yù)測鉆孔Z2，Z5的模擬效果也會產(chǎn)生一定影響。由上述結(jié)果也揭示了K為12時地層識別結(jié)果相對較優(yōu)。

圖10 預(yù)測鉆孔不同K值的誤差率

其中，當(dāng)K=12時，表5列舉了鉆孔Z2，Z5處20，40，60 m深度下土體單元分別為不同土體類型的概率值。

由表5可知，該方法將概率最大的土體類型作為該單元的土體類型。同時，根據(jù)K=12時的地層模擬結(jié)果，表6列舉了部分單元位置相應(yīng)的土體類型，其可作為虛擬鉆孔用于工程地質(zhì)資料分析。

表5 土體單元概率

表6 單元土體類型

3.3 結(jié)果驗證

利用二維馬爾可夫鏈模擬得到估計鉆孔Z2，Z5的地層信息；同時采用鉆孔Z1，Z3，Z4，Z6進行簡單線性插值得到鉆孔Z2，Z5的地層信息，進而對比分析兩種方法識別結(jié)果的誤差率；圖11給出了當(dāng)K=12時，鉆孔Z2，Z5模擬結(jié)果與實測數(shù)據(jù)的地層分布情況。

由圖11可知，鉆孔Z2，Z5利用線性插值法的誤差率分別為15%，10%；而利用二維馬爾可夫鏈模擬的誤差率為8.3%，6.6%；相比之下，二維馬爾可夫鏈模擬有效降低了地層識別結(jié)果的誤差率，使得地層識別結(jié)果更加接近于真實情況，從而驗證了該方法的可靠性。

4 結(jié) 論

本文利用二維馬爾可夫鏈理論模型，以武漢長江公鐵隧道右線部分鉆孔資料為例，進行了有效的地層識別，分析得到以下結(jié)論：

(1)在該工程中，基于有限鉆孔資料，運用二維馬爾可夫鏈模型進行有效的地層識別，并較簡單插值法的誤差率有所降低，但各鉆孔之間長度不宜過長，否則會對地層識別效果產(chǎn)生影響。

(2)采用二維馬爾可夫鏈模型進行地層識別，充分考慮了地層變異性對地層識別的影響，可為實際工程提供較為準確的地層識別結(jié)果，特別是可為彌補地質(zhì)勘查孔位不足時提供較為詳細的工程地質(zhì)資料。