基于蒙特卡洛模擬計(jì)算連通率不同方法比較研究

蘭志廣 徐 源 聶細(xì)江 郭 猛 王海宇

（天津市勘察設(shè)計(jì)院集團(tuán)有限公司，天津 300191）

0 引言

結(jié)構(gòu)面連通率是綜合描述巖體內(nèi)部結(jié)構(gòu)面發(fā)育尺寸及貫通程度的重要指標(biāo)，在評(píng)價(jià)巖體的強(qiáng)度及穩(wěn)定性等方面都發(fā)揮著積極作用[1]。因此，國(guó)內(nèi)外學(xué)者針對(duì)連通率的計(jì)算方法進(jìn)行了許多有益探討與研究?？傮w來(lái)看，根據(jù)計(jì)算結(jié)構(gòu)面連通率所需數(shù)據(jù)來(lái)源（實(shí)測(cè)或者模擬）不同，連通率計(jì)算方法可分為兩大類，即基于實(shí)測(cè)的概率模型估算法和基于蒙特卡洛模擬計(jì)算法[2]。如黃建安[3]從結(jié)構(gòu)面組合的幾何關(guān)系出發(fā)研究計(jì)算了連通率；黃潤(rùn)秋等[4]通過(guò)結(jié)構(gòu)面全跡長(zhǎng)調(diào)查估算了連通率；黃國(guó)明等[5]在實(shí)測(cè)結(jié)構(gòu)面平均跡長(zhǎng)以及巖橋平均長(zhǎng)度的基礎(chǔ)上推導(dǎo)了H-H 連通率估算公式；汪小剛等[6]基于蒙特卡洛模擬考慮結(jié)構(gòu)面與巖橋相互作用的破壞機(jī)理并運(yùn)用動(dòng)態(tài)規(guī)劃原理計(jì)算了連通率；陳劍平等[7]基于蒙特卡洛模擬運(yùn)用理想截面投影法計(jì)算了連通率；李啟源[8]基于H-H 連通率估算公式提出了EHH 法估算了連通率；Zhang 等[9]基于蒙特卡洛模擬運(yùn)用最短路徑法計(jì)算了連通率。

基于實(shí)測(cè)的概率模型估算法估算連通率時(shí)，需要追蹤調(diào)查結(jié)構(gòu)面沿走向方向的延伸情況或者調(diào)查統(tǒng)計(jì)特定尺寸窗口范圍內(nèi)所有結(jié)構(gòu)面。由于現(xiàn)場(chǎng)結(jié)構(gòu)面往往隱伏于地表之下，即使出露其規(guī)模尺寸也有限，因此，現(xiàn)場(chǎng)結(jié)構(gòu)面調(diào)查難度較大，耗時(shí)耗力。相較而言，基于蒙特卡洛模擬計(jì)算法，僅依靠有限露頭面上調(diào)查所得結(jié)構(gòu)面，通過(guò)蒙特卡洛模擬生成結(jié)構(gòu)面網(wǎng)絡(luò)模型即可計(jì)算獲得連通率。因此，考慮到工程實(shí)踐的經(jīng)濟(jì)性，基于蒙特卡洛模擬計(jì)算法計(jì)算連通率更有優(yōu)勢(shì)。本文具體介紹了三種基于蒙特卡洛模擬計(jì)算連通率的方法，并通過(guò)舉例對(duì)比分析了不同方法結(jié)果的差異性，為以后工程實(shí)踐提供理論指導(dǎo)。

1 方法介紹

基于蒙特卡洛模擬計(jì)算法計(jì)算連通率的基本思路是，調(diào)查有限露頭面上結(jié)構(gòu)面數(shù)據(jù)并利用蒙特卡洛模擬生成結(jié)構(gòu)面網(wǎng)絡(luò)，然后以模擬所得結(jié)構(gòu)面網(wǎng)絡(luò)為數(shù)據(jù)依據(jù)計(jì)算連通率。由于結(jié)構(gòu)面網(wǎng)絡(luò)模擬可分為二維結(jié)構(gòu)面網(wǎng)絡(luò)模擬和三維結(jié)構(gòu)面網(wǎng)絡(luò)模擬，相應(yīng)的，基于蒙特卡洛模擬計(jì)算法可分為二維連通率計(jì)算方法和三維連通率計(jì)算方法。本文主要針對(duì)二維連通率計(jì)算方法進(jìn)行比較研究。目前，基于蒙特卡洛模擬計(jì)算二維連通率的方法可歸納為三大類，分別是帶寬投影法、動(dòng)態(tài)規(guī)劃法以及最短路徑法，介紹如下：

1.1 帶寬投影法

帶寬投影法是一種比較傳統(tǒng)的計(jì)算連通率的方法，其基本思路是：根據(jù)工程實(shí)際在結(jié)構(gòu)面網(wǎng)絡(luò)中布置特定方向的測(cè)量基線并設(shè)定帶寬，然后向該測(cè)量基線投影設(shè)定帶寬范圍內(nèi)所有的結(jié)構(gòu)面跡線（見(jiàn)圖1），連通率則可通過(guò)計(jì)算測(cè)量基線上投影所得結(jié)構(gòu)面跡長(zhǎng)的總和與測(cè)量基線長(zhǎng)度的比值獲得[10]，計(jì)算公式為：

圖1 帶寬投影法簡(jiǎn)圖

式中：L為測(cè)量基線的長(zhǎng)度；D為帶寬；n為結(jié)構(gòu)面跡線的條數(shù)；θ 為結(jié)構(gòu)面跡線與測(cè)量基線的夾角。

帶寬投影法的優(yōu)點(diǎn)是計(jì)算簡(jiǎn)單快捷，但其值的可靠度嚴(yán)重依賴于帶寬的取值。文獻(xiàn)[11]研究表明，隨著帶寬增加，采用帶寬投影法計(jì)算所得連通率值會(huì)相應(yīng)增加直至100%。因此，合理選擇帶寬是保證帶寬投影法可靠度的重中之重，而目前為止，帶寬的取值尚沒(méi)有統(tǒng)一定論，多依賴于在實(shí)際調(diào)查基礎(chǔ)上進(jìn)行經(jīng)驗(yàn)分析取值[2]。

1.2 動(dòng)態(tài)規(guī)劃法

動(dòng)態(tài)規(guī)劃法將連通率由幾何意義拓展到力學(xué)意義，其基本思路是：考慮結(jié)構(gòu)面與巖橋相互作用的破壞機(jī)理，并運(yùn)用動(dòng)態(tài)規(guī)劃原理在結(jié)構(gòu)面網(wǎng)絡(luò)中搜索沿特定方向且具有最小抗剪強(qiáng)度的結(jié)構(gòu)面-巖橋組合而成的潛在破壞路徑，這條潛在破壞路徑上結(jié)構(gòu)面跡長(zhǎng)總和與潛在破壞路徑長(zhǎng)度的比值即為連通率[6]，即：

圖2 結(jié)構(gòu)面與巖橋相互作用破壞機(jī)理

動(dòng)態(tài)規(guī)劃法考慮了巖體破壞時(shí)結(jié)構(gòu)面與巖橋之間相互作用的破壞機(jī)制，所得連通率值更具有理論價(jià)值和實(shí)際意義。文獻(xiàn)[11]研究表明，動(dòng)態(tài)規(guī)劃法所得連通率值與帶寬、模擬剖面尺寸以及正應(yīng)力有關(guān)：在一定范圍內(nèi)，連通率值會(huì)隨帶寬增大而增大，當(dāng)帶寬增大到一定程度后，連通率會(huì)趨于穩(wěn)定；模擬剖面尺寸對(duì)連通率值的影響規(guī)律同帶寬一致；正應(yīng)力取值對(duì)連通率值有一定影響，且在低應(yīng)力條件下，這種影響更明顯。

1.3 最短路徑法

最短路徑法與動(dòng)態(tài)規(guī)劃法思想類似，也是以結(jié)構(gòu)面網(wǎng)絡(luò)上搜索沿特定方向的潛在破壞路徑為目的，然后以潛在破壞路徑上結(jié)構(gòu)面跡長(zhǎng)總和與潛在破壞路徑長(zhǎng)度的比值作為連通率，即：

區(qū)別于動(dòng)態(tài)規(guī)劃法，最短路徑法搜索潛在破壞路徑時(shí)并沒(méi)有考慮結(jié)構(gòu)面與巖橋相互作用的破壞機(jī)理。但其合理性在于：由于結(jié)構(gòu)面的力學(xué)強(qiáng)度較完整巖石小幾個(gè)數(shù)量級(jí)，在巖體破壞時(shí)，往往有沿結(jié)構(gòu)面滑動(dòng)的趨勢(shì)。而巖體破壞路徑經(jīng)過(guò)的結(jié)構(gòu)面越多，其抗滑力就越小。假設(shè)完整巖石之間兩點(diǎn)的距離為d，相同距離的結(jié)構(gòu)面跡線中的兩點(diǎn)距離為d′，可視d′ ＜＜d。這樣，在限定潛在破壞路徑的入口與出口時(shí)，潛在破壞路徑的總長(zhǎng)度越小，經(jīng)過(guò)的結(jié)構(gòu)面跡長(zhǎng)就越大，潛在破壞路徑最短時(shí)，則最大可能的經(jīng)過(guò)了巖體中的結(jié)構(gòu)面，沿此路徑破壞時(shí)所受的抗剪力最小[13]。因此，巖體潛在破壞路徑可視為沿入口與出口之間的最短路徑。

目前，最短路徑法常用到的算法有Dijkstra 算法以及Floyd 算法。本文采用Floyd 算法，算法介紹見(jiàn)文獻(xiàn)[14]。以圖3所示節(jié)點(diǎn)為例，應(yīng)用Floyd 算法計(jì)算節(jié)點(diǎn)a 與節(jié)點(diǎn)d 之間的最短路徑，通過(guò)更新矩陣D和矩陣P，可得到a→d 之間的最短路徑為a→c→d，最短距離為3+1=4。

圖3 Floyd 算法計(jì)算最短路徑示例

最短路徑法計(jì)算思想簡(jiǎn)單，但值得注意的是，最短路徑算法針對(duì)節(jié)點(diǎn)進(jìn)行搜索，故在實(shí)際搜索時(shí)，需將結(jié)構(gòu)面網(wǎng)絡(luò)中的結(jié)構(gòu)面跡線進(jìn)行離散化。文獻(xiàn)[9,13,14]研究表明，結(jié)構(gòu)面跡線離散化程度越高，搜索結(jié)果越精確，但計(jì)算量也越大；此外，對(duì)于指定不同搜索方向，最短路徑法計(jì)算所得連通率值也不相同，體現(xiàn)了巖體結(jié)構(gòu)的各向異性。

2 舉例計(jì)算

2.1 舉例模型

某一水平巖體表面隨機(jī)發(fā)育有大量結(jié)構(gòu)面，統(tǒng)計(jì)分析其幾何參數(shù)，采用蒙特卡洛模擬生成一個(gè)尺寸為100 m (X軸) × 100 m (Y軸) × 100 m (Z軸)的三維結(jié)構(gòu)面網(wǎng)絡(luò)模型，關(guān)鍵步驟如下[1,14,15]：

（a）結(jié)構(gòu)面產(chǎn)狀模擬：采用基于施密特下半球等面積投影網(wǎng)的概率統(tǒng)計(jì)法對(duì)各組結(jié)構(gòu)面產(chǎn)狀頻率進(jìn)行校正[16]，校正公式見(jiàn)式（4）：

式中：Nα為結(jié)構(gòu)面產(chǎn)狀頻率的校正值；α為結(jié)構(gòu)面的平均傾角；Nα

′為觀測(cè)頻率。

（b）結(jié)構(gòu)面尺寸模擬：基于結(jié)構(gòu)面形狀為圓盤(pán)的假設(shè)，結(jié)構(gòu)面尺寸可由結(jié)構(gòu)面圓盤(pán)直徑來(lái)表征，且結(jié)構(gòu)面圓盤(pán)直徑與跡長(zhǎng)存在空間理論關(guān)系[16]：

式中：ch為出露半跡長(zhǎng)；r為圓盤(pán)半徑。

（c）結(jié)構(gòu)面密度模擬：采用張量法[17]來(lái)計(jì)算用于三維結(jié)構(gòu)面網(wǎng)絡(luò)模擬的裂隙密度。具體計(jì)算公式見(jiàn)式（6）：

2)為結(jié)構(gòu)面直徑平方的均值。

（d）蒙特卡洛模擬：通過(guò)上述步驟確定的結(jié)構(gòu)面幾何參數(shù)（見(jiàn)表1），利用MATLAB 自編代碼采用蒙特卡洛模擬法，即可生成一個(gè)尺寸為100 m (X軸) ×100 m (Y軸) × 100 m (Z軸)的三維結(jié)構(gòu)面網(wǎng)絡(luò)模型（見(jiàn)圖4）。

表1 結(jié)構(gòu)面主要幾何參數(shù)統(tǒng)計(jì)

圖4 三維結(jié)構(gòu)面網(wǎng)絡(luò)模型

選擇某剖面與三維結(jié)構(gòu)面網(wǎng)絡(luò)模型相截切，即可得到該剖面上的二維結(jié)構(gòu)面網(wǎng)絡(luò)[1]。本文以走向?yàn)?0°，且過(guò)三維結(jié)構(gòu)面網(wǎng)絡(luò)模型體心的剖面為例，模擬所得結(jié)構(gòu)面網(wǎng)絡(luò)見(jiàn)圖5。

圖5 二維結(jié)構(gòu)面網(wǎng)絡(luò)模擬圖

2.2 帶寬投影法計(jì)算

采用帶寬投影法計(jì)算連通率時(shí)，本文考慮影響計(jì)算結(jié)果的因素有投影基線方向和帶寬。在與二維結(jié)構(gòu)面網(wǎng)絡(luò)XY軸夾角為10°～180°的方向上每隔10°布置一條測(cè)量基線（見(jiàn)圖6），分別計(jì)算其帶寬為1～15 m（間隔1 m）時(shí)的連通率，計(jì)算結(jié)果見(jiàn)圖7。

圖6 測(cè)量基線布置圖

由圖7可以看出，不同投影基線方向計(jì)算所得連通率值是不同的，以投影帶寬1 m 為例，投影基線10°時(shí)連通率值最小，為19.97%；投影基線100°時(shí)連通率值最大，為46.09%。對(duì)比不同投影帶寬時(shí)連通率計(jì)算結(jié)果，雖然不同投影基線方向計(jì)算所得連通率值沒(méi)有明顯變化規(guī)律，但整體而言，均是在投影基線100°附近連通率值達(dá)到最大。此外，投影帶寬對(duì)連通率值的影響較為明顯，變化規(guī)律較為明確，隨著投影帶寬的增大，連通率值也相應(yīng)增加，以投影基線10°為例，隨著投影帶寬由1 m 增大到15 m，連通率值也由19.97%增大到83.34%。而在投影基線100°、投影帶寬15 m 時(shí)，連通率值最大，為96.53%。

圖7 帶寬投影法計(jì)算結(jié)果統(tǒng)計(jì)圖

2.3 動(dòng)態(tài)規(guī)劃法計(jì)算

采用動(dòng)態(tài)規(guī)劃法計(jì)算連通率時(shí)，本文考慮影響計(jì)算結(jié)果的因素有剪切方向和帶寬，正應(yīng)力對(duì)連通率值的影響暫不做考慮。假定巖體摩擦系數(shù)為0.95、黏聚力為0.85 MPa，結(jié)構(gòu)面摩擦系數(shù)為0.5、黏聚力為0.15 MPa，在與二維結(jié)構(gòu)面網(wǎng)絡(luò)XY軸夾角為10°～180°的方向上每隔10°布置一個(gè)剪切方向（見(jiàn)圖6），分別計(jì)算其帶寬為1～15 m（間隔1 m）時(shí)，正應(yīng)力為3 MPa 時(shí)的連通率，計(jì)算結(jié)果見(jiàn)圖8。

圖8 動(dòng)態(tài)規(guī)劃法計(jì)算結(jié)果統(tǒng)計(jì)圖

由圖8可以看出，不同剪切方向計(jì)算所得連通率值也是不同的，以帶寬1m 為例，剪切方向10°時(shí)連通率值最小，為29.61%；剪切方向100°時(shí)連通率值最大，為51.45%。對(duì)比不同帶寬時(shí)連通率計(jì)算結(jié)果，不同剪切方向計(jì)算所得連通率值除在剪切方向100°附近連通率值達(dá)到最大外，也沒(méi)有其他明顯變化規(guī)律。此外，帶寬對(duì)連通率值的影響較為明顯，變化規(guī)律較為明確，在一定范圍內(nèi)，隨著帶寬的增大，連通率值相應(yīng)增加，但當(dāng)帶寬增大到一定值后，連通率值趨于穩(wěn)定，以剪切方向10°為例，在帶寬1～9 m，連通率值由29.61%增大到64.41%，在帶寬9 m 之后，連通率值在65.59%與66.49%之間波動(dòng)。

2.4 最短路徑法計(jì)算

采用最短路徑法計(jì)算連通率時(shí)，本文考慮影響計(jì)算結(jié)果的因素有路徑搜索方向和帶寬，結(jié)構(gòu)面跡線離散化程度暫不做考慮。在與二維結(jié)構(gòu)面網(wǎng)絡(luò)XY軸夾角為10°～180°的方向上每隔10°布置一個(gè)搜索路徑方向（見(jiàn)圖6），分別計(jì)算其帶寬為1～15 m（間隔1 m）時(shí)，跡線離散為6 個(gè)點(diǎn)時(shí)的連通率。計(jì)算結(jié)果見(jiàn)圖9。

由圖9可以看出，不同搜索方向計(jì)算所得連通率值也是不同的，以帶寬1 m 為例，搜索方向10°時(shí)連通率值最小，為30.27%；搜索方向100°時(shí)連通率值最大，為50.77%。對(duì)比不同帶寬時(shí)連通率計(jì)算結(jié)果，不同搜索方向計(jì)算所得連通率值除在搜索方向100°附近連通率值達(dá)到最大外，也沒(méi)有其他明顯變化規(guī)律。此外，同動(dòng)態(tài)規(guī)劃法一樣，帶寬對(duì)連通率值的影響較為明顯，變化規(guī)律較為明確，在一定范圍內(nèi)，隨著帶寬的增大，連通率值相應(yīng)增加，但當(dāng)帶寬增大到一定值后，連通率值趨于穩(wěn)定，以搜索方向10°為例，在帶寬1～9 m，連通率值由30.27%增大到64.04%，在帶寬9 m 之后，連通率值在65.22%與66.70%波動(dòng)。

圖9 最短路徑法計(jì)算結(jié)果統(tǒng)計(jì)圖

3 計(jì)算結(jié)果對(duì)比分析

分析對(duì)比上述三種計(jì)算方法所得連通率結(jié)果，可以發(fā)現(xiàn)三種計(jì)算方法結(jié)果均受限于投影基線（或剪切、或搜索）方向和帶寬的選擇。其中，投影基線（或剪切、或搜索）方向?qū)B通率值的影響規(guī)律尚不明確，但從三種方法計(jì)算所得連通率值的結(jié)果來(lái)看，不同投影基線（或剪切、或搜索）方向計(jì)算所得連通率值不同，說(shuō)明三種方法計(jì)算所得連通率值均是呈現(xiàn)各向異性的；在不同帶寬時(shí)，三種方法均是在投影基線（或剪切、或搜索）方向100°時(shí)連通率值達(dá)到最大，說(shuō)明這三種方法在尋找連通率值最大（即巖體最容易發(fā)生破壞）方向是一致的。

對(duì)比帶寬對(duì)三種計(jì)算方法所得連通率值的結(jié)果（見(jiàn)圖10），可以明顯看出，在不同投影基線（或剪切、或搜索）方向時(shí)，帶寬投影法計(jì)算所得連通率值會(huì)隨著帶寬的增大持續(xù)增大，而動(dòng)態(tài)規(guī)劃法和最短路徑法計(jì)算所得連通率值均是在帶寬增大達(dá)到一定程度后趨于穩(wěn)定；此外，動(dòng)態(tài)規(guī)劃法和最短路徑法計(jì)算所得連通率值在不同帶寬時(shí)均較為一致，而帶寬投影法計(jì)算所得連通率值只有在帶寬2～3 m 范圍內(nèi)與其他兩種方法保持一致。

圖10 不同方法計(jì)算結(jié)果對(duì)比圖（以基線/剪切/搜索方向50°、100°、160°為例）

綜上分析，認(rèn)為在工程實(shí)踐中需求一個(gè)經(jīng)濟(jì)安全的連通率值時(shí)，帶寬投影法由于其不能提供一個(gè)穩(wěn)定的值，顯然是不合適的；相較而言，動(dòng)態(tài)規(guī)劃法和最短路徑法能滿足要求。但是比較動(dòng)態(tài)規(guī)劃法和最短路徑法，動(dòng)態(tài)規(guī)劃法需要明確巖體和結(jié)構(gòu)面強(qiáng)度參數(shù)及其所受正應(yīng)力，在初始條件不足時(shí)，動(dòng)態(tài)規(guī)劃法應(yīng)用也將受限，因此，本文推薦使用最短路徑法。在工程實(shí)踐中，建議采用最短路徑法首先粗略計(jì)算不同搜索方向連通率值，待確定最大連通率值所在方向后，重點(diǎn)計(jì)算該方向不同帶寬時(shí)的連通率值，并以趨于穩(wěn)定后的不同帶寬連通率平均值作為連通率建議值。如本文案例在帶寬5 m 之后連通率值于78.38%～81.15%波動(dòng)，建議取帶寬5～15 m 連通率均值，即79.75%作為連通率建議值。

4 結(jié)論

本文通過(guò)舉例分析對(duì)比帶寬投影法、動(dòng)態(tài)規(guī)劃法以及最短路徑法計(jì)算所得連通率結(jié)果，得到如下結(jié)論：

（1）就尋找最大連通率值所在方向，帶寬投影法、動(dòng)態(tài)規(guī)劃法及最短路徑法結(jié)果較為一致，且三種方法計(jì)算所得連通率值均表現(xiàn)為各向異性。

（2）投影基線（或剪切、或搜索）方向一定，帶寬投影法計(jì)算所得連通率值會(huì)隨著帶寬的增大持續(xù)增大，而動(dòng)態(tài)規(guī)劃法和最短路徑法計(jì)算所得連通率值均是在帶寬增大達(dá)到一定程度后趨于穩(wěn)定；且動(dòng)態(tài)規(guī)劃法和最短路徑法計(jì)算所得連通率值在不同帶寬時(shí)均較為一致。

（3）綜合分析考慮，建議在工程實(shí)踐中采用最短路徑法，首先粗略計(jì)算不同搜索方向連通率值，待確定最大連通率值所在方向后，重點(diǎn)計(jì)算該方向不同帶寬時(shí)的連通率值，并以趨于穩(wěn)定后的不同帶寬連通率平均值作為連通率建議值。