地下水封洞庫圍巖質(zhì)量的空間變異性分析與三維表征

楊 凱，晏鄂川，陳 武

(1.中國地質(zhì)大學(xué)(武漢)工程學(xué)院，湖北 武漢 430074；2.上海市政工程設(shè)計研究總院(集團)有限公司，上海 200092)

地下水封洞庫是一種利用“地下水封存”的方式將油(氣)品儲存在地下水水位以下適當(dāng)深度的堅硬完整巖石洞室中的石油儲備形式。因其具有安全、環(huán)保、建造成本低等特點，已成為石油儲備的一種重要形式[1]。巖體質(zhì)量為工程穩(wěn)定性分析、巖體合理利用以及巖體等效力學(xué)參數(shù)的選取等提供了可靠的依據(jù)，是地下水封洞庫設(shè)計、施工中需重點關(guān)注的因素。而受到沉積、構(gòu)造、風(fēng)化等長期地質(zhì)作用的影響，作為綜合反映巖體地質(zhì)力學(xué)性質(zhì)的指標(biāo)——巖體質(zhì)量不可避免地具有空間變異性。巖體質(zhì)量的空間變異性表現(xiàn)在空間分布的不均勻性、自相關(guān)性和各向異性，以及由于實踐認(rèn)知的有限導(dǎo)致的不確定性。目前在洞庫設(shè)計中普遍采用巖體質(zhì)量粗略的統(tǒng)計值來評價庫區(qū)巖體條件和估算巖體等效力學(xué)參數(shù)，忽略了巖體質(zhì)量的空間變異性，事實上由于普遍采用的高邊墻、大跨度的洞室形式，巖體質(zhì)量的空間變異性將使得相關(guān)分析和結(jié)論存在較大的誤差。為了更準(zhǔn)確地預(yù)判圍巖在洞室開挖及運行過程中的穩(wěn)定性狀況，有必要開展地下水封洞庫圍巖質(zhì)量的空間變異性研究。

相對于巖體強度參數(shù)、變形參數(shù)和水力學(xué)參數(shù)等空間變異性及其表征方法的研究[2-8]，對綜合反映巖體地質(zhì)力學(xué)性質(zhì)的指標(biāo)——巖體質(zhì)量空間變異性的研究較少。Stavropoulou等[9]在分析某隧道巖體質(zhì)量RMR分類空間變異性的基礎(chǔ)上，采用普通克里格的方法對隧道圍巖RMR值進行了插值計算，得到巖體質(zhì)量三維非均質(zhì)場；賈明濤等[10]、馮武等[11]均借助區(qū)域化變量理論研究了巖體質(zhì)量RMR分類指標(biāo)RQD、節(jié)理間距、摩擦角和抗壓強度的空間變異性，并采用克里格插值技術(shù)建立了反映全局變化的礦巖質(zhì)量RMR分類三維評價模型；Esfahani等[12]研究了巖體質(zhì)量RQD的空間變異性，通過序貫高斯模擬的方法得到RQD值的空間分布模型；Egana等[13]采用序貫高斯模擬的地質(zhì)統(tǒng)計學(xué)方法表征了巖體質(zhì)量RMR值的不確定性；Ferrari等[14]分別在二維和三維的基礎(chǔ)上分析了巖體質(zhì)量RMR值的變差函數(shù)，并采用普通克里格和條件模擬的方法估計了未采樣的RMR值；Mayer等[15]分別對某露天礦山邊坡場區(qū)地質(zhì)強度指標(biāo)和單軸抗壓強度指標(biāo)的空間變異性進行了分析，并采用序貫高斯模擬的方法表征了邊坡巖體質(zhì)量的空間變異性。

序貫高斯模擬(Sequential Gaussian Simulation,SGS)是地質(zhì)統(tǒng)計學(xué)中的一種條件模擬方法，由于其能還原實測點數(shù)據(jù)并給出區(qū)域化變量的多種可能實現(xiàn)，能較好地表征地質(zhì)介質(zhì)固有的空間變異性，該方法已在巖體強度參數(shù)、巖體滲透系數(shù)等空間變異性的模擬中得到了廣泛應(yīng)用，并取得了可喜的研究成果。本文依托某地下水封洞庫項目，以圍巖質(zhì)量的空間變異性為研究對象，在采用基于區(qū)域化變量理論對場區(qū)巖體質(zhì)量Q值的空間變異性特征進行分析的基礎(chǔ)上，使用序貫高斯模擬來表征巖體質(zhì)量Q值的空間變異性，并構(gòu)建巖體質(zhì)量Q值的三維隨機場,通過對多次SGS隨機實現(xiàn)結(jié)果的統(tǒng)計分析，以概率的手段預(yù)測了研究區(qū)巖體質(zhì)量的空間分布狀況。

1 研究區(qū)概況

擬建的寧波LPG地下水封洞庫位于寧波大榭島東北部，洞庫區(qū)范圍見圖1中“本文模擬區(qū)域”。該洞庫區(qū)處于低山丘陵地貌，地層主要為上侏羅統(tǒng)西山頭組含角礫熔結(jié)凝灰?guī)r和高塢組花崗斑巖，巖體中局部穿插有煌斑巖脈、安山巖脈和輝綠巖脈。根據(jù)鉆探和音頻大地電磁法顯示，場區(qū)發(fā)育有F1斷層以及P1、P2破碎帶(見圖1)，走向分別為NE、NNE和NE方向；區(qū)域性節(jié)理廣泛發(fā)育，地表節(jié)理分為兩組，走向分別為60°～90°和340°～10°，地下節(jié)理以緩傾角為主，優(yōu)勢方向不明顯，相對優(yōu)勢走向為60°～110°和320°～10°。洞庫區(qū)地下水類型主要為第四系孔隙潛水和基巖構(gòu)造裂隙水，其中基巖構(gòu)造裂隙水主要受大氣降水補給，排泄途徑受構(gòu)造發(fā)育特征影響。寧波地區(qū)經(jīng)歷了復(fù)雜的區(qū)域地質(zhì)構(gòu)造過程，造成了區(qū)域內(nèi)錯綜復(fù)雜的地質(zhì)構(gòu)造形跡。從場區(qū)節(jié)理、斷層和破碎帶以及區(qū)域地質(zhì)構(gòu)造演化來看，場區(qū)巖體在成巖和演化過程中主要受到華夏系北東向斷裂以及新華夏系北北東向斷裂的影響，形成的斷層、破碎帶以及區(qū)域性節(jié)理共同造成了巖體性質(zhì)在空間不同位置上的差異以及在不同方向上連續(xù)性的差異，從而導(dǎo)致了場區(qū)巖體質(zhì)量固有的空間變異性。

擬建洞庫規(guī)模為200萬m3，主洞室軸向方向為NW285°，東西長約1 000 m，南北長約800 m，由2組平行的地下洞罐組成：一組V2401號洞罐容量為120×104m3，由4條平行的主洞室(V2401-1～V2401-4)組成；另一組V2301號洞罐容量為80×104m3，由3條平行的主洞室(V2301-1～V2301-3)組成。7條主洞室洞頂?shù)脑O(shè)計標(biāo)高均為-125 m，截面形狀呈馬蹄形(圓形拱頂、圓形邊墻、水平底板)，最大寬度為20 m，高度為26 m，凈斷面積為430 m2；7條主洞室的長度由南至北依次為576 m、655m、660 m、707 m、750 m、720 m和690 m，北面V2401號洞罐和南面V2301號洞罐凈間距為95 m，各洞罐內(nèi)主洞室凈間距為48.5 m。該地下水封洞庫結(jié)構(gòu)見圖2。

圖1 研究區(qū)平面概況圖Fig.1 Plane profile of the study area

圖2 某地下水封洞庫結(jié)構(gòu)示意圖Fig.2 Schematic diagram of the structure of an underground water-sealed cavern

2 研究區(qū)巖體質(zhì)量Q值的空間變異性分析

2.1 分析方法

基于區(qū)域化變量理論的地質(zhì)統(tǒng)計學(xué)是研究巖體質(zhì)量空間變異性的重要理論。從地質(zhì)學(xué)的角度，區(qū)域化變量可以反映地質(zhì)變量的以下特征：①局部性。區(qū)域化變量只限于一定的空間內(nèi)；②連續(xù)性。不同的區(qū)域化變量具有不同的連續(xù)性；③異向性。區(qū)域化變量在各個方向上如果性質(zhì)相同，則稱為各向同性；若在各個方向上性質(zhì)不同，則稱為各向異性；④可遷性。區(qū)域化變量在一定范圍內(nèi)具有明顯的空間相關(guān)性，但超過這一范圍，相關(guān)關(guān)系變得很弱，甚至消失。

變差函數(shù)能較好地反映區(qū)域化變量的這些特殊性質(zhì)。假設(shè)空間點x只在一維的x軸上變化，把區(qū)域化變量Z(x)在x、x+h兩點處的值之差的方差的一半定義為Z(x)在x軸方向上的變差函數(shù)，記為γ(x,h)，即：

(1)

當(dāng)γ(x,h)與x的取值無關(guān)時，γ(x,h)只依賴于h(滯后、間隔、步長)，則可將γ(x,h)寫成γ(h)，此時以h為橫坐標(biāo)、γ(h)為縱坐標(biāo)作出的圖形叫變差函數(shù)圖，見圖3。圖中：a為變程，反映了變量的影響范圍，在影響范圍內(nèi)，變量的自相關(guān)性通常表現(xiàn)為隨著距離的增大而減弱，當(dāng)兩點間距離超過這一范圍，空間相關(guān)性消失;C0為塊金值，指的是變差函數(shù)在原點處的間斷性，反映區(qū)域化變量即使在很短的距離內(nèi)也存在較大的差異性;C為基臺值，指的是超出變程后變差函數(shù)的穩(wěn)定值，其大小反映了區(qū)域化變量的變化幅度大小;C0/C值反映了變量連續(xù)性的好壞，其值越小，表明連續(xù)性越好。

圖3 變差函數(shù)示意圖Fig.3 Diagram of variation function

在實際應(yīng)用中，變差函數(shù)往往是通過觀察值來作出估計的。根據(jù)數(shù)理統(tǒng)計知識，可以通過求[Z(x)-Z(x+h)]2和[Z(x)-Z(x+h)]的平均數(shù)來估計其數(shù)學(xué)期望，進而得到變差函數(shù)，因此必須要有Z(x)和Z(x+h)的若干次觀測值，而實際上空間同一點處只能得到一個觀測值，由此必須對區(qū)域化變量Z(x)提出一些假設(shè)。

(1) 平穩(wěn)假設(shè)。嚴(yán)格的平穩(wěn)假設(shè)指的是區(qū)域化變量Z(x)的任意n維分布函數(shù)不因空間點x發(fā)生位移h而改變。這種要求是Z(x)的各階矩存在，且平穩(wěn)，這在實際中不能滿足，且不好驗證，故提出了二階平穩(wěn)假設(shè)。即區(qū)域化變量Z(x)滿足以下兩個條件時，則稱Z(x)滿足二階平穩(wěn)(弱平穩(wěn))假設(shè)：

①整個研究區(qū)內(nèi)，Z(x)的數(shù)學(xué)期望存在，且等于常數(shù)：

E[Z(x)]=const,?x

(2)

②整個研究區(qū)內(nèi)，Z(x)的協(xié)方差函數(shù)存在且平穩(wěn)(即只依賴于滯后h，而與x無關(guān))：

Cov[Z(x),Z(x+h)]=E[Z(x)Z(x+h)]-E[Z(x)]E[Z(x+h)]=E[Z(x)Z(x+h)]-m2=C(h),?x,?h

(3)

上式中：E[Z(x)]=m。

特殊地,當(dāng)h=0時，

Var[Z(x)]=C(0)

(4)

當(dāng)上述條件仍不能滿足時，條件進一步放寬，導(dǎo)致本征假設(shè)。

(2) 本征假設(shè)。區(qū)域化變量Z(x)的增量[Z(x)-Z(x+h)]滿足下列兩個條件：

①在整個研究區(qū)內(nèi)有：

E[Z(x)-Z(x+h)]=0,?x,?h

(5)

②增量[Z(x)-Z(x+h)]的方差函數(shù)存在且平穩(wěn)(不依賴于x)，即：

Var[Z(x)-Z(x+h)]=E[Z(x)-Z(x+h)]2-{E[Z(x)-Z(x+h)]}2=E[Z(x)-Z(x+h)]2，?x,?h

(6)

在研究區(qū)域化變量的空間變異性時，通常做平穩(wěn)假設(shè)或本征假設(shè)，此時有E[Z(x)-Z(x+h)]=0，?h,且與x無關(guān)，則變差函數(shù)公式變?yōu)?/p>

(7)

實驗變差函數(shù)就是根據(jù)觀測數(shù)據(jù)構(gòu)造變差函數(shù)γ(h)的估計值γ*(h)。在平穩(wěn)假設(shè)或本征假設(shè)下，可根據(jù)空間上距離向量為h的點對xi和xi+h上的觀測值[Z(xi),Z(xi+h)]，其中i=1,2,…,N(h),用求[(Z(xi)-Z(xi+h)]2的算術(shù)平均值的方法來計算γ*(h)，有：

(8)

式中:N(h)為數(shù)據(jù)對[Z(xi),Z(xi+h)]的對數(shù)。

為了對區(qū)域化變量的未知值作出估計，還需要將實驗變差函數(shù)擬合成相應(yīng)的理論變差函數(shù)模型，這些模型將直接參與后續(xù)估值或條件模擬中。

2.2 基本統(tǒng)計分析

研究區(qū)內(nèi)共完成16個鉆孔，其平面位置見圖1，鉆孔總長為3 476.94 m，其中有5個豎直孔、11個斜孔，斜孔偏斜角度在30°以內(nèi)。本文選擇標(biāo)高在-230～-45 m之間的鉆孔巖芯數(shù)據(jù)，巖芯總長為2 297 m，最小標(biāo)高為-225 m，以每米為一個樣品，共2 297個樣品，得到研究區(qū)樣本巖體質(zhì)量Q值直方圖，見圖4。在統(tǒng)計的研究區(qū)2 297個樣品巖體質(zhì)量Q值中，Q值最大值為200，最小值為0.01，均值為38.23，中值為23.4，方差為1 414.67。

圖4 研究區(qū)樣本巖體質(zhì)量Q值直方圖Fig.4 Histogram of Q-value of rock mass of the samples of the study area

由圖4可見，研究區(qū)樣品巖體質(zhì)量Q值不符合高斯分布，由于Q值的計算為乘法型，且多數(shù)評價指標(biāo)的取值為離散值而無連續(xù)性，所以造成樣品巖體質(zhì)量Q值在分布上存在多處集塊。

巖體質(zhì)量Q值的非高斯分布反映了Q值的比例效應(yīng)或異方差性[16]。比例效應(yīng)從一種角度來看，指的是研究區(qū)內(nèi)各子域計算的局部均值與局部方差之間存在的大尺度關(guān)系；從另一角度來看，指的是局部不確定分布的均值與變異性存在的相關(guān)性。由于變差函數(shù)的計算依賴于區(qū)域化變量的量值，所以在比例效應(yīng)下，很難得到準(zhǔn)確的變差函數(shù)圖，因此在研究區(qū)域化變量結(jié)構(gòu)性前必須消除其比例效應(yīng)。本文采用高斯變換的方法消除區(qū)域化變量Q值的比例效應(yīng)，將非高斯分布空間變量轉(zhuǎn)換成高斯分布空間變量。

2.3 正態(tài)得分變換

為契合序貫高斯模擬，采用正態(tài)得分變換的方法將非高斯分布Q值變換為標(biāo)準(zhǔn)高斯分布。設(shè)Q值變換后為QNor，則QNor的計算公式可表示為

(9)

式中：FQ(Q)為Q值非高斯分布的分布函數(shù)；FN(QNor)為Q值高斯分布的分布函數(shù)。

考慮在已從小到大排列好的總量為L的樣本Qi(i=1，2，…，L)中，有n個樣本Q值相同，即：Qi=Qi+1=Qi+2=…=Qi+n-1(i+n≤L+1)，此時并不采用FQ(Qi+n-1)的值來等同地賦給這n個樣本的累積分布函數(shù)值，而是采用i/L，(i+1)/L，…，(i+n-1)/L對這n個樣本的累積分布函數(shù)值進行隨機賦值，這個方法相當(dāng)于對每一個Q值附加了一個很小的噪聲數(shù)據(jù)，因此結(jié)果會產(chǎn)生相同Q值而存在微小差異的QNor值。本文將研究區(qū)源數(shù)據(jù)變換為均值和方差分別為0和1的標(biāo)準(zhǔn)高斯分布,見圖5。

圖5 QNor值直方圖Fig.5 Histogram of QNor value

2.4 變差函數(shù)分析

本文QNor值的實驗變差函數(shù)計算是在大型三維數(shù)字化礦山軟件Surpac中進行的，該軟件通過設(shè)置搜索窗口(見圖6)來選擇計算半方差的點對。在計算實驗變差函數(shù)時，首先要指定一個滯后距Δh，則計算的滯后區(qū)間為：(0,Δh],(Δh,2Δh],(2Δh,3Δh]，再選擇計算方向(方位角和傾伏角)，并設(shè)定角度允許誤差和允許誤差限制，以確定搜索窗口；然后在每一個樣本點處通過搜索窗口確定分別落在各個滯后區(qū)間內(nèi)的點，并將這些點對分別添加到這些滯后區(qū)間的點對堆中；最后在每一個滯后區(qū)間內(nèi)分別計算落在該區(qū)間內(nèi)點對的半方差值的平均值以及點對間距的平均值，據(jù)此可得到實驗變差函數(shù)圖。

圖6 Surpac軟件設(shè)置的搜索窗口Fig.6 Search window set in the Surpac software

由于斜孔孔斜不大，且研究區(qū)內(nèi)鉆孔密度較小，故本文不考慮鉆孔的偏斜，將16個鉆孔共2 297個樣本正態(tài)得分變換后得到的QNor值輸入Surpac軟件中計算各個方向的變差函數(shù)。

從以往的研究成果來看，巖體質(zhì)量變異性的主軸方向多為水平方向和豎直方向，也有不少研究對巖體質(zhì)量進行平面或空間估值或模擬時，直接考慮其主方向在水平方向和豎直方向上，這些研究也表明巖體質(zhì)量在豎直方向上的變異性遠(yuǎn)比水平方向上的變異性大得多，因此各個方向上變差函數(shù)的計算可簡化為分別計算豎直方向上的變差函數(shù)和水平面上多個方向上的變差函數(shù)。

為了分析滯后距對所計算的實驗變差函數(shù)的影響，本文對滯后距分別為4 m、8 m、16 m和32 m時研究區(qū)樣品QNor值在豎直方向上的實驗變差函數(shù)圖進行了比較分析，見圖7。

圖7 不同滯后距下研究區(qū)樣本QNor值在豎直方向上的實驗變差函數(shù)對比Fig.7 Comparison of the experimental variogram of QNor values of the samples of the study area in the vertical direction under different lag distances

由圖7可見，從總體上來看，不同滯后距下QNor值在豎直方向上的實驗變差函數(shù)圖均存在一段趨勢段和平穩(wěn)段(需要說明的是，當(dāng)間距增大時，用于計算半方差平均值的有效點對減少，使得到的半方差均值的可靠性變差)，且各實驗變差函數(shù)圖之間并無大的差異，表明區(qū)域化變量在豎直方向上具有良好的結(jié)構(gòu)性；而從局部來看，在滯后距從4 m變化到32 m的過程中，研究區(qū)樣本QNor值的各實驗變差函數(shù)圖在小間距段的結(jié)構(gòu)特征發(fā)生了改變，在曲線變得平滑的同時，塊金值逐漸變大，可見在大滯后距下的實驗變差函數(shù)圖喪失了區(qū)域化變量的局部特征，而在小滯后距下的實驗變差函數(shù)圖中，這種特征能很好地得以顯現(xiàn)。

在水平面上，選擇滯后距為40 m，角度允許誤差和允許誤差限制取20°和80 m，分別計算了研究區(qū)樣品QNor值在水平面308°、353°、38°、83°四個方向上的實驗變差函數(shù)圖，見圖8。

圖8 研究區(qū)樣本QNor值在水平面各方向上的實驗變差函數(shù)對比Fig.8 Comparison of the experimental variogram of QNor values of the samples of the study area in each direction on the horizontal plane

由圖8可見，計算得到的研究區(qū)樣本QNor值在水平方向上的實驗變差函數(shù)圖顯得較為曲折，相比豎直方向上，QNor值的結(jié)構(gòu)性稍差，這一方面是鉆孔數(shù)據(jù)平面上的稀疏不規(guī)則造成的，另一方面是區(qū)域上復(fù)雜的地質(zhì)構(gòu)造歷史影響了地質(zhì)屬性的結(jié)構(gòu)性；盡管得到的QNor值在水平面上4個方向的實驗變差函數(shù)圖質(zhì)量并不好，但從總體上仍然可以看到趨勢段和平穩(wěn)段，這也驗證了QNor值為一個區(qū)域化變量。

對比研究區(qū)樣品QNor值在豎直方向和水平面不同方向上的實驗變差函數(shù)圖中可以看出，樣品QNor值在豎直方向上的實驗變差函數(shù)圖趨勢段最短，即自相關(guān)性最差；而樣品QNor值在水平面上4個方向中計算方向為38°時的實驗變差函數(shù)圖的趨勢段最長，計算方向為308°時的實驗變差函數(shù)圖的趨勢段相對最短，這表明QNor值在水平面上38°方向的自相關(guān)性最好，這也說明了QNor值存在各向異性。

2.5 各向異性分析

各向異性按其性質(zhì)可分為幾何各向異性和帶狀各向異性，一般以幾何各向異性來研究地質(zhì)介質(zhì)的各向異性。幾何各向異性指的是變差函數(shù)在不同方向上具有相同的基臺值(設(shè)塊金值為0)和不同的變程，即在相同距離上不同方向的兩點間平均變異程度不同。巖土體地質(zhì)力學(xué)特性在地質(zhì)構(gòu)造、沉積條件等相關(guān)地質(zhì)因素的作用下，其空間變異性通常表現(xiàn)出各向異性的特征。

為了分析區(qū)域化變量QNor值的幾何各向異性，需要在同等基臺值下對QNor值在各方向上的實驗變差函數(shù)進行擬合。對于基臺值的取值一直是一個極具爭議的話題，不少學(xué)者認(rèn)為當(dāng)距離足夠大時，區(qū)域化變量自相關(guān)性將消失，因此基臺值應(yīng)該等于樣本方差；也有研究表明當(dāng)研究域比較小時，實驗變差函數(shù)圖的基臺值應(yīng)稍大于樣本方差，而研究域足夠大時(大于3倍變程)，可采用樣本方差來估計基臺值[17]。從計算得到的研究區(qū)樣本QNor值在各個方向上的實驗變差函數(shù)圖來看，其多在樣本方差附近趨于平穩(wěn)，因此本文以樣本方差作為變差函數(shù)的基臺值。

Surpac軟件的擬合變差函數(shù)功能為用戶提供了一個人機交互式界面，用戶可以通過反復(fù)改變變差函數(shù)模型及參數(shù)并通過直觀對比判斷，確定最優(yōu)的擬合方案。在多次試探后發(fā)現(xiàn)，當(dāng)塊金值為0.3、采用球狀模型時對變差函數(shù)的擬合效果最好。據(jù)此得到研究區(qū)樣本QNor值在水平面308°、353°、38°、83°四個方向上變程分別為340 m、390 m、530 m、370 m時，其各向異性橢圓的擬合結(jié)果，見圖9。

圖9 研究區(qū)樣本QNor值在水平面不同方向上各向異性橢圓的擬合結(jié)果Fig.9 Fitting of the anisotropy ellipse of QNor values of the samples of the study area in each direction on the horizontal plane

表1給出了研究區(qū)樣本QNor值在各向異性橢球體三個主軸方向上變差函數(shù)模型的參數(shù)。

由表1可知，研究區(qū)樣本QNor值在各向異性橢球體主軸上最大連續(xù)性方向為NE35°，其與區(qū)域上北北東、北東向斷裂體系有很好的對應(yīng)關(guān)系，印證了前述新華夏系北北東向斷裂體系和華夏系北東向斷裂體系主控了場區(qū)巖體地質(zhì)力學(xué)性質(zhì)的空間變異性。

表1 研究區(qū)樣本QNor值在各向異性橢球體三個主軸方向上的變差函數(shù)結(jié)構(gòu)Table 1 Variogram structure of QNor values of the study area of the samples in the three main axial directions of the anisotropy ellipse

3 研究區(qū)巖體質(zhì)量Q值的序貫高斯模擬

3.1 序貫高斯模擬方法

序貫高斯模擬(SGS)是在克里格方法的基礎(chǔ)上產(chǎn)生的，下面將介紹克里格方法和序貫高斯模擬算法。

3.1.1 克里格方法

對于滿足二階平穩(wěn)假設(shè)(或本征假設(shè))的區(qū)域化變量Z(x)，在空間上得到一組離散的樣本信息數(shù)據(jù)Z(xi)(i=1，2，…，n)，采用線性估計計算未采樣點Z(x0)的估計量Z*(x0)為

(10)

(11)

其中,Cov(xi,xj)=C(0)-γ(|xi-xj|)

(12)

(13)

3.1.2 序貫高斯模擬算法

高斯隨機函數(shù)是最經(jīng)典的隨機函數(shù)，可由均值向量和協(xié)方差矩陣完全描述，并且其所有的條件分布都是高斯型的。高斯隨機函數(shù)擁有非常特殊且重要的空間結(jié)構(gòu)與分布規(guī)律，其中值在空間中擁有最大的相關(guān)性，而極值間的相關(guān)性則逐漸變小，這也要求所研究的區(qū)域化變量需要具備這樣的特征，當(dāng)區(qū)域化變量不服從高斯分布時，可采用前述的正態(tài)得分變換方法予以解決。

由前文可知，通過克里格方法能得到待估點屬性值的均值和方差，在區(qū)域化變量服從高斯分布的條件下，即可用均值和方差來完全描述待估點屬性值的累積概率分布函數(shù)。但簡單地通過所構(gòu)建的累積概率分布函數(shù)對空間網(wǎng)格上各待估點屬性值進行隨機取值的這種方法，顯然沒能遵從空間點數(shù)據(jù)之間存在特定依賴性這一事實。因此，對于隨機模擬所有待估點屬性值的累積概率分布函數(shù)，實際上為一個多變量的聯(lián)合條件累積概率分布函數(shù)[18]：

Prob{Z(x)≤z,x∈Ω|n(x)}

(14)

式中：z為自變量;Ω表示空間上的離散點構(gòu)成的點集;n(x)表示空間上各離散點處的屬性值組成的一組數(shù)據(jù)。

SGS的原理是將這種相互依賴的多個變量形成的聯(lián)合條件累積概率分布函數(shù)分解為N個單變量的條件累積概率分布函數(shù)：

(15)

式中：M表示原始采樣點屬性值構(gòu)成的數(shù)據(jù)集；N表示所有待估值點構(gòu)成的隨機路徑的順序序列。

簡單來說，SGS就是將計算得到的待估點處的屬性值依次添加到原有采樣數(shù)據(jù)集(條件數(shù)據(jù)集)，從而得到新的條件數(shù)據(jù)集，對下一個待估點屬性值的累積概率分布函數(shù)是以新條件數(shù)據(jù)集作為條件來確定的。通過SGS，多變量的聯(lián)合模擬可以得到解決，而不斷增多的條件數(shù)據(jù)集造成的問題則可通過僅保留最接近(或最相關(guān))的幾個條件數(shù)據(jù)來解決。

SGS是將高斯概率理論與貫序模擬算法相結(jié)合，在進行模擬時，首先隨機生成一條能訪問所有未知點的隨機路徑，然后依次訪問每一個點，并在融入估值數(shù)據(jù)后的條件數(shù)據(jù)集中搜索與其臨近的若干個點，依據(jù)克里格方法計算的均值和方差構(gòu)建估值點處的累積概率分布函數(shù)，并從中抽取一個值作為其隨機實現(xiàn)值，訪問完這條隨機路徑即得到了一次SGS隨機實現(xiàn)，若要進行下一次實現(xiàn)，則只需重新選擇一條隨機路徑。

3.2 表征區(qū)域與模型構(gòu)建

主洞室是用于儲存丙烷的地下倉庫，是地下水封洞庫工程的核心部分。為了分析寧波地下封洞庫主洞室圍巖質(zhì)量的空間分布，本文采用ABAQUS有限元軟件建立了該洞庫區(qū)7條主洞室的SGS隨機模擬網(wǎng)格數(shù)值模型，見圖10。該數(shù)值模型以洞室軸向方向為Y軸，將主洞室的空間布置進行了一定的簡化，將各洞罐內(nèi)主洞室長度設(shè)為相同值，V2301洞罐取624 m，以主洞室南端點開始向北延伸，V2401洞罐取720 m，以主洞室北端點開始向南延伸；模型X軸方向以主洞室3倍洞跨(60 m)處為邊界，其長度為596 m，Y軸方向最大長度為832 m；Z軸方向以標(biāo)高-230 m和-45 m為邊界，長度為185 m，模擬區(qū)域范圍見圖1。

圖10 某地下水封洞庫主洞室的SGS隨機模擬網(wǎng)格數(shù)值模型圖Fig.10 Numerical mesh model of SGS of the main cavern of an underground water-sealed cavern

地質(zhì)統(tǒng)計學(xué)估值和模擬是建立在巖體表征單元體(REV)在數(shù)值離散化范圍內(nèi)存在這一事實上，巖體表征單元體的基本概念表明存在一個范圍，個體非均質(zhì)性和離散特征能被忽略，從而可采用均一化過程得到巖體表征單元體內(nèi)的等效連續(xù)性質(zhì)，采用連續(xù)介質(zhì)分析方法(如有限元法)解決巖體力學(xué)問題?，F(xiàn)有不少對巖體力學(xué)參數(shù)REV大小的研究，如周創(chuàng)兵等[19]推導(dǎo)了估算裂隙巖體REV的數(shù)學(xué)表達式；向文飛[20]、張婷婷[21]、安玉華等[22]分別采用有限元法、離散元法和基于三維裂隙網(wǎng)格體積節(jié)理數(shù)的方法給出了裂隙巖體力學(xué)參數(shù)REV的尺寸約為5 m×5 m～8 m×8 m(二維)。

考慮到裂隙巖體REV尺寸大小，本次將網(wǎng)格共劃分為145 817個節(jié)點、133 611個單元，單元平均尺寸為8.26 m(洞室周邊單元較小)。

3.3 Q場序貫高斯模擬

SGeMS為地質(zhì)統(tǒng)計學(xué)建模軟件，該軟件是斯坦福大學(xué)所開發(fā)的一套用于解決空間相關(guān)變量中所涉及問題的計算機軟件包，為地質(zhì)統(tǒng)計學(xué)從業(yè)者提供了一個用戶界面友好、交互式3D可視化和可供廣泛選擇的多種地質(zhì)統(tǒng)計學(xué)算法[18]。在建立用于隨機模擬的網(wǎng)格或點集時，SGeMS軟件不僅可采用其自帶的笛卡爾網(wǎng)格，還支持通過導(dǎo)入二進制格式的點坐標(biāo)文件來建立點集，這個點集可以是不規(guī)則的。

對于ABAQUS有限元軟件建立的網(wǎng)格模型，可導(dǎo)出每個單元的質(zhì)心坐標(biāo)點，并改寫為SGeMS支持格式后，導(dǎo)入SGeMS中，即生成模擬點集。根據(jù)研究區(qū)樣本數(shù)據(jù)以及巖體質(zhì)量Q值空間變異性的分析結(jié)果對模擬點集的巖體質(zhì)量Q值空間變異性進行了序貫高斯模擬。

3.3.1 序貫高斯模擬結(jié)果與分析

前已敘及，SGS進行待估點估值計算的基礎(chǔ)為克里格方程組，其是依據(jù)已知點與待估點的空間關(guān)系來建立的，當(dāng)已知點數(shù)據(jù)較大時，若在進行每個待估點的估值時均利用所有已知點建立克里格方程組，計算量將會極大，所以在計算實現(xiàn)時必須使用“滑動鄰域”的方法，其能限制計算每個待估點參數(shù)值所考慮的已知點點數(shù)，本文選取距離待估點距離最近的20個已知點來計算待估點的Q值。利用原始數(shù)據(jù)及前述獲得的各主方向上的變差函數(shù)通過SGeMS軟件計算得到網(wǎng)格各單元質(zhì)心Q值，以此代表各單元的Q值。每一次SGS均能得到一個不同的結(jié)果，圖11展示了SGS一次隨機實現(xiàn)的結(jié)果云圖(圍巖模型中坐標(biāo)軸單位均為m，下同)，并給出了X≤138.5 m和X≤389.5 m[見圖11(b)]、Y=300 m和Y=600 m[見圖11(a)]、Z=-120 m[見圖11(e)]和Z=-150 m[見圖11(f)]6個平面的Q場云圖，以更加全面地揭示Q值在空間上的分布及其變化特征。其中，X≤138.5 m和X≤389.5 m的平面Q場云圖[圖11(b)]分別表示移除洞室?guī)r體后X≤138.5 m和X≤389.5 m部分模型的右視圖，圖中黑線為洞庫頂?shù)變蓷l輪廓線，下同)。

由圖11可見，研究區(qū)巖體質(zhì)量Q場云圖呈星星點點分布，Q值多集中在4～40區(qū)間，根據(jù)巖體質(zhì)量Q分類與通用五級分類法之間的對照關(guān)系可知，研究區(qū)巖體多為Ⅱ級和Ⅲ級巖體，總體上巖體質(zhì)量較好；從三維圖[見圖11(c)、(d)]以及X≤138.5 m和X≤389.5 m[見圖11(b)]、Y=300 m和Y=600 m[見圖11(a)]平面等Q場云圖均可以看出，水平方向上Q值的連續(xù)性比豎直方向上好；從三維圖[見圖11(c)、(d)]以及Z=-120 m[見圖11(e)]和Z=-150 m[見圖11(f)]平面Q場云圖可以看出，Q值具有一定方向的叢聚現(xiàn)象，高質(zhì)量巖體(紅色區(qū)域)或低質(zhì)量巖體(藍色區(qū)域)均在某個方向上具有較強的連續(xù)性，而在與其正交的方向上連續(xù)性明顯較差，通過觀察估計其連續(xù)性較強的方向在-20°～-50°之間(實際為NE35°～NE65°)，與各向異性分析得到的NE35°水平方向較為接近，可見研究區(qū)巖體質(zhì)量Q值的各向異性得到了較好的還原，也能較為直觀地反映北北東向斷裂、北東向斷裂造成的場區(qū)巖體的各向異性特征；另外無論是三維視圖還是剖面圖均能看到Q值在較小的范圍內(nèi)能產(chǎn)生較大的突變，以上說明SGS很好地表征了研究區(qū)巖體質(zhì)量Q值的空間變異性。

圖11 研究區(qū)巖體質(zhì)量Q值的SGS一次隨機實現(xiàn)結(jié)果云圖Fig.11 Cloud map of the result of one random implementation of the rock mass Q-value of the study area by SGS

此外，為了與SGS隨機實現(xiàn)的結(jié)果進行比較，本文還采用克里格方法(簡單克里格，下文中均為簡單克里格)對網(wǎng)格質(zhì)心點進行了估值，該方法能提供基于原始數(shù)據(jù)下區(qū)域化變量的最優(yōu)線性無偏估計，其所依據(jù)的原始數(shù)據(jù)、變差函數(shù)以及估值所考慮的點數(shù)等均與SGS相同。需要指出的是，由于前述巖體質(zhì)量空間變異性分析中所求得的是正態(tài)得分QNor值的變差函數(shù)，因此采用克里格方法估值前需要把原始數(shù)據(jù)Q值轉(zhuǎn)換成正態(tài)得分QNor值，在估值后再變換成原始數(shù)據(jù)分布型的數(shù)值。

本文選取X≤321 m和Y=500 m兩個平面的Q場云圖對SGS和克里格方法的優(yōu)劣進行了比較，見圖12。

由圖12可見，克里格(Krige)估值所得到的Q

圖12 研究區(qū)巖體質(zhì)量Q值的克里格插值與SGS隨機實現(xiàn)結(jié)果比較圖Fig.12 Comparison between the results of Kriging interpolation and SGS random implementation of the rock mass Q-value of the study area

場云圖較為連續(xù)，呈水平線條狀分布，Q值大多集中在10～40區(qū)間內(nèi)，Q值最大值為100左右，遠(yuǎn)小于原始數(shù)據(jù)最大值200，而Q值最小值也遠(yuǎn)大于原始數(shù)據(jù)最小值0.01，這是由于克里格所采用的線性插值方法會造成估值結(jié)果不可避免地存在平滑效應(yīng)，具體表現(xiàn)為較小值被夸大而較大值被低估，因此，克里格估值不能反映出巖體質(zhì)量Q值在局部范圍內(nèi)的真實變化特征，不能較好地表征巖體質(zhì)量Q值的空間變異性；與克里格估值所得的水平線條狀分布Q場云圖不同，SGS的兩次隨機實現(xiàn)(Realization 1和Realization 2)結(jié)果云圖均呈現(xiàn)星星點點的分布特點，且其所得到的Q值最高值和最低值均與原始數(shù)據(jù)保持一致，因此SGS隨機模擬方法更加遵從原始數(shù)據(jù)，且能較好地反映巖體質(zhì)量的空間變異特性。另外，無論是從SGS隨機實現(xiàn)還是克里格估值的Q場云圖中均能看到，Q值在豎直方向上相比水平方向上具有更大的變異性。

3.3.2 結(jié)果驗證

由于在進行SGS隨機模擬時，將原始數(shù)據(jù)作為已知值輸入到未知點估值計算中，因此在已知點處的Q值毋庸置疑等于原始數(shù)據(jù)的Q值，即SGS完全再現(xiàn)了原始數(shù)據(jù)點的值。此外，本文對SGS的5次隨機實現(xiàn)(Realization 1～Realization 5)和克里格(Krige)估值結(jié)果做了基本的統(tǒng)計分析，并與原始數(shù)據(jù)進行了比較，得到的統(tǒng)計參數(shù)見表2。

此外，本文還對5次SGS隨機實現(xiàn)和克里格估值結(jié)果與原始數(shù)據(jù)直方圖進行了比較，并統(tǒng)計了各分布在相同累積分布概率時的分位數(shù)值，同時以原始數(shù)據(jù)的分位數(shù)值作為橫坐標(biāo)，以估值(克里格估值和SGS隨機實現(xiàn))得到的結(jié)果數(shù)據(jù)的同累積分布概率的分位數(shù)值作為縱坐標(biāo)，繪制了等累積分布概率的分位數(shù)比較圖(Q-Q圖)，見圖13。

表2 研究區(qū)巖體質(zhì)量Q值的5次SGS隨機實現(xiàn)和克里格估值結(jié)果與原始數(shù)據(jù)統(tǒng)計參數(shù)比較Table 2 Comparison between the statistical parameters of the evaluation results of five random implementation by SGS and Krige method and the original data of the rock mass Q-value of the study area

注：偏度是統(tǒng)計數(shù)據(jù)分布偏斜方向和程度的度量，其值大于0時為右偏態(tài)，小于0時為左偏態(tài)，絕對值越大非對稱程度越高；峰度表征概率密度分布曲線在平均值處峰值高低的特征數(shù)，其值越大曲線峰部越尖銳，正態(tài)分布曲線的峰度為3。

圖13 研究區(qū)巖體質(zhì)量Q值5次SGS隨機實現(xiàn)和克里格估值結(jié)果與原始數(shù)據(jù)直方圖對比Fig.13 Comparison between the histogram of the evalua- tion results of five random implementation by SGS and Krige method and the original data of the rock mass Q-value of the study area

由表2和圖13可以看出：

(1) SGS對于Q場的各次隨機實現(xiàn)的統(tǒng)計參數(shù)均與原始數(shù)據(jù)保持良好的一致性，僅有極小的差異，說明了SGS隨機模擬很好地還原了巖體質(zhì)量Q值的直方圖特征；細(xì)觀地比較SGS不同隨機實現(xiàn)與原始數(shù)據(jù)的統(tǒng)計參數(shù)可以看到，SGS隨機實現(xiàn)的3個分位數(shù)值(最小值、最大值和中值)與原始數(shù)據(jù)幾乎相等，而其均值和變異系數(shù)均略小于原始數(shù)據(jù)，且其偏度和峰度均略大于原始數(shù)據(jù)，從Q-Q圖中也可以看到相似現(xiàn)象，即其散點有偏向于原始數(shù)據(jù)一側(cè)的態(tài)勢，說明了SGS隨機實現(xiàn)得到的Q值總體上相比原始數(shù)據(jù)小，且密度分布略顯尖銳和偏斜。造成這種誤差的主要原因是原始數(shù)據(jù)的集塊特征和正態(tài)得分變換對其產(chǎn)生的微小噪聲，從前述正態(tài)得分變換原理來看，集塊處的相同值將變換為存在微小差異的不同正態(tài)得分值，而這些值均小于或等于無噪聲影響時原始數(shù)據(jù)的得分值，從而在反變換時不可避免地還原成小于原始數(shù)據(jù)的數(shù)值，此外還可以分析得到集塊程度的大小影響著SGS隨機實現(xiàn)結(jié)果偏離原始數(shù)據(jù)程度的結(jié)論，這里不再贅述。

(2) 相比SGS的各次隨機實現(xiàn)，克里格插值結(jié)果得到的統(tǒng)計參量、直方圖和Q-Q圖并不令人滿意，其統(tǒng)計分析得到的偏度、峰度均遠(yuǎn)大于原始數(shù)據(jù)，而最小值和最大值均不同程度地達不到原始數(shù)據(jù)相應(yīng)的數(shù)值，可見克里格的平滑效應(yīng)使得估值結(jié)果在某一區(qū)間內(nèi)過分地集中，遠(yuǎn)遠(yuǎn)偏離原始數(shù)據(jù)的直方圖，未能再現(xiàn)原始數(shù)據(jù)的統(tǒng)計特征。

對于Q場空間變異性表征效果的檢驗，除上述直方圖的檢驗外，對于變差函數(shù)的檢驗也是不可或缺的。為了能與前述原始數(shù)據(jù)的變差函數(shù)在同等情況下進行對比，需要將SGeMS中得到的Q值經(jīng)過正態(tài)得分變換后再輸入到Surpac軟件中計算Q值在豎直方向上和水平面各方向上的實驗變差函數(shù)。通過Surpac軟件分析得到Q場的第一次SGS隨機實現(xiàn)估值結(jié)果(Realization 1)在豎直方向和水平面各方向上的實驗變差函數(shù)，見圖14和圖15，作為對比，圖中還展現(xiàn)了原始數(shù)據(jù)和克里格估值結(jié)果的實驗變差函數(shù)。

圖14 研究區(qū)巖體質(zhì)量Q值在豎直方向上的實驗變差函數(shù)對比Fig.14 Comparison of experimental variogram of Q value in the vertical direction

圖15 研究區(qū)巖體質(zhì)量Q值在水平面各方向上的實驗變差函數(shù)對比圖Fig.15 Comparison of experimental variogram of the rock mass Q-value of the study area in all horizontal directions

由圖14和圖15可見，盡管SGS隨機實現(xiàn)的實驗變差函數(shù)曲線與原始數(shù)據(jù)的實驗變差函數(shù)圖存在一定的差異，但SGS隨機實現(xiàn)結(jié)果相比原始數(shù)據(jù)表現(xiàn)出相對較強的變異性，且在水平方向上的變異性增幅較豎直方向上的略大；從總體上來看，SGS隨機實現(xiàn)的實驗變差函數(shù)曲線在一定程度上再現(xiàn)了原始數(shù)據(jù)的實驗變差函數(shù)圖，且相比克里格估值結(jié)果的實驗變差函數(shù)圖，SGS隨機實現(xiàn)的實驗變差函數(shù)曲線對原始數(shù)據(jù)變差圖的還原度明顯更高。

綜合以上關(guān)于SGS隨機實現(xiàn)直方圖和實驗變差函數(shù)圖的分析以及其與克里格插值、原始數(shù)據(jù)的對比，可以看出本文對于巖體質(zhì)量Q值空間變異性的三維表征效果是較成功的。

4 研究區(qū)巖體質(zhì)量概率分析

在基于SGS隨機模擬得到巖體質(zhì)量Q場的多次隨機實現(xiàn)結(jié)果后(本文共得到巖體質(zhì)量Q場的100次不同的隨機實現(xiàn))，可以采用概率統(tǒng)計的方法評價空間上各單元的巖體質(zhì)量好壞。經(jīng)100次的SGS隨機實現(xiàn)后，空間上每一個單元均產(chǎn)生100個不同的隨機值，通過統(tǒng)計在每一個單元的巖體質(zhì)量Q值>10(對應(yīng)巖體質(zhì)量等級Ⅰ級和Ⅱ級)所發(fā)生的頻率，以此來代表其可能發(fā)生的概率P(Q值>10)(若需提高精度，可進行更多次的隨機實現(xiàn))，可以得到Q值>10的概率空間分布，圖16展示了Y=300 m、Y=400 m、Y=500 m和Y=600 m 4個剖面處的概率P(Q值>10)分布圖。

圖16 研究區(qū)不同剖面處的概率P(Q值>10)分布圖Fig.16 Distribution of the probability P(Q-value above 10) at different sections of the study area

由圖16可見，研究區(qū)絕大多數(shù)巖體有很高的可能性(P>0.75)表現(xiàn)為Ⅰ級或Ⅱ級巖體，而截面中所有巖體均有不小于0.55的可能性為Ⅰ級或Ⅱ級巖體，這也說明了研究區(qū)的巖體條件良好；比較不同剖面處的概率分布圖可知，剖面Y=500 m和Y=600 m相比剖面Y=300 m和Y=400 m巖體質(zhì)量為Ⅰ級或Ⅱ級的可能性更高，且高質(zhì)量巖體在Y=500 m和Y=600 m兩剖面西南部相比北東部出現(xiàn)的可能性更高，根據(jù)其空間位置關(guān)系，推測可能為破碎帶P2的影響，根據(jù)鉆孔資料和物探測試成果，破碎帶P2總體上為NE走向，其延展性由西南向北東增大，使得在Y=500 m剖面西南部和Y=600 m剖面北東部附近巖體質(zhì)量總體上相對偏低。

通過100次的SGS隨機實現(xiàn)，還可以統(tǒng)計得到各個單元的等分位數(shù)，進而得到研究區(qū)各個單元所有隨機實現(xiàn)構(gòu)成的直方圖的置信水平為20%的單側(cè)置信上限Q20%值云圖，即得到各個單元的巖體質(zhì)量Q值有80%的置信水平大于云圖中各單元所承載巖體質(zhì)量的數(shù)值，為了避免贅述，這里僅給出了研究區(qū)Y=300 m和Y=600 m兩個剖面處的Q20%值云圖，見圖17。

圖17 研究區(qū)Y=300 m和Y=600 m剖面處的Q20%值云圖Fig.17 Cloud maps of Q20% value at different sections of the study area

由圖17可見，研究區(qū)兩個剖面所有單元均有80%的置信水平為Ⅲ級及以上巖體，且Y=600 m剖面的巖體質(zhì)量總體上比Y=300 m剖面的巖體質(zhì)量好。

通過圖16和圖17，可以預(yù)測洞室在開挖過程中的相對不穩(wěn)定巖體，并指導(dǎo)安全施工，但其精度卻有賴于原始數(shù)據(jù)的好壞以及SGS隨機實現(xiàn)的次數(shù)。例如根據(jù)圖17可以預(yù)測，Y=300 m剖面處V2301號洞罐中間洞室的圍巖質(zhì)量保證率相對較低，其邊墻可能會發(fā)生失穩(wěn)，并有針對性地提出預(yù)防措施。

對于地下工程，對圍巖開挖后穩(wěn)定性影響最大的往往是洞室開挖線周邊最鄰近的巖體，淺層巖體的失穩(wěn)破壞將可能引發(fā)更大范圍的圍巖失穩(wěn)破壞，嚴(yán)重影響洞室的穩(wěn)定性，圖18(a)展示了V2301號和V2401號兩個洞罐各洞室淺層圍巖單元網(wǎng)格剖分橫截面，單元總數(shù)為8 598個。對SGS的各次隨機實現(xiàn)均可以統(tǒng)計在V2301號和V2401號洞罐圍巖質(zhì)量Q值小于或等于4的巖體單元個數(shù)N(Q值≤4)，由于Q值≤4時為Ⅳ級或Ⅴ級巖體，這類巖體在洞室開挖時會產(chǎn)生較大的變形，巖體失穩(wěn)的可能性較大，因此可認(rèn)為其為不穩(wěn)定單元，經(jīng)多次SGS隨機實現(xiàn)后可得到場區(qū)洞庫不穩(wěn)定單元個數(shù)N(Q值≤4)的直方圖，見圖18(b)。該直方圖表示多次SGS隨機實現(xiàn)中出現(xiàn)各個數(shù)量的不穩(wěn)定單元時的次數(shù)占比(頻數(shù))。

圖18 場區(qū)洞庫“淺層圍巖”及其不穩(wěn)定單元統(tǒng)計直方圖Fig.18 Statistical histogram of shallow surrounding rock of the cavern and its unstable units in the field area

由圖18(b)可見，場區(qū)洞庫不穩(wěn)定單元個數(shù)N(Q值≤4)大致呈正態(tài)分布，峰部在600～800之間，相對來說不穩(wěn)定單元較少的情況(<700)比較多的情況(>700)出現(xiàn)的可能性大，不穩(wěn)定單元平均約占8%，可以推測洞室開挖時巖體穩(wěn)定性較好。

5 結(jié) 論

(1) 某地下水封洞庫場區(qū)巖體質(zhì)量Q值的最大值為200，最小值為0.01，均值為38.23，其直方圖總體上呈正偏斜分布；變差函數(shù)在大滯后距下會喪失局部特征，建議以小滯后距來建立變差函數(shù)；QNor值在豎直方向上展現(xiàn)出良好的結(jié)構(gòu)性，而QNor值在水平面4個方向上的結(jié)構(gòu)性稍差，另一方面QNor值在水平面各個方向上的相關(guān)性存在差異。場區(qū)巖體質(zhì)量各向異性的分析表明，相關(guān)性最好的方向為NE35°水平方向，這與區(qū)域上北北東、北東向地質(zhì)構(gòu)造相對應(yīng)，印證了前述場區(qū)巖體質(zhì)量所表現(xiàn)出的空間變異性主要受控于新華夏構(gòu)造體系和華夏構(gòu)造體系的結(jié)論。

(2) 通過SGS隨機實現(xiàn)的平、縱、橫剖面以及三維視圖的展現(xiàn)，較為全面地揭示了研究區(qū)巖體質(zhì)量Q值的空間分布特征，并將其與克里格插值的方法得到的Q值分布圖進行了對比，結(jié)果表明：SGS隨機模擬較好地表征了研究區(qū)巖體質(zhì)量的空間自相關(guān)性、變異性、各向異性等特征，也直觀地反映了北北東向、北東向斷裂對場區(qū)巖體質(zhì)量的影響形跡。此外，通過SGS隨機實現(xiàn)結(jié)果、克里格插值結(jié)果與原始數(shù)據(jù)統(tǒng)計參量、直方圖、Q-Q圖以及變差函數(shù)圖的比較，說明SGS隨機模擬對場區(qū)巖體質(zhì)量空間變異性的表征效果是令人滿意的。

(3) 通過對100次SGS隨機實現(xiàn)結(jié)果的統(tǒng)計分析，以概率統(tǒng)計的手段預(yù)測了研究區(qū)巖體質(zhì)量的空間分布狀況。第一，通過統(tǒng)計模型中各單元Q值大于10(Ⅰ級和Ⅱ級巖體)所發(fā)生的頻率，以此代替其發(fā)生概率P(Q值>10)，結(jié)果表明研究區(qū)絕大多數(shù)巖體有大于75%的概率為Ⅰ級或Ⅱ級巖體；第二，通過統(tǒng)計得到各巖體單元Q值置信水平為20%的單側(cè)置信上限Q20%值，結(jié)果表明研究區(qū)所展示剖面處所有單元均有80%的置信水平為Ⅲ級及以上巖體；第三，通過統(tǒng)計場區(qū)7個洞室“淺層圍巖”的單元個數(shù)N(Q值≤4)所發(fā)生的頻率，可預(yù)測洞室開挖后不穩(wěn)定巖體的占比，結(jié)果表明不穩(wěn)定巖體平均約占8%，可以推測洞室開挖時圍巖的穩(wěn)定性較好。