CN112200420A - 基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法 - Google Patents

Abstract

一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，属于采空区稳定性评价领域，其步骤包括：(1)结构面三维激光扫描快速获取；(2)结构面聚类分析；(3)基于BQ指标的岩体质量计算；(4)岩体三维裂隙网络模型生成和剖切；(5)RQD_t各向异性图绘制；(6)基于BQ反演的最佳阈值t求解方法；(7)RQD_t各向异性求解方法；(8)Mathews稳定图法的改进方法；(9)改进Mathews稳定图评价方法；(10)围岩稳定性分析方法；(11)空区稳定性评价方法。本发明实现了RQD_t的最佳阈值t、RQD_t各向异性的求解和Mathews稳定图法的改进，实现了改进Mathews稳定图法和数值模拟相结合的空区围岩稳定性评价。本发明方法明确，适用于采空区围岩稳定性评价。

Description

基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价 方法

技术领域

本发明涉及一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，特别的是本发明基于三维激光扫描、近邻传播算法、BQ指标、裂隙网络模型和广义RQD理论，反演计算出RQD_t最佳阈值t的范围和最佳阈值t值，给出了RQD_t各向异性的求解方法，并基于RQD_t的各向异性，给出了Mathews 稳定图法的改进方法和改进后的Mathews稳定图评价方法，再结合数值模拟方法，提供了一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，属于采空区稳定性评价领域。

背景技术

金属矿山经过长时间地下开采，会形成大规模的采空区，采空区在地下相互贯通、重叠和集中，彼此间应力扰动，相互制约，形成了复杂的群体空间关系，岩体稳定状态极为复杂。

随着我国深部矿产资源开发力度的加大，采空区的数量不断增多，规模不断加大，开采条件不断恶化，安全生产形势更加严峻，事故起数与伤亡人数一直居高不下。特别是自20世纪80年代以来，由于开采秩序较为混乱，非法、无规划的乱采滥挖遗留了大量不明采空区，成为影响矿山安全生产的主要危险源之一。这些采空区，由于存在时间久，经过多次破坏，岩体结构已比较危险，随着开采深度的不断增加，造成岩体失稳、地表塌陷、岩爆、巷道破坏等灾害的可能性及不确定逐渐增加，加剧了后续矿体开采的难度。因此，开展采空区稳定性的分析评价具有重要意义。

Mathews稳定图法是采空区围岩稳定性分析方法的一种，是在1980年由Mathews等人通过分析大量的工程实例，提出的基于Q系统的稳定图法。Mathews稳定图法在评价过程中主要考虑采场特征、地质条件、节理产状等条件，通过计算岩体稳定数N和水力半径R，并将这两个因子绘制在划分为预测稳定区、潜在不稳定区和崩落区的图上。

在Mathews稳定图法的改进方面，1995年Stewart和Forsyth通过收集不同采矿深度的数据，验证和修正了稳定图法，将其划分为四个分区。1998年Potvin基于数据实例对稳定图法做了改进，在图中增加了节理方位修正系数和重力调整系数。2000年Pakalnis等在RMR的基础上总结出了临界跨度图表法。Trueman 等采用对数回归的方法，对稳定区和严重破坏区进行了重新定义。2004年Mawdeskey通过对405个实例进行逻辑回归分析，给出了稳定区、破坏与严重破坏区的等概率图。

但是岩体中存在着大量的结构面和裂隙，显著地影响岩体的力学特性、破碎程度和稳定性，是影响矿房稳定性的重要因素之一。岩体中的结构面和裂隙，是造成岩体质量各向异性的重要因素，而岩体质量的各向异性特征，在Mathews稳定图法中，尚未被考虑。

岩体的外观和质量强烈的依赖于岩体的各向异性程度，岩体的稳定性是由岩体质量所决定的。RQD是表征岩体质量的重要指标，RQD具有各向异性。如采用钻孔的方式获得RQD，不同的岩体部位得到的结果完全不同，RQD值结果将取决于方向，当钻孔方向与主要节理组平行时，得到较高的RQD值，当钻孔方向垂直于主要节理组方向时，会得到较低的RQD值。RQD的各向异性，从钻孔RQD到广义RQD的计算公式的拓展也可以看出。Deer在1964年提出了钻孔RQD的概念，由于钻孔RQD在应用时存在着如下两个缺点：对于不同工程规模岩体，100mm的阈值选取是否合理；钻孔钻探方向具有局限性，获得的RQD 只能反映局部岩体情况，不能反映出RQD各向异性。因此，部分学者引入阈值t，提出了广义RQD的概念，即对任一间距阈值t，把沿某一测线方向大于t的间距之和与测线总长之比的百分数定义为广义RQD，用RQD_t表示。广义RQD的引入，让RQD各向异性的求解成为了可能。

由RQD_t概念可以知道，阈值t是广义RQD的一项特别重要的参数，是广义RQD能否真实反映出岩体质量的关键，但是阈值t具有任意性，不具有唯一性，因此寻找到能够表征出岩体质量的最佳阈值t，具有非常重要的科学意义和研究意义。

RQD_t的阈值t受测线方向、结构面形态和分布特征的影响，而岩体中广泛发育的结构面，是破坏岩体连续性、完整性，控制岩体力学特征和稳定性的重要因素，对岩体质量起着控制性作用。结构面的发育模式和分布形态非常复杂，但同时，不同的节理组之间、节理和断层之间，又存在着一定的成生关系，构成某种特定的组合，表现出一定的规律性。因此，对结构面进行准确和有效的描述分析，获得结构面的产状和分布特征，是研究岩体质量的基础，也是最佳阈值t的研究基础。三维激光扫描技术作为一种高效的三维空间信息获取手段，在获取结构面方位和规模信息方面具有很大优势。近邻传播聚类算法在结构面产状聚类分析中具有较好效果，它具有不受初始聚类中心影响、计算效率高的优点，已经广泛应用于很多领域。因此，基于三维激光扫描技术和近邻传播算法，可以实现结构面的快速识别获取和聚类分析，为RQD_t的最佳阈值t的研究，提供数据基础。

国内外学者针对阈值t的研究，主要体现在以下两个方面：不同阈值t下的RQD_t计算和最佳阈值t的研究。在不同阈值t下的RQD_t计算方面，已有的研究主要是探讨阈值t的变化对RQD_t值的影响，主要是通过模拟岩体裂隙网络模型，布置虚拟钻孔，以求解不同RQD_t值，研究的载体是裂隙网络模型。但在这个方面的研究中，主要是为了研究RQD_t的各向异性特征，对最佳阈值t的研究鲜有涉及。

在最佳阈值t的研究方面，有些学者也曾对最佳阈值t展开过探讨。如有的学者基于三维结构面网络模拟技术，运用分形理论分别计算其结构面分布的分形维数，通过不断改变广义RQD的阈值，得到不同阈值下的RQD_t值，将各剖面的分形维数与广义RQD值对比分析，为准确选取广义RQD的最佳阈值提供依据。有的学者基于修改后的块度指数MBi，建立了35种假想的三维裂隙网络模型，测量了具有不同阈值的广义RQD值，确定广义RQD的最佳阈值。这两个方面的最佳阈值t的研究，适用条件都有一定的局限性，都是在特定理论和背景下的一种最佳阈值t求解方法，当背景或者模型变化了之后，其最佳阈值t 将不再具有最佳性。而且，由于分形维数或者块度指数本身并不具备表征岩体质量的功能，通过分形维数或者块度指数获得的最佳阈值t，能否反映真实的岩体质量，还有值得商榷的地方。

因此，找到并提出一种最能反映出岩体质量的RQD_t最佳阈值t的求解方法，具有非常重要的科学和研究意义，也是RQD_t阈值t研究中亟需解决的问题。

RQD的各向异性直接影响岩体质量，RQD各向异性对岩体质量的影响机理，目前还尚未探索清楚。由于在阈值t方面，还尚未有学者给出最佳阈值t的求解计算方法，因此也没有得到基于最佳阈值t的RQD_t各向异性计算公式，也尚未获得最能反映岩体质量的RQD_t各向异性求解方法。而获得RQD_t各向异性求解方法是研究岩体质量各向异性的基础和前提。

鉴于此，本发明提出了一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法。

发明内容

为了在采空区稳定性评价中考虑岩体质量的各向异性特征和岩体的破坏情况，本发明提供了一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法。本发明基于三维激光扫描、结构面聚类分析、BQ指标、裂隙网络模型和广义RQD理论，反演计算出RQD_t最佳阈值t的范围和最佳阈值t值，给出了RQD_t各向异性的求解方法，并基于RQD_t的各向异性，给出了Mathews稳定图法的改进方法和改进后的 Mathews稳定图评价方法，再结合数值模拟方法，提供了一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法。

为了解决上述技术问题，本发明提供如下的技术方案：

一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，所述方法包括以下步骤：

(1)结构面三维激光扫描快速获取，过程如下：

1.1：根据扫描目标和场地条件，选择扫描机位点，架设三脚架，架设中要确保仪器按照一定的扫描路线可以完整的获取边坡岩体的三维空间点云信息，同时要尽可能保证三脚架台面水平，并放置控制靶；

1.2：放置扫描仪主机在三脚架台面，固定旋钮，通过粗调脚架及微调扫描仪底座使主机气泡居中，设置扫描仪端口参数；

1.3：启动扫描控制软件，配置扫描仪相关参数，进入扫描仪控制界面，规划扫描角度，根据扫描目标设置扫描范围，调整相机配置参数，获取扫描目标图像；

1.4：固定扫描范围，获取扫描间距，设定采样间距，开始数据获取，并实时查看扫描点云数据及彩色信息情况，根据扫描成果随时调整扫描参数设定；

1.5：导出结构面点云数据；

(2)结构面聚类分析；

(3)基于BQ指标的岩体质量计算，过程如下：

3.1：根据结构面参数计算岩体完整性系数，公式如下：

式中：Jv为岩体体积节理数，单位条/m³；

3.2：Jv计算公式如下：

式中：L1，L2，...，Ln为垂直于结构面测线长度；N1，N2，...，Nn为同组结构面数目；

3.3：根据岩石单轴抗压强度值和岩体完整性系数值，计算BQ值：

BQ＝90+3R_c+250K_v

式中：Rc是岩石单轴抗压强度；Kv为岩体完整性系数；

3.4：在应用BQ计算公式过程中，遵循以下条件：

当Rc＞90Kv+30时，以Rc＝90Kv+30和Kv代入计算BQ值；

当Kv＞0.04Rc+0.4时，以Kv＝0.04Rc+0.4和Rc代入计算BQ值；

3.5：根据地下水、软弱结构面产状和天然应力影响对BQ进行修正，修正公式如下：

[BQ]＝BQ-100(K₁+K₂+K₃)

式中：K1为地下水影响修正系数；K2为软件结构面产状影响修正系数；K3为天然应力影响修正系数；

3.6：得到修正后的BQ岩体分级结果；

(4)岩体三维裂隙网络模型生成和剖切；

(5)RQD_t各向异性图绘制；

(6)基于BQ反演的最佳阈值t求解方法；

(7)RQD_t各向异性求解方法；

(8)Mathews稳定图法的改进方法；

(9)改进Mathews稳定图评价方法；

(10)围岩稳定性分析方法，过程如下：

10.1：根据矿房平面图和剖面图，利用CAD软件，建立矿房三维模型；

10.2：将三维模型导入到有限元分析软件中，赋予三维模型材料参数，包括岩体强度、密度、弹性模量、泊松比、粘聚力、摩擦角和容重；

10.3：赋予三维模型边界条件和初始应力；

10.4：进行矿房回采的数值模拟，首先进行未回采状态的数值模拟，得到初始应力场，然后进行矿房回采的数值模拟，得到采空区的应力场；

10.5：对矿房回采后的模拟结果进行分析，分析区域为矿房顶板、上下盘围岩和矿柱；

10.6：岩体剪切破坏的判断准则选用以主应力形式表示的Mohr-Coulomb屈服准则，当f_s＞0时，岩体发生剪切破坏，公式如下：

式中：σ₁、σ₃分别为最大主应力和最小主应力；c为粘聚力；

为摩擦角；

10.7：岩体张拉破坏的判断准则，根据岩体所受拉应力σ_t与抗拉强度R_t之间的关系确定，当σ_t≥R_t时， F≥0，即认为岩体发生张拉破坏，公式如下：

F＝σ_t-R_t

式中：σ_t为岩体所受拉应力；R_t为岩体抗拉强度；

10.8：根据分析结果，得到矿房顶板、上下盘围岩和矿柱的破坏模式；

(11)空区稳定性评价方法。

进一步，所述步骤(11)中，空区稳定性评价方法的过程如下：

11.1：以改进Mathews稳定图法得到的围岩稳定性评价结果作为采空区稳定性的初始评价结果；

11.2：利用数值模拟得到的围岩稳定性分析结果，对初始评价结果进行修正，过程如下：

11.2.1：当改进Mathews稳定图法评价的采空区围岩处于崩落区时，则该采空区处于崩落状态；

11.2.2：当改进Mathews稳定图法评价的采空区围岩处于破坏区时，则该采空区处于破坏状态；

11.2.3：当改进Mathews稳定图法评价的采空区围岩处于稳定区时，若数值模拟得到的采空区顶板、间柱、上下盘围岩未发生破坏，则该采空区处于稳定状态，若数值模拟得到的采空区顶板、间柱、上下盘围岩发生破坏，则该采空区有可能处于破坏状态；

11.3：得到采空区的稳定性评价结果。

进一步，所述步骤(9)中，改进Mathews稳定图评价方法的过程如下：

9.1：选取用于开展上下盘围岩稳定性评价的矿房；

9.2：测量出矿房上下盘围岩宽度X和长度Y，计算出岩石应力系数A；节理产状调整系数B和重力调整系数C；

9.3：测量并求解出矿房上下盘的节理组数J_n；节理粗糙度系数J_r；节理蚀变系数J_a；节理水折减系数J_w和SRF为应力折减系数；

9.4：求解出最佳阈值t下的RQD_tθ、RQD_tmin和

计算出相关拟合系数a和b，求解出考虑RQD_t各向异性特征的Q′值；

9.5：计算出矿房稳定数N和水力半径R；

9.6：计算出稳定区-破坏边界的曲线公式和破坏-崩落边界的曲线公式，并绘制到改进后的Mathews 稳定性图上；

9.7：将计算出的矿房稳定数N和水力半径R绘制到改进后的Mathews稳定图上，根据绘制的散点所处的区域，给出矿房上下盘围岩的稳定性评价结果，若散点处于稳定区，则矿房围岩评价结果为稳定，若散点处于破坏或主要破坏区，需再根据矿房围岩破坏和崩落的概率，给出矿房围岩的破坏概率和崩落概率，若散点处于崩落区，则矿房围岩评价结果为崩落。

进一步，所述步骤(8)中，Mathews稳定图法的改进方法的过程如下：

8.1：Mathews稳定图法的设计过程，过程如下：

8.1.1：Mathews稳定图方法的设计过程是以两个因子-稳定数N和水力半径R的计算为基础，然后将这两个因子绘制在划分为预测稳定区、潜在不稳定区和崩落区的图上，稳定数N代表岩体在给定应力条件下维持稳定的能力，水力半径R反映了采空区尺寸和形状，Mathews稳定图法计算公式为：

N＝Q′ABC

式中：X、Y为矿房上下盘围岩宽度和长度；Q′为根据勘测图或钻孔岩芯记录计算出的结果；A为岩石应力系数；B为节理产状调整系数；C为重力调整系数；

8.1.2：岩石应力系数A由完整岩石单轴抗压强度与采场中线采矿产生的地应力的比值确定，公式如下：

当

当

当

8.1.3：节理产状调整系数B的值是通过采场面倾角与主要节理组的倾角之差来度量；

8.1.4：重力调整系数C反映了在重力的影响下采场面产状对采场矿岩稳定性的影响，其大小取决于采场顶板暴露表面的崩落、滑落以及边帮的滑落等，重力调整系数C和采场表面倾角α的关系由下式确定：

C＝8-7cosα

8.2：Mathews稳定图法的改进方法，过程如下：

8.2.1：在求解Q′值时，用RQD_t值替代钻孔岩芯记录的RQD值，并将RQD_t的各向异性特征纳入考虑范畴，根据RQD_t各向异性求解公式，得到考虑RQD_t各向异性特征的Q′值，公式如下：

式中：RQD_tθ为角度θ下的RQD_t值；RQD_tmin为阈值t下的RQD最小值；

为RQD_t的均值；a，b为相关拟合系数；J_n为节理组数；J_r为节理粗糙度系数；J_a为节理蚀变系数；J_w为节理水折减系数；SRF 为应力折减系数；

8.2.2：结合等概率轮廓图对Mathews稳定图进行重新绘制，稳定区-破坏边界的曲线公式如下：

logN＝1.8206logR+1.618

破坏-崩落边界的曲线公式如下：

logN＝1.8076logR-3

8.2.3：改进后Mathews稳定性图分为三个区域：稳定区、破坏或主要破坏区和崩落区，处在稳定-破坏边界线上的工程，采场57％概率稳定，43％概率破坏，0％概率崩落；处在崩落-破坏边界线上的工程，采场0％概率稳定，5％概率破坏，95％概率崩落；处在稳定区的工程是稳定的；处于破坏或主要破坏区的工程，57％～0％的概率稳定，43％～5％的概率破坏，0％～95％的概率崩落；处在崩落区的工程，将持续发生崩落。

进一步，所述步骤(7)中，RQD_t各向异性求解方法的过程如下：

7.1：基于最佳阈值t和剖切的三个二维裂隙网络模型，求解出最佳阈值t下的RQD_t值，每个裂隙网络模型求解出36个RQD_t值；

7.2：计算出每个裂隙网络模型上RQD_t的最大值RQD_tmax、最小值RQD_tmin和均值

7.3：根据研究方向和RQD_tmin方向位置关系，结合方差进行修正，提出各向异性条件下RQD_t的计算公式如下：

式中：

为RQD_t的均值，RQD′_t为RQD_t的修正系数；

7.4：修正系数RQD′_t的求解考虑研究方向在某一具体方位角下的RQD_tθ值和方差D，按如下方法修正，当RQD_tθ＝RQD_tmin时，RQD′_t＝-D，D为方差；当RQD_tθ＝RQD_tmax时，RQD′_t＝D；当

时，RQD′_t＝0；

7.5：根据修正系数RQD′_t的求解过程，提出RQD′_t修正公式如下：

式中：RQD_tθ为θ角下的RQD_t值，RQD_tmin为阈值t下的RQD_t最小值，a，b为相关系数；

7.6：提出RQD_t各向异性计算公式如下：

7.7：相关系数a和b的求解，过程如下：

7.7.1：计算出每个裂隙网络模型上RQD_t的最大值RQD_tmax、最小值RQD_tmin和均值

7.7.2：根据RQD_tθ与RQD_tmax、RQD_tmin和

以及方差的关系，求解出三组RQD′_t值，并基于三组RQD′_t值，绘制出散点图；

7.7.3：根据散点图做出拟合曲线，曲线截距即为a值，斜率即为b值，曲线方程即为RQD′_t的修正公式；

7.8：将求解出的a、b值，带入到修正系数公式及RQD_t各向异性计算公式中，得到与角度θ有关的RQD_t各向异性公式；

7.9：根据RQD_t各向异性公式，求解出任意角度的RQD_t值。

进一步，所述步骤(6)中，基于BQ反演的最佳阈值t求解方法的过程如下：

6.1：基于BQ指标反演RQD_t范围，过程如下：

结合BQ分级计算出的岩体质量级别，查找《岩石质量指标》表，反演确定出该岩体级别下，RQD_t范围值；

6.2：最佳阈值t求解方法，过程如下：

6.2.1：在三维裂隙网络模型上，过中心点O，以任意角度剖切三个剖面，得到三个二维裂隙网络模型，导出二维裂隙网络模型及数据；

6.2.2：针对每一个二维裂隙网络模型，设置不同的阈值t，求解出不同阈值t下的RQD_t值；

6.2.3：导出不同阈值t下的RQD_t值，计算出RQD_t均值；

6.2.4：以阈值t为横坐标，以RQD_t的均值为纵坐标，绘制RQD_t随阈值t变化的散点图；

6.2.5：根据散点图，设置拟合方程，拟合RQD_t随阈值t变化的曲线图；

6.2.6：将反演出的RQD_t范围值，带入到拟合出的RQD_t随阈值t变化的曲线图中，结合函数方程和曲线图，求解出在该RQD_t范围内阈值t的范围，共得到三组阈值t的范围；

6.2.7：针对三组阈值t的范围，取其范围的交集，作为最佳阈值t的范围；

6.2.8：以最佳阈值t范围的中点值作为最佳阈值t值，得到最佳阈值t；

6.2.9：输出最佳阈值t的范围和最佳阈值t值。

进一步，所述步骤(5)中，RQD_t各向异性图绘制的过程如下：

5.1：RQD_t求解计算，过程如下：

5.1.1：RQD_t理论公式如下：

式中：x_i表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，RQD_t代表对应阈值t的岩石质量指标，即阈值t下的RQD_t值；

5.1.2确定二维裂隙网络模型的剖面中心点O、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，求解出中心点O坐标(X₀，Y₀)和边界方程；

5.1.3：以O为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度L等于O点到裂隙网络模型边界的距离，用L0～L35表示，求解出测线方程；

5.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(x_a，y_a)，将测线方程与边界方程依次连立，判断测线是否与边界相交；

5.1.5：求出测线与裂隙网络模型的边界方程交点(X_Z，α，Y_Z，α)；

5.1.6：确定测线所在区间，原理如下：

5.1.7：根据裂隙网络模型中每条节理的起点坐标(X_b，Y_b)和终点坐标(X_c，Y_c)，建立相应的解析方程；

5.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(X_j，Y_j)范围，如果交点符合a＜X_j＜c且b＜Y_j＜d，则记录该交点，遍历所有节理方程求出所有m个交点；

5.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；

5.1.10：计算相邻交点的距离d_i；

5.1.11：输入一个阈值t；

5.1.12：循环比较d_i和t的大小，设初始l_t＝0，规则如下：

l_t＝l_t+d_i，若d_i＞t

5.1.13：求解每条测线对应的RQD_t值，用m_α表示，以及对应的测线起点与终点的距离l_α；

5.1.14：循环求出每条测线所对应的m_α，获得36个方向上的RQD_t值；

5.2：RQD_t各向异性图绘制，过程如下：

5.2.1：将36个RQD_t值，按角度顺序依次排序；

5.2.2：以O点为圆心，以l为半径画圆，在射线角度为α上找到与圆心距离为

的点标出来；

5.2.3：依次连接36组射线的端点，若某一条射线上的RQD_t为0，则取圆心；

5.2.4：绘制出RQD_t的各向异性图；

5.2.5：输出RQD_t各向异性图；

5.3：不同阈值t下RQD_t各向异性图绘制，过程如下：

5.3.1：输入不同的阈值t，求解出对应的RQD_t值；

5.3.2：将不同阈值t下的RQD_t值，以角度为变量，以RQD_t值为函数，绘制到同个坐标系下；

5.3.3：得到不同阈值t下的RQD_t各向异性图；

5.3.4：输出不同阈值t下的RQD_t各向异性图；

5.4：输出RQD_t的值。

进一步，所述步骤(4)中，岩体三维裂隙网络模型生成和剖切的过程如下：

4.1：随机数求解，过程如下：

4.1.1：产生随机数的数学方法应满足以下条件：产生的随机数列应均匀分布在(0，1)区间；序列之间应无相关性；随机序列有足够长的重复周期，在计算机上产生的速度快，占有的内存空间小，具有完全可重复性；

4.1.2：利用Monte Carlo方法是根据确立的结构面几何概率模型，再现出服从这种模型的结构面网络模型，在(0，1)区间上生成均匀分布随机变量，利用这些均匀随机变量产生服从其他分布的随机数；

4.1.3：节理几何参数的密度函数有正态分布、对数正态分布、负指数分布、均匀分布四种；

4.1.4：根据求得的随机数，确定用于生成节理的基本几何参数；

4.2：岩体三维裂隙网络模型生成，过程如下：

4.2.1：根据结构面数据自动统计结果和求得的随机数，将每组结构面的数据保存到一个文本文件中，用st.dat表示；

4.2.2：st.dat数据内容格式依次为：每条结构面圆盘中心点坐标(x，y，z)，圆盘半径D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ和节理分组；

4.2.3：为区分不同组别的结构面，对相同组的结构面圆盘赋予相同的颜色，用数列GID表示；

4.2.4：利用Matlab软件，编写程序，读取结构面数据文件st.dat，运行后自动生成岩体三维裂隙网络模型；

4.2.5：得到岩体三维裂隙网络模型；

4.3：二维裂隙网络模型切割，过程如下：

4.3.1：在三维裂隙网络模型上，结合Matlab软件编程工具，以三维裂隙网络模型中心点为中心，实现任意角度的剖面切割功能；

4.3.2：得到穿过中心点的任意角度的二维裂隙网络模型；

4.3.3：在三维裂隙网络模型上，结合Matlab软件编程工具，在三维裂隙网络模型任意位置上，实现任意角度和方位的剖面切割功能；

4.3.4：得到任意角度和方位的二维裂隙网络模型；

4.3.5：将切割剖面上的数据，保存到st1.dat文件中，此时剖面处于三维坐标系下，文件中数据格式自左向右依次为：节理中心点坐标(x，y，z)，节理长度D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ；

4.3.6：将三维坐标系转化成二维坐标系，并将二维剖面数据保存到st2.dat文件中，文件中数据格式自左向右依次为：节理中心点坐标(x，y)，节理长度D，倾角DA，倾向DD，厚度thin，法向方向NX，NY，NZ；

4.4：输出三维裂隙网络模型的数据和任意二维裂隙网络模型的数据。

进一步，所述步骤(2)中，结构面聚类分析的过程如下：

2.1：点云数据处理，过程如下：

2.1.1：导入激光扫描获得的结构面点云数据；

2.1.2：计算拓扑构造后的点云中当前点与相邻点的距离与距离均值，通过距离阈值对点云数据中噪声点进行识别、剔除；

2.1.3：根据三维激光扫描仪自身空间坐标位置和现场结构面产状方位，确定点云数据的空间三维坐标；

2.1.4：基于下半球等角度投影方法进行点云数据的转换，将以倾向和倾角表示的节理产状数据转换为以节理单位法向量表示的结构面产状数据；

2.1.5：得到以单位法向量表示的结构面数据；

2.2：近邻传播算法聚类分析，过程如下：

2.2.1：设结构面的实测样本数量为N，每个样本数据的倾向为X_i，倾角为Y_i，以每个样本数据的倾向X_i，倾角Y_i作为一个聚类，确定一个初始聚类中心，共得到N个初始聚类中心；

2.2.2：通过相似性度量准则，遍历所有样本数据，计算每个样本数据距离聚类中心的距离，并将每个样本数据分配到距离它最近的聚类中心，得到N组数据；

2.2.3：对于每组数据，通过特征模量分析方法，求解计算每组数据的聚类中心，假设某组内存在l个数据，首先，按如下公式计算矩阵S：

式中：(x_i，y_i，z_i)为任意结构面的单位法向量，i∈(1，l)；

2.2.4：其次，求解矩阵S的特征值(τ₁，τ₂，τ₃)和特征向量(ξ₁，ξ₂，ξ₃)，其中τ₁＜τ₂＜τ₃，最大特征值对应的特征向量ξ₃为组内l个向量的平均向量，将ξ₃作为新的聚类中心；

2.2.5：针对所有样本数据，重复计算每个样本数据距离聚类中心的距离、矩阵S以及特征值和特征向量，直到所有聚类中心的位置都固定，确定出结构面的分组；

2.2.6：将以单位法向量表示的结构面产状数据转换为以倾向、倾角表示的结构面产状数据；

2.2.7：对结构面产状数据进行统计分析，计算结构面倾角的平均值m与标准差σ，计算倾角数据的稳健区间[m-σ，m+σ]；

2.2.8：判断样本数据的初始聚类中心的倾向X_i和倾角Y_i是否落在稳健区间稳健区间[m-σ，m+σ]，若是，则聚类分析完成；若不是，则需要对样本数据重新聚类，直到初始聚类中心的倾向X_i和倾角Y_i均落在稳健区间内[m-σ，m+σ]；

2.3：极点图绘制，过程如下；

2.3.1：基于结构面法向产状数据，根据结构面空间赤平投影图的纵剖面原理，求解出所有结构面法线的赤平投影坐标点(x_n，y_n)；

2.3.3：绘制一条直径为单位长度的基圆，绘制出铅直和水平两条直径，并标出E、S、W、N；

2.3.4：将所有结构面的赤平投影坐标(x_n，y_n)，绘制在基圆图上；

2.3.5：实现结构面极点图的绘制；

2.4：结构面统计分析，过程如下：

2.4.1：确定样本分区区间m；

2.4.2：求解样本极差

2.4.3：计算每个分区区间M_m：

2.4.4：确定样本落在每个分区区间里的概率，先利用计算机循环语言统计落在每一个区间的样本个数 N_m，结合样本总数N，计算样本数概率P_m；

2.4.5：求解样本均值

2.4.6：求解样本方差S²，其中S为标准差；

2.4.7：根据概率P_m值，绘制出每组结构面的倾向、倾角、迹长、间距和断距的概率分布形态；

2.5：输出结构面产状的分组信息，包括每组结构面的倾向、倾角、迹长、间距和断距的均值和方差。

本发明具有以下有益效果：

1、本发明基于三维激光扫描、结构面聚类分析、BQ指标、裂隙网络模型和广义RQD理论，提供了一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法；

2、本发明实现了结构面快速三维激光扫描和近邻传播算法聚类分析；

3、本发明实现了岩体BQ指标分级；

4、本发明实现了三维裂隙网络模型生成和二维剖面模型生成；

5、本发明实现了RQD_t求解计算和RQD_t各向异性图的绘制；

6、本发明实现了RQD_t的最佳阈值t求解；

7、本发明实现了RQD_t各向异性的求解；

8、本发明实现了考虑RQD_t的各向异性特征的Mathews稳定图法的改进；

9、本发明实现了改进Mathews稳定图法的空区围岩稳定性评价；

10、本发明实现了基于数值模拟的围岩稳定性分析；

11、本发明方法明确，提供的空区稳定性评价方法实现了在评价空区围岩稳定性中考虑岩体质量的各向异性特征和岩体的破坏情况。

附图说明

图1是结构面极点图。

图2是二维裂隙网络图。

图3是不同阈值t下的RQD_t各向异性图。

图4是RQD_t随阈值t的拟合曲线。

图5是RQD′_t的拟合曲线。

图6是矿房围岩稳定性评价结果图。

具体实施方式

下面参照附图对本发明做进一步说明。

参照图1～图6，一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，包括以下步骤：

1)结构面三维激光扫描快速获取，过程如下：

1.1：根据扫描目标和场地条件，选择扫描机位点，架设三脚架，架设中要确保仪器按照一定的扫描路线可以完整的获取边坡岩体的三维空间点云信息，同时要尽可能保证三脚架台面水平，并放置控制靶；

1.2：放置扫描仪主机在三脚架台面，固定旋钮，通过粗调脚架及微调扫描仪底座使主机气泡居中，设置扫描仪端口参数；

1.3：启动扫描控制软件，配置扫描仪相关参数，进入扫描仪控制界面，规划扫描角度，根据扫描目标设置扫描范围，调整相机配置参数，获取扫描目标图像；

1.4：固定扫描范围，获取扫描间距，设定采样间距，开始数据获取，并实时查看扫描点云数据及彩色信息情况，根据扫描成果随时调整扫描参数设定；

1.5：导出结构面点云数据；

2)结构面聚类分析，过程如下：

2.1：点云数据处理，过程如下；

2.1.1：导入结构面点云数据；

2.1.2：计算拓扑构造后的点云中当前点与相邻点的距离与距离均值，通过距离阈值对点云数据中噪声点进行识别、剔除；

2.1.3：根据三维激光扫描仪自身空间坐标位置和现场结构面产状方位，确定点云数据的空间三维坐标；

2.1.4：基于下半球等角度投影方法进行点云数据的转换；

2.1.5：将以倾向α_d和倾角β_d表示的节理产状数据转换为以节理单位法向量表示的结构面产状数据，设α_n和β_n分别为结构面单位法向量的倾伏向和倾伏角，对于任意结构面的单位法向量表示为 X＝(x₁，x₂，x₃)，此时半球面上每个点都对应一个节理产状，公式为：

X＝(x₁，x₂，x₃) (1)

α_d∈(0，360)，β_d∈(0，90) (4)

2.1.6：得到以单位法向量表示的结构面数据；

2.2：近邻传播算法聚类分析，过程如下：

2.2.1：设结构面的实测样本数量为N，每个样本数据的倾向为X_i，倾角为Y_i，i∈(1，N)；以每个样本数据的倾向X_i，倾角Y_i作为一个聚类，确定一个初始聚类中心，共得到N个初始聚类中心；

2.2.2：通过相似性度量准则，遍历所有样本数据，计算每个样本数据距离聚类中心的距离，并将每个样本数据分配到距离它最近的聚类中心，得到N组数据；

2.2.3：对于每组数据，通过特征模量分析方法，求解计算每组数据的聚类中心，假设某组内存在l个数据，聚类中心按如下方法求解：

2.2.4：首先，按如下公式计算矩阵S

式中：(x_i，y_i，z_i)为任意结构面的单位法向量，i∈(1，l)；

2.2.5：其次，求解矩阵S的特征值(τ₁，τ₂，τ₃)和特征向量(ξ₁，ξ₂，ξ₃)，其中τ₁＜τ₂＜τ₃，最大特征值对应的特征向量ξ₃为组内l个向量的平均向量，将ξ₃作为新的聚类中心；

2.2.6：针对所有样本数据，重复计算每个样本数据距离聚类中心的距离、矩阵S以及特征值和特征向量，直到所有聚类中心的位置都固定，确定出结构面的分组；

2.2.7：将以单位法向量表示的结构面产状数据转换为以倾向、倾角表示的结构面产状数据；

2.2.8：对结构面产状数据进行统计分析，计算结构面倾角的平均值m与标准差σ，计算倾角数据的稳健区间[m-σ，m+σ]；

2.2.9：判断样本数据的初始聚类中心的倾向X_i和倾角Y_i是否落在稳健区间稳健区间[m-σ，m+σ]，若是，则聚类分析完成；若不是，则需要对样本数据重新聚类，直到初始聚类中心的倾向X_i和倾角Y_i均落在稳健区间内[m-σ，m+σ]；

2.3：极点图绘制，过程如下；

2.3.1：基于结构面法向产状数据，根据结构面空间赤平投影图的纵剖面原理，设A’点为该平面法线的赤面投影，结合赤平投影原理，计算出A’在赤平投影图上的坐标x_n和y_n，公式如下：

2.3.2：求解出所有结构面法线的赤平投影坐标点(x_n，y_n)；

2.3.3：绘制一条直径为单位长度的基圆，绘制出铅直和水平两条直径，并标出E、S、W、N；

2.3.4：将所有结构面的赤平投影坐标(x_n，y_n)，绘制在基圆图上；

2.3.5：实现结构面极点图的绘制，如图1所示；

2.4：结构面统计分析，过程如下：

2.4.1：确定样本分区区间m；

2.4.2：求解样本极差

2.4.3：计算每个分区区间M_m：

2.4.4：确定样本落在每个分区区间里的概率，先利用计算机循环语言统计落在每一个区间的样本个数 N_m，结合样本总数N，计算样本数概率P_m：

2.4.5：求解样本均值

2.4.6：求解样本方差S²，其中S为标准差：

2.4.7：根据概率P_m值，绘制出每组结构面的倾向、倾角、迹长、间距和断距的概率分布形态；

2.5：输出结构面产状的分组信息，包括每组结构面的倾向、倾角、迹长、间距和断距的均值和方差；

3)基于BQ指标的岩体质量计算，过程如下：

3.1：根据结构面参数计算岩体完整性系数，公式如下：

式中：Jv为岩体体积节理数，单位条/m³；

3.2：Jv计算公式如下：

式中：L1，L2，...，Ln为垂直于结构面测线长度；N1，N2，...，Nn为同组结构面数目；

3.3：根据岩石单轴抗压强度值和岩体完整性系数值，计算BQ值：

BQ＝90+3R_c+250K_v (15)

式中：Rc是岩石单轴抗压强度；Kv为岩体完整性系数；

3.4：在应用BQ计算公式过程中，遵循以下条件：

当Rc＞90Kv+30时，以Rc＝90Kv+30和Kv代入计算BQ值；

当Kv＞0.04Rc+0.4时，以Kv＝0.04Rc+0.4和Rc代入计算BQ值；

3.5：根据地下水、软弱结构面产状和天然应力影响对BQ进行修正，修正公式如下：

[BQ]＝BQ-100(K₁+K₂+K₃) (16)

式中：K1为地下水影响修正系数；K2为软件结构面产状影响修正系数；K3为天然应力影响修正系数；

3.6：得到修正后的BQ岩体分级结果，为III级岩体；

4)岩体三维裂隙网络模型生成和剖切，过程如下；

4.1：随机数求解，过程如下；

4.1.1：产生随机数的数学方法应满足以下条件：产生的随机数列应均匀分布在(0，1)区间；序列之间应无相关性；随机序列有足够长的重复周期，在计算机上产生的速度快，占有的内存空间小，具有完全可重复性；

4.1.2：利用Monte Carlo方法是根据确立的结构面几何概率模型，再现出服从这种模型的结构面网络模型，在(0，1)区间上生成均匀分布随机变量，利用这些均匀随机变量产生服从其他分布的随机数；

4.1.3：节理几何参数的密度函数有正态分布、对数正态分布、负指数分布、均匀分布四种；

4.1.4：根据求得的随机数，确定用于生成节理的基本几何参数；

4.2：岩体三维裂隙网络模型生成，过程如下；

4.2.1：根据结构面数据自动统计结果和求得的随机数，将每组结构面的数据保存到一个文本文件中，用st.dat表示；

4.2.2：st.dat数据内容格式依次为：每条结构面圆盘中心点坐标(x，y，z)，圆盘半径D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ和节理分组；

4.2.3：为区分不同组别的结构面，对相同组的结构面圆盘赋予相同的颜色，用数列GID表示；

4.2.4：利用Matlab软件，编写程序，读取结构面数据文件st.dat，运行后自动生成岩体三维裂隙网络模型；

4.2.5：得到岩体三维裂隙网络模型；

4.3：二维裂隙网络模型切割，过程如下；

4.3.1：在三维裂隙网络模型上，结合Matlab软件编程工具，以三维裂隙网络模型中心点为中心，实现任意角度的剖面切割功能；

4.3.2：得到穿过中心点的任意角度的二维裂隙网络模型；

4.3.3：在三维裂隙网络模型上，结合Matlab软件编程工具，在三维裂隙网络模型任意位置上，实现任意角度和方位的剖面切割功能；

4.3.4：得到任意角度和方位的二维裂隙网络模型，如图2所示；

4.3.5：将切割剖面上的数据，保存到st1.dat文件中，此时剖面处于三维坐标系下，文件中数据格式自左向右依次为：节理中心点坐标(x，y，z)，节理长度D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ；

4.3.6：将三维坐标系转化成二维坐标系，并将二维剖面数据保存到st2.dat文件中，文件中数据格式自左向右依次为：节理中心点坐标(x，y)，节理长度D，倾角DA，倾向DD，厚度thin，法向方向NX，NY，NZ；

4.4：输出三维裂隙网络模型的数据和任意二维裂隙网络模型的数据；

5)RQD_t各向异性图绘制，过程如下：

5.1：RQD_t求解计算，过程如下：

5.1.1：RQD_t理论公式如下：

式中：x_i表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，RQD_t代表对应阈值t的岩石质量指标，即阈值t下的RQD_t值；

5.1.2确定二维裂隙网络模型的剖面中心点O、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，则中心点O坐标为：

边界方程为：

5.1.3：以O为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度L等于O点到裂隙网络模型边界的距离，用L0～L35表示，则测线方程为：

式中：s表示测线，α表示角度；

5.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(x_a，y_a)，将测线方程与边界方程依次连立，判断测线是否与边界相交，原理如下：

5.1.5：求出测线与裂隙网络模型的边界方程交点(X_Z，α，Y_Z，α)；

5.1.6：确定测线所在区间，原理如下：

5.1.7：根据裂隙网络模型中每条节理的起点坐标(X_b，Y_b)和终点坐标(X_c，Y_c)，建立相应的解析方程，定义节理方程如下：

5.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(X_j，Y_j)范围，如果交点符合 a＜X_j＜c且b＜Y_j＜d，则记录该交点，遍历所有节理方程求出所有m个交点；

5.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；

5.1.10：计算相邻交点的距离，公式如下：

x₀＝a/2 (26)

y₀＝b/2 (27)

x_m+1＝X_Z，α (28)

y_m+1＝Y_Z，α (29)

5.1.11：输入一个阈值t；

5.1.12：循环比较d_i和t的大小，设初始l_t＝0，规则如下：

l_t＝l_t+d_i，若d_i＞t (30)

5.1.13：求解每条测线对应的RQD_t值，用m_α表示，以及对应的测线起点与终点的距离l_α，公式如下：

5.1.14：循环求出每条测线所对应的m_α，获得36个方向上的RQD_t值；

5.2：RQD_t各向异性图绘制，过程如下：

5.2.1：将36个RQD_t值，按角度顺序依次排序；

5.2.2：以O点为圆心，以l为半径画圆，在射线角度为α上找到与圆心距离为

的点标出来；

5.2.3：依次连接36组射线的端点，若某一条射线上的RQD_t为0，则取圆心；

5.2.4：绘制出RQD_t的各向异性图；

5.2.5：输出RQD_t各向异性图；

5.3：不同阈值t下RQD_t各向异性图绘制，过程如下：

5.3.1：输入不同的阈值t，求解出对应的RQD_t值；

5.3.2：将不同阈值t下的RQD_t值，以角度为变量，以RQD_t值为函数，绘制到同个坐标系下；

5.3.3：得到不同阈值t下的RQD_t各向异性图；

3.3.4：输出不同阈值t下的RQD_t各向异性图，如图3所示；

5.4：输出RQD_t的值；

6)基于BQ反演的最佳阈值t求解方法，过程如下：

6.1：基于BQ指标反演RQD_t范围，过程如下：

结合BQ分级计算出的岩体质量级别，岩体为III级岩体，查找《岩石质量指标》表，见表1，反演确定出该岩体级别下，RQD_t范围值，RQD_t范围在50％～75％

表1；

6.2：最佳阈值t求解方法，过程如下：

6.2.1：在三维裂隙网络模型上，过中心点O，以任意角度剖切三个剖面，得到三个二维裂隙网络模型，导出二维裂隙网络模型及数据；

6.2.2：针对每一个二维裂隙网络模型，设置不同的阈值t，求解出不同阈值t下的RQD_t值；

6.2.3：导出不同阈值t下的RQD_t值，计算出RQD_t均值；

6.2.4：以阈值t为横坐标，以RQD_t的均值为纵坐标，绘制RQD_t随阈值t变化的散点图；

6.2.5：根据散点图，设置拟合方程y＝aexp(bx)，拟合RQD_t随阈值t变化的曲线图，如图4所示；

6.2.6：将反演出的RQD_t范围值，带入到拟合出的RQD_t随阈值t变化的曲线图中，结合函数方程和曲线图，求解出在该RQD_t范围内，阈值t的范围，共得到三组阈值t的范围，见表2

表2；

6.2.7：针对三组阈值t的范围，取其范围的交集，以0.076m～0.124m作为最佳阈值t的范围；

6.2.8：以最佳阈值t范围的中点值作为最佳阈值t值，得到最佳阈值t，最佳阈值t＝0.1m；

6.2.9：输出最佳阈值t的范围和最佳阈值t值；

7)RQD_t各向异性求解方法，过程如下：

7.1：基于最佳阈值t和剖切的三个二维裂隙网络模型，求解出最佳阈值t下的RQD_t值，每个裂隙网络模型求解出36个RQD_t值；

7.2：计算出每个裂隙网络模型上RQD_t的最大值RQD_tmax、最小值RQD_tmin和均值

7.3：根据研究方向和RQD_tmin方向位置关系，结合方差进行修正，提出各向异性条件下RQD_t的计算公式如下：

式中：

为RQD_t的均值，RQD′_t为RQD_t的修正系数；

7.4：修正系数RQD′_t的求解考虑研究方向在某一具体方位角下的RQD_tθ值和方差D，按如下方法修正，当RQD_tθ＝RQD_tmin时，RQD′_t＝-D，D为方差；当RQD_tθ＝RQD_tmax时，RQD′_t＝D；当

时，RQD′_t＝0；

7.5：根据修正系数RQD′_t的求解过程，提出RQD′_t修正公式如下：

式中：RQD_tθ为θ角下的RQD_t值，RQD_tmin为阈值t下的RQD_t最小值，a，b为相关系数；

7.6：提出RQD_t各向异性计算公式如下：

7.7：相关系数a和b的求解，过程如下：

7.7.1：计算出每个裂隙网络模型上RQD_t的最大值RQD_tmax、最小值RQD_tmin和均值

7.7.2：根据RQD_tθ与RQD_tmax、RQD_tmin和

以及方差的关系，求解出三组RQD′_t值，并基于三组RQD′_t值，绘制出散点图；

7.7.3：根据散点图做出拟合曲线，如图5所示，曲线截距即为a值，斜率即为b值，曲线方程即为RQD′_t的修正公式，得到的RQD′_t公式如下：

7.8：将求解出的a、b值，带入到修正系数公式及RQD_t各向异性计算公式中，得到与角度θ有关的RQD_t各向异性公式：

7.9：根据RQD_t各向异性公式，求解出任意角度的RQD_t值；

8)Mathews稳定图法的改进方法，过程如下：

8.1：Mathews稳定图法的设计过程，过程如下：

8.1.1：Mathews稳定图方法的设计过程是以两个因子-稳定数N和水力半径R的计算为基础，然后将这两个因子绘制在划分为预测稳定区、潜在不稳定区和崩落区的图上，稳定数N代表岩体在给定应力条件下维持稳定的能力，水力半径R反映了采空区尺寸和形状，Mathews稳定图法计算公式为：

N＝Q′ABC (38)

式中：X、Y为矿房上下盘围岩宽度和长度；Q′为根据勘测图或钻孔岩芯记录计算出的结果；A为岩石应力系数；B为节理产状调整系数；C为重力调整系数；

8.1.2：岩石应力系数A由完整岩石单轴抗压强度与采场中线采矿产生的地应力的比值确定，公式如下：

当

当

当

8.1.3：节理产状调整系数B的值是通过采场面倾角与主要节理组的倾角之差来度量；

8.1.4：重力调整系数C反映了在重力的影响下采场面产状对采场矿岩稳定性的影响，其大小取决于采场顶板暴露表面的崩落、滑落以及边帮的滑落等，重力调整系数C和采场表面倾角α的关系由下式确定：

C＝8-7cosα (43)

8.2：Mathews稳定图法的改进方法，过程如下：

8.2.1：在求解Q′值时，用RQD_t值替代钻孔岩芯记录的RQD值，并将RQD_t的各向异性特征纳入考虑范畴，根据RQD_t各向异性求解公式，得到考虑RQD_t各向异性特征的Q′值，公式如下：

式中：RQD_tθ为角度θ下的RQD_t值；RQD_tmin为阈值t下的RQD最小值；

为RQD_t的均值；a，b为相关拟合系数；J_n为节理组数；J_r为节理粗糙度系数；J_a为节理蚀变系数；J_w为节理水折减系数；SRF 为应力折减系数；

8.2.2：结合等概率轮廓图对Mathews稳定图进行重新绘制，稳定区-破坏边界的曲线公式如下：

logN＝1.8206logR+1.618 (45)

破坏-崩落边界的曲线公式如下：

logN＝1.8076logR-3 (46)

8.2.3：改进后Mathews稳定性图分为三个区域：稳定区、破坏或主要破坏区和崩落区，处在稳定-破坏边界线上的工程，采场57％概率稳定，43％概率破坏，0％概率崩落；处在崩落-破坏边界线上的工程，采场0％概率稳定，5％概率破坏，95％概率崩落；处在稳定区的工程是稳定的；处于破坏或主要破坏区的工程，57％～0％的概率稳定，43％～5％的概率破坏，0％～95％的概率崩落；处在崩落区的工程，将持续发生崩落；

9)改进Mathews稳定图评价方法，过程如下：

9.1：选取用于开展上下盘围岩稳定性评价的矿房；

9.2：测量出矿房上下盘围岩宽度X和长度Y，计算出岩石应力系数A；节理产状调整系数B和重力调整系数C；

9.3：测量并求解出矿房上下盘的节理组数J_n；节理粗糙度系数J_r；节理蚀变系数J_a；节理水折减系数J_w和SRF为应力折减系数；

9.4：求解出最佳阈值t下的RQD_tθ、RQD_tmin和

计算出相关拟合系数a和b，求解出考虑RQD_t各向异性特征的Q′值；

9.5：计算出矿房稳定数N和水力半径R；

9.6：计算出稳定区-破坏边界的曲线公式和破坏-崩落边界的曲线公式，并绘制到改进后的Mathews 稳定性图上；

9.7：将计算出的矿房稳定数N和水力半径R绘制到改进后的Mathews稳定图上，如图6所示，根据绘制的散点所处的区域，给出矿房上下盘围岩的稳定性评价结果，3203矿房和3204矿房处于崩落区，矿房围岩评价结果为崩落，3205矿房处于破坏或主要破坏区，围岩破坏的概率为55％，围岩崩落的概率为 33％，围岩稳定的概率为12％；

10)围岩稳定性分析方法，过程如下：

10.1：根据矿房平面图和剖面图，利用CAD软件，建立矿房三维模型；

10.2：将三维模型导入到有限元分析软件中，赋予三维模型材料参数，包括岩体强度、密度、弹性模量、泊松比、粘聚力、摩擦角和容重；

10.3：赋予三维模型边界条件和初始应力；

10.4：进行矿房回采的数值模拟，首先进行未回采状态的数值模拟，得到初始应力场，然后进行矿房回采的数值模拟，得到采空区的应力场；

10.5：对矿房回采后的模拟结果进行分析，分析区域为矿房顶板、上下盘围岩和矿柱；

10.6：岩体剪切破坏的判断准则选用以主应力形式表示的Mohr-Coulomb屈服准则，当f_s＞0时，岩体发生剪切破坏，公式如下：

式中：σ₁、σ₃分别为最大主应力和最小主应力；c为粘聚力；

为摩擦角；

10.7：岩体张拉破坏的判断准则，根据岩体所受拉应力σ_t与抗拉强度R_t之间的关系确定，当σ_t≥R_t时， F≥0，即认为岩体发生张拉破坏，公式如下：

F＝σ_t-R_t (48)

式中：σ_t为岩体所受拉应力；R_t为岩体抗拉强度；

10.8：根据分析结果，得到矿房顶板、上下盘围岩和矿柱的破坏模式；

11)空区稳定性评价方法，过程如下：

11.1：以改进Mathews稳定图法得到的围岩稳定性评价结果作为采空区稳定性的初始评价结果；

11.2：利用数值模拟得到的围岩稳定性分析结果，对初始评价结果进行修正，过程如下：

11.2.1：当改进Mathews稳定图法评价的采空区围岩处于崩落区时，则该采空区处于崩落状态；

11.2.2：当改进Mathews稳定图法评价的采空区围岩处于破坏区时，则该采空区处于破坏状态；

11.2.3：当改进Mathews稳定图法评价的采空区围岩处于稳定区时，若数值模拟得到的采空区顶板、间柱、上下盘围岩未发生破坏，则该采空区处于稳定状态，若数值模拟得到的采空区顶板、间柱、上下盘围岩发生破坏，则该采空区有可能处于破坏状态；

11.3：得到采空区的稳定性评价结果。

Claims (9)

1.一种基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述方法包括以下步骤：

(1)结构面三维激光扫描快速获取，过程如下：

1.1：根据扫描目标和场地条件，选择扫描机位点，架设三脚架，架设中要确保仪器按照一定的扫描路线可以完整的获取边坡岩体的三维空间点云信息，同时要尽可能保证三脚架台面水平，并放置控制靶；

1.2：放置扫描仪主机在三脚架台面，固定旋钮，通过粗调脚架及微调扫描仪底座使主机气泡居中，设置扫描仪端口参数；

1.3：启动扫描控制软件，配置扫描仪相关参数，进入扫描仪控制界面，规划扫描角度，根据扫描目标设置扫描范围，调整相机配置参数，获取扫描目标图像；

1.4：固定扫描范围，获取扫描间距，设定采样间距，开始数据获取，并实时查看扫描点云数据及彩色信息情况，根据扫描成果随时调整扫描参数设定；

1.5：导出结构面点云数据；

(2)结构面聚类分析；

(3)基于BQ指标的岩体质量计算，过程如下：

3.1：根据结构面参数计算岩体完整性系数，公式如下：

式中：Jv为岩体体积节理数，单位条/m³；

3.2：Jv计算公式如下：

式中：L1，L2，...，Ln为垂直于结构面测线长度；N1，N2，...，Nn为同组结构面数目；

3.3：根据岩石单轴抗压强度值和岩体完整性系数值，计算BQ值：

BQ＝90+3R_c+250K_v

式中：Rc是岩石单轴抗压强度；Kv为岩体完整性系数；

3.4：在应用BQ计算公式过程中，遵循以下条件：

当Rc＞90Kv+30时，以Rc＝90Kv+30和Kv代入计算BQ值；

当Kv＞0.04Rc+0.4时，以Kv＝0.04Rc+0.4和Rc代入计算BQ值；

3.5：根据地下水、软弱结构面产状和天然应力影响对BQ进行修正，修正公式如下：

[BQ]＝BQ-100(K₁+K₂+K₃)

式中：K1为地下水影响修正系数；K2为软件结构面产状影响修正系数；K3为天然应力影响修正系数；

3.6：得到修正后的BQ岩体分级结果；

(4)岩体三维裂隙网络模型生成和剖切；

(5)RQD_t各向异性图绘制；

(6)基于BQ反演的最佳阈值t求解方法；

(7)RQD_t各向异性求解方法；

(8)Mathews稳定图法的改进方法；

(9)改进Mathews稳定图评价方法；

(10)围岩稳定性分析方法，过程如下：

10.1：根据矿房平面图和剖面图，利用CAD软件，建立矿房三维模型；

10.2：将三维模型导入到有限元分析软件中，赋予三维模型材料参数，包括岩体强度、密度、弹性模量、泊松比、粘聚力、摩擦角和容重；

10.3：赋予三维模型边界条件和初始应力；

10.4：进行矿房回采的数值模拟，首先进行未回采状态的数值模拟，得到初始应力场，然后进行矿房回采的数值模拟，得到采空区的应力场；

10.5：对矿房回采后的模拟结果进行分析，分析区域为矿房顶板、上下盘围岩和矿柱；

10.6：岩体剪切破坏的判断准则选用以主应力形式表示的Mohr-Coulomb屈服准则，当f_s＞0时，岩体发生剪切破坏，公式如下：

式中：σ₁、σ₃分别为最大主应力和最小主应力；c为粘聚力；

为摩擦角；

10.7：岩体张拉破坏的判断准则，根据岩体所受拉应力σ_t与抗拉强度R_t之间的关系确定，当σ_t≥R_t时，F≥0，即认为岩体发生张拉破坏，公式如下：

F＝σ_t-R_t

式中：σ_t为岩体所受拉应力；R_t为岩体抗拉强度；

10.8：根据分析结果，得到矿房顶板、上下盘围岩和矿柱的破坏模式；

(11)空区稳定性评价方法。

2.如权利要求1所述的基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述步骤(11)中，空区稳定性评价方法的过程如下：

11.1：以改进Mathews稳定图法得到的围岩稳定性评价结果作为采空区稳定性的初始评价结果；

11.2：利用数值模拟得到的围岩稳定性分析结果，对初始评价结果进行修正，过程如下：

11.2.1：当改进Mathews稳定图法评价的采空区围岩处于崩落区时，则该采空区处于崩落状态；

11.2.2：当改进Mathews稳定图法评价的采空区围岩处于破坏区时，则该采空区处于破坏状态；

11.2.3：当改进Mathews稳定图法评价的采空区围岩处于稳定区时，若数值模拟得到的采空区顶板、间柱、上下盘围岩未发生破坏，则该采空区处于稳定状态，若数值模拟得到的采空区顶板、间柱、上下盘围岩发生破坏，则该采空区有可能处于破坏状态；

11.3：得到采空区的稳定性评价结果。

3.如权利要求1所述的基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述步骤(9)中，改进Mathews稳定图评价方法的过程如下：

9.1：选取用于开展上下盘围岩稳定性评价的矿房；

9.2：测量出矿房上下盘围岩宽度X和长度Y，计算出岩石应力系数A；节理产状调整系数B和重力调整系数C；

9.3：测量并求解出矿房上下盘的节理组数J_n；节理粗糙度系数J_r；节理蚀变系数J_a；节理水折减系数J_w和SRF为应力折减系数；

9.4：求解出最佳阈值t下的RQD_tθ、RQD_tmin和

计算出相关拟合系数a和b，求解出考虑RQD_t各向异性特征的Q′值；

9.5：计算出矿房稳定数N和水力半径R；

9.6：计算出稳定区-破坏边界的曲线公式和破坏-崩落边界的曲线公式，并绘制到改进后的Mathews稳定性图上；

9.7：将计算出的矿房稳定数N和水力半径R绘制到改进后的Mathews稳定图上，根据绘制的散点所处的区域，给出矿房上下盘围岩的稳定性评价结果，若散点处于稳定区，则矿房围岩评价结果为稳定，若散点处于破坏或主要破坏区，需再根据矿房围岩破坏和崩落的概率，给出矿房围岩的破坏概率和崩落概率，若散点处于崩落区，则矿房围岩评价结果为崩落。

4.如权利要求1所述的基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述步骤(8)中，Mathews稳定图法的改进方法的过程如下：

8.1：Mathews稳定图法的设计过程，过程如下：

8.1.1：Mathews稳定图方法的设计过程是以两个因子-稳定数N和水力半径R的计算为基础，然后将这两个因子绘制在划分为预测稳定区、潜在不稳定区和崩落区的图上，稳定数N代表岩体在给定应力条件下维持稳定的能力，水力半径R反映了采空区尺寸和形状，Mathews稳定图法计算公式为：

N＝Q′ABC

式中：X、Y为矿房上下盘围岩宽度和长度；Q′为根据勘测图或钻孔岩芯记录计算出的结果；A为岩石应力系数；B为节理产状调整系数；C为重力调整系数；

8.1.2：岩石应力系数A由完整岩石单轴抗压强度与采场中线采矿产生的地应力的比值确定，公式如下：

当

A＝0.1

当

当

A＝1

8.1.3：节理产状调整系数B的值是通过采场面倾角与主要节理组的倾角之差来度量；

8.1.4：重力调整系数C反映了在重力的影响下采场面产状对采场矿岩稳定性的影响，其大小取决于采场顶板暴露表面的崩落、滑落以及边帮的滑落等，重力调整系数C和采场表面倾角α的关系由下式确定：

C＝8-7cosα

8.2：Mathews稳定图法的改进方法，过程如下：

8.2.1：在求解Q′值时，用RQD_t值替代钻孔岩芯记录的RQD值，并将RQD_t的各向异性特征纳入考虑范畴，根据RQD_t各向异性求解公式，得到考虑RQD_t各向异性特征的Q′值，公式如下：

式中：RQD_tθ为角度θ下的RQD_t值；RQD_tmin为阈值t下的RQD最小值；

为RQD_t的均值；a，b为相关拟合系数；J_n为节理组数；J_r为节理粗糙度系数；J_a为节理蚀变系数；J_w为节理水折减系数；SRF为应力折减系数；

8.2.2：结合等概率轮廓图对Mathews稳定图进行重新绘制，稳定区-破坏边界的曲线公式如下：

logN＝1.8206logR+1.618

破坏-崩落边界的曲线公式如下：

logN＝1.8076logR-3

8.2.3：改进后Mathews稳定性图分为三个区域：稳定区、破坏或主要破坏区和崩落区，处在稳定-破坏边界线上的工程，采场57％概率稳定，43％概率破坏，0％概率崩落；处在崩落-破坏边界线上的工程，采场0％概率稳定，5％概率破坏，95％概率崩落；处在稳定区的工程是稳定的；处于破坏或主要破坏区的工程，57％～0％的概率稳定，43％～5％的概率破坏，0％～95％的概率崩落；处在崩落区的工程，将持续发生崩落。

5.如权利要求1所述的基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述步骤(7)中，RQD_t各向异性求解方法的过程如下：

7.1：基于最佳阈值t和剖切的三个二维裂隙网络模型，求解出最佳阈值t下的RQD_t值，每个裂隙网络模型求解出36个RQD_t值；

7.2：计算出每个裂隙网络模型上RQD_t的最大值RQD_tmax、最小值RQD_tmin和均值

7.3：根据研究方向和RQD_tmin方向位置关系，结合方差进行修正，提出各向异性条件下RQD_t的计算公式如下：

式中：

为RQD_t的均值，RQD′_t为RQD_t的修正系数；

7.4：修正系数RQD′_t的求解考虑研究方向在某一具体方位角下的RQD_tθ值和方差D，按如下方法修正，当RQD_tθ＝RQD_tmin时，RQD′_t＝-D，D为方差；当RQD_tθ＝RQD_tmax时，RQD′_t＝D；当

时，RQD′_t＝0；

7.5：根据修正系数RQD′_t的求解过程，提出RQD′_t修正公式如下：

式中：RQD_tθ为θ角下的RQD_t值，RQD_tmin为阈值t下的RQD_t最小值，a，b为相关系数；

7.6：提出RQD_t各向异性计算公式如下：

7.7：相关系数a和b的求解，过程如下：

7.7.1：计算出每个裂隙网络模型上RQD_t的最大值RQD_tmax、最小值RQD_tmin和均值

7.7.2：根据RQD_tθ与RQD_tmax、RQD_tmin和

以及方差的关系，求解出三组RQD′_t值，并基于三组RQD′_t值，绘制出散点图；

7.7.3：根据散点图做出拟合曲线，曲线截距即为a值，斜率即为b值，曲线方程即为RQD′_t的修正公式；

7.8：将求解出的a、b值，带入到修正系数公式及RQD_t各向异性计算公式中，得到与角度θ有关的RQD_t各向异性公式；

7.9：根据RQD_t各向异性公式，求解出任意角度的RQD_t值。

6.如权利要求1所述的基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述步骤(6)中，基于BQ反演的最佳阈值t求解方法的过程如下：

6.1：基于BQ指标反演RQD_t范围，过程如下：

结合BQ分级计算出的岩体质量级别，查找《岩石质量指标》表，反演确定出该岩体级别下，RQD_t范围值；

6.2：最佳阈值t求解方法，过程如下：

6.2.1：在三维裂隙网络模型上，过中心点O，以任意角度剖切三个剖面，得到三个二维裂隙网络模型，导出二维裂隙网络模型及数据；

6.2.2：针对每一个二维裂隙网络模型，设置不同的阈值t，求解出不同阈值t下的RQD_t值；

6.2.3：导出不同阈值t下的RQD_t值，计算出RQD_t均值；

6.2.4：以阈值t为横坐标，以RQD_t的均值为纵坐标，绘制RQD_t随阈值t变化的散点图；

6.2.5：根据散点图，设置拟合方程，拟合RQD_t随阈值t变化的曲线图；

6.2.6：将反演出的RQD_t范围值，带入到拟合出的RQD_t随阈值t变化的曲线图中，结合函数方程和曲线图，求解出在该RQD_t范围内阈值t的范围，共得到三组阈值t的范围；

6.2.7：针对三组阈值t的范围，取其范围的交集，作为最佳阈值t的范围；

6.2.8：以最佳阈值t范围的中点值作为最佳阈值t值，得到最佳阈值t；

6.2.9：输出最佳阈值t的范围和最佳阈值t值。

7.如权利要求1所述的基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述步骤(5)中，RQD_t各向异性图绘制的过程如下：

5.1：RQD_t求解计算，过程如下：

5.1.1：RQD_t理论公式如下：

式中：x_i表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，RQD_t代表对应阈值t的岩石质量指标，即阈值t下的RQD_t值；

5.1.2确定二维裂隙网络模型的剖面中心点O、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，求解出中心点O坐标(X₀，Y₀)和边界方程；

5.1.3：以O为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度L等于O点到裂隙网络模型边界的距离，用L0～L35表示，求解出测线方程；

5.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(x_a，y_a)，将测线方程与边界方程依次连立，判断测线是否与边界相交；

5.1.5：求出测线与裂隙网络模型的边界方程交点(X_Z，α，Y_Z，α)；

5.1.6：确定测线所在区间，原理如下：

5.1.7：根据裂隙网络模型中每条节理的起点坐标(X_b，Y_b)和终点坐标(X_c，Y_c)，建立相应的解析方程；

5.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(X_j，Y_j)范围，如果交点符合a＜X_j＜c且b＜Y_j＜d，则记录该交点，遍历所有节理方程求出所有m个交点；

5.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；

5.1.10：计算相邻交点的距离d_i；

5.1.11：输入一个阈值t；

5.1.12：循环比较d_i和t的大小，设初始l_t＝0，规则如下：

l_t＝l_t+d_i，若d_i＞t

5.1.13：求解每条测线对应的RQD_t值，用mα表示，以及对应的测线起点与终点的距离l_α；

5.1.14：循环求出每条测线所对应的m_α，获得36个方向上的RQD_t值；

5.2：RQD_t各向异性图绘制，过程如下：

5.2.1：将36个RQD_t值，按角度顺序依次排序；

5.2.2：以O点为圆心，以l为半径画圆，在射线角度为α上找到与圆心距离为

的点标出来；

5.2.3：依次连接36组射线的端点，若某一条射线上的RQD_t为0，则取圆心；

5.2.4：绘制出RQD_t的各向异性图；

5.2.5：输出RQD_t各向异性图；

5.3：不同阈值t下RQD_t各向异性图绘制，过程如下：

5.3.1：输入不同的阈值t，求解出对应的RQD_t值；

5.3.2：将不同阈值t下的RQD_t值，以角度为变量，以RQD_t值为函数，绘制到同个坐标系下；

5.3.3：得到不同阈值t下的RQD_t各向异性图；

5.3.4：输出不同阈值t下的RQD_t各向异性图；

5.4：输出RQD_t的值。

8.如权利要求1所述的基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述步骤(4)中，岩体三维裂隙网络模型生成和剖切的过程如下：

4.1：随机数求解，过程如下：

4.1.1：产生随机数的数学方法应满足以下条件：产生的随机数列应均匀分布在(0，1)区间；序列之间应无相关性；随机序列有足够长的重复周期，在计算机上产生的速度快，占有的内存空间小，具有完全可重复性；

4.1.2：利用Monte Carlo方法是根据确立的结构面几何概率模型，再现出服从这种模型的结构面网络模型，在(0，1)区间上生成均匀分布随机变量，利用这些均匀随机变量产生服从其他分布的随机数；

4.1.3：节理几何参数的密度函数有正态分布、对数正态分布、负指数分布、均匀分布四种；

4.1.4：根据求得的随机数，确定用于生成节理的基本几何参数；

4.2：岩体三维裂隙网络模型生成，过程如下：

4.2.1：根据结构面数据自动统计结果和求得的随机数，将每组结构面的数据保存到一个文本文件中，用st.dat表示；

4.2.2：st.dat数据内容格式依次为：每条结构面圆盘中心点坐标(x，y，z)，圆盘半径D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ和节理分组；

4.2.3：为区分不同组别的结构面，对相同组的结构面圆盘赋予相同的颜色，用数列GID表示；

4.2.4：利用Matlab软件，编写程序，读取结构面数据文件st.dat，运行后自动生成岩体三维裂隙网络模型；

4.2.5：得到岩体三维裂隙网络模型；

4.3：二维裂隙网络模型切割，过程如下：

4.3.1：在三维裂隙网络模型上，结合Matlab软件编程工具，以三维裂隙网络模型中心点为中心，实现任意角度的剖面切割功能；

4.3.2：得到穿过中心点的任意角度的二维裂隙网络模型；

4.3.3：在三维裂隙网络模型上，结合Matlab软件编程工具，在三维裂隙网络模型任意位置上，实现任意角度和方位的剖面切割功能；

4.3.4：得到任意角度和方位的二维裂隙网络模型；

4.3.5：将切割剖面上的数据，保存到st1.dat文件中，此时剖面处于三维坐标系下，文件中数据格式自左向右依次为：节理中心点坐标(x，y，z)，节理长度D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ；

4.3.6：将三维坐标系转化成二维坐标系，并将二维剖面数据保存到st2.dat文件中，文件中数据格式自左向右依次为：节理中心点坐标(x，y)，节理长度D，倾角DA，倾向DD，厚度thin，法向方向NX，NY，NZ；

4.4：输出三维裂隙网络模型的数据和任意二维裂隙网络模型的数据。

9.如权利要求1所述的基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法，其特征在于，所述步骤(2)中，结构面聚类分析的过程如下：

2.1：点云数据处理，过程如下：

2.1.1：导入激光扫描获得的结构面点云数据；

2.1.2：计算拓扑构造后的点云中当前点与相邻点的距离与距离均值，通过距离阈值对点云数据中噪声点进行识别、剔除；

2.1.3：根据三维激光扫描仪自身空间坐标位置和现场结构面产状方位，确定点云数据的空间三维坐标；

2.1.4：基于下半球等角度投影方法进行点云数据的转换，将以倾向和倾角表示的节理产状数据转换为以节理单位法向量表示的结构面产状数据；

2.1.5：得到以单位法向量表示的结构面数据；

2.2：近邻传播算法聚类分析，过程如下：

2.2.1：设结构面的实测样本数量为N，每个样本数据的倾向为X_i，倾角为Y_i，以每个样本数据的倾向X_i，倾角Y_i作为一个聚类，确定一个初始聚类中心，共得到N个初始聚类中心；

2.2.2：通过相似性度量准则，遍历所有样本数据，计算每个样本数据距离聚类中心的距离，并将每个样本数据分配到距离它最近的聚类中心，得到N组数据；

2.2.3：对于每组数据，通过特征模量分析方法，求解计算每组数据的聚类中心，假设某组内存在l个数据，首先，按如下公式计算矩阵S：

式中：(x_i，y_i，z_i)为任意结构面的单位法向量，i∈(1，l)；

2.2.4：其次，求解矩阵S的特征值(τ₁，τ₂，τ₃)和特征向量(ξ₁，ξ₂，ξ₃)，其中τ₁＜τ₂＜τ₃，最大特征值对应的特征向量ξ₃为组内l个向量的平均向量，将ξ₃作为新的聚类中心；

2.2.5：针对所有样本数据，重复计算每个样本数据距离聚类中心的距离、矩阵S以及特征值和特征向量，直到所有聚类中心的位置都固定，确定出结构面的分组；

2.2.6：将以单位法向量表示的结构面产状数据转换为以倾向、倾角表示的结构面产状数据；

2.2.7：对结构面产状数据进行统计分析，计算结构面倾角的平均值m与标准差σ，计算倾角数据的稳健区间[m-σ，m+σ]；

2.2.8：判断样本数据的初始聚类中心的倾向X_i和倾角Y_i是否落在稳健区间稳健区间[m-σ，m+σ]，若是，则聚类分析完成；若不是，则需要对样本数据重新聚类，直到初始聚类中心的倾向X_i和倾角Y_i均落在稳健区间内[m-σ，m+σ]；

2.3：极点图绘制，过程如下；

2.3.1：基于结构面法向产状数据，根据结构面空间赤平投影图的纵剖面原理，求解出所有结构面法线的赤平投影坐标点(x_n，y_n)；

2.3.3：绘制一条直径为单位长度的基圆，绘制出铅直和水平两条直径，并标出E、S、W、N；

2.3.4：将所有结构面的赤平投影坐标(x_n，y_n)，绘制在基圆图上；

2.3.5：实现结构面极点图的绘制；

2.4：结构面统计分析，过程如下：

2.4.1：确定样本分区区间m；

2.4.2：求解样本极差

2.4.3：计算每个分区区间M_m：

2.4.4：确定样本落在每个分区区间里的概率，先利用计算机循环语言统计落在每一个区间的样本个数N_m，结合样本总数N，计算样本数概率P_m；

2.4.5：求解样本均值

2.4.6：求解样本方差S²，其中S为标准差；

2.4.7：根据概率P_m值，绘制出每组结构面的倾向、倾角、迹长、间距和断距的概率分布形态；

2.5：输出结构面产状的分组信息，包括每组结构面的倾向、倾角、迹长、间距和断距的均值和方差。

Images