基于Markov鏈的隧道圍巖破碎狀態模擬研究

葉 亭，劉春原，李金龍，宋健康，徐良玉

（1.河北工業大學 土木與交通學院，天津 300401；2.河北省土木工程技術研究中心，天津 300401）

0 引言

地層變異性是確定隧道施工前圍巖狀態，選擇開挖方式的主要影響因素.地質鉆孔數據提供的地質信息是有限的，鉆孔之間的地質體分布狀態存在較大的隨機不確定性，傳統方法無法反映巖體的隨機分布特征[1-2].Halda和Tang在1949年首次將Markov隨機過程應用于沉積地層垂向疊置樣式分析[3]，Rumbein和Schercr提出了采用固定步長法分析土層的轉移概率，Ioannou在1984年將地層隨機轉移狀態法用于垂直方向巖土層分析.劉振峰以Markov鏈為基礎從轉移計數矩陣入手，成功的實現沿工程縱斷面方向的二維拓展[4].Elfeki在文獻中討論了利用兩個離散的一維Markov鏈組成聯合概率來模擬地層剖面的二維分布模擬的方法[5]，謝式千等將劉振峰的二維Markov鏈模擬方法進一步改進并應用到實際工程.劉振峰、楊長春將模擬退火算法引入Markov鏈模型中，在得到巖相轉移概率矩陣的基礎上對模擬結果進行擾動，直到得出較為滿意的圖像[6-7].李軍在2010年將尺度數據融合的方法引入Markov預測模型中，在計算過程中協調了巖心數據和物探數據進行二維地層巖相儲層模擬分析[8].然而，針對Markov鏈地層預測模型的研究主要是對大范圍區域地層巖相組合進行模擬預測，在隧道圍巖破碎狀態的研究未見報道.

為此，提出了間隔鉆孔數據修正的模擬方法，對鉆孔數據模擬結果進行互相修正降低了地層變異不可控的風險，并將該Markov鏈地層預測模型應用于地質體風化破碎狀態的預測.為隧道工程的設計及施工提供了有效的數據支持.

1 Markov鏈模擬方法原理

Markov鏈是1906年數學家Markov提出的“無后效性”的隨機過程，即未來狀態的發展與先前狀態無關，只取決于現在狀態.其數學表達式為

式中：X0，Xn，Xn+1為Markov隨機過程中第1個，第n個以及第n+1個位置的狀態；h，j，k分別為狀態X0，Xn，Xn+1的取值；Pr為事件發生概率[15].

由此可以看出Markov鏈是有方向性的，這一特性和沉積過程的方向性相吻合，這使得利用Markov鏈預測地層狀態成為可能[4].根據得到的地質鉆孔巖芯數據，設狀態空間為S{sk∈S，k=1，2，…，n}，同時假定已知地層下邊界和左邊界的鉆孔巖心數據，右邊界作為條件數據.統計計算出其相應方向上不同巖性的轉移狀態分布情況，求出二維空間中某位置出現sk狀態的概率，見圖1.網格（i，j）處出現狀態Sk的概率計算公式見式（2）

式中：用上標表示模擬隨機域，1為小尺度模擬隨機域即鉆孔巖芯統計域，2為大尺度模擬隨機域即物探數據統計域；d為大尺度網格狀態（根據物探波速統計的巖性狀態）；h為垂直方向，v為水平方向；i，j分別為狀態模擬網格的位置；Nx為條件網格位置；sl，sm分別為水平向、垂向已知狀態；sq為條件狀態；C為歸一化參數[5]；Pr為概率.

在參數修正過程中可以利用高斯轉換模型進行統計域之間數據的轉換[8]，具體轉換公式見式（3）.

其中：φ（x），φ（x）代表不同尺度數據之間的概率融合轉換映射過程；ml代表一個待模擬大尺度場中的小網格數；ε～N（μ，σ2）為準確度控制函數，在地質統計學中一般用正態分布，其中均值μ及方差σ2根據待模擬區域物探統計數據得出；xl，xs為不同尺度網格取值.

圖1 條件二維Markov鏈網格模擬計算示例圖Fig.1 The example of two-dimensional Markov chain mesh simulation

2 Markov鏈圍巖剖面模型的建立

用Markov鏈模型進行模擬的過程中，利用已知鉆孔實測數據，將區域巖層的破碎狀態分布統計規律[14]轉化為隨機轉移狀態矩陣，采用間隔鉆孔信息修正的方法，進行模擬計算，具體步驟如下：

1）對待模擬區根據地形條件及鉆孔數量進行區域網格劃分；

2） 選取符合條件的3個相鄰鉆孔，編號Y1，Y2，Y3，經χ2檢驗符合Markov條件后，將Y1鉆孔巖芯破碎狀態數據轉化為區域隨機轉移狀態矩陣，Y3號孔作為條件數據，進行隨機模擬計算，如圖2所示；

圖2 間隔鉆孔數據修正Markov隧道圍巖破碎帶分布預測示意圖Fig.2 Sketch map of the distribution of the surrounding rock broken zone in the Markov tunnel with interval data

3）對計算結果應用Monte-Carlo法進行大量次的隨機模擬賦值，形成Y1～Y3區域多個Monte-Carlo模擬圖像；

4）提取模擬結果虛擬Y2孔數據，并與真實Y2孔數據進行對比分析統計吻合率，選取吻合度最高的實現結果作為修正圖像；

5）將Y2孔作為原始數據，下一鉆孔作為條件數據進行遞進模擬實現，得到Y2～Y4區域模擬圖像，通過對比計算虛擬Y3孔與真實Y3孔的模擬結果，選取模擬效果最高的圖像；

6）結合兩次最優模擬結果對重疊區進行數據修正，依次進行得到實驗區完整的模擬圖像.

通過間隔鉆孔數據依次遞進計算，互相修正，減小地層變異性對圍巖破碎帶轉移狀態的影響，得到更為真實可靠的圍巖破碎帶分布數據.

3 隧道圍巖破碎帶模擬分析

3.1 工程概況

張石高速某隧道工程項目隧址區位于河北省保定市太行山中部山區，勘測路線起于淶源縣孤山北至紫荊關北.除局部緩坡及溝谷處堆積第四系全新統坡崩積層及坡洪積層外，大部分地段基巖裸露.模擬區RK73+450～RK77+450里程段地層巖性為燕山期花崗侵入巖，風化破碎程度不均，節理較發育，未發現斷層.

3.2 轉移狀態矩陣建立

首先選擇標高為600～700 m，里程RK73+450～R K77+450為模擬區段，選取其中9個鉆孔（ZKSD12-3～ZKSD12-14號）巖芯數據.取垂直方向步長為1 m，水平方向步長為5 m，模擬剖面為100×200的矩形模擬系統進行模擬計算.按層序生長方向對鉆孔進行巖芯破碎狀態進行轉移記數統計，采用χ2檢驗法對Markov性進行檢驗.統計結果為：垂向向下χ2=98.437 5，垂向向上χ2=96.876 2，其值均大于置信度為0.014χ2分布的臨界值χ20.014=13.277.由此可知該區域圍巖破碎分布序列具有較強的Markov性[9]，可以帶入隨機預測模型進行圍巖破碎帶轉移狀態的模擬預測.

按照Markov鏈模型對水平向巖相的轉移記數矩陣進行計算.從地震剖面上測得了層序不同方向上的傾角θ為14.5°并計算該區域的側向的延伸度之比14.5°傾角余切為3.866 7，得到水平向轉移概率矩陣[7]，統計計算的不同方向巖芯破碎狀態序列的轉移記數矩陣和轉移概率矩陣如表1所示.

3.3 物探數據統計修正

統計物探數據，以地震剖面劃分不同波速區間根據測井、巖芯資料統計得到不同波速區間的隧道圍巖破碎帶出現的概率分布，其地震縱波統計概率分布見圖3.統計并計算地震剖面圖不同尺度融合的精度控制函數ε～N（0.45，1.42），將物探數據帶入高斯分布隨作為約束條件，加入到條件化二維Markov鏈模型中[8，12-13].使用Matlab數據處理軟件經過數學編程，進行里程K46+625～K47+225區段隧道圍巖破碎帶模擬建模.

3.4 圍巖破碎帶模擬結果分析

Markov模擬網格區域垂向選取標高600～700 m，縱向選取里程RK73+450～RK77+450進行網格劃分然后沿隧道走向分區段進行，結果見圖4.其中a） 為里程段RK73+450～RK74+450，b） 為里程段RK74+450～RK75+450，c） 為里程段 RK75+450～RK76+450，d） 為里程段 RK76+450～RK77+450，由模擬圖像所反映出的圍巖破碎帶分布的區域性及連續性特征與根據鉆孔資料和地質超前預報資料分析得出的隧道圍巖破碎帶分布規律是一致的.對比隧道勘察數據，該區模擬準確率80.7%～89.6%.

由模擬圖a） 看出，在隧道孔頂標高635 m，孔低標高623 m.隧道垂直孔徑12 m，隧道在進口里程RK73+450處穿越中風化圍巖區，在里程RK73+950處開始圍巖破碎程度出現較大改變，其中在里程RK73+900～RK74+030區段隧道圍巖孔頂上方標高626～635 m，預測出大片中風化較破碎花崗巖.此處按照設計勘察報告為II級圍巖，圍巖完整性較好.

圖3 不同波速（縱波）區間各巖石破碎帶出現的概率條形統計圖Fig.3 The probability bar graph of each rock fracture zone with different wave velocity

表1 不同方向圍巖破碎帶轉移記數矩陣及相應的轉移概率矩陣Tab.1 Different directions surrounding rock with transfer count matrix and corresponding transtion-probablity matrix

圖4 RK76+450～RK77+450區段隧道圍巖破碎區模擬結果Fig.4 Simulation results of RK76+450～RK77+450 section tunnel surrounding rock broken zone

結合施工日志及隧道工程變更通知，在里程RK73+998～RK74+018區段發生隧道洞頂碎石塌落，圍巖由原設計II級改為Ⅲ級，該模型更準確的實現了隧道圍巖狀態模擬預測預報.

由模擬圖c） 可以看出模擬里程RK75+450～RK76+450區段，在里程RK75+800至RK76+050處隧道開始穿越圍巖中風化破碎區，且該區圍巖在標高618～642 m區域多處出現水平連續性花崗巖破碎帶，圍巖極不穩定.該區地質勘察報告設計為II級微風化花崗巖，圍巖較穩定.這是由于傳統地質剖面在里程RK75+560處勘測結果存在盲區.工程實例證明，通過圍巖分布模擬預測模型實現了隧道圍巖破碎帶有效的預測，修正了初始隧道圍巖地質剖面圖，解決了現有隧道地質勘測方法存在盲區的問題.該模型為隧道綜合風險評價提供了數據支持，對于隧道安全建設具有重要意義.

4 結論

張石高速隧道工程圍巖模擬預測的應用實例表明，對于采用間隔鉆孔數據修正的Markov鏈模型進行隧道圍巖破碎帶的模擬預測是可行的，可以反映隧道圍巖展布特征.主要結論如下：

1）通過間隔鉆孔數據依次遞進計算，互相修正，可以減小地層變異性對地層巖性轉移狀態的影響，得到更為真實可靠的巖性分布數據.

2）采用間隔鉆孔數據修正的Markov鏈模型進行二維隧道圍巖破碎帶隨機模擬，結合隧道開挖掌子面施工日志資料，預測準確率達80.7%～89.6%.

3）應用改進的Markov鏈模型進行圍巖破碎帶模擬，可以修正隧道圍巖地質剖面圖，解決了現有隧道地質勘測方法存在盲區的問題.

[1] 王仁銥，胡光道.線形地質統計學[M].北京：地質出版社，1988.

[2]Yarus J M，Chambers R L.隨機建模和地質統計學—原理、方法和實例研究[M].穆龍新，陳亮，譯.北京：石油工業出版社，2000.

[3]Schwarzacher W.沉積模型和定量地層學[M].徐桂榮，譯.北京：地質出版社，1984.

[4] 劉振峰，郝天珧，楊長春.基于鏈模型的儲層巖相隨機模擬[J].地球物理學進展，2003，18（4）：666-669.

[5] Elfeki A，Dekking A.Markov chain model for subsurface characterization[J].Theory and applications of Mathematical Geology，2001，33（5）：568-589.

[6] 劉振峰，郝天珧，楊長春.沉積模型和儲層隨機建模[J].地球物理學進展，2003，18（3）：519-523.

[7] 劉振峰，郝天珧，方輝.用鏈模型隨機模擬儲層巖相空間展布[J].石油學報，2005，26（5）：57-60.

[8] 李軍，熊利平，方石，等.基于多尺度數據融合Markov鏈模型的巖性隨機模擬[J].石油學報，2010，31（1）：74-77.

[9] 劉振峰．基于MRF模型的油氣儲層屬性隨機建模方法研究[D].北京：中國科學院研究生院，2004.

[10]Haldar A，Tang W H.Uneorminty analysis of relative density[J].Journal of Gentechnical Engineering ASCE，1979，110（4）：525-530.

[11]Norberg T，Rosen L，Baran A.On modeling discrete geological structures as Markov random fields[J].Mathematical Geology，2002，34（1）：63-77.

[12]沈洪濤，郭乃川.地質統計學反演技術在超薄儲層預測中的應用[J].地球物理學進展，2017，32（1）：248-253.

[13]侯斌，陳波，薄永德，等.基于地質統計學反演的薄互砂巖儲層預測[J].復雜油氣藏，2016，9（4）：12-15.

[14]王宇.基于巖性分析的公路隧道圍巖動態分級研究[D].長沙：長沙理工大學，2015.

[15]孫安黎，楊陽，伍焓熙，等.基于馬爾可夫過程的項目風險評價方法的研究[J].理論與算法，2016，24（1）：36-37.