CN112200431A - 基于激光扫描、bq、数值模拟的空区稳定性动态评价方法 - Google Patents
基于激光扫描、bq、数值模拟的空区稳定性动态评价方法 Download PDFInfo
- Publication number
- CN112200431A CN112200431A CN202011011892.7A CN202011011892A CN112200431A CN 112200431 A CN112200431 A CN 112200431A CN 202011011892 A CN202011011892 A CN 202011011892A CN 112200431 A CN112200431 A CN 112200431A
- Authority
- CN
- China
- Prior art keywords
- rqd
- stability
- follows
- value
- rock
- 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
Links
- 238000011156 evaluation Methods 0.000 title claims abstract description 65
- 238000004088 simulation Methods 0.000 title claims abstract description 44
- 239000011435 rock Substances 0.000 claims abstract description 215
- 238000000034 method Methods 0.000 claims abstract description 170
- 238000010586 diagram Methods 0.000 claims abstract description 59
- 238000004458 analytical method Methods 0.000 claims abstract description 41
- 230000006872 improvement Effects 0.000 claims abstract description 6
- 230000006641 stabilisation Effects 0.000 claims abstract description 4
- 238000011105 stabilization Methods 0.000 claims abstract description 4
- 230000008569 process Effects 0.000 claims description 72
- 238000004364 calculation method Methods 0.000 claims description 34
- 238000012937 correction Methods 0.000 claims description 31
- 239000013598 vector Substances 0.000 claims description 29
- 230000006378 damage Effects 0.000 claims description 26
- 238000009826 distribution Methods 0.000 claims description 19
- 230000005484 gravity Effects 0.000 claims description 19
- 230000006870 function Effects 0.000 claims description 16
- 238000013461 design Methods 0.000 claims description 12
- 230000009467 reduction Effects 0.000 claims description 12
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 12
- 238000005192 partition Methods 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 9
- 238000007621 cluster analysis Methods 0.000 claims description 8
- 238000004422 calculation algorithm Methods 0.000 claims description 7
- 238000007619 statistical method Methods 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 5
- 238000005065 mining Methods 0.000 claims description 5
- 238000013507 mapping Methods 0.000 claims description 4
- 238000005259 measurement Methods 0.000 claims description 4
- 238000007476 Maximum Likelihood Methods 0.000 claims description 3
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 3
- 230000004075 alteration Effects 0.000 claims description 3
- 238000010924 continuous production Methods 0.000 claims description 3
- 238000012938 design process Methods 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 3
- 238000012067 mathematical method Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000011084 recovery Methods 0.000 claims description 3
- 230000001105 regulatory effect Effects 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000009827 uniform distribution Methods 0.000 claims description 3
- 238000011160 research Methods 0.000 description 8
- 238000005553 drilling Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 3
- 230000001276 controlling effect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000007477 logistic regression Methods 0.000 description 2
- 239000011295 pitch Substances 0.000 description 2
- 206010057175 Mass conditions Diseases 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001149 cognitive effect Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000012854 evaluation process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
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
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/24—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
- G01B11/25—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures by projecting a pattern, e.g. one or more lines, moiré fringes on the object
- G01B11/2518—Projection by scanning of the object
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N3/00—Investigating strength properties of solid materials by application of mechanical stress
- G01N3/08—Investigating strength properties of solid materials by application of mechanical stress by applying steady tensile or compressive forces
-
- 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
- G01N2203/02—Details not specific for a particular testing method
- G01N2203/025—Geometry of the test
- G01N2203/0252—Monoaxial, i.e. the forces being applied along a single axis of the specimen
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
- G01N2203/02—Details not specific for a particular testing method
- G01N2203/06—Indicating or recording means; Sensing means
- G01N2203/0641—Indicating or recording means; Sensing means using optical, X-ray, ultraviolet, infrared or similar detectors
- G01N2203/0647—Image analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Business, Economics & Management (AREA)
- Human Resources & Organizations (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Geometry (AREA)
- Educational Administration (AREA)
- Economics (AREA)
- Evolutionary Computation (AREA)
- Strategic Management (AREA)
- Entrepreneurship & Innovation (AREA)
- Development Economics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Computer Hardware Design (AREA)
- Artificial Intelligence (AREA)
- Operations Research (AREA)
- General Business, Economics & Management (AREA)
- Tourism & Hospitality (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Mathematical Analysis (AREA)
- Game Theory and Decision Science (AREA)
- Computational Mathematics (AREA)
- Marketing (AREA)
- Quality & Reliability (AREA)
- Chemical & Material Sciences (AREA)
- Structural Engineering (AREA)
- Civil Engineering (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Health & Medical Sciences (AREA)
- Architecture (AREA)
- Analytical Chemistry (AREA)
Abstract
一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，属于采空区稳定性评价领域，其步骤包括：(1)结构面三维激光扫描快速获取；(2)结构面聚类分析；(3)基于BQ指标的岩体质量计算；(4)岩体三维裂隙网络模型生成和剖切；(5)RQDt各向异性图绘制；(6)基于BQ反演的最佳阈值t求解方法；(7)RQDt各向异性求解方法；(8)Mathews稳定图法的改进方法；(9)改进Mathews稳定图评价方法；(10)围岩稳定性动态分析方法；(11)空区稳定性动态评价方法。本发明实现了RQDt的最佳阈值t、RQDt各向异性的求解和Mathews稳定图法的改进，实现了改进Mathews稳定图法和数值模拟相结合的空区围岩稳定性动态评价。本发明方法明确，适用于采空区围岩稳定性动态评价。
Description
技术领域
本发明涉及一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，特别的是本发明基于三维激光扫描、近邻传播算法、BQ指标、裂隙网络模型和广义RQD理论，反演计算出RQDt最佳阈值t的范围和最佳阈值t值，给出了RQDt各向异性的求解方法，并基于RQDt的各向异性，给出了Mathews稳定图法的改进方法和改进后的Mathews稳定图评价方法，再结合数值模拟方法，提供了一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，属于采空区稳定性评价领域。
背景技术
金属矿山经过长时间地下开采，会形成大规模的采空区群，采空区群在地下相互贯通、重叠和集中，彼此间应力扰动，相互制约，形成了复杂的群体空间关系，岩体稳定状态处于动态稳定平衡中。采空区群由于岩体结构的特殊性和复杂性，以及在采空区群回采过程中的随机性、不确定性、连续性和规模性，由采空区群失稳坍塌造成的地质灾害非常严重、破坏力强悍，造成的人员财产损失巨大。
采空区群的破坏是一个动态的破坏过程。首先，由于空区群岩体结构的特殊性，采空区的破坏往往是成片成规模的，当一个采空区发生失稳破坏后，由于岩体结构的改变，相邻的空区会受到影响而产生链式破坏，进而造成空区群的规模性破坏。其次，由于空区群在回采时的连续性和规模性，往往是多个采空区同时回采，多处扰动同时存在，加剧了空区群岩体结构的动态破坏过程，也增加了岩体失稳破坏的风险。由于采空区群的这种复杂性、多变性、突发性和规模性的特征，人们对其破坏发生发展的过程仍缺乏足够的认知了解，许多关键的科学问题尚未得到有效的解决，也没有形成一个有效的空区稳定性动态评价方法。
数值模拟方法在研究围岩稳定性及应力演化规律方面具有显著优势，通过对任意尺寸形态的地质模型的模拟研究，可以实现围岩稳定性的动态分析。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、数值模拟的空区稳定性动态评价方法。
发明内容
为了实现采空区稳定性的动态评价，本发明提供了一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法。本发明基于三维激光扫描、结构面聚类分析、BQ指标、裂隙网络模型和广义RQD理论，反演计算出RQDt最佳阈值t的范围和最佳阈值t值，给出了RQDt各向异性的求解方法，并基于RQDt的各向异性，给出了Mathews稳定图法的改进方法和改进后的Mathews稳定图评价方法，再结合数值模拟方法，提供了一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法。
为了解决上述技术问题，本发明提供如下的技术方案：
一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，所述方法包括以下步骤：
(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稳定图法的改进方法；
(9)改进Mathews稳定图评价方法；
(10)围岩稳定性动态分析方法，过程如下：
10.1：根据矿房平面图和剖面图，利用CAD软件，建立多个连续矿房的三维模型；
10.2：将三维模型导入到有限元分析软件中，赋予三维模型材料参数，包括岩体强度、密度、弹性模量、泊松比、粘聚力、摩擦角和容重；
10.3：赋予三维模型边界条件和初始应力；
10.4：开展矿房动态回采的数值模拟，过程如下：
10.4.1：开展矿房未回采状态的数值模拟，得到初始应力场；
10.4.2：开展矿房动态回采的数值模拟，得到矿房动态回采过程中，矿房及采空区的动态应力场；
10.4.3：矿房的动态回采过程是指矿房自第一个矿房回采开始至所有矿房回采结束的连续过程，矿房每动态回采一次，得到一组动态矿房回采的数值模拟结果，所有矿房回采完成后，得到采空区的应力场；
10.5：围岩稳定性的动态分析，过程如下：
10.5.1：分析区域为矿房顶板、上下盘围岩和矿柱；
10.5.2：岩体剪切破坏的判断准则选用以主应力形式表示的Mohr-Coulomb屈服准则，当fs＞0时，岩体发生剪切破坏，公式如下：
式中：σ1、σ3分别为最大主应力和最小主应力；c为粘聚力；
为摩擦角；
10.5.3：岩体张拉破坏的判断准则，根据岩体所受拉应力σt与抗拉强度Rt之间的关系确定，当σt≥Rt时，F≥0，即认为岩体发生张拉破坏，公式如下：
F＝σt-Rt
式中：σt为岩体所受拉应力；Rt为岩体抗拉强度；
10.5.4：根据分析结果，得到每次矿房动态回采后，矿房顶板、上下盘围岩和矿柱的破坏模式；
10.5.5：将每次矿房动态回采后得到的矿房顶板、上下盘围岩和矿柱的破坏模式汇总，得到围岩的稳定性动态分析结果；
(11)空区稳定性动态评价方法。
进一步，所述步骤(11)中，空区稳定性动态评价方法的过程如下：
11.1：以改进Mathews稳定图法得到的围岩稳定性评价结果作为采空区稳定性的初始评价结果；
11.2：利用数值模拟得到的围岩稳定性动态分析结果，对初始评价结果进行修正，过程如下：
11.2.1：当改进Mathews稳定图法评价的采空区围岩处于崩落区时，则该采空区处于崩落状态；
11.2.2：当改进Mathews稳定图法评价的采空区围岩处于破坏区时，则该采空区处于破坏状态；
11.2.3：当改进Mathews稳定图法评价的采空区围岩处于稳定区时，若在围岩稳定性动态分析中得到的采空区顶板、间柱、上下盘围岩未发生破坏，则该采空区处于稳定状态，若在围岩稳定性动态分析中得到的采空区顶板、间柱、上下盘围岩发生破坏，则该采空区有可能处于破坏状态；
11.3：得到采空区的稳定性动态评价结果。
进一步，所述步骤(9)中，改进Mathews稳定图评价方法的过程如下：
9.1：选取用于开展上下盘围岩稳定性评价的矿房；
9.2：测量出矿房上下盘围岩宽度X和长度Y，计算出岩石应力系数A；节理产状调整系数B和重力调整系数C；
9.3：测量并求解出矿房上下盘的节理组数Jn；节理粗糙度系数Jr；节理蚀变系数Ja；节理水折减系数Jw和SRF为应力折减系数；
9.4：求解出最佳阈值t下的RQDtθ、RQDtmin和
计算出相关拟合系数a和b，求解出考虑RQDt各向异性特征的Q′值；
9.5：计算出矿房稳定数N和水力半径R；
9.6：计算出稳定区-破坏边界的曲线公式和破坏-崩落边界的曲线公式，并绘制到改进后的Mathews稳定性图上；
9.7：将计算出的矿房稳定数N和水力半径R绘制到改进后的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由完整岩石单轴抗压强度与采场中线采矿产生的地应力的比值确定，公式如下：
当
A＝0.1
当
当
A＝1
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最小值；
为RQDt的均值；a，b为相关拟合系数；Jn为节理组数；Jr为节理粗糙度系数；Ja为节理蚀变系数；Jw为节理水折减系数；SRF为应力折减系数；
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的计算公式如下：
式中：
为RQDt的均值，RQD′t为RQDt的修正系数；
7.4：修正系数RQD′t的求解考虑研究方向在某一具体方位角下的RQDtθ值和方差D，按如下方法修正，当RQDtθ＝RQDtmin时，RQD′t＝-D，D为方差；当RQDtθ＝RQDtmax时，RQD′t＝D；当
时，RQD′t＝0；
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和
以及方差的关系，求解出三组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求解计算，过程如下：
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：利用Monte Carlo方法是根据确立的结构面几何概率模型，再现出服从这种模型的结构面网络模型，在(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理论，RQDt各向异性求解方法、改进Mathews稳定图方法和数值模拟方法，提供了一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法；
2、本发明实现了结构面快速三维激光扫描和近邻传播算法聚类分析；
3、本发明实现了岩体BQ指标分级；
4、本发明实现了三维裂隙网络模型生成和二维剖面模型生成；
5、本发明实现了RQDt求解计算和RQDt各向异性图的绘制；
6、本发明实现了RQDt的最佳阈值t求解；
7、本发明实现了RQDt各向异性的求解；
8、本发明实现了考虑RQDt的各向异性特征的Mathews稳定图法的改进；
9、本发明实现了改进Mathews稳定图法的空区围岩稳定性评价；
10、本发明实现了基于数值模拟的围岩稳定性动态分析；
11、本发明方法明确，提供的空区稳定性动态评价方法实现了采空区稳定性的动态评价。
附图说明
图1是结构面极点图。
图2是二维裂隙网络图。
图3是不同阈值t下的RQDt各向异性图。
图4是RQDt随阈值t的拟合曲线。
图5是RQD′t的拟合曲线。
图6是矿房围岩稳定性评价结果图。
具体实施方式
下面参照附图对本发明做进一步说明。
参照图1～图6，一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，包括以下步骤：
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：利用Monte Carlo方法是根据确立的结构面几何概率模型，再现出服从这种模型的结构面网络模型，在(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的计算公式如下：
式中：
为RQDt的均值，RQD′t为RQDt的修正系数；
7.4：修正系数RQD′t的求解考虑研究方向在某一具体方位角下的RQDtθ值和方差D，按如下方法修正，当RQDtθ＝RQDtmin时，RQD′t＝-D，D为方差；当RQDtθ＝RQDtmax时，RQD′t＝D；当
时，RQD′t＝0；
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和
以及方差的关系，求解出三组RQD′t值，并基于三组RQD′t值，绘制出散点图；
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由完整岩石单轴抗压强度与采场中线采矿产生的地应力的比值确定，公式如下：
当
A＝0.1 (40)
当
当
A＝1 (42)
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最小值；
为RQDt的均值；a，b为相关拟合系数；Jn为节理组数；Jr为节理粗糙度系数；Ja为节理蚀变系数；Jw为节理水折减系数；SRF为应力折减系数；
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％的概率崩落；处在崩落区的工程，将持续发生崩落；
9)改进Mathews稳定图评价方法，过程如下：
9.1：选取用于开展上下盘围岩稳定性评价的矿房；
9.2：测量出矿房上下盘围岩宽度X和长度Y，计算出岩石应力系数A；节理产状调整系数B和重力调整系数C；
9.3；测量并求解出矿房上下盘的节理组数Jn；节理粗糙度系数Jr；节理蚀变系数Ja；节理水折减系数Jw和SRF为应力折减系数；
9.4：求解出最佳阈值t下的RQDtθ、RQDtmin和
计算出相关拟合系数a和b，求解出考虑RQDt各向异性特征的Q′值；
9.5：计算出矿房稳定数N和水力半径R；
9.6：计算出稳定区-破坏边界的曲线公式和破坏-崩落边界的曲线公式，并绘制到改进后的Mathews稳定性图上；
9.7：将计算出的矿房稳定数N和水力半径R绘制到改进后的Mathews稳定图上，如图6所示，根据绘制的散点所处的区域，给出矿房上下盘围岩的稳定性评价结果，3203矿房和3204矿房处于崩落区，矿房围岩评价结果为崩落，3205矿房处于破坏或主要破坏区，围岩破坏的概率为55％，围岩崩落的概率为33％，围岩稳定的概率为12％；
10)围岩稳定性动态分析方法，过程如下：
10.1：根据矿房平面图和剖面图，利用CAD软件，建立多个连续矿房的三维模型；
10.2：将三维模型导入到有限元分析软件中，赋予三维模型材料参数，包括岩体强度、密度、弹性模量、泊松比、粘聚力、摩擦角和容重；
10.3：赋予三维模型边界条件和初始应力；
10.4：开展矿房动态回采的数值模拟，过程如下：
10.4.1：开展矿房未回采状态的数值模拟，得到初始应力场；
10.4.2：开展矿房动态回采的数值模拟，得到矿房动态回采过程中，矿房及采空区的动态应力场；
10.4.3：矿房的动态回采过程是指矿房自第一个矿房回采开始至所有矿房回采结束的连续过程，矿房每动态回采一次，得到一组动态矿房回采的数值模拟结果，所有矿房回采完成后，得到采空区的应力场；
10.5：围岩稳定性的动态分析，过程如下：
10.5.1：分析区域为矿房顶板、上下盘围岩和矿柱；
10.5.2：岩体剪切破坏的判断准则选用以主应力形式表示的Mohr-Coulomb屈服准则，当fs＞0时，岩体发生剪切破坏，公式如下：
式中：σ1、σ3分别为最大主应力和最小主应力；c为粘聚力；
为摩擦角；
10.5.3：岩体张拉破坏的判断准则，根据岩体所受拉应力σt与抗拉强度Rt之间的关系确定，当σt≥Rt时，F≥0，即认为岩体发生张拉破坏，公式如下：
F＝σt-Rt (48)
式中：σt为岩体所受拉应力；Rt为岩体抗拉强度；
10.5.4：根据分析结果，得到每次矿房动态回采后，矿房顶板、上下盘围岩和矿柱的破坏模式；
10.5.5：将每次矿房动态回采后得到的矿房顶板、上下盘围岩和矿柱的破坏模式汇总，得到围岩的稳定性动态分析结果；
11)空区稳定性动态评价方法，过程如下：
11.1：以改进Mathews稳定图法得到的围岩稳定性评价结果作为采空区稳定性的初始评价结果；
11.2：利用数值模拟得到的围岩稳定性动态分析结果，对初始评价结果进行修正，过程如下：
11.2.1：当改进Mathews稳定图法评价的采空区围岩处于崩落区时，则该采空区处于崩落状态；
11.2.2：当改进Mathews稳定图法评价的采空区围岩处于破坏区时，则该采空区处于破坏状态；
11.2.3：当改进Mathews稳定图法评价的采空区围岩处于稳定区时，若在围岩稳定性动态分析中得到的采空区顶板、间柱、上下盘围岩未发生破坏，则该采空区处于稳定状态，若在围岩稳定性动态分析中得到的采空区顶板、间柱、上下盘围岩发生破坏，则该采空区有可能处于破坏状态；
11.3：得到采空区的稳定性动态评价结果。
Claims (9)
1.一种基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述方法包括以下步骤：
(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稳定图法的改进方法；
(9)改进Mathews稳定图评价方法；
(10)围岩稳定性动态分析方法，过程如下：
10.1：根据矿房平面图和剖面图，利用CAD软件，建立多个连续矿房的三维模型；
10.2：将三维模型导入到有限元分析软件中，赋予三维模型材料参数，包括岩体强度、密度、弹性模量、泊松比、粘聚力、摩擦角和容重；
10.3：赋予三维模型边界条件和初始应力；
10.4：开展矿房动态回采的数值模拟，过程如下：
10.4.1：开展矿房未回采状态的数值模拟，得到初始应力场；
10.4.2：开展矿房动态回采的数值模拟，得到矿房动态回采过程中，矿房及采空区的动态应力场；
10.4.3：矿房的动态回采过程是指矿房自第一个矿房回采开始至所有矿房回采结束的连续过程，矿房每动态回采一次，得到一组动态矿房回采的数值模拟结果，所有矿房回采完成后，得到采空区的应力场；
10.5：围岩稳定性的动态分析，过程如下：
10.5.1：分析区域为矿房顶板、上下盘围岩和矿柱；
10.5.2：岩体剪切破坏的判断准则选用以主应力形式表示的Mohr-Coulomb屈服准则，当fs＞0时，岩体发生剪切破坏，公式如下：
式中：σ1、σ3分别为最大主应力和最小主应力；c为粘聚力；
为摩擦角；
10.5.3：岩体张拉破坏的判断准则，根据岩体所受拉应力σt与抗拉强度Rt之间的关系确定，当σt≥Rt时，F≥0，即认为岩体发生张拉破坏，公式如下：
F＝σt-Rt
式中：σt为岩体所受拉应力；Rt为岩体抗拉强度；
10.5.4：根据分析结果，得到每次矿房动态回采后，矿房顶板、上下盘围岩和矿柱的破坏模式；
10.5.5：将每次矿房动态回采后得到的矿房顶板、上下盘围岩和矿柱的破坏模式汇总，得到围岩的稳定性动态分析结果；
(11)空区稳定性动态评价方法。
2.如权利要求1所述的基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述步骤(11)中，空区稳定性动态评价方法的过程如下：
11.1：以改进Mathews稳定图法得到的围岩稳定性评价结果作为采空区稳定性的初始评价结果；
11.2：利用数值模拟得到的围岩稳定性动态分析结果，对初始评价结果进行修正，过程如下：
11.2.1：当改进Mathews稳定图法评价的采空区围岩处于崩落区时，则该采空区处于崩落状态；
11.2.2：当改进Mathews稳定图法评价的采空区围岩处于破坏区时，则该采空区处于破坏状态；
11.2.3：当改进Mathews稳定图法评价的采空区围岩处于稳定区时，若在围岩稳定性动态分析中得到的采空区顶板、间柱、上下盘围岩未发生破坏，则该采空区处于稳定状态，若在围岩稳定性动态分析中得到的采空区顶板、间柱、上下盘围岩发生破坏，则该采空区有可能处于破坏状态；
11.3：得到采空区的稳定性动态评价结果。
3.如权利要求1所述的基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述步骤(9)中，改进Mathews稳定图评价方法的过程如下：
9.1：选取用于开展上下盘围岩稳定性评价的矿房；
9.2：测量出矿房上下盘围岩宽度X和长度Y，计算出岩石应力系数A；节理产状调整系数B和重力调整系数C；
9.3：测量并求解出矿房上下盘的节理组数Jn；节理粗糙度系数Jr；节理蚀变系数Ja；节理水折减系数Jw和SRF为应力折减系数；
9.4：求解出最佳阈值t下的RQDtθ、RQDtmin和
计算出相关拟合系数a和b，求解出考虑RQDt各向异性特征的Q′值；
9.5：计算出矿房稳定数N和水力半径R；
9.6：计算出稳定区-破坏边界的曲线公式和破坏-崩落边界的曲线公式，并绘制到改进后的Mathews稳定性图上；
9.7：将计算出的矿房稳定数N和水力半径R绘制到改进后的Mathews稳定图上，根据绘制的散点所处的区域，给出矿房上下盘围岩的稳定性评价结果，若散点处于稳定区，则矿房围岩评价结果为稳定，若散点处于破坏或主要破坏区，需再根据矿房围岩破坏和崩落的概率，给出矿房围岩的破坏概率和崩落概率，若散点处于崩落区，则矿房围岩评价结果为崩落。
4.如权利要求1所述的基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述步骤(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由完整岩石单轴抗压强度与采场中线采矿产生的地应力的比值确定，公式如下：
当
A＝0.1
当
当
A＝1
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最小值；
为RQDt的均值；a，b为相关拟合系数；Jn为节理组数；Jr为节理粗糙度系数；Ja为节理蚀变系数；Jw为节理水折减系数；SRF为应力折减系数；
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％的概率崩落；处在崩落区的工程，将持续发生崩落。
5.如权利要求1所述的基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述步骤(7)中，RQDt各向异性求解方法的过程如下：
7.1：基于最佳阈值t和剖切的三个二维裂隙网络模型，求解出最佳阈值t下的RQDt值，每个裂隙网络模型求解出36个RQDt值；
7.2：计算出每个裂隙网络模型上RQDt的最大值RQDtmax、最小值RQDtmin和均值
7.3：根据研究方向和RQDtmin方向位置关系，结合方差进行修正，提出各向异性条件下RQDt的计算公式如下：
式中：
为RQDt的均值，RQD′t为RQDt的修正系数；
7.4：修正系数RQD′t的求解考虑研究方向在某一具体方位角下的RQDtθ值和方差D，按如下方法修正，当RQDtθ＝RQDtmin时，RQD′t＝-D，D为方差；当RQDtθ＝RQDtmax时，RQD′t＝D；当
时，RQD′t＝0；
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和
以及方差的关系，求解出三组RQD′t值，并基于三组RQD′t值，绘制出散点图；
7.7.3：根据散点图做出拟合曲线，曲线截距即为a值，斜率即为b值，曲线方程即为RQD′t的修正公式；
7.8：将求解出的a、b值，带入到修正系数公式及RQDt各向异性计算公式中，得到与角度θ有关的RQDt各向异性公式；
7.9：根据RQDt各向异性公式，求解出任意角度的RQDt值。
6.如权利要求1所述的基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述步骤(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值。
7.如权利要求1所述的基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述步骤(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，a)；
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α表示，以及对应的测线起点与终点的距离la；
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的值。
8.如权利要求1所述的基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述步骤(4)中，岩体三维裂隙网络模型生成和剖切的过程如下：
4.1：随机数求解，过程如下：
4.1.1：产生随机数的数学方法应满足以下条件：产生的随机数列应均匀分布在(0，1)区间；序列之间应无相关性；随机序列有足够长的重复周期，在计算机上产生的速度快，占有的内存空间小，具有完全可重复性；
4.1.2：利用Monte Carlo方法是根据确立的结构面几何概率模型，再现出服从这种模型的结构面网络模型，在(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：输出三维裂隙网络模型的数据和任意二维裂隙网络模型的数据。
9.如权利要求1所述的基于激光扫描、BQ、数值模拟的空区稳定性动态评价方法，其特征在于，所述步骤(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：输出结构面产状的分组信息，包括每组结构面的倾向、倾角、迹长、间距和断距的均值和方差。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011011892.7A CN112200431A (zh) | 2020-09-16 | 2020-09-16 | 基于激光扫描、bq、数值模拟的空区稳定性动态评价方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011011892.7A CN112200431A (zh) | 2020-09-16 | 2020-09-16 | 基于激光扫描、bq、数值模拟的空区稳定性动态评价方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112200431A true CN112200431A (zh) | 2021-01-08 |
Family
ID=74014568
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011011892.7A Withdrawn CN112200431A (zh) | 2020-09-16 | 2020-09-16 | 基于激光扫描、bq、数值模拟的空区稳定性动态评价方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112200431A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113420067A (zh) * | 2021-06-22 | 2021-09-21 | 北京房江湖科技有限公司 | 目标地点的位置可信度评估方法和装置 |
-
2020
- 2020-09-16 CN CN202011011892.7A patent/CN112200431A/zh not_active Withdrawn
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113420067A (zh) * | 2021-06-22 | 2021-09-21 | 北京房江湖科技有限公司 | 目标地点的位置可信度评估方法和装置 |
| CN113420067B (zh) * | 2021-06-22 | 2024-01-19 | 贝壳找房(北京)科技有限公司 | 目标地点的位置可信度评估方法和装置 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112150001A (zh) | 基于摄影测量、BQ、改进Mathews稳定图的围岩稳定性评价方法 | |
| CN112200426A (zh) | 基于激光扫描、bq、数值模拟的围岩稳定性动态评价方法 | |
| CN112200419A (zh) | 基于激光扫描、BQ、改进Mathews稳定图的围岩稳定性评价方法 | |
| CN112200423A (zh) | 一种基于bq、数值模拟的围岩稳定性动态评价方法 | |
| CN112200417A (zh) | 基于摄影测量、BQ、RQDt、地应力的改进Mathews稳定图评价方法 | |
| CN112200431A (zh) | 基于激光扫描、bq、数值模拟的空区稳定性动态评价方法 | |
| CN112132407A (zh) | 一种基于BQ反演最佳阈值t的空间RQDt求解方法 | |
| CN112149996A (zh) | 基于激光扫描、BQ、RQDt各向异性的改进Mathews稳定图评价方法 | |
| CN112200418A (zh) | 基于摄影测量、bq、数值模拟的围岩稳定性动态评价方法 | |
| CN112150002A (zh) | 基于激光扫描、BQ、RQDt、地应力的改进Mathews稳定图评价方法 | |
| CN112132411A (zh) | 基于激光扫描、BQ、RQDt各向异性的Q各向异性求解方法 | |
| CN112200420A (zh) | 基于激光扫描、BQ、改进Mathews稳定图的空区稳定性评价方法 | |
| CN112150005A (zh) | 基于激光扫描、BQ、RQDt各向异性的Mathews稳定图法的改进方法 | |
| CN112132408A (zh) | 一种基于激光扫描和BQ反演最佳阈值t的空间RQDt求解方法 | |
| CN112200428A (zh) | 基于摄影测量、bq、数值模拟的空区稳定性动态评价方法 | |
| CN112464514A (zh) | 一种基于摄影测量、RQD和RQDt的巷道开挖不利方位求解方法 | |
| CN112464516A (zh) | 一种基于激光扫描和RQD反演最佳阈值t的空间RQDt求解方法 | |
| CN112200432A (zh) | 基于摄影测量、BQ、改进Mathews稳定图的空区稳定性评价方法 | |
| CN112200421A (zh) | 基于激光扫描、BQ、RQDt、地应力的Mathews稳定图法的改进方法 | |
| CN112200416A (zh) | 一种基于bq、数值模拟的空区稳定性动态评价方法 | |
| CN112150003A (zh) | 基于摄影测量、BQ、RQDt各向异性的改进Mathews稳定图评价方法 | |
| CN112132403A (zh) | 一种基于摄影测量和BQ反演的RQDt最佳阈值t求解方法 | |
| CN112132410A (zh) | 一种基于摄影测量和BQ反演最佳阈值t的空间RQDt求解方法 | |
| CN112200425A (zh) | 一种基于BQ、RQDt各向异性的改进Mathews稳定图评价方法 | |
| CN112150000A (zh) | 基于摄影测量、BQ、RQDt各向异性的Mathews稳定图法的改进方法 |
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 |
Application publication date: 20210108 |
|
| WW01 | Invention patent application withdrawn after publication |