基于激光扫描、BQ、RQDt各向异性的Mathews稳定图法的改进方法与流程
本发明涉及一种基于激光扫描、bq、rqdt各向异性的mathews稳定图法的改进方法，特别的是本发明基于三维激光扫描、近邻传播算法、bq指标、裂隙网络模型和广义rqd理论，反演计算出rqdt最佳阈值t的范围和最佳阈值t值，给出了rqdt各向异性的求解方法，并基于rqdt的各向异性，提供了一种基于激光扫描、bq、rqdt各向异性的mathews稳定图法的改进方法，属于mathews稳定图法的方法改进领域。
背景技术：
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，用rqdt表示。广义rqd的引入，让rqd各向异性的求解成为了可能。
由rqdt概念可以知道，阈值t是广义rqd的一项特别重要的参数，是广义rqd能否真实反映出岩体质量的关键，但是阈值t具有任意性，不具有唯一性，因此寻找到能够表征出岩体质量的最佳阈值t，具有非常重要的科学意义和研究意义。
rqdt的阈值t受测线方向、结构面形态和分布特征的影响，而岩体中广泛发育的结构面，是破坏岩体连续性、完整性，控制岩体力学特征和稳定性的重要因素，对岩体质量起着控制性作用。结构面的发育模式和分布形态非常复杂，但同时，不同的节理组之间、节理和断层之间，又存在着一定的成生关系，构成某种特定的组合，表现出一定的规律性。因此，对结构面进行准确和有效的描述分析，获得结构面的产状和分布特征，是研究岩体质量的基础，也是最佳阈值t的研究基础。三维激光扫描技术作为一种高效的三维空间信息获取手段，在获取结构面方位和规模信息方面具有很大优势。近邻传播聚类算法在结构面产状聚类分析中具有较好效果，它具有不受初始聚类中心影响、计算效率高的优点，已经广泛应用于很多领域。因此，基于三维激光扫描技术和近邻传播算法，可以实现结构面的快速识别获取和聚类分析，为rqdt的最佳阈值t的研究，提供数据基础。
国内外学者针对阈值t的研究，主要体现在以下两个方面：不同阈值t下的rqdt计算和最佳阈值t的研究。在不同阈值t下的rqdt计算方面，已有的研究主要是探讨阈值t的变化对rqdt值的影响，主要是通过模拟岩体裂隙网络模型，布置虚拟钻孔，以求解不同rqdt值，研究的载体是裂隙网络模型。但在这个方面的研究中，主要是为了研究rqdt的各向异性特征，对最佳阈值t的研究鲜有涉及。
在最佳阈值t的研究方面，有些学者也曾对最佳阈值t展开过探讨。如有的学者基于三维结构面网络模拟技术，运用分形理论分别计算其结构面分布的分形维数，通过不断改变广义rqd的阈值，得到不同阈值下的rqdt值，将各剖面的分形维数与广义rqd值对比分析，为准确选取广义rqd的最佳阈值提供依据。有的学者基于修改后的块度指数mbi，建立了35种假想的三维裂隙网络模型，测量了具有不同阈值的广义rqd值，确定广义rqd的最佳阈值。这两个方面的最佳阈值t的研究，适用条件都有一定的局限性，都是在特定理论和背景下的一种最佳阈值t求解方法，当背景或者模型变化了之后，其最佳阈值t将不再具有最佳性。而且，由于分形维数或者块度指数本身并不具备表征岩体质量的功能，通过分形维数或者块度指数获得的最佳阈值t，能否反映真实的岩体质量，还有值得商榷的地方。
因此，找到并提出一种最能反映出岩体质量的rqdt最佳阈值t的求解方法，具有非常重要的科学和研究意义，也是rqdt阈值t研究中亟需解决的问题。
rqd的各向异性直接影响岩体质量，rqd各向异性对岩体质量的影响机理，目前还尚未探索清楚。由于在阈值t方面，还尚未有学者给出最佳阈值t的求解计算方法，因此也没有得到基于最佳阈值t的rqdt各向异性计算公式，也尚未获得最能反映岩体质量的rqdt各向异性求解方法。而获得rqdt各向异性求解方法是研究岩体质量各向异性的基础和前提。
鉴于此，本发明提出了一种基于激光扫描、bq、rqdt各向异性的mathews稳定图法的改进方法。
技术实现要素：
为了实现mathews稳定图法的方法改进，本发明提供了一种基于激光扫描、bq、rqdt各向异性的mathews稳定图法的改进方法。本发明基于三维激光扫描、近邻传播算法、bq指标、裂隙网络模型和广义rqd理论，反演计算出rqdt最佳阈值t的范围和最佳阈值t值，给出了rqdt各向异性的求解方法，并基于rqdt的各向异性，提供了一种基于激光扫描、bq、rqdt各向异性的mathews稳定图法的改进方法。
为了解决上述技术问题，本发明提供如下的技术方案：
一种基于激光扫描、bq、rqdt各向异性的mathews稳定图法的改进方法，所述方法包括以下步骤：
(1)结构面三维激光扫描快速获取，过程如下：
1.1：根据扫描目标和场地条件，选择扫描机位点，架设三脚架，架设中要确保仪器按照一定的扫描路线可以完整的获取边坡岩体的三维空间点云信息，同时要尽可能保证三脚架台面水平，并放置控制靶；
1.2：放置扫描仪主机在三脚架台面，固定旋钮，通过粗调脚架及微调扫描仪底座使主机气泡居中，设置扫描仪端口参数；
1.3：启动扫描控制软件，配置扫描仪相关参数，进入扫描仪控制界面，规划扫描角度，根据扫描目标设置扫描范围，调整相机配置参数，获取扫描目标图像；
1.4：固定扫描范围，获取扫描间距，设定采样间距，开始数据获取，并实时查看扫描点云数据及彩色信息情况，根据扫描成果随时调整扫描参数设定；
1.5：导出结构面点云数据；
(2)结构面聚类分析；
(3)基于bq指标的岩体质量计算，过程如下：
3.1：根据结构面参数计算岩体完整性系数，公式如下：
式中：jv为岩体体积节理数，单位条/m3；
3.2：jv计算公式如下：
式中：l1，l2，...，ln为垂直于结构面测线长度；n1，n2，...，nn为同组结构面数目；
3.3：根据岩石单轴抗压强度值和岩体完整性系数值，计算bq值：
bq＝90+3rc+250kv
式中：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(k1+k2+k3)
式中：k1为地下水影响修正系数；k2为软件结构面产状影响修正系数；k3为天然应力影响修正系数；
3.6：得到修正后的bq岩体分级结果；
(4)岩体三维裂隙网络模型生成和剖切；
(5)rqdt各向异性图绘制；
(6)基于bq反演的最佳阈值t求解方法；
(7)rqdt各向异性求解方法；
(8)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′值时，用rqdt值替代钻孔岩芯记录的rqd值，并将rqdt的各向异性特征纳入考虑范畴，根据rqdt各向异性求解公式，得到考虑rqdt各向异性特征的q′值，公式如下：
式中：rqdtθ为角度θ下的rqdt值；rqdtmin为阈值t下的rqd最小值；
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)中，rqdt各向异性求解方法的过程如下：
7.1：基于最佳阈值t和剖切的三个二维裂隙网络模型，求解出最佳阈值t下的rqdt值，每个裂隙网络模型求解出36个rqdt值；
7.2：计算出每个裂隙网络模型上rqdt的最大值rqdtmax、最小值rqdtmin和均值
7.3：根据研究方向和rqdtmin方向位置关系，结合方差进行修正，提出各向异性条件下rqdt的计算公式如下：
式中：
7.4：修正系数rqd′t的求解考虑研究方向在某一具体方位角下的rqdtθ值和方差d，按如下方法修正，当rqdtθ＝rqdtmin时，rqd′t＝-d，d为方差；当rqdtθ＝rqdtmax时，rqd′t＝d；当
7.5：根据修正系数rqd′t的求解过程，提出rqd′t修正公式如下：
式中：rqdtθ为θ角下的rqdt值，rqdtmin为阈值t下的rqdt最小值，a，b为相关系数；
7.6：提出rqdt各向异性计算公式如下：
7.7：相关系数a和b的求解，过程如下：
7.7.1：计算出每个裂隙网络模型上rqdt的最大值rqdtmax、最小值rqdtmin和均值
7.7.2：根据rqdtθ与rqdtmax、rqdtmin和
7.7.3：根据散点图做出拟合曲线，曲线截距即为a值，斜率即为b值，曲线方程即为rqd′t的修正公式；
7.8：将求解出的a、b值，带入到修正系数公式及rqdt各向异性计算公式中，得到与角度θ有关的rqdt各向异性公式；
7.9：根据rqdt各向异性公式，求解出任意角度的rqdt值。
进一步，所述步骤(6)中，基于bq反演的最佳阈值t求解方法的过程如下：
6.1：基于bq指标反演rqdt范围，过程如下：
结合bq分级计算出的岩体质量级别，查找《岩石质量指标》表，反演确定出该岩体级别下，rqdt范围值；
6.2：最佳阈值t求解方法，过程如下：
6.2.1：在三维裂隙网络模型上，过中心点o，以任意角度剖切三个剖面，得到三个二维裂隙网络模型，导出二维裂隙网络模型及数据；
6.2.2：针对每一个二维裂隙网络模型，设置不同的阈值t，求解出不同阈值t下的rqdt值；
6.2.3：导出不同阈值t下的rqdt值，计算出rqdt均值；
6.2.4：以阈值t为横坐标，以rqdt的均值为纵坐标，绘制rqdt随阈值t变化的散点图；
6.2.5：根据散点图，设置拟合方程，拟合rqdt随阈值t变化的曲线图；
6.2.6：将反演出的rqdt范围值，带入到拟合出的rqdt随阈值t变化的曲线图中，结合函数方程和曲线图，求解出在该rqdt范围内阈值t的范围，共得到三组阈值t的范围；
6.2.7：针对三组阈值t的范围，取其范围的交集，作为最佳阈值t的范围；
6.2.8：以最佳阈值t范围的中点值作为最佳阈值t值，得到最佳阈值t；
6.2.9：输出最佳阈值t的范围和最佳阈值t值。
进一步，所述步骤(5)中，rqdt各向异性图绘制的过程如下：
5.1：rqdt求解计算，过程如下：
5.1.1：rqdt理论公式如下：
式中：xi表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，rqdt代表对应阈值t的岩石质量指标，即阈值t下的rqdt值；
5.1.2确定二维裂隙网络模型的剖面中心点o、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，求解出中心点o坐标(x0，y0)和边界方程；
5.1.3：以o为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度l等于o点到裂隙网络模型边界的距离，用l0～l35表示，求解出测线方程；
5.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(xa，ya)，将测线方程与边界方程依次连立，判断测线是否与边界相交；
5.1.5：求出测线与裂隙网络模型的边界方程交点(xz，α，yz，α)；
5.1.6：确定测线所在区间，原理如下：
5.1.7：根据裂隙网络模型中每条节理的起点坐标(xb，yb)和终点坐标(xc，yc)，建立相应的解析方程；
5.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(xj，yj)范围，如果交点符合a＜xj＜c且b＜yj＜d，则记录该交点，遍历所有节理方程求出所有m个交点；
5.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；
5.1.10：计算相邻交点的距离di；
5.1.11：输入一个阈值t；
5.1.12：循环比较di和t的大小，设初始lt＝0，规则如下：
lt＝lt+di，若di＞t
5.1.13：求解每条测线对应的rqdt值，用mα表示，以及对应的测线起点与终点的距离lα；
5.1.14：循环求出每条测线所对应的mα，获得36个方向上的rqdt值；
5.2：rqdt各向异性图绘制，过程如下：
5.2.1：将36个rqdt值，按角度顺序依次排序；
5.2.2：以o点为圆心，以1为半径画圆，在射线角度为α上找到与圆心距离为
5.2.3：依次连接36组射线的端点，若某一条射线上的rqdt为0，则取圆心；
5.2.4：绘制出rqdt的各向异性图；
5.2.5：输出rqdt各向异性图；
5.3：不同阈值t下rqdt各向异性图绘制，过程如下：
5.3.1：输入不同的阈值t，求解出对应的rqdt值；
5.3.2：将不同阈值t下的rqdt值，以角度为变量，以rqdt值为函数，绘制到同个坐标系下；
5.3.3：得到不同阈值t下的rqdt各向异性图；
5.3.4：输出不同阈值t下的rqdt各向异性图；
5.4：输出rqdt的值。
进一步，所述步骤(4)中，岩体三维裂隙网络模型生成和剖切的过程如下：
4.1：随机数求解，过程如下：
4.1.1：产生随机数的数学方法应满足以下条件：产生的随机数列应均匀分布在(0，1)区间；序列之间应无相关性；随机序列有足够长的重复周期，在计算机上产生的速度快，占有的内存空间小，具有完全可重复性；
4.1.2：利用montecarlo方法是根据确立的结构面几何概率模型，再现出服从这种模型的结构面网络模型，在(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，每个样本数据的倾向为xi，倾角为yi，以每个样本数据的倾向xi，倾角yi作为一个聚类，确定一个初始聚类中心，共得到n个初始聚类中心；
2.2.2：通过相似性度量准则，遍历所有样本数据，计算每个样本数据距离聚类中心的距离，并将每个样本数据分配到距离它最近的聚类中心，得到n组数据；
2.2.3：对于每组数据，通过特征模量分析方法，求解计算每组数据的聚类中心，假设某组内存在l个数据，首先，按如下公式计算矩阵s：
式中：(xi，yi，zi)为任意结构面的单位法向量，i∈(1，l)；
2.2.4：其次，求解矩阵s的特征值(τ1，τ2，τ3)和特征向量(ξ1，ξ2，ξ3)，其中τ1＜τ2＜τ3，最大特征值对应的特征向量ξ3为组内l个向量的平均向量，将ξ3作为新的聚类中心；
2.2.5：针对所有样本数据，重复计算每个样本数据距离聚类中心的距离、矩阵s以及特征值和特征向量，直到所有聚类中心的位置都固定，确定出结构面的分组；
2.2.6：将以单位法向量表示的结构面产状数据转换为以倾向、倾角表示的结构面产状数据；
2.2.7：对结构面产状数据进行统计分析，计算结构面倾角的平均值m与标准差σ，计算倾角数据的稳健区间[m-σ，m+σ]；
2.2.8：判断样本数据的初始聚类中心的倾向xi和倾角yi是否落在稳健区间稳健区间[m-σ，m+σ]，若是，则聚类分析完成；若不是，则需要对样本数据重新聚类，直到初始聚类中心的倾向xi和倾角yi均落在稳健区间内[m-σ，m+σ]；
2.3：极点图绘制，过程如下；
2.3.1：基于结构面法向产状数据，根据结构面空间赤平投影图的纵剖面原理，求解出所有结构面法线的赤平投影坐标点(xn，yn)；
2.3.3：绘制一条直径为单位长度的基圆，绘制出铅直和水平两条直径，并标出e、s、w、n；
2.3.4：将所有结构面的赤平投影坐标(xn，yn)，绘制在基圆图上；
2.3.5：实现结构面极点图的绘制；
2.4：结构面统计分析，过程如下：
2.4.1：确定样本分区区间m；
2.4.2：求解样本极差
2.4.3：计算每个分区区间mm：
2.4.4：确定样本落在每个分区区间里的概率，先利用计算机循环语言统计落在每一个区间的样本个数nm，结合样本总数n，计算样本数概率pm；
2.4.5：求解样本均值
2.4.6：求解样本方差s2，其中s为标准差；
2.4.7：根据概率pm值，绘制出每组结构面的倾向、倾角、迹长、间距和断距的概率分布形态；
2.5：输出结构面产状的分组信息，包括每组结构面的倾向、倾角、迹长、间距和断距的均值和方差。
本发明具有以下有益效果：
1、本发明基于三维激光扫描、近邻传播算法、bq指标、裂隙网络模型和广义rqd理论，提供了一种基于激光扫描、bq、rqdt各向异性的mathews稳定图法的改进方法；
2、本发明实现了结构面快速三维激光扫描和近邻传播算法聚类分析；
3、本发明实现了岩体bq指标分级；
4、本发明实现了三维裂隙网络模型生成和二维剖面模型生成；
5、本发明实现了rqdt求解计算和rqdt各向异性图的绘制；
6、本发明实现了rqdt的最佳阈值t求解；
7、本发明实现了rqdt各向异性的求解；
8、本发明实现了mathews稳定图法的改进；
9、本发明方法明确，提供的mathews稳定图法的改进方法考虑了rqdt的各向异性特征。
附图说明
图1是结构面极点图。
图2是二维裂隙网络图。
图3是不同阈值t下的rqdt各向异性图。
图4是rqdt随阈值t的拟合曲线。
图5是rqd′t的拟合曲线。
具体实施方式
下面参照附图对本发明做进一步说明。
参照图1～图5，一种基于激光扫描、bq、rqdt各向异性的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＝(x1，x2，x3)，此时半球面上每个点都对应一个节理产状，公式为：
x＝(x1，x2，x3)(1)
αd∈(0，360)，βd∈(0，90)(4)
2.1.6：得到以单位法向量表示的结构面数据；
2.2：近邻传播算法聚类分析，过程如下：
2.2.1：设结构面的实测样本数量为n，每个样本数据的倾向为xi，倾角为yi，i∈(1，n)；以每个样本数据的倾向xi，倾角yi作为一个聚类，确定一个初始聚类中心，共得到n个初始聚类中心；
2.2.2：通过相似性度量准则，遍历所有样本数据，计算每个样本数据距离聚类中心的距离，并将每个样本数据分配到距离它最近的聚类中心，得到n组数据；
2.2.3：对于每组数据，通过特征模量分析方法，求解计算每组数据的聚类中心，假设某组内存在l个数据，聚类中心按如下方法求解：
2.2.4：首先，按如下公式计算矩阵s
式中：(xi，yi，zi)为任意结构面的单位法向量，i∈(1，l)；
2.2.5：其次，求解矩阵s的特征值(τ1，τ2，τ3)和特征向量(ξ1，ξ2，ξ3)，其中τ1＜τ2＜τ3，最大特征值对应的特征向量ξ3为组内l个向量的平均向量，将ξ3作为新的聚类中心；
2.2.6：针对所有样本数据，重复计算每个样本数据距离聚类中心的距离、矩阵s以及特征值和特征向量，直到所有聚类中心的位置都固定，确定出结构面的分组；
2.2.7：将以单位法向量表示的结构面产状数据转换为以倾向、倾角表示的结构面产状数据；
2.2.8：对结构面产状数据进行统计分析，计算结构面倾角的平均值m与标准差σ，计算倾角数据的稳健区间[m-σ，m+σ]；
2.2.9：判断样本数据的初始聚类中心的倾向xi和倾角yi是否落在稳健区间稳健区间[m-σ，m+σ]，若是，则聚类分析完成；若不是，则需要对样本数据重新聚类，直到初始聚类中心的倾向xi和倾角yi均落在稳健区间内[m-σ，m+σ]；
2.3：极点图绘制，过程如下；
2.3.1：基于结构面法向产状数据，根据结构面空间赤平投影图的纵剖面原理，设a’点为该平面法线的赤面投影，结合赤平投影原理，计算出a’在赤平投影图上的坐标xn和yn，公式如下：
2.3.2：求解出所有结构面法线的赤平投影坐标点(xn，yn)；
2.3.3：绘制一条直径为单位长度的基圆，绘制出铅直和水平两条直径，并标出e、s、w、n；
2.3.4：将所有结构面的赤平投影坐标(xn，yn)，绘制在基圆图上；
2.3.5：实现结构面极点图的绘制，如图1所示；
2.4：结构面统计分析，过程如下：
2.4.1：确定样本分区区间m；
2.4.2：求解样本极差
2.4.3：计算每个分区区间mm：
2.4.4：确定样本落在每个分区区间里的概率，先利用计算机循环语言统计落在每一个区间的样本个数nm，结合样本总数n，计算样本数概率pm：
2.4.5：求解样本均值
2.4.6：求解样本方差s2，其中s为标准差：
2.4.7：根据概率pm值，绘制出每组结构面的倾向、倾角、迹长、间距和断距的概率分布形态；
2.5：输出结构面产状的分组信息，包括每组结构面的倾向、倾角、迹长、间距和断距的均值和方差；
3)基于bq指标的岩体质量计算，过程如下：
3.1：根据结构面参数计算岩体完整性系数，公式如下：
式中：jv为岩体体积节理数，单位条/m3；
3.2：jv计算公式如下：
式中：l1，l2，...，ln为垂直于结构面测线长度；n1，n2，...，nn为同组结构面数目；
3.3：根据岩石单轴抗压强度值和岩体完整性系数值，计算bq值：
bq＝90+3rc+250kv(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(k1+k2+k3)(16)
式中：k1为地下水影响修正系数；k2为软件结构面产状影响修正系数；k3为天然应力影响修正系数；
3.6：得到修正后的bq岩体分级结果，为iii级岩体；
4)岩体三维裂隙网络模型生成和剖切，过程如下；
4.1：随机数求解，过程如下；
4.1.1：产生随机数的数学方法应满足以下条件：产生的随机数列应均匀分布在(0，1)区间；序列之间应无相关性；随机序列有足够长的重复周期，在计算机上产生的速度快，占有的内存空间小，具有完全可重复性；
4.1.2：利用montecarlo方法是根据确立的结构面几何概率模型，再现出服从这种模型的结构面网络模型，在(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)rqdt各向异性图绘制，过程如下：
5.1：rqdt求解计算，过程如下：
5.1.1：rqdt理论公式如下：
式中：xi表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，rqdt代表对应阈值t的岩石质量指标，即阈值t下的rqdt值；
5.1.2确定二维裂隙网络模型的剖面中心点o、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，则中心点o坐标为：
边界方程为：
5.1.3：以o为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度l等于o点到裂隙网络模型边界的距离，用l0～l35表示，则测线方程为：
式中：s表示测线，α表示角度；
5.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(xa，ya)，将测线方程与边界方程依次连立，判断测线是否与边界相交，原理如下：
5.1.5：求出测线与裂隙网络模型的边界方程交点(xz，α，yz，α)；
5.1.6：确定测线所在区间，原理如下：
5.1.7：根据裂隙网络模型中每条节理的起点坐标(xb，yb)和终点坐标(xc，yc)，建立相应的解析方程，定义节理方程如下：
5.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(xj，yj)范围，如果交点符合a＜xj＜c且b＜yj＜d，则记录该交点，遍历所有节理方程求出所有m个交点；
5.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；
5.1.10：计算相邻交点的距离，公式如下：
x0＝a/2(26)
y0＝b/2(27)
xm+1＝xz，α(28)
ym+1＝yz，α(29)
5.1.11：输入一个阈值t；
5.1.12：循环比较di和t的大小，设初始lt＝0，规则如下：
lt＝lt+di，若di＞t(30)
5.1.13：求解每条测线对应的rqdt值，用mα表示，以及对应的测线起点与终点的距离lα，公式如下：
5.1.14：循环求出每条测线所对应的mα，获得36个方向上的rqdt值；
5.2：rqdt各向异性图绘制，过程如下：
5.2.1：将36个rqdt值，按角度顺序依次排序；
5.2.2：以o点为圆心，以1为半径画圆，在射线角度为α上找到与圆心距离为
5.2.3：依次连接36组射线的端点，若某一条射线上的rqdt为0，则取圆心；
5.2.4：绘制出rqdt的各向异性图；
5.2.5：输出rqdt各向异性图；
5.3：不同阈值t下rqdt各向异性图绘制，过程如下：
5.3.1：输入不同的阈值t，求解出对应的rqdt值；
5.3.2：将不同阈值t下的rqdt值，以角度为变量，以rqdt值为函数，绘制到同个坐标系下；
5.3.3：得到不同阈值t下的rqdt各向异性图；
3.3.4：输出不同阈值t下的rqdt各向异性图，如图3所示；
5.4：输出rqdt的值；
6)基于bq反演的最佳阈值t求解方法，过程如下：
6.1：基于bq指标反演rqdt范围，过程如下：
结合bq分级计算出的岩体质量级别，岩体为iii级岩体，查找《岩石质量指标》表，见表1，反演确定出该岩体级别下，rqdt范围值，rqdt范围在50％～75％
表1；
6.2：最佳阈值t求解方法，过程如下：
6.2.1：在三维裂隙网络模型上，过中心点o，以任意角度剖切三个剖面，得到三个二维裂隙网络模型，导出二维裂隙网络模型及数据；
6.2.2：针对每一个二维裂隙网络模型，设置不同的阈值t，求解出不同阈值t下的rqdt值；
6.2.3：导出不同阈值t下的rqdt值，计算出rqdt均值；
6.2.4：以阈值t为横坐标，以rqdt的均值为纵坐标，绘制rqdt随阈值t变化的散点图；
6.2.5：根据散点图，设置拟合方程y＝aexp(bx)，拟合rqdt随阈值t变化的曲线图，如图4所示；
6.2.6：将反演出的rqdt范围值，带入到拟合出的rqdt随阈值t变化的曲线图中，结合函数方程和曲线图，求解出在该rqdt范围内，阈值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)rqdt各向异性求解方法，过程如下：
7.1：基于最佳阈值t和剖切的三个二维裂隙网络模型，求解出最佳阈值t下的rqdt值，每个裂隙网络模型求解出36个rqdt值；
7.2：计算出每个裂隙网络模型上rqdt的最大值rqdtmax、最小值rqdtmin和均值
7.3：根据研究方向和rqdtmin方向位置关系，结合方差进行修正，提出各向异性条件下rqdt的计算公式如下：
式中：
7.4：修正系数rqd′t的求解考虑研究方向在某一具体方位角下的rqdtθ值和方差d，按如下方法修正，当rqdtθ＝rqdtmin时，rqd′t＝-d，d为方差；当rqdtθ＝rqdtmax时，rqd′t＝d；当
7.5：根据修正系数rqd′t的求解过程，提出rqd′t修正公式如下：
式中：rqdtθ为θ角下的rqdt值，rqdtmin为阈值t下的rqdt最小值，a，b为相关系数；
7.6：提出rqdt各向异性计算公式如下：
7.7：相关系数a和b的求解，过程如下：
7.7.1：计算出每个裂隙网络模型上rqdt的最大值rqdtmax、最小值rqdtmin和均值
7.7.2：根据rqdtθ与rqdtmax、rqdtmin和
7.7.3：根据散点图做出拟合曲线，如图5所示，曲线截距即为a值，斜率即为b值，曲线方程即为rqd′t的修正公式，得到的rqd′t公式如下：
7.8：将求解出的a、b值，带入到修正系数公式及rqdt各向异性计算公式中，得到与角度θ有关的rqdt各向异性公式：
7.9：根据rqdt各向异性公式，求解出任意角度的rqdt值；
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′值时，用rqdt值替代钻孔岩芯记录的rqd值，并将rqdt的各向异性特征纳入考虑范畴，根据rqdt各向异性求解公式，得到考虑rqdt各向异性特征的q′值，公式如下：
式中：rqdtθ为角度θ下的rqdt值；rqdtmin为阈值t下的rqd最小值；
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％的概率崩落；处在崩落区的工程，将持续发生崩落。
- 还没有人留言评论。精彩留言会获得点赞！