Acoustic emission / micro-seismic positioning method in mined area-containing and non-uniform wave velocity scene

Abstract

The invention relates to an acoustic emission / micro-seismic positioning method under a mined-out area-containing and non-uniform wave velocity scene. The method comprises the following steps: step 1, carrying out block division on a micro-seismic monitoring area, calculating time required for vibration wave propagation between vertexes in each block, constructing a travel time network model, and adding a sensor into the travel time network model; 2, determining coordinates of an initial iteration point of a seismic source by adopting a least square algorithm; 3, when the coordinates of the initial iteration point of the seismic source are located in the goaf, correcting the coordinates of the initial iteration point of the seismic source by using a BFS algorithm; step 4, using a Dijkstra ray tracing algorithm to calculate a propagation duration and a propagation path from the seismic source to the sensor; and step 5, based on the time difference between the propagation duration from the seismic source to each sensor and the actually received signal, correcting the coordinates of the seismic source through a Geiger iterative algorithm, and when the coordinates of the iteration point of the seismic source are located in the goaf, correcting by using a BFS algorithm.

Description

Technical Field

[0001] The invention belongs to the technical field of acoustic emission / microseismic positioning methods in rock engineering, and relates to an acoustic emission / microseismic positioning method in a scene containing a void area and uneven wave velocity. Background Art

[0002] The rapid economic growth of my country has greatly stimulated the development of resources and energy, most of which involve rock engineering. The construction process of rock engineering will cause the formation of fracture surfaces inside the rock body and the propagation of stress waves. When cracks or damage occur inside the rock due to stress concentration, elastic energy will be released, generating tiny vibration waves. These vibration waves propagate outward at a certain speed. Multiple acoustic emission / microseismic sensors deployed in the monitoring area receive vibration signals and record the arrival time and waveform characteristics of the vibration waves. The spatial coordinates and occurrence time of the earthquake source can be calculated using the positioning algorithm.

[0003] At present, microseismic monitoring technology is widely used in the location and early warning of rock damage areas and potential rock bursts and rock bursts. Classical source location methods such as Geiger, simplex method, LocSAT method, and double difference location method have good positioning effects in uniform wave velocity scenarios, but do not consider the impact of voids and wave velocity inhomogeneity on the positioning effect. In view of the inhomogeneity of the geological body itself and the spatiotemporal evolution of the mining area, in order to achieve better monitoring and early warning effects, the scope of application of the acoustic emission / microseismic location algorithm has gradually evolved from the idealized scenario of uniform wave velocity to the complex scenario of non-uniform wave velocity and voids. In order to consider the impact of voids on the wave propagation path and arrival time, scholars began to combine the shortest travel time model with the positioning algorithm. For example, the A* algorithm and Dijkstra algorithm were combined with grid search to realize the source positioning in the uniform wave velocity field containing voids. However, the non-uniformity of the wave velocity field was not considered. In addition, the search time of the grid search algorithm in three-dimensional space is proportional to the cube of the grid size, while the positioning accuracy is inversely proportional to the grid size. It is impossible to achieve both positioning time and accuracy, which limits the engineering application effect of the algorithm. Summary of the invention

[0004] In order to solve the above technical problems, the purpose of the present invention is to provide an acoustic emission / microseismic positioning method in a scene with void areas and uneven wave velocity.

[0005] The present invention provides an acoustic emission / microseismic positioning method in a scene containing a void area and uneven wave velocity, comprising:

[0006] Step 1: Divide the microseismic monitoring area into blocks, calculate the time required for the seismic wave to propagate between the vertices in each block, build a travel time network model, and add the sensor to the travel time network model;

[0007] Step 2: Use the least squares algorithm to determine the coordinates of the initial iteration point of the earthquake source;

[0008] Step 3: When the coordinates of the initial iteration point of the earthquake source are located in the empty area, the coordinates of the initial iteration point of the earthquake source are corrected using the BFS algorithm;

[0009] Step 4: Use the Dijkstra ray tracing algorithm to calculate the propagation time and propagation path from the source to the sensor;

[0010] Step 5: Based on the time difference between the propagation time from the earthquake source to each sensor and the time when the signal is actually received, the earthquake source coordinates are corrected by the Geiger iterative algorithm. When the coordinates of the iterative point of the earthquake source are located in the empty area, the BFS algorithm is used for correction.

[0011] An acoustic emission / microseismic positioning method in a scene with a void area and uneven wave velocity has the following beneficial effects:

[0012] The positioning method of the present invention overcomes the difficult problem of low positioning accuracy under the influence of the uneven wave velocity field caused by the formation of voids due to rock excavation and lithology differences. The influence of the existence of voids and the difference in wave velocity spatial distribution on the positioning result is considered through the three-dimensional wave velocity field travel time network model and ray tracing. At the same time, ray tracing is combined with an iterative algorithm to overcome the problem that the positioning time consumption and accuracy cannot be achieved in the process of combining ray tracing with a grid search algorithm. At the same time, the introduced void escape algorithm can avoid the location of the source in the void. The positioning algorithm proposed by the present invention significantly reduces the source positioning error caused by the existence of voids and uneven wave velocity, and reduces the positioning time consumption, which is conducive to engineering applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 It is a flow chart of an acoustic emission / microseismic positioning method in a scene containing a void area and uneven wave velocity according to the present invention;

[0014] Figure 2 is a schematic diagram of a given microseismic monitoring area in an embodiment of the present invention;

[0015] Figure 3 It is a schematic diagram of cutting the microseismic monitoring area into a plurality of identical cubic blocks in an embodiment of the present invention;

[0016] Figure 4 A schematic diagram of establishing undirected edges between vertices of a cubic block in an embodiment of the present invention;

[0017] Figure 5 ; A schematic diagram of establishing undirected edges between sensors, seismic sources and vertices of the blocks to which they belong in an embodiment of the present invention; DETAILED DESCRIPTION

[0018] In response to the problems existing in the current acoustic emission / microseismic positioning algorithm, this paper proposes a positioning algorithm suitable for the scenario of non-uniform wave velocity field containing void areas based on the assumption that wave propagation has the shortest travel time, combined with the travel time network model considering void areas and non-uniform wave velocity field, the iterative algorithm with fast convergence, and the void area escape algorithm.

[0019] like Figure 1 As shown, the acoustic emission / microseismic positioning method of the present invention in a scene containing a void area and uneven wave velocity includes:

[0020] Step 1: Divide the microseismic monitoring area into blocks, calculate the time required for the seismic wave to propagate between the vertices in each block, build a travel time network model, and add sensors to the travel time network model. Specifically:

[0021] Step 1.1: Cut the microseismic monitoring area into multiple identical cubic blocks, determine the vertex coordinates of each block, and for each block, establish undirected edges between its vertices.

[0022] In specific implementation, given Figure 2 The microseismic monitoring area is shown in the figure. Different gray areas represent different lithologies. Figure 3 The microseismic monitoring area is divided into multiple identical cubic blocks. Figure 4 Schematic diagram of establishing undirected edges between vertices of a cubic block.

[0023] Step 1.2: Measure the wave velocities of different lithologies in the monitoring area. Set the wave velocity in the empty area to the speed of sound. Set the wave velocity of the undirected edge to the wave velocity of the lithology at the center of the cube. The weight of the undirected edge is the wave propagation time, i.e., travel time.

[0024] Step 1.3: Set up multiple sensors in the microseismic monitoring area and establish undirected edges between sensors, seismic sources and the vertices of the blocks to which they belong, such as Figure 5 As shown, in this embodiment, eight sensors S1-S8 are set in the microseismic monitoring area, and Source represents the earthquake source. Figure 5 Undirected edges are established between sensors, sources and each vertex of the block to which they belong.

[0025] Step 2: Use the least squares algorithm to determine the coordinates of the initial iteration point of the earthquake source, specifically:

[0026] Step 2.1: For N sensors receiving the shock waves released by the same earthquake source, establish the following shock wave propagation distance equation:

[0027] (x i -x) 2 +(y i -y) 2 +(z i -z)2 =v 2 (t i -t) 2

[0028] Where: i = 1, 2, ..., N, N> 4 represents the number of sensors; x i ,y i , z i represents the coordinates of the ith sensor; x, y, z represent the coordinates of the initial iteration point of the source to be solved; v represents the velocity of the shock wave, and the global average wave velocity is used here; t i represents the time when the shock wave reaches the i-th sensor; t represents the moment when the shock wave is generated.

[0029] Step 2.2: Select an equation and subtract it from the other equations to get:

[0030] d j x++b j y+c j z+d j t=e j

[0031] Where: j = 1, 2, ..., N-1, a j , b j , c j , d j , e j are the coefficients of each item, let:

[0032]

[0033] Then the above formula is rewritten as AX=B, and the least squares method is used to solve it:

[0034] x * =(A T A) -1 B

[0035] Among them, x * The coordinates of the initial iteration point representing the earthquake source.

[0036] Step 3: When the coordinates of the initial iteration point of the earthquake source are located in the empty area, the BFS algorithm is used to correct the coordinates of the initial iteration point of the earthquake source, specifically:

[0037] Step 3.1: Define G = (V, E) to represent a weighted graph, where V represents the set of vertices of all blocks and E represents the set of all undirected edges; assume that the coordinates of the initial iteration point s of the earthquake source are V s , given a queue Q storing the nodes to be traversed, a set M storing the traversed nodes, v air Indicates the wave velocity in the empty area.

[0038] Step 3.2: Initially, add the initial iteration point s of the source to the queue Q and the set M.

[0039] Step 3.3: Take the head q from the queue Q and determine whether the head q satisfies v q >v air , v q represents the wave velocity of the undirected edge where the team head q is located. If it satisfies, let s = q, and the corrected source coordinates are V q , exit step 3; otherwise, proceed to step 3.4.

[0040] Step 3.4: Traverse all adjacent nodes of q that are not included in the set M Will Add queue Q and set M and repeat step 3.3.

[0041] Step 4: Use the Dijkstra ray tracing algorithm to calculate the propagation time and propagation path from the source to the sensor, specifically:

[0042] Step 4.1: In the weighted graph G = (V, E), define d(i, j) to represent the shortest travel time from node i to node j, and w(i, j) to represent the weight of the undirected edge from node i to node j, i.e., the travel time; ensure that w(i, j) ≥ 0. If there is no edge between node i and node j, then w(i, j) is ∞.

[0043] Define the relaxation operation. For d(s,v) and node u, if d(s,v)>d(s,u)+w(u,v), then d(s,v)=d(s,u)+w(u,v) means that d(s,v) is relaxed by node u. Define prec[v] as the previous node of the shortest travel time path from the source s to the node v, that is, the predecessor node.

[0044] Step 4.2: Initialization: for all nodes v∈V, set the initial travel time d(s,v)=∞, indicating that the shortest travel time path from the source to each node is unknown; set the source d(s,s)=0, and mark all nodes as unvisited.

[0045] Step 4.3: Traverse the nodes, select the node u with the smallest d(s,u) from the unvisited nodes each time, mark u as visited, and perform relaxation operations on all its adjacent nodes v:

[0046] d(s,v)=min(d(s,v),d(s,u)+w(u,v))

[0047] prev[v]=u

[0048] Step 4.4: Repeat step 4.3 until all nodes are visited and go to step 5.

[0049] Step 5: Based on the propagation time from the source to each sensor and the time difference between the actual received signal, the source coordinates are corrected by the Geiger iterative algorithm. When the coordinates of the iterative point of the source are located in the empty area, BFS is used for correction, specifically:

[0050] Step 5.1: Assume that the shock wave propagates in a straight line, satisfying the equation:

[0051]

[0052] Where: x*, y*, z* are the corrected earthquake source coordinates in step 3; t* is the earthquake time; x i ,y i , z i represents the coordinates of the i-th sensor; t i Represents the shock wave reaching the i-th sensor s i time; i represents the average wave velocity along the shortest travel time path from the source to the sensor after correction in step 3; P(s,s i ) is the corrected source to the i-th sensor s in step 3 i The physical path length of all edges of d(s,s i ) is the corrected source to the i-th sensor s in step 3 i The travel time.

[0053] Step 5.2: Based on the propagation time from the source to each sensor and the time difference between the actual received signal, determine whether the accuracy requirement is met or the maximum number of iterations is reached. If so, the algorithm ends; otherwise, proceed to step 5.3.

[0054] Step 5.3: Perform a first-order Taylor expansion on the equation in step 5.1:

[0055]

[0056] Where: t oi is the arrival time of the shock wave detected by the i-th sensor; t ci is the time when the shock wave reaches the i-th sensor calculated from the source coordinates.

[0057] Step 5.4: Solve by least square method, let:

[0058]

[0059] Write the first-order Taylor expansion formula in matrix form:

[0060] CΔX=D

[0061] Solving by least squares method, we get:

[0062] ΔX=(C T C) -1 C T D

[0063] Step 5.5: Get the iteratively updated earthquake source coordinates:

[0064] X′=(X+ΔX)

[0065]

[0066] Where x′, y′, z′ are the new source coordinates, and t′ is the time when the new source signal is generated.

[0067] Step 5.6: Use the BFS algorithm to correct the iterative coordinates of the source located in the empty area, the same as step 3.

[0068] Step 5.7: Add the new source coordinates to the travel time network.

[0069] Step 5.8: Repeat step 4 until the accuracy requirement is met or the maximum number of iterations is reached.

[0070] The above description is only a preferred embodiment of the present invention and is not intended to limit the concept of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. An acoustic emission / microseismic positioning method in a scene with voids and uneven wave velocity, characterized in that: include: Step 1: Divide the microseismic monitoring area into blocks, calculate the time required for the seismic wave to propagate between the vertices in each block, build a travel time network model, and add the sensor to the travel time network model; Step 2: Use the least squares algorithm to determine the coordinates of the initial iteration point of the earthquake source; Step 3: When the coordinates of the initial iteration point of the earthquake source are located in the empty area, the coordinates of the initial iteration point of the earthquake source are corrected using the BFS algorithm; Step 4: Use the Dijkstra ray tracing algorithm to calculate the propagation time and propagation path from the source to the sensor; Step 5: Based on the time difference between the propagation time from the earthquake source to each sensor and the time when the signal is actually received, the earthquake source coordinates are corrected by the Geiger iterative algorithm. When the coordinates of the iterative point of the earthquake source are located in the empty area, the BFS algorithm is used for correction.

2. The acoustic emission / microseismic positioning method in a scenario with voids and uneven wave velocity as claimed in claim 1, characterized in that: The step 1 is specifically as follows: Step 1.1: Cut the microseismic monitoring area into multiple identical cubic blocks, determine the vertex coordinates of each block, and for each block, establish undirected edges between its vertices; Step 1.2: Measure the wave velocity of different lithologies in the monitoring area. Set the wave velocity in the empty area to the speed of sound. Set the wave velocity of the undirected edge to the wave velocity of the lithology at the center of the cube. The weight of the undirected edge is the wave propagation time, i.e., travel time. Step 1.3: Set up multiple sensors in the microseismic monitoring area and establish undirected edges between the sensors, seismic sources and the vertices of the blocks to which they belong.

3. The acoustic emission / microseismic positioning method in a scene with void area and uneven wave velocity as claimed in claim 2, characterized in that: The step 2 is specifically as follows: Step 2.1: For N sensors receiving the shock waves released by the same earthquake source, establish the following shock wave propagation distance equation: (x i -x) 2 +(y i -y) 2 +(z i -z) 2 =v 2 (t i -t) 2 Where: i = 1, 2, ..., N, N > 4 represents the number of sensors; x i ,y i , z i represents the coordinates of the ith sensor; x, y, z represent the coordinates of the initial iteration point of the source to be solved; v represents the velocity of the shock wave, and the global average wave velocity is used here; t i Represents the time when the shock wave reaches the i-th sensor; t represents the time when the shock wave is generated; Step 2.2: Select an equation and subtract it from the other equations to get: a j x+b j y+c j z+d j t=e j Where: j = 1, 2, ..., N-1, a j , b j , c j , d j , e j are the coefficients of each item, let: Then the above formula is rewritten as AX=B, and the least squares method is used to solve it: x * =(A T A) -1 B Among them, x * The coordinates of the initial iteration point representing the earthquake source.

4. The acoustic emission / microseismic positioning method in a scenario with voids and uneven wave velocity as claimed in claim 2, characterized in that: The step 3 is specifically as follows: Step 3.1: Define G = (V, E) to represent a weighted graph, where V represents the set of vertices of all blocks and E represents the set of all undirected edges; assume that the coordinates of the initial iteration point s of the earthquake source are V s , given a queue Q storing the nodes to be traversed, a set M storing the traversed nodes, v air It represents the wave speed in the empty zone; Step 3.2: Initially, add the initial iteration point s of the earthquake source to the queue Q and the set M; Step 3.3: Take the head q from the queue Q and determine whether the head q satisfies v q >v air , v q represents the wave velocity of the undirected edge where the team head q is located. If it satisfies, let s = q, and the corrected source coordinates are V q , exit step 3; Otherwise, proceed to step 3.4; Step 3.4: Traverse all adjacent nodes of q that are not included in the set M Will Add queue Q and set M and repeat step 3.

3.

5. The acoustic emission / microseismic positioning method in a scenario with voids and uneven wave velocity as claimed in claim 4, characterized in that: The step 4 is specifically as follows: Step 4.1: In the weighted graph G = (V, E), define d(i, j) to represent the shortest travel time from node i to node j, and w(i, j) to represent the weight of the undirected edge from node i to node j, i.e., the travel time; ensure that w(i, j) ≥ 0. If there is no edge between node i and node j, then w(i, j) is ∞; Define the relaxation operation. For d(s, v) and node u, if d(s, v)>d(s, u)+w(u, v), then d(s, v)=d(s, u)+w(u, v) means d(s, v) is relaxed by node u. Define prec[v] as the previous node of the shortest travel time path from the earthquake source s to the node v, that is, the predecessor node. Step 4.2: Initialization: for all nodes v∈V, set the initial travel time d(s, v) = ∞, indicating that the shortest travel time path from the source to each node is unknown; set the source d(s, s) = 0, and mark all nodes as unvisited; Step 4.3: Traverse the nodes, select the node u with the smallest d(s,u) from the unvisited nodes each time, mark u as visited, and perform relaxation operations on all its adjacent nodes v: d(s,v)=min(d(s,v),d(s,u)+w(u,v)) prev[v]=u Step 4.4: Repeat step 4.3 until all nodes are visited and go to step 5.

6. The acoustic emission / microseismic positioning method in a scenario with voids and uneven wave velocity as claimed in claim 1, characterized in that: The step 5 is specifically as follows: Step 5.1: Assume that the shock wave propagates in a straight line, satisfying the equation: Where: x*, y*, z* are the corrected earthquake source coordinates in step 3; t* is the earthquake time; x i ,y i , z i represents the coordinates of the i-th sensor; t i Represents the shock wave reaching the i-th sensor s i time; i represents the average wave velocity along the shortest travel time path from the source to the sensor after correction in step 3; P(s,s i ) is the corrected source to the i-th sensor s in step 3 i The physical path length of all edges of d(s,s i ) is the corrected source to the i-th sensor s in step 3 i Travel time; Step 5.2: Based on the propagation time from the source to each sensor and the time difference between the actual received signal, determine whether the accuracy requirement is met or the maximum number of iterations is reached. If so, the algorithm ends; otherwise, proceed to step 5.3; Step 5.3: Perform a first-order Taylor expansion on the equation in step 5.1: Where: t oi is the arrival time of the shock wave detected by the i-th sensor; t ci is the time when the shock wave reaches the i-th sensor calculated from the earthquake source coordinates; Step 5.4: Solve by least square method, let: Write the first-order Taylor expansion formula in matrix form: CΔX=D Solving by least squares method, we get: ΔX=(C T C) -1 C T D Step 5.5: Get the iteratively updated earthquake source coordinates: X'=(X+ΔX) In the formula, x', y', z' are the new source coordinates, and t' is the time when the new source signal is generated; Step 5.6: Use BFS to correct the iterative coordinates of the source located in the empty area, the same as step 3; Step 5.7: Add the new source coordinates to the travel time network; Step 5.8: Repeat step 4 until the accuracy requirement is met or the maximum number of iterations is reached.

Images