A microseismic positioning method for self-adapting wave velocity field optimization in complex rock mass

Abstract

This invention relates to a microseismic location method with adaptive wave velocity field optimization in complex rock masses, belonging to the field of microseismic monitoring. The method includes: delineating a monitoring area according to the monitoring task, deploying microseismic sensors within the monitoring area, and representing the coordinates of any three-dimensional grid node within the monitoring area as (…). x, y, z ), will 3D mesh nodes ( x, y, z ) imparts a wave velocity field V0 to the rock mass x, y, z Starting from each sensor, the arrival times of the first arrivals at each three-dimensional grid node are calculated. A residual function between the theoretical arrival time and the actual arrival time is established. The first n nodes corresponding to the smallest values ​​among all the calculated residual function results are selected, and then the coordinates of the fracture source point are calculated. This invention achieves the location of rock fracture in areas with cavities and complex rock mass wave velocity by calculating the residual values ​​between the theoretical arrival time and the actual acquired microseismic wave arrival times and searching for nodes that satisfy the minimum residual.

Description

Technical Field

[0001] This invention relates to the field of microseismic monitoring, and more particularly to a microseismic location method with adaptive wave velocity field optimization in complex rock masses. Background Technology

[0002] Human engineering activities, such as rock excavation and blasting, can cause the rock mass to fracture and release elastic waves into the outside world. Microseismic monitoring technology can collect this stress wave information, locate the fracture source, and ultimately conduct stability analysis and early warning of the rock mass. This technology has been widely used in many fields such as water conservancy and mining.

[0003] Microseismic location technology utilizes microseismic sensors deployed in the engineering monitoring area to collect stress waves generated by the fracturing of the engineering rock mass. Combined with on-site engineering geological conditions, such as rock mass wave velocity measurement results, the spatial coordinates of the microseismic events of the rock mass fracturing are inverted. Currently, microseismic location technology methods include linear solution location algorithms, such as the Geiger algorithm, relative location method, and optimization function joint inversion method, as well as microseismic location algorithms that adopt nonlinear solution approaches, such as Newton's method and simplex method. These methods simplify the engineering rock mass into rock masses with the same wave velocity and achieve accurate location of the rock mass as the source of the seismic fracture through an iterative process and optimization location. However, while existing microseismic location methods can identify the location of rupture sources in a relatively short time when applied to relatively small study areas or areas with homogeneous lithology, they are difficult to accurately measure when applied to engineering geological environments with complex rock mass wave velocities. Therefore, treating the surrounding rock as a rock mass with a uniform wave velocity leads to excessive errors in the analysis results. In addition, as engineering activities progress, human engineering activities will create large cavities in the surrounding rock environment. If a uniform equivalent wave velocity model is used, it will greatly interfere with the microseismic location results. Furthermore, the rock mass wave velocity information used in the location work is usually based on the analysis results of a limited number of measuring points. Using a simple rock mass wave velocity model will also face the problem of large location error. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the prior art and provide a microseismic location method with adaptive wave velocity field optimization in complex rock masses, thus solving the deficiencies of the prior art.

[0005] The objective of this invention is achieved through the following technical solution: a microseismic location method with adaptive wave velocity field optimization in complex rock masses, the microseismic location method comprising:

[0006] The monitoring area is delineated according to the monitoring task, microseismic sensors are deployed in the monitoring area, and the coordinates of any three-dimensional grid node in the monitoring area are represented as (x,y,z). The three-dimensional grid node (x,y,z) is assigned the rock mass wave velocity field V0(x,y,z).

[0007] Starting from each sensor, the initial arrival time to each 3D grid node is calculated to obtain the initial arrival time of the rock mass fracture elastic wave from the i-th microseismic sensor to the 3D grid node (x,y,z). When a microseismic event is transmitted in the monitoring area, the arrival time of the rock fracture stress wave collected on-site by the i-th microseismic sensor is: ;

[0008] Establish the residual function between theoretical arrival time and actual arrival time, select the first n nodes corresponding to the smallest value among all the results of the calculated residual function, and then calculate the coordinates of the rupture source point.

[0009] The calculation of the initial arrival time to each 3D mesh node, starting from each sensor, specifically includes the following:

[0010] A1. Equation for a three-dimensional wave velocity field , Let x0, y0, z0 represent the initial arrival / departure time when reaching the 3D mesh node (x, y, z). Then, the initial arrival / departure time for the starting point (x0, y0, z0) is... ;

[0011] A2. The initial arrival and arrival times of the three-dimensional mesh nodes (x, y, z) in the non-uniform wave velocity field are equivalently deformed as follows: , This represents the initial arrival and arrival times of the three-dimensional mesh nodes (x, y, z) obtained using the equivalent uniform wave velocity field model. It is represented as the initial arrival travel time optimization correction coefficient located at the 3D mesh node (x,y,z);

[0012] A3. For the previously known rock mass wave velocity field V0(x,y,z), the optimized rock mass wave velocity field is expressed as follows based on posterior or optimization information: V0(x,y,z) represents the previously known rock mass wave velocity field of the three-dimensional grid node (x,y,z). It is represented as the wave velocity field optimization correction coefficient located at the three-dimensional grid node (x,y,z);

[0013] A4. According to the formula and For the equation of function The transformation yields the equation used to fuse physical information and machine learning algorithms. ;

[0014] A5. Construct the input features of the machine learning network model, construct the intermediate and output layers of the network model, and calculate the results based on the loss function;

[0015] Repeat steps A1-A5. Stop repeating when the maximum number of repetitions is reached or the final error function result is less than the threshold. Substitute into the formula In the process, when the initial arrival time to each 3D mesh node is obtained, the obtained... Substitute into the formula The optimized wave velocity values ​​of each three-dimensional mesh node are obtained.

[0016] The A5 step specifically includes the following:

[0017] Input features for constructing a machine learning network model: The input to the machine learning network model consists of travel time and wave speed adjustment coefficients for each 3D grid node. Initially, let... ;

[0018] Constructing the intermediate layers of the network model: The network model has a total of p layers, with q neurons in each layer. The weights of the i-th neuron in the j-th layer are represented by W as linear and constant terms, respectively. pq With b pq ;

[0019] Construct the model output layer: The output result is represented as The loss function of the network model is expressed as: ;

[0020] Based on the loss function calculation results, the loss function is calculated on W using the steepest gradient descent method. pq , b pq , and The partial derivatives are calculated, and a learning rate of 0.0001 is used to update these four types of variables, completing one parameter update.

[0021] The residual function is:

[0022] ,in, Let and represent the initial arrival and arrival times of the stress wave generated by the rock fracture, respectively, acquired by the j-th microseismic sensor. N represents the initial arrival time of the elastic wave from the j-th microseismic sensor to the 3D grid node (x,y,z) of the rock mass fracture. k This represents the total number of microseismic sensors in the monitoring area.

[0023] The process of delineating the monitoring area according to the monitoring task, deploying microseismic sensors within the monitoring area, and representing the coordinates of any three-dimensional grid node within the monitoring area as (x,y,z), and assigning the three-dimensional grid node (x,y,z) to the rock mass wave velocity field V0(x,y,z), specifically includes the following:

[0024] A regular cuboid is selected as the monitoring area, with its length, width, and height defined as L meters, W meters, and H meters, respectively. A Cartesian coordinate system is established along the three directions, and the grid size is set to d meters. The monitoring area of ​​the cuboid is discretized into N... tol = (L / d)×(W / d)×(H / d) three-dimensional mesh nodes, and define the length, width and height directions as the X direction, Y direction and Z direction respectively. Then the coordinates of any three-dimensional mesh node can be represented as (x,y,z).

[0025] Microseismic sensors are deployed in the monitoring area, and the three-dimensional spatial coordinates of each sensor are recorded. For the i-th sensor, its spatial coordinates are recorded as (Sxi, Syi, Szi).

[0026] Based on the preliminary engineering geological survey data, the three-dimensional grid nodes (x,y,z) of the monitoring area are assigned the rock mass wave velocity V0(x,y,z).

[0027] This invention offers the following advantages: a microseismic location method based on adaptive wave velocity field optimization in complex rock masses can update the wave velocity field characteristics of the monitoring area based on existing data and information, and calculate the theoretical first arrival time of each node to different sensors. By calculating the residual values ​​between the theoretical first arrival time and the actual arrival time of the acquired microseismic waves, and searching for nodes that satisfy the minimum residual, the method can locate the rock mass fracture location in wave velocity regions containing cavities and complex rock masses. This method can better utilize microseismic monitoring technology to delineate the surrounding rock damage area and provide accurate real-time early warning, thus offering assistance. Attached Figure Description

[0028] Figure 1 This is a schematic diagram of the network model structure of the present invention;

[0029] Figure 2 A schematic diagram of the wave velocity characteristics of a rock mass on a hydropower slope;

[0030] Figure 3 This is a schematic diagram comparing the positioning results of the present invention with those of the classical method. Detailed Implementation

[0031] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. The components of the embodiments of this application described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the detailed description of the embodiments of this application provided below with reference to the accompanying drawings is not intended to limit the scope of protection of the claimed application, but merely represents selected embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without inventive effort are within the scope of protection of this application. The present invention will be further described below with reference to the accompanying drawings.

[0032] This invention specifically relates to a microseismic location method based on adaptive wave velocity field optimization in complex rock masses. The method employs a fusion of mesh expansion and a deep learning network based on physical information to solve the equation. Based on the wave velocity field of the monitoring area obtained previously, and combined with the physical constraints, boundary adjustments, and initial arrival times of the equation, the prior wave velocity conditions of the monitoring area are optimized to obtain an optimized wave velocity model of the engineering rock mass in the monitoring area, and the theoretical initial arrival times of stress wave propagation are calculated. Subsequently, based on the measured initial arrival times of stress waves from various sensors deployed at the engineering site, the residual between the theoretical and actual arrival times is calculated, and the point that minimizes the residual value is found, thereby achieving the location of the fracture source. Specifically, this includes the following:

[0033] (1) First, based on the monitoring task, the monitoring area is delineated. Usually, a regular cuboid is chosen as the monitoring area, with its length, width, and height determined to be L meters, W meters, and H meters, respectively. A rectangular coordinate system is established along the three directions, and the grid size is set to d meters. The monitoring area of ​​the cuboid is discretized into N... tol = (L / d)×(W / d)×(H / d) three-dimensional mesh nodes, and define the length, width and height directions as the X direction, Y direction and Z direction respectively. Then the coordinates of any node can be represented as (x,y,z).

[0034] (2) Microseismic sensors are arranged in the monitoring area. Their spatial distribution needs to be arranged in a mesh structure. The three-dimensional spatial coordinates of each sensor are recorded. For the i-th sensor, its spatial coordinates are recorded as (Sxi, Syi, Szi).

[0035] (3) Based on the preliminary engineering geological survey data, the three-dimensional grid nodes (x,y,z) of the monitoring area are assigned the rock mass wave velocity V0(x,y,z). For the excavated cavity area, an initial wave velocity of 1 m / s is assigned.

[0036] (4) Using the machine learning algorithm fused with physical information proposed in this invention, the arrival and departure times to each node are calculated, starting from each sensor. The calculation approach is as follows:

[0037] ① The equation for the three-dimensional wave velocity field:

[0038] (1)

[0039] in, Let represent the initial arrival / departure time at node (x, y, z). For the starting point (x0, y0, z0), the initial arrival / departure time satisfies:

[0040] (2)

[0041] To construct a suitable structure, the initial arrival and arrival times of nodes (x, y, z) in a non-uniform wave velocity field are equivalently deformed as follows:

[0042] (3)

[0043] In the formula, This represents the first arrival and arrival times of the node (x, y, z) calculated using an equivalent uniform wave velocity field model. It can also be replaced by the first arrival and arrival times obtained using other simplified models. This is represented as the initial arrival travel time optimization correction coefficient located at the 3D mesh node (x,y,z). and The combination of can characterize .

[0044] In addition, for the previously known rock mass wave velocity field Based on posterior information or other optimization information, the optimized rock mass wave velocity field can be expressed as follows using formula (4):

[0045] (4)

[0046] In the formula, Let the rock mass wave velocity at node (x, y, z) be the previously known wave velocity. It is represented as the wave velocity field optimization correction coefficient located at node (x,y,z).

[0047] Based on formulas (3) and (4), by making appropriate changes to formula (1), an equation that can be used to integrate physical information and machine learning algorithms is constructed:

[0048] (5)

[0049] ② Constructing the input features for machine learning network models, such as Figure 1As shown, the input to the machine learning network model consists of the travel time and wave speed adjustment coefficients for each node, with a total of L / d×W / d×H / d×2 input features. Initially, let... .

[0050] ③ Construct the intermediate layers of the network model, where the network has p layers, each with q neurons (recommended p>10, q>10). The linear and constant weights of the q-th neuron in the p-th layer are represented by W. pq With b pq .

[0051] ④ Construct the model output layer and output the results. Represented as:

[0052] (6)

[0053] The loss function of this model can then be expressed as:

[0054] (7)

[0055] In the equation, the second and third terms on the right-hand side represent the nodes for the known initial arrival and arrival correction coefficients and wave velocity correction coefficients, respectively, and their quantities are respectively... and .

[0056] ⑤ Based on the loss function calculation results, calculate the damage function for W using the steepest gradient descent method. pq , b pq F t (x,y,z), The partial derivatives are calculated, and a learning rate of 0.0001 is used to update these four types of variables, completing one parameter update.

[0057] ⑥ Repeat steps ①-⑤, continuously updating the four types of variables mentioned. Stop the process when the maximum number of repetitions is reached or the final error function result is less than the threshold. Then, obtain the F... t Substituting (x, y, z) into formula (3) yields the initial arrival and arrival times at each node. Substitute these values ​​into formula (4) to obtain the optimized wave velocity values ​​for each node.

[0058] (5) Using the method mentioned above, calculate the initial arrival and arrival times from the i-th microseismic sensor to the three-dimensional mesh node (x,y,z) in the region. .

[0059] (6) When a microseismic event occurs in the monitoring area, the arrival and arrival times of the rock mass fracture stress wave collected on-site by the i-th microseismic sensor are: .

[0060] (7) Establish the residual function between theoretical arrival time and actual arrival time:

[0061] (8)

[0062] in, and These represent the stress waves generated by rock fractures collected by the i-th and j-th microseismic sensors, respectively. and N represents the initial arrival and arrival times of the rock mass fracture elastic waves from the i-th and j-th microseismic sensors to the three-dimensional grid node (x,y,z) calculated in step (4); k This represents the total number of microseismic sensors in the monitoring area.

[0063] (8) The calculated results From all results, select the first n nodes corresponding to the smallest values, i.e., (x1, y1, z1), ..., (x n ,y n ,z n Typically, n=10 is chosen, and the coordinates (x, y) of the rupture origin point are calculated using formula (9). f ,y f ,z f ):

[0064] (9)

[0065] This invention uses a hydropower slope as an example to illustrate specific implementation methods.

[0066] (1) Based on the actual monitoring range, the study area to be monitored (300 m × 300 m × 264 m) was determined. The monitoring area was divided into 24,009,265 grid nodes according to a grid size of 1 m. Based on previous exploration data, specific rock masses in certain areas were assigned corresponding rock mass wave velocities, and the wave velocity of excavated cavities was set to 1 m / s. The wave velocity zoning is as follows: Figure 2 As shown in Table 1, the wave velocity in each region is as follows.

[0067] Table 1. Wave velocity table for each region

[0068] 1 <![CDATA[V1]]> 4500 2 <![CDATA[V2]]> 5100 3 <![CDATA[V3]]> 4100 4 <![CDATA[V4]]> 4800

[0069] (2) Microseismic monitoring sensors were deployed in the monitoring area. A total of 13 microseismic sensors were deployed, and their corresponding spatial coordinates are recorded in Table 2.

[0070] Table 2. Spatial Coordinates of Microseismic Sensing

[0071] 1 20 191 187 2 120 257 178 3 73 228 184 4 98 191 187 5 116 156 186 6 136 140 97 7 61 190 99 8 86 162 98 9 159 150 97 10 120 231 101 11 157 281 101 12 201 103 44 13 182 189 42

[0072] (3) Based on the initial set nodal wave velocity values, the first arrival time of the stress wave reaching each node is calculated using the method mentioned in this invention, taking the 13 microseismic sensors in the area as the starting point, and the optimized wave velocity field is obtained to establish a first arrival time database.

[0073] (4) A blasting test was conducted at the engineering site, and the initial arrival time of the stress wave to each sensor was extracted. Using the method of this invention, the residual function was calculated, and the grid node that satisfies the minimum residual was searched to locate the blasting position. The results were compared with the commonly used uniform velocity positioning model. The positioning results are shown in Table 3. The results show that the positioning accuracy was improved by 50%.

[0074] Table 3. Positioning Results Table

[0075] Uniform velocity model positioning 16.2m 17.3m 17.8m Considering voids and wave velocity stratification 9.3m 8.2m 9.2m Reduction amount (%) 42.59 52.60 48.31

[0076] (5) Based on the rock mass fracture microseismic information monitored in May of a certain year, the arrival times of the stress waves of the actual microseismic events were extracted. The residual function was calculated, and the grid nodes that minimized the residuals were searched to locate the microseismic positions. The results are shown in... Figure 3 The positioning results show that, compared with the commonly used uniform velocity positioning model, the positioning results of this invention are closer to the interlayer slip zone of the slope affected by engineering disturbances. Figure 3 (Right), and all of them use the velocity model, and their positioning results are further away from the interlayer fault zone.

[0077] The above description is merely a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the forms disclosed herein and should not be construed as excluding other embodiments. It can be used in various other combinations, modifications, and improvements, and can be altered within the scope of the concept described herein through the above teachings or related technologies or knowledge. Modifications and variations made by those skilled in the art that do not depart from the spirit and scope of the present invention should be within the protection scope of the appended claims.

Claims

1. A microseismic location method with adaptive wave velocity field optimization in complex rock masses, characterized in that: The microseismic location method includes: The monitoring area is delineated according to the monitoring task, microseismic sensors are deployed in the monitoring area, and the coordinates of any three-dimensional grid node in the monitoring area are represented as (x,y,z). The three-dimensional grid node (x,y,z) is assigned the rock mass wave velocity field V0(x,y,z). Starting from each sensor, the initial arrival and arrival times of the rock mass fracture elastic wave at each three-dimensional grid node are calculated, resulting in the initial arrival and arrival times of the wave at the three-dimensional grid node (x,y,z) of the i-th microseismic sensor. When a microseismic event is transmitted in the monitoring area, the arrival time of the rock fracture stress wave collected on-site by the i-th microseismic sensor is: ; Establish the residual function between theoretical arrival time and actual arrival time, select the first n nodes corresponding to the smallest value among all the results of the calculated residual function, and then calculate the coordinates of the rupture source point; The process of delineating the monitoring area according to the monitoring task, deploying microseismic sensors within the monitoring area, and representing the coordinates of any three-dimensional grid node within the monitoring area as (x,y,z), and assigning the three-dimensional grid node (x,y,z) to the rock mass wave velocity field V0(x,y,z), specifically includes the following: A regular cuboid is selected as the monitoring area, with its length, width, and height defined as L, W, and H, respectively. A rectangular coordinate system is established along the three directions, and the grid size is set to d meters. The monitoring area of ​​the cuboid is then discretized into N... tol = (L / d)×(W / d)×(H / d) three-dimensional mesh nodes, and define the length, width and height directions as the X direction, Y direction and Z direction respectively. Then the coordinates of any three-dimensional mesh node can be represented as (x,y,z). Microseismic sensors are deployed in the monitoring area, and the three-dimensional spatial coordinates of each sensor are recorded. For the i-th sensor, its spatial coordinates are recorded as (Sxi, Syi, Szi). Based on the preliminary engineering geological survey data, the three-dimensional grid nodes (x,y,z) of the monitoring area are assigned the rock mass wave velocity field V0(x,y,z).

2. The microseismic location method for adaptive wave velocity field optimization in complex rock masses according to claim 1, characterized in that: The calculation of the initial arrival time to each 3D mesh node, starting from each sensor, specifically includes the following: A1. Equation for a three-dimensional wave velocity field , Let x0, y0, z0 represent the initial arrival / departure time when reaching the 3D mesh node (x, y, z). Then, the initial arrival / departure time for the starting point (x0, y0, z0) is... ; A2. The initial arrival and arrival times of the three-dimensional mesh nodes (x, y, z) in the non-uniform wave velocity field are equivalently deformed as follows: , This represents the initial arrival and arrival times of the three-dimensional mesh nodes (x, y, z) obtained using the equivalent uniform wave velocity field model. It is represented as the initial arrival travel time optimization correction coefficient located at the 3D mesh node (x,y,z); A3. For the previously known rock mass wave velocity field V0(x,y,z), the optimized rock mass wave velocity field is expressed as follows based on posterior or optimization information: , The previously known rock mass wave velocity field is represented as a three-dimensional grid node (x,y,z). It is represented as the wave velocity field optimization correction coefficient located at the three-dimensional grid node (x,y,z); A4. According to the formula and For the equation of function The transformation yields the equation used to fuse physical information and machine learning algorithms. ; A5. Construct the input features of the machine learning network model, construct the intermediate and output layers of the network model, and calculate the results based on the loss function; Repeat steps A1-A5. Stop repeating when the maximum number of repetitions is reached or the final error function result is less than the threshold. Substitute into the formula In the process, when the initial arrival time to each 3D mesh node is obtained, the obtained... Substitute into the formula The optimized wave velocity values ​​of each three-dimensional mesh node are obtained.

3. The microseismic location method for adaptive wave velocity field optimization in complex rock masses according to claim 2, characterized in that: Step A5 specifically includes the following: Input features for constructing a machine learning network model: Input the travel time and wave speed adjustment coefficients of each 3D grid node into the machine learning network model. At the initial input, let... ; Constructing the intermediate layers of the network model: The network model has a total of p layers, with q neurons in each layer. The weights of the i-th neuron in the j-th layer are represented by W as linear and constant terms, respectively. pq With b pq ; Construct the model output layer: The output result is represented as The loss function of the network model is expressed as: ; Based on the loss function calculation results, the loss function is calculated on W using the steepest gradient descent method. pq b pq partial derivatives and The four types of variables are updated using a learning rate of 0.0001, thus completing one parameter update.

4. The microseismic location method for adaptive wave velocity field optimization in complex rock masses according to claim 1, characterized in that: The residual function is: ,in, Let represent the initial arrival and arrival times of the stress wave generated by the rock fracture, as collected by the j-th microseismic sensor. N represents the initial arrival time of the elastic wave from the j-th microseismic sensor to the 3D grid node (x,y,z) indicating the rock mass fracture. k This represents the total number of microseismic sensors in the monitoring area.

Images