CN112200426A - Dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation - Google Patents
Dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation Download PDFInfo
- Publication number
- CN112200426A CN112200426A CN202011011852.2A CN202011011852A CN112200426A CN 112200426 A CN112200426 A CN 112200426A CN 202011011852 A CN202011011852 A CN 202011011852A CN 112200426 A CN112200426 A CN 112200426A
- Authority
- CN
- China
- Prior art keywords
- rqd
- stability
- rock
- value
- follows
- 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
- 239000011435 rock Substances 0.000 title claims abstract description 243
- 238000011156 evaluation Methods 0.000 title claims abstract description 64
- 238000004088 simulation Methods 0.000 title claims abstract description 44
- 238000000034 method Methods 0.000 claims abstract description 192
- 238000010586 diagram Methods 0.000 claims abstract description 65
- 238000004458 analytical method Methods 0.000 claims abstract description 41
- 230000006641 stabilisation Effects 0.000 claims abstract description 3
- 238000011105 stabilization Methods 0.000 claims abstract description 3
- 230000008569 process Effects 0.000 claims description 79
- 238000004364 calculation method Methods 0.000 claims description 35
- 238000012937 correction Methods 0.000 claims description 31
- 239000013598 vector Substances 0.000 claims description 29
- 238000012360 testing method Methods 0.000 claims description 27
- 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
- 238000009933 burial Methods 0.000 claims description 15
- 238000013461 design Methods 0.000 claims description 12
- 230000009467 reduction Effects 0.000 claims description 12
- 238000005070 sampling Methods 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
- 238000005259 measurement Methods 0.000 claims description 8
- 238000004422 calculation algorithm Methods 0.000 claims description 7
- 238000005553 drilling Methods 0.000 claims description 7
- 230000006872 improvement Effects 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 6
- 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
- 230000000694 effects Effects 0.000 claims description 4
- 238000013507 mapping Methods 0.000 claims description 4
- 238000000691 measurement method 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
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 238000013028 emission testing Methods 0.000 claims description 3
- 238000002474 experimental method Methods 0.000 claims description 3
- 229910052500 inorganic mineral Inorganic materials 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 3
- 238000012067 mathematical method Methods 0.000 claims description 3
- 239000011707 mineral Substances 0.000 claims description 3
- 238000011084 recovery Methods 0.000 claims description 3
- 230000001105 regulatory effect Effects 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 9
- 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
- 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)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Strategic Management (AREA)
- Evolutionary Computation (AREA)
- Entrepreneurship & Innovation (AREA)
- Development Economics (AREA)
- Educational Administration (AREA)
- Economics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Computer Hardware Design (AREA)
- Life Sciences & Earth Sciences (AREA)
- Operations Research (AREA)
- Mathematical Analysis (AREA)
- Evolutionary Biology (AREA)
- Marketing (AREA)
- Bioinformatics & Computational Biology (AREA)
- Quality & Reliability (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Game Theory and Decision Science (AREA)
- Computational Mathematics (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Analytical Chemistry (AREA)
- Chemical & Material Sciences (AREA)
- Biochemistry (AREA)
- Software Systems (AREA)
Abstract
A dynamic evaluation method for surrounding rock stability based on laser scanning, BQ and numerical simulation belongs to the field of goaf stability evaluation, and comprises the following steps: (1) three-dimensional laser scanning of the structural surface is rapidly acquired; (2) structural surface clustering analysis; (3) calculating the rock mass based on the BQ index; (4) generating and sectioning a rock three-dimensional fracture network model; (5) RQDtDrawing an anisotropy map; (6) solving method of optimal threshold t based on BQ inversion; (7) RQDtAn anisotropy solution method; (8) a ground stress measuring method; (9) an improved method of Mathews graph stabilization; (10) improving a Mathews stability chart evaluation method; (11) a dynamic analysis method of the stability of the surrounding rock; (12) a dynamic evaluation method for stability of a dead zone. The invention realizes dynamic evaluation of the stability of the surrounding rock of the dead zone based on the combination of an improved Mathews stability diagram method and numerical simulation. The method is clear and suitable for dynamic evaluation of the stability of the surrounding rock of the goaf.
Description
Technical Field
The invention relates to a dynamic evaluation method for surrounding rock stability based on laser scanning, BQ and numerical simulation, in particular to an RQD (total Quadrature-resolution) inversion calculation method based on three-dimensional laser scanning, a neighbor propagation algorithm, a BQ index, a fracture network model and a generalized RQD theorytThe range of the optimal threshold t and the optimal threshold t value give the RQDtMethod for solving anisotropy and based on RQDtThe anisotropy and ground stress fitting method provides an improved method of a Mathews stability diagram method and an improved Mathews stability diagram evaluation method, and provides a dynamic evaluation method of the stability of the surrounding rock based on laser scanning, BQ and numerical simulation by combining a numerical simulation method, and belongs to the field of goaf stability evaluation.
Background
After long-time underground mining, a large-scale goaf group can be formed in the metal mine, the goaf group is mutually communicated, overlapped and concentrated underground, stress disturbance and mutual restriction are realized between the goaf group and the goaf group, a complex group space relationship is formed, and a stable rock mass state is in dynamic stable balance. Due to the particularity and complexity of a rock mass structure, randomness, uncertainty, continuity and scale in the goaf group stoping process, geological disasters caused by unstable collapse of the goaf group are serious, destructive power is strong, and personnel and property losses caused by the fact are huge.
The destruction of the goaf group is a dynamic destruction process. Firstly, due to the particularity of the rock mass structure of the empty area group, the damage of the goafs is often scaled, and after one goaf is destabilized and damaged, due to the change of the rock mass structure, the adjacent goafs can be affected to generate chain damage, so that the scale damage of the empty area group is further caused. Secondly, due to the continuity and the scale of the goaf group during stoping, a plurality of goafs are stoped at the same time, and a plurality of disturbances exist at the same time, so that the dynamic damage process of the rock mass structure of the goaf group is aggravated, and the risk of rock instability damage is increased. Due to the characteristics of complexity, variability, burstiness and scalability of the goaf group, people still lack sufficient cognitive understanding of the process of the damage occurrence and development, and many key scientific problems are not effectively solved, and an effective goaf stability dynamic evaluation method is not formed.
The numerical simulation method has remarkable advantages in the aspect of researching the stability and the stress evolution rule of the surrounding rock, and dynamic analysis of the stability of the surrounding rock can be realized through simulation research on geological models with any size and shape. The Mathews stability diagram method is one of goaf surrounding rock stability analysis methods, and stability evaluation results of the goaf can be obtained through stability evaluation of surrounding rocks of the upper and lower walls of the goaf. A dynamic evaluation method for the stability of the goaf can be provided and obtained by combining a Mathews stability diagram method and a numerical simulation method.
The Mathews-stabilized map method is a Q-system-based stabilized map method proposed in 1980 by Mathews et al by analyzing a large number of engineering examples. In the evaluation process, conditions such as stope characteristics, geological conditions, joint occurrence and the like are mainly considered by the Mathews stability diagram method, the rock mass stability number N and the hydraulic radius R are calculated, and the two factors are drawn on a diagram divided into a prediction stable area, a potential unstable area and a collapse area.
In terms of an improvement of the Mathews stabbing map method, Stewart and Forsyth verified and modified the stabbing map method in 1995 by collecting data at different mining depths, dividing it into four partitions. In 1998, Potvin improved the stable graph method based on the data example, and a joint orientation correction coefficient and a gravity adjustment coefficient are added in the graph. In 2000, Pakalnis et al summarized the critical span chart method based on RMR. Trueman et al redefined the plateau and severe disruption regions using a logistic regression approach. Mawdeskey in 2004 gave an equiprobable map of stable, damaged and severely damaged regions by logistic regression analysis of 405 examples.
However, a large number of structural surfaces and cracks exist in the rock mass, which significantly affect the mechanical properties, the crushing degree and the stability of the rock mass and are one of the important factors affecting the stability of the chamber. Structural planes and cracks in the rock mass are important factors causing the anisotropy of the mass of the rock mass, and the anisotropy characteristic of the mass of the rock mass is not considered in the Mathews stability diagram method.
The ground stress is the stress existing in the earth crust, is the acting force on the unit area inside the medium caused by rock deformation, the stress state of the ground stress at each underground point is different, the ground stress linearly increases along with the increase of the depth, the structure degree and the geographical position of each underground rock body are different, the gradient of the increase of the ground stress at each point is different, and the research on the relation of the ground stress changing along with the buried depth is relatively less in the Mathews stable diagram method.
The appearance and quality of the rock mass strongly depend on the degree of anisotropy of the rock mass, and the stability of the rock mass is determined by the quality of the rock mass. The RQD is an important index for representing the quality of rock mass, and has anisotropy. If the drilling mode is adopted to obtain the RQD, the results obtained from different rock body parts are completely different, the result of the RQD value depends on the direction, when the drilling direction is parallel to the main joint group, a higher RQD value is obtained, and when the drilling direction is perpendicular to the main joint group, a lower RQD value is obtained. The anisotropy of the RQD can also be seen from the expansion of the calculation formula from the borehole RQD to the generalized RQD. De er proposed the concept of a drill RQD in 1964, which has two disadvantages in its application: for rock masses of different engineering scales, whether the threshold value of 100mm is reasonable or not is selected; the drilling direction of the drill hole has limitation, and the obtained RQD can only reflect the local rock mass condition and cannot reflect the anisotropy of the RQD. Therefore, some scholars have introduced the threshold t and proposed the concept of generalized RQD, i.e. for any pitch threshold t, defining the percentage of the sum of pitches greater than t in a certain line direction to the total length of the line as the generalized RQD, using RQDtAnd (4) showing. The introduction of a generalized RQD makes it possible to solve the RQD anisotropy.
By RQDtThe concept can be known that the threshold t is a particularly important parameter of the generalized RQD and is a key for whether the generalized RQD can truly reflect the quality of the rock mass, but the threshold t has arbitrariness and no uniqueness, so that the optimal threshold t capable of representing the quality of the rock mass is found, and the method has very important scientific significance and research significance.
RQDtThe threshold value t is influenced by the measurement line direction, the structural plane form and the distribution characteristics, and the structural plane which is widely developed in the rock mass is an important factor for destroying the continuity and the integrity of the rock mass and controlling the mechanical characteristics and the stability of the rock mass, and plays a role in controlling the quality of the rock mass. The development mode and distribution form of the structural surface are very complex, but at the same time, certain generation relations exist among different joint groups, joints and faults, and a certain specific combination is formed to show certain regularity. Therefore, accurate and effective description analysis is carried out on the structural plane, and the occurrence and distribution characteristics of the structural plane are obtained, so that the method is the basis for researching the quality of the rock mass and the optimal threshold value t. The three-dimensional laser scanning technology is used as an efficient three-dimensional space information acquisition means, and has great advantages in the aspect of acquiring the azimuth and scale information of the structural surface. The neighbor propagation clustering algorithm has a good effect in structural plane occurrence clustering analysis, has the advantages of no influence of an initial clustering center and high calculation efficiency, and is widely applied to a plurality of fields. Therefore, based on the three-dimensional laser scanning technology and the neighbor propagation algorithm, the rapid identification acquisition and cluster analysis of the structural plane can be realized, and the structural plane is RQDtProvides a data base for the study of the optimal threshold t.
The study of scholars at home and abroad on the threshold t is mainly embodied in the following two aspects: RQD at different thresholds ttCalculation and study of the optimal threshold t. RQD at different thresholds ttIn terms of calculation, the existing research mainly discusses the change of the threshold t to the RQDtThe influence of the value is mainly realized by simulating a rock fracture network model and arranging virtual drilling holes to solve different RQDstThe vehicle studied is a fracture network model. However, in this study, it was mainly for the purpose of studying RQDtThe anisotropic characteristic of (2) is rarely involved in the research of the optimal threshold value t.
Some researchers have also studied the optimal threshold t in the research of the optimal threshold t. For example, some scholars respectively calculate the fractal dimension of the structural surface distribution by applying a fractal theory based on a three-dimensional structural surface network simulation technology, and by continuously changing the threshold of the generalized RQD,obtaining RQD under different thresholdstAnd comparing and analyzing the fractal dimension of each section with the generalized RQD value, and providing a basis for accurately selecting the optimal threshold of the generalized RQD. Some scholars establish 35 kinds of hypothetical three-dimensional fracture network models based on the modified block degree index MBi, measure generalized RQD values with different thresholds, and determine the optimal threshold of the generalized RQD. The research and the application conditions of the optimal threshold t in the two aspects have certain limitations, and the method is an optimal threshold t solving method under specific theories and backgrounds, and when the background or a model is changed, the optimal threshold t is not optimal any more. Moreover, since the fractal dimension or the block index does not have the function of characterizing the quality of the rock mass, whether the true quality of the rock mass can be reflected or not can be obtained through the optimal threshold t obtained through the fractal dimension or the block index is questionable.
Therefore, an RQD which can reflect the mass of the rock mass most is found and providedtThe method for solving the optimal threshold t has very important scientific and research significance and is also RQDtThe threshold t research needs to solve the problem.
The anisotropy of the RQD directly influences the quality of the rock mass, and the influence mechanism of the anisotropy of the RQD on the quality of the rock mass is not yet explored clearly. In terms of the threshold t, no one has yet given a solution calculation method for the optimal threshold t, and therefore no RQD based on the optimal threshold t has been obtained eithertThe anisotropy calculation formula does not obtain the RQD which can reflect the mass of the rock mass mosttAnd (3) an anisotropy solving method. To obtain RQDtThe anisotropy solving method is the basis and the premise for researching the mass anisotropy of the rock mass.
In view of the above, the invention provides a dynamic evaluation method for the stability of surrounding rock based on laser scanning, BQ and numerical simulation.
Disclosure of Invention
In order to realize dynamic evaluation of goaf stability, the invention provides a dynamic evaluation method of surrounding rock stability based on laser scanning, BQ and numerical simulation. The method is based on three-dimensional laser scanning, structural plane cluster analysis, BQ index, fracture network model and generalized RQD theory, and is used for invertingCalculate RQDtThe range of the optimal threshold t and the optimal threshold t value give the RQDtMethod for solving anisotropy and based on RQDtThe anisotropy and ground stress fitting method provides an improved method of a Mathews stability diagram method and an improved Mathews stability diagram evaluation method, and provides a dynamic evaluation method of the stability of the surrounding rock based on laser scanning, BQ and numerical simulation by combining a numerical simulation method.
In order to solve the technical problems, the invention provides the following technical scheme:
a dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation comprises the following steps:
(1) the three-dimensional laser scanning of the structural surface is rapidly obtained, and the process is as follows:
1.1: selecting a scanner site according to a scanning target and site conditions, erecting a tripod, ensuring that an instrument can completely acquire three-dimensional space point cloud information of a slope rock mass according to a certain scanning route during erection, ensuring that a tripod table top is horizontal as far as possible, and placing a control target;
1.2: placing a scanner host on the table surface of a tripod, fixing a knob, centering bubbles of the host through a coarse adjustment foot rest and a fine adjustment scanner base, and setting parameters of a scanner port;
1.3: starting scanning control software, configuring relevant parameters of a scanner, entering a scanner control interface, planning a scanning angle, setting a scanning range according to a scanning target, adjusting configuration parameters of a camera, and acquiring an image of the scanning target;
1.4: fixing a scanning range, acquiring a scanning interval, setting a sampling interval, starting data acquisition, checking scanning point cloud data and color information conditions in real time, and adjusting scanning parameter setting at any time according to scanning results;
1.5: exporting structural surface point cloud data;
(2) structural surface clustering analysis;
(3) the rock mass calculation based on the BQ index comprises the following steps:
3.1: and calculating the rock integrity coefficient according to the structural plane parameters, wherein the formula is as follows:
in the formula: jv is the volume regulating number of rock mass, unit bar/m3；
3.2: the Jv calculation is as follows:
in the formula: l1, L2.., Ln is the length of the measuring line perpendicular to the structure surface; n1, N2.., Nn is the number of structural surfaces in the same group;
3.3: calculating the BQ value according to the uniaxial compressive strength value and the integrity coefficient value of the rock mass:
BQ＝90+3Rc+250Kv
in the formula: rc is the uniaxial compressive strength of the rock; kv is the integrity coefficient of the rock mass;
3.4: in applying the BQ calculation formula, the following conditions are followed:
when Rc is greater than 90Kv +30, substituting Rc to 90Kv +30 and Kv to calculate a BQ value;
calculating the BQ value by substituting Kv of 0.04Rc +0.4 and Rc when Kv is more than 0.04Rc + 0.4;
3.5: correcting the BQ according to the occurrence of the underground water and the weak structural plane and the influence of natural stress, wherein the correction formula is as follows:
[BQ]＝BQ-100(K1+K2+K3)
in the formula: k1 is an underground water influence correction coefficient; k2 is a software structural plane attitude influence correction coefficient; k3 is a natural stress influence correction coefficient;
3.6: obtaining a corrected BQ rock mass grading result;
(4) generating and sectioning a rock three-dimensional fracture network model;
(5)RQDtdrawing an anisotropy map;
(6) solving method of optimal threshold t based on BQ inversion;
(7)RQDtan anisotropy solution method;
(8) the ground stress measuring method comprises the following steps:
8.1: the measurement of the ground stress is carried out by adopting an acoustic emission method in a direct measurement method;
8.2: in the vertical direction, sampling points are selected at certain depth distances;
8.3: drilling a hole in the original rock on site at a sampling point to extract a core sample, and processing the core sample into a standard cylindrical test piece, wherein the core taking process is as follows;
8.3.1: preparing test pieces in 6 different directions in the original rock at the position to be cored;
8.3.2: if the local coordinate system of the point is oxyz, 3 directions are set as coordinate axis directions, and the other 3 directions are selected as axial angle bisector directions in an oxy, oyz and ozx plane;
8.3.3: the sampling quantity of the test pieces in each direction is 20;
8.4: an indoor acoustic emission test is developed, and the process is as follows:
8.4.1: carrying out an experiment by adopting an acoustic emission testing system;
8.4.2: the test loading mode adopts displacement control loading, the loading rate is 0.01mm/min, the threshold value set by the acoustic emission system is 45dB, and data are automatically acquired and stored under the control of a microcomputer;
8.4.3: drawing a time-stress-acoustic emission counting relation curve of each test specimen according to the test data;
8.4.4: determining the Keyzing effect characteristic point of each test piece and the corresponding vertical principal stress of each test piece, and calculating the horizontal maximum principal stress, the minimum principal stress and the maximum principal stress direction;
8.5: the method for fitting the relation between the ground stress and the burial depth comprises the following steps:
8.5.1: drawing a scatter diagram of the vertical principal stress, the maximum principal stress and the minimum principal stress with the depth as an abscissa and the stress value as an ordinate;
8.5.2: selecting a proper fitting equation according to a scatter diagram rule, and fitting a relation curve of the vertical principal stress, the maximum principal stress and the minimum principal stress with the buried depth h;
8.5.3: obtaining values of vertical principal stress, maximum principal stress and minimum principal stress at any burial depth position according to the fitting relation curve;
(9) an improved method of Mathews graph stabilization;
(10) improving a Mathews stability chart evaluation method;
(11) the dynamic analysis method for the stability of the surrounding rock comprises the following steps:
11.1: establishing a three-dimensional model of a plurality of continuous rooms by using CAD software according to the plan view and the section view of the rooms;
11.2: importing the three-dimensional model into finite element analysis software, and endowing the three-dimensional model with material parameters including rock strength, density, elastic modulus, Poisson's ratio, cohesion, friction angle and volume weight;
11.3: endowing a three-dimensional model with boundary conditions and initial stress;
11.4: carrying out numerical simulation of dynamic stoping of the chamber, wherein the process is as follows:
11.4.1: carrying out numerical simulation of a non-stoping state of the chamber to obtain an initial stress field;
11.4.2: carrying out numerical simulation of dynamic stoping of the chamber to obtain a dynamic stress field of the chamber and the goaf in the dynamic stoping process of the chamber;
11.4.3: the dynamic stoping process of the chamber refers to a continuous process from the beginning of stoping of the first chamber to the end of stoping of all the chambers, a group of numerical simulation results of the stoping of the dynamic chamber are obtained every time the chamber is dynamically stoped, and a stress field of a goaf is obtained after the stoping of all the chambers is completed;
11.5: dynamic analysis of the stability of the surrounding rock, the process is as follows:
11.5.1: the analysis area comprises a roof of the chamber, surrounding rocks of an upper tray and a lower tray and a pillar;
11.5.2: the judgment criterion of rock mass shear failure is represented by Mohr-Coulomb yield criterion in a principal stress form, when fs is more than 0, the rock mass shear failure occurs, and the formula is as follows:
in the formula: sigma1、σ3Maximum principal stress and minimum principal stress, respectively; c is cohesion;
is a friction angle;
11.5.3: the criterion of the rock mass tension failure is determined according to the tensile stress sigma borne by the rock masstAnd tensile strength RtThe relationship between the two is determined when the sigma ist≥RtWhen F is larger than or equal to 0, the rock mass is considered to be subjected to tension failure, and the formula is as follows:
F＝σt-Rt
in the formula: sigmatThe tensile stress borne by the rock mass; rtThe tensile strength of the rock mass;
11.5.4: according to the analysis result, obtaining the destruction modes of the roof, the upper and lower wall surrounding rocks and the ore pillars of the chamber after the dynamic recovery of the chamber each time;
11.5.5: summarizing the failure modes of the roof, the upper and lower wall surrounding rocks and the ore pillars of the chamber obtained after the dynamic stoping of the chamber each time to obtain the dynamic stability analysis result of the surrounding rocks;
(12) a dynamic evaluation method for stability of a dead zone.
Further, in the step (12), the process of the dynamic evaluation method for the stability of the empty zone is as follows:
12.1: taking a surrounding rock stability evaluation result obtained by improving a Mathews stability diagram method as an initial evaluation result of the goaf stability;
12.2: and correcting the initial evaluation result by utilizing a dynamic analysis result of the stability of the surrounding rock obtained by numerical simulation, wherein the process is as follows:
12.2.1: when the surrounding rock of the goaf evaluated by the improved Mathews stability diagram method is in a caving zone, the goaf is in a caving state;
12.2.2: when the surrounding rock of the goaf evaluated by the improved Mathews stability diagram method is in a damaged area, the goaf is in a damaged state;
12.2.3: when the surrounding rock of the goaf evaluated by the Mathews stability diagram method is in a stable region, if the top plate, the stud, the surrounding rock of the upper and lower trays of the goaf obtained in the dynamic analysis of the stability of the surrounding rock are not damaged, the goaf is in a stable state, and if the top plate, the stud, the surrounding rock of the upper and lower trays of the goaf obtained in the dynamic analysis of the stability of the surrounding rock are damaged, the goaf is possibly in a damaged state;
12.3: and obtaining a dynamic evaluation result of the stability of the goaf.
Further, in the step (10), the process of improving the Mathews stability chart evaluation method is as follows:
10.1: selecting a chamber for evaluating the stability of the surrounding rock of the upper and lower trays;
10.2: measuring the width X and the length Y of surrounding rocks of an upper and a lower plate of the chamber, and calculating a joint attitude adjustment coefficient B and a gravity adjustment coefficient C;
10.3: calculating the ground stress value of the position of the chamber according to a ground stress fitting formula, and solving a rock stress coefficient A;
10.4: measuring and solving joint group number J of upper and lower plates of the chambern(ii) a Joint roughness coefficient Jr(ii) a Joint coefficient of change Ja(ii) a Joint water reduction coefficient JwAnd SRF is the stress reduction factor;
10.5: solving RQD under the optimal threshold ttθ、RQDtminAnd
calculating related fitting coefficients a and b, and solving to consider RQDtThe Q' value of the anisotropic feature;
10.6: calculating the stability number N and the hydraulic radius R of the chamber;
10.7: calculating a curve formula of a stable region-damage boundary and a curve formula of a damage-collapse boundary, and drawing the curve formulas on an improved Mathews stability graph;
10.8: and drawing the calculated stability number N and hydraulic radius R of the chamber on an improved Mathews stability diagram, giving out the stability evaluation result of the surrounding rocks of the upper and lower walls of the chamber according to the area where the drawn scattering point is located, if the scattering point is located in a stable area, the evaluation result of the surrounding rocks of the chamber is stable, if the scattering point is located in a damage or main damage area, the damage probability and the collapse probability of the surrounding rocks of the chamber are given out according to the damage and collapse probability of the surrounding rocks of the chamber, and if the scattering point is located in a collapse area, the evaluation result of the surrounding rocks of the chamber is collapsed.
Further, in the step (9), the process of the improved method of the Mathews stabilized mapping method is as follows:
9.1: the Mathews stable graph method comprises the following steps:
9.1.1: the design process of the Mathews stability diagram method is based on the calculation of two factors, namely a stability number N and a hydraulic radius R, and then the two factors are drawn on a diagram divided into a prediction stable area, a potential unstable area and a collapse area, wherein the stability number N represents the capability of a rock body to maintain stability under a given stress condition, the hydraulic radius R reflects the size and the shape of a goaf, and the Mathews stability diagram method has the following calculation formula:
N＝Q′ABC
in the formula: x, Y the width and length of the surrounding rock of the upper and lower walls of the chamber; q' is a result calculated according to the survey map or the borehole core record; a is the stress coefficient of the rock; b is a joint attitude adjustment coefficient; c is a gravity adjustment coefficient;
9.1.2: the rock stress coefficient A is determined by the ratio of the uniaxial compressive strength of the intact rock to the ground stress generated by stope centerline mining, and the formula is as follows:
when in use
When in use
When in use
9.1.3: the value of the joint attitude adjustment coefficient B is measured by the difference between the dip angle of the stope face and the dip angle of the main joint group;
9.1.4: the gravity adjustment coefficient C reflects the influence of the stope face attitude on the stope rock stability under the influence of gravity, the size of the gravity adjustment coefficient C depends on the caving, the sliding and the sliding of the side slope and the like of the exposed surface of the stope roof, and the relationship between the gravity adjustment coefficient C and the stope surface inclination angle alpha is determined by the following formula:
C＝8-7 cosα
9.2: the improvement method of Mathews stable graph method comprises the following steps:
9.2.1: when solving for the Q' value, use RQDtReplace RQD values recorded for the borehole core and apply the RQD valuestTaking into account the anisotropic character of (A), according to RQDtSolving the formula for anisotropy to obtain the RQDtThe Q' value of the anisotropic characteristic, the formula is as follows:
in the formula: RQDtθIs RQD at angle thetatA value; RQDtminIs the RQD minimum at threshold t;
is RQDtThe mean value of (a); a and b are related fitting coefficients; j. the design is a squarenThe number of joints is; j. the design is a squarerThe joint roughness coefficient; j. the design is a squareaIs the joint alteration coefficient; j. the design is a squarewTo save the water reduction coefficient; SRF is the stress reduction coefficient;
9.2.2: the ground stress value is obtained by measuring through an acoustic emission method, a relation curve of the ground stress and the burial depth h is obtained according to a fitting curve of measured data of measuring points at different depths, and the ground stress value of any burial depth position and any mineral house is obtained;
9.2.3: and (3) redrawing the Mathews stability graph by combining the equal probability contour graph, wherein the curve formula of the stability region-destruction boundary is as follows:
logN＝1.8206logR+1.618
the curve formula for the break-breakout boundary is as follows:
logN＝1.8076logR-3
9.2.4: the improved Mathews stability map is divided into three regions: a stable zone, a failure or major failure zone and a breakout zone, projects on a stable-failure boundary line, stopes with 57% probability stable, 43% probability failure, 0% probability breakout; engineering on the caving-damage boundary line, stope is stable with 0% probability, damaged with 5% probability and caving with 95% probability; the engineering in the stable region is stable; in the engineering of a damage or main damage area, the probability is stable at 57-0%, the probability is damaged at 43-5%, and the probability is collapsed at 0-95%; the project in the breakout area will continue to have breakout.
Further, in the step (7), RQDtThe process of the anisotropy solution method is as follows:
7.1: solving RQD under the optimal threshold t based on the optimal threshold t and three dissected two-dimensional fracture network modelstValue, each fracture network model solves 36 RQDstA value;
7.2: calculating RQD on each fracture network modeltMaximum value RQD oftmaxMinimum RQDtminSum mean value
7.3: according to study Direction and RQDtminCorrecting the direction position relation by combining the variance, and providing RQD under the anisotropic conditiontThe calculation formula of (a) is as follows:
in the formula:
is RQDtMean value of, RQD'tIs RQDtThe correction coefficient of (2);
7.4: correction coefficient RQD'tSolution of (a) takes into account the RQD of the study direction at a particular azimuthtθThe value and variance D, corrected as follows, when RQDtθ＝RQDtminThen, RQD't-D is variance; when RQD is reachedtθ＝RQDtmaxThen, RQD'tD; when in use
Then, RQD't＝0；
7.5: according to a correction coefficient RQD'tThe solution process of (2) provides RQD'tThe correction formula is as follows:
in the formula: RQDtθRQD at angle θtValue, RQDtminIs RQD below threshold ttMinimum value, a, b are correlation coefficients;
7.6: propose RQDtThe anisotropy calculation formula is as follows:
7.7: solving the correlation coefficients a and b, wherein the process is as follows:
7.7.1: calculating RQD on each fracture network modeltMaximum value RQD oftmaxMinimum RQDtminSum mean value
7.7.2: according to RQDtθWith RQDtmax、RQDtminAnd
and the variance are obtained to obtain three groups of RQD'tValues, and based on three groups RQD'tDrawing a scatter diagram;
7.7.3: a fitting curve is made according to the scatter diagram, the intercept of the curve is a value a, the slope is b, and the equation of the curve is RQD'tThe correction formula of (2);
7.8: the solved values of a and b are substituted into a correction coefficient formula and RQDtIn the anisotropy calculation formula, the RQD related to the angle theta is obtainedtAn anisotropy formula;
7.9: according to RQDtThe anisotropy formula is used for solving the RQD of any angletThe value is obtained.
Further, in the step (6), the process of the optimal threshold t solving method based on BQ inversion is as follows:
6.1: RQD inversion based on BQ indextThe range, process is as follows:
and (3) searching a rock quality index table by combining the rock quality grades calculated by BQ grading, and performing inversion to determine RQD (maximum-degree-of-saturation) at the rock quality gradetA range value;
6.2: the optimal threshold value t solving method comprises the following steps:
6.2.1: cutting three sections at any angle on the three-dimensional fracture network model through a central point O to obtain three two-dimensional fracture network models, and deriving the two-dimensional fracture network models and data;
6.2.2: setting different thresholds t aiming at each two-dimensional fracture network model, and solving RQD (maximum likelihood ratio) under different thresholds ttA value;
6.2.3: deriving RQD at different thresholds ttValue, calculating RQDtMean value;
6.2.4: the threshold value t is used as the abscissa, and RQD is usedtIs the ordinate, draws RQDtScatter plots that vary with threshold t;
6.2.5: according to the scatter diagram, setting a fitting equation and fitting RQDtA plot of variation with threshold t;
6.2.6: RQD to be invertedtRange values, brought into fitted RQDtIn the curve graph changing along with the threshold t, the function equation and the curve graph are combined to solve the RQDtObtaining three groups of threshold value t ranges in the range of the threshold value t;
6.2.7: taking the intersection of the ranges of the three groups of threshold values t as the range of the optimal threshold value t;
6.2.8: taking a midpoint value in the range of the optimal threshold t as the optimal threshold t to obtain the optimal threshold t;
6.2.9: and outputting the range of the optimal threshold value t and the value of the optimal threshold value t.
Further, in the step (5), RQDtThe process of anisotropy mapping is as follows:
5.1：RQDtsolving calculation, the process is as follows:
5.1.1：RQDtthe theoretical formula is as follows:
in the formula: x is the number ofiRepresenting the length of the whole rock or interval, RQD, of the ith greater than a given threshold t, in a certain line directiontRepresenting the rock quality index corresponding to the threshold t, i.e. RQD at the threshold ttA value;
5.1.2: determining a section center point O, a length a and a width b of the two-dimensional fracture network model, solving a center point O coordinate (X) by taking the lower left corner of the model as a coordinate origin, the horizontal right side as an X axis and the vertical upward side as a y axis0，Y0) And a boundary equation;
5.1.3: taking O as a starting point, drawing a measuring line at an angle of every alpha of 10 degrees, intersecting with the fracture network model, drawing 36 measuring lines in total, wherein the length L of the measuring line is equal to the distance from the O point to the boundary of the fracture network model, and is represented by L0-L35, and solving a measuring line equation;
5.1.4: judging the intersection point of the measuring line and the boundary, and setting the intersection point of the measuring line equation and the boundary equation as (x)a，ya) Sequentially connecting a measuring line equation and a boundary equation, and judging whether the measuring line is intersected with the boundary;
5.1.5: solving the intersection point (X) of the boundary equation of the survey line and the fracture network modelZ，α，YZ，α)；
5.1.6: determining the interval of the measuring line, wherein the principle is as follows:
5.1.7: from the coordinates (X) of the start of each joint in the fracture network modelb，Yb) And endpoint coordinate (X)c，Yc) Establishing a corresponding analytical equation;
5.1.8: solving the intersection points of the first measuring line and each joint equation, and circularly judging each intersection point (X)j，Yj) Range, if the intersection point is coincident with a < XjC and b < YjIf d, recording the intersection point, and traversing all the joint equations to obtain all m intersection points;
5.1.9: sequencing the recorded m intersection points and the starting point coordinates and the ending point coordinates of the measuring lines from small to large according to the abscissa or the ordinate;
5.1.10: calculating the distance d between adjacent intersectionsi；
5.1.11: inputting a threshold value t;
5.1.12: circular comparison diAnd t is the size of the initial l t0, the rule is as follows:
lt＝lt+diif d isi＞t
5.1.13: solving RQD corresponding to each measuring linetValue of mαIndicates, and the corresponding distance l between the start and end of the survey lineα；
5.1.14: circularly solving m corresponding to each measuring lineαObtaining RQD in 36 directionstA value;
5.2：RQDtthe anisotropy plot was plotted as follows:
5.2.1: mixing 36 RQDstValues, ordered in order of angle;
5.2.2: drawing a circle by taking the O point as the center of the circle and 1 as the radius, and finding the distance from the center of the circle by taking the ray angle as alpha
The points of (a) are marked;
5.2.3: the end points of 36 groups of rays are connected in sequence, if the RQD on one ray istIf the value is 0, the center of the circle is taken;
5.2.4: plotting RQDtAn anisotropy map of (a);
5.2.5: output RQDtAn anisotropy map;
5.3: RQD at different thresholds ttThe anisotropy plot was plotted as follows:
5.3.1: inputting different threshold values t and solving corresponding RQDtA value;
5.3.2: RQD at different thresholds ttValue, in degrees, in RQDtThe values are functions and are drawn under the same coordinate system;
5.3.3: obtaining RQD under different threshold values ttAn anisotropy map;
5.3.4: output RQD under different thresholds ttAn anisotropy map;
5.4: output RQDtThe value of (c).
Further, in the step (4), the generation and sectioning processes of the rock three-dimensional fracture network model are as follows:
4.1: solving the random number, wherein the process is as follows:
4.1.1: the mathematical method for generating random numbers should satisfy the following conditions: the generated random number sequence is uniformly distributed in the (0, 1) interval; there should be no correlation between sequences; the random sequence has a long enough repetition period, the generation speed on a computer is high, the occupied memory space is small, and the repeatability is complete;
4.1.2: according to the Monte Carlo method, a structural surface network model obeying the model is reproduced according to the established structural surface geometric probability model, uniformly distributed random variables are generated in a (0, 1) interval, and random numbers obeying other distributions are generated by using the uniformly distributed random variables;
4.1.3: the density function of the joint geometric parameters comprises normal distribution, log-normal distribution, negative exponential distribution and uniform distribution;
4.1.4: determining basic geometric parameters for generating joints according to the obtained random numbers;
4.2: generating a rock three-dimensional fracture network model, wherein the process is as follows:
4.2.1: according to the automatic statistical result of the structural surface data and the obtained random number, storing the data of each group of structural surfaces into a text file, and expressing the data by st.dat;
4.2.2: and the content format of dat data is as follows: the coordinates (x, y and z) of the center point of the disc of each structural plane, the radius D of the disc, the inclination angle DA, the inclination DD, the trend SD, the thickness thin, the normal directions NX, NY and NZ and joints are grouped;
4.2.3: in order to distinguish the structural surfaces of different groups, the structural surface discs of the same group are endowed with the same color and are represented by a number array GID;
4.2.4: writing a program by using Matlab software, reading a structural plane data file st.dat, and automatically generating a rock three-dimensional fracture network model after running;
4.2.5: obtaining a rock three-dimensional fracture network model;
4.3: cutting a two-dimensional fracture network model, wherein the process is as follows:
4.3.1: combining a Matlab software programming tool on the three-dimensional fracture network model, and taking the central point of the three-dimensional fracture network model as a center to realize the section cutting function at any angle;
4.3.2: obtaining a two-dimensional fracture network model of any angle passing through a central point;
4.3.3: combining a Matlab software programming tool on the three-dimensional fracture network model, and realizing the section cutting function of any angle and any direction on any position of the three-dimensional fracture network model;
4.3.4: obtaining a two-dimensional fracture network model at any angle and any direction;
4.3.5: and storing the data on the cut section into an st1.dat file, wherein the section is in a three-dimensional coordinate system, and the data formats in the file are as follows from left to right: joint center point coordinates (x, y, z), joint length D, inclination DA, inclination DD, strike SD, thickness thin, normal directions NX, NY, NZ;
4.3.6: converting the three-dimensional coordinate system into a two-dimensional coordinate system, and storing the two-dimensional profile data into a st2.dat file, wherein the data formats in the file are as follows from left to right: joint center point coordinates (x, y), joint length D, inclination DA, inclination DD, thickness thin, normal directions NX, NY, NZ;
4.4: and outputting the data of the three-dimensional fracture network model and the data of any two-dimensional fracture network model.
Further, in the step (2), the process of structural plane clustering analysis is as follows:
2.1: processing point cloud data, wherein the process is as follows:
2.1.1: importing structural surface point cloud data obtained by laser scanning;
2.1.2: calculating the distance between the current point and the adjacent point in the point cloud after the topological structure and the distance mean value, and identifying and eliminating noise points in the point cloud data through a distance threshold;
2.1.3: determining a spatial three-dimensional coordinate of the point cloud data according to the self spatial coordinate position of the three-dimensional laser scanner and the occurrence orientation of the field structure surface;
2.1.4: converting point cloud data based on a lower hemisphere equal angle projection method, and converting joint attitude data represented by inclination and dip angles into structural surface attitude data represented by joint unit normal vectors;
2.1.5: obtaining structural surface data expressed by unit normal vectors;
2.2: the process of cluster analysis by a neighbor propagation algorithm is as follows:
2.2.1: let the number of actually measured samples of the structural surface be N, and the tendency of each sample data be XiThe angle of inclination is YiWith a tendency X of each sample dataiAngle of inclination YiDetermining an initial clustering center as a cluster to obtain N initial clustering centers;
2.2.2: traversing all sample data according to a similarity measurement criterion, calculating the distance between each sample data and a clustering center, and distributing each sample data to the nearest clustering center to obtain N groups of data;
2.2.3: for each group of data, solving and calculating the clustering center of each group of data by a characteristic modulus analysis method, and firstly, calculating a matrix S according to the following formula under the assumption that l data exist in a certain group:
in the formula: (x)i，yi，zi) Is a unit normal vector of any structural surface, i belongs to (1, l);
2.2.4: next, the eigenvalues (τ) of the matrix S are solved1，τ2，τ3) And a feature vector ([ xi ])1，ξ2，ξ3) In which τ is1＜τ2＜τ3Feature vector xi corresponding to the maximum feature value3For the average vector of l vectors in the group, ξ3As a new cluster center;
2.2.5: for all sample data, repeatedly calculating the distance from each sample data to the clustering center, the matrix S, the characteristic value and the characteristic vector until the positions of all the clustering centers are fixed, and determining the grouping of the structural plane;
2.2.6: converting the structural surface attitude data expressed by unit normal vector into structural surface attitude data expressed by inclination and dip angle;
2.2.7: carrying out statistical analysis on structural surface attitude data, calculating an average value m and a standard deviation sigma of the inclination angle of the structural surface, and calculating a steady interval [ m-sigma, m + sigma ] of the inclination angle data;
2.2.8: judging tendency X of initial clustering center of sample dataiAnd the inclination angle YiWhether it falls within the robust interval [ m-sigma, m + sigma ]]If yes, finishing clustering analysis; if not, re-clustering the sample data until the tendency X of the initial clustering centeriAnd the inclination angle YiAll fall within the robust interval [ m-sigma, m + sigma]；
2.3: drawing a polar diagram, wherein the process is as follows;
2.3.1: based on normal attitude data of structural plane and according to spatial plano-projection diagram of structural planeSolving the coordinate point (x) of the declination projection of all the normal lines of the structural plane based on the longitudinal section principlen，yn)；
2.3.3: drawing a base circle with the diameter as the unit length, drawing two diameters of vertical and horizontal, and marking E, S, W, N;
2.3.4: the coordinate (x) of the declination plane of all the structural surfaces is measuredn，yn) Plotted on a base circle graph;
2.3.5: drawing a polar point diagram of the structural surface;
2.4: structural surface statistical analysis, the process is as follows:
2.4.1: determining a sample partition interval m;
2.4.2: solving sample range
2.4.3: calculate each partition interval Mm：
2.4.4: determining the probability of the sample falling in each partition interval, and counting the number N of the samples falling in each interval by using a computer circulation languagemCalculating the probability P of the number of samples by combining the total number N of samplesm；
2.4.5: solving sample mean
2.4.6: solving for sample variance S2Wherein S is the standard deviation;
2.4.7: according to the probability PmDrawing the probability distribution forms of the inclination, the dip angle, the trace length, the spacing and the fault distance of each group of structural surfaces;
2.5: and outputting grouping information of the occurrence of the structural planes, wherein the grouping information comprises the mean value and the variance of the inclination, the dip angle, the trace length, the spacing and the fault distance of each group of structural planes.
The invention has the following beneficial effects:
1. the method is based on three-dimensional laser scanning, structural plane cluster analysis, BQ indexes, fracture network models, generalized RQD theory and RQDtAnisotropy solving method and ground stressThe method comprises a force measurement method, an improved Mathews stability diagram method and a numerical simulation method, and provides a dynamic evaluation method for the stability of the surrounding rock based on laser scanning, BQ and numerical simulation;
2. the invention realizes the rapid three-dimensional laser scanning of the structural surface and the cluster analysis of the nearest neighbor propagation algorithm;
3. the invention realizes the BQ index classification of the rock mass;
4. the method realizes the generation of a three-dimensional fracture network model and the generation of a two-dimensional profile model;
5. the invention realizes RQDtSolving the calculation sum RQDtDrawing an anisotropy map;
6. the invention realizes RQDtSolving the optimal threshold t;
7. the invention realizes RQDtSolving anisotropy;
8. the invention realizes the measurement of the ground stress and the fitting of the relation between the ground stress and the buried depth;
9. the invention realizes the RQD-basedtThe improvement of the method of Mathews stability mapping of anisotropic features and ground stress fitting of (1);
10. the invention realizes the evaluation of the stability of the surrounding rock of the dead zone by improving the Mathews stability diagram method;
11. the method realizes dynamic analysis of the stability of the surrounding rock based on numerical simulation;
12. the method is clear, and the dynamic evaluation method for the goaf stability realizes the dynamic evaluation of the goaf stability.
Drawings
Fig. 1 is a structural plane pole view diagram.
FIG. 2 is a two-dimensional fracture network diagram.
FIG. 3 is RQD at different thresholds ttAn anisotropy map.
FIG. 4 is RQDtFitting a curve with threshold t.
FIG. 5 is RQD'tThe fitted curve of (1).
Fig. 6 is a fitted curve of maximum principal stress versus depth of burial.
Fig. 7 is a graph of the results of the stability evaluation of the surrounding rock of the chamber.
Detailed Description
The invention will be further explained with reference to the drawings.
Referring to fig. 1 to 7, a method for dynamically evaluating the stability of a surrounding rock based on laser scanning, BQ and numerical simulation includes the following steps:
1) the three-dimensional laser scanning of the structural surface is rapidly obtained, and the process is as follows:
1.1: selecting a scanner site according to a scanning target and site conditions, erecting a tripod, ensuring that an instrument can completely acquire three-dimensional space point cloud information of a slope rock mass according to a certain scanning route during erection, ensuring that a tripod table top is horizontal as far as possible, and placing a control target;
1.2: placing a scanner host on the table surface of a tripod, fixing a knob, centering bubbles of the host through a coarse adjustment foot rest and a fine adjustment scanner base, and setting parameters of a scanner port;
1.3: starting scanning control software, configuring relevant parameters of a scanner, entering a scanner control interface, planning a scanning angle, setting a scanning range according to a scanning target, adjusting configuration parameters of a camera, and acquiring an image of the scanning target;
1.4: fixing a scanning range, acquiring a scanning interval, setting a sampling interval, starting data acquisition, checking scanning point cloud data and color information conditions in real time, and adjusting scanning parameter setting at any time according to scanning results;
1.5: exporting structural surface point cloud data;
2) structural surface clustering analysis, the process is as follows:
2.1: processing point cloud data, wherein the process is as follows;
2.1.1: importing structural surface point cloud data;
2.1.2: calculating the distance between the current point and the adjacent point in the point cloud after the topological structure and the distance mean value, and identifying and eliminating noise points in the point cloud data through a distance threshold;
2.1.3: determining a spatial three-dimensional coordinate of the point cloud data according to the self spatial coordinate position of the three-dimensional laser scanner and the occurrence orientation of the field structure surface;
2.1.4: converting point cloud data based on a lower hemisphere equal-angle projection method;
2.1.5: will tend to be alphadAnd angle of inclination betadThe represented joint attitude data is converted into structural plane attitude data expressed by normal vector of joint unit, and alpha is setnAnd betanThe inclination direction and the inclination angle are respectively the unit normal vector of the structural plane, and the unit normal vector of any structural plane is expressed as X ═ X (X1，x2，x3) At this time, each point on the hemispherical surface corresponds to a joint occurrence form, and the formula is as follows:
X＝(x1，x2，x3) (1)
αd∈(0，360)，βd∈(0，90) (4)
2.1.6: obtaining structural surface data expressed by unit normal vectors;
2.2: the process of cluster analysis by a neighbor propagation algorithm is as follows:
2.2.1: let the number of actually measured samples of the structural surface be N, and the tendency of each sample data be XiThe angle of inclination is YiI ∈ (1, N); by inclination X of each sample dataiAngle of inclination YiDetermining an initial clustering center as a cluster to obtain N initial clustering centers;
2.2.2: traversing all sample data according to a similarity measurement criterion, calculating the distance between each sample data and a clustering center, and distributing each sample data to the nearest clustering center to obtain N groups of data;
2.2.3: for each group of data, solving and calculating the clustering center of each group of data by a characteristic modulus analysis method, and assuming that l data exist in a certain group, solving the clustering center according to the following method:
2.2.4: first, a matrix S is calculated as follows
In the formula: (x)i，yi，zi) Is a unit normal vector of any structural surface, i belongs to (1, l);
2.2.5: next, the eigenvalues (τ) of the matrix S are solved1，τ2，τ3) And a feature vector ([ xi ])1，ξ2，ξ3) In which τ is1＜τ2＜τ3Feature vector xi corresponding to the maximum feature value3For the average vector of l vectors in the group, ξ3As a new cluster center;
2.2.6: for all sample data, repeatedly calculating the distance from each sample data to the clustering center, the matrix S, the characteristic value and the characteristic vector until the positions of all the clustering centers are fixed, and determining the grouping of the structural plane;
2.2.7: converting the structural surface attitude data expressed by unit normal vector into structural surface attitude data expressed by inclination and dip angle;
2.2.8: carrying out statistical analysis on structural surface attitude data, calculating an average value m and a standard deviation sigma of the inclination angle of the structural surface, and calculating a steady interval [ m-sigma, m + sigma ] of the inclination angle data;
2.2.9: judging tendency X of initial clustering center of sample dataiAnd the inclination angle YiWhether it falls within the robust interval [ m-sigma, m + sigma ]]If yes, finishing clustering analysis; if not, re-clustering the sample data until the tendency X of the initial clustering centeriAnd the inclination angle YiAll fall within the robust interval [ m-sigma, m + sigma]；
2.3: drawing a polar diagram, wherein the process is as follows;
2.3.1: based on normal attitude data of the structural plane and based on the longitudinal section principle of the spatial plano projection diagram of the structural plane, the point A' is set as the planeThe red surface projection of the surface normal is combined with the red plane projection principle to calculate the coordinate x of A' on the red plane projection diagramnAnd ynThe formula is as follows:
2.3.2: solving the coordinate point (x) of the declination projection of all the structural surface normalsn，yn)；
2.3.3: drawing a base circle with the diameter as the unit length, drawing two diameters of vertical and horizontal, and marking E, S, W, N;
2.3.4: the coordinate (x) of the declination plane of all the structural surfaces is measuredn，yn) Plotted on a base circle graph;
2.3.5: drawing a polar point diagram of the structural plane is realized, as shown in FIG. 1;
2.4: structural surface statistical analysis, the process is as follows:
2.4.1: determining a sample partition interval m:
2.4.2: solving sample range
2.4.3: calculate each partition interval Mm：
2.4.4: determining the probability of the sample falling in each partition interval, and counting the number N of the samples falling in each interval by using a computer circulation languagemBinding to the sampleNumber N, calculating sample number probability Pm：
2.4.5: solving sample mean
2.4.6: solving for sample variance S2Wherein S is the standard deviation:
2.4.7: according to the probability PmDrawing the probability distribution forms of the inclination, the dip angle, the trace length, the spacing and the fault distance of each group of structural surfaces;
2.5: outputting grouping information of the occurrence of the structural surfaces, wherein the grouping information comprises the mean value and the variance of the inclination, the inclination angle, the trace length, the spacing and the fault distance of each group of structural surfaces;
3) the rock mass calculation based on the BQ index comprises the following steps:
3.1: and calculating the rock integrity coefficient according to the structural plane parameters, wherein the formula is as follows:
in the formula: jv is the volume regulating number of rock mass, unit bar/m3；
3.2: the Jv calculation is as follows:
in the formula: l1, L2.., Ln is the length of the measuring line perpendicular to the structure surface; n1, N2.., Nn is the number of structural surfaces in the same group;
3.3: calculating the BQ value according to the uniaxial compressive strength value and the integrity coefficient value of the rock mass:
BQ＝90+3Rc+250Kv (15)
in the formula: rc is the uniaxial compressive strength of the rock; kv is the integrity coefficient of the rock mass;
3.4: in applying the BQ calculation formula, the following conditions are followed:
when Rc is greater than 90Kv +30, substituting Rc to 90Kv +30 and Kv to calculate a BQ value;
calculating the BQ value by substituting Kv of 0.04Rc +0.4 and Rc when Kv is more than 0.04Rc + 0.4;
3.5: correcting the BQ according to the occurrence of the underground water and the weak structural plane and the influence of natural stress, wherein the correction formula is as follows:
[BQ]＝BQ-100(K1+K2+K3) (16)
in the formula: k1 is an underground water influence correction coefficient; k2 is a software structural plane attitude influence correction coefficient; k3 is a natural stress influence correction coefficient;
3.6: obtaining a corrected BQ rock mass grading result which is a grade III rock mass;
4) generating and sectioning a rock three-dimensional fracture network model, and carrying out the following process;
4.1: solving the random number, wherein the process is as follows;
4.1.1: the mathematical method for generating random numbers should satisfy the following conditions: the generated random number sequence is uniformly distributed in the (0, 1) interval; there should be no correlation between sequences; the random sequence has a long enough repetition period, the generation speed on a computer is high, the occupied memory space is small, and the repeatability is complete;
4.1.2: according to the Monte Carlo method, a structural surface network model obeying the model is reproduced according to the established structural surface geometric probability model, uniformly distributed random variables are generated in a (0, 1) interval, and random numbers obeying other distributions are generated by using the uniformly distributed random variables;
4.1.3: the density function of the joint geometric parameters comprises normal distribution, log-normal distribution, negative exponential distribution and uniform distribution;
4.1.4: determining basic geometric parameters for generating joints according to the obtained random numbers;
4.2: generating a rock three-dimensional fracture network model, and the process is as follows;
4.2.1: according to the automatic statistical result of the structural surface data and the obtained random number, storing the data of each group of structural surfaces into a text file, and expressing the data by st.dat;
4.2.2: and the content format of dat data is as follows: the coordinates (x, y and z) of the center point of the disc of each structural plane, the radius D of the disc, the inclination angle DA, the inclination DD, the trend SD, the thickness thin, the normal directions NX, NY and NZ and joints are grouped;
4.2.3: in order to distinguish the structural surfaces of different groups, the structural surface discs of the same group are endowed with the same color and are represented by a number array GID;
4.2.4: writing a program by using Matlab software, reading a structural plane data file st.dat, and automatically generating a rock three-dimensional fracture network model after running;
4.2.5: obtaining a rock three-dimensional fracture network model;
4.3: cutting a two-dimensional fracture network model, wherein the process is as follows;
4.3.1: combining a Matlab software programming tool on the three-dimensional fracture network model, and taking the central point of the three-dimensional fracture network model as a center to realize the section cutting function at any angle;
4.3.2: obtaining a two-dimensional fracture network model of any angle passing through a central point;
4.3.3: combining a Matlab software programming tool on the three-dimensional fracture network model, and realizing the section cutting function of any angle and any direction on any position of the three-dimensional fracture network model;
4.3.4: obtaining a two-dimensional fracture network model of any angle and any direction, as shown in FIG. 2;
4.3.5: and storing the data on the cut section into an st1.dat file, wherein the section is in a three-dimensional coordinate system, and the data formats in the file are as follows from left to right: joint center point coordinates (x, y, z), joint length D, inclination DA, inclination DD, strike SD, thickness thin, normal directions NX, NY, NZ;
4.3.6: converting the three-dimensional coordinate system into a two-dimensional coordinate system, and storing the two-dimensional profile data into a st2.dat file, wherein the data formats in the file are as follows from left to right: joint center point coordinates (x, y), joint length D, inclination DA, inclination DD, thickness thin, normal directions NX, NY, NZ;
4.4: outputting data of the three-dimensional fracture network model and data of any two-dimensional fracture network model;
5)RQDtthe anisotropy plot was plotted as follows:
5.1：RQDtsolving calculation, the process is as follows:
5.1.1：RQDtthe theoretical formula is as follows:
in the formula: x is the number ofiRepresenting the length of the whole rock or interval, RQD, of the ith greater than a given threshold t, in a certain line directiontRepresenting the rock quality index corresponding to the threshold t, i.e. RQD at the threshold ttA value;
5.1.2: determining a section center point O, a length a and a width b of the two-dimensional fracture network model, taking the lower left corner of the model as a coordinate origin, the horizontal right direction as an x axis and the vertical upward direction as a y axis, and then the coordinate of the center point O is as follows:
the boundary equation is:
5.1.3: taking O as a starting point, drawing a measuring line every other angle alpha of 10 degrees, intersecting with the fracture network model, drawing 36 measuring lines in total, wherein the length L of the measuring line is equal to the distance from the O point to the boundary of the fracture network model and is represented by L0-L35, and then the equation of the measuring line is as follows:
in the formula: s represents a survey line, and α represents an angle;
5.1.4: judging the intersection point of the measuring line and the boundary, and setting the intersection point of the measuring line equation and the boundary equation as (x)a，ya) And sequentially connecting a measuring line equation and a boundary equation, and judging whether the measuring line intersects with the boundary, wherein the principle is as follows:
5.1.5: solving the intersection point (X) of the boundary equation of the survey line and the fracture network modelZ，α，YZ，α)；
5.1.6: determining the interval of the measuring line, wherein the principle is as follows:
5.1.7: from the coordinates (X) of the start of each joint in the fracture network modelb，Yb) And endpoint coordinate (X)c，Yc) Establishing a corresponding analytical equation, and defining a joint equation as follows:
5.1.8: solving the intersection points of the first measuring line and each joint equation, and circularly judging each intersection point (X)j，Yj) Range, if the intersection point is coincident with a < XjC and b < YjIf d, recording the intersection point, and traversing all the joint equations to obtain all m intersection points;
5.1.9: sequencing the recorded m intersection points and the starting point coordinates and the ending point coordinates of the measuring lines from small to large according to the abscissa or the ordinate;
5.1.10: calculating the distance between adjacent intersection points, and the formula is as follows:
x0＝a/2 (26)
y0＝b/2 (27)
xm+1＝XZ，α (28)
ym+1＝YZ，α (29)
5.1.11: inputting a threshold value t;
5.1.12: circular comparison diAnd t is the size of the initial l t0, the rule is as follows:
lt＝lt+diif d isi＞t (30)
5.1.13: solving RQD corresponding to each measuring linetValue of mαIndicates, and the corresponding distance l between the start and end of the survey lineαThe formula is as follows:
5.1.14: circularly solving m corresponding to each measuring lineαObtaining RQD in 36 directionstA value;
5.2：RQDtthe anisotropy plot was plotted as follows:
5.2.1: mixing 36 RQDstValues, ordered in order of angle;
5.2.2: drawing a circle by taking the O point as the center of the circle and 1 as the radius, and finding the distance from the center of the circle by taking the ray angle as alpha
The points of (a) are marked;
5.2.3: the end points of 36 groups of rays are connected in sequence, if the RQD on one ray istIf the value is 0, the center of the circle is taken;
5.2.4: plotting RQDtAn anisotropy map of (a);
5.2.5: output RQDtAn anisotropy map;
5.3: RQD at different thresholds ttThe anisotropy plot was plotted as follows:
5.3.1: inputting different threshold values t and solving corresponding RQDtA value;
5.3.2: RQD at different thresholds ttValue, in degrees, in RQDtThe values are functions and are drawn under the same coordinate system;
5.3.3: obtaining RQD under different threshold values ttAn anisotropy map;
3.3.4: output RQD under different thresholds ttAnisotropy map, as shown in FIG. 3;
5.4: output RQDtA value of (d);
6) the optimal threshold t solving method based on BQ inversion comprises the following steps:
6.1: RQD inversion based on BQ indextThe range, process is as follows:
combining the rock mass quality grade calculated by BQ grading, wherein the rock mass is a grade III rock mass, searching a rock quality index table, and determining RQD under the rock mass grade through inversion according to the table 1tRange value, RQDtThe range is 50 to 75 percent
Table 1;
6.2: the optimal threshold value t solving method comprises the following steps:
6.2.1: cutting three sections at any angle on the three-dimensional fracture network model through a central point O to obtain three two-dimensional fracture network models, and deriving the two-dimensional fracture network models and data;
6.2.2: setting different thresholds t aiming at each two-dimensional fracture network model, and solving RQD (maximum likelihood ratio) under different thresholds ttA value;
6.2.3: deriving RQD at different thresholds ttValue, calculating RQDtMean value;
6.2.4: the threshold value t is used as the abscissa, and RQD is usedtIs the ordinate, draws RQDtScatter plots that vary with threshold t;
6.2.5: according to the scatter diagram, a fitting equation y ═ aexp (bx) is set, and an RQD is fittedtA graph of the variation with threshold t, as shown in fig. 4;
6.2.6: RQD to be invertedtRange values, brought into fitted RQDtIn the curve graph changing along with the threshold t, the function equation and the curve graph are combined to solve the RQDtWithin the range, threshold t ranges, resulting in three sets of threshold t ranges, see Table 2
Table 2;
6.2.7: aiming at the range of the three groups of threshold values t, taking the intersection of the ranges, and taking 0.076 m-0.124 m as the range of the optimal threshold value t;
6.2.8: taking a midpoint value in the range of the optimal threshold t as the optimal threshold t to obtain the optimal threshold t, wherein the optimal threshold t is 0.1 m;
6.2.9: outputting the range of the optimal threshold t and the value of the optimal threshold t;
7)RQDtthe anisotropy solving method comprises the following steps:
7.1: solving RQD under the optimal threshold t based on the optimal threshold t and three dissected two-dimensional fracture network modelstValue, each fracture network model solves 36 RQDstA value;
7.2: calculating RQD on each fracture network modeltMaximum value RQD oftmaxMinimum RQDtminSum mean value
7.3: according to study Direction and RQDtminCorrecting the direction position relation by combining the variance, and providing RQD under the anisotropic conditiontThe calculation formula of (a) is as follows:
in the formula:
is RQDtMean value of, RQD'tIs RQDtThe correction coefficient of (2);
7.4: correction coefficient RQD'tSolution of (a) takes into account the RQD of the study direction at a particular azimuthtθThe value and variance D, corrected as follows, when RQDtθ＝RQDtminThen, RQD't-D is variance; when RQD is reachedtθ＝RQDtmaxThen, RQD'tD; when in use
Then, RQD't＝0；
7.5: according to a correction coefficient RQD'tThe solution process of (2) provides RQD'tThe correction formula is as follows:
in the formula: RQDtθRQD at angle θtValue, RQDtminIs RQD below threshold ttMinimum value, a, b are correlation coefficients;
7.6: propose RQDtThe anisotropy calculation formula is as follows:
7.7: solving the correlation coefficients a and b, wherein the process is as follows:
7.7.1: calculating RQD on each fracture network modeltMaximum value RQD oftmaxMinimum RQDtminSum mean value
7.7.2: according to RQDtθWith RQDtmax、RQDtminAnd
and the variance are obtained to obtain three groups of RQD'tValues, and based on three groups RQD'tDrawing a scatter diagram;
7.7.3: a fitting curve is made according to the scatter diagram, as shown in FIG. 5, the intercept of the curve is a value a, the slope is b value, and the equation of the curve is RQD'tRQD 'obtained by the correction formula'tThe formula is as follows:
7.8: the solved values of a and b are substituted into a correction coefficient formula and RQDtIn the anisotropy calculation formula, the RQD related to the angle theta is obtainedtThe formula of anisotropy:
7.9: according to RQDtThe anisotropy formula is used for solving the RQD of any angletA value;
8) the ground stress measuring method comprises the following steps:
8.1: the measurement of the ground stress is carried out by adopting an acoustic emission method in a direct measurement method;
8.2: in the vertical direction, sampling points are selected at certain depth distances;
8.3: drilling a hole in the original rock on site at a sampling point to extract a core sample, and processing the core sample into a standard cylindrical test piece, wherein the core taking process is as follows;
8.3.1: preparing test pieces in 6 different directions in the original rock at the position to be cored;
8.3.2: if the local coordinate system of the point is oxyz, 3 directions are set as coordinate axis directions, and the other 3 directions are selected as axial angle bisector directions in an oxy, oyz and ozx plane;
8.3.3: the sampling quantity of the test pieces in each direction is 20;
8.4: an indoor acoustic emission test is developed, and the process is as follows:
8.4.1: carrying out an experiment by adopting an acoustic emission testing system;
8.4.2: the test loading mode adopts displacement control loading, the loading rate is 0.01mm/min, the threshold value set by the acoustic emission system is 45dB, and data are automatically acquired and stored under the control of a microcomputer;
8.4.3: drawing a time-stress-acoustic emission counting relation curve of each test specimen according to the test data;
8.4.4: determining the Keyzing effect characteristic point of each test piece and the corresponding vertical principal stress of each test piece, and calculating the horizontal maximum principal stress, the minimum principal stress and the maximum principal stress direction;
8.5: the method for fitting the relation between the ground stress and the burial depth comprises the following steps:
8.5.1: drawing a scatter diagram of the vertical principal stress, the maximum principal stress and the minimum principal stress with the depth as an abscissa and the stress value as an ordinate;
8.5.2: according to the scatter diagram rule, selecting a proper fitting equation, fitting a relation curve of the vertical principal stress, the maximum principal stress and the minimum principal stress and the burial depth h, wherein the fitting curve of the maximum principal stress and the burial depth is shown in FIG. 6;
8.5.3: obtaining values of vertical principal stress, maximum principal stress and minimum principal stress at any burial depth position according to the fitting relation curve;
9) the improvement method of Mathews stable graph method comprises the following steps:
9.1: the Mathews stable graph method comprises the following steps:
9.1.1: the design process of the Mathews stability diagram method is based on the calculation of two factors, namely a stability number N and a hydraulic radius R, and then the two factors are drawn on a diagram divided into a prediction stable area, a potential unstable area and a collapse area, wherein the stability number N represents the capability of a rock body to maintain stability under a given stress condition, the hydraulic radius R reflects the size and the shape of a goaf, and the Mathews stability diagram method has the following calculation formula:
N＝Q′ABC (38)
in the formula: x, Y the width and length of the surrounding rock of the upper and lower walls of the chamber; q' is a result calculated according to the survey map or the borehole core record; a is the stress coefficient of the rock; b is a joint attitude adjustment coefficient; c is a gravity adjustment coefficient;
9.1.2: the rock stress coefficient A is determined by the ratio of the uniaxial compressive strength of the intact rock to the ground stress generated by stope centerline mining, and the formula is as follows:
when in use
When in use
When in use
9.1.3: the value of the joint attitude adjustment coefficient B is measured by the difference between the dip angle of the stope face and the dip angle of the main joint group;
9.1.4: the gravity adjustment coefficient C reflects the influence of the stope face attitude on the stope rock stability under the influence of gravity, the size of the gravity adjustment coefficient C depends on the caving, the sliding and the sliding of the side slope and the like of the exposed surface of the stope roof, and the relationship between the gravity adjustment coefficient C and the stope surface inclination angle alpha is determined by the following formula:
C＝8-7 cosα (43)
9.2: the improvement method of Mathews stable graph method comprises the following steps:
9.2.1: when solving for the Q' value, use RQDtReplace RQD values recorded for the borehole core and apply the RQD valuestTaking into account the anisotropic character of (A), according to RQDtSolving the formula for anisotropy to obtain the RQDtThe Q' value of the anisotropic characteristic, the formula is as follows:
in the formula: RQDtθIs RQD at angle thetatA value; RQDtminIs the RQD minimum at threshold t;
is RQDtThe mean value of (a); a and b are related fitting coefficients; j. the design is a squarenThe number of joints is; j. the design is a squarerThe joint roughness coefficient; j. the design is a squareaIs the joint alteration coefficient; j. the design is a squarewTo save the water reduction coefficient; SRF is the stress reduction coefficient;
9.2.2: the ground stress value is obtained by measuring through an acoustic emission method, a relation curve of the ground stress and the burial depth h is obtained according to a fitting curve of measured data of measuring points at different depths, and the ground stress value of any burial depth position and any mineral house is obtained;
9.2.3: and (3) redrawing the Mathews stability graph by combining the equal probability contour graph, wherein the curve formula of the stability region-destruction boundary is as follows:
logN＝1.8206logR+1.618 (45)
the curve formula for the break-breakout boundary is as follows:
logN＝1.8076logR-3 (46)
9.2.4: the improved Mathews stability map is divided into three regions: a stable zone, a failure or major failure zone and a breakout zone, projects on a stable-failure boundary line, stopes with 57% probability stable, 43% probability failure, 0% probability breakout; engineering on the caving-damage boundary line, stope is stable with 0% probability, damaged with 5% probability and caving with 95% probability; the engineering in the stable region is stable; in the engineering of a damage or main damage area, the probability is stable at 57-0%, the probability is damaged at 43-5%, and the probability is collapsed at 0-95%; in the engineering in the collapse area, the collapse can continuously occur;
10) the method for evaluating the Mathews stability diagram is improved and comprises the following steps:
10.1: selecting a chamber for evaluating the stability of the surrounding rock of the upper and lower trays;
10.2: measuring the width X and the length Y of surrounding rocks of an upper and a lower plate of the chamber, and calculating a joint attitude adjustment coefficient B and a gravity adjustment coefficient C;
10.3: calculating the ground stress value of the position of the chamber according to a ground stress fitting formula, and solving a rock stress coefficient A;
10.4: measuring and solving joint group number J of upper and lower plates of the chambern(ii) a Joint roughness coefficient Jr(ii) a Joint coefficient of change Ja(ii) a Joint water reduction coefficient JwAnd SRF is the stress reduction factor;
10.5: solving RQD under the optimal threshold ttθ、RQDtminAnd
calculating related fitting coefficients a and b, and solving to consider RQDtThe Q' value of the anisotropic feature;
10.6: calculating the stability number N and the hydraulic radius R of the chamber;
10.7: calculating a curve formula of a stable region-damage boundary and a curve formula of a damage-collapse boundary, and drawing the curve formulas on an improved Mathews stability graph;
10.8: drawing the calculated stability number N and hydraulic radius R of the chamber on an improved Mathews stability diagram, as shown in FIG. 7, according to the area where the drawn scattering points are located, giving out the stability evaluation results of the surrounding rocks of the upper and lower walls of the chamber, wherein the 3203 chamber and 3204 chamber are in a caving area, the evaluation results of the surrounding rocks of the chamber are caving, the 3205 chamber is in a damaged or main damaged area, the probability of the damage of the surrounding rocks is 55%, the probability of the collapse of the surrounding rocks is 33%, and the probability of the stability of the surrounding rocks is 12%;
11) the dynamic analysis method for the stability of the surrounding rock comprises the following steps:
11.1: establishing a three-dimensional model of a plurality of continuous rooms by using CAD software according to the plan view and the section view of the rooms;
11.2: importing the three-dimensional model into finite element analysis software, and endowing the three-dimensional model with material parameters including rock strength, density, elastic modulus, Poisson's ratio, cohesion, friction angle and volume weight;
11.3: endowing a three-dimensional model with boundary conditions and initial stress;
11.4: carrying out numerical simulation of dynamic stoping of the chamber, wherein the process is as follows:
11.4.1: carrying out numerical simulation of a non-stoping state of the chamber to obtain an initial stress field;
11.4.2: carrying out numerical simulation of dynamic stoping of the chamber to obtain a dynamic stress field of the chamber and the goaf in the dynamic stoping process of the chamber;
11.4.3: the dynamic stoping process of the chamber refers to a continuous process from the beginning of stoping of the first chamber to the end of stoping of all the chambers, a group of numerical simulation results of the stoping of the dynamic chamber are obtained every time the chamber is dynamically stoped, and a stress field of a goaf is obtained after the stoping of all the chambers is completed;
11.5: dynamic analysis of the stability of the surrounding rock, the process is as follows:
11.5.1: the analysis area comprises a roof of the chamber, surrounding rocks of an upper tray and a lower tray and a pillar;
11.5.2: the judgment criterion of rock mass shear failure is represented by Mohr-Coulomb yield criterion expressed in the form of principal stress, when f issWhen the pressure is higher than 0, the rock mass is subjected to shear failure, and the formula is as follows:
in the formula: sigma1、σ3Maximum principal stress and minimum principal stress, respectively; c is cohesion;
is a friction angle;
11.5.3: the criterion of the rock mass tension failure is determined according to the tensile stress sigma borne by the rock masstAnd tensile strength RtThe relationship between the two is determined when the sigma ist≥RtWhen F is larger than or equal to 0, the rock mass is considered to be subjected to tension failure, and the formula is as follows:
F＝σt-Rt (48)
in the formula: sigmatThe tensile stress borne by the rock mass; rtThe tensile strength of the rock mass;
11.5.4: according to the analysis result, obtaining the destruction modes of the roof, the upper and lower wall surrounding rocks and the ore pillars of the chamber after the dynamic recovery of the chamber each time;
11.5.5: summarizing the failure modes of the roof, the upper and lower wall surrounding rocks and the ore pillars of the chamber obtained after the dynamic stoping of the chamber each time to obtain the dynamic stability analysis result of the surrounding rocks;
12) the dynamic evaluation method of the stability of the dead zone comprises the following steps:
12.1: taking a surrounding rock stability evaluation result obtained by improving a Mathews stability diagram method as an initial evaluation result of the goaf stability;
12.2: and correcting the initial evaluation result by utilizing a dynamic analysis result of the stability of the surrounding rock obtained by numerical simulation, wherein the process is as follows:
12.2.1: when the surrounding rock of the goaf evaluated by the improved Mathews stability diagram method is in a caving zone, the goaf is in a caving state;
12.2.2: when the surrounding rock of the goaf evaluated by the improved Mathews stability diagram method is in a damaged area, the goaf is in a damaged state;
12.2.3: when the surrounding rock of the goaf evaluated by the Mathews stability diagram method is in a stable region, if the top plate, the stud, the surrounding rock of the upper and lower trays of the goaf obtained in the dynamic analysis of the stability of the surrounding rock are not damaged, the goaf is in a stable state, and if the top plate, the stud, the surrounding rock of the upper and lower trays of the goaf obtained in the dynamic analysis of the stability of the surrounding rock are damaged, the goaf is possibly in a damaged state;
12.3: and obtaining a dynamic evaluation result of the stability of the goaf.
Claims (9)
1. A dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation is characterized by comprising the following steps:
(1) the three-dimensional laser scanning of the structural surface is rapidly obtained, and the process is as follows:
1.1: selecting a scanner site according to a scanning target and site conditions, erecting a tripod, ensuring that an instrument can completely acquire three-dimensional space point cloud information of a slope rock mass according to a certain scanning route during erection, ensuring that a tripod table top is horizontal as far as possible, and placing a control target;
1.2: placing a scanner host on the table surface of a tripod, fixing a knob, centering bubbles of the host through a coarse adjustment foot rest and a fine adjustment scanner base, and setting parameters of a scanner port;
1.3: starting scanning control software, configuring relevant parameters of a scanner, entering a scanner control interface, planning a scanning angle, setting a scanning range according to a scanning target, adjusting configuration parameters of a camera, and acquiring an image of the scanning target;
1.4: fixing a scanning range, acquiring a scanning interval, setting a sampling interval, starting data acquisition, checking scanning point cloud data and color information conditions in real time, and adjusting scanning parameter setting at any time according to scanning results;
1.5: exporting structural surface point cloud data;
(2) structural surface clustering analysis;
(3) the rock mass calculation based on the BQ index comprises the following steps:
3.1: and calculating the rock integrity coefficient according to the structural plane parameters, wherein the formula is as follows:
in the formula: jv is the volume regulating number of rock mass, unit bar/m3；
3.2: the Jv calculation is as follows:
in the formula: l1, L2.., Ln is the length of the measuring line perpendicular to the structure surface; n1, N2.., Nn is the number of structural surfaces in the same group;
3.3: calculating the BQ value according to the uniaxial compressive strength value and the integrity coefficient value of the rock mass:
BQ＝90+3Rc+250Kv
in the formula: rc is the uniaxial compressive strength of the rock; kv is the integrity coefficient of the rock mass;
3.4: in applying the BQ calculation formula, the following conditions are followed:
when Rc is greater than 90Kv +30, substituting Rc to 90Kv +30 and Kv to calculate a BQ value;
calculating the BQ value by substituting Kv of 0.04Rc +0.4 and Rc when Kv is more than 0.04Rc + 0.4;
3.5: correcting the BQ according to the occurrence of the underground water and the weak structural plane and the influence of natural stress, wherein the correction formula is as follows:
[BQ]＝BQ-100(K1+K2+K3)
in the formula: k1 is an underground water influence correction coefficient; k2 is a software structural plane attitude influence correction coefficient; k3 is a natural stress influence correction coefficient;
3.6: obtaining a corrected BQ rock mass grading result;
(4) generating and sectioning a rock three-dimensional fracture network model;
(5)RQDtdrawing an anisotropy map;
(6) solving method of optimal threshold t based on BQ inversion;
(7)RQDtan anisotropy solution method;
(8) the ground stress measuring method comprises the following steps:
8.1: the measurement of the ground stress is carried out by adopting an acoustic emission method in a direct measurement method;
8.2: in the vertical direction, sampling points are selected at certain depth distances;
8.3: drilling a hole in the original rock on site at a sampling point to extract a core sample, and processing the core sample into a standard cylindrical test piece, wherein the core taking process is as follows;
8.3.1: preparing test pieces in 6 different directions in the original rock at the position to be cored;
8.3.2: if the local coordinate system of the point is oxyz, 3 directions are set as coordinate axis directions, and the other 3 directions are selected as axial angle bisector directions in an oxy, oyz and ozx plane;
8.3.3: the sampling quantity of the test pieces in each direction is 20;
8.4: an indoor acoustic emission test is developed, and the process is as follows:
8.4.1: carrying out an experiment by adopting an acoustic emission testing system;
8.4.2: the test loading mode adopts displacement control loading, the loading rate is 0.01mm/min, the threshold value set by the acoustic emission system is 45dB, and data are automatically acquired and stored under the control of a microcomputer;
8.4.3: drawing a time-stress-acoustic emission counting relation curve of each test specimen according to the test data;
8.4.4: determining the Keyzing effect characteristic point of each test piece and the corresponding vertical principal stress of each test piece, and calculating the horizontal maximum principal stress, the minimum principal stress and the maximum principal stress direction;
8.5: the method for fitting the relation between the ground stress and the burial depth comprises the following steps:
8.5.1: drawing a scatter diagram of the vertical principal stress, the maximum principal stress and the minimum principal stress with the depth as an abscissa and the stress value as an ordinate;
8.5.2: selecting a proper fitting equation according to a scatter diagram rule, and fitting a relation curve of the vertical principal stress, the maximum principal stress and the minimum principal stress with the buried depth h;
8.5.3: obtaining values of vertical principal stress, maximum principal stress and minimum principal stress at any burial depth position according to the fitting relation curve;
(9) an improved method of Mathews graph stabilization;
(10) improving a Mathews stability chart evaluation method;
(11) the dynamic analysis method for the stability of the surrounding rock comprises the following steps:
11.1: establishing a three-dimensional model of a plurality of continuous rooms by using CAD software according to the plan view and the section view of the rooms;
11.2: importing the three-dimensional model into finite element analysis software, and endowing the three-dimensional model with material parameters including rock strength, density, elastic modulus, Poisson's ratio, cohesion, friction angle and volume weight;
11.3: endowing a three-dimensional model with boundary conditions and initial stress;
11.4: carrying out numerical simulation of dynamic stoping of the chamber, wherein the process is as follows:
11.4.1: carrying out numerical simulation of a non-stoping state of the chamber to obtain an initial stress field;
11.4.2: carrying out numerical simulation of dynamic stoping of the chamber to obtain a dynamic stress field of the chamber and the goaf in the dynamic stoping process of the chamber;
11.4.3: the dynamic stoping process of the chamber refers to a continuous process from the beginning of stoping of the first chamber to the end of stoping of all the chambers, a group of numerical simulation results of the stoping of the dynamic chamber are obtained every time the chamber is dynamically stoped, and a stress field of a goaf is obtained after the stoping of all the chambers is completed;
11.5: dynamic analysis of the stability of the surrounding rock, the process is as follows:
11.5.1: the analysis area comprises a roof of the chamber, surrounding rocks of an upper tray and a lower tray and a pillar;
11.5.2: the judgment criterion of rock mass shear failure is represented by Mohr-Coulomb yield criterion expressed in the form of principal stress, when f issWhen the pressure is higher than 0, the rock mass is subjected to shear failure, and the formula is as follows:
in the formula: sigma1、σ3Maximum principal stress and minimum principal stress, respectively; c is cohesion;
is a friction angle;
11.5.3: the criterion of the rock mass tension failure is determined according to the tensile stress sigma borne by the rock masstAnd tensile strength RtThe relationship between the two is determined when the sigma ist≥RtWhen F is larger than or equal to 0, the rock mass is considered to be subjected to tension failure, and the formula is as follows:
F＝σt-Rt
in the formula: sigmatThe tensile stress borne by the rock mass; rtThe tensile strength of the rock mass;
11.5.4: according to the analysis result, obtaining the destruction modes of the roof, the upper and lower wall surrounding rocks and the ore pillars of the chamber after the dynamic recovery of the chamber each time;
11.5.5: summarizing the failure modes of the roof, the upper and lower wall surrounding rocks and the ore pillars of the chamber obtained after the dynamic stoping of the chamber each time to obtain the dynamic stability analysis result of the surrounding rocks;
(12) a dynamic evaluation method for stability of a dead zone.
2. The laser scanning, BQ and numerical simulation based surrounding rock stability dynamic evaluation method as claimed in claim 1, wherein in the step (12), the process of the dead zone stability dynamic evaluation method is as follows:
12.1: taking a surrounding rock stability evaluation result obtained by improving a Mathews stability diagram method as an initial evaluation result of the goaf stability;
12.2: and correcting the initial evaluation result by utilizing a dynamic analysis result of the stability of the surrounding rock obtained by numerical simulation, wherein the process is as follows:
12.2.1: when the surrounding rock of the goaf evaluated by the improved Mathews stability diagram method is in a caving zone, the goaf is in a caving state;
12.2.2: when the surrounding rock of the goaf evaluated by the improved Mathews stability diagram method is in a damaged area, the goaf is in a damaged state;
12.2.3: when the surrounding rock of the goaf evaluated by the Mathews stability diagram method is in a stable region, if the top plate, the stud, the surrounding rock of the upper and lower trays of the goaf obtained in the dynamic analysis of the stability of the surrounding rock are not damaged, the goaf is in a stable state, and if the top plate, the stud, the surrounding rock of the upper and lower trays of the goaf obtained in the dynamic analysis of the stability of the surrounding rock are damaged, the goaf is possibly in a damaged state;
12.3: and obtaining a dynamic evaluation result of the stability of the goaf.
3. The method for dynamically evaluating the stability of the surrounding rock based on laser scanning, BQ and numerical simulation as claimed in claim 1, wherein in the step (10), the process of improving the Mathews stability chart evaluation method is as follows:
10.1: selecting a chamber for evaluating the stability of the surrounding rock of the upper and lower trays;
10.2: measuring the width X and the length Y of surrounding rocks of an upper and a lower plate of the chamber, and calculating a joint attitude adjustment coefficient B and a gravity adjustment coefficient C;
10.3: calculating the ground stress value of the position of the chamber according to a ground stress fitting formula, and solving a rock stress coefficient A;
10.4: measuring and solving joint group number J of upper and lower plates of the chambern(ii) a Joint roughness coefficient Jr(ii) a Joint coefficient of change Ja(ii) a Joint water reduction coefficient JwAnd SRF is the stress reduction factor;
10.5: solving RQD under the optimal threshold ttθ、RQDtminAnd
calculating related fitting coefficients a and b, and solving to consider RQDtThe Q' value of the anisotropic feature;
10.6: calculating the stability number N and the hydraulic radius R of the chamber;
10.7: calculating a curve formula of a stable region-damage boundary and a curve formula of a damage-collapse boundary, and drawing the curve formulas on an improved Mathews stability graph;
10.8: and drawing the calculated stability number N and hydraulic radius R of the chamber on an improved Mathews stability diagram, giving out the stability evaluation result of the surrounding rocks of the upper and lower walls of the chamber according to the area where the drawn scattering point is located, if the scattering point is located in a stable area, the evaluation result of the surrounding rocks of the chamber is stable, if the scattering point is located in a damage or main damage area, the damage probability and the collapse probability of the surrounding rocks of the chamber are given out according to the damage and collapse probability of the surrounding rocks of the chamber, and if the scattering point is located in a collapse area, the evaluation result of the surrounding rocks of the chamber is collapsed.
4. The dynamic evaluation method for the stability of the surrounding rock based on laser scanning, BQ and numerical simulation as claimed in claim 1, wherein in the step (9), the process of the improved Mathews stability chart method is as follows:
9.1: the Mathews stable graph method comprises the following steps:
9.1.1: the design process of the Mathews stability diagram method is based on the calculation of two factors, namely a stability number N and a hydraulic radius R, and then the two factors are drawn on a diagram divided into a prediction stable area, a potential unstable area and a collapse area, wherein the stability number N represents the capability of a rock body to maintain stability under a given stress condition, the hydraulic radius R reflects the size and the shape of a goaf, and the Mathews stability diagram method has the following calculation formula:
N＝Q′ABC
in the formula: x, Y the width and length of the surrounding rock of the upper and lower walls of the chamber; q' is a result calculated according to the survey map or the borehole core record; a is the stress coefficient of the rock; b is a joint attitude adjustment coefficient; c is a gravity adjustment coefficient;
9.1.2: the rock stress coefficient A is determined by the ratio of the uniaxial compressive strength of the intact rock to the ground stress generated by stope centerline mining, and the formula is as follows:
when in use
A＝0.1
When in use
When in use
A＝1
9.1.3: the value of the joint attitude adjustment coefficient B is measured by the difference between the dip angle of the stope face and the dip angle of the main joint group;
9.1.4: the gravity adjustment coefficient C reflects the influence of the stope face attitude on the stope rock stability under the influence of gravity, the size of the gravity adjustment coefficient C depends on the caving, the sliding and the sliding of the side slope and the like of the exposed surface of the stope roof, and the relationship between the gravity adjustment coefficient C and the stope surface inclination angle alpha is determined by the following formula:
C＝8-7cosα
9.2: the improvement method of Mathews stable graph method comprises the following steps:
9.2.1: when solving for the Q' value, use RQDtReplace RQD values recorded for the borehole core and apply the RQD valuestTaking into account the anisotropic character of (A), according to RQDtSolving the formula for anisotropy to obtain the RQDtThe Q' value of the anisotropic characteristic, the formula is as follows:
in the formula: RQDtθIs RQD at angle thetatA value; RQDtminIs the RQD minimum at threshold t;
is RQDtThe mean value of (a); a and b are related fitting coefficients; j. the design is a squarenThe number of joints is; j. the design is a squarerThe joint roughness coefficient; j. the design is a squareaIs the joint alteration coefficient; j. the design is a squarewTo save the water reduction coefficient; SRF is the stress reduction coefficient;
9.2.2: the ground stress value is obtained by measuring through an acoustic emission method, a relation curve of the ground stress and the burial depth h is obtained according to a fitting curve of measured data of measuring points at different depths, and the ground stress value of any burial depth position and any mineral house is obtained;
9.2.3: and (3) redrawing the Mathews stability graph by combining the equal probability contour graph, wherein the curve formula of the stability region-destruction boundary is as follows:
logN＝1.8206logR+1.618
the curve formula for the break-breakout boundary is as follows:
logN＝1.8076logR-3
9.2.4: the improved Mathews stability map is divided into three regions: a stable zone, a failure or major failure zone and a breakout zone, projects on a stable-failure boundary line, stopes with 57% probability stable, 43% probability failure, 0% probability breakout; engineering on the caving-damage boundary line, stope is stable with 0% probability, damaged with 5% probability and caving with 95% probability; the engineering in the stable region is stable; in the engineering of a damage or main damage area, the probability is stable at 57-0%, the probability is damaged at 43-5%, and the probability is collapsed at 0-95%; the project in the breakout area will continue to have breakout.
5. The dynamic evaluation method for the stability of the surrounding rock based on laser scanning, BQ and numerical simulation as claimed in claim 1, wherein in the step (7), RQD is usedtThe process of the anisotropy solution method is as follows:
7.1: solving RQD under the optimal threshold t based on the optimal threshold t and three dissected two-dimensional fracture network modelstValue, each fracture network model solves 36 RQDstA value;
7.2: calculating RQD on each fracture network modeltMaximum value RQD oftmaxMinimum RQDtminSum mean value
7.3: according to study Direction and RQDtminDirection of orientationCorrecting the relationship by combining the variance, and providing RQD under the anisotropic conditiontThe calculation formula of (a) is as follows:
in the formula:
is RQDtMean value of, RQD'tIs RQDtThe correction coefficient of (2);
7.4: correction coefficient RQD'tSolution of (a) takes into account the RQD of the study direction at a particular azimuthtθThe value and variance D, corrected as follows, when RQDtθ＝RQDtminThen, RQD't-D is variance; when RQD is reachedtθ＝RQDtmaxThen, RQD'tD; when in use
Then, RQD't＝0；
7.5: according to a correction coefficient RQD'tThe solution process of (2) provides RQD'tThe correction formula is as follows:
in the formula: RQDtθRQD at angle θtValue, RQDtminIs RQD below threshold ttMinimum value, a, b are correlation coefficients;
7.6: propose RQDtThe anisotropy calculation formula is as follows:
7.7: solving the correlation coefficients a and b, wherein the process is as follows:
7.7.1: calculating RQD on each fracture network modeltMaximum value RQD oftmaxMinimum RQDtminSum mean value
7.7.2: according to RQDtθWith RQDtmax、RQDtminAnd
and the variance are obtained to obtain three groups of RQD'tValues, and based on three groups RQD'tDrawing a scatter diagram;
7.7.3: a fitting curve is made according to the scatter diagram, the intercept of the curve is a value a, the slope is b, and the equation of the curve is RQD'tThe correction formula of (2);
7.8: the solved values of a and b are substituted into a correction coefficient formula and RQDtIn the anisotropy calculation formula, the RQD related to the angle theta is obtainedtAn anisotropy formula;
7.9: according to RQDtThe anisotropy formula is used for solving the RQD of any angletThe value is obtained.
6. The dynamic evaluation method for the stability of the surrounding rock based on laser scanning, BQ and numerical simulation as claimed in claim 1, wherein in the step (6), the process of the optimal threshold t solution method based on BQ inversion is as follows:
6.1: RQD inversion based on BQ indextThe range, process is as follows:
and (3) searching a rock quality index table by combining the rock quality grades calculated by BQ grading, and performing inversion to determine RQD (maximum-degree-of-saturation) at the rock quality gradetA range value;
6.2: the optimal threshold value t solving method comprises the following steps:
6.2.1: cutting three sections at any angle on the three-dimensional fracture network model through a central point O to obtain three two-dimensional fracture network models, and deriving the two-dimensional fracture network models and data;
6.2.2: setting different thresholds t aiming at each two-dimensional fracture network model, and solving RQD (maximum likelihood ratio) under different thresholds ttA value;
6.2.3: deriving RQD at different thresholds ttValue, calculating RQDtMean value;
6.2.4: the threshold value t is used as the abscissa, and RQD is usedtIs the ordinate, draws RQDtScatter plots that vary with threshold t;
6.2.5: according to the scatter diagram, setting a fitting equation and fitting RQDtA plot of variation with threshold t;
6.2.6: RQD to be invertedtRange values, brought into fitted RQDtIn the curve graph changing along with the threshold t, the function equation and the curve graph are combined to solve the RQDtObtaining three groups of threshold value t ranges in the range of the threshold value t;
6.2.7: taking the intersection of the ranges of the three groups of threshold values t as the range of the optimal threshold value t;
6.2.8: taking a midpoint value in the range of the optimal threshold t as the optimal threshold t to obtain the optimal threshold t;
6.2.9: and outputting the range of the optimal threshold value t and the value of the optimal threshold value t.
7. The dynamic evaluation method for the stability of the surrounding rock based on laser scanning, BQ and numerical simulation as claimed in claim 1, wherein in the step (5), RQD is usedtThe process of anisotropy mapping is as follows:
5.1：RQDtsolving calculation, the process is as follows:
5.1.1：RQDtthe theoretical formula is as follows:
in the formula: x is the number ofiRepresenting the length of the whole rock or interval, RQD, of the ith greater than a given threshold t, in a certain line directiontRepresenting the rock quality index corresponding to the threshold t, i.e. RQD at the threshold ttA value;
5.1.2: determining a section center point O, a length a and a width b of the two-dimensional fracture network model, solving a center point O coordinate (X) by taking the lower left corner of the model as a coordinate origin, the horizontal right side as an X axis and the vertical upward side as a y axis0，Y0) And a boundary equation;
5.1.3: taking O as a starting point, drawing a measuring line at an angle of every alpha of 10 degrees, intersecting with the fracture network model, drawing 36 measuring lines in total, wherein the length L of the measuring line is equal to the distance from the O point to the boundary of the fracture network model, and is represented by L0-L35, and solving a measuring line equation;
5.1.4: judging the intersection point of the measuring line and the boundary, and setting the intersection point of the measuring line equation and the boundary equation as (x)a，ya) Sequentially connecting a measuring line equation and a boundary equation, and judging whether the measuring line is intersected with the boundary;
5.1.5: solving the intersection point (X) of the boundary equation of the survey line and the fracture network modelZ，α，YZ，α)；
5.1.6: determining the interval of the measuring line, wherein the principle is as follows:
5.1.7: from the coordinates (X) of the start of each joint in the fracture network modelb，Yb) And endpoint coordinate (X)c，Yc) Establishing a corresponding analytical equation;
5.1.8: solving the intersection points of the first measuring line and each joint equation, and circularly judging each intersection point (X)j，Yj) Range, if the intersection point is coincident with a < XjC and b < YjIf d, recording the intersection point, and traversing all the joint equations to obtain all m intersection points;
5.1.9: sequencing the recorded m intersection points and the starting point coordinates and the ending point coordinates of the measuring lines from small to large according to the abscissa or the ordinate;
5.1.10: calculating the distance d between adjacent intersectionsi；
5.1.11: inputting a threshold value t;
5.1.12: circular comparison diAnd t is the size of the initial lt0, the rule is as follows:
lt＝lt+diif d isi＞t
5.1.13: solving RQD corresponding to each measuring linetValue of mαIndicates, and the corresponding distance l between the start and end of the survey lineα；
5.1.14: circularly solving m corresponding to each measuring lineαObtaining RQD in 36 directionstA value;
5.2：RQDtthe anisotropy plot was plotted as follows:
5.2.1: mixing 36 RQDstValues, ordered in order of angle;
5.2.2: drawing a circle by taking the O point as the center of the circle and 1 as the radius, and finding the distance from the center of the circle by taking the ray angle as alpha
The points of (a) are marked;
5.2.3: the end points of 36 groups of rays are connected in sequence, if the RQD on one ray istIf the value is 0, the center of the circle is taken;
5.2.4: plotting RQDtAn anisotropy map of (a);
5.2.5: output RQDtAn anisotropy map;
5.3: RQD at different thresholds ttThe anisotropy plot was plotted as follows:
5.3.1: inputting different threshold values t and solving corresponding RQDtA value;
5.3.2: RQD at different thresholds ttValue, in degrees, in RQDtThe values are functions and are drawn under the same coordinate system;
5.3.3: obtaining RQD under different threshold values ttAn anisotropy map;
5.3.4: output RQD under different thresholds ttAn anisotropy map;
5.4: output RQDtThe value of (c).
8. The dynamic evaluation method for the stability of the surrounding rock based on the laser scanning, BQ and numerical simulation as claimed in claim 1, wherein in the step (4), the generation and sectioning processes of the three-dimensional fracture network model of the rock body are as follows:
4.1: solving the random number, wherein the process is as follows:
4.1.1: the mathematical method for generating random numbers should satisfy the following conditions: the generated random number sequence is uniformly distributed in the (0, 1) interval; there should be no correlation between sequences; the random sequence has a long enough repetition period, the generation speed on a computer is high, the occupied memory space is small, and the repeatability is complete;
4.1.2: according to the Monte Carlo method, a structural surface network model obeying the model is reproduced according to the established structural surface geometric probability model, uniformly distributed random variables are generated in a (0, 1) interval, and random numbers obeying other distributions are generated by using the uniformly distributed random variables;
4.1.3: the density function of the joint geometric parameters comprises normal distribution, log-normal distribution, negative exponential distribution and uniform distribution;
4.1.4: determining basic geometric parameters for generating joints according to the obtained random numbers;
4.2: generating a rock three-dimensional fracture network model, wherein the process is as follows:
4.2.1: according to the automatic statistical result of the structural surface data and the obtained random number, storing the data of each group of structural surfaces into a text file, and expressing the data by st.dat;
4.2.2: and the content format of dat data is as follows: the coordinates (x, y and z) of the center point of the disc of each structural plane, the radius D of the disc, the inclination angle DA, the inclination DD, the trend SD, the thickness thin, the normal directions NX, NY and NZ and joints are grouped;
4.2.3: in order to distinguish the structural surfaces of different groups, the structural surface discs of the same group are endowed with the same color and are represented by a number array GID;
4.2.4: writing a program by using Matlab software, reading a structural plane data file st.dat, and automatically generating a rock three-dimensional fracture network model after running;
4.2.5: obtaining a rock three-dimensional fracture network model;
4.3: cutting a two-dimensional fracture network model, wherein the process is as follows:
4.3.1: combining a Matlab software programming tool on the three-dimensional fracture network model, and taking the central point of the three-dimensional fracture network model as a center to realize the section cutting function at any angle;
4.3.2: obtaining a two-dimensional fracture network model of any angle passing through a central point;
4.3.3: combining a Matlab software programming tool on the three-dimensional fracture network model, and realizing the section cutting function of any angle and any direction on any position of the three-dimensional fracture network model;
4.3.4: obtaining a two-dimensional fracture network model at any angle and any direction;
4.3.5: and storing the data on the cut section into an st1.dat file, wherein the section is in a three-dimensional coordinate system, and the data formats in the file are as follows from left to right: joint center point coordinates (x, y, z), joint length D, inclination DA, inclination DD, strike SD, thickness thin, normal directions NX, NY, NZ;
4.3.6: converting the three-dimensional coordinate system into a two-dimensional coordinate system, and storing the two-dimensional profile data into a st2.dat file, wherein the data formats in the file are as follows from left to right: joint center point coordinates (x, y), joint length D, inclination DA, inclination DD, thickness thin, normal directions NX, NY, NZ;
4.4: and outputting the data of the three-dimensional fracture network model and the data of any two-dimensional fracture network model.
9. The dynamic surrounding rock stability evaluation method based on laser scanning, BQ and numerical simulation as claimed in claim 1, wherein in the step (2), the process of structural plane cluster analysis is as follows:
2.1: processing point cloud data, wherein the process is as follows:
2.1.1: importing structural surface point cloud data obtained by laser scanning;
2.1.2: calculating the distance between the current point and the adjacent point in the point cloud after the topological structure and the distance mean value, and identifying and eliminating noise points in the point cloud data through a distance threshold;
2.1.3: determining a spatial three-dimensional coordinate of the point cloud data according to the self spatial coordinate position of the three-dimensional laser scanner and the occurrence orientation of the field structure surface;
2.1.4: converting point cloud data based on a lower hemisphere equal angle projection method, and converting joint attitude data represented by inclination and dip angles into structural surface attitude data represented by joint unit normal vectors;
2.1.5: obtaining structural surface data expressed by unit normal vectors;
2.2: the process of cluster analysis by a neighbor propagation algorithm is as follows:
2.2.1: let the number of actually measured samples of the structural surface be N, and the tendency of each sample data be XiThe angle of inclination is YiWith a tendency X of each sample dataiAngle of inclination YiDetermining an initial clustering center as a cluster to obtain N initial clustering centers;
2.2.2: traversing all sample data according to a similarity measurement criterion, calculating the distance between each sample data and a clustering center, and distributing each sample data to the nearest clustering center to obtain N groups of data;
2.2.3: for each group of data, solving and calculating the clustering center of each group of data by a characteristic modulus analysis method, and firstly, calculating a matrix S according to the following formula under the assumption that l data exist in a certain group:
in the formula: (x)i，yi，zi) Is a unit normal vector of any structural surface, i belongs to (1, l);
2.2.4: next, the eigenvalues (τ) of the matrix S are solved1，τ2，τ3) And a feature vector ([ xi ])1，ξ2，ξ3) In which τ is1＜τ2＜τ3Feature vector xi corresponding to the maximum feature value3For the average vector of l vectors in the group, ξ3As new clustersA center;
2.2.5: for all sample data, repeatedly calculating the distance from each sample data to the clustering center, the matrix S, the characteristic value and the characteristic vector until the positions of all the clustering centers are fixed, and determining the grouping of the structural plane;
2.2.6: converting the structural surface attitude data expressed by unit normal vector into structural surface attitude data expressed by inclination and dip angle;
2.2.7: carrying out statistical analysis on structural surface attitude data, calculating an average value m and a standard deviation sigma of the inclination angle of the structural surface, and calculating a steady interval [ m-sigma, m + sigma ] of the inclination angle data;
2.2.8: judging tendency X of initial clustering center of sample dataiAnd the inclination angle YiWhether it falls within the robust interval [ m-sigma, m + sigma ]]If yes, finishing clustering analysis; if not, re-clustering the sample data until the tendency X of the initial clustering centeriAnd the inclination angle YiAll fall within the robust interval [ m-sigma, m + sigma]；
2.3: drawing a polar diagram, wherein the process is as follows;
2.3.1: solving the orthographic projection coordinate points (x) of all the structural surface normals according to the longitudinal section principle of the structural surface space orthographic projection diagram based on the structural surface normal attitude datan，yn)；
2.3.3: drawing a base circle with the diameter as the unit length, drawing two diameters of vertical and horizontal, and marking E, S, W, N;
2.3.4: the coordinate (x) of the declination plane of all the structural surfaces is measuredn，yn) Plotted on a base circle graph;
2.3.5: drawing a polar point diagram of the structural surface;
2.4: structural surface statistical analysis, the process is as follows:
2.4.1: determining a sample partition interval m;
2.4.2: solving sample range
2.4.3: calculate each partition interval Mm：
2.4.4: determining the probability of the sample falling in each partition interval, and counting the number N of the samples falling in each interval by using a computer circulation languagemCalculating the probability P of the number of samples by combining the total number N of samplesm；
2.4.5: solving sample mean
2.4.6: solving for sample variance S2Wherein S is the standard deviation;
2.4.7: according to the probability PmDrawing the probability distribution forms of the inclination, the dip angle, the trace length, the spacing and the fault distance of each group of structural surfaces;
2.5: and outputting grouping information of the occurrence of the structural planes, wherein the grouping information comprises the mean value and the variance of the inclination, the dip angle, the trace length, the spacing and the fault distance of each group of structural planes.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011011852.2A CN112200426A (en) | 2020-09-16 | 2020-09-16 | Dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011011852.2A CN112200426A (en) | 2020-09-16 | 2020-09-16 | Dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112200426A true CN112200426A (en) | 2021-01-08 |
Family
ID=74014495
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011011852.2A Withdrawn CN112200426A (en) | 2020-09-16 | 2020-09-16 | Dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112200426A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112903565A (en) * | 2021-02-01 | 2021-06-04 | 核工业北京地质研究院 | Permeability determination method considering internal geometric characteristics of rock fracture |
| CN113343343A (en) * | 2021-07-16 | 2021-09-03 | 龙欢 | Rock slope stability evaluation method |
-
2020
- 2020-09-16 CN CN202011011852.2A patent/CN112200426A/en not_active Withdrawn
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112903565A (en) * | 2021-02-01 | 2021-06-04 | 核工业北京地质研究院 | Permeability determination method considering internal geometric characteristics of rock fracture |
| CN112903565B (en) * | 2021-02-01 | 2022-10-18 | 核工业北京地质研究院 | Permeability determination method considering internal geometric characteristics of rock fracture |
| CN113343343A (en) * | 2021-07-16 | 2021-09-03 | 龙欢 | Rock slope stability evaluation method |
| CN113343343B (en) * | 2021-07-16 | 2024-02-02 | 龙欢 | Rock slope stability evaluation method |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112150001A (en) | Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram | |
| CN112200426A (en) | Dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation | |
| CN112200419A (en) | Surrounding rock stability evaluation method based on laser scanning, BQ and improved Mathews stability diagram | |
| CN112200423A (en) | Surrounding rock stability dynamic evaluation method based on BQ and numerical simulation | |
| CN112200417A (en) | Based on photogrammetry, BQ,RQDtImproved Mathews stability chart evaluation method for ground stress | |
| CN112200431A (en) | Dynamic evaluation method for stability of empty zone based on laser scanning, BQ and numerical simulation | |
| CN112150002A (en) | Based on laser scanning, BQ and RQDtImproved Mathews stability chart evaluation method for ground stress | |
| CN112200418A (en) | Dynamic evaluation method for stability of surrounding rock based on photogrammetry, BQ and numerical simulation | |
| CN112149996A (en) | Based on laser scanning, BQ and RQDtImproved Mathews stability chart evaluation method for anisotropy | |
| CN112132407A (en) | Space RQD based on BQ inversion optimal threshold ttSolving method | |
| CN112132411A (en) | Based on laser scanning, BQ and RQDtMethod for solving Q anisotropy of anisotropy | |
| CN112200420A (en) | Method for evaluating stability of empty area based on laser scanning, BQ and improved Mathews stability diagram | |
| CN112200421A (en) | Based on laser scanning, BQ and RQDtImprovement method of Mathews stability diagram method of ground stress | |
| CN112464514A (en) | Based on photogrammetry, RQD and RQDtMethod for solving unfavorable position of roadway excavation | |
| CN112132408A (en) | Space RQD based on laser scanning and BQ inversion optimal threshold ttSolving method | |
| CN112150005A (en) | Based on laser scanning, BQ and RQDtImproved method of anisotropic Mathews stabilization mapping | |
| CN112464516A (en) | Space RQD based on laser scanning and RQD inversion optimal threshold ttSolving method | |
| CN112200428A (en) | Dynamic evaluation method for stability of empty zone based on photogrammetry, BQ and numerical simulation | |
| CN112200432A (en) | Dead zone stability evaluation method based on photogrammetry, BQ and improved Mathews stable diagram | |
| CN112200429A (en) | Based on BQ, RQDtImproved Mathews stability chart evaluation method for ground stress | |
| CN112150004A (en) | Based on photogrammetry, BQ, RQDtImprovement method of Mathews stability diagram method of ground stress | |
| CN112200430A (en) | Based on BQ, RQDtImprovement method of Mathews stability diagram method of ground stress | |
| CN112200416A (en) | Dynamic evaluation method for stability of empty zone based on BQ and numerical simulation | |
| CN112200427A (en) | Surrounding rock stability evaluation method based on BQ and improved Mathews stability diagram method | |
| CN112150003A (en) | Based on photogrammetry, BQ, RQDtImproved Mathews stability chart evaluation method for anisotropy |
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: 20210108 |