A three-dimensional positioning method of slope microseismic events considering complex stratum structure and terrain effect

Abstract

The present application provides a kind of complex stratum structure and topographic effect considering slope microseismic event 3D positioning method, the technical scheme includes S1: the theoretical travel time model suitable for complex slope environment is built;S2: the actual travel time difference is extracted from slope monitoring data;S3: based on the 3D positioning of theoretical and observed travel time difference matching, the ambient noise surface wave tomography technology is applied to slope velocity modeling, without artificial seismic source, high resolution velocity structure can be obtained;Realize 3D ray tracing under the constraint of digital elevation model, ensure that wave propagation path meets physical law, positioning accuracy reaches 3 m, can reach 1.5 m in dense array coverage area, compared with traditional method, it is improved by one order of magnitude;Successfully applied to complex environment with more than 700m topographic relief and extreme velocity variation, can be applied to microseismic monitoring in various complex geological environments.

Description

Technical Field

[0001] This invention relates to the field of geophysical monitoring and geological engineering technology, specifically to a three-dimensional localization method for microseismic events that takes into account complex slope strata structure and topographic effects. It is particularly suitable for microseismic monitoring of mountain slopes with strong geological structures, extreme topographic undulations, and high heterogeneity. Technical Background

[0002] Slope stability monitoring is a crucial means of disaster prevention and mitigation. Microseismic monitoring technology, by monitoring and analyzing the elastic waves generated during rock mass deformation, enables real-time monitoring of internal slope deformation. Accurate source location is the core of microseismic monitoring, directly determining whether potential slip surfaces can be correctly identified and slope stability assessed. However, the mountainous slope environment presents significant challenges to traditional microseismic location methods.

[0003] First, mountain slopes possess extremely complex geological structures. From the surface downwards, they typically include loose deposits, weathered zones, fractured zones, and finally relatively intact bedrock. Elastic wave velocities can vary by an order of magnitude within a range of tens of meters, from 200 m / s to 2500 m / s. Traditional methods, employing uniform or simple layered velocity models, cannot reflect this strong velocity non-uniformity, leading to significant deviations between calculated theoretical travel times and actual values.

[0004] Secondly, the undulating terrain of mountain slopes significantly affects the propagation of microwaves. Since the detectors are deployed on the surface while the seismic source is underground, elastic waves must propagate under the constraints of the terrain and cannot travel through the air. In areas with complex terrain, the actual propagation path deviates considerably from the straight-line assumption; ignoring this effect would lead to significant positioning errors.

[0005] Existing microseismic location methods, whether based on travel time-based Geiger methods or waveform-based methods, are all based on simplified velocity models and the assumption of straight-line rays. In slope environments with large topographic relief and complex geological structures, the location errors of these methods typically reach 20-30 m or even greater, failing to meet the accuracy requirements of engineering applications. Therefore, developing a high-precision location method that can accurately consider the influence of complex strata and topography is of great significance. Summary of the Invention

[0006] The present invention aims to provide a high-precision three-dimensional positioning method for slope microseismic analysis, which significantly improves positioning accuracy by accurately depicting complex geological structures and taking into account topographic effects.

[0007] The technical solution adopted in this invention includes three core steps: theoretical travel time model construction, microseismic event identification and travel time extraction, and three-dimensional mesh search and positioning. The specific steps are as follows:

[0008] S1: Construct a theoretical travel-time model adapted to complex slope environments, including:

[0009] S1.1: A high-resolution three-dimensional velocity model reflecting the complex geological structure of the slope is constructed using environmental noise surface wave tomography with a dense seismic array. Specifically, a dense seismic array (e.g., with a spacing of 50 m, determined based on actual site conditions) is deployed in the slope deformation area to collect environmental noise data for more than 72 hours. The noise data is preprocessed, including instrument response removal, time-domain normalization, and spectral whitening. The cross-correlation function between all station pairs is calculated, and the empirical Green's function is extracted by phase-weighted superposition. The Rayleigh wave dispersion curve is extracted from the cross-correlation function using the MASW multichannel surface wave analysis method combined with high-resolution Radon transform. The dispersion curve is converted into an S-wave velocity-depth profile through damped least squares inversion. The one-dimensional velocity profiles of each path are constructed into a three-dimensional velocity model using natural neighborhood interpolation.

[0010] S1.2: Based on the constraints of the Digital Elevation Model (DEM), ray tracing is performed using the fast travel method. Specifically, this includes: establishing a three-dimensional computational domain that includes the slope terrain; using the DEM data as the upper boundary of the computational domain and excluding grid nodes above the terrain surface; applying the fast travel method in the terrain-constrained three-dimensional velocity field to calculate the travel time field from each grid node to all stations; and tracing backward along the negative travel gradient direction to determine the minimum travel time ray path that follows the terrain constraints and velocity structure.

[0011] S1.3: Integrating the velocity model and terrain-constrained paths, a theoretical travel time database from each potential seismic source location to the monitoring station is established. The formula for calculating the theoretical travel time is as follows:

[0012]

[0013] Where T(x,y,z,S) is the theoretical travel time from the assumed source location (x,y,z) to station S; P is the shortest path (i.e., the path with the minimum travel time); (i,j) are a pair of adjacent grid nodes on the shortest path; d i,j The Euclidean distance between adjacent nodes i and j; v i and v j Let be the wave velocity values ​​at nodes i and j. Let Π be the set of all possible propagation paths; |∇v| (i,j) β is the velocity gradient magnitude between nodes i and j; β is the velocity gradient correction coefficient, used to account for the impact of velocity changes on the path; α is the weighting factor, adjusting the influence of the velocity difference term; cosθ ij It is the cosine of the angle between the ray direction and the velocity transition direction, used to correct the propagation effect in non-uniform media.

[0014] S2: Extract the actual travel time difference from the slope monitoring data, including:

[0015] S2.1: Using both signal-to-noise ratio (SNR) and power spectral density (PSD) as criteria, effective microseismic events generated within the slope are identified from continuous monitoring data. The microseismic event identification criteria are as follows:

[0016] (1) Signal-to-noise ratio criterion:

[0017]

[0018] Where N S x represents the count of data locations within the signal search window. Sk x represents the amplitude of the cross-correlation function at position k within the signal search window. tkm and x lkm Let N represent the cross-correlation function amplitudes of the leading noise window and the trailing noise window corresponding to the signal window at position k, respectively, while N... n This indicates the number of data points contained in the leading or trailing noise window, with a duration that is n times the size of the signal window. `max(k=1,…Ns)` identifies the peak amplitude ratio by examining all positions within the entire signal search window for a specific noise window duration setting, while `min(n=1,…30)` extracts the minimum value from these peak amplitude ratios calculated for different noise window durations.

[0019] (2) Power spectral density criterion: The average PSD value of the signal window exceeds the PSD value of the noise window by more than 1.5 times;

[0020] (3) Spatial consistency criterion: it must be detected at least at 3 stations and have a time difference consistent with the source-station distance.

[0021] S2.2: Accurately extract the observation travel time difference of the identified event between each station pair through waveform cross-correlation analysis;

[0022] S3: 3D positioning based on matching theoretical and observed travel time differences. A grid search algorithm is used to find the location in 3D space that minimizes the residual between the theoretical and observed travel time differences, thus determining the precise 3D coordinates of the microseismic event. The grid search positioning employs a weighted objective residual function:

[0023]

[0024] Where, ΔT theo (x,y,z,S i ,S j ) represents the theoretical connection between grid nodes (x, y, z) and stations in the time difference database (S). i ,S j Theoretical time difference, ΔT obs (S i ,S jThe weight w represents the actual time difference observed at this station, obtained through cross-correlation analysis. ij From the following formula, we can obtain:

[0025]

[0026] Among them, GAP ij SNR is the azimuth difference between the station and the relative source. ij CC represents the average signal-to-noise ratio of the station pair. ij The maximum cross-correlation coefficient is given by α, β, and γ, which are weighting coefficients whose optimal values ​​are determined through experience or experimentation.

[0027] A two-stage optimization strategy is adopted: The first stage is to divide the entire monitoring area into coarse grids with a large size (about 10 m) and perform a global search to identify the global minimum value region of the target residual function; The second stage is to divide the minimum value region into fine grids with a small size (about 2 m) and perform a local fine search. When the target residual function value is less than 0.1 seconds and the target residual function values ​​of the 8 surrounding adjacent nodes are all greater than the center point, it is determined as the final positioning result.

[0028] In constructing the theoretical travel time model, this invention innovatively employs environmental noise surface wave tomography to obtain high-resolution velocity structures. By deploying a dense array of temporary seismic stations in the slope area, surface wave signals are extracted using the cross-correlation function of environmental noise recorded by the stations, and then the three-dimensional distribution of underground S-wave velocity is obtained through inversion. This method eliminates the need for artificial seismic sources, is economical and practical, and can obtain a high-resolution, detailed velocity model. Simultaneously, this invention introduces a digital elevation model into ray tracing calculations, employing a fast travel method to calculate the elastic wave propagation path in a terrain-constrained three-dimensional space, ensuring that the obtained ray path conforms to physical laws and accurately reflects the influence of terrain on wave propagation.

[0029] In microseismic event identification and travel time extraction, this invention employs a dual criterion to ensure the reliability of the identification. By comprehensively analyzing the signal-to-noise ratio and power spectral density characteristics, and requiring consistency detection across multiple stations, various interference signals are effectively eliminated. For identified valid events, waveform cross-correlation technology is used to accurately measure the travel time difference between stations, with measurement accuracy reaching the millisecond level.

[0030] In 3D positioning, this invention defines a weighted residual function that comprehensively considers multiple factors and employs a two-stage grid search strategy to improve positioning efficiency. A coarse search quickly determines the target area, and a fine search obtains the precise location, achieving a balance between computational efficiency and positioning accuracy.

[0031] This invention applies environmental noise surface wave tomography to slope velocity modeling, enabling the acquisition of high-resolution velocity structures without the need for artificial seismic sources; it achieves three-dimensional ray tracing under the constraints of a digital elevation model, ensuring that wave propagation paths conform to physical laws; and it establishes a complete technical process, systematically solving the problem of microseismic location in complex slope environments.

[0032] The invention has significant advantages: the positioning accuracy reaches within 3 m, and can reach 1.5 m in the dense array coverage area, which is an order of magnitude higher than the traditional method; it has been successfully applied to complex environments with terrain undulations of more than 700 m and extreme velocity changes, proving the reliability of the method; it reduces monitoring costs by eliminating the need for artificial seismic sources such as blasting based on environmental noise; and the method has good universality and can be extended to microseismic monitoring in various complex geological environments. Attached Figure Description

[0033] Figure 1 This is a flowchart illustrating the overall process of the present invention, showing the complete flow from data acquisition to final positioning;

[0034] Figure 2 The diagram shows the spatial distribution of short-term dense arrays of stations and long-term monitoring stations.

[0035] Figure 3 The flowchart shows the surface wave tomography processing flow, including (a) the original environmental noise record; (b) the empirical Green's function extracted from the cross-correlation; (c) the dispersion curve extraction; and (d) the velocity-depth profile inversion results.

[0036] Figure 4 This is a diagram illustrating the principle of terrain-constrained ray tracing, showcasing the ray path characteristics under complex strata and terrain conditions.

[0037] Figure 5 The flowchart for microseismic event identification and travel time difference extraction includes: (a) waveform and SNR calculation from multiple stations; (b) PSD analysis of time spectrum plots; and (c) cross-correlation travel time difference extraction.

[0038] Figure 6 The spatial distribution map of the microseismic event location results in the study area shows the spatial distribution characteristics of 1470 microseismic events during the one-year monitoring period. Detailed Implementation

[0039] In this embodiment, the Shangquesuo (SQS) and Mindu (MD) creep slopes in the Jinsha River basin of the Hengduan Mountains in eastern Tibet Autonomous Region were selected as the test site. The terrain of this area has an elevation change of more than 700 m and a complex geological structure, making it an ideal location to test the method of the present invention.

[0040] This embodiment was implemented on a slope in the reservoir area of ​​a large hydropower station in eastern Tibet Autonomous Region. The slope has an elevation difference of approximately 700 m, multiple deformable bodies, and complex geological conditions, making it an ideal site for verifying the effectiveness of the method of this invention.

[0041] During the velocity model construction phase, 150 temporary seismic stations were deployed on two main deformable bodies, with station spacing ranging from 50 to 100 meters depending on the terrain conditions. Each station was equipped with a three-component seismic detector, sampling at 250 Hz, and continuously recorded 96 hours of environmental noise data. Data processing began with conventional preprocessing, including instrument response removal, mean removal, linear trend removal, time-domain normalization, and spectral whitening. Then, the cross-correlation function of all station pairs was calculated, and a stable empirical Green's function was obtained by superimposing 72 hours of data. Rayleigh wave dispersion curves, with a frequency range of 1-15 Hz, were extracted from the cross-correlation function using multichannel surface wave analysis. The dispersion curves were converted into S-wave velocity-depth profiles through nonlinear inversion, reaching a depth of 150 meters. Finally, the one-dimensional velocity profiles of all paths were combined into a three-dimensional velocity model using natural neighborhood interpolation. The inversion results show that the velocity of the loose surface deposits is 200-500 m / s, the velocity of the weathered and fractured zone is 500-1500 m / s, and the velocity of the relatively intact bedrock exceeds 2000 m / s, which is in good agreement with the stratigraphic structure revealed by the borehole.

[0042] In the ray tracing calculations, digital elevation model data with a resolution of 12.5 m was first acquired and matched to a 2 m computational grid using bilinear interpolation. The established three-dimensional computational domain has a horizontal range of approximately 3 × 6 km and a vertical depth of 150 m. The fast travel method was used to calculate the travel time from each grid node to all monitoring stations, with strict constraints that the ray path remain below the terrain surface. Comparison revealed a significant difference between the ray path and the straight path after considering terrain constraints, especially near valleys and ridges, where the path deviation could reach tens of meters, corresponding to travel time differences exceeding 0.1 seconds.

[0043] During the long-term monitoring phase, 15 long-term stations were deployed, forming a good spatial envelope around the main deformable body. Over the one-year monitoring period, the system automatically identified and located 1470 microseismic events. Event identification employed a short-window averaging to long-window averaging ratio method for initial triggering, followed by confirmation through signal-to-noise ratio and power spectral density analysis. For travel time extraction, a bandpass filter was designed based on the dominant frequency characteristics of the events, the envelope function was extracted using Hilbert transform, and the accurate travel time was obtained through normalized cross-correlation calculation.

[0044] The location calculation employed a two-stage strategy. The first stage involved searching the entire monitoring area with a 10-meter grid spacing, calculating the residual function value for each grid point to quickly determine the global minimum region of the residual function. The second stage involved a finer search within a 100×100×50-meter area at 2-meter intervals. Weighting coefficients were optimized through 10 artificial blasting experiments, and the final weights used were: azimuth coverage weight 0.4, signal-to-noise ratio weight 0.3, and cross-correlation coefficient weight 0.3. The final epicenter location was determined when the target residual function value was less than 0.1 seconds and was a local minimum.

[0045] Verified through manual blasting, the average positioning error of this method is 3 m, and within the core area covered by the dense array, the positioning error is less than 1.5 m. The positioning results show that microseismic events are mainly distributed within a depth range of 0-43 m, spatially highly concentrated in the surface deformation zones identified by InSAR, confirming the reliability of the positioning results. Compared with traditional methods, the positioning accuracy is improved by an order of magnitude, fully meeting the requirements for fine-grained slope stability evaluation.

[0046] This invention achieves high-precision three-dimensional positioning of microseismic events in complex slope environments by accurately considering complex geological structures and topographic effects, providing strong technical support for slope safety monitoring. This method is not only applicable to slopes in hydropower projects, but can also be extended to slope monitoring in mining, transportation, and land resources sectors, showing broad application prospects.

[0047] When applying the method described in this invention, the following should be noted: the dense array should cover the main area of ​​interest; the environmental noise collection time should be no less than 48 hours; long-term stations should form a good spatial envelope; and the system should be calibrated regularly through blasting experiments.

[0048] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A three-dimensional localization method for slope microseismic events considering complex geological structures and topographic effects, characterized in that, The steps include the following: S1: Construct a theoretical travel-time model adapted to complex slope environments, including: S1.1: Construct a high-resolution three-dimensional velocity model reflecting the complex geological structure of the slope using environmental noise surface wave tomography with a dense seismic array; S1.2: Based on the constraints of the digital elevation model (DEM), the fast travel method is used for ray tracing to calculate the minimum travel time path of seismic waves considering topographic effects; S1.3: Integrating velocity models and terrain-constrained paths, a theoretical travel time database from each potential seismic source location to the monitoring station is established; S2: Extract the actual travel time difference from the slope monitoring data, including: S2.1: The effective microseismic events generated inside the slope are identified from continuous monitoring data using both signal-to-noise ratio (SNR) and power spectral density (PSD) criteria. S2.2: Accurately extract the observation travel time difference of the identified event between each station pair through waveform cross-correlation analysis; S3: Based on the matching of theoretical and observed travel time differences, the three-dimensional positioning method uses a grid search algorithm to find the position in three-dimensional space that minimizes the residual between the theoretical and observed travel time differences, thereby determining the precise three-dimensional coordinates of the microseismic event.

2. The three-dimensional location method for slope microseismic events considering complex geological structures and topographic effects according to claim 1, characterized in that, The construction of the high-resolution three-dimensional velocity model in step S1.1 specifically includes: A dense array of seismic stations was deployed in the slope deformation area to collect more than 72 hours of environmental noise data. Preprocessing of noisy data includes instrument response removal, time-domain normalization, and spectral whitening; Calculate the cross-correlation function between all station pairs, and extract the empirical Green's function by phase-weighted superposition; The Rayleigh wave dispersion curve was extracted from the cross-correlation function by using the multichannel surface wave analysis method MASW combined with high-resolution Radon transform. The dispersion curve is converted into an S-wave velocity-depth profile by damped least squares inversion. Natural neighborhood interpolation is used to construct a three-dimensional velocity model from the one-dimensional velocity profiles of each path.

3. The three-dimensional location method for slope microseismic events considering complex geological structures and topographic effects according to claim 1, characterized in that, The ray tracing in step S1.2 specifically includes: Establish a three-dimensional computational domain that includes the slope topography; Use DEM data as the upper boundary of the computational domain and exclude grid nodes above the terrain surface; In a terrain-constrained 3D velocity field, the fast travel method is applied to calculate the travel time field from each grid node to all stations; the travel time ray path is determined by tracing in reverse along the negative travel time gradient direction, thus determining the minimum travel time ray path that follows the terrain constraints and velocity structure.

4. The three-dimensional location method for slope microseismic events considering complex geological structures and topographic effects according to claim 1, characterized in that, The theoretical timekeeping calculation in step S1.3 uses the following formula: ; Where T(x,y,z,S) is the theoretical travel time from the assumed source location (x,y,z) to station S; P is the shortest path; (i,j) are a pair of adjacent grid nodes on the shortest path; d i,j v is the Euclidean distance between adjacent nodes i and j; i and v j Let be the wave velocity values ​​at nodes i and j; Π be the set of all possible propagation paths; |∇ v| (i,j) β is the velocity gradient magnitude between nodes i and j; β is the velocity gradient correction coefficient, used to account for the impact of velocity changes on the path; α is the weighting factor, adjusting the influence of the velocity difference term; cosθ ij It is the cosine of the angle between the ray direction and the velocity transition direction, used to correct the propagation effect in non-uniform media.

5. The three-dimensional location method for slope microseismic events considering complex geological structures and topographic effects according to claim 1, characterized in that, The microseismic event identification criteria in step S2.1 are as follows: Signal-to-noise ratio criterion: ; in, This indicates the number of data locations present within the signal search window. This represents the amplitude of the cross-correlation function at position k within the signal search window. and Let represent the amplitudes of the cross-correlation functions of the leading noise window and the trailing noise window corresponding to the signal window at position k, respectively. This indicates the number of data points contained in the leading or trailing noise window, with a duration that is n times the size of the signal window. The maximum value is the peak amplitude ratio that the process identifies by examining all positions in the entire signal search window under a specific noise window duration setting, while the minimum value is the minimum value extracted from these peak amplitude ratios calculated by the process for different noise window durations.

6. The three-dimensional location method for slope microseismic events considering complex geological structures and topographic effects according to claim 5, characterized in that, The power spectral density criterion is that the average PSD value of the signal window exceeds the PSD value of the noise window by 1.5 times.

7. The three-dimensional localization method for slope microseismic events considering complex geological structures and topographic effects according to claim 1, characterized in that, The grid search localization in step S3 uses a weighted target residual function: ; in, For the theoretical to time difference database grid node (x,y,z) to station pair The theory of time difference, The weights represent the actual time difference observed at this station, obtained through cross-correlation analysis. From the following formula, we can obtain: ; Among them, GAP pq SNR is the azimuth difference between the station and the relative source. pq CC represents the average signal-to-noise ratio of the station pair. pq To maximize the cross-correlation coefficient, δ, θ, γ These are weighting coefficients, and their optimal values ​​are determined through experience or experimentation.

8. The three-dimensional location method for slope microseismic events considering complex geological structures and topographic effects according to claim 7, characterized in that, A two-stage optimization strategy is adopted, including Phase 1: Divide the entire monitoring area into large-scale coarse grids, perform a global search, and identify the global minimum region of the target residual function; The second stage involves dividing the minimum value region into a fine grid with a small size and performing a local fine search. When the target residual function value is less than 0.1 seconds and the target residual function values ​​of the eight surrounding adjacent nodes are all greater than the center point, the result is determined as the final location result. At this time, the result has clear spatial convergence characteristics.

Images