裂隙網(wǎng)絡的三維隨機模擬方法研究

沈瑞文

(安徽理工大學 地球與環(huán)境學院，安徽 淮南 232001)

裂隙網(wǎng)絡的模擬研究在地質(zhì)學領(lǐng)域的諸多方面都有十分重要的應用，裂隙網(wǎng)絡的準確刻畫和模擬是地質(zhì)學領(lǐng)域具有挑戰(zhàn)性的問題之一[1]。利用隨機方法模擬裂隙網(wǎng)絡最早于20世紀80年代末由Neuman提出，該方法早期多應用于模等效連續(xù)介質(zhì)模型，90年代后期逐步應用于裂隙網(wǎng)絡的隨機模擬[2]。目前，裂隙網(wǎng)絡模擬技術(shù)已逐漸成熟，其成果已在工程上得到了廣泛的應用。陶凱等利用加權(quán)Voronoi圖法生成二維裂隙網(wǎng)絡模型，并應用在膨脹土的力學性質(zhì)和滲流特性[3]；王晉麗等利用Monte-Carlo隨機模擬技術(shù)生成二維裂隙網(wǎng)絡滲流模型，研究了裂隙網(wǎng)絡中水流在不同邊界條件下的流動特征[4]；敖雪菲等在研究壩基裂隙巖體灌漿時，基于Monte-Carlo方法建立了裂隙巖體三維網(wǎng)絡模型，結(jié)果表明該法能有效反應裂隙真實特征并獲得貼近實際的灌漿結(jié)果[5]；張弛等利用布爾模擬和多點統(tǒng)計學方法進行了裂隙網(wǎng)絡的三維模擬，并對比兩種模型的精度[6]；李瑞金等基于Monte-Carlo方法實現(xiàn)裂隙的三維隨機模擬，并應用于壩基多孔分序灌漿數(shù)值模擬[7]。

近年來，隨著裂隙網(wǎng)絡模擬研究在采礦設(shè)計以及地下工程施工中應用的增加，對裂隙網(wǎng)絡模擬的精度要求越來越高。以往根據(jù)裂隙的定量數(shù)據(jù)對裂隙進行二維模擬的方法僅僅對模擬對象進行平面展布，沒有考慮模擬對象方位角度的變化，與實際裂隙網(wǎng)絡之間的差異較大，難以達到高精度地質(zhì)模型的要求[8]。三維地質(zhì)隨機模擬技術(shù)相較于二維模擬具有諸多優(yōu)點，它在原二維模擬的基礎(chǔ)上考慮了裂隙的位置和方向(走向和傾角)，有利于對地質(zhì)變量各屬性間的關(guān)系進行充分研究，確定裂隙的變化趨勢，建立更加符合實際情況的裂隙網(wǎng)絡，為采礦設(shè)計、地質(zhì)資源預測、地下工程施工、地質(zhì)災害預報等領(lǐng)域的研究和工作提供科學的定量依據(jù)[9]。本文在裂隙變量數(shù)據(jù)的基礎(chǔ)上，采用Monte-Carlo和地質(zhì)統(tǒng)計學相結(jié)合的三維隨機模擬方法模擬出裂隙元，并同時考慮裂隙的距離和角度，按照一定的標準連接成裂隙面，形成高精度的、與實際裂隙吻合程度較高的裂隙網(wǎng)絡。

1 三維隨機模擬原理

近年來，在不同領(lǐng)域有許多方法用來模擬裂隙網(wǎng)絡的分布，常見的方法有加權(quán)Voronoi圖法、序貫指示模擬、布爾模擬、多點統(tǒng)計學模擬、拉丁超立方抽樣模擬(LHS)、Monte-Carlo等等[10]。Monte-Carlo模擬是通過一定的隨機數(shù)生成方法，生成服從某一隨機變量的概率分布形式的隨機數(shù)序列[11]。裂隙具有方向、位置、形狀等多重屬性，它在巖體中的分布是隨機的、復雜的[12]，要對裂隙網(wǎng)絡進行模擬具有一定的難度，但是大量實踐證明裂隙的定量數(shù)據(jù)(密度)往往服從一定的概率分布，因此，可以運用Monte-Carlo來研究裂隙分布特征以及模擬裂隙網(wǎng)絡分布[13]。地質(zhì)統(tǒng)計學的基本分析方法是克里金法，克里金是一種確定性內(nèi)插方法，通過它計算得到的實現(xiàn)期望偏差相對最小。克里金方法是建立在隨機地質(zhì)介質(zhì)模型基礎(chǔ)之上(即變差函數(shù)基礎(chǔ)上)的確定性求解，可以求出與未知的客觀實際值之間偏差最小的唯一解。將Monte-Carlo和地質(zhì)統(tǒng)計學結(jié)合在一起，流程如圖1所示，對裂隙網(wǎng)絡進行三維模擬研究，得出的結(jié)果與實際裂隙之間的吻合性較高，同時可以更好地為后續(xù)的采礦設(shè)計以及地下工程施工提供精度較高的依據(jù)[14]。

圖1 裂隙網(wǎng)絡模擬流程

2 裂隙網(wǎng)絡非定量數(shù)據(jù)的模擬方法

根據(jù)上文所述的原理，歸納出裂隙網(wǎng)絡模擬實現(xiàn)的步驟：

1)裂隙位置模擬。在裂隙網(wǎng)絡的三維隨機模擬過程中，裂隙密度的計算占有首要位置，它反映了裂隙在研究區(qū)域中發(fā)育的密集程度,密度是定量數(shù)據(jù)，以往通常采用線密度或面密度等二維形式進行計算,這里對裂隙密度用三維形式——體密度進行計算，體密度是指研究區(qū)域內(nèi)單位體積分布的裂隙條數(shù)[15]，其計算公式為

(1)

式中：n為裂隙總數(shù)目；ri為裂隙平均直徑；V為取樣體積。

對于研究區(qū)域內(nèi)的每一個空間單元格，可以計算落入其中的裂隙中心點的數(shù)目，從而根據(jù)單元的大小得到各個單元格的裂隙密度。

Monte-Carlo隨機模擬方法的基本原理是利用[0，1]區(qū)間標準均分隨機數(shù)，根據(jù)由已知裂隙幾何參數(shù)的密度分布函數(shù)求得的抽樣公式來獲得服從給定裂隙幾何參數(shù)分布形式的隨機變量的模擬模型，其抽樣方法中直接抽樣方法最為常用和有效，具體實施方案：設(shè)隨機變量x服從分布函數(shù)y(x)和累積分布函數(shù)Y(x)，Y(x)的值域為[0，1]，ti作為Y(xi)的函數(shù)值，則x和t之間的對應關(guān)系為[16]

(2)

其反函數(shù)為

x=Y-1(t).

(3)

通過計算機程序產(chǎn)生一系列均勻分布隨機數(shù)x1，x2，…，xn，將這些均勻分布隨機數(shù)代入式(2)，即得到隨機變量X的模擬模型[17]

Xi=Y-1(ti)，i=1,2,…,n.

(4)

式中：n為試驗次數(shù)，實踐表明，試驗次數(shù)越多，X的頻率分布越接近其真實的概率分布，在實際中n一般取頻率分布收斂時所對應的試驗次數(shù)，在三維隨機模擬研究過程中，當模擬的次數(shù)足夠多時，就可以得到相對比較精確的概率值[18]。

裂隙位置的模擬是根據(jù)已經(jīng)計算出的具有正態(tài)分布特征的密度值，運用Monte-Carlo和地質(zhì)統(tǒng)計學相結(jié)合的方法，在密度的半變異函數(shù)模型的基礎(chǔ)上，應用普通克里格對裂隙密度的空間展布進行模擬，并隨機模擬出單位體積內(nèi)的裂隙條數(shù)，最后得到裂隙位置空間分布圖。

2)裂隙方向模擬。走向和傾角在裂隙幾何特征中占有重要位置。裂隙在構(gòu)造應力的作用下會產(chǎn)生優(yōu)勢方向，根據(jù)走向玫瑰花圖、傾角分布的直方圖對其產(chǎn)狀分布特征進行分析，用n個互斥的方向組來表示裂隙的方向，裂隙方向用(t1t2…tn)形式表達。若裂隙方向在該方向組內(nèi)出現(xiàn)，則這一組的值為1，同時，其他組因裂隙方向不在這些組內(nèi)，則賦值為0[19]，從而地質(zhì)表示的裂隙方向即可以實現(xiàn)向指示形式表示的轉(zhuǎn)換。

主成分分析實質(zhì)上是一種數(shù)學變換的方法，它把一組相關(guān)變量通過線性變換轉(zhuǎn)換成另一組不相關(guān)的變量，將這些新的不相關(guān)的變量按照方差依次遞減的順序排列，就形成所謂的主成分，使得第一主成分具有最大的方差，第二主成分的方差次之，并且和第一主成分不相關(guān)，依次類推，形成p個主成分[20]。

運用主成分分析法對裂隙的指示值進行計算，得出裂隙方向組的主成分，并計算主成分的實驗變異函數(shù)，再用普通克里格法模擬出計算得到的主成分的空間分布，用n個0和1組合成的指示形式將預估點處的主成分值反演出來，其中，最大值對應的方向組作為該點裂隙方向所屬的組, 根據(jù)該方向組內(nèi)樣本裂隙走向和傾角的累積分布函數(shù), 用隨機函數(shù)產(chǎn)生該裂隙的方向[21]，如圖2所示。

圖2 裂隙方向模擬流程

3)裂隙面的建立。對已經(jīng)模擬出的裂隙元，以一定的距離和角度方向為標準，視空間分布位置、方向相同或相近的裂隙同屬于一個相同的裂隙面，將其以一定的準則連接為裂隙面[22]。

令裂隙中心點之間的距離為Lc，裂隙與其中心點連線的夾角為α1，α2，如圖3所示，設(shè)L為裂隙元連接的最大距離，一般規(guī)定為2～3個單元格大?。沪翞榱严对B接的最大角度，一般規(guī)定為5°～10°。將符合式(5)的裂隙連接成同一個裂隙面[23]

圖3 裂隙元連接

Lc

(5)

3 研究實例

由于朱仙莊礦區(qū)煤層主要位于二疊系地層內(nèi)部，本文進行裂隙模擬的地層即為朱仙莊煤系地層——二疊系。對朱仙莊礦區(qū)內(nèi)的實際裂隙數(shù)據(jù)進行搜集，對搜集到的樣本數(shù)據(jù)進行分析，以10 m為單位，統(tǒng)計單元格內(nèi)的裂隙數(shù)目，并以此計算出裂隙密度，如圖4所示,得出的數(shù)據(jù)符合正態(tài)分布，并根據(jù)實際裂隙的展布情況，做出走向玫瑰花圖，如圖5所示，可以看出裂隙主要分布在NW、NE等方向。

圖4 裂隙密度頻率

圖5 走向玫瑰花圖

根據(jù)裂隙密度，將Monte-Carlo模擬和地質(zhì)統(tǒng)計學方法相結(jié)合，在密度的半變異函數(shù)模型基礎(chǔ)上，應用普通克里格模擬出裂隙空間分布，用普通克里格對裂隙密度的空間展布進行模擬，并隨機模擬出單位體積內(nèi)的裂隙條數(shù)，最后得到裂隙位置空間分布圖，如圖6所示。

圖6 裂隙密度的半變異函數(shù)(A)、裂隙密度空間分布(B)和裂隙位置空間分布(C)

將裂隙方向在走向a∈[0°,360°]和傾角b∈[0°,360°]的范圍內(nèi)分為6個方向組：

T1{a∈[0°,60°),b∈(0°,90°]},

T2{a∈[60°,120°),b∈(0°,90°]},

T3{a∈[120°,180°),b∈(0°,90°]},

T4{a∈[180°,240°),b∈(0°,90°],

T5{a∈[240°,300°),b∈(0°,90°]},

T6{a∈[300°,360°),b∈(0°,90°]}.

對這6個方向組進行處理，將其用二值數(shù)字0和1組成的指示形式表示出來。利用主成分分析這6個指示變量，求出3個主成分(PC1，PC2，PC3)，對其進行計算并分別模擬出半變異函數(shù)模型，如圖7所示。基于半變異函數(shù)模型應用普通克里格模擬出PC1，PC2，PC3的空間分布，如圖8所示。按照6個0和1構(gòu)成的變量反算出待估裂隙的主成分值，則該裂隙方向即為最大值對應的方向組，根據(jù)模擬得到的裂隙方向組，應用Monte-Carlo和地質(zhì)統(tǒng)計學相結(jié)合，模擬出裂隙方向。

圖8 3個主成分對應的空間分布(Ⅰ、 Ⅱ、 Ⅲ)

圖7 3個主成分對應的半變異函數(shù)(Ⅰ、 Ⅱ、 Ⅲ)

根據(jù)三維隨機模擬得到的裂隙元，在綜合考慮其距離和角度的基礎(chǔ)上，將屬于同一裂隙面的裂隙元連接起來，形成裂隙面，如圖9所示?？梢钥闯觯芯繀^(qū)域內(nèi)的大裂隙面主要分布在中部、東南部，展布較為均勻，裂隙的發(fā)育方向主要在NW和NE。在應用Monte-Carlo進行模擬時，模擬次數(shù)越多結(jié)果越接近實際情況，而在實際模擬過程中，模擬次數(shù)不可能趨于無窮多，同時在收集實際裂隙數(shù)據(jù)時存在一定的誤差，因此，模擬結(jié)果與實際情況具有一定的偏差。綜合考慮多次模擬結(jié)果，得到與實際裂隙的走向趨勢和空間分布具有較好的吻合性的模型。在模擬出的裂隙三維網(wǎng)絡的基礎(chǔ)上，可以清晰地看出研究區(qū)域內(nèi)斷裂走向變化的趨勢，從而為進一步的采礦設(shè)計、地下工程施工等提供依據(jù)。

圖9 實際裂隙分布(A)和模擬的裂隙網(wǎng)絡空間分布(B)

4 結(jié)束語

運用Monte-Carlo和地質(zhì)統(tǒng)計學相結(jié)合的方法對裂隙網(wǎng)絡進行三維隨機模擬，在Monte-Carlo模擬方法的實現(xiàn)過程中，會產(chǎn)生多個等概率的反映裂隙網(wǎng)絡空間結(jié)構(gòu)的模擬結(jié)果，各個結(jié)果間的差別反映了模擬次數(shù)、累積概率以及資料不足等原因而引起的不確定性， 綜合考慮這些模擬結(jié)果，可以得到與實際裂隙展布情況較一致的裂隙網(wǎng)絡。在對朱仙莊含煤系地層的裂隙網(wǎng)絡模擬研究上，可以清楚地看出模擬出的裂隙與原實際裂隙在空間位置分布和方向變化趨勢上都具有較好的吻合性，進一步說明在運用Monte-Carlo模擬方法的基礎(chǔ)上，結(jié)合地質(zhì)統(tǒng)計學，綜合考慮裂隙的方向，有利于實現(xiàn)高精度裂隙網(wǎng)絡模型的模擬，對合理評價礦區(qū)內(nèi)裂隙結(jié)構(gòu)特征具有重要意義，同時也能夠為采礦設(shè)計和地下工程施工提供精度較高的依據(jù)。