CN112132414A - 一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法 - Google Patents

一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法 Download PDF

Info

Publication number
CN112132414A
CN112132414A CN202010893122.3A CN202010893122A CN112132414A CN 112132414 A CN112132414 A CN 112132414A CN 202010893122 A CN202010893122 A CN 202010893122A CN 112132414 A CN112132414 A CN 112132414A
Authority
CN
China
Prior art keywords
rqd
value
follows
network model
fracture network
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Withdrawn
Application number
CN202010893122.3A
Other languages
English (en)
Inventor
胡高建
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Shaoxing
Original Assignee
University of Shaoxing
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Shaoxing filed Critical University of Shaoxing
Priority to CN202010893122.3A priority Critical patent/CN112132414A/zh
Publication of CN112132414A publication Critical patent/CN112132414A/zh
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
    • G01C11/04—Interpretation of pictures
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/84—Systems specially adapted for particular applications
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F7/00—Methods or arrangements for processing data by operating upon the order or content of the data handled
    • G06F7/58—Random or pseudo-random number generators
    • G06F7/588—Random number generators, i.e. based on natural stochastic processes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06395—Quality analysis or management
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/20—Drawing from basic elements, e.g. lines or circles
    • G06T11/203—Drawing of straight lines or curves
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/20—Drawing from basic elements, e.g. lines or circles
    • G06T11/206—Drawing of charts or graphs
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Evolutionary Computation (AREA)
  • Strategic Management (AREA)
  • Economics (AREA)
  • Development Economics (AREA)
  • Educational Administration (AREA)
  • Computer Hardware Design (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Tourism & Hospitality (AREA)
  • Artificial Intelligence (AREA)
  • Game Theory and Decision Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Business, Economics & Management (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Biochemistry (AREA)
  • Software Systems (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Pathology (AREA)
  • Immunology (AREA)
  • Multimedia (AREA)
  • General Health & Medical Sciences (AREA)

Abstract

一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法，属于RQDt各向异性求解领域，其步骤包括：(1)结构面数字摄影测量快速获取；(2)结构面极点图和走向玫瑰花图绘制；(3)基于BQ指标的岩体质量计算；(4)岩体三维裂隙网络模型生成和剖切；(5)RQDt各向异性图绘制；(6)基于BQ反演的最佳阈值t求解方法；(7)RQDt各向异性求解方法。本发明将摄影测量、模糊等价分析、BQ理论、裂隙网络模型、广义RQD理论和反演计算相结合，实现了RQDt的最佳阈值t和各向异性的求解，提出了一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法。本发明方法明确，适用于岩体RQDt的各向异性求解。

Description

一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法
技术领域
本发明涉及一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法，特别的是本发明基于摄影测量、模糊等价分析、BQ指标、裂隙网络模型和广义RQD理论，反演计算出RQDt最佳阈值t的范围和最佳阈值t值，给出了RQDt各向异性的求解方法，提供了一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法，属于RQDt各向异性求解领域。
背景技术
在自然界中，各向异性无处不在。岩体表现出非常明显的各向异性，其性质随观察角度和测量方向的变化而变化。岩体各向异性变化的起源可归结于地质起源，如沉积岩中独特的层理、玄武岩中独具的天然孔洞、片麻岩中的片麻状构造和板岩中的劈理，以及贯穿这些岩体的断层等。除了固有的地质起源造成的岩体各向异性外，一些其它因素也会影响岩体的各向异性，如岩体中节理组的产状和分布形态，岩体中的节理会造成岩体的应力各向异性，影响岩体的质量。
岩体的外观和质量强烈的依赖于岩体的各向异性程度，岩体的稳定性是由岩体质量所决定的。一些地质现象也已经证明了这一点。例如，在巷道或高速公路隧道开挖过程中，会沿着主要节理组倾角出现非垂直性的变形破坏，或者隧道中的某一块岩体出现楔形块体掉落，进而造成整个岩体的失稳或者巷道的坍塌，这种失稳或坍塌对巷道开挖是非常不利的，不仅影响工程施工进度，更重要的是会造成严重的人员伤害或者财产安全。因此，研究并获得岩体质量的各向异性特征，具有重要的意义。
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值，研究的载体是裂隙网络模型。但在这个方面的研究中，主要是为了研究RQD的各向异性特征，对最佳阈值t的研究鲜有涉及。
在最佳阈值t的研究方面，有些学者也曾对最佳阈值t展开过探讨。如有的学者基于三维结构面网络模拟技术，运用分形理论分别计算其结构面分布的分形维数，通过不断改变广义RQD的阈值，得到不同阈值下的RQD值，将各剖面的分形维数与广义RQD值对比分析，为准确选取广义RQD的最佳阈值提供依据。有的学者基于修改后的块度指数MBi，建立了35种假想的三维裂隙网络模型，测量了具有不同阈值的广义RQD值，确定广义RQD的最佳阈值。这两个方面的最佳阈值t的研究，适用条件都有一定的局限性，都是在特定理论和背景下的一种最佳阈值t求解方法，当背景或者模型变化了之后，其最佳阈值t将不再具有最佳性。而且，由于分形维数或者块度指数本身并不具备表征岩体质量的功能，通过分形维数或者块度指数获得的最佳阈值t，能否反映真实的岩体质量，还有值得商榷的地方。
因此，找到并提出一种最能反映出岩体质量的RQDt最佳阈值t的求解方法，具有非常重要的科学和研究意义，也是RQDt阈值t研究中亟需解决的问题。
RQD的各向异性直接影响岩体质量，RQD各向异性对岩体质量的影响机理，目前还尚未探索清楚。由于在阈值t方面，还尚未有学者给出最佳阈值t的求解计算方法，因此也没有得到基于最佳阈值t的RQD各向异性计算公式，也尚未获得最能反映岩体质量的RQD各向异性求解方法。而获得RQD各向异性求解方法是研究岩体质量各向异性的基础和前提。
鉴于此，本发明提出了一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法。
发明内容
为了实现岩体RQDt各向异性的求解计算，本发明提供了一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法。本发明基于摄影测量、模糊等价分析、BQ指标、裂隙网络模型和广义RQD理论，反演计算出RQDt最佳阈值t的范围和最佳阈值t值，给出了RQDt各向异性的求解方法，提供了一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法。
为了解决上述技术问题，本发明提供如下的技术方案：
一种基于摄影测量、BQ和最佳阈值的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：根据结构面参数计算岩体完整性系数，公式如下：
Figure BSA0000218223910000031
式中：Jv为岩体体积节理数，单位条/m3；
3.2：Jv计算公式如下：
Figure BSA0000218223910000041
式中：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：基于最佳阈值t和剖切的三个二维裂隙网络模型，求解出最佳阈值t下的RQDt值，每个裂隙网络模型求解出36个RQDt值；
7.2：计算出每个裂隙网络模型上RQDt的最大值RQDtmax、最小值RQDtmin和均值
Figure BSA0000218223910000042
7.3：根据研究方向和RQDtmin方向位置关系，结合方差进行修正，提出各向异性条件下RQDt的计算公式如下：
Figure BSA0000218223910000044
式中：
Figure BSA0000218223910000043
为RQDt的均值，RQD′t为RQDt的修正系数；
7.4：修正系数RQD′t的求解考虑研究方向在某一具体方位角下的RQDtθ值和方差D，按如下方法修正，当RQDtθ＝RQDtmin时，RQD′t＝-D，D为方差；当RQDtθ＝RQDtmax时，RQD′t＝D；当
Figure BSA0000218223910000051
时，RQD′t＝0；
7.5：根据修正系数RQD′t的求解过程，提出RQD′t修正公式如下：
Figure BSA0000218223910000052
式中：RQDtθ为θ角下的RQDt值，RQDtmin为阈值t下的RQDt最小值，a，b为相关系数；
7.6：提出RQDt各向异性计算公式如下：
Figure BSA0000218223910000053
7.7：相关系数a和b的求解，过程如下：
7.7.1：计算出每个裂隙网络模型上RQDt的最大值RQDtmax、最小值RQDtmin和均值
Figure BSA0000218223910000054
7.7.2：根据RQDtθ与RQDtmax、RQDtmin和
Figure BSA0000218223910000055
以及方差的关系，求解出三组RQD′t值，并基于三组RQD′t值，绘制出散点图；
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求解计算
基于RQDt理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的RQDt值，过程如下：
5.1.1：RQDt理论公式如下：
Figure BSA0000218223910000061
式中：xi表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，RQDt代表对应阈值t的岩石质量指标，即阈值t下的RQDt值；
5.1.2确定二维裂隙网络模型的剖面中心点O、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，则中心点O坐标为：
Figure BSA0000218223910000062
Figure BSA0000218223910000063
边界方程为：
Figure BSA0000218223910000064
5.1.3：以O为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度L等于O点到裂隙网络模型边界的距离，用L0～L35表示，则测线方程为：
Figure BSA0000218223910000065
式中：s表示测线，α表示角度；
5.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(xa，ya)，将测线方程与边界方程依次连立，判断测线是否与边界相交，原理如下：
Figure BSA0000218223910000071
5.1.5：求出测线与裂隙网络模型的边界方程交点(XZ，α，YZ，α)；
5.1.6：确定测线所在区间，原理如下：
Figure BSA0000218223910000072
5.1.7：根据裂隙网络模型中每条节理的起点坐标(Xb，Yb)和终点坐标(Xc，Yc)，建立相应的解析方程，定义节理方程如下：
Figure BSA0000218223910000073
5.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(Xj，Yj)范围，如果交点符合a＜Xj＜c且b＜Yj＜d，则记录该交点，遍历所有节理方程求出所有m个交点；
5.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；
5.1.10：计算相邻交点的距离，公式如下：
Figure BSA0000218223910000074
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α，公式如下：
Figure BSA0000218223910000075
Figure BSA0000218223910000081
5.1.14：循环求出每条测线所对应的mα，获得36个方向上的RQDt值；
5.2：RQDt各向异性图绘制
根据求解出的36个RQDt值，绘制出RQDt的各向异性图，过程如下：
5.2.1：将36个RQDt值，按角度顺序依次排序；
5.2.2：以O点为圆心，以1为半径画圆，在射线角度为α上找到与圆心距离为
Figure BSA0000218223910000082
的点标出来；
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：实现三维裂隙网络模型的数据输出和任意二维裂隙网络模型的数据输出。
进一步，所述步骤(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：统计每组节理的数目
Figure BSA0000218223910000101
和每组节理的平均走向
Figure BSA0000218223910000102
2.6.3：根据作图的大小和各组节理数目，选取一定长度的线段代表一组节理，确定线段的比例尺LT；
2.6.4：以等于按比例尺LT表示的、数目最多的一组节理
Figure BSA0000218223910000103
的线段长度
Figure BSA0000218223910000104
为半径，作半圆，过圆心作南北线及东西线，在圆周上标明方位角；
2.6.5：对每组节理Ti，按平均走向
Figure BSA0000218223910000105
为方位角，在半圆上做出记号，自圆心向圆周记号点半径方向，按组内数目
Figure BSA0000218223910000106
和比例尺LT定出一点
Figure BSA0000218223910000107
该点即表示该组节理平均走向和节理数目；
2.6.6：顺次连接
Figure BSA0000218223910000108
和
Figure BSA0000218223910000109
如其中某组节理数为零，则连线回到圆心，再由圆心引出与下一组相连；
2.6.7：绘制出节理走向玫瑰花图。
本发明具有以下有益效果：
1、本发明基于摄影测量、模糊等价分析、BQ指标、裂隙网络模型和广义RQD理论，提供了一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法；
2、本发明实现了结构面快速摄影测量和模糊等价聚类分析；
3、本发明实现了岩体BQ指标分级；
4、本发明实现了三维裂隙网络模型生成和二维剖面模型生成；
5、本发明实现了RQDt求解计算和RQDt各向异性图的绘制；
6、本发明实现了RQDt的最佳阈值t求解；
7、本发明实现了RQDt各向异性的求解；
8、本发明方法明确，得到的最佳阈值t能够真实反映出岩体质量，得到的RQDt各向异性求解方法适用于求解RQDt的各向异性。
附图说明
图1是结构面极点图。
图2是二维裂隙网络图。
图3是不同阈值t下的RQDt各向异性图。
图4是RQDt随阈值t的拟合曲线。
图5是RQD′t的拟合曲线。
具体实施方式
下面参照附图对本发明做进一步说明。
参照图1～图5，一种基于摄影测量、BQ和最佳阈值的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为：
Figure BSA0000218223910000121
矩阵中的元素rij为第i个样本与第j个样本间的相似系数，表征其两者间的相似程度；rij越大表示样本i与样本j越相似；
2.2.2：计算相似系数rij：
Figure BSA0000218223910000122
式中，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中的元素
Figure BSA0000218223910000123
n级模糊关系矩阵R，即为n个R连乘；即
Figure BSA0000218223910000124
当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)
Figure BSA0000218223910000131
Figure BSA0000218223910000132
αd∈(0，360)，βd∈(0，90) (11)
2.3.3：得到以单位法向量表示的结构面产状数据；
2.3.4：基于结构面法向产状数据，根据结构面空间赤平投影图的纵剖面原理图，A’点为该平面法线的赤面投影；结合赤平投影原理图，计算出A’在赤平投影图上的坐标xn和yn，公式如下：
Figure BSA0000218223910000133
Figure BSA0000218223910000134
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：求解样本极差
Figure BSA0000218223910000141
Figure BSA0000218223910000142
2.4.3：计算每个分区区间Mm：
Figure BSA0000218223910000143
2.4.4：确定样本落在每个分区区间里的概率，先利用计算机循环语言统计落在每一个区间的样本个数Nm，结合样本总数N，计算样本数概率Pm：
Figure BSA0000218223910000144
2.4.5：求解样本均值
Figure BSA0000218223910000145
Figure BSA0000218223910000146
2.4.6：求解样本方差S2，其中S为标准差：
Figure BSA0000218223910000147
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：统计每组节理的数目
Figure BSA0000218223910000151
和每组节理的平均走向
Figure BSA0000218223910000152
Figure BSA0000218223910000153
i∈(1，10)∪(27，36) (23)
2.6.3：根据作图的大小和各组节理数目，选取一定长度的线段代表一组节理，确定线段的比例尺LT；
2.6.4：以等于按比例尺LT表示的、数目最多的一组节理
Figure BSA0000218223910000154
的线段长度
Figure BSA0000218223910000155
为半径，作半圆，过圆心作南北线及东西线，在圆周上标明方位角；
Figure BSA0000218223910000156
2.6.5：对每组节理Ti，按平均走向
Figure BSA0000218223910000157
为方位角，在半圆上做出记号，自圆心向圆周记号点半径方向，按组内数目
Figure BSA0000218223910000158
和比例尺LT定出一点
Figure BSA0000218223910000159
该点即表示该组节理平均走向和节理数目；
Figure BSA00002182239100001510
2.6.6：顺次连接
Figure BSA00002182239100001511
和
Figure BSA00002182239100001512
如其中某组节理数为零，则连线回到圆心，再由圆心引出与下一组相连；
2.6.7：绘制出节理走向玫瑰花图；
3)基于BQ指标的岩体质量计算，过程如下：
3.1：根据结构面参数计算岩体完整性系数，公式如下：
Figure BSA00002182239100001513
式中：Jv为岩体体积节理数，单位条/m3；
3.2：Jv计算公式如下：
Figure BSA00002182239100001514
式中：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.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)RQDt各向异性图绘制，过程如下：
5.1：RQDt求解计算
基于RQDt理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的RQDt值，过程如下：
5.1.1：RQDt理论公式如下：
Figure BSA0000218223910000171
式中：xi表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，RQDt代表对应阈值t的岩石质量指标，即阈值t下的RQDt值；
5.1.2确定二维裂隙网络模型的剖面中心点O、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，则中心点O坐标为：
Figure BSA0000218223910000181
Figure BSA0000218223910000182
边界方程为：
Figure BSA0000218223910000183
5.1.3：以O为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度L等于O点到裂隙网络模型边界的距离，用L0～L35表示，则测线方程为：
Figure BSA0000218223910000184
式中：s表示测线，α表示角度；
5.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(xa，ya)，将测线方程与边界方程依次连立，判断测线是否与边界相交，原理如下：
Figure BSA0000218223910000185
5.1.5：求出测线与裂隙网络模型的边界方程交点(XZ，α，YZ，α)；
5.1.6：确定测线所在区间，原理如下：
Figure BSA0000218223910000186
5.1.7：根据裂隙网络模型中每条节理的起点坐标(Xb，Yb)和终点坐标(Xc，Yc)，建立相应的解析方程，定义节理方程如下：
Figure BSA0000218223910000187
5.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(Xj，Yj)范围，如果交点符合a＜Xj＜c且b＜Yj＜d，则记录该交点，遍历所有节理方程求出所有m个交点；
5.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；
5.1.10：计算相邻交点的距离，公式如下：
Figure BSA0000218223910000191
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α，公式如下：
Figure BSA0000218223910000192
Figure BSA0000218223910000193
5.1.14：循环求出每条测线所对应的mα，获得36个方向上的RQDt值；
5.2：RQDt各向异性图绘制
根据求解出的36个RQDt值，绘制出RQDt的各向异性图，过程如下：
5.2.1：将36个RQDt值，按角度顺序依次排序；
5.2.2：以O点为圆心，以1为半径画圆，在射线角度为α上找到与圆心距离为
Figure BSA0000218223910000194
的点标出来；
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％
Figure BSA0000218223910000201
表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
Figure BSA0000218223910000202
表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和均值
Figure BSA0000218223910000215
7.3：根据研究方向和RQDtmin方向位置关系，结合方差进行修正，提出各向异性条件下RQDt的计算公式如下：
Figure BSA0000218223910000216
式中：
Figure BSA0000218223910000217
为RQDt的均值，RQD′t为RQDt的修正系数；
7.4：修正系数RQD′t的求解考虑研究方向在某一具体方位角下的RQDtθ值和方差D，按如下方法修正，当RQDtθ＝RQDtmin时，RQD′t＝-D，D为方差；当RQDtθ＝RQDtmax时，RQD′t＝D；当
Figure BSA0000218223910000218
时，RQD′t＝0；
7.5：根据修正系数RQD′t的求解过程，提出RQD′t修正公式如下：
Figure BSA0000218223910000211
式中：RQDtθ为θ角下的RQDt值，RQDtmin为阈值t下的RQDt最小值，a，b为相关系数；
7.6：提出RQDt各向异性计算公式如下：
Figure BSA0000218223910000212
7.7：相关系数a和b的求解，过程如下：
7.7.1：计算出每个裂隙网络模型上RQDt的最大值RQDtmax、最小值RQDtmin和均值
Figure BSA0000218223910000219
7.7.2：根据RQDtθ与RQDtmax、RQDtmin和
Figure BSA00002182239100002110
以及方差的关系，求解出三组RQD′t值，并基于三组RQD′t值，绘制出散点图；
7.7.3：根据散点图做出拟合曲线，如图5所示，曲线截距即为a值，斜率即为b值，曲线方程即为RQD′t的修正公式，得到的RQD′t公式如下：
Figure BSA0000218223910000213
7.8：将求解出的a、b值，带入到修正系数公式及RQDt各向异性计算公式中，得到与角度θ有关的RQDt各向异性公式：
Figure BSA0000218223910000214
7.9：根据RQDt各向异性公式，求解出任意角度的RQDt值。

Claims (6)

1.一种基于摄影测量、BQ和最佳阈值的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：根据结构面参数计算岩体完整性系数，公式如下：
Figure FSA0000218223900000011
式中：Jv为岩体体积节理数，单位条/m3；
3.2：Jv计算公式如下：
Figure FSA0000218223900000012
式中：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各向异性求解方法。
2.如权利要求1所述的一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法，其特征在于，所述步骤(7)中，RQDt各向异性求解方法的过程如下：
7.1：基于最佳阈值t和剖切的三个二维裂隙网络模型，求解出最佳阈值t下的RQDt值，每个裂隙网络模型求解出36个RQDt值；
7.2：计算出每个裂隙网络模型上RQDt的最大值RQDtmax、最小值RQDtmin和均值
Figure FSA0000218223900000025
7.3：根据研究方向和RQDtmin方向位置关系，结合方差进行修正，提出各向异性条件下RQDt的计算公式如下：
Figure FSA0000218223900000026
式中：
Figure FSA0000218223900000021
为RQDt的均值，RQD′t为RQDt的修正系数；
7.4：修正系数RQD′t的求解考虑研究方向在某一具体方位角下的RQDtθ值和方差D，按如下方法修正，当RQDtθ＝RQDtmin时，RQD′t＝-D，D为方差；当RQDtθ＝RQDtmax时，RQD′t＝D；当
Figure FSA0000218223900000022
时，RQD′t＝0；
7.5：根据修正系数RQD′t的求解过程，提出RQD′t修正公式如下：
Figure FSA0000218223900000023
式中：RQDtθ为θ角下的RQDt值，RQDtmin为阈值t下的RQDt最小值，a，b为相关系数；
7.6：提出RQDt各向异性计算公式如下：
Figure FSA0000218223900000024
7.7：相关系数a和b的求解，过程如下：
7.7.1：计算出每个裂隙网络模型上RQDt的最大值RQDtmax、最小值RQDtmin和均值
Figure FSA0000218223900000031
7.7.2：根据RQDtθ与RQDtmax、RQDtmin和
Figure FSA0000218223900000032
以及方差的关系，求解出三组RQD′t值，并基于三组RQD′t值，绘制出散点图；
7.7.3：根据散点图做出拟合曲线，曲线截距即为a值，斜率即为b值，曲线方程即为RQD′t的修正公式；
7.8：将求解出的a、b值，带入到修正系数公式及RQDt各向异性计算公式中，得到与角度θ有关的RQDt各向异性公式；
7.9：根据RQDt各向异性公式，求解出任意角度的RQDt值。
3.如权利要求1所述的一种基于摄影测量、BQ和最佳阈值的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值。
4.如权利要求1所述的一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法，其特征在于，所述步骤(5)中，RQDt各向异性图绘制的过程如下：
5.1：RQDt求解计算
基于RQDt理论，确定出二维裂隙网络模型的中心点O，以每隔10°角度做射线，穿过岩体裂隙网络模型，求解出岩体36个方位的RQDt值，过程如下：
5.1.1：RQDt理论公式如下：
Figure FSA0000218223900000041
式中：xi表示沿某一测线方向的第i个大于给定阈值t的整段岩石或间距长度，RQDt代表对应阈值t的岩石质量指标，即阈值t下的RQDt值；
5.1.2确定二维裂隙网络模型的剖面中心点O、长度a和宽度b，以模型的左下角为坐标原点，水平向右为x轴，垂直向上为y轴，则中心点O坐标为：
Figure FSA0000218223900000042
Figure FSA0000218223900000043
边界方程为：
Figure FSA0000218223900000044
5.1.3：以O为起点，每隔α＝10°角绘制一条测线，与裂隙网络模型相交，共绘制36条测线，测线长度L等于O点到裂隙网络模型边界的距离，用L0～L35表示，则测线方程为：
Figure FSA0000218223900000045
式中：s表示测线，α表示角度；
5.1.4：判断测线与边界的交点，设测线方程与边界方程的交点为(xa，ya)，将测线方程与边界方程依次连立，判断测线是否与边界相交，原理如下：
Figure FSA0000218223900000046
5.1.5：求出测线与裂隙网络模型的边界方程交点(XZ，α，YZ，α)；
5.1.6：确定测线所在区间，原理如下：
Figure FSA0000218223900000047
5.1.7：根据裂隙网络模型中每条节理的起点坐标(Xb，Yb)和终点坐标(Xc，Yc)，建立相应的解析方程，定义节理方程如下：
Figure FSA0000218223900000051
5.1.8：求解第一条测线与各节理方程的交点，循环判断每个交点(Xj，Yj)范围，如果交点符合a＜Xj＜c且b＜Yj＜d，则记录该交点，遍历所有节理方程求出所有m个交点；
5.1.9：将记录的m个交点及测线起点坐标和终点坐标，按横坐标或纵坐标从小到大排序；
5.1.10：计算相邻交点的距离，公式如下：
Figure FSA0000218223900000052
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α，公式如下：
Figure FSA0000218223900000053
Figure FSA0000218223900000054
5.1.14：循环求出每条测线所对应的mα，获得36个方向上的RQDt值；
5.2：RQDt各向异性图绘制
根据求解出的36个RQDt值，绘制出RQDt的各向异性图，过程如下：
5.2.1：将36个RQDt值，按角度顺序依次排序；
5.2.2：以O点为圆心，以l为半径画圆，在射线角度为α上找到与圆心距离为
Figure FSA0000218223900000055
的点标出来；
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的值。
5.如权利要求1所述的一种基于摄影测量、BQ和最佳阈值的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：实现三维裂隙网络模型的数据输出和任意二维裂隙网络模型的数据输出。
6.如权利要求1所述的一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法，其特征在于，所述步骤(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：统计每组节理的数目
Figure FSA0000218223900000081
和每组节理的平均走向
Figure FSA0000218223900000082
2.6.3：根据作图的大小和各组节理数目，选取一定长度的线段代表一组节理，确定线段的比例尺LT；
2.6.4：以等于按比例尺LT表示的、数目最多的一组节理
Figure FSA0000218223900000083
的线段长度
Figure FSA0000218223900000084
为半径，作半圆，过圆心作南北线及东西线，在圆周上标明方位角；
2.6.5：对每组节理Ti，按平均走向
Figure FSA0000218223900000085
为方位角，在半圆上做出记号，自圆心向圆周记号点半径方向，按组内数目
Figure FSA0000218223900000086
和比例尺LT定出一点
Figure FSA0000218223900000087
该点即表示该组节理平均走向和节理数目；
2.6.6：顺次连接
Figure FSA0000218223900000088
和
Figure FSA0000218223900000089
如其中某组节理数为零，则连线回到圆心，再由圆心引出与下一组相连；
2.6.7：绘制出节理走向玫瑰花图。
CN202010893122.3A 2020-08-25 2020-08-25 一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法 Withdrawn CN112132414A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010893122.3A CN112132414A (zh) 2020-08-25 2020-08-25 一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010893122.3A CN112132414A (zh) 2020-08-25 2020-08-25 一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法

Publications (1)

Publication Number Publication Date
CN112132414A true CN112132414A (zh) 2020-12-25

Family

ID=73848693

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010893122.3A Withdrawn CN112132414A (zh) 2020-08-25 2020-08-25 一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法

Country Status (1)

Country Link
CN (1) CN112132414A (zh)

Similar Documents

Publication Publication Date Title
CN112150001A (zh) 基于摄影测量、BQ、改进Mathews稳定图的围岩稳定性评价方法
CN112464513A (zh) 一种基于摄影测量和RQD反演的RQDt最佳阈值t求解方法
CN112200417A (zh) 基于摄影测量、BQ、RQDt、地应力的改进Mathews稳定图评价方法
CN112200418A (zh) 基于摄影测量、bq、数值模拟的围岩稳定性动态评价方法
CN112132407A (zh) 一种基于BQ反演最佳阈值t的空间RQDt求解方法
CN114088015A (zh) 一种岩体三维裂隙网络模型快速智能生成和剖切方法
CN112200426A (zh) 基于激光扫描、bq、数值模拟的围岩稳定性动态评价方法
CN112464514A (zh) 一种基于摄影测量、RQD和RQDt的巷道开挖不利方位求解方法
CN112132403A (zh) 一种基于摄影测量和BQ反演的RQDt最佳阈值t求解方法
CN112132410A (zh) 一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法
CN112132411A (zh) 基于激光扫描、BQ、RQDt各向异性的Q各向异性求解方法
CN112464516A (zh) 一种基于激光扫描和RQD反演最佳阈值t的空间RQDt求解方法
CN112132406A (zh) 一种基于激光扫描、BQ和RQDt的巷道开挖不利方位求解方法
CN112132408A (zh) 一种基于激光扫描和BQ反演最佳阈值t的空间RQDt求解方法
CN112132414A (zh) 一种基于摄影测量、BQ和最佳阈值的RQDt各向异性求解方法
CN112464506A (zh) 一种基于摄影测量、RQD和最佳阈值的RQDt各向异性求解方法
CN112149996A (zh) 基于激光扫描、BQ、RQDt各向异性的改进Mathews稳定图评价方法
CN112132412A (zh) 一种基于摄影测量、BQ和RQDt的巷道开挖不利方位求解方法
CN112132404A (zh) 基于摄影测量、BQ、RQDt各向异性的Q各向异性求解方法
CN112150002A (zh) 基于激光扫描、BQ、RQDt、地应力的改进Mathews稳定图评价方法
CN112464507A (zh) 基于摄影测量、RQD、RQDt各向异性的Q各向异性求解方法
CN113324513A (zh) 一种基于摄影测量和RQD反演最佳阈值t的空间RQDt求解方法
CN112200431A (zh) 基于激光扫描、bq、数值模拟的空区稳定性动态评价方法
CN112132405A (zh) 一种基于激光扫描、BQ和最佳阈值的RQDt各向异性求解方法
CN112464515A (zh) 一种基于激光扫描、RQD和最佳阈值的RQDt各向异性求解方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
WW01 Invention patent application withdrawn after publication
WW01 Invention patent application withdrawn after publication

Application publication date: 20201225