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

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

0 引言

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

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

1 Markov鏈模擬方法原理

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

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

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

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

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

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

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

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

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

1）對待模擬區(qū)根據(jù)地形條件及鉆孔數(shù)量進行區(qū)域網(wǎng)格劃分；

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

圖2 間隔鉆孔數(shù)據(jù)修正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區(qū)域多個Monte-Carlo模擬圖像；

4）提取模擬結果虛擬Y2孔數(shù)據(jù)，并與真實Y2孔數(shù)據(jù)進行對比分析統(tǒng)計吻合率，選取吻合度最高的實現(xiàn)結果作為修正圖像；

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

6）結合兩次最優(yōu)模擬結果對重疊區(qū)進行數(shù)據(jù)修正，依次進行得到實驗區(qū)完整的模擬圖像.

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

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

3.1 工程概況

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

3.2 轉移狀態(tài)矩陣建立

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

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

3.3 物探數(shù)據(jù)統(tǒng)計修正

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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

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

[2]Yarus J M，Chambers R L.隨機建模和地質統(tǒng)計學—原理、方法和實例研究[M].穆龍新，陳亮，譯.北京：石油工業(yè)出版社，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] 李軍，熊利平，方石，等.基于多尺度數(shù)據(jù)融合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]沈洪濤，郭乃川.地質統(tǒng)計學反演技術在超薄儲層預測中的應用[J].地球物理學進展，2017，32（1）：248-253.

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

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

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