地下水封洞库围岩质量的空间变异性分析与三维表征

杨 凯，晏鄂川，陈 武

(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%，可以推测洞室开挖时围岩的稳定性较好。