t求解方法附图说明" src="/d/file/p/2020/12-30/39771991fb8c7e48c5086966a3b03cc8.gif" />
一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法
技术领域
本发明涉及一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法，特别的是本发明基于摄影测量、模糊等价分析、BQ指标、裂隙网络模型、广义RQD理论和圆锥曲线方程，反演计算出RQDt最佳阈值t的范围和最佳阈值t值，得到了最佳阈值t下的RQDt椭圆曲线和各向异性椭球体，提供了一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法，属于RQDt各向异性求解领域。
背景技术
在自然界中，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的求解计算方法，因此众多针对RQDt空间效应的研究，主要还是集中在三维RQDt的计算上，如开展不同阈值t下的三维RQDt计算，或者在此基础上研究不同方向不同位置上的RQDt值，或者建立三维空间云图。很少有学者针对最佳阈值t，开展最能反映岩体质量的RQDt空间效应研究。因此，得到最佳阈值t下的空间RQDt的求解方法具有重要的理论意义和研究意义。
鉴于此，本发明提出了一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法。
发明内容
为了实现岩体RQDt空间各向异性的求解计算，本发明提供了一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法。本发明基于摄影测量、模糊等价分析、BQ指标、裂隙网络模型、广义RQD理论和圆锥曲线方程，反演计算出RQDt最佳阈值t的范围和最佳阈值t值，得到了最佳阈值t下的RQDt椭圆曲线和各向异性椭球体，提供了一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法。
为了解决上述技术问题，本发明提供如下的技术方案：
一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法，所述方法包括以下步骤：
(1)结构面数字摄影测量快速获取，过程如下：
1.1：根据观测区域岩体范围与空间位置，选取表面节理发育较好且无障碍物的岩体作为摄影测量区域，将标杆垂直的立于测量区域一侧，用于标定最后生成的三维图像上任意两点间的距离；
1.2：在岩体表面选取一块露头的、区域比较大、比较平滑的结构面作为标定点，用罗盘测量出倾向及倾角，并做标记，用于后处理时图像的方位真实化；
1.3：使用高分辨率相机，在选定区域正前方左、右位置，对岩体依次拍照，两次拍摄时，镜头离所测岩体的距离D及两次成像位置之间的距离B满足关系B＝D/8～D/5；
1.4：测点数据采集完成后，取回标杆，返回室内作进一步后处理操作；
1.5：将野外摄影测量所获取的左、右视图导入到软件分析系统，采用基准标定、像素点匹配、图像变形偏正纠正对左、右视图中的像素点进行匹配，合成岩体表面三维实体模型；
1.6：根据标杆尺寸以及罗盘量出的标定点产状，进行三维实体模型方位、尺寸和距离的真实化；
1.7：基于真实化的实体模型，识别定位出每条结构面，并导出结构面数据信息；
(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求解方法。
进一步，所述步骤(7)中，空间RQDt求解方法的过程如下：
7.1：平面RQDt求解方法
根据最佳阈值t下的36个方位的RQDt值，求解出任意角度的平面RQDt值，过程如下：
7.1.1：求解出最佳阈值t下的36个RQDt值；
7.1.2：以RQDt值沿x方向的分量RQDx作为横坐标，沿y轴方向的分量RQDy作为纵坐标，设定拟合圆锥曲线方程为：
式中：p1，p2，p3为拟合系数；
7.1.3：基于最小二乘法，应用Matlab软件编写拟合程序；
7.1.4：拟合得到RQDt的椭圆曲线；
7.1.5：通过该拟合椭圆曲线，可以得到平面上任意角度的RQDt值；
7.2：空间RQDt求解方法
根据拟合得到的RQDt椭圆曲线，求解出空间内任意方位任意角度的RQDt值，过程如下：
7.2.1：根据拟合得到的RQDt椭圆曲线，求解出椭圆的长轴所在位置；
7.2.2：以椭圆曲线的长轴作为旋转轴，利用Matlab软件编写程序，将椭圆曲线沿长轴旋转360°，得到RQDt的各向异性椭球体；
7.2.3：得到空间内任意方位任意角度的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求解计算
基于RQDt理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的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
y0＝b/2
xm+1＝XZ，α
ym+1＝YZ，α
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各向异性图绘制
根据求解出的36个RQDt值，绘制出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各向异性图绘制
用于绘制不同阈值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：给定分布下随机变量数值的产生，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.2：基于模糊等价聚类算法的结构面聚类分析，过程如下：
2.2.1：设结构面的实测样本数为N，第i个样本表示为(xi1，xi2)，xi1为结构面倾向，xi2为结构面倾角；
2.2.2：计算模糊关系矩阵R和相似系数rij；
2.2.3：求解闭包t(R)；
2.2.4：进行结构面分组判断；
2.3：结构面极点图绘制，过程如下：
2.3.1：采用下半球等角度投影方法绘制极点图；
2.3.2：将以倾向和倾角表示的节理产状数据转换为以节理单位法向量表示的结构面产状数据，得到以单位法向量表示的结构面产状数据；
2.3.3：求解出所有结构面法线的赤平投影坐标点；
2.3.4：绘制一条直径为单位长度的基圆，绘制出铅直和水平两条直径，并标出E、S、W、N；
2.3.5：将所有结构面的赤平投影坐标，绘制在基圆图上，实现结构面极点图的绘制；
2.4：结构面统计分析，过程如下：
2.4.1：确定样本分区区间m；
2.4.2：求解样本极差；
2.4.3：计算每个分区区间Mm；
2.4.4：确定样本落在每个分区区间里的概率，先利用计算机循环语言统计落在每一个区间的样本个数，结合样本总数N，计算样本数概率；
2.4.5：求解样本均值和样本方差；
2.4.6：根据概率值，绘制出每组结构面的倾向、倾角、迹长、间距和断距的概率分布形态；
2.5：输出结构面数据和概率分布形态；
2.6：走向玫瑰花图绘制，过程如下：
2.6.1：将测量所得的节理走向数据，换算成北东和北西方向，按节理走向方位角大小依次排序，每隔α＝10°角度进行分组，每组命名Ti；
2.6.2：统计每组节理的数目
2.6.3：根据作图的大小和各组节理数目，选取一定长度的线段代表一组节理，确定线段的比例尺LT；
2.6.4：以等于按比例尺LT表示的、数目最多的一组节理
2.6.5：对每组节理Ti，按平均走向
2.6.6：顺次连接
2.6.7：绘制出节理走向玫瑰花图。
本发明具有以下有益效果：
1、本发明基于摄影测量、模糊等价分析、BQ指标、裂隙网络模型、广义RQD理论和圆锥曲线方程，提供了一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法；
2、本发明实现了结构面快速摄影测量和模糊等价聚类分析；
3、本发明实现了岩体BQ指标分级；
4、本发明实现了三维裂隙网络模型生成和二维剖面模型生成；
5、本发明实现了RQDt求解计算和RQDt各向异性图的绘制；
6、本发明实现了RQDt的最佳阈值t求解；
7、本发明实现了RQDt椭圆曲线和各向异性椭球体的求解绘制；
8、本发明方法明确，所得最佳阈值t能够反映出真实岩体质量，所得RQDt椭圆曲线和各向异性椭球体能够反映出真实岩体的RQDt各向异性特征，适用于岩体RQDt的空间各向异性求解。
附图说明
图1是结构面极点图。
图2是二维裂隙网络图。
图3是不同阈值t下的RQDt各向异性图。
图4是RQDt随阈值t的拟合曲线。
图5是RQDt椭圆曲线。
具体实施方式
下面参照附图对本发明做进一步说明。
参照图1～图5，一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法，包括以下步骤：
1)结构面数字摄影测量快速获取，过程如下：
1.1：根据观测区域岩体范围与空间位置，选取表面节理发育较好且无障碍物的岩体作为摄影测量区域，将标杆垂直的立于测量区域一侧，用于标定最后生成的三维图像上任意两点间的距离；
1.2：在岩体表面选取一块露头的、区域比较大、比较平滑的结构面作为标定点，用罗盘测量出倾向及倾角，并做标记，用于后处理时图像的方位真实化；
1.3：使用高分辨率相机，在选定区域正前方左、右位置，对岩体依次拍照，两次拍摄时，镜头离所测岩体的距离D及两次成像位置之间的距离B满足关系B＝D/8～D/5；
1.4：测点数据采集完成后，取回标杆，返回室内作进一步后处理操作；
1.5：将野外摄影测量所获取的左、右视图导入到软件分析系统，采用基准标定、像素点匹配、图像变形偏正纠正对左、右视图中的像素点进行匹配，合成岩体表面三维实体模型；
1.6：根据标杆尺寸以及罗盘量出的标定点产状，进行三维实体模型方位、尺寸和距离的真实化；
1.7：基于真实化的实体模型，识别定位出每条结构面，并导出结构面数据信息；
2)结构面极点图和走向玫瑰花图绘制，过程如下：
2.1：识别数字摄影测量获得的结构面数据；
2.2：基于模糊等价聚类算法的结构面聚类分析
基于模糊等价聚类算法，实现结构面产状的模糊等价聚类，过程如下：
2.2.1：设结构面的实测样本数为N，第i个样本表示为(xi1，xi2)，xi1为结构面倾向，xi2为结构面倾角，其模糊关系矩阵R为：
矩阵中的元素rij为第i个样本与第j个样本间的相似系数，表征其两者间的相似程度；rij越大表示样本i与样本j越相似；
2.2.2：计算相似系数rij：
式中，i＝1，2...N；j＝1，2...N；c为计算参数(0≤c≤1)，适当选取c值使得rij在[0，1]中分散开来；
2.2.3：求解闭包t(R)：
R2＝RR
R4＝R2R2
…(3)
2.2.4：进行结构面分组判断：模糊矩阵乘法运算步骤与普通矩阵乘法相似，不同的是，并非先两项相乘后再相加，而是先取小后取大；若C＝AB，则C中的元素
当Rn＝Rn+1＝Rn+2＝…时(5)
有，模糊等价矩阵t(R)＝Rn (6)
取定截集水平λ∈[0，1]，若t(R)中rij≥λ，则结构面i和j属于同一类；即
rij≥λ (7)
2.3：极点图绘制
根据聚类结果和结构面分组结果，绘制出结构面产状极点图，过程如下；
2.3.1：采用下半球等角度投影方法绘制极点图；
2.3.2：将以倾向αd和倾角βd表示的节理产状数据转换为以节理单位法向量表示的结构面产状数据，设αn和βn分别为结构面单位法向量的倾伏向和倾伏角，对于任意结构面的单位法向量表示为X＝(x1，x2，x3)，此时半球面上每个点都对应一个节理产状，公式为：
X＝(x1，x2，x3)(8)
αd∈(0，360)，βd∈(0，90) (11)
2.3.3：得到以单位法向量表示的结构面产状数据；
2.3.4：基于结构面法向产状数据，根据结构面空间赤平投影图的纵剖面原理图，A’点为该平面法线的赤面投影；结合赤平投影原理图，计算出A’在赤平投影图上的坐标xn和yn，公式如下：
2.3.5：求解出所有结构面法线的赤平投影坐标点(xn，yn)；
2.3.6：绘制一条直径为单位长度的基圆，绘制出铅直和水平两条直径，并标出E、S、W、N；
2.3.7：将所有结构面的赤平投影坐标(xn，yn)，绘制在基圆图上；
2.3.8：实现结构面极点图的绘制，如图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：数据输出
输出结构面产状的分类信息，包括每组结构面的倾向、倾角、迹长、间距和断距的均值和方差；
2.6：走向玫瑰花图绘制
根据走向玫瑰花图绘制方法，绘制出结构面走向玫瑰花图，过程如下；
2.6.1：将测量所得的节理走向数据，换算成北东和北西方向，按节理走向方位角大小依次排序，每隔α＝10°角度进行分组，每组命名Ti；
Ti＝{α，α+9°}(19)
α＝10(i-1) (20)
i∈(1，10)∪(27，36) (21)
2.6.2：统计每组节理的数目
i∈(1，10)∪(27，36) (23)
2.6.3：根据作图的大小和各组节理数目，选取一定长度的线段代表一组节理，确定线段的比例尺LT；
2.6.4：以等于按比例尺LT表示的、数目最多的一组节理
2.6.5：对每组节理Ti，按平均走向
2.6.6：顺次连接
2.6.7：绘制出节理走向玫瑰花图；
3)基于BQ指标的岩体质量计算，过程如下：
3.1：根据结构面参数计算岩体完整性系数，公式如下：
式中：Jv为岩体体积节理数，单位条/m3；
3.2：Jv计算公式如下：
式中：L1，L2，...，Ln为垂直于结构面测线长度；N1，N2，...，Nn为同组结构面数目；
3.3：根据岩石单轴抗压强度值和岩体完整性系数值，计算BQ值：
BQ＝90+3Rc+250Kv(28)
式中：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) (29)
式中：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.dar文件中，此时剖面处于三维坐标系下，文件中数据格式自左向右依次为：节理中心点坐标(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.dar文件数据信息；
4.4.3：输出st2.dat文件数据信息；
5)RQDt各向异性图绘制，过程如下：
5.1：RQDt求解计算
基于RQDt理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的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 (39)
y0＝b/2(40)
xm+1＝XZ，α(41)
ym+1＝YZ，α(42)
5.1.11：输入一个阈值t；
5.1.12：循环比较di和t的大小，设初始lt＝0，规则如下：
lt＝lt+di，若di＞t(43)
5.1.13：求解每条测线对应的RQDt值，用mα表示，以及对应的测线起点与终点的距离lα，公式如下：
5.1.14：循环求出每条测线所对应的mα，获得36个方向上的RQDt值；
5.2：RQDt各向异性图绘制
根据求解出的36个RQDt值，绘制出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各向异性图绘制
用于绘制不同阈值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：平面RQDt求解方法
根据最佳阈值t下的36个方位的RQDt值，求解出任意角度的平面RQDt值，过程如下：
7.1.1：求解出最佳阈值t下的36个RQDt值；
7.1.2：以RQDt值沿x方向的分量RQDx作为横坐标，沿y轴方向的分量RQDy作为纵坐标，设定拟合圆锥曲线方程为：
式中：p1，p2，p3为拟合系数；
7.1.3：基于最小二乘法，应用Matlab软件编写拟合程序；
7.1.4：拟合得到RQDt的椭圆曲线，如图5所示；
7.1.5：通过该拟合椭圆曲线，可以得到平面上任意角度的RQDt值；
7.2：空间RQDt求解方法
根据拟合得到的RQDt椭圆曲线，求解出空间内任意方位任意角度的RQDt值，过程如下：
7.2.1：根据拟合得到的RQDt椭圆曲线，求解出椭圆的长轴所在位置；
7.2.2：以椭圆曲线的长轴作为旋转轴，利用Matlab软件编写程序，将椭圆曲线沿长轴旋转360°，得到RQDt的各向异性椭球体；
7.2.3：得到空间内任意方位任意角度的RQDt值。
