An acoustic emission / microseismic positioning method in a scene with voids and uneven wave velocity

Abstract

The present invention provides an acoustic emission / microseismic positioning method in a scenario containing a void area and uneven wave velocity, comprising: step 1: dividing the microseismic monitoring area into blocks, calculating the time required for vibration wave propagation between vertices within each block, constructing a travel time network model, and adding sensors to the travel time network model; step 2: using a least squares algorithm to determine the coordinates of an initial iteration point of the earthquake source; step 3: when the coordinates of the initial iteration point of the earthquake source are located in a void area, using a BFS algorithm to correct the coordinates of the initial iteration point of the earthquake source; step 4: using a Dijkstra ray tracing algorithm to calculate the propagation time and propagation path from the earthquake source to the sensor; step 5: based on the propagation time from the earthquake source to each sensor and the time difference between the actual received signal, correcting the earthquake source coordinates using a Geiger iterative algorithm, and when the coordinates of the iteration point of the earthquake source are located in a void area, correction is performed using the BFS algorithm.

Description

Technical Field

[0001] The present 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 voids and uneven wave velocity. Background Art

[0002] my country's rapid economic growth has greatly stimulated the development of resources and energy, much of which involves rock engineering. The construction process of rock engineering projects causes the formation of fracture surfaces within the rock mass, accompanied by the propagation of stress waves. When stress concentration forms cracks or damage occurs within the rock, elastic energy is released, generating tiny vibration waves that propagate outward at a constant speed. Multiple acoustic emission / microseismic sensors deployed in the monitoring area receive the vibration signals and record the arrival time and waveform characteristics of the vibration waves. Using positioning algorithms, the spatial coordinates and occurrence time of the earthquake source can be calculated.

[0003] Currently, microseismic monitoring technology is widely used to locate and provide early warning for rock damage areas and potential rockbursts and rockbursts. Classic source location methods such as the Geiger method, simplex method, LocSAT method, and double-difference method offer good positioning results in scenarios with uniform velocity, but they fail to consider the impact of voids and velocity heterogeneity on positioning. Given the inherent heterogeneity of geological bodies and the spatiotemporal evolution of mining areas, to achieve better monitoring and early warning results, the scope of application of acoustic emission / microseismic location algorithms has gradually evolved from idealized scenarios with uniform velocity to complex scenarios with heterogeneous 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 a uniform wave velocity field containing voids. However, the non-uniformity of the wave velocity field was not considered. Moreover, 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. The positioning time and accuracy cannot be achieved at the same time, 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 voids 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 vertices within each block, build a travel-time network model, and add sensors 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 difference between the propagation time from the source to each sensor and the actual time of the received signal, the source coordinates are corrected using the Geiger iterative algorithm. When the coordinates of the iterative point of the source are in the empty area, the BFS algorithm is used for correction.

[0011] An acoustic emission / microseismic positioning method in a scene with voids 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 due to the formation of voids due to rock excavation and the influence of uneven wave velocity fields caused by differences in lithology. Through the three-dimensional wave velocity field travel time network model and ray tracing, the influence of the existence of voids and the differences in the spatial distribution of wave velocities on the positioning results is taken into account. 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 simultaneously when ray tracing is used in conjunction with a grid search algorithm. At the same time, the void escape algorithm introduced can prevent the earthquake source from being located in a void. The positioning algorithm proposed by the present invention significantly reduces the earthquake source positioning error caused by the existence of voids and uneven wave velocity, and reduces the positioning time consumption, which is beneficial to engineering applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 This is a flow chart of an acoustic emission / microseismic positioning method in a scenario with voids 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 Schematic diagram of cutting the microseismic monitoring area into multiple identical cubic blocks in an embodiment of the present invention;

[0016] Figure 4 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 respective vertices of the blocks to which they belong in an embodiment of the present invention; DETAILED DESCRIPTION

[0018] To address the problems existing in the current acoustic emission / microseismic positioning algorithms, this paper proposes a positioning algorithm suitable for scenarios with non-uniform wave velocity fields containing voids, based on the assumption that wave propagation has the shortest travel time, combined with a travel time network model that considers voids and non-uniform wave velocity fields, a rapidly convergent iterative algorithm, and a void escape algorithm.

[0019] like Figure 1 As shown, the present invention provides an acoustic emission / microseismic positioning method in a scenario with a void area and uneven wave velocity, comprising:

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

[0021] Step 1.1: Divide 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 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 cube.

[0023] 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., the travel time.

[0024] Step 1.3: Set up multiple sensors in the microseismic monitoring area and establish undirected edges between the sensors, the 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 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 i-th 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 For each coefficient, let:

[0032]

[0033] The above formula is rewritten as AX=B and solved using the least squares method:

[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 It represents the wave velocity of the undirected edge where the team head q is located. If it satisfies the requirement, 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 earthquake source s to the node v, that is, the predecessor node.

[0044] Step 4.2: Initialize. 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 difference between the propagation time from the source to each sensor and the actual time of receiving the signal, the source coordinates are corrected using the Geiger iterative algorithm. When the coordinates of the iteration point of the source are in the empty area, the 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 after 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; v 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 time difference between the propagation time from the source to each sensor and the actual time of the 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 the problem by the 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 earthquake 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 principles of the present invention should be included in the scope of protection 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 vertices within each block, build a travel-time network model, and add sensors 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 difference between the propagation time from the source to each sensor and the actual time of the received signal, the source coordinates are corrected using the Geiger iterative algorithm. When the coordinates of the iterative point of the source are 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 according to claim 1, characterized in that: The step 1 is specifically as follows: Step 1.1: Divide 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., the travel time. Step 1.3: Place multiple sensors in the microseismic monitoring area and establish undirected edges between the sensors, the seismic sources, and the vertices of the blocks to which they belong.

3. The acoustic emission / microseismic positioning method in a void area and non-uniform wave velocity scenario according to claim 2, characterized in that: The step 2 is specifically as follows: Step 2.1: For N sensors receiving 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 i-th 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; 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 For each coefficient, let: The above formula is rewritten as AX=B and solved using the least squares method: 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 void area and non-uniform wave velocity scenario according to 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 represents the wave velocity in the null region; 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 It represents the wave velocity of the undirected edge where the team head q is located. If it satisfies the requirement, 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 according to 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 nodes i and 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), indicating that d(s, v) is relaxed by node u. Define prec[v] as the previous node of the shortest travel time path from source s to node v, i.e., the predecessor node. Step 4.2: Initialize. 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 void area and non-uniform wave velocity scenario according to 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 after 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; v 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 source coordinates; Step 5.4: Solve the problem by the 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) Where 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