t最佳阈值t求解方法附图说明" src="/d/file/p/2020/12-30/aa5ee9e213d17ed719c25bfea32c4544.gif" />
一种基于BQ反演的RQDt最佳阈值t求解方法
技术领域
本发明涉及一种基于BQ反演的RQDt最佳阈值t求解方法，特别的是本发明基于BQ指标、裂隙网络模型和广义RQD理论，反演计算出RQDt最佳阈值t的范围和最佳阈值t值，提供了一种基于BQ反演的RQDt最佳阈值t求解方法，属于RQDt最佳阈值t求解领域。
背景技术
在自然界中，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，将具有非常重要的科学意义和研究意义。
国内外学者针对阈值t的研究，主要体现在以下两个方面：不同阈值t下的RQDt计算和最佳阈值t的研究。
在不同阈值t下的RQDt计算方面，已有的研究主要是探讨阈值t的变化对RQDt值的影响，主要是通过模拟岩体裂隙网络模型，布置虚拟钻孔，以求解不同RQDt值，研究的载体是裂隙网络模型。但在这个方面的研究中，主要是为了研究RQD的各向异性特征，对最佳阈值t的研究鲜有涉及。
在最佳阈值t的研究方面，有些学者也曾对最佳阈值t展开过探讨。如有的学者基于三维结构面网络模拟技术，运用分形理论分别计算其结构面分布的分形维数，通过不断改变广义RQD的阈值，得到不同阈值下的RQD值，将各剖面的分形维数与广义RQD值对比分析，为准确选取广义RQD的最佳阈值提供依据。有的学者基于修改后的块度指数MBi，建立了35种假想的三维裂隙网络模型，测量了具有不同阈值的广义RQD值，确定广义RQD的最佳阈值。这两个方面的最佳阈值t的研究，适用条件都有一定的局限性，都是在特定理论和背景下的一种最佳阈值t求解方法，当背景或者模型变化了之后，其最佳阈值t将不再具有最佳性。而且，由于分形维数或者块度指数本身并不具备表征岩体质量的功能，通过分形维数或者块度指数获得的最佳阈值t，能否反映真实的岩体质量，还有值得商榷的地方。
因此，找到并提出一种最能反映出岩体质量的RQDt最佳阈值t的求解方法，具有非常重要的科学和研究意义，也是RQDt阈值t研究中亟需解决的问题。
鉴于此，本发明提出了一种基于BQ反演的RQDt最佳阈值t求解方法。
发明内容
为了实现岩体RQDt最佳阈值t的求解计算，本发明提供了一种基于BQ反演的RQDt最佳阈值t求解方法。本发明基于BQ指标、裂隙网络模型和广义RQD理论，反演计算出RQDt最佳阈值t的范围和最佳阈值t值，提供了一种基于BQ反演的RQDt最佳阈值t求解方法。
为了解决上述技术问题，本发明提供如下的技术方案：
一种基于BQ反演的RQDt最佳阈值t求解方法，所述方法包括以下步骤：
(1)基于BQ指标的岩体质量计算，过程如下：
1.1：根据结构面参数计算岩体完整性系数，公式如下：
式中：Jv为岩体体积节理数，单位条/m3；
1.2：Jv计算公式如下：
式中：L1，L2，...，Ln为垂直于结构面测线长度；N1，N2，...，Nn为同组结构面数目；
1.3：根据岩石单轴抗压强度值和岩体完整性系数值，计算BQ值：
BQ＝90+3Rc+250Kv
式中：Rc是岩石单轴抗压强度；Kv为岩体完整性系数；
1.4：在应用BQ计算公式过程中，遵循以下条件：
当Rc＞90Kv+30时，以Rc＝90Kv+30和Kv代入计算BQ值；
当Kv＞0.04Rc+0.4时，以Kv＝0.04Rc+0.4和Rc代入计算BQ值；
1.5：根据地下水、软弱结构面产状和天然应力影响对BQ进行修正，修正公式如下：
[BQ]＝BQ-100(K1+K2+K3)
式中：K1为地下水影响修正系数；K2为软件结构面产状影响修正系数；K3为天然应力影响修正系数；
1.6：得到修正后的BQ岩体分级结果；
(2)岩体三维裂隙网络模型生成和剖切；
(3)RQDt各向异性图绘制；
(4)基于BQ反演的最佳阈值t求解方法。
进一步，所述步骤(4)中，基于BQ反演的最佳阈值t求解方法过程如下：
4.1：基于BQ指标反演RQDt范围
结合BQ分级计算出的岩体质量级别，查找《岩石质量指标》表，反演确定出该岩体级别下，RQDt范围值；
4.2：最佳阈值t求解方法
4.2.1：在三维裂隙网络模型上，过中心点O，以任意角度剖切三个剖面，得到三个二维裂隙网络模型，导出二维裂隙网络模型及数据；
4.2.2：针对每一个二维裂隙网络模型，设置不同的阈值t，求解出不同阈值t下的RQDt值；
4.2.3：导出不同阈值t下的RQDt值，计算出RQDt均值；
4.2.4：以阈值t为横坐标，以RQDt的均值为纵坐标，绘制RQDt随阈值t变化的散点图；
4.2.5：根据散点图，设置拟合方程，拟合RQDt随阈值t变化的曲线图；
4.2.6：将反演出的RQDt范围值，带入到拟合出的RQDt随阈值t变化的曲线图中，结合函数方程和曲线图，求解出在该RQDt范围内阈值t的范围，共得到三组阈值t的范围；
4.2.7：针对三组阈值t的范围，取其范围的交集，作为最佳阈值t的范围；
4.2.8：以最佳阈值t范围的中点值作为最佳阈值t值，得到最佳阈值t；
4.2.9：输出最佳阈值t的范围和最佳阈值t值。
进一步，所述步骤(3)中，RQDt各向异性图绘制过程如下：
3.1：RQDt求解计算
基于RQDt理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的RQDt值，过程如下：
3.1.1：RQDt理论公式如下：
式中：xi表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，RQDt代表对应阈值t的岩石质量指标，即阈值t下的RQDt值；
3.1.2确定二维裂隙网络模型的剖面中心点O、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，则中心点O坐标为：
边界方程为：
3.1.3：以O为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度L等于O点到裂隙网络模型边界的距离，用L0～L35表示，则测线方程为：
式中：s表示测线，α表示角度；
3.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(xa，ya)，将测线方程与边界方程依次连立，判断测线是否与边界相交，原理如下：
3.1.5：求出测线与裂隙网络模型的边界方程交点(XZ，α，YZ，α)；
3.1.6：确定测线所在区间，原理如下：
3.1.7：根据裂隙网络模型中每条节理的起点坐标(Xb，Yb)和终点坐标(Xc，Yc)，建立相应的解析方程，定义节理方程如下：
3.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(Xj，Yj)范围，如果交点符合a＜Xj＜c且b＜Yj＜d，则记录该交点，遍历所有节理方程求出所有m个交点；
3.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；
3.1.10：计算相邻交点的距离，公式如下：
x0＝a/2
y0＝b/2
xm+1＝XZ，α
ym+1＝YZ，α
3.1.11：输入一个阈值t；
3.1.12：循环比较di和t的大小，设初始lt＝0，规则如下：
lt＝lt+di，若di＞t
3.1.13：求解每条测线对应的RQDt值，用mα表示，以及对应的测线起点与终点的距离lα，公式如下：
3.1.14：循环求出每条测线所对应的mα，获得36个方向上的RQDt值；
3.2：RQDt各向异性图绘制
根据求解出的36个RQDt值，绘制出RQDt的各向异性图，过程如下：
3.2.1：将36个RQDt值，按角度顺序依次排序；
3.2.2：以O点为圆心，以1为半径画圆，在射线角度为α上找到与圆心距离为
3.2.3：依次连接36组射线的端点，若某一条射线上的RQDt为0，则取圆心；
3.2.4：绘制出RQDt的各向异性图；
3.2.5：输出RQDt各向异性图；
3.3：不同阈值t下RQDt各向异性图绘制
用于绘制不同阈值t下的岩体RQDt各向异性图，过程如下：
3.3.1：输入不同的阈值t，求解出对应的RQDt值；
3.3.2：将不同阈值t下的RQDt值，以角度为变量，以RQDt值为函数，绘制到同个坐标系下；
3.3.3：得到不同阈值t下的RQDt各向异性图；
3.3.4：输出不同阈值t下的RQDt各向异性图；
3.4：输出RQDt的值。
再进一步，所述步骤(2)中，岩体三维裂隙网络模型生成和剖切过程如下：
2.1：随机数求解
基于蒙特卡洛方法和节理分布形态，求解随机数，过程如下；
2.1.1：伪随机数的产生，产生随机数的数学方法应满足以下条件：产生的随机数列应均匀分布在(0，1)区间；序列之间应无相关性；随机序列有足够长的重复周期，在计算机上产生的速度快，占有的内存空间小，具有完全可重复性；
2.1.2：给定分布下随机变量数值的产生，Monte Carlo方法是根据确立的结构面几何概率模型，再现服从这种模型的结构面网络模型，该方法的核心是抽样数的随机性，高品质随机数可以得到良好的模拟结果，即在(0，1)区间上生成均匀分布随机变量ri，利用这些均匀随机变量产生服从其他分布的随机数；
2.1.3：节理几何参数的密度函数有正态分布、对数正态分布、负指数分布、均匀分布四种；
2.1.4：根据求得的随机数，确定用于生成节理的基本几何参数；
2.2：岩体三维裂隙网络模型生成
根据确立的结构面几何概率模型，利用求得的随机数，确定出用于生成节理的基本几何参数，再现出服从这种模型的结构面网络模型，过程如下；
2.2.1：根据结构面数据自动统计结果和求得的随机数，将每组结构面的数据保存到一个文本文件中，用st.dat表示；
2.2.2：st.dat数据内容格式依次为：每条结构面圆盘中心点坐标(x，y，z)，圆盘半径D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ和节理分组；
2.2.3：为区分不同组别的结构面，对相同组的结构面圆盘赋予相同的颜色，用数列GID表示；
2.2.4：利用Matlab软件，编写程序，读取结构面数据文件st.dat，运行后自动生成岩体三维裂隙网络模型；
2.2.5：得到岩体三维裂隙网络模型；
2.3：二维裂隙网络模型切割
在三维裂隙网络模型的基础上，实现任意角度任意方位的二维裂隙网络模型切割，过程如下；
2.3.1：在三维裂隙网络模型上，结合Matlab软件编程工具，以三维裂隙网络模型中心点为中心，实现任意角度的剖面切割功能；
2.3.2：得到穿过中心点的任意角度的二维裂隙网络模型；
2.3.3：在三维裂隙网络模型上，结合Matlab软件编程工具，在三维裂隙网络模型任意位置上，实现任意角度和方位的剖面切割功能；
2.3.4：得到任意角度和方位的二维裂隙网络模型；
2.3.5：将切割剖面上的数据，保存到st1.dat文件中，此时剖面处于三维坐标系下，文件中数据格式自左向右依次为：节理中心点坐标(x，y，z)，节理长度D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ；
2.3.6：将三维坐标系转化成二维坐标系，并将二维剖面数据保存到st2.dat文件中，文件中数据格式自左向右依次为：节理中心点坐标(x，y)，节理长度D，倾角DA，倾向DD，厚度thin，法向方向NX，NY，NZ；
2.4：实现三维裂隙网络模型的数据输出和任意二维裂隙网络模型的数据输出。
本发明具有以下有益效果：
1、本发明基于BQ指标、裂隙网络模型和广义RQD理论，提出了一种基于BQ反演的RQDt最佳阈值t求解方法；
2、本发明实现了岩体BQ指标分级；
3、本发明实现了三维裂隙网络模型生成和二维剖面模型生成；
4、本发明实现了RQDt求解计算和RQDt各向异性图的绘制；
5、本发明实现了RQDt的最佳阈值t求解；
6、本发明方法明确，所得阈值t能够反映出真实岩体质量，适用于岩体RQDt的最佳阈值t求解。
附图说明
图1是二维裂隙网络图。
图2是不同阈值t下的RQDt各向异性图。
图3是RQDt随阈值t的拟合曲线。
具体实施方式
下面参照附图对本发明做进一步说明。
参照图1～图3，一种基于BQ反演的RQDt最佳阈值t求解方法，包括以下步骤：
1)基于BQ指标的岩体质量计算，过程如下：
1.1：根据结构面参数计算岩体完整性系数，公式如下：
式中：Jv为岩体体积节理数，单位条/m3；
1.2：Jv计算公式如下：
式中：L1，L2，...，Ln为垂直于结构面测线长度；N1，N2，...，Nn为同组结构面数目；
1.3：根据岩石单轴抗压强度值和岩体完整性系数值，计算BQ值：
BQ＝90+3Rc+250Kv(3)
式中：Rc是岩石单轴抗压强度；Kv为岩体完整性系数；
1.4：在应用BQ计算公式过程中，遵循以下条件：
当Rc＞90Kv+30时，以Rc＝90Kv+30和Kv代入计算BQ值；
当Kv＞0.04Rc+0.4时，以Kv＝0.04Rc+0.4和Rc代入计算BQ值；
1.5：根据地下水、软弱结构面产状和天然应力影响对BQ进行修正，修正公式如下：
[BQ]＝BQ-100(K1+K2+K3)(4)
式中：K1为地下水影响修正系数；K2为软件结构面产状影响修正系数；K3为天然应力影响修正系数；
1.6：得到修正后的BQ岩体分级结果，岩体为III级岩体；
2)岩体三维裂隙网络模型生成和剖切，过程如下；
2.1：随机数求解
基于蒙特卡洛方法和节理分布形态，求解随机数，过程如下；
2.1.1：伪随机数的产生，产生随机数的数学方法应满足以下条件：产生的随机数列应均匀分布在(0，1)区间；序列之间应无相关性；随机序列有足够长的重复周期，在计算机上产生的速度快，占有的内存空间小，具有完全可重复性；
2.1.2：给定分布下随机变量数值的产生，Monte Carlo方法是根据确立的结构面几何概率模型，再现服从这种模型的结构面网络模型，该方法的核心是抽样数的随机性，高品质随机数可以得到良好的模拟结果，即在(0，1)区间上生成均匀分布随机变量ri，利用这些均匀随机变量产生服从其他分布的随机数；
2.1.3：节理几何参数的密度函数有正态分布、对数正态分布、负指数分布、均匀分布四种；
2.1.4：根据求得的随机数，确定用于生成节理的基本几何参数；
2.2：岩体三维裂隙网络模型生成
根据确立的结构面几何概率模型，利用求得的随机数，确定出用于生成节理的基本几何参数，再现出服从这种模型的结构面网络模型，过程如下；
2.2.1：根据结构面数据自动统计结果和求得的随机数，将每组结构面的数据保存到一个文本文件中，用st.dat表示；
2.2.2：st.dat数据内容格式依次为：每条结构面圆盘中心点坐标(x，y，z)，圆盘半径D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ和节理分组；
2.2.3：为区分不同组别的结构面，对相同组的结构面圆盘赋予相同的颜色，用数列GID表示；
2.2.4：利用Matlab软件，编写程序，读取结构面数据文件st.dat，运行后自动生成岩体三维裂隙网络模型；
2.2.5：得到岩体三维裂隙网络模型；
2.3：二维裂隙网络模型切割
在三维裂隙网络模型的基础上，实现任意角度任意方位的二维裂隙网络模型切割，过程如下；
2.3.1：在三维裂隙网络模型上，结合Matlab软件编程工具，以三维裂隙网络模型中心点为中心，实现任意角度的剖面切割功能；
2.3.2：得到穿过中心点的任意角度的二维裂隙网络模型；
2.3.3：在三维裂隙网络模型上，结合Matlab软件编程工具，在三维裂隙网络模型任意位置上，实现任意角度和方位的剖面切割功能；
2.3.4：得到任意角度和方位的二维裂隙网络模型，如图1所示；
2.3.5：将切割剖面上的数据，保存到st1.dat文件中，此时剖面处于三维坐标系下，文件中数据格式自左向右依次为：节理中心点坐标(x，y，z)，节理长度D，倾角DA，倾向DD，走向SD，厚度thin，法向方向NX，NY，NZ；
2.3.6：将三维坐标系转化成二维坐标系，并将二维剖面数据保存到st2.dat文件中，文件中数据格式自左向右依次为：节理中心点坐标(x，y)，节理长度D，倾角DA，倾向DD，厚度thin，法向方向NX，NY，NZ；
2.4：裂隙网络模型数据输出
实现三维裂隙网络模型的数据输出和任意二维裂隙网络模型的数据输出，过程如下；
2.4.1：输出st.dat文件数据信息；
2.4.2：输出st1.dat文件数据信息；
2.4.3：输出st2.dat文件数据信息；
3)RQDt各向异性图绘制，过程如下：
3.1：RQDt求解计算
基于RQDt理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的RQDt值，过程如下：
3.1.1：RQDt理论公式如下：
式中：xi表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，RQDt代表对应阈值t的岩石质量指标，即阈值t下的RQDt值；
3.1.2确定二维裂隙网络模型的剖面中心点O、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，则中心点O坐标为：
边界方程为：
3.1.3：以O为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度L等于O点到裂隙网络模型边界的距离，用L0～L35表示，则测线方程为：
式中：s表示测线，α表示角度；
3.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(xa，ya)，将测线方程与边界方程依次连立，判断测线是否与边界相交，原理如下：
3.1.5：求出测线与裂隙网络模型的边界方程交点(XZ，α，YZ，α)；
3.1.6：确定测线所在区间，原理如下：
3.1.7：根据裂隙网络模型中每条节理的起点坐标(Xb，Yb)和终点坐标(Xc，Yc)，建立相应的解析方程，定义节理方程如下：
3.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(Xj，Yj)范围，如果交点符合a＜Xj＜c且b＜Yj＜d，则记录该交点，遍历所有节理方程求出所有m个交点；
3.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；
3.1.10：计算相邻交点的距离，公式如下：
x0＝a/2 (14)
y0＝b/2 (15)
xm+1＝XZ，α (16)
ym+1＝YZ，α (17)
3.1.11：输入一个阈值t；
3.1.12：循环比较di和t的大小，设初始lt＝0，规则如下：
lt＝lt+di，若di＞t (18)
3.1.13：求解每条测线对应的RQDt值，用mα表示，以及对应的测线起点与终点的距离lα，公式如下：
3.1.14：循环求出每条测线所对应的mα，获得36个方向上的RQDt值；
3.2：RQDt各向异性图绘制
根据求解出的36个RQDt值，绘制出RQDt的各向异性图，过程如下：
3.2.1：将36个RQDt值，按角度顺序依次排序；
3.2.2：以O点为圆心，以1为半径画圆，在射线角度为α上找到与圆心距离为
3.2.3：依次连接36组射线的端点，若某一条射线上的RQDt为0，则取圆心；
3.2.4：绘制出RQDt的各向异性图；
3.2.5：输出RQDt各向异性图；
3.3：不同阈值t下RQDt各向异性图绘制
用于绘制不同阈值t下的岩体RQDt各向异性图，过程如下：
3.3.1：输入不同的阈值t，求解出对应的RQDt值；
3.3.2：将不同阈值t下的RQDt值，以角度为变量，以RQDt值为函数，绘制到同个坐标系下；
3.3.3：得到不同阈值t下的RQDt各向异性图；
3.3.4：输出不同阈值t下的RQDt各向异性图，如图2所示；
3.4：输出RQDt的值；
4)基于BQ反演的最佳阈值t求解方法，过程如下：
4.1：基于BQ指标反演RQDt范围
结合BQ分级计算出的岩体质量级别，岩体为III级岩体，查找《岩石质量指标》表，见表1，反演确定出该岩体级别下，RQDt范围值，RQDt范围在50％～75％
表1；
4.2：最佳阈值t求解方法
4.2.1：在三维裂隙网络模型上，过中心点O，以任意角度剖切三个剖面，得到三个二维裂隙网络模型，导出二维裂隙网络模型及数据；
4.2.2：针对每一个二维裂隙网络模型，设置不同的阈值t，求解出不同阈值t下的RQDt值；
4.2.3：导出不同阈值t下的RQDt值，计算出RQDt均值；
4.2.4：以阈值t为横坐标，以RQDt的均值为纵坐标，绘制RQDt随阈值t变化的散点图；
4.2.5：根据散点图，设置拟合方程y＝aexp(bx)，拟合RQDt随阈值t变化的曲线图，如图3所示；
4.2.6：将反演出的RQDt范围值，带入到拟合出的RQDt随阈值t变化的曲线图中，结合函数方程和曲线图，求解出在该RQDt范围内，阈值t的范围，共得到三组阈值t的范围，见表2
表2；
4.2.7：针对三组阈值t的范围，取其范围的交集，以0.076m～0.124m作为最佳阈值t的范围；
4.2.8：以最佳阈值t范围的中点值作为最佳阈值t值，得到最佳阈值t，最佳阈值t＝0.1m；
4.2.9：输出最佳阈值t的范围和最佳阈值t值。
