基于二維馬爾可夫鏈的武漢長(zhǎng)江公鐵隧道地層識(shí)別

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

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

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

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

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

1 工程基本概況

1.1 工程概況

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

圖1 武漢長(zhǎng)江公鐵隧道地理位置

1.2 水文地質(zhì)條件

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

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

1.3 地質(zhì)鉆孔信息

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

圖3 鉆孔相對(duì)位置

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

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

2 二維馬爾可夫鏈模型

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

地質(zhì)沉積過(guò)程可以看作馬爾可夫過(guò)程[12]，利用馬爾可夫鏈的前提條件是土體轉(zhuǎn)移狀態(tài)符合一階馬爾可夫性；巖土體在橫向和縱向上的沉積不同，假設(shè)分別由兩個(gè)不同的一維馬爾可夫鏈組成，然后將其按Elfeki等[13]的方法進(jìn)行耦合，其在二維空間的分布如圖5所示，左右兩側(cè)單元可以由鉆孔資料獲取，按行從左到右依次進(jìn)行模擬得到未知單元的狀態(tài)，其中具體步驟見(jiàn)文獻(xiàn)[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)移計(jì)數(shù)矩陣可由鉆孔資料統(tǒng)計(jì)獲得，記為Cv：

(2)

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

(3)

(4)

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

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

(5)

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

(6)

(7)

2.1.4K值估計(jì)

地層在橫向與縱向延伸長(zhǎng)度之比K可以利用地質(zhì)鉆孔進(jìn)行估計(jì)獲得，其基本思路類似于極大似然估計(jì)，不同K值得到的土體類型不同，將各鉆孔區(qū)間的轉(zhuǎn)移似然度相乘即為總的土體轉(zhuǎn)移似然度Q；當(dāng)Q取得最大值時(shí)，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個(gè)地質(zhì)鉆孔剖面區(qū)間內(nèi)對(duì)應(yīng)的土體類型；Zo,…,Zp為第z+1個(gè)地質(zhì)鉆孔剖面區(qū)間內(nèi)對(duì)應(yīng)的土體類型；Qz為鉆孔z與鉆孔z-1兩側(cè)鉆孔之間的土體類型轉(zhuǎn)移發(fā)生概率的似然度；Q為鉆孔1到鉆孔N所有相鄰兩鉆孔之間轉(zhuǎn)移發(fā)生概率的似然度相乘之積。圖6為各地質(zhì)鉆孔剖面。

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

2.2 二維馬爾可夫鏈的計(jì)算流程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表3 橫向轉(zhuǎn)移計(jì)數(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時(shí)，似然度Q最大。

3.2 地層識(shí)別結(jié)果

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

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

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

(10)

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

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

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

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

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

表5 土體單元概率

表6 單元土體類型

3.3 結(jié)果驗(yàn)證

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

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

4 結(jié) 論

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

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

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