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

楊 凱，晏鄂川，陳 武

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

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

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

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

1 研究區概況

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

擬建洞庫規模為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條主洞室洞頂的設計標高均為-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，各洞罐內主洞室凈間距為48.5 m。該地下水封洞庫結構見圖2。

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

圖2 某地下水封洞庫結構示意圖Fig.2 Schematic diagram of the structure of an underground water-sealed cavern

2 研究區巖體質量Q值的空間變異性分析

2.1 分析方法

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

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

(1)

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

圖3 變差函數示意圖Fig.3 Diagram of variation function

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

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

①整個研究區內，Z(x)的數學期望存在，且等于常數：

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

(2)

②整個研究區內，Z(x)的協方差函數存在且平穩(即只依賴于滯后h，而與x無關)：

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。

特殊地,當h=0時，

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

(4)

當上述條件仍不能滿足時，條件進一步放寬，導致本征假設。

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

①在整個研究區內有：

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

(5)

②增量[Z(x)-Z(x+h)]的方差函數存在且平穩(不依賴于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)

在研究區域化變量的空間變異性時，通常做平穩假設或本征假設，此時有E[Z(x)-Z(x+h)]=0，?h,且與x無關，則變差函數公式變為

(7)

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

(8)

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

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

2.2 基本統計分析

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

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

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

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

2.3 正態得分變換

為契合序貫高斯模擬，采用正態得分變換的方法將非高斯分布Q值變換為標準高斯分布。設Q值變換后為QNor，則QNor的計算公式可表示為

(9)

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

考慮在已從小到大排列好的總量為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個樣本的累積分布函數值，而是采用i/L，(i+1)/L，…，(i+n-1)/L對這n個樣本的累積分布函數值進行隨機賦值，這個方法相當于對每一個Q值附加了一個很小的噪聲數據，因此結果會產生相同Q值而存在微小差異的QNor值。本文將研究區源數據變換為均值和方差分別為0和1的標準高斯分布,見圖5。

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

2.4 變差函數分析

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

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

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

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

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

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

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

圖8 研究區樣本QNor值在水平面各方向上的實驗變差函數對比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可見，計算得到的研究區樣本QNor值在水平方向上的實驗變差函數圖顯得較為曲折，相比豎直方向上，QNor值的結構性稍差，這一方面是鉆孔數據平面上的稀疏不規則造成的，另一方面是區域上復雜的地質構造歷史影響了地質屬性的結構性；盡管得到的QNor值在水平面上4個方向的實驗變差函數圖質量并不好，但從總體上仍然可以看到趨勢段和平穩段，這也驗證了QNor值為一個區域化變量。

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

2.5 各向異性分析

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

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

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

圖9 研究區樣本QNor值在水平面不同方向上各向異性橢圓的擬合結果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給出了研究區樣本QNor值在各向異性橢球體三個主軸方向上變差函數模型的參數。

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

表1 研究區樣本QNor值在各向異性橢球體三個主軸方向上的變差函數結構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值的序貫高斯模擬

3.1 序貫高斯模擬方法

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

3.1.1 克里格方法

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

(10)

(11)

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

(12)

(13)

3.1.2 序貫高斯模擬算法

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

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

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

(14)

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

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

(15)

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

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

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

3.2 表征區域與模型構建

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

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

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

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

3.3 Q場序貫高斯模擬

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

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

3.3.1 序貫高斯模擬結果與分析

前已敘及，SGS進行待估點估值計算的基礎為克里格方程組，其是依據已知點與待估點的空間關系來建立的，當已知點數據較大時，若在進行每個待估點的估值時均利用所有已知點建立克里格方程組，計算量將會極大，所以在計算實現時必須使用“滑動鄰域”的方法，其能限制計算每個待估點參數值所考慮的已知點點數，本文選取距離待估點距離最近的20個已知點來計算待估點的Q值。利用原始數據及前述獲得的各主方向上的變差函數通過SGeMS軟件計算得到網格各單元質心Q值，以此代表各單元的Q值。每一次SGS均能得到一個不同的結果，圖11展示了SGS一次隨機實現的結果云圖(圍巖模型中坐標軸單位均為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)]分別表示移除洞室巖體后X≤138.5 m和X≤389.5 m部分模型的右視圖，圖中黑線為洞庫頂底兩條輪廓線，下同)。

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

圖11 研究區巖體質量Q值的SGS一次隨機實現結果云圖Fig.11 Cloud map of the result of one random implementation of the rock mass Q-value of the study area by SGS

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

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

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

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

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

3.3.2 結果驗證

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

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

表2 研究區巖體質量Q值的5次SGS隨機實現和克里格估值結果與原始數據統計參數比較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

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

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

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

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

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

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

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

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

4 研究區巖體質量概率分析

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

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

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

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

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

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

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

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

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

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

5 結 論

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

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

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