CN116466401A - Wave velocity anisotropic rock mass fracture positioning method for cross communication area of underground engineering - Google Patents

Wave velocity anisotropic rock mass fracture positioning method for cross communication area of underground engineering Download PDF

Info

Publication number
CN116466401A
CN116466401A CN202310423762.1A CN202310423762A CN116466401A CN 116466401 A CN116466401 A CN 116466401A CN 202310423762 A CN202310423762 A CN 202310423762A CN 116466401 A CN116466401 A CN 116466401A
Authority
CN
China
Prior art keywords
grid
point
adjacent
wave velocity
value
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.)
Pending
Application number
CN202310423762.1A
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.)
Sichuan University
Changan University
Original Assignee
Sichuan University
Changan University
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 Sichuan University, Changan University filed Critical Sichuan University
Priority to CN202310423762.1A priority Critical patent/CN116466401A/en
Publication of CN116466401A publication Critical patent/CN116466401A/en
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/626—Physical property of subsurface with anisotropy
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Abstract

The invention discloses a wave velocity anisotropic rock mass fracture positioning method for an underground engineering cross communication area, which comprises the steps of obtaining coordinate point positions of all sensors and coordinate point positions of actual vibration source points; traversing and calculating the travel values from each sensor to all grid nodes in the pre-built three-dimensional Cartesian coordinate grid model, and performing difference calculation to establish a theoretical arrival time difference database; calculating the arrival time difference of each sensor based on the acquisition initial time of the P waves, and establishing an actual arrival time difference database; comparing the theoretical time difference database with the actual time difference database to obtain h function values and sequencing the h function values according to the order from small to large; and selecting the objective h function values, marking the objective grid nodes corresponding to the objective h function values, and calculating and outputting the space coordinates of the seismic source points. The method can accurately monitor the space position of the occurrence of the microseism event, and can improve the microseism source positioning accuracy of the regional rock mass with anisotropic characteristics.

Description

Wave velocity anisotropic rock mass fracture positioning method for cross communication area of underground engineering
Technical Field
The invention relates to the field of engineering safety technology monitoring, in particular to a wave velocity anisotropic rock mass fracture positioning method for a cross communication area of an underground engineering.
Background
In recent years, development and utilization of deep resources have been increasingly emphasized, and development of underground engineering and tunnel engineering has been promoted. Many underground projects must be constructed in layered rock mass because layered rock mass is commonly distributed throughout the earth's surface. The layered rock mass is formed from pre-existing virgin material on the crust. During the entire diagenetic process of the layered rock mass, geological formations in the continuity of the rock mass may cause delamination of the rock and the wave velocity of the seismic waves exhibit anisotropy. However, complex geological environments result in numerous geological disasters occurring during the process of engineering excavation. Therefore, in engineering activities, the safety monitoring of underground engineering and the disaster prevention and reduction monitoring of engineering are particularly necessary. Microseismic and geophysical-based microseismic monitoring technologies have been developed, and have been well applied to the fields of underground engineering excavation, slope stability monitoring and the like. The accurate earthquake focus positioning algorithm is used, and the acquisition of accurate cracking occurrence time and spatial position is an important point of a microseismic monitoring technology.
The microseismic monitoring technology is a dynamic, continuous, real-time monitoring technology that utilizes elastic wave signals generated during the breaking process of a rock mass to study and evaluate the stability of the rock mass. The seismic source positioning is a process that a seismic source positioning system utilizes the acquired microseism signals to invert the space position and the earthquake starting moment of a microseism event in the monitoring area range through a seismic source positioning algorithm, and evaluates the reliability of a seismic source positioning result. Currently, common source positioning methods include time-of-day source positioning and three-axis sensor-based source positioning. The seismic source positioning method based on the arrival time theory is various, wherein the algorithms such as the classical Geiger algorithm, the Nelder simplex method and the like are widely applied, and good positioning effect is obtained under the research and improvement of a plurality of students. The three-axis sensor-based seismic source positioning method can relatively stably acquire positioning results in consideration of arrival time of P waves and S waves, but has limitation in a small-scale microseism monitoring range. These traditional localization methods are mostly based on the assumption of a single wave velocity model, however, the geological conditions in actual engineering have complexity. The cross arrangement of the underground chambers in the cross-connected areas of the underground works complicates the construction process. With the progress of the underground chamber excavation process, a cavity area is formed in the underground engineering. The wave velocity of the cavity area is greatly different from the wave velocity of surrounding rock mass, and the rock mass presents obvious anisotropic characteristics due to multi-layered or segmented distribution. Therefore, in the area, the wave velocity values of the shock waves have larger difference, and a single wave velocity model can cause larger error on a positioning result, so that the real space position of the rock mass fracture is difficult to accurately position. For the positioning of anisotropic rock mass, the traditional seismic source positioning method cannot accurately and effectively calculate the seismic source position of the fracture of the layered rock mass by taking the difference of the P-wave velocity values and the change of the P-wave ray paths of the elastic waves passing through different rock layers into consideration.
Disclosure of Invention
The invention aims to provide a method for locating wave velocity anisotropic rock mass fracture in a cross-connected area of underground engineering, which aims to solve the problem that the traditional method for locating the seismic source cannot accurately and effectively calculate the position of the seismic source of lamellar rock mass fracture.
In order to solve the technical problems, the aim of the invention is realized by the following technical scheme: the utility model provides a wave velocity anisotropic rock mass fracture positioning method of a cross communication area of underground engineering, which comprises the following steps:
s101, acquiring coordinate point positions of all sensors and coordinate point positions of actual vibration source points;
s102, traversing and calculating the running values of all the grid nodes from each sensor to a pre-built three-dimensional Cartesian coordinate grid model, and performing difference calculation on the running values of all the grid nodes from each sensor to each grid node in the three-dimensional Cartesian coordinate grid model to obtain a calculation result and establish a theoretical time difference database Th, wherein the three-dimensional Cartesian coordinate grid model is obtained by discretizing a pre-divided monitoring area;
s103, acquiring the acquisition initial time of the P wave, calculating the arrival time difference of each sensor based on the acquisition initial time of the P wave, and establishing an actual arrival time difference database Ac;
S104, using a first matrix T ij Representing the theoretical time difference database Th and using a second matrix A ij Representing an actual time difference database Ac, wherein the matrix A ij Element psi in (a) ij Expressed by the following formula:
ψ cd ＝p c -p d ＝(p in +Δp c )-(p in +Δp d )＝Δp c -Δp d (1≤c≤m,1≤d
≤m)
wherein m represents the number of all sensors in the monitored area, p c And p d Indicating the time values of arrival of the c-th sensor and the d-th sensor, p in Indicating the time of vibration, Δp c And Δp d The travel values of the seismic source to the c-th sensor and the d-th sensor are respectively represented;
s105, comparing the theoretical arrival time difference database Th with the actual arrival time difference database Ac according to the following formula to obtain an h function value:
wherein T is ij k Representing a theoretical arrival time difference matrix at a kth mesh node;
s106, after all the h function values are obtained, sequencing all the h function values according to the order from small to large;
s107, selecting the first q objective h function values, and labeling the objective grid nodes corresponding to the objective h function values, wherein the label format is (X q ,Y q ,Z q )；
S108, calculating and outputting the space coordinates (X) of the vibration source points according to the following formula o ,Y o ,Z o )：
Wherein, (X p ,Y p ,Z p ) And the coordinate value of any grid node in all target grid nodes is represented.
In addition, the technical problem to be solved by the invention is to provide a wave velocity anisotropic rock mass fracture positioning device for a cross communication area of underground engineering, which comprises the following components:
The acquisition unit is used for acquiring coordinate point positions of all the sensors and coordinate point positions of actual vibration source points;
the theoretical to time difference database building unit is used for traversing and calculating the running values of all grid nodes from each sensor to a pre-built three-dimensional Cartesian coordinate grid model, and performing difference calculation on the running values of all grid nodes from each sensor to each grid node in the three-dimensional Cartesian coordinate grid model to obtain a calculation result and building a theoretical to time difference database Th, wherein the three-dimensional Cartesian coordinate grid model is obtained by discretizing a pre-divided monitoring area;
the actual arrival time difference database establishing unit is used for acquiring the acquisition initial time of the P wave, calculating the arrival time difference of each sensor based on the acquisition initial time of the P wave and establishing an actual arrival time difference database Ac;
a representation unit for using the first matrix T ij Representing the theoretical time difference database Th and using a second matrix A ij Representing an actual time difference database Ac, wherein the matrix A ij Element psi in (a) ij Expressed by the following formula:
ψ cd ＝p c -p d ＝(p in +Δp c )-(p in +Δp d )＝Δp c -Δp d (1≤c≤m,1≤d
≤m)
wherein m represents the number of all sensors in the monitored area, p c And p d Indicating the time values of arrival of the c-th sensor and the d-th sensor, p in Indicating time of vibration，Δp c And Δp d The travel values of the seismic source to the c-th sensor and the d-th sensor are respectively represented;
the comparison unit is used for comparing the theoretical arrival time difference database Th with the actual arrival time difference database Ac according to the following formula to obtain the h function value:
wherein T is ij k Representing a theoretical arrival time difference matrix at a kth mesh node;
the sorting unit is used for sorting all the h function values according to the order from small to large after obtaining all the h function values;
a labeling unit for selecting the first q objective h function values and labeling the objective grid nodes corresponding to the objective h function values, wherein the labeling format is (X q ,Y q ,Z q )；
An output unit for calculating and outputting the spatial coordinates (X) of the source points according to the following formula o ,Y o ,Z o )：
Wherein, (X p ,Y p ,Z p ) And the coordinate value of any grid node in all target grid nodes is represented.
In addition, the embodiment of the invention also provides computer equipment, which comprises a memory, a processor and a computer program stored on the memory and capable of running on the processor, wherein the processor realizes the method for locating the wave velocity anisotropic rock mass fracture of the cross communication area of the underground engineering in the first aspect when executing the computer program.
In addition, an embodiment of the present invention further provides a computer readable storage medium, where the computer readable storage medium stores a computer program, where the computer program when executed by a processor causes the processor to execute the method for locating a wave velocity anisotropic rock mass fracture in a cross-connected area of underground engineering according to the first aspect.
The embodiment of the invention discloses a method for locating the fracture of a wave velocity anisotropic rock mass in a cross-connected area of underground engineering, which can accurately monitor the space position of a microseism event and can improve the microseism source locating precision of the rock mass in the area with anisotropic characteristics.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for the description of the embodiments will be briefly described below, and it is obvious that the drawings in the following description are some embodiments of the present invention, and other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a schematic flow chart of a method for locating wave velocity anisotropic rock mass fracture in a cross communication area of underground engineering, which is provided by the embodiment of the invention;
FIG. 2 is a diagram of a monitoring area sensor profile provided by an embodiment of the present invention;
FIG. 3 is a schematic diagram of an exemplary embodiment of an elongated domain distribution diagram;
FIG. 4 is a grid template vector diagram provided by an embodiment of the present invention;
FIG. 5 is a front view of the fracture positioning results provided by an embodiment of the present invention;
FIG. 6 is a top view of a fracture positioning result provided by an embodiment of the present invention;
FIG. 7 is a schematic block diagram of a device for locating wave velocity anisotropic rock mass fracture in a cross-connected region of an underground engineering provided by an embodiment of the invention;
fig. 8 is a schematic block diagram of a computer device according to an embodiment of the present invention.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and fully with reference to the accompanying drawings, in which it is evident that the embodiments described are some, but not all embodiments of the invention. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
It should be understood that the terms "comprises" and "comprising," when used in this specification and the appended claims, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof.
It is also to be understood that the terminology used in the description of the invention herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used in this specification and the appended claims, the singular forms "a," "an," and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise.
It should be further understood that the term "and/or" as used in the present specification and the appended claims refers to any and all possible combinations of one or more of the associated listed items, and includes such combinations.
Referring to fig. 1, fig. 1 is a schematic flow chart of a method for locating wave velocity anisotropic rock mass fracture in a cross-connected area of an underground engineering according to an embodiment of the invention;
as shown in fig. 1, the method includes steps S101 to S108.
S101, acquiring coordinate point positions of all sensors and coordinate point positions of actual vibration source points;
s102, traversing and calculating the running values of all the grid nodes from each sensor to a pre-built three-dimensional Cartesian coordinate grid model, and performing difference calculation on the running values of all the grid nodes from each sensor to each grid node in the three-dimensional Cartesian coordinate grid model to obtain a calculation result and establish a theoretical time difference database Th, wherein the three-dimensional Cartesian coordinate grid model is obtained by discretizing a pre-divided monitoring area;
It should be noted that, the monitoring area is defined according to the occurrence range of the microseismic event, and the monitoring area is made to be a cuboid, and is discretized into a three-dimensional cartesian coordinate grid model, wherein the length direction, the width direction and the height direction of the cuboid are respectively defined as an x-axis, a y-axis and a z-axis, and meanwhile, all the sensors in step S101 are all the sensors within the monitoring area.
Specifically, in this embodiment, the "the three-dimensional cartesian coordinate grid model is obtained by discretizing a pre-divided monitoring area" in step S102 includes:
s20, discretizing a pre-divided monitoring area into a three-dimensional Cartesian coordinate grid model;
s21, acquiring the microseismic longitudinal wave velocity of the monitoring area, and decomposing the microseismic longitudinal wave velocity into a first wave velocity v in the direction of a rock layer p And a second wave velocity v perpendicular to the formation direction v And characterizing the wave velocity of the cavity excavated in the monitoring area as a third wave velocity v k ；
S22, respectively establishing a first matrix V based on all the first wave velocity, the second wave velocity and the third wave velocity P Second matrix V V Third matrix V k Wherein the number of network grid points of the three-dimensional Cartesian coordinate grid model is determined based on the number of three-dimensional rows, columns and pages of the first matrix;
It should be added that since there is no other medium in the cavity, it can be considered that there is only air, and the wave velocity of the air medium is 340m/s, so that the wave velocity of all cavity areas is 340m/s, i.e. v k ＝340m/s。
S23, calculating the length Lx, the width Ly and the height Lz of the three-dimensional Cartesian coordinate grid model according to the following formula:
wherein ux, vy and wz characterize the grid sizes of the three-dimensional cartesian coordinate grid model in the x-axis, y-axis and z-axis directions, respectively.
In the present embodiment, pressIllumination V P The number of rows, columns and pages in three dimensions of the matrix determines the number of grid points of the monitoring area. The grid numbers u, V and w in the x-axis, y-axis and z-axis directions of the rectangular parallelepiped region are respectively equal to V P Row, column, page number of the matrix. Thus, each mesh node within the monitoring area is assigned a different wave velocity value.
In the monitoring area, according to the calculation condition of the travel time of the grid nodes, the grid nodes are divided into the following three grid node types: the points are calculated, the points are not calculated, and the points are to be fixed. In the three-dimensional Cartesian grid model, 6 neighbors (e.g., upper neighbor A 1 Next adjacent point A 2 Left adjacent point A 3 Right adjacent point A 4 Front adjacent point A 5 Post adjacent point A 6 ) The points for which the calculation of the time-lapse times of (a) are completed are called calculated points; the point where the running value is infinite is called the uncomputed point; all nodes except the calculated and uncalculated points are referred to as pending points. The elongated areas of points to be fixed are called elongated domains.
Before "walk time values of each sensor to all grid nodes in the three-dimensional cartesian coordinate grid model" in the step S102, the steps include:
s30, acquiring a coordinate point position of a current sensor, taking the current sensor as an initial point, and taking a grid node adjacent to the current sensor as an adjacent point;
in the present embodiment, a sensor S is used 1 (x 1 ,y 1 ,z 1 ) For example, where x 1 、y 1 And z 1 Respectively represent the sensor S 1 In the three-dimensional cartesian grid model, the coordinate positions in the X-axis, Y-axis and Z-axis are not necessarily adjacent to each other, and therefore, the sensors are not necessarily adjacent to each other, but may be common grid nodes.
S31, calculating adjacent point coordinates (i, j, k) corresponding to all adjacent points according to the following formula based on the coordinate point positions of the initial points to obtain corresponding long and narrow domains:
wherein h is x 、h y 、h z Steps in the x-axis, y-axis and z-axis directions, respectively;
s32, when the values of i, j and k in the adjacent point coordinates are all larger than 0 and the grid node type corresponding to the adjacent point is an uncomputed point, setting the step length according to the following formula:
h x ＝h y ＝h z ＝h
the last h in the formula is a preset value, and can be set according to the conditions of the size of the monitoring area, the accuracy requirement of the positioning result and the like.
S33, judging whether the initial point and the adjacent point are located on a plane formed by an x axis and a y axis in the three-dimensional Cartesian coordinate grid model, if so, calculating a transverse travel value of the adjacent point according to a formula (1), and if not, calculating a longitudinal travel value of the adjacent point according to a formula (2).
Wherein p (i, j, k) represents the running value of the grid node (i, j, k), and p (x, y, z) represents the running value of the grid node (x, y, z);
the meaning of h in the formula (2) is the same as that of h in the step S32, and the steps are all step sizes.
S34, acquiring a historical travel value t (i, j, k) of the adjacent point, and comparing the historical travel value with the calculated current travel value according to the following formula:
in this embodiment, the travel time must be a minimum travel time value, and equation (3) is a formula for finding the minimum travel time value, where the minimum travel time value can make the seismic wave propagation path obtained by the travel time calculation more conform to the actual propagation path, so that the positioning result is more accurate.
And S35, after the travel values of all the adjacent points in the long and narrow domain are obtained, selecting the adjacent point with the smallest current travel value in the long and narrow domain as an expansion point, and marking and storing the expansion point as a calculated point.
In the present embodiment, for the upper adjacent point A 1 Next adjacent point A 2 Left adjacent point A 3 Right adjacent point A 4 Front adjacent point A 5 Post adjacent point A 6 After the calculation of the running values of the adjacent points A is completed, the adjacent point A is assumed 1 If the running value of (a) is minimum, the upper adjacent point A 1 Marked as calculated point, its coordinate value is stored, and the calculation is performed as expansion point for the next partial walk, and the other points in the long and narrow domain are marked as A 2 、A 3 、A 4 、A 5 、A 6 。
In the first embodiment, the step S102 includes:
s40, acquiring adjacent points of the expansion points, judging the grid node type of the adjacent points, if the grid node type of the adjacent points is a calculated point, not updating the running value of the adjacent points, and continuously judging the next adjacent points;
s41, if the grid node type of the adjacent point is a to-be-determined point, judging whether the initial point and the adjacent point are located on a plane formed by an x axis and a y axis in the three-dimensional Cartesian coordinate grid model, if the initial point and the adjacent point are located on a plane formed by the x axis and the y axis in the three-dimensional Cartesian coordinate grid model, calculating a transverse travel value of the adjacent point according to a formula (4), and updating the travel value of the adjacent point according to a formula (3); if the two adjacent points are not located on the plane formed by the x axis and the y axis in the same three-dimensional Cartesian coordinate grid model, calculating the longitudinal travel values of the adjacent points according to a formula (2), and updating the travel values of the adjacent points according to a formula (3):
If the grid node type of the adjacent point is an uncomputed point, calculating the travel value of the adjacent point according to a formula (4) or a formula (2), and updating the grid node type of the adjacent point to be fixed;
s42, after updating the grid node types of all adjacent points of the expansion point, the new undetermined point forms a new long and narrow domain;
s43, acquiring the grid node with the smallest travel time value in the new long and narrow domain, performing adjacent point travel time calculation based on the grid node with the smallest travel time value, updating the node type of the adjacent point, and determining the new long and narrow domain until all the grid nodes are all calculated points.
After the calculation of the running values of the current sensor is completed, steps S30 to S43 are repeated, and the running values from the initial point to all grid nodes in the calculation area are calculated by using the remaining (m-1) sensors as the initial points.
The long and narrow domain refers to the long and narrow domain updated each time in the extension process of the long and narrow domain updated in S42.
In another embodiment, the step S102 includes:
s50, acquiring adjacent points of the initial point and adjacent points of the expansion point, judging the grid node type of the adjacent points, and if the grid node type of the adjacent points is a calculated point, continuing to judge the next adjacent point;
S51, if the grid node type of the adjacent point is a point to be fixed, calculating the travel time value of the adjacent point according to a formula (5), updating the travel time value of the adjacent point according to a formula (3), and if the grid node type of the adjacent point is an uncomputed point, updating the travel time value of the adjacent point according to a formula (5), and updating the grid node type of the adjacent point to be the point to be fixed:
wherein n represents O - With O + The order of P a,b,c Representing the running value, Q, of a grid node (a, b, c) a,b,c Representing the inverse speed of the grid node (a, b, c), wherein the inverse speed Q a,b,c Obtained by the following formula:
wherein V is a,b,c Representing an equivalent wave velocity at a grid node (a, b, c), the equivalent wave velocity at the grid node (a, b, c) being obtained by:
wherein, beta represents an included angle between the rock stratum and the propagation direction of the seismic wave;
it should be noted that, because the propagation direction of the seismic wave is not constant when it passes through the rock layer, the angle between the direction of the rock layer and the propagation direction of the seismic wave is denoted by β, and the angle β at each grid node is not necessarily the same.
Wherein in the x-axis direction, O - P and O + The first, second and third order forms of P are represented by the following formulas, respectively:
wherein the running value P at the grid node (a, b, c) a,b,c Calculated by the above formula (5).
The above formula is for explaining O in formula (5) - P term and O + The meaning of the term P is brought into equation (5) to obtain a variant of equation (5), and then Px, py, pz, pp, pq, pr, ps and the inverse speed are brought into the variant of equation (5), so that the final result can be solved, and the equation (5) is different due to the different ordersAlso different, but calculated based on equation (5).
Note that, the formula (5) calculates travel time by using the upper, lower, left, right, front, rear, adjacent points of the adjacent points (a, b, c) and the adjacent points in the diagonal direction (i.e., travel time by using nodes around the adjacent points), and meanwhile, the following explanation is given to the formula (5):
the running value is atIn the equivalent item, P a,b,c Is the running value of the grid node (a, b, c) to be calculated. N in formula (5) represents +.>The order of the equal terms is adjusted by +.>The order of the terms can change the accuracy and running speed of the algorithm, so O in the form of first order, second order and third order is exemplified in the application - P and O + And P, an appropriate order can be selected according to the requirement.
Taking n=1 as an example,that is, the first term in the first max function of equation (5), such that for O, respectively - P and O + P carry over (other O - P and O + P in P a-1,b,c Becomes P in the formula (6) and the formula (7) a+1,b,c 、P a,b,c-1 、P a,b-1,c 、P a,b+1,c Etc.), where P represents the travel value and the angle label is the node coordinates. (a-1, b, c) represents the adjacent point of (a, b, c) in the negative x direction, and the other points are the same. Δx represents the step length in the x direction, O - P and O + The sign of P and x, y, z represent positive and negative directions along the x, y, and z axes, respectively; in this way, only the running value of the point (a, b, c) to be solved is unknown, Q a,b,c Can also be according to->And (5) obtaining.
In this embodiment, in the step S51: "the running value of the adjacent point according to formula (5)" includes:
s60, if O - P and O + The order of P is 1, O in the x-axis direction is obtained as follows - P and O + P：
Wherein P is x 、P y 、P z Representing the minimum running values of the grid nodes (a, b, c) in the forward direction and the backward direction of the adjacent nodes in the x, y and z directions respectively, and calculating P according to the following formulas x 、P y 、P z ：
It should be noted that (a-1, b, c) in the above formula represents the neighboring points of the node (a, b, c) in the negative x-axis direction, because the grid node in the monitored area is divided into calculated points, uncomputed points, to-be-fixed points, the running values of the calculated points are known and remain unchanged, the running values of the to-be-fixed points are known but continue to be updated, and only the running values of the uncomputed points are unknown, so that the running times of the uncomputed points can be initially set to infinity, i.e. there is no point P in the above formula a-1,b,c Presence), i.e., the presence of the present application may also be expressed in the meaning of "known".
S61, based on the historical travel time information corresponding to the grid nodes on the diagonal positions of the grid nodes (a, b and c), selecting travel time values of the grid nodes (a, b and c) adjacent to 8 diagonal nodes, comparing the travel time of the corresponding nodes, and using P p 、P q 、P r 、P s Representing the minimum running values of grid nodes (a, b, c) in the forward and backward directions of two left-falling diagonals and the nodes adjacent to the two left-falling diagonals respectively, and representing P by the following formula p 、P q 、P r 、P s ：
In this embodiment, the min function is used to compare the travel times of the corresponding nodes, for example, in one of the diagonal directions, the travel values of the adjacent nodes of the nodes in the forward and backward directions are compared, and a smaller value is taken.
S62, P x 、P y 、P z 、P p 、P q 、P r 、P s Bringing into (5) the running value P of the grid node (a, b, c) is obtained a,b,c The method comprises the steps of carrying out a first treatment on the surface of the If the grid node type of the grid nodes (a, b and c) is the point to be fixed, adopting a formula (3) to update the running value; if the grid node (a, b, c) is an uncomputed point, updating the grid node type to be the point to be fixed.
In another embodiment, in the step S51: "the running value of the adjacent point according to formula (5)" includes:
S70, if O - P and O + The order of P is 2, and the running value calculation equation of the grid node is expressed by the following formula:
wherein the running value of the grid node (d, e) is P d,e 。Q d,e Representing the inverse velocity at the grid node (d, e), in two dimensions, the equivalent wave velocity at the grid node (d, e) is represented by:
V d,e ＝v p (d,e)·cosω+v v (d,e)·sinω
wherein ω is the angle between the formation and the direction of propagation of the seismic wave;
the grid nodes (d, e) are grid nodes in a two-dimensional space, and the grid nodes (a, b, c) are grid nodes in a three-dimensional space.
S71, the O - P and O + P is represented by the following formula:
note that Δx in the formula represents a step in the x-axis direction, Δy represents a step in the y-axis direction, and the formula in step S72 may be understood in correspondence with Px, py, and the like in S60 and S61. Pn (1) and Pn (2) correspond to travel time minimum expressions of Px and Py (x-axis and y-axis directions), and Pn (3) and Pn (4) correspond to travel time minimum expressions of Pp and Pq (diagonal directions). To facilitate the expression of S74 and S75 "t=1, 2" "" t=3, 4", the formal expression of Pn (t) is employed.
S72, use P n (1)、P n (2) Representing the minimum running values of grid nodes (d, e) adjacent to the grid nodes in the directions of the x axis and the y axis in the forward direction and the backward direction respectively, and using P n (3)、P n (4) Representing minimum running values of grid nodes (d, e) in forward and backward directions of a left-falling diagonal line and grid nodes adjacent to the left-falling diagonal line, respectively, and determining P by the following formula n (t)(t＝1,2,3,4)：
Condition 1: p (P) d-2,e <P d-1,e And P is d-1,e Presence;
condition 2: p (P) d+2,e <P d+1,e And P is d+1,e Presence;
condition 3: p (P) d,e-2 <P d,e-1 And P is d,e-1 Presence;
condition 4: p (P) d,e+2 <P d,e+1 And P is d,e+1 Presence;
condition 5: p (P) d-2,e+2 <P d-1,e+1 And P is d-1,e+1 Presence;
condition 6: p (P) d+2,e-2 <P d+1,e-1 And P is d+1,e-1 Presence;
condition 7: p (P) d-2,e-2 <P d-1,e-1 And P is d-1,e-1 Presence;
condition 8: p (P) d+2,e+2 <P d+1,e+1 And P is d+1,e+1 Presence;
condition 1 and condition 2 are simultaneously established
Condition 3 and condition 4 are simultaneously satisfied
Condition 5 and condition 6 are simultaneously satisfied
Condition 7 and condition 8 are simultaneously satisfied
S73, the formula (8) is arranged into a unitary quadratic equation, coefficients are D, E, F respectively, root calculation is carried out on the formula (8) by adopting the following formula based on D noteq 0, and a calculation result is represented by Pp (t) (t=1, 2,3, 4):
s74, based on the input conditional instruction, determining whether there is travel time information provided by the diagonal position, and if there is no travel time information provided by the diagonal position, calculating the coefficient D, E, F according to the following formula:
calculating the result Pp at the time of the calculation 1 (t)＝max(Pp(1),Pp(2))；
S75, if the travel time information provided by the diagonal position is provided, the coefficient D, E, F is calculated according to the following formula:
Pp 2 (t) =max (Pp (3), pp (4)), if Pp 2 (t) is present, the final calculation result pp=min (Pp 1 (t),Pp 2 (t))。
It should be emphasized that the conventional method generally uses the longitudinal wave velocity to calculate the running value directly, and the method of steps S30-S35 of the present applicationThe method considers v p And v v Namely, the wave velocity along the formation direction and the wave velocity perpendicular to the formation direction, and it is emphasized that S60 to S75 do not directly use the longitudinal wave velocity of the conventional method, but use the equivalent wave velocity; the method in steps S60-S62 considers the 8 diagonal running values, and the method in steps S70-75 considers the 4 diagonal running values, so that the calculated running values are more accurate, and the final calculation result is improved in accuracy.
S103, acquiring the acquisition initial time of the P wave, calculating the arrival time difference of each sensor based on the acquisition initial time of the P wave, and establishing an actual arrival time difference database Ac;
s104, using a first matrix T ij Representing the theoretical time difference database Th and using a second matrix A ij Representing an actual time difference database Ac, wherein the matrix A ij Element psi in (a) ij Expressed by the following formula:
ψ cd ＝p c -p d ＝(p in +Δp c )-(p in +Δp d )＝Δp c -Δp d (1≤c≤m,1≤d
≤m)
wherein m represents the number of all sensors in the monitored area, p c And p d Indicating the time values of arrival of the c-th sensor and the d-th sensor, p in Indicating the time of vibration, Δp c And Δp d The travel values of the seismic source to the c-th sensor and the d-th sensor are respectively represented;
s105, comparing the theoretical arrival time difference database Th with the actual arrival time difference database Ac according to the following formula to obtain an h function value:
wherein T is ij k Representing a theoretical arrival time difference matrix at a kth mesh node;
s106, after all the h function values are obtained, sequencing all the h function values according to the order from small to large;
in this embodiment, the calculated h function values are sorted from small to large, the first q h function values are taken as the smallest, and the corresponding nodes are labeled (X 1 ,Y 1 ,Z 1 ),(X 2 ,Y 2 ,Z 2 ),(X 3 ,Y 3 ,Z 3 )…(X q ,Y q ,Z q )。
S107, selecting the first q objective h function values, and labeling the objective grid nodes corresponding to the objective h function values, wherein the label format is (X q ,Y q ,Z q )；
S108, calculating and outputting the space coordinates (X) of the vibration source points according to the following formula o ,Y o ,Z o )：
Wherein, (X p ,Y p ,Z p ) And the coordinate value of any grid node in all target grid nodes is represented.
The term (X) o ,Y o ,Z o ) The method is the optimal result of the focus positioning obtained by the method.
In the embodiment, the method is adopted to perform microseismic source positioning on a certain underground chamber group with anisotropic characteristics, surrounding rocks at the position of the underground chamber group are mainly limestone, dolomite and quartz rocks, and class A and class B breccia are distributed locally. According to the geological survey data of the area, the wave speed range of the rock mass in the field area is 4500 m/s-6000 m/s. The wave velocity of the excavated part of the underground chamber group is 340m/s.
And defining a cuboid monitoring area as shown in the figure, and positioning a microseismic source in the anisotropic rock mass area. The length, width and height directions of the cuboid are respectively an x axis, a y axis and a z axis. Establishing a three-dimensional square grid coordinate model, wherein the grid sizes of the monitoring area in the x-axis direction, the y-axis direction and the z-axis direction are ux=vy=wz=1m, and the length, the width and the height of the three-dimensional square grid coordinate model are lx=200 m, ly=174 m and lz=96 m respectively. Corner O (1081.2, -65.4,754.7) is the origin of coordinates, and the following grid nodes are generated in the monitoring domain:
the grid nodes are positively coded along the x-axis, the y-axis and the z-axis, the coding of the original point position O (1081.2-65.4,754.7) is (1, 1), the x-axis, the y-axis and the z-axis of the cuboid (namely the three-dimensional square grid coordinate model) are respectively regarded as rows, columns and layers, the coding of the nodes positioned on the a-th row, the B-th column and the C-th layer is (a, B, C), the corresponding node space sitting marks are (A, B, C), wherein the value range of a is 1-201, the value range of B is 1-175, and the value range of C is 1-97.
In the application, 18 microseismic monitoring sensors are arranged in a monitoring area and are numbered S i (i=1 to 18), and the coordinates are denoted as (x) i ,y i ,z i ) Three-dimensional space coordinates of each sensor are measured and recorded, when the codes of the sensors are determined, the sensors are classified into adjacent nearest grid nodes according to the coordinate positions of the sensors, and the code positions (a i ,b i ,c i ). The spatial coordinates and codes of each sensor are shown in the following table:
TABLE 1
The distribution of the sensors in the underground chamber group is shown in figure 2.
In a certain time, a microseismic event occurs in the monitored area due to engineering disturbance, and the position coordinates of the source points are assumed to be (x) o ,y o ,z o ) Sensor S 1 ～S 18 When the first arrival of the micro-seismic wave is received, the first arrival time of the ith sensor is recorded as p i (1.ltoreq.i.ltoreq.18), in the present embodiment, S 1 ～S 18 The received P wave first arrival time is shown in the following table:
TABLE 2
Through geological exploration, rock bodies in a monitoring area are distributed in a layered manner, all rock layers can be divided into 4 wave velocity categories according to rock lithology differences, and the wave velocity value V of the rock layers is recorded m (1.ltoreq.m.ltoreq.4). The wave velocity value of the cavity part artificially excavated in the monitoring area is recorded as V k ，V k =340 m/s. Since the anisotropic rock mass has a large velocity difference between the rock formation and the direction perpendicular to the rock formation, V is m Decomposing into wave velocities V parallel to the formation direction p And velocity V perpendicular to the formation direction v And reconstructing the velocity model with the equivalent wave velocity v. V is represented by the following formula:
where β is the angle between the formation and the direction of propagation of the seismic wave.
The velocity values of the wave velocity layers are shown in the following table:
TABLE 3 Table 3
Wave velocity V to be parallel to the formation direction p And wave velocity V perpendicular to the formation direction v And (3) endowing grid nodes of the corresponding area, and calculating the equivalent speed at each grid node according to the included angle between the rock stratum and the seismic wave propagation direction. For the hole area, the speed value of the grid node where the hole area is located is assigned to 340m/s. Based on the wave velocity values at each node, a velocity matrix v of size 201×175×97 is established.
Each sensor (x i ,y i ,z i ) The running value Δp to each grid node (a, b, c) in the monitoring area i (a,b,c)。
In the calculation process, at a certain moment, the distribution of the long and narrow domains and the three types of nodes in the monitoring area is shown in fig. 3.
The nodes in the long and narrow domain are ordered, the position with the minimum travel time is determined, and the adjacent point travel time calculation is performed based on the position. In the three-dimensional area, a grid template vector diagram for calculating node travel time by using the method of the invention is shown in fig. 4.
And repeating the extension process of the long narrow domain until the types of the updated points are all calculated points, namely, when the number of nodes in the long narrow domain is 0, and calculating the seismic wave travel values of all grid nodes in the monitoring area.
Travel time Δp for each sensor to each grid node i (a, b, c) working out differences two by two, building a theoretical to time difference database Th, and using a first matrix T ij And (3) representing.
Based on the time-of-day data p received by the 18 sensors i Calculating the arrival time difference between the sensors, establishing an actual arrival time difference database Ac, and using a second matrix A ij And (3) representing. Because of the large amount of data calculated, the matrix T is not listed here ij Matrix A ij The values of the elements in (a).
Setting a time residual function h, comparing the databases Th and Ac, and calculating a time difference function value h at each node a,b,c . H is represented by a,b,c ：
Where n represents the number of sensors and i and j are the labels of the sensors. T (T) ij a,b,c Representing a theoretical to moveout matrix at the node (a, b, c).
H is calculated to be a,b,c Sequencing from small to large, taking the first 15 h a,b,c A grid node corresponding to the function value of (a).
TABLE 4 Table 4
The arithmetic mean of these 15 grid nodes is calculated as the final positioning result of each source point using the following formula, the source coordinates are calculated using (x o ,y o ,z o ) And (3) representing.
Wherein (X) p ,Y p ,Z p ) Coordinate values representing any one of the 15 nodes, q being the number of nodes, q=15.
Calculated, the source coordinates are (x o ,y o ,z o )＝(1146.3,54.5,779.6)。
And calculating the position coordinates of other microseismic events in the monitored area by adopting the same method. A total of 68 rock mass breaking events occurred during a certain period of time. The results of the source localization using the method of the present invention are shown in fig. 5 and 6.
In this embodiment, the actual spatial coordinates of the source in the field are (X o ,Y o ,Z o ) = (1142.9,49.4,777.5). True source location (X) o ,Y o ,Z o ) And calculating the source coordinates (x o ,y o ,z o ) The three-dimensional space distance between the two is regarded as the positioning error of the microseismic positioning. The positioning error is calculated by the following formula:
the positioning error was calculated to be 6.5m.
Comparative example
The comparative example aims at discussing the difference between the seismic source positioning method and the single wave velocity model positioning method in the aspect of microseism positioning results. The single wave velocity model treats the rock mass medium in the monitoring area as homogeneous and single. According to the wave velocity of the rock mass in the area, the wave velocity value is uniformly set to 5300m/s. In performing the travel time calculation, v=5300 m/s is assigned to each node within the area.
The sensor layout is the same as in example 1. And carrying out objective function solution by adopting an optimization algorithm Newton iteration method, wherein the objective function of the example is shown as follows:
wherein p is a And p b Indicating the time delta P when the a-th and b-th sensors receive the first arrival of the P wave a And Δp b Representing the travel time of the P-wave from the source location to the a-th and b-th sensors.
Establishing a single speed model, and adopting a traditional optimization algorithm Newton method to calculate an objective function minimum value asThe spatial coordinates corresponding to the position are (1153.3,37.8,791.4). The positioning error at this time is also calculated by the formula (4), that is:
As can be seen from the positioning errors of the embodiment and the comparative example, the method provided by the invention has the advantages that the positioning result of the seismic source is more accurate and is closer to the actual seismic source position. Therefore, the method can accurately monitor the space position of the microseism event, and can improve the microseism source positioning accuracy of the regional rock mass with anisotropic characteristics.
The embodiment of the invention also provides a device for positioning the wave velocity anisotropic rock mass fracture of the cross communication area of the underground engineering, which is used for executing any embodiment of the method for positioning the wave velocity anisotropic rock mass fracture of the cross communication area of the underground engineering. Specifically, referring to fig. 7, fig. 7 is a schematic block diagram of an apparatus for locating a fracture of an anisotropic rock mass in a wave velocity in a cross-connected area of an underground engineering according to an embodiment of the present invention.
As shown in fig. 7, the wave velocity anisotropic rock mass fracture positioning device 500 for the cross-connected region of the underground engineering comprises:
an obtaining unit 501, configured to obtain coordinate point positions of all sensors and coordinate point positions of actual seismic source points;
the theoretical-time difference database establishing unit 502 is configured to traverse and calculate running values from each sensor to all grid nodes in a pre-built three-dimensional cartesian coordinate grid model, and perform difference calculation on the running values from each sensor to each grid node in the three-dimensional cartesian coordinate grid model to obtain a calculation result and establish a theoretical-time difference database Th, where the three-dimensional cartesian coordinate grid model is obtained by discretizing a pre-divided monitoring area;
An actual arrival time difference database establishing unit 503, configured to obtain an initial acquisition time of a P-wave, calculate an arrival time difference of each sensor based on the initial acquisition time of the P-wave, and establish an actual arrival time difference database Ac;
a representation unit 504 for using a first matrix T ij Representing the theoretical time difference database Th and using a second matrix A ij Representing an actual time difference database Ac, wherein the matrix A ij Element psi in (a) ij Expressed by the following formula:
ψ cd ＝p c -p d ＝(p in +Δp c )-(p in +Δp d )＝Δp c -Δp d (1≤c≤m,1≤d≤m)
wherein m represents the number of all sensors in the monitored area, p c And p d Indicating the time values of arrival of the c-th sensor and the d-th sensor, p in Indicating the time of vibration, Δp c And Δp d The travel values of the seismic source to the c-th sensor and the d-th sensor are respectively represented;
the comparison unit 505 is configured to compare the theoretical arrival time difference database Th with the actual arrival time difference database Ac according to the following formula, to obtain an h function value:
wherein T is ij k Representing a theoretical arrival time difference matrix at a kth mesh node;
a sorting unit 506, configured to sort all the h function values according to a descending order after obtaining all the h function values;
a labeling unit 507 for selecting the first q objective h function values and labeling the objective grid nodes corresponding to the objective h function values, wherein the labeling format is (X q ,Y q ,Z q )；
An output unit 508 for calculating and outputting the spatial coordinates (X) of the source points according to the following formula o ,Y o ,Z o )：
Wherein, (X p ,Y p ,Z p ) And the coordinate value of any grid node in all target grid nodes is represented.
It will be clear to those skilled in the art that, for convenience and brevity of description, specific working procedures of the apparatus and units described above may refer to corresponding procedures in the foregoing method embodiments, which are not described herein again.
The above-described wave velocity anisotropic rock mass fracture positioning apparatus for an underground engineering cross-connect region may be implemented in the form of a computer program which may be run on a computer device as shown in fig. 8.
Referring to fig. 8, fig. 8 is a schematic block diagram of a computer device according to an embodiment of the present invention. The computer device 1100 is a server, and the server may be a stand-alone server or a server cluster formed by a plurality of servers.
With reference to FIG. 8, the computer device 1100 includes a processor 1102, memory, and a network interface 1105 connected through a system bus 1101, wherein the memory may include a non-volatile storage medium 1103 and an internal memory 1104.
The non-volatile storage medium 1103 may store an operating system 11031 and computer programs 11032. The computer program 11032, when executed, causes the processor 1102 to perform a method for locating a wave velocity anisotropic rock mass fracture in a cross-connect area of an underground works.
The processor 1102 is operable to provide computing and control capabilities to support the operation of the overall computer device 1100.
The internal memory 1104 provides an environment for the execution of a computer program 11032 in the non-volatile storage medium 1103, which computer program 11032, when executed by the processor 1102, causes the processor 1102 to perform a method for locating wave velocity anisotropic rock mass fracture in a cross-over connected area of an underground project.
The network interface 1105 is used for network communication such as providing transmission of data information, etc. It will be appreciated by those skilled in the art that the architecture shown in fig. 8 is merely a block diagram of some of the architecture relevant to the present inventive arrangements and is not limiting of the computer device 1100 to which the present inventive arrangements may be implemented, and that a particular computer device 1100 may include more or fewer components than shown, or may combine some of the components, or have a different arrangement of components.
Those skilled in the art will appreciate that the embodiment of the computer device shown in fig. 8 is not limiting of the specific construction of the computer device, and in other embodiments, the computer device may include more or less components than those shown, or certain components may be combined, or a different arrangement of components. For example, in some embodiments, the computer device may include only a memory and a processor, and in such embodiments, the structure and function of the memory and the processor are consistent with the embodiment shown in fig. 8, and will not be described again.
It should be appreciated that in embodiments of the invention, the processor 1102 may be a central processing unit (Central Processing Unit, CPU), the processor 1102 may also be other general purpose processors, digital signal processors (Digital Signal Processor, DSPs), application specific integrated circuits (Application Specific Integrated Circuit, ASICs), off-the-shelf programmable gate arrays (Field-Programmable Gate Array, FPGAs) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, or the like. Wherein the general purpose processor may be a microprocessor or the processor may be any conventional processor or the like.
In another embodiment of the invention, a computer-readable storage medium is provided. The computer readable storage medium may be a non-volatile computer readable storage medium. The computer readable storage medium stores a computer program, wherein the computer program when executed by a processor implements the wave velocity anisotropic rock mass fracture positioning method of the cross-connected region of the underground engineering according to the embodiment of the invention.
The storage medium is a physical, non-transitory storage medium, and may be, for example, a U-disk, a removable hard disk, a Read-Only Memory (ROM), a magnetic disk, or an optical disk.
It will be clearly understood by those skilled in the art that, for convenience and brevity of description, specific working procedures of the apparatus, device and unit described above may refer to corresponding procedures in the foregoing method embodiments, which are not repeated herein.
While the invention has been described with reference to certain preferred embodiments, it will be understood by those skilled in the art that various changes and substitutions of equivalents may be made and equivalents will be apparent to those skilled in the art without departing from the scope of the invention. Therefore, the protection scope of the invention is subject to the protection scope of the claims.

Claims (10)

1. A wave velocity anisotropic rock mass fracture positioning method for a cross communication area of an underground engineering is characterized by comprising the following steps:
s101, acquiring coordinate point positions of all sensors and coordinate point positions of actual vibration source points;
s102, traversing and calculating the running values of all the grid nodes from each sensor to a pre-built three-dimensional Cartesian coordinate grid model, and performing difference calculation on the running values of all the grid nodes from each sensor to each grid node in the three-dimensional Cartesian coordinate grid model to obtain a calculation result and establish a theoretical time difference database Th, wherein the three-dimensional Cartesian coordinate grid model is obtained by discretizing a pre-divided monitoring area;
S103, acquiring the acquisition initial time of the P wave, calculating the arrival time difference of each sensor based on the acquisition initial time of the P wave, and establishing an actual arrival time difference database Ac;
s104, using a first matrix T ij Representing the theoretical time difference database Th and using a second matrix A ij Representing an actual time difference database Ac, wherein the matrix A ij Element psi in (a) ij Expressed by the following formula:
ψ cd ＝p c -p d ＝(p in +Δp c )-(p in +Δp d )＝Δp c -Δp d (1≤c≤m,1≤d
≤m)
wherein m represents the number of all sensors in the monitored area, p c And p d Indicating the time values of arrival of the c-th sensor and the d-th sensor, p in Indicating the time of vibration, Δp c And Δp d The travel values of the seismic source to the c-th sensor and the d-th sensor are respectively represented;
s105, comparing the theoretical arrival time difference database Th with the actual arrival time difference database Ac according to the following formula to obtain an h function value:
wherein T is ij k Representing a theoretical arrival time difference matrix at a kth mesh node;
s106, after all the h function values are obtained, sequencing all the h function values according to the order from small to large;
s107, select the first qThe objective h function values are marked, and the objective grid nodes corresponding to the objective h function values are marked, wherein the marking format is (X q ,Y q ,Z q )；
S108, calculating and outputting the space coordinates (X) of the vibration source points according to the following formula o ,Y o ,Z o )：
Wherein, (X p ,Y p ,Z p ) And the coordinate value of any grid node in all target grid nodes is represented.
2. The method for locating wave velocity anisotropic rock mass fracture in cross-over communication area of underground engineering according to claim 1, wherein the three-dimensional cartesian coordinate grid model is obtained by discretizing a pre-divided monitoring area in step S102, comprising:
s20, discretizing a pre-divided monitoring area into a three-dimensional Cartesian coordinate grid model;
s21, acquiring the microseismic longitudinal wave velocity of the monitoring area, and decomposing the microseismic longitudinal wave velocity into a first wave velocity v in the direction of a rock layer p And a second wave velocity v perpendicular to the formation direction v And characterizing the wave velocity of the cavity excavated in the monitoring area as a third wave velocity v k ；
S22, respectively establishing a first matrix V based on all the first wave velocity, the second wave velocity and the third wave velocity P Second matrix V V Third matrix V k Wherein the number of network grid points of the three-dimensional Cartesian coordinate grid model is determined based on the number of three-dimensional rows, columns and pages of the first matrix;
s23, calculating the length Lx, the width Ly and the height Lz of the three-dimensional Cartesian coordinate grid model according to the following formula:
wherein ux, vy and wz characterize the grid sizes of the three-dimensional cartesian coordinate grid model in the x-axis, y-axis and z-axis directions, respectively.
3. The method for locating a wave velocity anisotropic rock mass fracture in a cross-connected region of an underground engineering according to claim 2, wherein before "the traversing and calculating the travel values of each sensor to all grid nodes in the three-dimensional cartesian coordinate grid model" in step S102, the method comprises:
s30, acquiring a coordinate point position of a current sensor, taking the current sensor as an initial point, and taking a grid node adjacent to the current sensor as an adjacent point;
s31, calculating adjacent point coordinates (i, j, k) corresponding to all adjacent points according to the following formula based on the coordinate point positions of the initial points to obtain corresponding long and narrow domains:
wherein h is x 、h y 、h z Steps in the x-axis, y-axis and z-axis directions, respectively;
s32, when the values of i, j and k in the adjacent point coordinates are all larger than 0 and the grid node type corresponding to the adjacent point is an uncomputed point, setting the step length according to the following formula:
h x ＝h y ＝h z ＝h
s33, judging whether the initial point and the adjacent point are located on a plane formed by an x axis and a y axis in the three-dimensional Cartesian coordinate grid model, if so, calculating a transverse travel value of the adjacent point according to a formula (1), and if not, calculating a longitudinal travel value of the adjacent point according to a formula (2).
Wherein p (i, j, k) represents the running value of the grid node (i, j, k), and p (x, y, z) represents the running value of the grid node (x, y, z);
s34, acquiring a historical travel value t (i, j, k) of the adjacent point, and comparing the historical travel value with the calculated current travel value according to the following formula:
and S35, after the travel values of all the adjacent points in the long and narrow domain are obtained, selecting the adjacent point with the smallest current travel value in the long and narrow domain as an expansion point, and marking and storing the expansion point as a calculated point.
4. The method for locating wave velocity anisotropic rock mass fracture in cross-over communication area of underground works according to claim 3, wherein said step S102 comprises:
s40, acquiring adjacent points of the expansion points, judging the grid node type of the adjacent points, if the grid node type of the adjacent points is a calculated point, not updating the running value of the adjacent points, and continuously judging the next adjacent points;
s41, if the grid node type of the adjacent point is a to-be-determined point, judging whether the initial point and the adjacent point are located on a plane formed by an x axis and a y axis in the three-dimensional Cartesian coordinate grid model, if the initial point and the adjacent point are located on a plane formed by the x axis and the y axis in the three-dimensional Cartesian coordinate grid model, calculating a transverse travel value of the adjacent point according to a formula (4), and updating the travel value of the adjacent point according to a formula (3); if the two adjacent points are not located on the plane formed by the x axis and the y axis in the same three-dimensional Cartesian coordinate grid model, calculating the longitudinal travel values of the adjacent points according to a formula (2), and updating the travel values of the adjacent points according to a formula (3):
If the grid node type of the adjacent point is an uncomputed point, calculating the travel value of the adjacent point according to a formula (4) or a formula (2), and updating the grid node type of the adjacent point to be fixed;
s42, after updating the grid node types of all adjacent points of the expansion point, the new undetermined point forms a new long and narrow domain;
s43, acquiring the grid node with the smallest travel time value in the new long and narrow domain, performing adjacent point travel time calculation based on the grid node with the smallest travel time value, updating the node type of the adjacent point, and determining the new long and narrow domain until all the grid nodes are all calculated points.
5. The method for locating wave velocity anisotropic rock mass fracture in cross-over communication area of underground works according to claim 4, wherein said step S102 comprises:
s50, acquiring adjacent points of the initial point and adjacent points of the expansion point, judging the grid node type of the adjacent points, and if the grid node type of the adjacent points is a calculated point, continuing to judge the next adjacent point;
s51, if the grid node type of the adjacent point is a point to be fixed, calculating the travel time value of the adjacent point according to a formula (5), updating the travel time value of the adjacent point according to a formula (3), and if the grid node type of the adjacent point is an uncomputed point, updating the travel time value of the adjacent point according to a formula (5), and updating the grid node type of the adjacent point to be the point to be fixed:
Wherein n represents O - With O + The order of P a,b,c Representing the running value, Q, of a grid node (a, b, c) a,b,c Representing the inverse speed of the grid node (a, b, c), wherein the inverse speed Q a,b,c Obtained by the following formula:
wherein V is a,b,c Representing an equivalent wave velocity at a grid node (a, b, c), the equivalent wave velocity at the grid node (a, b, c) being obtained by:
wherein, beta represents an included angle between the rock stratum and the propagation direction of the seismic wave;
wherein in the x-axis direction, O - P and O + The first, second and third order forms of P are represented by the following formulas, respectively:
wherein the running value P at the grid node (a, b, c) a,b,c Calculated by the above formula (5).
6. The method for locating wave velocity anisotropic rock mass fracture in cross-over communication area of underground works according to claim 5, wherein in step S51: "the running value of the adjacent point according to formula (5)" includes:
s60, if O - P and O + The order of P is 1, O in the x-axis direction is obtained as follows - P and O + P：
Wherein P is x 、P y 、P z Representing the minimum running values of the grid nodes (a, b, c) in the forward direction and the backward direction of the adjacent nodes in the x, y and z directions respectively, and calculating P according to the following formulas x 、P y 、P z ：
S61, based on the historical travel time information corresponding to the grid nodes on the diagonal positions of the grid nodes (a, b and c), selecting travel time values of the grid nodes (a, b and c) adjacent to 8 diagonal nodes, comparing the travel time of the corresponding nodes, and using P p 、P q 、P r 、P s Representing the minimum running values of grid nodes (a, b, c) in the forward and backward directions of two left-falling diagonals and the nodes adjacent to the two left-falling diagonals respectively, and representing P by the following formula p 、P q 、P r 、P s ：
S62, P x 、P y 、P z 、P p 、P q 、P r 、P s Bringing into (5) the running value P of the grid node (a, b, c) is obtained a,b,c The method comprises the steps of carrying out a first treatment on the surface of the If the grid node type of the grid nodes (a, b and c) is the point to be fixed, adopting a formula (3) to update the running value; if the grid node (a, b, c) is an uncomputed point, updating the grid node type to be the point to be fixed.
7. The method for locating wave velocity anisotropic rock mass fracture in cross-over communication area of underground works according to any one of claims 5 to 6, wherein in step S51: "the running value of the adjacent point according to formula (5)" includes:
s70, if O - P and O + The order of P is 2, and the running value calculation equation of the grid node is expressed by the following formula:
wherein the running value of the grid node (d, e) is P d,e 。Q d,e Representing the inverse velocity at the grid node (d, e), in two dimensions, the equivalent wave velocity at the grid node (d, e) is represented by:
V d,e ＝v p (d,e)·cosω+v v (d,e)·sinω
wherein ω is the angle between the formation and the direction of propagation of the seismic wave;
s71, the O - P and O + P is represented by the following formula:
s72, use P n (1)、P n (2) Representing the minimum running values of grid nodes (d, e) adjacent to the grid nodes in the directions of the x axis and the y axis in the forward direction and the backward direction respectively, and using P n (3)、P n (4) Representing minimum running values of grid nodes (d, e) in forward and backward directions of a left-falling diagonal line and grid nodes adjacent to the left-falling diagonal line, respectively, and determining P by the following formula n (t)(t＝1,2,3,4)：
Condition 1: p (P) d-2,e <P d-1,e And P is d-1,e Presence;
condition 2: p (P) d+2,e <P d+1,e And P is d+1,e Presence;
condition 3: p (P) d,e-2 <P d,e-1 And P is d,e-1 Presence;
condition 4: p (P) d,e+2 <P d,e+1 And P is d,e+1 Presence;
condition 5: p (P) d-2,e+2 <P d-1,e+1 And P is d-1,e+1 Presence;
condition 6: p (P) d+2,e-2 <P d+1,e-1 And P is d+1,e-1 Presence;
condition 7: p (P) d-2,e-2 <P d-1,e-1 And P is d-1,e-1 Presence;
condition 8: p (P) d+2,e+2 <P d+1,e+1 And P is d+1,e+1 Presence;
s73, the formula (8) is arranged into a unitary quadratic equation, coefficients are D, E, F respectively, root calculation is carried out on the formula (8) by adopting the following formula based on D noteq 0, and a calculation result is represented by Pp (t) (t=1, 2,3, 4):
s74, based on the input conditional instruction, judging whether the travel time information provided by the diagonal position exists, if the travel time information provided by the diagonal position does not exist, calculating the coefficient D, E, F according to the following formula to obtain a travel time calculation result Pp 1 (t)＝max(Pp(1),Pp(2))：
S75, if the travel time information provided by the diagonal position is provided, the coefficient D, E, F is calculated according to the following formula:
Pp 2 (t) =max (Pp (3), pp (4)), if Pp 2 (t) is present, the final calculation result pp=min (Pp 1 (t),Pp 2 (t))。
8. The utility model provides a wave velocity anisotropic rock mass fracture positioner of underground works cross communication region which characterized in that includes:
The acquisition unit is used for acquiring coordinate point positions of all the sensors and coordinate point positions of actual vibration source points;
the theoretical to time difference database building unit is used for traversing and calculating the running values of all grid nodes from each sensor to a pre-built three-dimensional Cartesian coordinate grid model, and performing difference calculation on the running values of all grid nodes from each sensor to each grid node in the three-dimensional Cartesian coordinate grid model to obtain a calculation result and building a theoretical to time difference database Th, wherein the three-dimensional Cartesian coordinate grid model is obtained by discretizing a pre-divided monitoring area;
the actual arrival time difference database establishing unit is used for acquiring the acquisition initial time of the P wave, calculating the arrival time difference of each sensor based on the acquisition initial time of the P wave and establishing an actual arrival time difference database Ac;
a representation unit for using the first matrix T ij Representing the theoretical time difference database Th and using a second matrix A ij Representing an actual time difference database Ac, wherein the matrix A ij Element psi in (a) ij Expressed by the following formula:
ψ cd ＝p c -p d ＝(p in +Δp c )-(p in +Δp d )＝Δp c -Δp d (1≤c≤m,1≤d
≤m)
wherein m represents the number of all sensors in the monitored area, p c And p d Indicating the time values of arrival of the c-th sensor and the d-th sensor, p in Indicating the time of vibration, Δp c And Δp d The travel values of the seismic source to the c-th sensor and the d-th sensor are respectively represented;
the comparison unit is used for comparing the theoretical arrival time difference database Th with the actual arrival time difference database Ac according to the following formula to obtain the h function value:
wherein T is ij k Representing a theoretical arrival time difference matrix at a kth mesh node;
the sorting unit is used for sorting all the h function values according to the order from small to large after obtaining all the h function values;
a labeling unit for selecting the first q objective h function values and labeling the objective grid nodes corresponding to the objective h function values, wherein the labeling format is (X q ,Y q ,Z q )；
An output unit for calculating and outputting the spatial coordinates (X) of the source points according to the following formula o ,Y o ,Z o )：
Wherein, (X p ,Y p ,Z p ) And the coordinate value of any grid node in all target grid nodes is represented.
9. A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, wherein execution of the computer program by the processor implements the method for locating a wave velocity anisotropic rock mass fracture in a cross-connected region of an underground works as claimed in any one of claims 1 to 7.
10. A computer readable storage medium, characterized in that the computer readable storage medium stores a computer program which, when executed by a processor, causes the processor to perform the method of locating a wave velocity anisotropic rock mass fracture in a cross-connected area of an underground works as claimed in any one of claims 1 to 7.
CN202310423762.1A 2023-04-19 2023-04-19 Wave velocity anisotropic rock mass fracture positioning method for cross communication area of underground engineering Pending CN116466401A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310423762.1A CN116466401A (en) 2023-04-19 2023-04-19 Wave velocity anisotropic rock mass fracture positioning method for cross communication area of underground engineering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310423762.1A CN116466401A (en) 2023-04-19 2023-04-19 Wave velocity anisotropic rock mass fracture positioning method for cross communication area of underground engineering

Publications (1)

Publication Number Publication Date
CN116466401A true CN116466401A (en) 2023-07-21

Family

ID=87174865

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310423762.1A Pending CN116466401A (en) 2023-04-19 2023-04-19 Wave velocity anisotropic rock mass fracture positioning method for cross communication area of underground engineering

Country Status (1)

Country Link
CN (1) CN116466401A (en)

Similar Documents

Publication Publication Date Title
Lu Systematic identification of polyhedral rock blocks with arbitrary joints and faults
Shi et al. Automatic generation of road network map from massive GPS, vehicle trajectories
US10795053B2 (en) Systems and methods of multi-scale meshing for geologic time modeling
BRPI0820174B1 (en) METHOD FOR FORMING A GEOLOGICAL MODEL OF AN EAR REGION AND METHOD FOR DRILLING A HOLE IN AN EARTH REGION
CA2680021A1 (en) Model-based time-preserving tomography
CA2897304A1 (en) Method to invert for fault activity and tectonic stress
Cameron et al. Seismic velocity estimation from time migration
Zhang Advances in three-dimensional block cutting analysis and its applications
CN105301639A (en) Speed field inversion method and device based on VSP double-weight travel time tomography
KR20080023946A (en) 3-d gravity inversion method of underground cavities using euler deconvolution and 3-d imaging method using it
CN109459787B (en) coal mine underground structure imaging method and system based on seismic channel wave full-waveform inversion
CN116774292B (en) Seismic wave travel time determining method, system, electronic equipment and storage medium
CN110221342A (en) Seismic source location method, apparatus and storage medium based on three-dimensional velocity structure
CN112035939A (en) Rock-soil body parameter random field modeling method for double-side-wall pilot tunnel
CN114943178A (en) Three-dimensional geological model modeling method and device and computer equipment
CN103698810A (en) Hybrid network minimum travel time ray tracing tomography method
CN107817524B (en) The method and apparatus of three-dimensional seismic tomography
CN117331119A (en) Rapid earthquake wave travel time calculation method for tunnel detection
CN108828669A (en) A kind of two-dimensional intersection survey line static corrections processing method, apparatus and system
Fu et al. Spatial topology identification of three-dimensional complex block system of rock masses
CN116466401A (en) Wave velocity anisotropic rock mass fracture positioning method for cross communication area of underground engineering
CN106950601A (en) Static correcting method and device
Ban et al. Critical gates identification for fault-tolerant design in math circuits
Potter et al. Ordered line integral methods for solving the eikonal equation
US6505147B1 (en) Method for process simulation

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