CN112150001A - Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram - Google Patents

Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram Download PDF

Info

Publication number
CN112150001A
CN112150001A CN202011011774.6A CN202011011774A CN112150001A CN 112150001 A CN112150001 A CN 112150001A CN 202011011774 A CN202011011774 A CN 202011011774A CN 112150001 A CN112150001 A CN 112150001A
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
Application number
CN202011011774.6A
Other languages
Chinese (zh)
Inventor
胡高建
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Shaoxing
Original Assignee
University of Shaoxing
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Shaoxing filed Critical University of Shaoxing
Priority to CN202011011774.6A priority Critical patent/CN112150001A/en
Publication of CN112150001A publication Critical patent/CN112150001A/en
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
    • G01C11/04—Interpretation of pictures
    • G—PHYSICS
    • 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
    • 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)
  • Evolutionary Computation (AREA)
  • Strategic Management (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Economics (AREA)
  • Educational Administration (AREA)
  • Development Economics (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Optimization (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Business, Economics & Management (AREA)
  • Pure & Applied Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Tourism & Hospitality (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Game Theory and Decision Science (AREA)
  • Quality & Reliability (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Architecture (AREA)
  • Structural Engineering (AREA)
  • Civil Engineering (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Multimedia (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)

Abstract

A surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability chart belongs to goaf stabilityThe field of qualitative evaluation, comprising the steps of: (1) fast acquisition of structural surface digital photogrammetry; (2) structural surface fuzzy equivalent 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 surrounding rock stability analysis method; (12) a method for evaluating stability of a dead zone. The invention realizes the consideration of RQDtThe improvement of the Mathews stability diagram method of anisotropy and ground stress fitting realizes the evaluation of the stability of the surrounding rock based on the improved Mathews stability diagram method and numerical simulation. The method is clear and suitable for evaluating the stability of the surrounding rock of the goaf.

Description

Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram
Technical Field
The invention relates to a method based on photographyThe invention relates to a method for evaluating the stability of surrounding rocks of measurement, BQ and improved Mathews stable graphs, in particular to a method for calculating RQD by inversion based on photogrammetry, fuzzy equivalent analysis, BQ indexes, fracture network models and 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 surrounding rock stability evaluation method based on photogrammetry, BQ and the improved Mathews stability diagram 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 can be formed in the metal mine, the goafs are mutually communicated, overlapped and concentrated underground, stress disturbance and mutual restriction are caused, a complex group space relationship is formed, and the stable state of a rock body is extremely complex.
Along with the increase of the development of deep mineral resources in China, the number of goafs is continuously increased, the scale is continuously increased, the mining conditions are continuously worsened, the safety production situation is more severe, and the accident starting number and the casualties are always high. Particularly since the 80 s of the 20 th century, because the mining order is disordered, a large amount of blind goafs are left by illegal and unplanned disorderly mining and abusing, and the method becomes one of the main dangerous sources influencing the safety production of mines. In the goaf, due to the fact that the goaf exists for a long time, the rock mass structure is dangerous after being damaged for many times, the possibility and uncertainty of disasters such as rock mass instability, surface subsidence, rock burst, roadway damage and the like are gradually increased along with the continuous increase of the mining depth, and the difficulty of subsequent ore body mining is aggravated. Therefore, the analysis and evaluation of the goaf stability are significant.
The Mathews stability diagram method is one of goaf surrounding rock stability analysis methods, and is a stability diagram method based on a Q system, which is proposed in 1980 by Mathews et al through 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 RQD is obtained by adopting a drilling mode, the results obtained by different rock body parts are completely different, the RQD value result 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 parallel to the main joint groupPerpendicular to the main joint group direction, lower RQD values are 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 photogrammetry technology is a brand-new method for quickly, efficiently, accurately and comprehensively acquiring the random rock mass structural plane information, and is particularly advanced in the aspect of solving the structural plane direction and scale information. Cluster analysis is a method for statistically studying the classification problem, and its task is to distribute all the sample data to several samplesDry clusters, such that sample data of the same cluster is clustered around the cluster center, relatively close in distance, while sample data of different clusters are relatively far apart. Therefore, through the digital photogrammetry technology and fuzzy equivalent cluster analysis, the rapid acquisition and rapid processing analysis of structural plane data can be realized, and the structural plane data 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. Some scholars respectively calculate fractal dimension of structural surface distribution by applying fractal theory based on three-dimensional structural surface network simulation technology, and obtain RQDs under different thresholds by continuously changing thresholds of generalized RQDstAnd 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 surrounding rock stability evaluation method based on photogrammetry, BQ and an improved Mathews stability diagram.
Disclosure of Invention
In order to consider the anisotropy characteristic of the rock mass, the change characteristic of the ground stress along with the burial depth and the destruction characteristic of the rock mass in the goaf stability evaluation, the invention provides a surrounding rock stability evaluation method based on photogrammetry, BQ and an improved Mathews stability diagram. The method is based on photogrammetry, fuzzy equivalence analysis, BQ indexes, a fracture network model and a generalized RQD theory, and RQD is calculated through inversiontThe 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 surrounding rock stability evaluation method based on photogrammetry, BQ and the improved Mathews stability diagram by combining a numerical simulation method.
In order to solve the technical problems, the invention provides the following technical scheme:
a method for evaluating the stability of surrounding rocks based on photogrammetry, BQ, and modified Mathews stability charts, the method comprising the steps of:
(1) the method comprises the following steps of (1) quickly obtaining the digital photogrammetry of the structural surface:
1.1: selecting a rock mass with good surface joint development and no obstacles as a photogrammetric region according to the rock mass range and the spatial position of the observation region, and vertically erecting a marker post on one side of the measurement region for calibrating the distance between any two points on the finally generated three-dimensional image;
1.2: selecting a structural surface which is exposed, has a large area and is smooth on the surface of the rock mass as a calibration point, measuring the inclination and the dip angle by using a compass, and marking the structural surface for orientation reality of an image during post-processing;
1.3: sequentially photographing a rock mass at the left and right positions right in front of the selected region by using a high-resolution camera, wherein the distance D between a lens and the measured rock mass and the distance B between two imaging positions satisfy the relation B of D/8-D/5 when the two times of photographing are carried out;
1.4: after the measurement point data are collected, the marker post is taken back, and the marker post returns to the room for further post-processing operation;
1.5: importing the left and right views obtained by field photogrammetry into a software analysis system, matching the pixels in the left and right views by adopting reference calibration, pixel matching and image deformation correction to synthesize a three-dimensional solid model of the rock surface;
1.6: according to the size of the marker post and the appearance of the calibration point measured by the compass, the orientation, the size and the distance of the three-dimensional solid model are realized;
1.7: identifying and positioning each structural surface based on a realistic entity model, and deriving structural surface data information;
(2) structural surface fuzzy equivalent 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:
Figure BSA0000220400580000041
in the formula: jv is the rock volume rational number in unitsStrip/m3；
3.2: the Jv calculation is as follows:
Figure BSA0000220400580000042
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 method for analyzing the stability of the surrounding rock comprises the following steps:
11.1: establishing a three-dimensional model of the mineral house by using CAD software according to the plan view and the section view of the mineral house;
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 stoping of the chamber, firstly carrying out numerical simulation of a non-stoping state to obtain an initial stress field, and then carrying out numerical simulation of stoping of the chamber to obtain a stress field of a goaf;
11.5: analyzing the simulation result after stoping of the chamber, wherein the analysis area comprises a roof of the chamber, surrounding rocks of an upper tray and a lower tray and a pillar;
11.6: 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:
Figure BSA0000220400580000061
in the formula: sigma1、σ3Maximum principal stress and minimum principal stress, respectively; c is cohesion;
Figure BSA0000220400580000062
is a friction angle;
11.7: 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.8: according to the analysis result, obtaining the damage modes of the roof, the surrounding rocks of the upper and lower trays and the ore pillars of the chamber;
(12) a method for evaluating stability of a dead zone.
Further, in the step (12), the procedure of the method for evaluating 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: correcting the initial evaluation result by using a surrounding rock stability analysis result 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 rocks of the goaf evaluated by the Mathews stability diagram method are in a stable region, if the top plate, the stud, the surrounding rocks of the upper and lower discs of the goaf obtained by numerical simulation are not damaged, the goaf is in a stable state, and if the top plate, the stud, the surrounding rocks of the upper and lower discs of the goaf obtained by numerical simulation are damaged, the goaf is possibly in a damaged state;
12.3: and obtaining the stability evaluation result 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: measuringMeasuring and solving joint group number J of upper and lower plates of a mineral housen(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
Figure BSA0000220400580000072
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
Figure BSA0000220400580000071
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
Figure BSA0000220400580000081
When in use
Figure BSA0000220400580000082
When in use
Figure BSA0000220400580000083
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:
Figure BSA0000220400580000084
in the formula: RQDtθIs RQD at angle thetatA value; RQDtminIs the RQD minimum at threshold t;
Figure BSA0000220400580000085
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, per crackThe network model solves 36 RQDstA value;
7.2: calculating RQD on each fracture network modeltMaximum value RQD oftmaxMinimum RQDtminSum mean value
Figure BSA0000220400580000093
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:
Figure BSA0000220400580000094
in the formula:
Figure BSA0000220400580000095
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
Figure BSA0000220400580000096
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:
Figure BSA0000220400580000091
in the formula: RQDtθRQD at angle θtValue, RQDtminIs RQD below threshold ttMinimum value, a, b are correlation coefficients;
7.6: lifting deviceDischarge RQDtThe anisotropy calculation formula is as follows:
Figure BSA0000220400580000092
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
Figure BSA0000220400580000097
7.7.2: according to RQDtθWith RQDtmax、RQDtminAnd
Figure BSA0000220400580000098
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 (4) searching a rock quality index table by combining the rock quality grades calculated by BQ grading to determine RQD (total green weight) of the rock mass 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:
Figure BSA0000220400580000101
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 value t, i.e.RQD at 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:
Figure BSA0000220400580000102
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
Figure BSA0000220400580000111
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 fuzzy equivalent clustering analysis of the structural plane is as follows:
2.1: identifying structural surface data obtained by digital photogrammetry;
2.2: structural plane clustering analysis based on a fuzzy equivalent clustering algorithm comprises the following steps:
2.2.1: let the number of actual measurement samples of the structural surface be N, and the ith sample be (x)i1，xi2)，xi1For structural orientation, xi2The inclination angle of the structural surface is shown;
2.2.2: calculating fuzzy relation matrix R and similarity coefficient Rij；
2.2.3: solving the closure t (R);
2.2.4: carrying out grouping judgment on the structural surface;
2.3: drawing a pole point diagram of the structural surface, wherein the process is as follows:
2.3.1: drawing a polar point diagram by adopting a lower hemisphere equal-angle projection method;
2.3.2: converting the joint attitude data expressed by the inclination and the dip into structural plane attitude data expressed by a joint unit normal vector to obtain structural plane attitude data expressed by the unit normal vector;
2.3.3: solving the erythroplanic projection coordinate points of all the structural surface normals;
2.3.4: 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.5: drawing the stereographic projection coordinates of all the structural surfaces on a base circle diagram to realize the drawing of a polar point diagram of the structural surfaces;
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, firstly counting the number of the samples falling in each interval by utilizing a computer circulation language, and calculating the probability of the number of the samples by combining the total number N of the samples;
2.4.5: solving a sample mean value and a sample variance;
2.4.6: drawing probability distribution forms of the inclination, the inclination angle, the trace length, the spacing and the fault distance of each group of structural surfaces according to the probability values;
2.5: and outputting structural surface data and probability distribution morphology.
The invention has the following beneficial effects:
1. the method is based on photogrammetry, fuzzy equivalence analysis, BQ indexes, a fracture network model, a generalized RQD theory and RQDtAn anisotropy solving method, a ground stress measuring method, a method for improving Mathews stability diagram and a numerical simulation method are providedA surrounding rock stability evaluation method based on photogrammetry, BQ and an improved Mathews stability diagram;
2. the invention realizes the rapid photogrammetry and fuzzy equivalent clustering analysis of the structural surface;
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 consideration of RQDtThe 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 the stability analysis of the surrounding rock based on numerical simulation;
12. the method provided by the invention realizes the consideration of the anisotropic characteristic of rock mass quality, the change characteristic of ground stress along with the burial depth and the destruction characteristic of the rock mass in the evaluation of the stability of the surrounding rock of the dead zone.
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 evaluating the stability of surrounding rocks based on photogrammetry, BQ and an improved Mathews stability map comprises the following steps:
1) the method comprises the following steps of (1) quickly obtaining the digital photogrammetry of the structural surface:
1.1: selecting a rock mass with good surface joint development and no obstacles as a photogrammetric region according to the rock mass range and the spatial position of the observation region, and vertically erecting a marker post on one side of the measurement region for calibrating the distance between any two points on the finally generated three-dimensional image;
1.2: selecting a structural surface which is exposed, has a large area and is smooth on the surface of the rock mass as a calibration point, measuring the inclination and the dip angle by using a compass, and marking the structural surface for orientation reality of an image during post-processing;
1.3: sequentially photographing a rock mass at the left and right positions right in front of the selected region by using a high-resolution camera, wherein the distance D between a lens and the measured rock mass and the distance B between two imaging positions satisfy the relation B of D/8-D/5 when the two times of photographing are carried out;
1.4: after the measurement point data are collected, the marker post is taken back, and the marker post returns to the room for further post-processing operation;
1.5: importing the left and right views obtained by field photogrammetry into a software analysis system, matching the pixels in the left and right views by adopting reference calibration, pixel matching and image deformation correction to synthesize a three-dimensional solid model of the rock surface;
1.6: according to the size of the marker post and the appearance of the calibration point measured by the compass, the orientation, the size and the distance of the three-dimensional solid model are realized;
1.7: identifying and positioning each structural surface based on a realistic entity model, and deriving structural surface data information;
2) structural surface fuzzy equivalent clustering analysis, the process is as follows:
2.1: identifying structural surface data obtained by digital photogrammetry;
2.2: structural plane clustering analysis based on a fuzzy equivalent clustering algorithm comprises the following steps:
2.2.1: let the number of actual measurement samples of the structural surface be N, and the ith sample be (x)i1，xi2)，xi1For structural orientation, xi2The structural plane inclination angle is shown as the fuzzy relation matrix R:
Figure BSA0000220400580000151
element r in the matrixijRepresenting the similarity degree between the ith sample and the jth sample for the similarity coefficient between the ith sample and the jth sample; r isijLarger indicates that sample i is more similar to sample j;
2.2.2: calculating a similarity coefficient rij：
Figure BSA0000220400580000152
In the formula: n ═ 1, 2.. N; n ═ 1, 2.. N; c is a calculation parameter (c is more than or equal to 0 and less than or equal to 1), and the value of c is properly selected to ensure that r isijIn [0, 1 ]]Disperse from the middle;
2.2.3: solving the closure t (R):
R2＝RR
R4＝R2R2
… (3)
2.2.4: and (3) carrying out structural surface grouping judgment: the fuzzy matrix multiplication steps are similar to the common matrix multiplication, and the difference is that the multiplication of two terms is not carried out first and then the addition is carried out, but the multiplication is carried out first and then the multiplication is carried out; if C ═ AB, then the elements in C
Figure BSA0000220400580000153
n-level fuzzy relation matrix R is n R continuous multiplication; namely, it is
Figure BSA0000220400580000154
When R isn＝Rn+1＝Rn+2When … (5)
There, the fuzzy equivalence matrix t (R) ═ Rn (6)
Taking the definite intercept set level lambda belongs to [0, 1 ]]If r in t (R)ijIf the structural plane i and the structural plane j belong to the same class, the structural plane i and the structural plane j belong to the same class; namely, it is
rij≥λ (7)
2.3: drawing a polar diagram, wherein the process is as follows;
2.3.1: drawing a polar point diagram by adopting a lower hemisphere equal-angle projection method;
2.3.2: 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) (8)
Figure BSA0000220400580000161
Figure BSA0000220400580000162
αd∈(0，360)，βd∈(0，90) (11)
2.3.3: obtaining structural surface attitude data expressed by unit normal vectors;
2.3.4: based on normal attitude data of the structural plane and according to a longitudinal section schematic diagram of a spatial red-flat projection diagram of the structural plane, the point A' is the red-flat projection of the plane normal; combining with the schematic diagram of the red projection, the coordinate x of A' on the schematic diagram of the red projection is calculatednAnd ynThe formula is as follows:
Figure BSA0000220400580000163
Figure BSA0000220400580000164
2.3.5: solving the coordinate point (x) of the declination projection of all the structural surface normalsn，yn)；
2.3.6: 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.7: the coordinate (x) of the declination plane of all the structural surfaces is measuredn，yn) Plotted on a base circle graph;
2.3.8: 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
Figure BSA0000220400580000171
Figure BSA0000220400580000172
2.4.3: calculate each partition interval Mm：
Figure BSA0000220400580000173
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：
Figure BSA0000220400580000174
2.4.5: solving sample mean
Figure BSA0000220400580000175
Figure BSA0000220400580000176
2.4.6: solving for sample variance S2Wherein S is the standard deviation:
Figure BSA0000220400580000177
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 classification information of the occurrence of the structural surfaces, including the mean value and 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:
Figure BSA0000220400580000178
in the formula: jv is the volume regulating number of rock mass, unit bar/m3；
3.2: the Jv calculation is as follows:
Figure BSA0000220400580000179
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 (21)
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) (22)
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:
Figure BSA0000220400580000191
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:
Figure BSA0000220400580000192
Figure BSA0000220400580000193
the boundary equation is:
Figure BSA0000220400580000194
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:
Figure BSA0000220400580000201
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:
Figure BSA0000220400580000202
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:
Figure BSA0000220400580000203
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:
Figure BSA0000220400580000204
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:
Figure BSA0000220400580000205
x0＝a/2 (32)
y0＝b/2 (33)
xm+1＝XZ，α (34)
ym+1＝YZ，α (35)
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 (36)
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:
Figure BSA0000220400580000211
Figure BSA0000220400580000212
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
Figure BSA0000220400580000213
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, see table 1, and determining RQD under the rock mass gradetRange value, RQDtThe range is 50 to 75 percent
Figure BSA0000220400580000214
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
Figure BSA0000220400580000221
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
Figure BSA0000220400580000223
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:
Figure BSA0000220400580000224
in the formula:
Figure BSA0000220400580000225
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
Figure BSA0000220400580000226
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:
Figure BSA0000220400580000222
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:
Figure BSA0000220400580000231
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
Figure BSA0000220400580000234
7.7.2: according to RQDtθWith RQDtmax、RQDtminAnd
Figure BSA0000220400580000235
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:
Figure BSA0000220400580000232
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:
Figure BSA0000220400580000233
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 (44)
Figure BSA0000220400580000241
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
Figure BSA0000220400580000242
When in use
Figure BSA0000220400580000243
When in use
Figure BSA0000220400580000244
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α (49)
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:
Figure BSA0000220400580000251
in the formula: RQDtθIs RQD at angle thetatA value; RQDtminIs the RQD minimum at threshold t;
Figure BSA0000220400580000252
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 (51)
the curve formula for the break-breakout boundary is as follows:
logN＝1.8076logR-3 (52)
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
Figure BSA0000220400580000253
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 surrounding rocks on the upper and lower walls of the chamber, wherein 3200 chambers and 3202 chambers are located in stable areas, the evaluation results of the surrounding rocks of the chamber are stable, 3201 chambers are located in a damaged or main damaged area, the probability of the damage of the surrounding rocks is 46%, the probability of the collapse of the surrounding rocks is 46%, and the probability of the stability of the surrounding rocks is 8%;
11) the method for analyzing the stability of the surrounding rock comprises the following steps:
11.1: establishing a three-dimensional model of the mineral house by using CAD software according to the plan view and the section view of the mineral house;
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 stoping of the chamber, firstly carrying out numerical simulation of a non-stoping state to obtain an initial stress field, and then carrying out numerical simulation of stoping of the chamber to obtain a stress field of a goaf;
11.5: analyzing the simulation result after stoping of the chamber, wherein the analysis area comprises a roof of the chamber, surrounding rocks of an upper tray and a lower tray and a pillar;
11.6: 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:
Figure BSA0000220400580000261
in the formula: sigma1、σ3Maximum principal stress and minimum principal stress, respectively; c is cohesion;
Figure BSA0000220400580000262
is a friction angle;
11.7: 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 (54)
in the formula: sigmatIs a rock mass placeA tensile stress; rtThe tensile strength of the rock mass;
11.8: according to the analysis result, obtaining the damage modes of the roof, the surrounding rocks of the upper and lower trays and the ore pillars of the chamber;
12) the method for evaluating 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: correcting the initial evaluation result by using a surrounding rock stability analysis result 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 rocks of the goaf evaluated by the Mathews stability diagram method are in a stable region, if the top plate, the stud, the surrounding rocks of the upper and lower discs of the goaf obtained by numerical simulation are not damaged, the goaf is in a stable state, and if the top plate, the stud, the surrounding rocks of the upper and lower discs of the goaf obtained by numerical simulation are damaged, the goaf is possibly in a damaged state;
12.3: and obtaining the stability evaluation result of the goaf.

Claims (9)

1. A method for evaluating the stability of surrounding rocks based on photogrammetry, BQ and an improved Mathews stability chart is characterized by comprising the following steps:
(1) the method comprises the following steps of (1) quickly obtaining the digital photogrammetry of the structural surface:
1.1: selecting a rock mass with good surface joint development and no obstacles as a photogrammetric region according to the rock mass range and the spatial position of the observation region, and vertically erecting a marker post on one side of the measurement region for calibrating the distance between any two points on the finally generated three-dimensional image;
1.2: selecting a structural surface which is exposed, has a large area and is smooth on the surface of the rock mass as a calibration point, measuring the inclination and the dip angle by using a compass, and marking the structural surface for orientation reality of an image during post-processing;
1.3: sequentially photographing a rock mass at the left and right positions right in front of the selected region by using a high-resolution camera, wherein the distance D between a lens and the measured rock mass and the distance B between two imaging positions satisfy the relation B of D/8-D/5 when the two times of photographing are carried out;
1.4: after the measurement point data are collected, the marker post is taken back, and the marker post returns to the room for further post-processing operation;
1.5: importing the left and right views obtained by field photogrammetry into a software analysis system, matching the pixels in the left and right views by adopting reference calibration, pixel matching and image deformation correction to synthesize a three-dimensional solid model of the rock surface;
1.6: according to the size of the marker post and the appearance of the calibration point measured by the compass, the orientation, the size and the distance of the three-dimensional solid model are realized;
1.7: identifying and positioning each structural surface based on a realistic entity model, and deriving structural surface data information;
(2) structural surface fuzzy equivalent 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:
Figure FSA0000220400570000011
in the formula: jv is the volume regulating number of rock mass, unit bar/m3；
3.2: the Jv calculation is as follows:
Figure FSA0000220400570000012
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 method for analyzing the stability of the surrounding rock comprises the following steps:
11.1: establishing a three-dimensional model of the mineral house by using CAD software according to the plan view and the section view of the mineral house;
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 stoping of the chamber, firstly carrying out numerical simulation of a non-stoping state to obtain an initial stress field, and then carrying out numerical simulation of stoping of the chamber to obtain a stress field of a goaf;
11.5: analyzing the simulation result after stoping of the chamber, wherein the analysis area comprises a roof of the chamber, surrounding rocks of an upper tray and a lower tray and a pillar;
11.6: 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:
Figure FSA0000220400570000031
in the formula: sigma1、σ3Maximum principal stress and minimum principal stress, respectively; c is cohesion;
Figure FSA0000220400570000032
is a friction angle;
11.7: 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.8: according to the analysis result, obtaining the damage modes of the roof, the surrounding rocks of the upper and lower trays and the ore pillars of the chamber;
(12) a method for evaluating stability of a dead zone.
2. The method for evaluating the stability of surrounding rocks based on photogrammetry, BQ, modified Mathews stability charts of claim 1, wherein in the step (12), the procedure of the dead zone stability 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: correcting the initial evaluation result by using a surrounding rock stability analysis result 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 rocks of the goaf evaluated by the Mathews stability diagram method are in a stable region, if the top plate, the stud, the surrounding rocks of the upper and lower discs of the goaf obtained by numerical simulation are not damaged, the goaf is in a stable state, and if the top plate, the stud, the surrounding rocks of the upper and lower discs of the goaf obtained by numerical simulation are damaged, the goaf is possibly in a damaged state;
12.3: and obtaining the stability evaluation result of the goaf.
3. The method for evaluating the stability of surrounding rocks based on photogrammetry, BQ, modified Mathews charts of claim 1, wherein in the step (10), the process of the modified Mathews chart evaluating 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 isA stress reduction factor;
10.5: solving RQD under the optimal threshold ttθ、RQDtminAnd
Figure FSA0000220400570000041
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 method for evaluating the stability of surrounding rocks based on photogrammetry, BQ, modified Mathews's stability map as set forth in claim 1, wherein the modified Mathews' stability map method in the step (9) is performed 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
Figure FSA0000220400570000051
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
Figure FSA0000220400570000052
A＝0.1
When in use
Figure FSA0000220400570000053
When in use
Figure FSA0000220400570000054
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:
Figure FSA0000220400570000055
in the formula: RQDtθIs RQD at angle thetatA value; RQDtminIs the RQD minimum at threshold t;
Figure FSA0000220400570000056
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 method for evaluating the stability of surrounding rock based on photogrammetry, BQ, improved Mathews stability charts of claim 1, wherein in 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
Figure FSA0000220400570000066
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:
Figure FSA0000220400570000061
in the formula:
Figure FSA0000220400570000062
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
Figure FSA0000220400570000063
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:
Figure FSA0000220400570000064
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:
Figure FSA0000220400570000065
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
Figure FSA0000220400570000067
7.7.2: according to RQDtθWith RQDtmax、RQDtminAnd
Figure FSA0000220400570000068
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 method for evaluating the stability of surrounding rocks based on photogrammetry, BQ and improved Mathews stability map 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 (4) searching a rock quality index table by combining the rock quality grades calculated by BQ grading to determine RQD (total green weight) of the rock mass 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 method for evaluating the stability of surrounding rock based on photogrammetry, BQ, improved Mathews stability charts of claim 1, wherein in 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:
Figure FSA0000220400570000071
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:
Figure FSA0000220400570000081
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 corresponding analytic partyA process;
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+dtif 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
Figure FSA0000220400570000082
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 method for evaluating the stability of the surrounding rock based on photogrammetry, BQ and Mathews-improved stability diagram in the step (4), wherein the generation and sectioning 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 method for evaluating the stability of surrounding rocks based on photogrammetry, BQ and improved Mathews stability charts in claim 1, wherein in the step (2), the fuzzy equivalent clustering analysis of the structural planes is performed as follows:
2.1: identifying structural surface data obtained by digital photogrammetry;
2.2: structural plane clustering analysis based on a fuzzy equivalent clustering algorithm, wherein the process is as follows;
2.2.1: let the number of actual measurement samples of the structural surface be N, and the ith sample be (x)i1，xi2)，xi1For structural orientation, xi2The inclination angle of the structural surface is shown;
2.2.2: calculating fuzzy relation matrix R and similarity coefficient Rij；
2.2.3: solving the closure t (R);
2.2.4: carrying out grouping judgment on the structural surface;
2.3: drawing a pole point diagram of the structural surface, wherein the process is as follows:
2.3.1: drawing a polar point diagram by adopting a lower hemisphere equal-angle projection method;
2.3.2: converting the joint attitude data expressed by the inclination and the dip into structural plane attitude data expressed by a joint unit normal vector to obtain structural plane attitude data expressed by the unit normal vector;
2.3.3: solving the erythroplanic projection coordinate points of all the structural surface normals;
2.3.4: 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.5: drawing the stereographic projection coordinates of all the structural surfaces on a base circle diagram to realize the drawing of a polar point diagram of the structural surfaces;
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, firstly counting the number of the samples falling in each interval by utilizing a computer circulation language, and calculating the probability of the number of the samples by combining the total number N of the samples;
2.4.5: solving a sample mean value and a sample variance;
2.4.6: drawing probability distribution forms of the inclination, the inclination angle, the trace length, the spacing and the fault distance of each group of structural surfaces according to the probability values;
2.5: and outputting structural surface data and probability distribution morphology.
CN202011011774.6A 2020-09-16 2020-09-16 Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram Withdrawn CN112150001A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011011774.6A CN112150001A (en) 2020-09-16 2020-09-16 Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011011774.6A CN112150001A (en) 2020-09-16 2020-09-16 Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram

Publications (1)

Publication Number Publication Date
CN112150001A true CN112150001A (en) 2020-12-29

Family

ID=73896222

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011011774.6A Withdrawn CN112150001A (en) 2020-09-16 2020-09-16 Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram

Country Status (1)

Country Link
CN (1) CN112150001A (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112989481A (en) * 2021-05-17 2021-06-18 中国科学院武汉岩土力学研究所 Method for processing stable visual image data of complex geological tunnel construction surrounding rock
CN113280951A (en) * 2021-07-22 2021-08-20 中国科学院地质与地球物理研究所 Method for establishing stress field distribution of sloping field in canyon region
CN114119508A (en) * 2021-11-10 2022-03-01 广东粤海珠三角供水有限公司 Shield tunnel surrounding rock quality judgment method based on monitoring video

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112989481A (en) * 2021-05-17 2021-06-18 中国科学院武汉岩土力学研究所 Method for processing stable visual image data of complex geological tunnel construction surrounding rock
CN113280951A (en) * 2021-07-22 2021-08-20 中国科学院地质与地球物理研究所 Method for establishing stress field distribution of sloping field in canyon region
CN113280951B (en) * 2021-07-22 2021-10-08 中国科学院地质与地球物理研究所 Method for establishing stress field distribution of sloping field in canyon region
CN114119508A (en) * 2021-11-10 2022-03-01 广东粤海珠三角供水有限公司 Shield tunnel surrounding rock quality judgment method based on monitoring video

Similar Documents

Publication Publication Date Title
CN112150001A (en) Surrounding rock stability evaluation method based on photogrammetry, BQ and improved Mathews stability diagram
CN112200419A (en) Surrounding rock stability evaluation method based on laser scanning, BQ and improved Mathews stability diagram
CN112200417A (en) Based on photogrammetry, BQ,RQDtImproved Mathews stability chart evaluation method for ground stress
CN112200426A (en) Dynamic evaluation method for stability of surrounding rock based on laser scanning, BQ and numerical simulation
CN112200423A (en) Surrounding rock stability dynamic evaluation method based on BQ and numerical simulation
CN112464513A (en) RQD based on photogrammetry and RQD inversiontOptimal threshold t solving method
CN112200418A (en) Dynamic evaluation method for stability of surrounding rock based on photogrammetry, BQ and numerical simulation
CN112132407A (en) Space RQD based on BQ inversion optimal threshold ttSolving method
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
CN112149996A (en) Based on laser scanning, BQ and RQDtImproved Mathews stability chart evaluation method for anisotropy
CN112464514A (en) Based on photogrammetry, RQD and RQDtMethod for solving unfavorable position of roadway excavation
CN112200432A (en) Dead zone stability evaluation method based on photogrammetry, BQ and improved Mathews stable diagram
CN112200428A (en) Dynamic evaluation method for stability of empty zone based on photogrammetry, BQ and numerical simulation
CN112132411A (en) Based on laser scanning, BQ and RQDtMethod for solving Q anisotropy of anisotropy
CN112150004A (en) Based on photogrammetry, BQ, RQDtImprovement method of Mathews stability diagram method of ground stress
CN112132403A (en) RQD based on photogrammetry and BQ inversiontOptimal threshold t solving method
CN112200420A (en) Method for evaluating stability of empty area based on laser scanning, BQ and improved Mathews stability diagram
CN112132410A (en) Space RQD based on photogrammetry and BQ inversion optimal threshold ttSolving method
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
CN112150000A (en) Based on photogrammetry, BQ, RQDtImproved method of anisotropic Mathews stabilization mapping
CN112200429A (en) Based on BQ, RQDtImproved Mathews stability chart evaluation method for ground stress
CN112200430A (en) Based on BQ, RQDtImprovement method of Mathews stability diagram method of ground stress
CN112200421A (en) Based on laser scanning, BQ and RQDtImprovement method of Mathews stability diagram method of ground stress

Legal Events

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

Application publication date: 20201229

WW01 Invention patent application withdrawn after publication