Tailing pond monitoring method based on passive source earthquake joint inversion

Abstract

The invention discloses a tailing pond monitoring method based on passive source seismic joint inversion. The method comprises the following steps: step 1, carrying out passive source data observation; 2, performing signal demodulation on the passive source data to obtain time-phased long-time continuous background noise and micro-seismic surface wave data, and performing preprocessing; 3, after preprocessing, multi-channel cross-correlation and frequency dispersion extraction are carried out, then two-step surface wave inversion imaging is carried out, and underground space shear wave velocity structure models in different periods are obtained; 4, performing micro-seismic effective event inversion by using a double-difference positioning algorithm to obtain spatial distribution of micro-seismic events in the underground space rock mass; 5, time-frequency analysis is carried out on the background noise and the micro-seismic surface wave data, and a time corresponding relation between time-frequency parameters and underground space internal evolution is established; and step 6, forming a tailing pond underground space three-dimensional dynamic reconstruction sensing system, and realizing the functions of real-time monitoring of potential geological catastrophe in the whole life cycle of the tailing pond and risk early warning.

Description

Technical Field

[0001] The present invention relates to the field of geophysical monitoring technology, in particular to a tailings pond monitoring method based on passive source seismic joint inversion. Background Art

[0002] Traditional underground monitoring technologies (such as electromagnetic methods and borehole exploration) have problems such as high cost, low resolution, and difficulty in long-term dynamic monitoring. Although existing DAS technology can achieve high-density observation, it lacks the ability to conduct multi-parameter joint inversion and real-time early warning for complex scenarios. For example, the identification of small-scale anomalies under strong noise interference in urban underground spaces and the lack of accuracy in the rapid positioning of microseismic events in mines have not been effectively solved. Microseismic events refer to earthquake events with extremely small magnitudes, usually between -3 and 2 on the Richter scale. The energy is so weak that humans cannot perceive it directly and must be detected with the help of highly sensitive instruments (such as seismographs, downhole detectors, etc.). Its occurrence may be caused by natural factors (such as tiny fractures in the crust, fluid migration) or human activities (such as hydraulic fracturing, mining, and geothermal development). Summary of the Invention

[0003] Purpose of the invention: The technical problem to be solved by the present invention is to provide a tailings pond monitoring method based on passive source seismic joint inversion in response to the shortcomings of the existing technology.

[0004] In order to solve the above technical problems, the present invention discloses a tailings pond monitoring method based on passive source seismic joint inversion, comprising the following steps:

[0005] Step 1: Conduct long-term continuous passive source data observation of the tailings pond using distributed fiber optic seismic sensors;

[0006] Step 2: demodulate the passive source data to obtain long-term continuous background noise and microseismic surface wave data in different time periods and perform preprocessing;

[0007] Step 3: Perform multi-channel cross-correlation and dispersion extraction on the pre-processed background noise surface wave data, and then perform two-step surface wave inversion imaging to obtain the shear wave velocity structure model of the underground space at different periods;

[0008] Step 4: Based on this velocity structure model, the double-difference positioning algorithm is used to perform microseismic effective event inversion to obtain the spatial distribution of microseismic events inside the underground rock mass;

[0009] Step 5: Use the generalized S transform to perform time-frequency analysis on the background noise and microseismic surface wave data obtained in step 2, and establish the time correspondence between the time-frequency parameters and the internal evolution of the underground space;

[0010] Step 6: Based on the shear wave velocity structure model obtained in step 3, the spatial distribution of microseismic events obtained in step 4, and the time correspondence obtained in step 5, a three-dimensional dynamic reconstruction perception system of the underground space of the tailings pond based on multi-source information is formed to realize the real-time monitoring and risk warning functions of potential geological disasters throughout the life cycle of the tailings pond.

[0011] The dynamic monitoring system constructed based on the three-dimensional dynamic reconstruction perception system of the tailings pond underground space includes: optical fiber seismic sensors and distributed optical fiber acoustic wave sensor demodulators buried on the surface of the tailings pond; optical fiber seismic sensors are buried on the surface of the tailings pond to conduct long-term passive source seismic data acquisition, extract high-density observation background noise and microseismic information respectively, and jointly invert to obtain the three-dimensional velocity structure evolution and spatiotemporal distribution law of microseismic events in the underground space at different periods, and at the same time apply generalized S transform to obtain signal time-frequency parameters, and thereby construct a more refined multi-parameter digital transparent display system for underground space information).

[0012] The fiber-optic seismic sensors were deployed shallowly in the ground with soil backfill, using a three-dimensional arrangement with a specific trace spacing, gauge length, and sampling rate. The optical fiber was shallowly buried in a trench approximately 30 cm long and backfilled with soil to lay out a 500-meter fiber-optic seismic sensor. During data acquisition, to meet the requirements of high-density small-scale detection, the trace spacing was set to 0.5 m, the gauge length was set to 0.5 m, and the sampling rate was 2000 Hz. The three-dimensional arrangement allowed for long-term continuous observation, with a single imaging data acquisition period of at least 5 hours. The MS-DAS2000 distributed fiber-optic acoustic sensor demodulator was planned to be used.

[0013] In step 2, the preprocessing of the background noise surface wave data includes: removing instrument response, removing median, bandpass filtering, continuous record cutting, time domain regularization, and spectrum whitening;

[0014] The preprocessing of the microseismic data includes filtering, denoising and energy recovery.

[0015] Step 3 is as follows:

[0016] Step 3-1, perform multi-channel cross-correlation processing on the pre-processed background noise surface wave data to obtain the cross-correlation dispersion energy spectrum function:

[0017] Assuming that the noise sources are randomly and uniformly distributed and the noise field is discrete, in a horizontally layered isotropic medium, the frequency domain cross-correlation function C(r,ω) of continuous background noise records of two x1 and x2 radial components at any distance r is proportional to the Green's function G:

[0018]

[0019] Where A is the amplitude constant, which is used to adjust the dimension or match the amplitude of the actual observed data. Im[·] is the imaginary part of the complex number, which extracts the part of the complex function related to energy dissipation, phase delay or attenuation. ω represents the angular frequency in radians per second, which represents the oscillation frequency of the signal. represents the radial component of the Green's function;

[0020] Perform FJ transformation on the frequency domain cross-correlation function to obtain the cross-correlation dispersion energy spectrum I(ω,k):

[0021]

[0022] Where J0(kr) is the zero-order Bessel function of the first kind, and k represents the wave number;

[0023] Substitute the frequency domain cross-correlation function C(r,ω) and the cross-correlation dispersion energy spectrum I(ω,k) into the radial Green's function

[0024]

[0025] Where g r (z=0,κ,ω) is the kernel function, which is independent of the distance r, z represents the depth, and κ represents the integral variable wave number; and according to the orthogonality of the zero-order Bessel function, the orthogonalized cross-correlation dispersion energy spectrum is obtained:

[0026]

[0027] Where δ(k-κ) is the Dirac delta function. The terms with k=κ in the integral are screened so that the integral result only retains the contribution of the target wavenumber k. When k=κ, the properties of the Dirichlet function and the relationship between the wavenumber k and the phase velocity c are as follows: The orthogonalized cross-correlation dispersion energy spectrum formula is further simplified as:

[0028] I(ω,k)=A·Im{g r (z = 0,κ,ω)};

[0029] Step 3-2, extracting the dispersion curve from the function obtained by multi-channel cross-correlation;

[0030] In step 3-3, the surface wave dispersion curve extracted in step 3-2 is used to perform a 1-D S-wave velocity structure inversion using the Computer Programs in Seismology (CPS330) program. A single survey line is interpolated to obtain a 2D velocity profile. Furthermore, a 3D velocity model, or shear wave velocity structure model of the subsurface space, is derived from the 2D velocity profiles of multiple survey lines. The survey line is a pre-designed straight path along which instruments (such as our fiber optic seismic sensors) are deployed to collect subsurface structure data.

[0031] The microseismic effective event inversion in step 4 is specifically as follows:

[0032] Step 4-1: Detect and identify effective microseismic events:

[0033] The short-to-long time average ratio method (STA / LTA) and the artificial intelligence event automatic detection method (PhaseNet) are used to pick the first arrival time of the S / P (transverse / longitudinal) wave group of the microseismic event (the first arrival time refers to the time when the seismic wave (such as P wave or S wave) first arrives at a certain observation point);

[0034] Step 4-2: Based on the underground space shear wave velocity structure model obtained by the background noise surface wave imaging obtained in step 3, the HypoDD double-difference positioning program is used to perform microseismic positioning inversion on effective events to obtain the spatial distribution of microseismic events inside the underground rock mass. The microseismic event inversion positioning tracks microseismic events induced by rock mass fractures such as karst collapse and delamination cracks.

[0035] The spatial distribution of microseismic events includes microseismic event density, which reflects the intensity and spatial distribution characteristics of underground rock fractures or fault activities. The microseismic event density is calculated by time density (the number of microseismic events occurring per unit time (such as the number of events / hour)), space density (the number of events per unit volume or area (such as the number of events / cubic meter or the number of events / square kilometer)), and joint space-time density; the density of the comprehensive space-time distribution (such as the number of events / (km 3 ·day)))measure.

[0036] Step 4-2 is as follows:

[0037] Step 4-2-1: For microseismic events i and j recorded at the same station (fiber optic seismic sensor location) k, the double difference is defined as

[0038]

[0039] Where t is the travel time of the micro-earthquake, is the observed value of microseismic travel time difference, is the theoretical value of microseismic travel time difference;

[0040] Step 4-2-2, the double difference equations of all events recorded by all stations are combined to obtain the reliable locations and occurrence times of all microseismic sources:

[0041] WGm=Wd

[0042] Where G is an M×4N partial differential matrix, M is the number of double-difference observations, and N is the number of microseismic events; m is the change in source parameters, d is the double-difference vector, and W is the diagonal weighted matrix. Considering the varying quality of data received at each receiving point and the fact that events have less influence on each other as they are farther apart, the diagonal weighted matrix W is added. By solving for the change in source parameters m and iteratively updating m, we gradually approach the optimal solution and ultimately obtain the reliable locations and onset times of all microseismic sources. This describes the spatiotemporal evolution of microfractures and fissures within the rock mass during the initial stages of underground anomaly formation and its development.

[0043] When the distance between the two sources is small enough compared to the epicentral distance and the scale of velocity inhomogeneity, the double difference in step 4-1 is expressed as:

[0044]

[0045] Among them, the source parameter disturbance transformation amount Δm, which is specifically (Δx, Δy, Δz, Δτ), Δx represents the position correction amount of the source in the horizontal direction (east-west direction) (unit: usually kilometers or meters), a positive value indicates correction to the east, and a negative value indicates correction to the west (depending on the coordinate system definition); Δy represents the position correction amount of the source in the horizontal direction (north-south direction) (unit: kilometers or meters), a positive value indicates correction to the north, and a negative value indicates correction to the south; Δz represents the position correction amount of the source in the vertical direction (depth) (unit: kilometers or meters), a positive value indicates correction to deeper crust, and a negative value indicates correction to shallower part; Δτ represents the correction amount of the origin time of the source (unit: seconds), a positive value indicates a delay in the origin time, and a negative value indicates an advance.

[0046] Then the above double difference formula is expanded into:

[0047]

[0048] The source parameter perturbation transformation, Δm, represents the correction to the initial source parameters of the earthquake event. Its purpose is to minimize the residual between theoretical and observed travel times through iterative optimization. The double-difference positioning method uses the residual between the observed and theoretical microseismic travel time differences to invert the relative position of each event in the microseismic cluster to the cluster center. This method can constrain the mutual microseismic events within the rock mass and reduce errors caused by inaccurate velocity models.

[0049] The transformation in step 5, generalized S transform (GST), is a non-stationary signal analysis and processing method between short-time Fourier transform (STFT) and wavelet transform (WT), and the formula is as follows:

[0050]

[0051] Where S(τ,f) is the time-frequency distribution function, h(t) is the original time domain signal, is the normalization coefficient, is the Gaussian window function, f is the analysis frequency, τ is the analysis time point, λ is the window function width control parameter, p is the frequency power adjustment parameter, e -i2πft is the basis function of Fourier transform.

[0052] Step 6 is as follows:

[0053] Based on the velocity change gradient obtained from the underground space shear wave velocity structure model in step 3, the microseismic event density in step 4 and the spectrum resonance warning threshold interval matrix in step 5 (a tool for monitoring and warning system state changes based on spectrum analysis theory (signal processing, converting time domain signals into frequency domain signals, identifying the energy distribution of different frequency components) and resonance phenomenon), the possibility and time of instability of geological anomalies are determined to realize the function of real-time monitoring and risk warning of potential geological disasters throughout the life cycle of the tailings pond.

[0054] This velocity structure model provides important parameters for identifying and analyzing underground spatial distribution patterns, constructing digital, transparent displays of geological information, and real-time dynamic evolution monitoring. It also provides a reliable initial velocity model for microseismic positioning. (Distributed fiber-optic seismic sensing uses background noise surface waves to invert the S-wave velocity structure of the underground space, identifying and analyzing underground spatial distribution patterns, and enabling digital, transparent displays of underground geological information.)

[0055] The microseismic event density reveals the spatiotemporal evolution of cracks and internal deformation at various stages of disasters, such as karst collapse and delamination fissures, providing a scientific basis for revealing the mechanisms of multi-field coupled geological disasters in underground spaces and for monitoring and forecasting. (Distributed fiber-optic seismic sensing microseismic monitoring locates microfracture events induced by deformation of surrounding rock masses in underground spaces in real time, characterizing the spatiotemporal evolution of micro-, fine-, and macro-stratum damage and healing, and building a system for identifying and forecasting geological disaster precursor information.)

[0056] The time correspondence is achieved by identifying the dominant resonant frequency from the time-frequency parameters of the noise data and analyzing the time-frequency characteristics of rock microfracture signals at different deformation stages. This establishes a temporal relationship between the time-frequency parameters and the damage to underground structures, providing a reference for analyzing underground stability. (Time-frequency analysis of distributed fiber-optic seismic sensing signals establishes a temporal relationship between the time-frequency parameters and the damage to underground structures.)

[0057] Beneficial effects: Distributed fiber optic seismic sensing technology has the advantage of high spatial resolution. It can be deployed over long distances on the surface in complex environments, achieve long-term ultra-high-density observations at a relatively low cost, realize fine structural imaging of underground space, and reveal small-scale structural anomalies and disaster mechanisms.

[0058] Microseismic positioning monitoring can obtain real-time information on microscopic fractures of rock masses inside underground spaces, and directly monitor the initial form of abnormal structures such as karst and ground fissures and the entire process of internal deformation evolution.

[0059] Establish a spatial ultra-high-density passive source seismic joint imaging process using distributed fiber-optic seismic sensing technology to meet the needs of fine structure surveys at smaller scales underground.

[0060] Universal high-precision underground space structure evolution perception and real-time monitoring technology can achieve three-dimensional transparent display of multi-parameter geological information throughout the entire life cycle of underground space, and characterize the spatiotemporal evolution characteristics of internal deformation of underground space, providing technical support for multi-scale complex underground space resource utilization evaluation, geological disaster formation mechanism and monitoring and forecasting.

[0061] Distributed fiber-optic seismic sensing technology itself has many advantages. ① It uses the optical fiber itself as the sensing unit, which can continuously sample the seismic wave field along the optical fiber. ② The optical fiber itself is resistant to extreme temperatures and strong electromagnetic interference, which greatly reduces operation and maintenance costs and is also conducive to construction in harsh environments. ③ Not only can specially laid special sensor cables be used, but also existing sensor cables can be used, thereby reducing field deployment costs. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 Flowchart of the present invention.

[0063] Figure 2 is the background noise Rayleigh wave phase velocity dispersion energy spectrum

[0064] Figure 3 The shear wave velocity model is inverted from the dispersion curve.

[0065] Figure 4 Time-frequency analysis of passive source earthquake records.

[0066] Figure 5 (a) is the three-dimensional velocity structure of the underground surrounding rock mass fracture, and (b) is a slice of the spatiotemporal evolution of the microfracture event. DETAILED DESCRIPTION

[0067] A tailings pond monitoring method based on passive source seismic joint inversion includes the following steps:

[0068] Step 1: Conduct long-term continuous passive source data observation of the tailings pond using distributed fiber optic seismic sensors;

[0069] Step 2: demodulate the passive source data to obtain long-term continuous background noise and microseismic surface wave data in different time periods and perform preprocessing;

[0070] Step 3: Perform multi-channel cross-correlation and dispersion extraction on the pre-processed background noise surface wave data, and then perform two-step surface wave inversion imaging to obtain the shear wave velocity structure model of the underground space at different periods;

[0071] Step 4: Based on this velocity structure model, the double-difference positioning algorithm is used to perform microseismic effective event inversion to obtain the spatial distribution of microseismic events inside the underground rock mass;

[0072] Step 5: Use the generalized S transform to perform time-frequency analysis on the background noise and microseismic surface wave data obtained in step 2, and establish the time correspondence between the time-frequency parameters and the internal evolution of the underground space;

[0073] Step 6: Based on the shear wave velocity structure model obtained in step 3, the spatial distribution of microseismic events obtained in step 4, and the time correspondence obtained in step 5, a three-dimensional dynamic reconstruction perception system of the underground space of the tailings pond based on multi-source information is formed to realize the real-time monitoring and risk warning functions of potential geological disasters throughout the life cycle of the tailings pond.

[0074] In step 2, the preprocessing of the background noise surface wave data includes: removing instrument response, removing median, bandpass filtering, continuous record cutting, time domain regularization, and spectrum whitening;

[0075] The preprocessing of the microseismic data includes filtering, denoising and energy recovery.

[0076] Step 3 is as follows:

[0077] Step 3-1, perform multi-channel cross-correlation processing on the pre-processed background noise surface wave data to obtain the cross-correlation dispersion energy spectrum function:

[0078] Assuming that the noise sources are randomly and uniformly distributed and the noise field is discrete, in a horizontally layered isotropic medium, the frequency domain cross-correlation function C(r,ω) of continuous background noise records of two x1 and x2 radial components at any distance r is proportional to the Green's function G:

[0079]

[0080] Where A is the amplitude constant, Im[·] is the imaginary part of the complex number, which extracts the part of the complex function related to energy dissipation, phase delay or attenuation, and ω represents the angular frequency in radians per second, which represents the oscillation frequency of the signal. represents the radial component of the Green's function;

[0081] Perform FJ transformation on the frequency domain cross-correlation function to obtain the cross-correlation dispersion energy spectrum I(ω,k):

[0082]

[0083] Where J0(kr) is the zero-order Bessel function of the first kind, and k represents the wave number;

[0084] Substitute the frequency domain cross-correlation function C(r,ω) and the cross-correlation dispersion energy spectrum I(ω,k) into the radial Green's function

[0085]

[0086] Where g r (z=0,κ,ω) is the kernel function, which is independent of the distance r, z represents the depth, and κ represents the integral variable wave number; and according to the orthogonality of the zero-order Bessel function, the orthogonalized cross-correlation dispersion energy spectrum is obtained:

[0087]

[0088] Where δ(k-κ) is the Dirac delta function. The terms with k=κ in the integral are screened so that the integral result only retains the contribution of the target wavenumber k. When k=κ, the properties of the Dirichlet function and the relationship between the wavenumber k and the phase velocity c are as follows: The orthogonalized cross-correlation dispersion energy spectrum formula is further simplified as:

[0089] I(ω,k)=A·Im{g r (z = 0,κ,ω)};

[0090] Step 3-2, extracting the dispersion curve from the function obtained by multi-channel cross-correlation;

[0091] In step 3-3, a 1-D S-wave velocity structure inversion is performed on the surface wave dispersion curve extracted in step 3-2, and a 2D velocity profile is obtained by interpolation of a single survey line. Then, a 3D velocity model, i.e., a shear wave velocity structure model of the underground space, is obtained from the 2D velocity profiles of multiple survey lines.

[0092] The microseismic effective event inversion in step 4 is specifically as follows:

[0093] Step 4-1: Detect and identify effective microseismic events:

[0094] The short-to-long time average ratio method and artificial intelligence event automatic detection method are used to pick the first arrival time of the S / P wave group of microseismic events;

[0095] Step 4-2: Based on the underground space shear wave velocity structure model obtained by background noise surface wave imaging obtained in step 3, the HypoDD double-difference positioning program is used to perform microseismic positioning inversion on effective events to obtain the spatial distribution of microseismic events inside the underground space rock mass.

[0096] The spatial distribution of microseismic events includes the density of microseismic events, which is measured by time density, space density and joint space-time density.

[0097] Step 4-2 is as follows:

[0098] Step 4-2-1: For microseismic events i and j recorded at the same station k, the double difference is defined as

[0099]

[0100] Where t is the travel time of the micro-earthquake, is the observed value of microseismic travel time difference, is the theoretical value of microseismic travel time difference;

[0101] Step 4-2-2, the double difference equations of all events recorded by all stations are combined to obtain the reliable locations and occurrence times of all microseismic sources:

[0102] WGm=Wd

[0103] Where G is an M×4N partial differential matrix, M is the number of double-difference observations, and N is the number of microseismic events. m is the change in source parameters, d is the double-difference vector, and W is the diagonal weighted matrix. By solving the change in source parameters m and iteratively updating m, we gradually approach the optimal solution and finally obtain the reliable locations and occurrence times of all microseismic sources.

[0104] When the distance between the two sources is small enough compared to the epicentral distance and the scale of velocity inhomogeneity, the double difference in step 4-1 is expressed as:

[0105]

[0106] Among them, the source parameter disturbance transformation Δm is specifically (Δx, Δy, Δz, Δτ), and the above double difference formula is expanded into:

[0107]

[0108] The generalized S transform in step 5 is a method for analyzing and processing non-stationary signals between the short-time Fourier transform and the wavelet transform. The formula is as follows:

[0109]

[0110] Where S(τ,f) is the time-frequency distribution function, h(t) is the original time domain signal, is the normalization coefficient, is the Gaussian window function, f is the analysis frequency, τ is the analysis time point, λ is the window function width control parameter, p is the frequency power adjustment parameter, e -i2πft is the basis function of Fourier transform.

[0111] Step 6 is as follows:

[0112] Based on the velocity gradient obtained from the underground space shear wave velocity structure model in step 3, the microseismic event density in step 4, and the spectrum resonance warning threshold interval matrix in step 5, the possibility and time of geological anomaly instability are determined to achieve the function of real-time monitoring and risk warning of potential geological disasters throughout the life cycle of the tailings pond.

[0113] This example corresponds to a real tailings pond monitoring:

[0114] The dynamic monitoring system constructed based on the three-dimensional dynamic reconstruction perception system of the tailings pond underground space includes: fiber optic seismic sensors buried on the tailings pond surface and distributed fiber optic acoustic wave sensor demodulators. (The fiber optic seismic sensors buried on the tailings pond surface perform long-term passive source seismic data acquisition, extract high-density observation background noise and microseismic information, and jointly invert to obtain the three-dimensional velocity structure evolution and spatiotemporal distribution of microseismic events in the underground space at different periods. At the same time, the generalized S transform is applied to obtain the signal time-frequency parameters, thereby constructing a more refined multi-parameter digital transparent display system for underground space information.)

[0115] The fiber optic seismic sensors described in step 1 were deployed in three dimensions using shallow ground burial and backfilling, with a specific trace spacing, gauge length, and sampling rate. The optical fiber was shallowly buried in a trench approximately 30 cm long and backfilled with soil to deploy a 500-meter fiber optic seismic sensor. During data acquisition, to meet the requirements of high-density, small-scale detection, the trace spacing was set to 0.5 m, the gauge length was 0.5 m, and the sampling rate was 2000 Hz. Long-term continuous observation was performed, with a single imaging data acquisition period of at least 5 hours. The MS-DAS2000 distributed fiber optic acoustic sensor demodulator was used.

[0116] like Figure 2 After preprocessing in step 2, the cross-correlation dispersion energy spectrum is calculated, and the surface wave dispersion curve is picked up from the energy convergence point in the figure. The horizontal axis in the figure is frequency, the vertical axis is velocity, and the red area is the spectrum energy.

[0117] like Figure 3The horizontal axis is shear wave velocity, the vertical axis is depth, initial is the initial value, inversion is the inversion value, and mean is the mean. Based on the dispersion curve obtained in step 3, the Computer Programs in Seismology (CPS330) program is used to perform 1-D S-wave velocity structure inversion. A 2D velocity profile is obtained by interpolating the single survey line. A 3D velocity model, i.e., the shear wave velocity structure model of the underground space, is then obtained from the 2D velocity profiles of multiple survey lines. The survey line is a pre-designed straight path along which instruments (such as our fiber optic seismic sensors) are deployed to collect underground structure data.

[0118] In step 4, according to the formula WGm = Wd, taking into account the varying data quality at each receiving point and the fact that events have less influence on each other as they are farther apart, a diagonal weighted matrix W is added. By solving for the change in the source parameters m and iteratively updating m, the optimal solution is gradually approached, ultimately yielding the reliable locations and onset times of all microseismic sources. This describes the spatiotemporal evolution of microfractures and fissures within the rock mass during the initial stages of the underground anomaly's catastrophic formation and its development.

[0119] like Figure 4 The left figure shows the waveform of part of the earthquake record, with the horizontal axis being time and the vertical axis being frequency; the right figure shows the time-frequency analysis, with the horizontal axis being time and the vertical axis being Hertz. The time-frequency analysis results of the waveform recorded by a single station are obtained from step 5. The frequency and energy of the earthquake wave change with the deformation of the underground space. According to the formula The generalized S transform not only has multiscale focusing properties but also maintains the absolute phase of the frequency. The basic transformation function does not need to meet the admissibility condition, which compensates for the shortcomings of the short-time Fourier transform and wavelet transform, and has variable time-frequency resolution. The generalized S transform is used to perform time-frequency analysis on each continuous recording, and its corresponding relationship with the changing patterns of the underground space at different stages is statistically analyzed.

[0120] From step 6, the temporal and spatial distribution of the underground S-wave (transverse wave) velocity structure and microseismic positioning at each stage of the underground space of the tailings pond is obtained, and a three-dimensional transparent display of geological information and dynamic evolution monitoring are established. Figure 5 , which are part of the results of the previous physical simulation experiments on rock fracture in underground space. Figure 5 (a) The coordinates are the three-dimensional spatial position, red is low speed, blue is high speed, and the change of shear wave velocity can be observed. Figure 5 The three figures in (b) correspond to the microseismic density and shear wave velocity slices at the same location at different times. The coordinates are the three-dimensional spatial positions, the red and blue represent the differential changes in velocity, and the black dots are the density of microseismic events. The density changes of microseismic events can be observed through Figure 5 (a) Observe the change of speed and Figure 5 (b) The velocity differential rate of change over different periods enables real-time monitoring and risk warning of tailings pond geological disasters. A yellow alert is activated when the velocity differential rate exceeds 2–3 standard deviations of the baseline value (usually based on historical data or monitoring data within a specific time period (such as several days, weeks, or the initial stable stage of operations, calculated using statistical methods (such as mean, median, moving average, etc.). It represents the "normal state" of the system when there are no abnormal disturbances) or the microseismic event density increases by 30–50% (lasting ≥4 hours). A red alert is activated when the velocity differential rate exceeds 5 standard deviations of the baseline value or changes suddenly by 50%, or the microseismic event density increases by more than 80% (or has significant spatial clustering). A yellow alert is activated when the velocity differential rate exceeds 2 standard deviations of the baseline value and the microseismic event density increases by 40%. A red alert is activated when the velocity differential rate exceeds 3 standard deviations of the baseline value and the microseismic event density increases by 60%.

[0121] The present invention provides a tailings pond monitoring method based on passive source seismic joint inversion. There are many methods and approaches to implement this technical solution. The above is only a preferred embodiment of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention. These improvements and modifications should also be considered within the scope of protection of the present invention. Any components not specified in this embodiment can be implemented using existing technologies.

Claims

1. A tailings pond monitoring method based on passive source seismic joint inversion, characterized in that: The following steps are involved: Step 1: Conduct long-term continuous passive source data observation of the tailings pond using distributed fiber optic seismic sensors; Step 2: demodulate the passive source data to obtain long-term continuous background noise and microseismic surface wave data in different time periods and perform preprocessing; Step 3: Perform multi-channel cross-correlation and dispersion extraction on the pre-processed background noise surface wave data, and then perform two-step surface wave inversion imaging to obtain the shear wave velocity structure model of the underground space at different periods; Step 4: Based on this velocity structure model, the double-difference positioning algorithm is used to perform microseismic effective event inversion to obtain the spatial distribution of microseismic events inside the underground rock mass; Step 5: Use the generalized S transform to perform time-frequency analysis on the background noise and microseismic surface wave data obtained in step 2, and establish the time correspondence between the time-frequency parameters and the internal evolution of the underground space; Step 6: Based on the shear wave velocity structure model obtained in step 3, the spatial distribution of microseismic events obtained in step 4, and the time correspondence obtained in step 5, a three-dimensional dynamic reconstruction perception system of the underground space of the tailings pond based on multi-source information is formed to realize the real-time monitoring and risk warning functions of potential geological disasters throughout the life cycle of the tailings pond.

2. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 1 is characterized in that: The dynamic monitoring system constructed based on the three-dimensional dynamic reconstruction perception system of the underground space of the tailings pond includes: optical fiber seismic sensors buried on the surface of the tailings pond and distributed optical fiber acoustic wave sensor demodulators; the optical fiber seismic sensors are laid out in a shallow manner and backfilled with soil, and are deployed in three dimensions according to a certain track spacing, gauge length and sampling rate.

3. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 1 is characterized in that: In step 2, the preprocessing of the background noise surface wave data includes: removing instrument response, removing median, bandpass filtering, continuous record cutting, time domain regularization, and spectrum whitening; The preprocessing of the microseismic data includes filtering, denoising and energy recovery.

4. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 1 is characterized in that: Step 3 is as follows: Step 3-1, perform multi-channel cross-correlation processing on the pre-processed background noise surface wave data to obtain the cross-correlation dispersion energy spectrum function: Assuming that the noise sources are randomly and uniformly distributed and the noise field is discrete, in a horizontally layered isotropic medium, the frequency domain cross-correlation function C(r,ω) of continuous background noise records of two x1 and x2 radial components at any distance r is proportional to the Green's function G: Where A is the amplitude constant, Im[·] is the imaginary part of the complex number, which extracts the part of the complex function related to energy dissipation, phase delay or attenuation, and ω represents the angular frequency in radians per second, which represents the oscillation frequency of the signal. represents the radial component of the Green's function; Perform FJ transformation on the frequency domain cross-correlation function to obtain the cross-correlation dispersion energy spectrum I(ω,k): Where J0(kr) is the zero-order Bessel function of the first kind, and k represents the wave number; Substitute the frequency domain cross-correlation function C(r,ω) and the cross-correlation dispersion energy spectrum I(ω,k) into the radial Green's function Where g r (z=0,κ,ω) is the kernel function, which is independent of the distance r, z represents the depth, and κ represents the integral variable wave number; and according to the orthogonality of the zero-order Bessel function, the orthogonalized cross-correlation dispersion energy spectrum is obtained: Where δ(k-κ) is the Dirac delta function. The terms with k=κ in the integral are screened so that the integral result only retains the contribution of the target wavenumber k. When k=κ, the properties of the Dirichlet function and the relationship between the wavenumber k and the phase velocity c are as follows: The orthogonalized cross-correlation dispersion energy spectrum formula is further simplified as: I(ω,k)=A·Im{g r (z=0,κ,ω)}; Step 3-2, extracting the dispersion curve from the function obtained by multi-channel cross-correlation; In step 3-3, a 1-D S-wave velocity structure inversion is performed on the surface wave dispersion curve extracted in step 3-2, and a 2D velocity profile is obtained by interpolation of a single survey line. Then, a 3D velocity model, i.e., a shear wave velocity structure model of the underground space, is obtained from the 2D velocity profiles of multiple survey lines.

5. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 1 is characterized in that: The microseismic effective event inversion in step 4 is specifically as follows: Step 4-1: Detect and identify effective microseismic events: The short-to-long time average ratio method and artificial intelligence event automatic detection method are used to pick the first arrival time of the S / P wave group of microseismic events; Step 4-2: Based on the underground space shear wave velocity structure model obtained by background noise surface wave imaging obtained in step 3, the HypoDD double-difference positioning program is used to perform microseismic positioning inversion on effective events to obtain the spatial distribution of microseismic events inside the underground space rock mass.

6. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 5 is characterized in that: The spatial distribution of microseismic events includes microseismic event density, which is measured by time density, space density, and joint space-time density.

7. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 1 is characterized in that: Step 4-2 is as follows: Step 4-2-1: For microseismic events i and j recorded at the same station k, the double difference is defined as Where t is the travel time of the micro-earthquake, is the observed value of microseismic travel time difference, is the theoretical value of microseismic travel time difference; Step 4-2-2, the double difference equations of all events recorded by all stations are combined to obtain the reliable locations and occurrence times of all microseismic sources: WGm=Wd Where G is an M×4N partial differential matrix, M is the number of double-difference observations, and N is the number of microseismic events. m is the change in source parameters, d is the double-difference vector, and W is the diagonal weighted matrix. By solving the change in source parameters m and iteratively updating m, we gradually approach the optimal solution and finally obtain the reliable locations and occurrence times of all microseismic sources.

8. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 7 is characterized in that: When the distance between the two sources is small enough compared to the epicentral distance and the scale of velocity inhomogeneity, the double difference in step 4-1 is expressed as: Among them, the source parameter disturbance transformation Δm is specifically (Δx, Δy, Δz, Δτ), and the above double difference formula is expanded into:

9. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 1 is characterized in that: The generalized S transform in step 5 is a method for analyzing and processing non-stationary signals between the short-time Fourier transform and the wavelet transform. The formula is as follows: Where S(τ,f) is the time-frequency distribution function, h(t) is the original time domain signal, is the normalization coefficient, is the Gaussian window function, f is the analysis frequency, τ is the analysis time point, λ is the window function width control parameter, p is the frequency power adjustment parameter, e -i2πft is the basis function of Fourier transform.

10. The tailings pond monitoring method based on passive source seismic joint inversion according to claim 8, characterized in that: Step 6 is as follows: Based on the velocity change gradient obtained from the underground space shear wave velocity structure model in step 3, the microseismic event density in step 4, and the spectrum resonance warning threshold interval matrix in step 5, the possibility and time of instability of the geological anomaly body are determined, realizing the function of real-time monitoring and risk warning of potential geological disasters throughout the life cycle of the tailings pond.

Images