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

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

0 引言

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

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

1 Markov鏈模擬方法原理

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

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

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

式中：用上標(biāo)表示模擬隨機(jī)域，1為小尺度模擬隨機(jī)域即鉆孔巖芯統(tǒng)計域，2為大尺度模擬隨機(jī)域即物探數(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ù)修正過程中可以利用高斯轉(zhuǎn)換模型進(jìn)行統(tǒng)計域之間數(shù)據(jù)的轉(zhuǎn)換[8]，具體轉(zhuǎn)換公式見式（3）.

其中：φ（x），φ（x）代表不同尺度數(shù)據(jù)之間的概率融合轉(zhuǎn)換映射過程；ml代表一個待模擬大尺度場中的小網(wǎng)格數(shù)；ε～N（μ，σ2）為準(zhǔn)確度控制函數(shù)，在地質(zhì)統(tǒng)計學(xué)中一般用正態(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鏈模型進(jìn)行模擬的過程中，利用已知鉆孔實測數(shù)據(jù)，將區(qū)域巖層的破碎狀態(tài)分布統(tǒng)計規(guī)律[14]轉(zhuǎn)化為隨機(jī)轉(zhuǎn)移狀態(tài)矩陣，采用間隔鉆孔信息修正的方法，進(jìn)行模擬計算，具體步驟如下：

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

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

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

3）對計算結(jié)果應(yīng)用Monte-Carlo法進(jìn)行大量次的隨機(jī)模擬賦值，形成Y1～Y3區(qū)域多個Monte-Carlo模擬圖像；

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

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

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

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

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

3.1 工程概況

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

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

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

按照Markov鏈模型對水平向巖相的轉(zhuǎn)移記數(shù)矩陣進(jìn)行計算.從地震剖面上測得了層序不同方向上的傾角θ為14.5°并計算該區(qū)域的側(cè)向的延伸度之比14.5°傾角余切為3.866 7，得到水平向轉(zhuǎn)移概率矩陣[7]，統(tǒng)計計算的不同方向巖芯破碎狀態(tài)序列的轉(zhuǎn)移記數(shù)矩陣和轉(zhuǎn)移概率矩陣如表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ù)學(xué)編程，進(jìn)行里程K46+625～K47+225區(qū)段隧道圍巖破碎帶模擬建模.

3.4 圍巖破碎帶模擬結(jié)果分析

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

[2]Yarus J M，Chambers R L.隨機(jī)建模和地質(zhì)統(tǒng)計學(xué)—原理、方法和實例研究[M].穆龍新，陳亮，譯.北京：石油工業(yè)出版社，2000.

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

[4] 劉振峰，郝天珧，楊長春.基于鏈模型的儲層巖相隨機(jī)模擬[J].地球物理學(xué)進(jìn)展，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ī)建模[J].地球物理學(xué)進(jìn)展，2003，18（3）：519-523.

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

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

[9] 劉振峰．基于MRF模型的油氣儲層屬性隨機(jī)建模方法研究[D].北京：中國科學(xué)院研究生院，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]沈洪濤，郭乃川.地質(zhì)統(tǒng)計學(xué)反演技術(shù)在超薄儲層預(yù)測中的應(yīng)用[J].地球物理學(xué)進(jìn)展，2017，32（1）：248-253.

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

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

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