北山沙枣园花岗岩体随机结构面的空间分布

霍 亮，王贵宾，李亚伟，魏 翔，刘桓兑

1) 东华理工大学核资源与环境国家重点实验室，南昌 330013

2) 中国科学院武汉岩土力学研究所岩土力学与工程国家重点实验室，武汉430071

3) 核工业北京地质研究院中核集团高放废物地质处置评价技术重点实验室，北京 100029

4) 重庆工商大学，重庆 400067

国际公认的高放废物最安全可行的方法是深地质处置法，将高放废物埋藏在深部地质处置库[1]（通常指500～1000 m）中，对放射性核素形成包容、阻滞作用，使它与人类的生存环境永久隔离.由于高放射性核废物漫长的半衰期，深地质处置库的安全评价期限至少要达到1×104a[1]，对高放废物的隔离最终依赖围岩完成.我国的高放废物处置工程的重点围岩是花岗岩[2].花岗岩中广泛分布的随机结构面是核素迁移的天然通道，是控制岩体变形和渗透特性的关键因素.因此针对花岗岩体中随机结构面的空间分布开展研究，对于我国高放废物处置工程的长久有效密闭至关重要.

在随机结构面的数据获取方面，技术人员主要采取的测量手段有两种：露头节理调查和钻孔裂隙编录.对于北山花岗岩体中发育良好的天然露头节理，采用综合节理测线法[3]测量随机结构面的倾向、倾角、迹长和间距，采用圆形窗口法[3-4]计算露头节理的中点面密度、平均迹长，采用统计方法分析结构面几何特征满足的概率分布[5-7]，采用Miller 方法开展岩体均质区的划分[6]，采用基于粒子群的K-means 聚类算法[7]划分地表露头发育的节理优势组，并进一步开展节理三维网络模拟和渗透张量分析[8].

对于北山花岗岩体中钻孔裂隙的研究，运用超声波钻孔电视技术识别钻孔裂隙，采用WellCAD软件解译裂隙的深度、倾向和倾角值[9-10]，统计分析钻孔裂隙中线密度和产状随深度的变化以及隙宽分布[10]，并开展钻孔岩体均质区的划分和三维网络模拟[11-12].部分学者结合地表露头节理和钻孔裂隙的产状数据，采用Fracman 软件[13]和3DEC 软件[14]开展北山岩体随机结构面的三维网络模拟.

总结上述北山花岗岩体随机结构面的研究现状：由于地表露头和钻孔开挖面的限制，以及天然岩体经历构造运动的复杂性，当前尚未实现对岩体内全部结构面几何参数的系统性测量，进而无法开展岩体结构面整体的完备性描述工作.研究人员为了克服调查方法的局限性，依据区域岩体中部分随机结构面的测量数据，推测区域岩体中随机结构面整体的发育特征，从而拓展随机结构面的认识边界.当前研究主要存在两方面问题：（1）多数研究采用单一的调查手段开展天然岩体中随机结构面的研究与模拟，钻孔裂隙与露头节理的调查研究相互割裂：在认识岩体结构面整体时，通常假设地表与地下的随机结构面分布规律相同，或者地下的随机结构面具有水平分层发育的特征，存在以偏概全的误区；（2）研究中过于注重随机结构面的统计特征，忽略结构面的空间属性：岩体随机结构面整体并非纯粹随机生成，而是具有空间组织特征，仅基于概率统计分布的三维网络模拟忽视了随机结构面的空间复杂性，导致其网络模拟结果难以指导工程实践.

针对上述两个问题，本文选择我国高放废物处置北山沙枣园备选场址中尺寸为600 m×600 m×600 m 的岩体为研究区域，以其中的单一钻孔裂隙以及地表的露头节理为研究对象.首先开展结构面的产状统计分析，比较岩体钻孔裂隙与地表露头节理的产状分布，并采用聚类分析开展结构面优势组划分；其次开展随机结构面的位置分析，分别恢复钻孔裂隙和露头节理产状聚类分析的深度分布和平面分布，全面认识北山沙枣园目标花岗岩体中随机结构面的空间分布特征.

1 岩体随机结构面

北山沙枣园岩体位于我国甘肃省酒泉市金塔县西北部，如图1（a）所示，区域构造背景在文献[15]中有详细介绍，本文不再赘述；沙枣园岩体的海拔高度为1200～1500 m，岩性以花岗闪长岩为主，岩体地表风化严重，由浑圆山丘组成，坡度较小，山谷低平，相对高差几十米，砾石较普遍，呈戈壁荒漠地貌，如图1（b）所示；对沙枣园目标岩体的随机结构面调研以垂直钻孔为中心轴开展，示意如图1（c）所示；采用超声波钻孔电视开展600 米深钻孔的裂隙识别工作，获取68 条裂隙的深度和产状数据，裂隙编录的技术细节参考文献[9-10]，钻孔裂隙的三维重建如图1（d）所示；采用综合测线法开展岩体地表36 万m2的露头节理测量，技术细节参考文献[7]，钻孔位于地表中心位置，地表分布10 个露头，测得204 条节理的产状数据，将节理简化为线段，设置钻孔为坐标零点，开展节理的平面展布恢复，如图1（e）所示，将典型露头O36放大呈现，如图1（f）所示.

图1 沙枣园区域600 m×600 m×600 m 岩体随机结构面的调查示意图.(a)研究区域位置;(b)研究区域地貌;(c)研究区域随机结构面的调查方法;(d)钻孔裂隙建模;(e)区域露头节理平面展布;(f)O36 露头平面展布Fig.1 Investigation map of the stochastic discontinuities of 600 m × 600 m × 600 m rock mass in the Shazaoyuan: (a) location of the study area;(b) geomorphology of the study area;(c) survey methods for stochastic discontinuities;(d) 3D modeling of cracks in the borehole;(e) joint distribution on the surface of rock mass;(f) joint distribution of the O36 outcrop

图1（c）中钻孔岩体的岩体质量指标值(Rock quality designation，RQD)为94，图1（e）中岩体地表露头测量信息以及据测线计算的岩体RQD 值，如表1 所示.

表1 岩体地表露头的测量信息及RQD 值Table 1 Measuring information and RQD of outcrops

表1 中地表露头测量面积总和为1487 m2，占据区域岩体地表36 万m2调查面积的0.413%，测量面积有限.露头测线的RQD 值均在97 以上，钻孔的RQD 值为94，两者的RQD 值均大于90，显示该研究区域岩体完整，岩体质量“很好”.下面分别针对随机结构面的产状和位置开展详细分析.

2 随机结构面的产状分析

倾角和倾向是随机结构面的空间产出状态.首先绘制露头节理与钻孔裂隙倾角的频数柱状分布图；其次绘制结构面的极点等密度图；最后参考极点汇集区域数量，设置优势组的聚类数目，运用基于粒子群的K-means 算法[16]对露头节理与钻孔裂隙的产状分别进行聚类分析，并对该聚类算法获取的计算聚类中心进行优化，选择邻近该中心的真实结构面产状作为聚类中心，如图2 所示.

图2 随机结构面的产状统计图.(a)露头节理倾角柱状分布图;(b)露头节理极点等密度图;(c)露头节理聚类分析图;(d)钻孔裂隙倾角柱状分布图;(e)钻孔裂隙极点等密度图;(f)钻孔裂隙聚类分析图Fig.2 Statistical analysis of the occurrences of stochastic discontinuities: (a) dip distribution of outcrop joints;(b) pole equal density map of outcrop joints;(c) cluster analysis of outcrop joints;(d) dip distribution of borehole cracks;(e) pole equal density map of borehole cracks;(f) cluster analysis of borehole cracks

对比图2（a）与（d），区间[65°，85°]的倾角数据分别占露头节理与钻孔裂隙总量的94%与58%；区间[0°，45°]的倾角数据分别占露头节理与钻孔裂隙总量的1.5%与11.5%.两种调查方法中倾角的集中区间相似，沙枣园岩体主要发育倾角大于65°的陡倾角结构面.

由于随机结构面的倾角越大，与近水平的露头面相交概率越大，与垂直钻孔相交的概率越小.图2（a）与（d）中，区间[75°，85°]的倾角数据分别占露头节理与钻孔裂隙的52%与4%，体现出不同调查方法对随机结构面采样的偏差.因此，对于主要发育陡倾角结构面的沙枣园岩体，露头调查的取样偏差相对较小.岩体随机结构面真实倾角的峰值区间在两种调查方法峰值之间，且偏向于露头节理.

观察图2（b）与（e），将其中走向相近，倾向相差180°的极点划分为一个汇集区域.图2（b）中露头节理的极点汇集在四个区域，依次为：NNE 向陡倾角节理、NE 向陡倾角节理、NW 向陡倾角节理，最后一个极点汇集区域没有对称区域.图2（e）中钻孔裂隙的极点同样汇集在四个区域，依次发育：北北东（North-northeast，NNE）向陡倾角裂隙、近东西（Eastwest，EW）向陡倾角裂隙、北西（Northwest，NW）向陡倾角裂隙、北西西（Northwest by west，NWbW）向缓倾角裂隙.

根据图2（b）与（e）中极点汇集区域的数量，设置聚类数目为4，对比图2（c）、（f）与图2（b）、（e），聚类中心位置与极点汇集区域位置基本吻合，钻孔裂隙与露头节理优势组的产状及数据占比，如表2所示.

表2 随机结构面的优势产状分组Table 2 Dominant sets of stochastic discontinuities

表2 中钻孔中发育的EW 向（189.5°∠71.2°）和NW 向（235.6°∠69.2°）优势组与露头节理中发育的EW 向（180°∠77°）和NW 向（231°∠67°）优势组分别对应，走向分别相差9.5°和4.6°，倾角分别相差5.8°和2.2°，数据占比分别相差4.4%和9.2%.

钻孔中发育最显著的NNE 向优势组（295.4°∠68.3°），在露头节理的优势组中对应为两组优势产状：NNE 向（282°∠80°）和NE 向（320°∠79°）优势组，走向相差13.4°和24.6°，倾角相差11.7°和10.7°.观察图2(e)，两处优势产状的极点汇集区域相连，并在相连处凹陷，呈现倾角较小的特征，结合图2（a）与（d）中[75°，80°]区间倾角占比4%与40%的差异性，在该优势组的钻孔采样受到钻孔位置及采样窗口局限的影响，明显偏少.同样由于取样偏差，钻孔中观察到的NWbW 向缓倾角优势组66.0°∠33.0°，在露头节理的聚类分析中没有对应优势组.

因此，钻孔裂隙与露头节理中的产状分布以相似性为主，其中NW 向、EW 向结构面呈现明显的相似性，而NNE 向结构面则呈现一定差异性.

3 随机结构面的位置分析

结构面的产状决定了其空间形态，位置决定了结构面的空间存在.

3.1 露头节理的平面分析

基于调查露头的空间位置和数量，将美国气象学家Thiessen[17]提出的Voronoi 多边形法，拓展至区域岩体地表划分中.用测线起点坐标表征露头，连接各个露头点，形成Delaunay 三角形，绘制三角形各边的垂直平分线，若干平分线围成Voronoi多边形，将研究区域的岩体表面划分为10 个区块.

区块中包含唯一露头，采用露头名称指代区块，如图3（a）所示，以O36露头为例，采用相切圆的圆形窗口法[18]计算单一露头中节理的平均迹长值L和中点面密度值P20，将O36露头的属性值赋予O36区块，计算公式如式（1）所示，计算值如表3.

表3 岩体地表完整露头及拆分露头的平均迹长和中点面密度计算表Table 3 Average trace lengths and midpoint plane densities of complete outcrops and split outcrops of the rock mass surface

图3 完整露头节理的平面特征分布图.(a)岩体地表的Voronoi 剖分;(b)平均迹长分布;(c)中点面密度分布Fig.3 Plane characteristic distribution map of complete outcrops: (a) Voronoi dissection of the rock surface;(b) average trace length distribution map of outcrops;(c) midpoint plane density distribution map of outcrops

式中，c为圆形窗口半径；圆形窗口与窗口内取样节理有三类相互关系：贯穿型，相交型和包容型，数目依次为N0、N1、N2.

将上述方法推广至地表其他区块，采用完整露头的平均迹长和中点面密度值充填区块，直观呈现地表不同区块中随机结构面的平面特征差异，如图3（b）、（c）所示.

图3（b）、（c）中颜色深浅表征节理平均迹长值和中点面密度的相对大小.图3（b）东北部区块的平均迹长小于西南部，参照表3 中完整露头节理的参数值，东北部O18、O29区块的平均迹长值在5 m左右，而钻孔附近的O1、O19区块的平均迹长值在8.5 m 左右，西南部区块的平均迹长值在11 m 以上.

图3（c）中点面密度的平面分布规律不同于平均迹长.钻孔附近O1区块达到最大值0.12 m-2，区域北部和东南部的中点面密度较小，其中O18、O29区块的中点面密度分别为0.042 m-2和0.067 m-2，呈现中间大，周围小的分布规律.

由于岩体的各向异性，在经历构造运动过程中，地表不同区块露头中发育的节理组不同，单一优势组在不同区块的发育强度不同.以图2（c）中节理产状优势组划分为依据，以露头O36为例，保留图3（a）中露头节理的计算窗口，对单一露头进行节理拆分，如图4（a）所示；同理，放大呈现地表其他区块中的露头节理走向拆分图，如图4（b）所示.

图4 露头节理的拆分计算图.(a) 单一露头节理的优势组拆分;(b) 岩体地表露头节理的优势组拆分Fig.4 Split calculation diagram of outcrop joints: (a) split joints of O36;(b) spatial distribution of joint dominant sets

由图4（a），露头O36中仅有三组优势节理，参照式（1），可得拆分露头中不同优势组的平均迹长Li和中点面密度参数值，如式（2）：

式（2）中，i代表优势组的划分类型，c为相切圆的半径，据i类优势组节理与相切圆窗口的关系划分为：贯穿型，相交型与包容型，三种类型相加之和为Ni.

据Mardia[19]提出的方向统计理论，提出拆分露头中不同优势组节理的走向平均值 αi，如式（3）：

式（3）中，m表征单一露头中归属优势组i的节理数目，表征单一露头中归属优势组i的任意节理j的走向值，该值以优势组内节理相邻为取值范围.

由图4（b）可得，发育最显著的黄色NNE 向节理，除在O1区块没有发育外，其他区块均发育，不同区块发育程度不同，在位于钻孔东北与西南一侧的O18、O29、O2、O37、O43区块中是最优产状；发育次显著的蓝色NE 向节理，在O43区块没有发育，在其他区块中均有发育，在位于钻孔东南侧的O1、O11、O19、O38区块中是最优产状；橙色EW 向节理，在O43区块没有发育，在位于钻孔西北侧的O36区块是最优产状，在O38、O1、O19区块中是次优产状；紫色NW 向节理，区块中均非最优产状，在O43、O29、O11区块中为次优产状，在钻孔西侧的O36、O37、O38区块中没有发育.

图4（b）定性呈现了地表不同优势组的分布差异，下面开展定量分析.据式（1）计算完整露头的平面特征参数，以及区域节理的平均迹长值和中点面密度值平均值；据式（2）、（3）计算拆分露头中不同优势组的平面特征参数，并计算区域不同优势组的平均迹长值、平均中点面密度值和走向平均值，以及上述数值的标准差，其中方向数据的标准差，计算如式（4）所示，并采用离散系数衡量数据的均匀性，计算值如表3 所示.

式（5）中K为随机结构面的聚类数目.由式（2）可得，拆分露头中不同优势组的Li值和值与该优势组中节理的数目无关，而与圆形窗口中的不同类型节理的数目相关，并且值取决于相交型和包容型节理的数目相加.值得注意的是，如果圆形窗口内的节理类型为贯穿型节理，即均为零，则值为零，Li值无意义.因此图5（b）中露头O2发育1 条EW 向优势组节理，但在表（3）中该露头i=2(EW 向优势组)中并没有数值.

图5 拆分露头节理的平面分析矩阵图.(a) 1 组优势节理的平均走向图;(b) 2 组优势节理的平均走向图;(c)3 组优势节理的平均走向图;(d) 4 组优势节理的平均走向图;(e) 1 组优势节理的平均迹长分布图;(f) 2 组优势节理的平均迹长分布图;(g) 3 组优势节理的平均迹长分布图;(h) 4 组优势节理的平均迹长分布图;(i) 1 组优势节理的中点面密度分布图;(j) 2 组优势节理的中点面密度分布图;(k) 3 组优势节理的中点面密度分布图;(l) 4 组优势节理的中点面密度分布图Fig.5 Matrix of plane analysis of the split outcrops: (a) average strike distribution map of cluster 1 joints;(b) average strike distribution map of cluster 2 joints;(c) average strike distribution map of cluster 3 joints;(d) average strike distribution map of cluster 4 joints;(e) average trace length distribution map of cluster 1 joints;(f) average trace length distribution map of cluster 2 joints;(g) average trace length distribution map of cluster 3 joints;(h) average trace length distribution map of cluster 4 joints;(i) midpoint surface density distribution map of cluster 1 joints;(j) midpoint surface density distribution map of cluster 2 joints;(k) midpoint surface density distribution map of cluster 3 joints;(l) midpoint surface density distribution map of cluster 4 joints

表3 拆分露头节理中，相对于优势组i=4，优势组i=1、2、3 时的区域走向平均值的标准差均在15°以下，说明区域发育的NE 向、EW 向与NNE向优势组结构面走向更集中.

表3 中相对于优势组i=2、4，优势组i=1、3 时，区域拆分露头节理的迹长平均值均在8 m 以上，中点面密度值均在0.024 m-2以上，离散系数均小于0.7，说明相对于EW 向和NW 向结构面，区域的NE 向与NNE 向结构面发育规模更大、更密集；并且与NNE 向结构面比较，NE 向结构面的平均迹长较大，中点面密度值较小.

以图4（b）中岩体的Voronoi 剖分为基础，以节理聚类优势组i为横轴，以走向平均值 αi、平均迹长值Li和中点面密度值为纵轴，将表3 中拆分露头的平面特征值赋予区块，直观呈现区域特征值的平面差异性，如图5 所示.

图5 中空白区域表示该区块没有数值，颜色深浅表征平面特征数值大小.以平均走向 αi为横轴，观察图5（a）、（b）、（c）、（d）中不同优势组的平均走向分布.在i=1 中，相对于其他区块，东北部区块偏北；在i=3 中，东北部区块北偏东，西南部区块则北偏西；其余两个优势组的走向分布并没有明显规律.

其次以平均迹长Li为横轴，与图3（b）对照分析.图3（b）中平均迹长L较小的区块位于东北部，图5（e）、（f）、（g）、（h）中位于东北部的区块，其平均迹长Li同样较小，呈现相似性特征；图3（b）中平均迹长L较大的区块位于西南部，而图5（e）、（f）、（g）、（h）中的不同优势组平均迹长发育特征不同：例如i=1（NE 向）优势组节理在O36、O37、O38区块的平均迹长值均大于9 m，而i=4（NW 向）优势组节理在同样的三个区块中缺失；i=3（NNE 向）优势组节理在O2、O37、O43区块中发育的平均迹长值大于8 m，而i=2（EW 向）向优势组节理缺失.

3.2 钻孔裂隙的深度分析

钻孔的深度区间为[0，-600 m]，以深度为纵轴，以10 m 为间隔，绘制不同深度的裂隙频数柱状图，如图6（a）所示；其次依据图2（f）中钻孔裂隙的优势组划分，将裂隙中不同优势组的走向、倾角在深度位置上展开，如图6（b）、（c）所示：

图6 钻孔裂隙优势组的深度分布图.(a) 10 m 间隔裂隙数量；(b)走向分布；(c) 倾角分布Fig.6 Depth distribution map of dominant sets of borehole cracks: (a) counts of borehole cracks at 10-m intervals;(b) strike distribution map of borehole cracks;(c) dip distribution map of borehole cracks

图6（a）中，裂隙数量随深度呈现不均匀分布特征.钻孔裂隙的线密度平均值为0.11 m-1；在[-190 m，-310 m]区域和[-460 m，-600 m]区域，裂隙发育相对稀疏，线密度均值分别为0.05 m-1和0.036 m-1；在[0 m，-190 m]区域和[-310 m，-460 m]区域，裂隙发育相对密集，裂隙线密度均值分别为0.14 m-1和0.2 m-1，线密度在[-150 m，-160 m]区域达到最大值0.60 m-1.

图6（b）可得，裂隙走向随深度呈现不均匀分布特征.钻孔裂隙中发育最显著的蓝色NNE 向优势组295.4°∠68.3°，仅在[0 m，-380 m]范围内发育，并在[-310 m，-370 m]范围内密集发育，如图中蓝色小框区域内裂隙，其平均走向为33.21°，走向随深度增加顺时针偏转；橙色近EW 向优势组189.5°∠71.2°的深度发育区间为[-10 m，-580 m]，其中[-170，-380 m]区间并没有发育，但在与蓝色小框相邻的[-380 m，-400 m]区间内密集发育，走向均值为98.31°，同样随深度增加顺时针偏转.同一优势组内部及不同优势组之间随深度呈现不均匀分布特征.

对比图6（a）与图6（b），图6（a）中裂隙线密度最大的[-150 m，-160 m]区间，在图6（b）中由三个优势组的结构面构成，并且相邻深度的裂隙并不发育，裂隙走向呈现分散特征；而图6（a）裂隙相对密集发育的[-310 m，-370 m]和[-380 m，-400 m]区间，在图6（b）中仅由NNE 向裂隙和EW 向裂隙构成.裂隙密度与裂隙走向之间呈现复杂性.

由图6（c）可得，裂隙倾角以陡倾角为主，而小于45°的缓倾角裂隙集中在[-150，-160 m]，对照图3（b），缓倾角裂隙主要为NWbW 向.

4 讨论

本文以北山沙枣园花岗岩体中分布的随机结构面为研究对象，目标岩体尺寸为600 m×600 m×600 m，以露头节理测量和钻孔裂隙编录的随机结构面数据为来源，对结构面倾角、倾向开展地质统计与聚类分析.在随机结构面的产状分析中，露头节理与钻孔裂隙的产状呈现相似性特征：钻孔裂隙中发育的四组优势产状，发育数量最多的三组优势产状在露头节理的优势产状中均有对应发育，数据占比排序相似，而没有对应的缓倾角优势产状则源于露头取样方法的局限性；而在随机结构面的位置分析中，露头节理与钻孔裂隙则充分体现了分布的不均匀性：完整露头节理的平均迹长和中点面密度呈现不均匀性，而相对于完整露头节理，拆分露头节理中不同优势组节理的平均迹长和面密度在平面分布上更不均匀；对钻孔裂隙而言，裂隙在深度分布上密度不均匀，不同优势组裂隙的深度分布范围并不一致、密度分布不均一，并随深度增加出现走向偏转的现象.

我国西北地区构造地应力的整体环境为纵向压缩，由NNE 向至NEE 向顺时针旋转[20-21].北山属于天山海西褶皱带的东延部分，主要构造线为近EW 向，采用“梯子-剪切”模式[7]解释沙枣园岩体地表多尺度结构面组合形态的形成过程.由于不同期次的构造运动对岩体深度与地表区域的影响范围不同，随机结构面在位置分析中呈现的不均匀性，考虑风化作用对岩体的影响深度达到20～35 m[22]，以及地表露头与钻孔岩体有限的揭露范围，两者的测量结果易受到测量位置的影响呈现统计差异，然而深部钻孔裂隙与岩体地表露头节理产状仍然呈现显著的相似性特征.由上述现象，得到两点启发：（1）对于地表风化严重的沙枣园岩体而言，地表露头节理调查可能反映岩体深部随机结构面的产状分布特征；（2）据地表多尺度结构面组合形态提出的“梯子-剪切”模式可能同样适用于地下岩体.

由图1 可知，研究区域内随机结构面的测量面积非常有限，这体现了野外地质调查的局限性：只能就地表发育良好、可测的露头开展节理调查.调查露头的数量多少、面积大小以及在调查区域中分布位置，均受到岩体地表的风化状态以及勘测技术的限制.从调查方法的角度，由于不同倾角结构面与钻孔和露头取样窗口相交的概率不同，Terzaghi[23]提出了不同倾角结构面取样数量的矫正概率模型.考虑沙枣园天然岩体主要发育陡倾角裂隙，而倾角与倾向是结构面不可分割的产状属性，为避免对后续结构面的产状聚类分析造成偏差，本文并没有根据倾角矫正结构面数量.

岩体随机结构面的产状纷繁复杂，产状分析的关键是确定结构面的优势分组.本文选用基于粒子群算法的K-means 方法[16]搜索结构面的优势产状.当前研究人员开发了人工蜂群[24]、遗传算法[25]等生物启发式算法来优化K-means 算法的迭代方向，使结构面的聚类结果达到全局最优.该类算法有两个缺陷：（1）基于距离衡量结构面的产状相似程度，而不是基于概率分布，不能根据局部区域的产状聚类划分推测结构面整体的分布规律；（2）算法将位于不同平面位置、产状属性相似的节理归为一类，却忽略了结构面的空间属性，以节理产状的聚类为例，将单一露头中相近产状节理拆分为两类.后续应基于von Mises-Fisher 分布[26]与结构面空间属性，开展随机结构面聚类分析研究，为后续沙枣园岩体随机结构面的三维离散裂隙网络模拟奠定理论基础.

在露头节理的平面分析中，本文考虑区域岩体的平面均质性和异质性，沿用经典的Voronoi 多边形法[17]，直观呈现了岩体地表的节理分布特征.需要进行区分的是，张开雨等[27]基于岩石的细观晶粒结构与Voronoi 多边形的几何近似，将Voronoi多边形运用在花岗岩试样破裂模拟中，这一点不同于Voronoi 多边形法的起源.

本文仅针对随机结构面的真实空间分布开展研究，没有就岩体随机结构面产状分析的确定性与位置分析的不均匀性原因开展深入分析.假设目标岩体是一个天然的完整试验试样，视其经历的地质构造运动为试验加载过程，岩体内部的随机结构面就是试样在加载过程中产生的微裂纹.李世愚等[28]发现由Mohr-Coulomb 准则所确定的破裂面裂纹角与用Griffith 准则所确定的最优裂纹角一致.Griffith 准则解析单一裂纹张拉破坏行为，而Mohr-Coulomb 准则描述岩石剪切破坏时大量微裂纹行为共同作用产生的统计结果.理论准则与经验准则的统一说明：从岩体破坏孕育的角度，加载过程中天然岩体的非均质性会产生应力集中现象，从而导致微裂纹的空间分布不均匀，而该过程服从于岩体宏观破裂面的形成过程.上述发现将随机结构面产状统计的确定性与位置分析的不均匀性统一于岩体随机结构面的建造过程，需要后续进一步开展室内力学实验进行验证.

5 结论

本文以600 m×600 m×600 m 天然岩体中孕育的随机结构面为研究对象，目标岩体的钻孔和露头RQD 值均在94 以上，岩体质量完整，针对岩体地表的露头节理数据和岩体深部的钻孔裂隙编录数据开展产状统计、聚类分析以及钻孔裂隙深度分析，并开发出基于Voronoi 多边形法与结构面优势组的区域露头节理平面分析技术，得到如下结论：

(1) 沙枣园岩体调查区域的露头节理与钻孔裂隙的优势组呈现相似性特征.钻孔裂隙的优势组为295.4°∠68.3°、189.5°∠71.2°、235.6°∠69.2°和66.0°∠33.0°，数据占比分别为41.2%、23.5%、20.6%、14.7%，露头节理的优势组为282°∠80°、320°∠79°、180°∠77°和231°∠67°，数据占比分别为37.3%、31.9%、19.1%和11.8%.

(2) 沙枣园岩体调查区域中完整露头平均迹长和中点面密度值呈现平面不均匀分布特征，区域东北部、中部、西南部的平均迹长值依次在5、8.5、11 m 左右，呈现增大趋势；中点面密度值呈现中间大，周边小的分布，最大值达0.12 m-2；节理拆分后的NE 向和NNE 向区域迹长平均值为9.52 m 和8.18 m、区域中点面密度均值为0.0246 m-2和0.0348 m-2均大于EW 向和NW 向优势组.

(3) 沙枣园岩体调查区域的钻孔裂隙的数量与走向呈现深度不均匀分布，以10 m 为间隔，裂隙平均线密度为0.11 m-1，最大值为0.60 m-1，在[-310 m，-370 m] 线密度达到0.22 m-1，发育平均走向33.21°的NNE 向裂隙，在[-380，-400 m]区间线密度达到0.4 m-1，发育平均走向为98.31°的EW 向裂隙.

本文是将随机结构面统计分析与空间分布相结合的初步研究成果，为后续开展沙枣园岩体随机结构面的分区域三维网络模拟提供数据支撑.