Microseismic source end-to-end intelligent positioning method for complex geological conditions

Abstract

The application discloses a microseismic source end-to-end intelligent positioning method for complex geological conditions and belongs to the technical field of microseismic positioning. In view of the industry pain points that the traditional microseismic source positioning method has limited accuracy and relies on initial values and a large amount of labeled data, the application adopts a ray tracing method to construct a high-fidelity P-wave first arrival time simulation dataset based on actual geological data, replaces high-cost numerical simulation and field labeled data, designs an adaptive residual neural network, effectively solves the problems of gradient disappearance of a deep network and unstable training through residual connection, batch normalization and an adaptive matrix estimation Adam optimizer, and realizes end-to-end mapping from monitoring point coordinates and P-wave first arrival time to source three-dimensional coordinates. Test results show that the method has excellent positioning accuracy and anti-interference robustness and is suitable for microseismic high-precision real-time monitoring under complex geological conditions.

Description

Technical Field

[0001] This invention belongs to the field of microseismic source localization technology, specifically relating to an end-to-end intelligent localization method for microseismic sources under complex geological conditions. Background Technology

[0002] Microseismic monitoring technology, as an effective means of obtaining the dynamic response of underground rock masses, has been widely used in fields such as coal mining, oil extraction, and water conservancy and hydropower projects. In deep underground engineering, the synergistic effect of in-situ stress and mining activities can easily induce rock mass disturbance and crack propagation, ultimately leading to geological disasters such as rockbursts and roof collapses. Microseismic source location is the core link of microseismic monitoring technology. It captures the elastic waves excited by rock mass fractures through deployed sensors and inverts the spatial location of the seismic source, which is the core foundation for early warning and prevention of deep geological disasters.

[0003] Currently, mainstream microseismic source location methods are mainly divided into two categories: one is the traditional model-driven location method, mainly represented by the Geiger method, linear inversion algorithm, Newton's iteration method, and grid search method. Among them, Newton's iteration method relies on good initial values ​​and is prone to getting trapped in local extrema, resulting in extremely large location errors in layered non-homogeneous media; although the grid search method avoids the local extrema problem, it is limited by the grid resolution, making it difficult to balance location accuracy and computational efficiency. At the same time, most traditional methods are based on uniform wave velocity models, simplifying the elastic wave propagation path to a straight line. However, actual strata are mostly layered non-homogeneous media, and elastic waves will be reflected, transmitted, and refracted at the interlayer interfaces. Such simplification introduces a large systematic error, which cannot meet the requirements of high-precision location.

[0004] Another category is data-driven deep learning-based localization methods, which have shown great potential in the field of microseismic localization due to their excellent nonlinear fitting capabilities. Existing research has proposed microseismic localization models based on CNNs, U-nets, and graph neural networks. However, these methods all rely on a large amount of measured microseismic data with source location labels for training. Acquiring on-site measured labeled data is extremely costly, and much of it contains unlabeled and invalid data, severely limiting the engineering application of the models. Although existing research has proposed using synthetic data to replace measured data for training networks, synthetic data is mostly generated based on numerical simulation methods. Simulating a single model is time-consuming and extremely inefficient. Furthermore, existing network architectures are mostly designed for waveform data and have poor adaptability to the low-dimensional features of P-wave first arrival, easily leading to gradient vanishing and training instability.

[0005] Therefore, given the geological characteristics of layered and non-homogeneous media in coal mines, developing a microseismic source location method that requires no large amount of measured and labeled data, has high positioning accuracy, high training efficiency, and strong anti-interference ability is an urgent technical problem to be solved in the field of safety monitoring. Summary of the Invention

[0006] To address the problems existing in the prior art, this invention provides an end-to-end intelligent positioning method for microseismic sources under complex geological conditions. This method overcomes the shortcomings of existing technologies, such as low microseismic positioning accuracy, reliance on initial values ​​and large amounts of high-cost labeled data, and unstable network training. It achieves high-precision, highly adaptable, and low-cost end-to-end positioning of microseismic events.

[0007] To achieve the above objectives, the technical solution adopted by this invention is: an end-to-end intelligent positioning method for microseismic sources under complex geological conditions, comprising the following steps: S1. Construction of a three-dimensional model of layered media: Based on the actual geological survey data of the target area, a three-dimensional layered media model matching the on-site strata distribution is constructed. According to the physical and mechanical parameters such as the thickness, density, elastic modulus, and Poisson's ratio of each rock layer, the P-wave velocity of the corresponding strata is calculated using the elasticity formula.

[0008] S2. Construction of a high-fidelity simulation dataset based on ray tracing: Based on Snell's law, a three-dimensional ray tracing is performed using a trial-and-error method combined with a bisection iterative strategy. The initial emission angle of the ray is iteratively optimized based on Fermat's principle to determine the minimum travel time propagation path. The critical angle condition is monitored in real time, and the travel times of the direct elastic wave and the first wave at the critical angle are calculated in parallel. The minimum value is taken as the first arrival time of the P-wave. Source points are set up in the three-dimensional model space at a preset step size. The first arrival time of the P-wave at the target monitoring point under different sources is obtained through the ray tracing algorithm. A simulation dataset containing the spatial coordinates of the monitoring point, the first arrival time of the P-wave, and the corresponding three-dimensional coordinates of the source is constructed.

[0009] S3. Adaptive Residual Neural Network Model Construction: An adaptive residual neural network model with fully connected layers as the main body is constructed. The model input is a 16-dimensional feature vector composed of the three-dimensional coordinates of four monitoring points and the corresponding P-wave first arrival times. The output is the three-dimensional coordinates of the seismic source. The model adopts a 10-layer progressively reduced dimensionality fully connected hidden layer architecture, introduces residual connections between key layers, sets a batch normalization layer after each hidden layer, uses LeakyReLU and ReLU as activation functions, and sets a gradient clipping strategy to avoid gradient vanishing and gradient exploding.

[0010] S4. Training and Optimization of the Adaptive Residual Neural Network Model: The Z-score normalization method is used to preprocess the input features and output target of the dataset. The preprocessed dataset is randomly shuffled and divided into training and test sets in a 9:1 ratio. The Adam adaptive moment estimation algorithm is used as the optimizer, combined with the learning rate warm-up and step decay strategy. The mean squared error with L2 regularization is used as the loss function to complete the iterative training of the adaptive residual neural network model. The model's localization accuracy is verified through the test set.

[0011] S5. Microseismic Event Location: For microseismic events to be located by on-site monitoring, the spatial coordinates and P-wave arrival times collected by the target monitoring point are extracted. After standardized preprocessing, the data are input into the trained adaptive residual neural network model. The model directly outputs the three-dimensional spatial coordinates of the source of the microseismic event, thus completing the microseismic location.

[0012] Furthermore, in step S1, the P-wave velocity is calculated using elastic mechanical parameters, and the calculation formula is as follows: in, For P-wave velocity, The elastic modulus of the formation. The Poisson's ratio of the strata. The density is the medium.

[0013] Furthermore, in step S2, the expression for Snell's law in three-dimensional space is: in, For the first The angle of incidence in the layer, For the corresponding layer P-wave velocity, These are the ray parameters.

[0014] During ray tracing, the first The formulas for calculating the propagation path length and travel time of elastic waves within a layer are as follows: in, For the first Vertical layer thickness This is the propagation path length within this layer. This is the travel time propagated within this layer.

[0015] Total time The sum of the travel times of the wave across all the layers it passes through: in, This represents the total number of strata traversed by the ray.

[0016] Furthermore, in step S2, when the ray enters the high-speed layer from the low-velocity layer, the critical angle... The calculation formula is: in, This refers to the propagation of waves from a low-speed medium to a high-speed medium.

[0017] The formula for calculating the first wave time is: in, The total horizontal distance between the seismic source and the receiver. , These represent the vertical distances from the seismic source and receiver to the interface, respectively.

[0018] During ray tracing, the critical angle condition is monitored in real time, and the travel times of the direct wave and the first wave are calculated in parallel. The minimum value is taken as the final P-wave arrival time to ensure the physical correctness of the travel time calculation.

[0019] Furthermore, in step S3, the adaptive residual neural network model includes 10 fully connected hidden layers and 1 output layer, adopting a progressive dimensionality reduction architecture: the first hidden layer is configured with 2048 neurons, the second hidden layer with 1024 neurons, the third hidden layer with 512 neurons, the fourth hidden layer with 256 neurons, the fifth hidden layer with 128 neurons, the sixth hidden layer with 64 neurons, the seventh hidden layer with 32 neurons, the eighth hidden layer with 16 neurons, the ninth hidden layer with 8 neurons, and the tenth hidden layer with 4 neurons, progressively refining and compressing feature information; the output layer is configured with 3 neurons, corresponding to the x, y, and z coordinate values ​​of the earthquake source in three-dimensional space, respectively.

[0020] Furthermore, in step S3, the residual connection directly transmits shallow features to the deep network via a shortcut path. layer to the first During the layer propagation process, residual paths are constructed, and shallow features are projected onto the same dimensional space as deep features through the residual projection matrix, thereby achieving cross-layer feature reuse.

[0021] The formula for batch normalization is: in, , The mean and variance of the current batch. , For learnable transformation parameters, is the numerical stability constant.

[0022] The formula for calculating the LeakyReLU activation function is: in, The slope coefficient is set to 0.01. This design maintains the advantages of ReLU sparse activation while ensuring that the gradient in the negative interval is not zero, avoiding the problem of neuron "death" and promoting a more stable training process.

[0023] To prevent gradient explosion during training, a gradient pruning strategy is implemented, limiting the gradients of all parameters to the range [-1, 1]. This technique ensures the numerical stability of the training process, which is particularly important in deep network architectures. The specific formula is as follows: Where g is the gradient of the parameters before clipping, g c This represents the gradient of the clipped parameters.

[0024] Furthermore, in step S4, the input and output data are preprocessed using the Z-score normalization method, and the calculation formula is as follows: in, , The input features are the mean and standard deviation. , To output the mean and standard deviation of the target.

[0025] Furthermore, in step S4, the input and output data are preprocessed using the Z-score normalization method, and the calculation formula is as follows: in, , The input features are the mean and standard deviation. , To output the mean and standard deviation of the target.

[0026] Furthermore, in step S4, the Adam optimizer updates its parameters by calculating the first-order moment estimate and the second-order moment estimate of the gradient. The calculation formula is as follows: in, and These represent the first-order moment estimate and the second-order moment estimate of the gradient, respectively. For the current gradient, , This is the attenuation rate hyperparameter.

[0027] The formula for calculating the estimated value after deviation correction is: The final parameter update formula is: in, For learning rate, To prevent small positive values ​​from being divided by zero, These are the parameters to be optimized at the current time step.

[0028] Furthermore, in step S4, the learning rate strategy is as follows: a warm-up mechanism is adopted in the early stage of training, and the learning rate is linearly increased from 0 to the initial value of 0.02 in the first 200 iterations to avoid gradient instability in the early stage of training; after the warm-up, a step decay strategy is adopted, and the current learning rate is multiplied by a decay factor of 0.9 every 100 training iterations to ensure fine adjustment of parameters and smooth convergence in the later stage of training.

[0029] Furthermore, in step S4, the loss function is calculated using the following formula: Where J is the original loss function for mean squared error. This is the regularization penalty coefficient. The total number of network layers. Let F be the weight matrix of the i-th layer of the neural network, and F be the Frobenius norm.

[0030] Compared with the prior art, the present invention has the following advantages: 1. Based on actual geological data and Snell's law, this invention uses ray tracing to construct a high-fidelity simulation dataset of P-wave first arrival times. It eliminates the need for costly field annotation data and time-consuming numerical simulations. The calculation of a single travel time only takes milliseconds, and the maximum error of P-wave arrival time calculation is extremely small. While ensuring the accuracy of the dataset, it significantly reduces the data cost of model training.

[0031] 2. This invention addresses the 16-dimensional low-dimensional spatiotemporal features of microseismic localization in target areas by custom-designing a 10-layer progressive dimensionality reduction adaptive residual neural network model. Through residual connections, batch normalization, LeakyReLU and ReLU activation functions, and gradient pruning strategies, it effectively solves the problems of gradient vanishing and training instability in deep networks. Combined with the Adam optimizer and a refined learning rate strategy, the model converges quickly and has strong training stability. Compared with standard gradient descent optimized networks, the average localization error is significantly reduced.

[0032] 3. The method of the present invention is adapted to the geological characteristics of layered non-homogeneous media, can be directly adapted to different geological models, has strong generalization ability, does not require complex manual parameter adjustment, has high computational efficiency, can realize real-time location of microseismic events, provides reliable technical support for early warning and prevention of geological disasters, and has important engineering application value. Attached Figure Description

[0033] Figure 1 This is a diagram showing the architecture design of the adaptive residual neural network model in this invention.

[0034] Figure 2 This is a distribution diagram of the source location error of the adaptive residual neural network model in this invention.

[0035] Figure 3 This is a heatmap of the positioning error of the network positioning results in this invention. Detailed Implementation

[0036] The present invention will be further described below.

[0037] The intelligent microseismic source localization method based on ray tracing and adaptive residual neural network provided by this invention has the following specific implementation steps: S1. Construction of a three-dimensional model of layered medium: Collect geological survey data of the target area, clarify the distribution of strata, the thickness, density, elastic modulus, Poisson's ratio and other physical and mechanical parameters of each rock layer, construct a three-dimensional layered medium model that is completely matched with the strata on site, calculate the P-wave velocity of each rock layer through the P-wave velocity calculation formula, and provide a basic wave velocity model for subsequent ray tracing.

[0038] S2. Construction of a high-fidelity simulation dataset based on ray tracing: S2.1 Determine the layout scheme of the monitoring points of the on-site microseismic monitoring system, select 4 core monitoring points, record their three-dimensional spatial coordinates, and map them to the constructed three-dimensional layered medium model.

[0039] S2.2 Within the monitoring space of the three-dimensional model, seismic source points are evenly distributed along the x, y, and z directions at preset step lengths to cover the entire monitoring area.

[0040] S2.3. Based on Snell's law, a three-dimensional ray tracing method is adopted using a trial firing method combined with a bisection iterative strategy: the coordinates of the seismic source and the monitoring point are initialized, and the horizontal and vertical distances are calculated. If the seismic source and the monitoring point are in the same stratum, the travel time of the straight wave is directly calculated; if they are in different strata, the initial emission angle of the ray is optimized by iterative optimization using the bisection method, and the incident angle, propagation path length and travel time of each layer are calculated. The total travel time of the straight wave is obtained by summing them up.

[0041] S2.4 Real-time monitoring of the incident angle of the ray. When the critical angle is reached, the travel time of the first wave is calculated. The travel time of the direct wave and the travel time of the first wave are compared in parallel. The minimum value is taken as the arrival time of the P wave at the monitoring point corresponding to the source.

[0042] S2.5. Traverse all deployed seismic source points, complete the calculation of the first arrival time of the P-wave for all sources, and construct a dataset. Each dataset contains the three-dimensional coordinates of four monitoring points, the corresponding first arrival time of the P-wave, and the three-dimensional coordinate labels of the seismic source corresponding to that dataset.

[0043] S3. Adaptive Residual Neural Network Model Construction: An adaptive residual neural network model was constructed. The main body of the model consists of 10 fully connected hidden layers and 1 output layer. The input layer has a dimension of 16, corresponding to the x, y, and z coordinates of the 4 monitoring points and the first arrival time of the P wave, totaling 16 dimensions of features. The first hidden layer is configured with 2048 neurons, the second hidden layer with 1024 neurons, the third hidden layer with 512 neurons, the fourth hidden layer with 256 neurons, the fifth hidden layer with 128 neurons, the sixth hidden layer with 64 neurons, the seventh hidden layer with 32 neurons, the eighth hidden layer with 16 neurons, the ninth hidden layer with 8 neurons, the tenth hidden layer with 4 neurons, and the output layer with 3 neurons, corresponding to the x, y, and z coordinates of the earthquake source.

[0044] Residual connections are set between the 1st hidden layer and the 3rd hidden layer, the 3rd hidden layer and the 5th hidden layer and the 7th hidden layer, respectively, and the feature dimension is matched by the projection matrix. A batch normalization layer is set after each hidden layer, and the activation functions are LeakyReLU and ReLU. The gradient clipping threshold is set to [-1,1].

[0045] S4. Training optimization of the adaptive residual neural network model: S4.1. The Z-score standardization method is used to standardize the input features and output labels of the dataset to eliminate the influence of different physical dimensions.

[0046] S4.2 Randomly shuffle the standardized dataset and divide it into a training set and a test set in a 9:1 ratio. The training set is used to update model parameters, and the test set is used to verify the model's localization accuracy.

[0047] S4.3 Initialize model parameters, set the initial learning rate of the Adam optimizer to 0.02. =0.9, =0.999; set the learning rate warm-up steps to 200, the step decay step size to 50, and the decay factor to 0.9; set the L2 regularization coefficient of the loss function. =10 -4 The total number of iterations was 2000.

[0048] S4.4 Start model training. In each iteration, the training set data is input into the model in batches. The predicted source coordinates are obtained through forward propagation, the loss function value is calculated, the model parameters are updated through backpropagation, and the learning rate adjustment and gradient clipping are performed simultaneously.

[0049] S4.5 After training, use test set data to verify the model performance, calculate the average positioning error in three-dimensional space and the positioning error in each direction. If the accuracy meets the engineering requirements, save the trained model weights.

[0050] S5. Microseismic Event Location: For microseismic events acquired by the on-site microseismic monitoring system, the first arrival time of the P-wave at four target monitoring points is picked up using the AIC algorithm. Combined with the three-dimensional spatial coordinates of the monitoring points, a 16-dimensional input feature vector is formed. After standardization preprocessing consistent with the training set, it is input into the trained adaptive residual neural network model. The model directly outputs the three-dimensional spatial coordinates of the source of the microseismic event, thus completing the microseismic location.

[0051] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. An end-to-end intelligent positioning method for microseismic sources under complex geological conditions, characterized in that, Includes the following steps: S1. Construction of a three-dimensional model of layered medium: Based on the actual geological survey data of the target area, a three-dimensional layered medium model matching the on-site strata distribution is constructed, and the P-wave velocity of the corresponding strata is calculated according to the physical and mechanical parameters of each rock layer. S2. Construction of a high-fidelity simulation dataset based on ray tracing: Based on Snell's law, a three-dimensional ray tracing is performed using a trial-and-error method combined with a bisection iterative strategy. The travel times of the direct elastic wave and the critical angle first wave are calculated in parallel. The minimum value is taken as the arrival time of the P-wave. A simulation dataset containing the spatial coordinates of the monitoring points, the arrival time of the P-wave, and the three-dimensional coordinates of the corresponding seismic source is constructed. S3. Construction of Adaptive Residual Neural Network Model: Build an adaptive residual neural network model with fully connected layers as the main body. The hidden layers of the model introduce residual connections and batch normalization layers. LeakyReLU and ReLU are used as activation functions, and gradient pruning strategy is set. S4. Training and optimization of the adaptive residual neural network model: After standardizing and preprocessing the dataset, the training set and test set are divided. The Adam adaptive moment estimation algorithm is used as the optimizer, combined with the learning rate warm-up and step decay strategy. The mean squared error with L2 regularization is used as the loss function to complete the iterative training of the adaptive residual neural network model. S5. Microseismic Event Location: Input the spatial coordinates of the monitoring points corresponding to the target microseismic event and the P-wave first arrival time data into the trained adaptive residual neural network model, and directly output the three-dimensional spatial coordinates of the source of the microseismic event to complete the microseismic source location.

2. The method according to claim 1, characterized in that, In step S1, the P-wave velocity is calculated using elastic mechanical parameters, and the calculation formula is as follows: in, For P-wave velocity, The elastic modulus of the formation. The Poisson's ratio of the strata. The density is the medium.

3. The method according to claim 1, characterized in that, In step S2, the expression for Snell's law in three-dimensional space is: in, For the first The angle of incidence in the layer, For the corresponding layer P-wave velocity, For ray parameters; During ray tracing, the first The formulas for calculating the propagation path length and travel time of elastic waves within a layer are as follows: in, For the first Vertical layer thickness This is the propagation path length within this layer. This refers to the travel time propagated within this layer; Total time The sum of the travel times of the wave across all the layers it passes through: in, This represents the total number of strata traversed by the ray.

4. The method according to claim 1, characterized in that, In step S2, when the ray enters the high-speed layer from the low-velocity layer, the critical angle... The calculation formula is: in, That is, the wave propagates from a low-speed medium to a high-speed medium; The formula for calculating the first wave time is: in, The total horizontal distance between the seismic source and the receiver. , These are the vertical distances from the seismic source and receiver to the interface, respectively. During ray tracing, the critical angle condition is monitored in real time, and the travel times of the direct wave and the first wave are calculated in parallel. The minimum value is taken as the final P-wave arrival time.

5. The method according to claim 1, characterized in that, In step S3, the input to the adaptive residual neural network model is a 16-dimensional feature vector composed of the three-dimensional coordinates of the four monitoring points and the corresponding P-wave arrival times, and the output is the three-dimensional coordinates of the earthquake source. The model contains 10 fully connected hidden layers and 1 output layer, and adopts a progressive dimensionality reduction architecture: the first hidden layer is configured with 2048 neurons, the second hidden layer with 1024 neurons, the third hidden layer with 512 neurons, the fourth hidden layer with 256 neurons, the fifth hidden layer with 128 neurons, the sixth hidden layer with 64 neurons, the seventh hidden layer with 32 neurons, the eighth hidden layer with 16 neurons, the ninth hidden layer with 8 neurons, the tenth hidden layer with 4 neurons, and the output layer with 3 neurons, corresponding to the x, y, and z coordinate values ​​of the earthquake source in three-dimensional space, respectively.

6. The method according to claim 1, characterized in that, In step S3, the formula for batch normalization is: in, , The mean and variance of the current batch. , For learnable transformation parameters, It is the numerical stability constant; The formula for calculating the LeakyReLU activation function is: in, It is a negative slope coefficient; The gradient pruning strategy restricts the gradients of all parameters to the range [-1, 1]. The specific formula is as follows: Where g is the gradient of the parameters before clipping, g c This represents the gradient of the clipped parameters.

7. The method according to claim 1, characterized in that, In step S4, the input and output data are preprocessed using the Z-score normalization method, and the calculation formula is as follows: in, , The input features are the mean and standard deviation. , To output the mean and standard deviation of the target.

8. The method according to claim 1, characterized in that, In step S4, the Adam optimizer updates its parameters by calculating the first-order moment estimate and the second-order moment estimate of the gradient. The calculation formula is as follows: in, and These represent the first-order moment estimate and the second-order moment estimate of the gradient, respectively. For the current gradient, , This is the attenuation rate hyperparameter; The formula for calculating the estimated value after deviation correction is: The final parameter update formula is: in, For learning rate, To prevent small positive values ​​from being divided by zero, These are the parameters to be optimized at the current time step.

9. The method according to claim 1, characterized in that, In step S4, the learning rate strategy is as follows: a warm-up mechanism is adopted in the early stage of training, and the linear learning rate is increased from 0 to the initial value of 0.02 in the first 200 iterations; after the warm-up, a step decay strategy is adopted, and the current learning rate is multiplied by a decay factor of 0.9 every 100 training iterations.

10. The method according to claim 1, characterized in that, In step S4, the loss function is calculated using the following formula: Where J is the original loss function for mean squared error. This is the regularization penalty coefficient. The total number of network layers. Let F be the weight matrix of the i-th layer of the neural network, and F be the Frobenius norm.

Images