CN112132406A - 一种基于激光扫描、BQ和RQDt的巷道开挖不利方位求解方法 - Google Patents

Abstract

一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，属于巷道开挖方位领域，其步骤包括：(1)结构面三维激光扫描快速获取；(2)结构面极点图和倾向玫瑰花图绘制；(3)基于BQ指标的岩体质量计算；(4)岩体三维裂隙网络模型生成和剖切；(5)RQD_t各向异性图绘制；(6)基于BQ反演的最佳阈值t求解方法；(7)空间RQD_t求解方法；(8)巷道开挖不利方位求解方法。本发明实现了RQD_t的最佳阈值t求解、RQD_t椭圆曲线和各向异性椭球体绘制，实现了最佳阈值t下的巷道开挖不利方位求解。本发明方法明确，所得最佳阈值t、椭圆曲线和各向异性椭球体能够准确反映出岩体的RQD_t各向异性特征，得到的巷道开挖不利方位可用于指导巷道开挖施工应用。

Description

一种基于激光扫描、BQ和RQDt的巷道开挖不利方位求解方法

技术领域

本发明涉及一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，特别的是本发明基于三维激光扫描、近邻传播算法、BQ指标、裂隙网络模型、广义RQD理论和圆锥曲线方程，反演计算出RQD_t最佳阈值t的范围和最佳阈值t值，得到了最佳阈值t下的RQD_t椭圆曲线和各向异性椭球体，得到了最佳阈值t下的巷道开挖不利方位，提供了一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，属于巷道开挖方位领域。

背景技术

在自然界中，各向异性无处不在。岩体表现出非常明显的各向异性，其性质随观察角度和测量方向的变化而变化。岩体各向异性变化的起源可归结于地质起源，如沉积岩中独特的层理、玄武岩中独具的天然孔洞、片麻岩中的片麻状构造和板岩中的劈理，以及贯穿这些岩体的断层等。除了固有的地质起源造成的岩体各向异性外，一些其它因素也会影响岩体的各向异性，如岩体中节理组的产状和分布形态，岩体中的节理会造成岩体的应力各向异性，影响岩体的质量。

岩体的外观和质量强烈的依赖于岩体的各向异性程度，岩体的稳定性是由岩体质量所决定的。一些地质现象也已经证明了这一点。例如，在巷道或高速公路隧道开挖过程中，会沿着主要节理组倾角出现非垂直性的变形破坏，或者隧道中的某一块岩体出现楔形块体掉落，进而造成整个岩体的失稳或者巷道的坍塌，这种失稳或坍塌对巷道开挖是非常不利的，不仅影响工程施工进度，更重要的是会造成严重的人员伤害或者财产安全。因此，在巷道开采的前期，如何准确的获取岩体质量的各向异性程度，以及如何依据岩体质量的各向异性程度，选择巷道的开采方向和开挖走向，具有重要的研究价值和现实指导意义。而目前国内外的研究中，从岩体质量各向异性角度考虑巷道开挖方向的研究，尚不多见。

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的求解计算方法，因此众多针对RQD_t空间效应的研究，主要还是集中在三维RQD_t的计算上，如开展不同阈值t下的三维RQD_t计算，或者在此基础上研究不同方向不同位置上的RQD_t值，或者建立三维空间云图。很少有学者针对最佳阈值t，开展最能反映岩体质量的RQD_t空间效应研究。因此，得到最佳阈值t下的空间RQD_t的求解方法具有重要的理论意义和研究意义。

鉴于此，本发明提出了一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法。

发明内容

为了实现巷道开挖的不利方位求解，本发明提供了一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法。本发明基于三维激光扫描、近邻传播算法、BQ指标、裂隙网络模型、广义RQD理论和圆锥曲线方程，反演计算出RQD_t最佳阈值t的范围和最佳阈值t值，得到了最佳阈值t下的RQD_t椭圆曲线和各向异性椭球体，得到了最佳阈值t下的巷道开挖不利方位，提供了一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法。

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

一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，所述方法包括以下步骤：

(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)巷道开挖不利方位求解方法。

进一步，所述步骤(8)中，巷道开挖不利方位求解方法的过程如下：

8.1：平面上巷道开挖不利方位求解方法，过程如下：

8.1.1：根据RQD_t椭圆曲线，求解出曲线上最大值RQD_tmax和最小值RQD_tmin的值；

8.1.2：求解出RQD_t的各向异性比值，各向异性比计算公式如下：

式中：RQD_tmax为椭圆曲线上的RQD_t最大值，RQD_tmin为椭圆曲线上的RQD_t最小值；

8.1.3：各向异性比越大，岩体的不利方位越明显；

8.1.4：判断巷道开挖方向与RQD_tmin方向的关系，当巷道轴线与RQD_tmin方向一致时，对围岩稳定非常不利；当巷道轴线与RQD_tmin方向垂直时，对围岩稳定有利；

8.1.5：得到巷道开挖的不利方位和最优方位；

8.2：空间上巷道开挖不利方位求解方法，过程如下：

8.2.1：根据RQD_t各向异性椭球体，求解出椭球体上最大值RQD_tmax和最小值RQD_tmin的值；

8.2.2：求解出RQD_t的空间各向异性比值，空间各向异性比计算公式如下：

式中：RQD_tmax为椭球体上的RQD_t最大值，RQD_tmin为椭球体上的RQD_t最小值；

8.2.3：各向异性比越大，岩体的不利方位越明显；

8.2.4：判断巷道开挖方向与RQD_tmin方向的关系，当巷道轴线与RQD_tmin方向一致时，对围岩稳定非常不利；当巷道轴线与RQD_tmin方向垂直时，对围岩稳定有利；

8.2.5：得到巷道开挖的空间不利方位和最优方位。

进一步，所述步骤(7)中，空间RQD_t求解方法的过程如下：

7.1：平面RQD_t求解方法

根据最佳阈值t下的36个方位的RQD_t值，求解出任意角度的平面RQD_t值，过程如下：

7.1.1：求解出最佳阈值t下的36个RQD_t值；

7.1.2：以RQD_t值沿x方向的分量RQD_x作为横坐标，沿y轴方向的分量RQD_y作为纵坐标，设定拟合圆锥曲线方程为：

式中：p₁，p₂，p₃为拟合系数；

7.1.3：基于最小二乘法，应用Matlab软件编写拟合程序；

7.1.4：拟合得到RQD_t的椭圆曲线；

7.1.5：通过该拟合椭圆曲线，可以得到平面上任意角度的RQD_t值；

7.2：空间RQD_t求解方法

根据拟合得到的RQD_t椭圆曲线，求解出空间内任意方位任意角度的RQD_t值，过程如下：

7.2.1：根据拟合得到的RQD_t椭圆曲线，求解出椭圆的长轴所在位置；

7.2.2：以椭圆曲线的长轴作为旋转轴，利用Matlab软件编写程序，将椭圆曲线沿长轴旋转360°，得到RQD_t的各向异性椭球体；

7.2.3：得到空间内任意方位任意角度的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求解计算

基于RQD_t理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的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

y₀＝b/2

x_m+1＝X_Z，α

y_m+1＝Y_Z，α

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各向异性图绘制

根据求解出的36个RQD_t值，绘制出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各向异性图绘制

用于绘制不同阈值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)区间上生成均匀分布随机变量ri，利用这些均匀随机变量产生服从其他分布的随机数；

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：裂隙网络模型数据输出

实现三维裂隙网络模型的数据输出和任意二维裂隙网络模型的数据输出，过程如下；

4.4.1：输出st.dat文件数据信息；

4.4.2：输出st1.dat文件数据信息；

4.4.3：输出st2.dat文件数据信息。

再进一步，所述步骤(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：数据输出

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

2.6：倾向玫瑰花图绘制，过程如下：

2.6.1：将节理倾向数据，按节理倾向方位角大小依次排序，每隔θ＝10°角度进行分组，每组命名D_j；

2.6.2：统计每组节理的数目

和每组节理的平均倾向

2.6.3：根据作图的大小和各组节理数目，选取一定长度的线段代表一组节理，确定线段的比例尺L_D；

2.6.4：以等于按比例尺L_D表示的、数目最多的一组节理

的线段长度

为半径作圆，过圆心作南北线及东西线，在圆周上标明方位角；

2.6.5：对每组节理D_j，按平均倾向

为方位角，在半圆上做出记号，自圆心向圆周记号点半径方向，按组内数目

和比例尺L_D定出一点

该点即表示该组节理平均倾向和节理数目；

2.6.6：顺次连接

和

如其中某组节理数为零，则连线回到圆心，再由圆心引出与下一组相连；

2.6.7：绘制出结构面倾向玫瑰花图。

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

1、本发明基于三维激光扫描、近邻传播算法、BQ指标、裂隙网络模型、广义RQD理论和圆锥曲线方程，提供了一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法；

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

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

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

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

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

7、本发明实现了RQD_t椭圆曲线和各向异性椭球体的求解绘制；

8、本发明实现了巷道开挖的不利方位求解和最优方位求解；

9、本发明方法明确，得到的最佳阈值t能够准确反映出岩体质量，得到的RQD_t椭圆曲线和各向异性椭球体能够准确反映出岩体的RQD_t各向异性特征，得到的巷道开挖不利方位和最优方位可用于指导巷道开挖施工应用。

附图说明

图1是结构面极点图。

图2是二维裂隙网络图。

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

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

图5是RQD_t椭圆曲线。

具体实施方式

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

参照图1～图5，一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，包括以下步骤：

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：数据输出

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

2.6：倾向玫瑰花图绘制

根据倾向玫瑰花图绘制方法，绘制出结构面倾向玫瑰花图，过程如下；

2.6.1：将节理倾向数据，按节理倾向方位角大小依次排序，每隔θ＝10°角度进行分组，每组命名D_j；

D_j＝{θ，θ+9°} (13)

θ＝10(j-1) (14)

j∈(1，36) (15)

2.6.2：统计每组节理的数目

和每组节理的平均倾向

j∈(1，36) (17)

2.6.3：根据作图的大小和各组节理数目，选取一定长度的线段代表一组节理，确定线段的比例尺L_D；

2.6.4：以等于按比例尺L_D表示的、数目最多的一组节理

的线段长度

为半径作圆，过圆心作南北线及东西线，在圆周上标明方位角；

N_D为常数 (18)

2.6.5：对每组节理D_j，按平均倾向

为方位角，在半圆上做出记号，自圆心向圆周记号点半径方向，按组内数目

和比例尺L_D定出一点

该点即表示该组节理平均倾向和节理数目；

2.6.6：顺次连接

和

如其中某组节理数为零，则连线回到圆心，再由圆心引出与下一组相连；

2.6.7：绘制出结构面倾向玫瑰花图；

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

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

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

3.2：Jv计算公式如下：

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

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

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

式中：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₃) (23)

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

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

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

4.1：随机数求解

基于蒙特卡洛方法和节理分布形态，求解随机数，过程如下；

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

4.1.2：给定分布下随机变量数值的产生，Monte Carlo方法是根据确立的结构面几何概率模型，再现服从这种模型的结构面网络模型，该方法的核心是抽样数的随机性，高品质随机数可以得到良好的模拟结果，即在(0，1)区间上生成均匀分布随机变量ri，利用这些均匀随机变量产生服从其他分布的随机数；

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：裂隙网络模型数据输出

实现三维裂隙网络模型的数据输出和任意二维裂隙网络模型的数据输出，过程如下；

4.4.1：输出st.dat文件数据信息；

4.4.2：输出st1.dat文件数据信息；

4.4.3：输出st2.dat文件数据信息；

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

5.1：RQD_t求解计算

基于RQD_t理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的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 (33)

y₀＝b/2 (34)

x_m+1＝X_Z，α (35)

y_m+1＝Y_Z，α (36)

5.1.11：输入一个阈值t；

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

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

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

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

5.2：RQD_t各向异性图绘制

根据求解出的36个RQD_t值，绘制出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各向异性图绘制

用于绘制不同阈值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范围在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：平面RQD_t求解方法

根据最佳阈值t下的36个方位的RQD_t值，求解出任意角度的平面RQD_t值，过程如下：

7.1.1：求解出最佳阈值t下的36个RQD_t值；

7.1.2：以RQD_t值沿x方向的分量RQD_x作为横坐标，沿y轴方向的分量RQD_y作为纵坐标，设定拟合圆锥曲线方程为：

式中：p₁，p₂，p₃为拟合系数；

7.1.3：基于最小二乘法，应用Matlab软件编写拟合程序；

7.1.4：拟合得到RQD_t的椭圆曲线，如图5所示；

7.1.5：通过该拟合椭圆曲线，可以得到平面上任意角度的RQD_t值；

7.2：空间RQD_t求解方法

根据拟合得到的RQD_t椭圆曲线，求解出空间内任意方位任意角度的RQD_t值，过程如下：

7.2.1：根据拟合得到的RQD_t椭圆曲线，求解出椭圆的长轴所在位置；

7.2.2：以椭圆曲线的长轴作为旋转轴，利用Matlab软件编写程序，将椭圆曲线沿长轴旋转360°，得到RQD_t的各向异性椭球体；

7.2.3：得到空间内任意方位任意角度的RQD_t值；

8)巷道开挖不利方位求解方法，过程如下：

8.1：平面上巷道开挖不利方位求解方法，过程如下：

8.1.1：根据RQD_t椭圆曲线，求解出曲线上最大值RQD_tmax和最小值RQD_tmin的值；

8.1.2：求解出RQD_t的各向异性比值，各向异性比计算公式如下：

式中：RQD_tmax为椭圆曲线上的RQD_t最大值，RQD_tmin为椭圆曲线上的RQD_t最小值；

8.1.3：各向异性比越大，岩体的不利方位越明显，计算得到的各向异性比为3.33；

8.1.4：判断巷道开挖方向与RQD_tmin方向的关系，当巷道轴线与RQD_tmin方向一致时，对围岩稳定非常不利；当巷道轴线与RQD_tmin方向垂直时，对围岩稳定有利；计算得到RQD_tmin的方向为230°，即当巷道开挖方向为230°时，为巷道开挖的不利方位，围岩稳定性处于最不利的状态，当巷道开挖方向为140°时，为巷道开挖的最优方位，围岩稳定性处于最优状态；

8.1.5：得到巷道开挖的不利方位和最优方位，巷道开挖的不利方位为230°，最优方位为140°；

8.2：空间上巷道开挖不利方位求解方法，过程如下：

8.2.1：根据RQD_t各向异性椭球体，求解出椭球体上最大值RQD_tmax和最小值RQD_tmin的值；

8.2.2：求解出RQD_t的空间各向异性比值，空间各向异性比计算公式如下：

式中：RQD_tmax为椭球体上的RQD_t最大值，RQD_tmin为椭球体上的RQD_t最小值；

8.2.3：各向异性比越大，岩体的不利方位越明显；

8.2.4：判断巷道开挖方向与RQD_tmin方向的关系，当巷道轴线与RQD_tmin方向一致时，对围岩稳定非常不利；当巷道轴线与RQD_tmin方向垂直时，对围岩稳定有利；

8.2.5：得到巷道开挖的空间不利方位和最优方位。

Claims (7)

1.一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，其特征在于，所述方法包括以下步骤：

(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)巷道开挖不利方位求解方法。

2.如权利要求1所述的一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，其特征在于，所述步骤(8)中，巷道开挖不利方位求解方法过程如下：

8.1：平面上巷道开挖不利方位求解方法，过程如下：

8.1.1：根据RQD_t椭圆曲线，求解出曲线上最大值RQD_tmax和最小值RQD_tmin的值；

8.1.2：求解出RQD_t的各向异性比值，各向异性比计算公式如下：

式中：RQD_tmax为椭圆曲线上的RQD_t最大值，RQD_tmin为椭圆曲线上的RQD_t最小值；

8.1.3：各向异性比越大，岩体的不利方位越明显；

8.1.4：判断巷道开挖方向与RQD_tmin方向的关系，当巷道轴线与RQD_tmin方向一致时，对围岩稳定非常不利；当巷道轴线与RQD_tmin方向垂直时，对围岩稳定有利；

8.1.5：得到巷道开挖的不利方位和最优方位；

8.2：空间上巷道开挖不利方位求解方法，过程如下：

8.2.1：根据RQD_t各向异性椭球体，求解出椭球体上最大值RQD_tmax和最小值RQD_tmin的值；

8.2.2：求解出RQD_t的空间各向异性比值，空间各向异性比计算公式如下：

式中：RQD_tmax为椭球体上的RQD_t最大值，RQD_tmin为椭球体上的RQD_t最小值；

8.2.3：各向异性比越大，岩体的不利方位越明显；

8.2.4：判断巷道开挖方向与RQD_tmin方向的关系，当巷道轴线与RQD_tmin方向一致时，对围岩稳定非常不利；当巷道轴线与RQD_tmin方向垂直时，对围岩稳定有利；

8.2.5：得到巷道开挖的空间不利方位和最优方位。

3.如权利要求1所述的一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，其特征在于，所述步骤(7)中，空间RQD_t求解方法过程如下：

7.1：平面RQD_t求解方法

根据最佳阈值t下的36个方位的RQD_t值，求解出任意角度的平面RQD_t值，过程如下：

7.1.1：求解出最佳阈值t下的36个RQD_t值；

7.1.2：以RQD_t值沿x方向的分量RQD_x作为横坐标，沿y轴方向的分量RQD_y作为纵坐标，设定拟合圆锥曲线方程为：

式中：p₁，p₂，p₃为拟合系数；

7.1.3：基于最小二乘法，应用Matlab软件编写拟合程序；

7.1.4：拟合得到RQD_t的椭圆曲线；

7.1.5：通过该拟合椭圆曲线，可以得到平面上任意角度的RQD_t值；

7.2：空间RQD_t求解方法

根据拟合得到的RQD_t椭圆曲线，求解出空间内任意方位任意角度的RQD_t值，过程如下：

7.2.1：根据拟合得到的RQD_t椭圆曲线，求解出椭圆的长轴所在位置；

7.2.2：以椭圆曲线的长轴作为旋转轴，利用Matlab软件编写程序，将椭圆曲线沿长轴旋转360°，得到RQD_t的各向异性椭球体；

7.2.3：得到空间内任意方位任意角度的RQD_t值。

4.如权利要求1所述的一种基于激光扫描、BQ和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.如权利要求1所述的一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，其特征在于，所述步骤(5)中，RQD_t各向异性图绘制过程如下：

5.1：RQD_t求解计算

基于RQD_t理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的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

y₀＝b/2

x_m+1＝X_Z，α

y_m+1＝Y_Z，α

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各向异性图绘制

根据求解出的36个RQD_t值，绘制出RQD_t的各向异性图，过程如下：

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

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

的点标出来；

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

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

5.2.5：输出RQD_t各向异性图；

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

用于绘制不同阈值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的值。

6.如权利要求1所述的一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，其特征在于，所述步骤(4)中，岩体三维裂隙网络模型生成和剖切过程如下：

4.1：随机数求解

基于蒙特卡洛方法和节理分布形态，求解随机数，过程如下；

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

4.1.2：给定分布下随机变量数值的产生，Monte Carlo方法是根据确立的结构面几何概率模型，再现服从这种模型的结构面网络模型，该方法的核心是抽样数的随机性，高品质随机数可以得到良好的模拟结果，即在(0，1)区间上生成均匀分布随机变量ri，利用这些均匀随机变量产生服从其他分布的随机数；

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：裂隙网络模型数据输出

实现三维裂隙网络模型的数据输出和任意二维裂隙网络模型的数据输出，过程如下；

4.4.1：输出st.dat文件数据信息；

4.4.2：输出st1.dat文件数据信息；

4.4.3：输出st2.dat文件数据信息。

7.如权利要求1所述的一种基于激光扫描、BQ和RQD_t的巷道开挖不利方位求解方法，其特征在于，所述步骤(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的特征值(τ₁，τ₂，t₃)和特征向量(ξ₁，ξ₂，ξ₃)，其中τ₁＜τ₂＜τ₃，最大特征值对应的特征向量ξ₃为组内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：数据输出

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

2.6：倾向玫瑰花图绘制，过程如下：

2.6.1：将节理倾向数据，按节理倾向方位角大小依次排序，每隔θ＝10°角度进行分组，每组命名D_j；

2.6.2：统计每组节理的数目

和每组节理的平均倾向

2.6.3：根据作图的大小和各组节理数目，选取一定长度的线段代表一组节理，确定线段的比例尺L_D；

2.6.4：以等于按比例尺L_D表示的、数目最多的一组节理

的线段长度

为半径作圆，过圆心作南北线及东西线，在圆周上标明方位角；

2.6.5：对每组节理D_j，按平均倾向

为方位角，在半圆上做出记号，自圆心向圆周记号点半径方向，按组内数目

和比例尺L_D定出一点

该点即表示该组节理平均倾向和节理数目；

2.6.6：顺次连接

和

如其中某组节理数为零，则连线回到圆心，再由圆心引出与下一组相连；

2.6.7：绘制出结构面倾向玫瑰花图。

Images