Intelligent rockburst dynamic risk prediction method fusing real-time microseismic data

Abstract

This invention relates to the field of safety monitoring and disaster early warning technology for deep underground engineering, and discloses an intelligent rockburst dynamic risk prediction method integrating real-time microseismic data. This method first constructs a high-precision medium model reflecting the real-time stress and attenuation characteristics of the rock mass based on background noise interferometric imaging and the adjoint state method inversion. Then, it reconstructs a four-dimensional energy density field, defines and calculates the waveguide locking index to accurately screen key energy localization regions. Subsequently, it extracts the energy residual sequence of the key regions, calculates the lag first-order autocorrelation coefficient and time-domain variance, and uses the critical slowing theory to determine the precursors of system dynamic instability and trigger graded early warnings. This invention corrects imaging biases caused by energy dissipation through physics-driven medium parameter inversion, and combined with data-driven nonlinear dynamic analysis, achieves accurate spatiotemporal positioning of the rockburst incubation core area and early identification of disaster risks.

Description

Technical Field

[0001] This invention relates to the field of safety monitoring and disaster early warning technology for deep underground engineering, specifically to an intelligent rockburst dynamic risk prediction method that integrates real-time microseismic data. Background Technology

[0002] As underground engineering projects extend into deeper strata, rockburst hazards induced by high ground stress environments have become a major bottleneck restricting engineering safety. Microseismic monitoring technology, as the mainstream method for assessing rock mass stability, captures elastic wave signals generated by rock mass fracturing to invert source parameters, providing data support for rockburst risk assessment.

[0003] However, existing microseismic monitoring technologies suffer from model applicability limitations when applied to complex deep geological environments. Traditional source location and wavefield analysis largely rely on pre-defined static velocity models and ideal elastic medium assumptions, neglecting the real-time dynamic changes in wave velocity fields caused by excavation disturbances and the inherent viscous attenuation effect of deep media. This deviation in the physical model leads to reduced resolution in reverse time migration imaging and makes it difficult to accurately reconstruct energy dissipation during wavefield propagation, resulting in distortions in source location accuracy and energy inversion results.

[0004] Furthermore, in terms of risk identification and early warning mechanisms, existing methods mainly rely on the spatial distribution density and energy statistical thresholds of discrete microseismic events for assessment. This approach struggles to quantify the continuous flow and local accumulation of energy within the rock mass, and is insufficient to effectively identify the energy waveguide structures that lead to disasters. Simultaneously, relying solely on threshold warnings based on numerical exceedances often exhibits a lag, failing to capture critical slowdown precursors of reduced system dynamic stability before macroscopic damage occurs in the rock mass, resulting in insufficient on-site evacuation response time. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides an intelligent rockburst dynamic risk prediction method that integrates real-time microseismic data. This method solves the problems of low imaging inversion accuracy, difficulty in accurately identifying waveguide structures with energy accumulation, and lag in traditional threshold early warning mechanisms caused by the reliance on static elastic models in existing microseismic monitoring technologies.

[0006] The first aspect of this invention provides an intelligent rockburst dynamic risk prediction method that integrates real-time microseismic data. It is applied to a system that includes a microseismic monitoring sensor array, a data transmission unit, and a central processing server. The aim is to achieve spatiotemporal quantification of the energy accumulation state of deep rock masses and early identification of disaster precursors through physical-driven medium parameter inversion and data-driven critical slowing characteristic analysis.

[0007] First, data splitting and medium velocity modeling are performed. The short-to-long-time average ratio algorithm is used to separate the acquired signal into a background noise data stream and a microseismic event data stream. Based on the principles of seismic interferometry, the empirical Green's function is recovered from the background noise data stream. Combined with the equation of process and tomographic algorithms, a three-dimensional heterogeneous velocity model reflecting the real-time stress state of the rock mass is inverted and constructed.

[0008] Building upon this, this invention employs the adjoint state method to invert attenuation parameters. A viscous acoustic wave equation is constructed, using a three-dimensional heterogeneous velocity model as initial values, and inverse-time wavefield extension is performed using microseismic event data streams. Wavefield focusing entropy is introduced as the objective function for imaging quality evaluation, and its gradient relative to the attenuation parameter model is calculated and iteratively updated until the optimal quality factor model is obtained, thereby correcting the energy dissipation effect caused by medium viscosity.

[0009] Based on the aforementioned dual model, this invention implements energy field reconstruction and waveguide structure screening. A full-waveform forward simulation is performed on the microseismic data, and kinetic and potential energy densities are superimposed to generate a four-dimensional energy density field. To quantitatively identify waveguide structures with energy localization characteristics, this invention defines and calculates the waveguide locking index, determined by the product of the time-integrated energy and the peak energy at the spatial grid points.

[0010] The grids with waveguide locking exponents exceeding a preset threshold are selected to form key regions, and the average energy density time series of the key regions formed by the grids is extracted.

[0011] For key regions, this invention further utilizes critical slowing theory to extract dynamic precursors. A detrending process is performed on the energy density sequence to obtain the residual sequence, and the sliding window technique is used to calculate the hysteresis first-order autocorrelation coefficient reflecting the system's memory effect and the time-domain variance reflecting the intensity of fluctuations.

[0012] Ultimately, this invention determines risk through a comprehensive early warning decision module. The aforementioned indicators are normalized and weighted to generate a comprehensive instability probability index. When the comprehensive instability probability index exceeds a preset yellow or red warning threshold and shows a monotonically increasing trend, it is determined that the rock mass system's resilience has decreased and it has entered a critical slowdown state, triggering a graded risk warning signal.

[0013] This invention provides an intelligent rockburst dynamic risk prediction method that integrates real-time microseismic data. It has the following beneficial effects: 1. This invention utilizes background noise interferometry and the adjoint state method to construct a real-time attenuation model, which not only avoids the risk of artificial seismic source excitation in deep engineering, but also corrects the wave field distortion caused by medium viscosity attenuation, thereby improving the reconstruction accuracy of microseismic energy reverse time migration imaging.

[0014] 2. This invention constructs a four-dimensional energy density field and waveguide locking index. By quantifying the long-term accumulation of energy and peak intensity characteristics, it effectively eliminates transient wave interference and accurately identifies the rockburst-generating core waveguide structure with energy trapping effect.

[0015] 3. Based on the critical slowing theory, this invention monitors the lag first-order autocorrelation coefficient and time-domain variance of the energy residual sequence, and keenly captures the precursors of dynamic instability of the system before macroscopic failure of the rock mass, thus solving the lag problem of traditional single energy threshold monitoring. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the intelligent rockburst dynamic risk prediction system of the present invention; Figure 2 This is a schematic diagram of the comprehensive early warning decision-making process of the present invention.

[0017] Among them, 10 is the microseismic monitoring sensor array; 20 is the data transmission unit; 30 is the central processing server; 301 is the data diversion and preprocessing module; 302 is the medium velocity model construction module; 303 is the attenuation parameter iterative inversion module; 304 is the energy field reconstruction and screening module; 305 is the dynamic characteristic analysis module; and 306 is the comprehensive early warning decision module. Detailed Implementation

[0018] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] See attached document Figure 1 , Figure 1 This is a schematic diagram of the structure of an intelligent rockburst dynamic risk prediction system that integrates real-time microseismic data according to an embodiment of the present invention. The system provided in this embodiment includes: a microseismic monitoring sensor array 10, a data transmission unit 20, and a central processing server 30.

[0020] The microseismic monitoring sensor array 10 is deployed inside or on the surface of the rock mass in the monitoring area. It consists of multiple spatially distributed acceleration or velocity sensors and is used to collect continuous mechanical vibration signals. The data transmission unit 20 connects the sensor array 10 to the central processing server 30 and is responsible for the analog-to-digital conversion and transmission of the signals.

[0021] The central processing server 30 includes a data offloading and preprocessing module 301, a medium velocity model construction module 302, an attenuation parameter iterative inversion module 303, an energy field reconstruction and screening module 304, a dynamic characteristic analysis module 305, and a comprehensive early warning decision module 306. These modules work collaboratively to perform data processing, inversion, and risk identification.

[0022] The intelligent rockburst dynamic risk prediction method is based on the above system and includes the following steps: Step S100, data acquisition and transmission. The microseismic monitoring sensor array 10 acquires continuous vibration signals from the monitored area in real time, and transmits them to the central processing server 30 via the data transmission unit 20.

[0023] Step S200: Data splitting and preprocessing. The data splitting and preprocessing module 301 receives continuous waveform data and separates it into a background noise data stream and a micro-vibration event data stream based on the signal amplitude characteristics. Simultaneously, it performs mean removal, trend removal, and bandpass filtering on the original signal to eliminate instrument drift and low-frequency environmental interference.

[0024] Step S300: Construct a real-time velocity model. The medium velocity model construction module 302 utilizes the background noise data stream, calculates the cross-correlation function between sensor nodes based on the principle of seismic interferometry, and extracts the Green's function. The volume wave travel time information in the Green's function is processed using a tomographic imaging algorithm to invert and construct a three-dimensional P-wave velocity model and a three-dimensional S-wave velocity model reflecting the current rock mass stress state.

[0025] Step S400: Invert the attenuation parameter model. The attenuation parameter iterative inversion module 303, based on the viscous acoustic wave equation containing the attenuation term, uses the real-time three-dimensional velocity model generated in step S300 as the initial condition to perform reverse time wavefield extension on the microseismic event data stream, calculates the wavefield focusing entropy of the imaging results to quantify the degree of focusing. If the focusing entropy does not meet the convergence criterion, the gradient of the focusing entropy relative to the quality factor is calculated using the adjoint state method, and the quality factor model is iteratively corrected until the optimal quality factor model is obtained.

[0026] Step S500: Reconstruct the energy field and screen waveguide regions. The energy field reconstruction and screening module 304 uses a real-time three-dimensional velocity model and an optimal quality factor model to perform full-waveform reverse time migration on the microseismic data, generating a four-dimensional energy density field containing spatiotemporal evolution information. The waveguide locking index is calculated, and effective energy accumulation regions are screened based on the spatial geometric overlap between the wave velocity gradient field and the energy density field. The average energy density time series of the effective energy accumulation regions is then extracted.

[0027] Step S600: Analyze nonlinear dynamic characteristics. The dynamic characteristic analysis module 305 uses a Gaussian kernel function to smooth the extracted energy density time series to separate the trend term and obtain the residual series. Within a set sliding time window, the lag first-order autocorrelation coefficient and time-domain variance of the residual series are calculated to quantify the critical slowing characteristics of the rock mass system.

[0028] Step S700, Comprehensive Early Warning Decision. The comprehensive early warning decision module 306 monitors the evolution of dynamic indicators. When the first-order autocorrelation coefficient of the lag increases monotonically and approaches the critical value, and the time-domain variance shows a divergent trend, while the wavefield focusing entropy is lower than the effectiveness threshold, the system is determined to have entered a critical slowdown state, triggering a rockburst dynamic risk early warning signal.

[0029] During the background noise data extraction and preprocessing process performed by the data splitting preprocessing module 301, in order to ensure that the subsequent medium velocity model construction module 302 can accurately extract the Green's function based on the high signal-to-noise ratio noise signal, it is necessary to perform rigorous signal screening and conditioning on the original continuous waveform data. This involves not only conventional denoising but also spectral shaping and time-domain normalization processing tailored to the characteristics of seismic interferometry. The specific implementation steps are as follows: In step S201, the data splitting and preprocessing module 301 receives the raw continuous digital signal stream transmitted by the microseismic monitoring sensor array 10. To accurately separate data segments containing effective background noise from the continuous signal, the Short-Long-Time Average Ratio (STA / LTA) algorithm is used for event-triggered detection. For the discrete-time series of any sensor channel, the short-time-window average energy is calculated respectively. and long window average energy : ; ; in, The number of sampling points for the short time window. The number of sampling points for the long time window. This is the index of the sampling point at the current moment. For the summation variable, For a discrete time series at time t The data splitting preprocessing module 301 calculates the ratio of the two values. The system has a preset noise threshold. and event determination threshold When the ratio Continuously below the noise threshold When the ratio is [value missing], the data for the corresponding time period is determined to be background noise data; when [value missing] Above the event determination threshold When the data is identified as a microseismic event, the data splitting and preprocessing module 301 segments the time series according to the identification result and stores it into the background noise buffer queue and the microseismic event buffer queue respectively.

[0030] In step S202, the data splitting preprocessing module 301 performs mean and trend removal processing on the extracted background noise data segments. Because the acquisition instrument may experience zero-point drift or long-period temperature drift, directly calculating the cross-correlation will lead to a DC component deviation in the results. First, the arithmetic mean of the data segments is calculated and subtracted to eliminate the DC component. Then, the least squares method is used to fit the linear trend term in the data segments. And subtract the linear trend term from the original data, where The slope This is the intercept.

[0031] In step S203, to eliminate spectral aliasing and retain low-frequency components effective for deep imaging, the data splitting preprocessing module 301 applies a zero-phase bandpass filter to the detrended background noise data. The passband frequency range of the filter is set to... The passband frequency range is determined based on the superior frequency characteristics of the rock mass in the monitoring area and the frequency response characteristics of the microseismic monitoring sensor array 10. Zero-phase filtering is used to avoid introducing additional time phase delay during the filtering process, ensuring the accuracy of the subsequently extracted Green's function travel time information.

[0032] Step S204: The data splitting preprocessing module 301 performs time-domain normalization on the background noise data. Due to the complexity of the monitoring environment, the background noise often contains non-stationary transient high-energy interference. This high-energy interference dominates the results in cross-correlation calculations, masking the true background noise wavefield characteristics. Therefore, a moving absolute average normalization method is used to suppress high-energy transient interference. For time... data points Its normalized value The calculation is as follows: ; ; in, To half the width of the sliding normalized window, The average absolute amplitude weight within the window. For the summation index within the window, This is a small regularization constant. Through step S204, the phase information of the background noise is preserved, while the influence of abnormal amplitudes is suppressed.

[0033] In step S205, the data splitting preprocessing module 301 further performs spectral whitening processing on the data. Background noise sources are typically unevenly distributed across the spectrum, and the Green's function spectrum obtained through direct cross-correlation will be modulated by the noise source spectrum, leading to reduced resolution. The data splitting preprocessing module 301 transforms the time-domain normalized data to the frequency domain using a Fourier transform to obtain the spectrum. The envelope spectrum is obtained by smoothing the spectral amplitude. And construct the whitened spectrum : ,in, This is the regularization parameter. Subsequently, an inverse Fourier transform is performed... The data is then converted back to the time domain. Spectral whitening expands the effective bandwidth of the background noise data, resulting in a sharper waveform for the subsequently calculated cross-correlation function, which improves the accuracy of travel time acquisition. The data processed through steps S201 to S205 serves as high-quality input data for constructing the medium velocity model.

[0034] After the data splitting and preprocessing module 301 completes the cleaning and standardization of the background noise data, the medium velocity model construction module 302 performs cross-correlation calculations to recover the empirical Green's function between the sensors. This process is based on the principle of seismic interferometry, where the expected cross-correlation of the random wavefields recorded by two receiving points corresponds to the Green's function of the data excited at one point and received at the other. The medium velocity model construction module 302 specifically performs the following steps: Step S301: The medium velocity model construction module 302 pairs and combines the sensor nodes in the microseismic monitoring sensor array 10. Assume the microseismic monitoring sensor array 10 includes... One sensor, construct all possible non-repeating sensor pairs. ,in and For sensor number index, and For each pair of sensors, the medium velocity model building module 302 obtains the preprocessed background noise data stream within the synchronization time period from the data splitting and preprocessing module 301.

[0035] In step S302, the medium velocity model construction module 302 performs segmented cross-correlation calculations on the long-term background noise data stream. To reduce memory consumption and avoid the effects of clock drift in long-term sequences, the sensor... and sensors The data stream is divided into multiple lengths. The time segment. Let the first... Within a time segment, the sensor The data sequence is ,sensor The data sequence is ,in This is the index of the discrete-time sampling points within the segment, with values ​​ranging from 1 to... , For the corresponding The total number of sampling points. The medium velocity model construction module 302 calculates the cross-correlation function of the time segment. : , in the formula, The discrete-time delay variable has a range of values. . For the maximum cross-correlation delay time, The maximum span of the monitoring area and the estimated minimum wave velocity are used to determine the propagation time of the direct wave between the two points covered.

[0036] In step S303, the medium velocity model construction module 302 performs linear superposition processing of the cross-correlation function. Because the signal-to-noise ratio of a single, short-time cross-correlation result is low, it cannot clearly represent the characteristics of the Green's function. Therefore, for all... The cross-correlation functions of each time segment are superimposed to generate the superimposed cross-correlation function. : .

[0037] The superposition operation utilizes the statistical stability of the background noise source under long-term observation, enhances the coherent wave (surface wave or volume wave) signal propagating between two points through constructive interference, and suppresses incoherent random noise through destructive interference.

[0038] Step S304: The medium velocity model construction module 302 recovers the empirical Green's function based on the superimposed cross-correlation function. According to the relationship between cross-correlation and the Green's function in wave theory, the time derivative of the cross-correlation function approximates the Green's function. The medium velocity model construction module 302... By performing time differentiation and taking the negative value, the sensor is obtained. and sensors Empirical Green's function between : .

[0039] In the discrete domain implementation, the medium velocity model construction module 302 uses the central difference method to approximate the above differential operation. This is to eliminate the causal part of the Green's function caused by the asymmetric distribution of noise sources. ) and non-causal part ( To address the asymmetry of the empirical Green's function, a symmetric averaging process is applied to extract its symmetric components. : ; This empirical Green's function symmetric component Includes sensors and sensors The kinematic information of elastic wave propagation between them is obtained, where the wave packet with the largest amplitude corresponds to the direct wave or surface wave signal, and its travel time accurately reflects the wave velocity structure characteristics of the medium. The medium velocity model construction module 302 will calculate the symmetric components of the empirical Green's function for all sensors. The data set serves as the basic input for subsequent tomographic inversion.

[0040] After the medium velocity model construction module 302 obtains the empirical Green's function between each sensor pair, it further performs accurate acquisition of travel time information and tomographic inversion of the three-dimensional velocity structure to construct an initial velocity model that can reflect the heterogeneous characteristics of the rock mass. The specific implementation steps are as follows: Step S305, the medium velocity model construction module 302 constructs the symmetric components of the empirical Green's function obtained in step S304. Envelope analysis is performed to extract the arrival time of the stable signal. Since the Green's function waveform recovered from background noise is complex, directly picking the phase jump point introduces errors; therefore, the Hilbert transform is used to calculate the instantaneous amplitude envelope of the signal. For each sensor pair... Calculate the Hilbert transform of the symmetric components of the empirical Green's function. : Subsequently, the medium velocity model construction module 302 constructs the instantaneous amplitude envelope function. : .

[0041] The medium velocity model construction module 302 searches for the envelope function within a preset physically reasonable time window. The maximum point, and the time delay corresponding to the maximum point. Determined to be a sensor and sensors Observation time between Observation of time travel It represents the actual propagation time of an elastic wave between two points.

[0042] Step S306: The medium velocity model construction module 302 discretizes the rock mass space of the monitoring area into a grid. The monitoring area is defined as a three-dimensional spatial domain. , three-dimensional spatial domain Divided into There are cuboid discrete mesh elements, with a total number of meshes of [number]. Define each grid cell. The slowness (the reciprocal of speed) is ,in For grid indices, the value range is 1 to... The slowness of all grid cells constitutes the parameter vector of the model to be inverted. The medium velocity model construction module 302 sets the initial uniform slow speed model. The value is the reciprocal of the average wave velocity in the monitored area.

[0043] Step S307: The medium velocity model construction module 302 constructs a forward modeling computation based on the Eikonal Equation. To establish a mathematical mapping between the slowness model and travel time, the FastSweeping Method is used to solve the Eikonal Equation. Under the assumption of an isotropic medium, the Eikonal Equation is described as follows: ,in, Spatial grid points The arrival time of the wavefront at that location, Spatial grid points The slowness function. For the th The slow model of the next iteration (in (For iteration number index), medium velocity model building module 302 with sensor Using the location of the source point as the reference point, the above equations are solved numerically to calculate the wavefront propagation to the sensor. Theoretical timekeeping of location .

[0044] Step S308: The medium velocity model construction module 302 establishes a system of linear equations between the travel time residual and the slowness correction. It defines a vector representing the observation travel times of all sensor pairs. The corresponding theoretical time-form vector The residual vectors of the two are According to Fermat's principle, the propagation path of a wave is approximately invariant under small slow-motion perturbations, and the travel time residual and the slow-motion correction are related. They satisfy a linear relationship: In the formula, For a Jacobian matrix, the elements are... Indicates the first The ray path passes through the first The length of each grid cell, where This provides an index for ray paths. Due to the non-uniformity of ray path coverage, linear equation systems are often ill-conditioned. Module 302, which constructs the medium velocity model, introduces Tikhonov regularization constraints to construct the objective functional. : ,in, As the regularization damping factor, This is the Laplace smoothing operator matrix, used to constrain the spatial smoothness of the model and prevent non-physical oscillations in the inversion results.

[0045] In step S309, the medium velocity model construction module 302 uses the least squares QR decomposition (LSQR) algorithm or the joint iterative reconstruction technique (SIRT) to minimize the objective functional and obtain the slowness correction. And update the current slowness model: .

[0046] The medium velocity model construction module 302 determines whether the root mean square error of the travel time residual is less than a preset convergence threshold. If convergence fails, the updated slow-speed model will be used. Repeat the ray tracing and forward modeling calculations in step S307 to reconstruct the Jacobian matrix. Then proceed to the next iteration; if convergence has been achieved, output the final slow-speed model. And take the slow speed model The reciprocal yields the three-dimensional nonhomogeneous velocity model. Three-dimensional heterogeneous velocity model It contains real-time spatial wave velocity distribution information of the rock mass medium in the monitoring area, which is used as the initial velocity field input for the subsequent attenuation parameter iterative inversion module 303 to perform reverse time migration.

[0047] The medium velocity model construction module 302 generates a three-dimensional heterogeneous velocity model containing real-time spatial wave velocity distribution information of the rock mass medium in the monitoring area. Subsequently, the attenuation parameter iterative inversion module 303 further incorporates the attenuation characteristics of the medium to construct a viscous acoustic wave equation to simulate the propagation process of seismic waves in inelastic rock mass. A physical and mathematical model for the propagation process is established for subsequent reverse-time migration imaging and focusing entropy calculation. The specific implementation steps are as follows: Step S401: The attenuation parameter iterative inversion module 303 establishes a viscoelastic constitutive model of the rock mass. To characterize the absorption and dispersion effects of seismic wave energy by the rock mass, the attenuation parameter iterative inversion module 303 uses the Generalized Standard Linear Solid (GSLS) model to approximate the constant... The model, the relationship between stress and strain, is derived through the relaxation modulus function. describe: ,in, Represents the physical time variable. It is the relaxation modulus (i.e., the low-frequency limiting modulus). The number of relaxation mechanisms. For the relaxation mechanism index, For the first The stress relaxation time constant under a relaxation mechanism For the first The strain relaxation time constant under a relaxation mechanism It is a unit step function. The attenuation parameter iterative inversion module 303 determines the attenuation parameter based on the bandwidth requirements of the monitoring frequency band. The value is taken in this embodiment. To balance computational efficiency and simulation accuracy, it corresponds to the standard linear solid model.

[0048] Step S402: The attenuation parameter iterative inversion module 303 derives and constructs the first-order viscous acoustic wave equations in terms of velocity and stress. Under the assumption of an isotropic medium, the attenuation contribution of shear waves is ignored, and the focus is on the attenuation characteristics of longitudinal waves. The wave field state variables are defined to include the particle vibration velocity vector. Sound pressure scalar and memory variables The governing equations consist of equations of motion, state equations, and evolution equations for memory variables: ; ; ; In the formula, It is a spatial position coordinate vector. For the density of the rock mass, For Hamiltonian operators, For the focal force term, This represents the bulk modulus of the rock mass. (Memory variable) It is used to characterize the effect of historical strain on current stress, thereby realizing a mathematical description of energy dissipation in the time domain.

[0049] Step S403: The attenuation parameter iterative inversion module 303 establishes the attenuation parameter model. Quality factors With relaxation time parameter stress relaxation time constant strain relaxation time constant The mapping relationship between them. In order to accurately describe the dissipation characteristics of the medium during the inversion process, a decay parameter model is defined. Quality factor The reciprocal of, that is Set the reference angular frequency (Taking the dominant frequency of the microseismic signal). Based on viscoelastic theory, the attenuation parameter model... The relaxation time satisfies the following constraint equation: ; In the single relaxation mechanism ( Under these conditions, the attenuation parameter iterative inversion module 303 solves the above equations analytically to invert the spatial attenuation parameter model to be retrieved. Converted into specific micro-relaxation parameters: ; ; Through the above mapping, the attenuation parameter iterative inversion module 303 can transform the macroscopic attenuation parameter model. The wave equation in step S402 is directly coupled into the wave field extension process, which enables the wave field extension process to reflect the absorption and attenuation characteristics of the non-uniform spatial distribution of the rock mass medium.

[0050] In step S404, the attenuation parameter iterative inversion module 303 uses the staggered-grid finite difference method to discretize and solve the above viscous acoustic wave equations. The attenuation parameter iterative inversion module 303 approximates the spatial derivative using a higher-order difference operator in the spatial domain and performs time stepping using a second-order difference operator in the time domain. To eliminate the interference of artificial boundary reflections on the wavefield simulation, the attenuation parameter iterative inversion module 303 loads a perfectly matched layer (PML) absorbing boundary condition at the boundary of the computational domain. The discretization solver serves as the forward modeling engine, used in subsequent steps to calculate the predicted wavefield and the backpropagated residual wavefield.

[0051] After constructing the viscous acoustic wave equation and numerical solver, the attenuation parameter iterative inversion module 303 reconstructs the microseismic wavefield using reverse time migration (RTM) technology. The reconstruction process includes the forward extension of the source wavefield and the reverse extension of the receiver wavefield, aiming to recover the true propagation path of the wavefield in the spatiotemporal domain. The specific implementation steps are as follows: In step S405, the attenuation parameter iterative inversion module 303 performs a forward continuation calculation of the source wavefield, based on the microseismic source location determined in step S305. As the excitation point, load the source time function. The three-dimensional heterogeneous velocity model output in step S309 is used. and the decay parameter model under the current iteration step Through the high-order staggered mesh finite difference solver in step S404, along the positive time direction ( to ,in To simulate the total duration, solve the viscous acoustic wave equation.

[0052] During this process, the attenuation parameter iterative inversion module 303 calculates and stores each simulation time. Source pressure wave field The evolution of the source wavefield follows the attenuation characteristics of the medium; that is, as the propagation distance increases, high-frequency energy is gradually attenuated by the medium. The absorption dissipation simulates the real seismic wave propagation process.

[0053] In step S406, the attenuation parameter iterative inversion module 303 performs the reverse continuation calculation of the receiver point wavefield. The attenuation parameter iterative inversion module 303 converts the multi-channel microseismic observation data received by the microseismic monitoring sensor array 10 into... Perform time reversal processing to obtain reversed data. .Will As a boundary source term, it is loaded at the corresponding sensor location, using the same three-dimensional heterogeneous velocity model. and attenuation parameter model Perform wave field simulation.

[0054] To compensate for energy loss during wave propagation in imaging, the attenuation parameter iterative inversion module 303 flips the sign of the dissipation term in the viscous acoustic wave equation during reverse continuation. Specifically, the sign of the memory variable evolution equation in step S402 is adjusted to represent a negative dissipation (i.e., energy gain) process, thereby numerically recovering the amplitude attenuation and dispersion caused by the viscosity of the medium, in reverse time (from... to Solve the equations to obtain the results for each simulation time step. Detector point pressure wave field .

[0055] Step S407: The attenuation parameter iterative inversion module 303 constructs the microseismic source imaging results using cross-correlation imaging conditions. Based on the reverse time migration principle, when the forward-propagating source wavefield and the backward-propagating receiver wavefield coincide at the same time and spatial location, the location is the reflection point or the source point. The attenuation parameter iterative inversion module 303 iterates the source wavefield at each spatial grid point. and detector wave field Perform zero-delay cross-correlation integration to generate spatial imaging values. : .

[0056] To eliminate imaging amplitude distortion caused by uneven source illumination, the attenuation parameter iterative inversion module 303 further normalizes the imaging results using source illumination to obtain normalized imaging results. : ,in, To prevent the use of a stability factor with a denominator of zero, normalized imaging results are used. It intuitively reflects the microseismic energy focusing state under the current velocity and attenuation model, serving as the basic data for subsequent calculation of wavefield focusing entropy.

[0057] In obtaining normalized imaging results Subsequently, the attenuation parameter iterative inversion module 303 requires a quantitative index to evaluate the focusing degree of the normalized imaging results. This invention introduces wavefield focusing entropy as the objective function. The physical meaning is that when the medium model (velocity and attenuation) is accurate, the microseismic energy, after reverse time migration, will ideally focus on the true source location, resulting in a concentrated imaging energy distribution and a minimum entropy value. Conversely, if the medium model is inaccurate, the energy is defocused and its distribution is chaotic, leading to a larger entropy value. Therefore, by minimizing the wavefield focusing entropy, the attenuation parameter model can be iteratively optimized.

[0058] Step S408, the attenuation parameter iterative inversion module 303 processes the normalized imaging result obtained in step S407. The image is then processed to convert it into a normalized probability density function. The attenuation parameter iterative inversion module 303 first performs energy squared normalization on the imaging results to obtain the energy distribution function. : ,in, Represents the entire imaging space domain. This is the integration variable. Through processing, The integral over the entire spatial domain is 1, satisfying the fundamental property of the probability density function. The numerical value reflects the spatial location The relative probability of observing source energy at a location.

[0059] Step S409, the attenuation parameter iterative inversion module 303 is based on the energy distribution function Calculate the current attenuation parameter model Corresponding wave field focusing entropy This invention uses Shannon entropy to define focusing entropy, and the specific calculation formula is as follows: For a discretized three-dimensional mesh space, the above integral form is transformed into a summation form. Let the imaging space be divided into... There are 10 grid cells, and the index of each cell is 1. Then the discretized wavefield focusing entropy The calculation is as follows: ,in, It is the first Energy distribution function values ​​on each grid cell, It is the first The volume of each grid cell.

[0060] In step S410, the attenuation parameter iterative inversion module 303 focuses the wavefield entropy. As the objective function, it is used in the subsequent iterative optimization process. In each iteration, the current wavefield focusing entropy is calculated, and it is determined whether the wavefield focusing entropy has converged. The convergence criterion can be: the current wavefield focusing entropy is less than a preset convergence threshold. Alternatively, the rate of change of wavefield focusing entropy in two consecutive iterations may be less than another preset threshold. If convergence fails, the gradient of the focusing entropy relative to the attenuation parameter model needs to be calculated to guide the model update direction. When the focusing entropy reaches its minimum, it indicates that the energy focusing degree of the reverse time migration imaging result is the highest, and the corresponding attenuation parameter model... It is considered the current optimal solution.

[0061] To achieve an inversion process aimed at minimizing the wavefield focusing entropy, the attenuation parameter iterative inversion module 303 employs the adjoint state method to efficiently calculate the gradient of the objective function relative to the attenuation parameter model. By constructing an adjoint wavefield, the adjoint state method avoids calculating partial derivatives by performing forward perturbations on each model parameter, thereby reducing the cost of gradient calculation to a level comparable to two forward simulations.

[0062] In step S411, the attenuation parameter iterative inversion module 303 constructs the adjoint source term. The adjoint source is a virtual source driving the back propagation of the adjoint wavefield; its physical meaning is the partial derivative of the objective function (wavefield focusing entropy) with respect to the wavefield variables. According to the definition of wavefield focusing entropy in step S409, the attenuation parameter iterative inversion module 303 calculates the wavefield focusing entropy by applying the chain rule. Normalized imaging results The partial derivatives, and further traced back to the pressure wave field of the seismic source. and detector pressure wave field The partial derivatives of . The results of this series of calculations constitute the adjoint source term distributed throughout the entire spatiotemporal domain, denoted as . and The derivation of the adjoint source term can be performed mathematically by those skilled in the art based on the wavefield focusing entropy form and imaging conditions used. This is a conventional technique in the field, and its specific expression will not be elaborated here.

[0063] In step S412, the attenuation parameter iterative inversion module 303 performs backpropagation of the adjoint wave field, introducing two adjoint wave field variables. and , respectively, as the pressure wave field corresponding to the seismic source. and detector pressure wave field The Lagrange multipliers. These two accompanying wavefields follow the same viscous acoustic wave equations as the original wavefield, but are driven by the accompanying source terms calculated in step S411, and in reverse time (from...). to The solution is obtained by performing a numerical inversion of the attenuation parameter iteration inversion module 303, which utilizes the stored source pressure wave field. and detector pressure wave field Combined with the adjoint source term, the adjoint wave field in the entire spatiotemporal domain is calculated using the finite difference method. and .

[0064] Step S413: The attenuation parameter iterative inversion module 303 calculates the gradient of the objective function with respect to the attenuation parameter model. Utilizing the core conclusion of the adjoint state method, the wavefield focusing entropy... For attenuation parameter model gradient It can be obtained through the cross-correlation of the forward wave field and the adjoint wave field. The formula for calculating the gradient in the attenuation parameter iterative inversion module 303 is as follows: ,in, Represents the operator of the viscous acoustic wave equation. Represents the operator pair decay parameter model The partial derivatives. Specifically, the partial derivatives act on the wave field variables, and the non-zero terms mainly originate from the attenuation-related parts of the state equation of the viscous acoustic wave equation and the memory variable evolution equation in step S402. Through the gradient formula, the attenuation parameter iterative inversion module 303 can obtain a gradient field of the same size as the model parameter space with low computational cost. .

[0065] In step S414, the attenuation parameter iterative inversion module 303 updates the attenuation parameter model using a nonlinear optimization algorithm. After obtaining the gradient, the attenuation parameter iterative inversion module 303 uses an iterative optimization strategy, such as the conjugate gradient method (CG) or the L-BFGS algorithm, to update the attenuation parameter model. Let the current step be the... In this iteration, the model update process mainly includes: First, the decay parameter iterative inversion module 303 calculates an effective descent search direction based on the current gradient and historical iteration information. To improve the uneven illumination problem of the gradient, preconditioning operators, such as a diagonal matrix based on ray illumination, can be applied to correct the gradient and accelerate convergence.

[0066] Then, the attenuation parameter iterative inversion module 303 determines an optimal step size along the calculated search direction using a linear search method (such as parabolic fitting or backtracking search). This process aims to find a step size value that minimizes the objective function (wavefield focusing entropy) in the search direction.

[0067] Finally, the decay parameter iterative inversion module 303 moves the current decay parameter model along the search direction by the optimal step size to obtain the updated model. To ensure the physical validity of the model, the decay parameter iterative inversion module 303 imposes constraints on the updated model, such as through projection or truncation, to ensure that the decay parameter value is always positive.

[0068] Step S415, the attenuation parameter iterative inversion module 303 updates the attenuation parameter model. As new input, the process returns to step S405 to begin a new round of reverse-time migration wavefield extension, focusing entropy calculation, and gradient solving until the convergence condition set in step S410 is met. When the iteration terminates, the output is the optimal quality factor model (obtained through...). (Originated from conversion).

[0069] To obtain the optimal three-dimensional heterogeneous velocity model and quality factor Subsequently, the waveguide effect analysis module 304 initiates the reconstruction process of the spatiotemporal evolution characteristics of the wavefield energy of the microseismic event. The core of this process lies in generating a data volume that can characterize the flow and accumulation of energy in the underground medium, namely a four-dimensional energy density field.

[0070] In step S501, the waveguide effect analysis module 304 performs the final positive wave field simulation, using the microseismic source location determined by the microseismic data preprocessing and location module 302. and source time function As an excitation condition, the key is that the medium model used in this forward wave field simulation is the optimal three-dimensional heterogeneous velocity model obtained through iterative inversion. and optimal decay parameter model The energy spatiotemporal evolution analysis module 304 calls the viscous acoustic wave numerical solver from step S404, along the positive time direction ( to Solve the viscous acoustic wave equation to calculate the evolution of the source wave field over time in the entire three-dimensional spatial domain under the optimal medium model.

[0071] In step S502, the waveguide effect analysis module 304 synchronously extracts the wavefield dynamic parameters. This occurs at each simulation moment in the forward simulation. The waveguide effect analysis module 304 extracts physical quantities at each spatial grid point within the computational region in real time, including pressure wave fields in scalar form. and particle vibration velocity in vector form Among them, the particle vibration velocity It contains components in three orthogonal directions. This process generates physical field data that evolves continuously over time, providing the fundamental variables for subsequent energy density calculations.

[0072] Step S503, waveguide effect analysis module 304 calculates the instantaneous total energy density in the entire space. Seismic waves, as a type of elastic wave (or pressure wave in the acoustic approximation), have a total energy consisting of kinetic and potential energy as they propagate through a medium. For each spatial grid point... and every simulation moment The module performs the following calculations: Regarding the kinetic energy component, i.e., instantaneous kinetic energy density The value depends on the inertial properties of the medium and the motion state of the particles. The waveguide effect analysis module 304 analyzes the waveguide effect based on a preset medium density model. and the extracted particle vibration velocity The calculations are performed, and the numerical value is equal to the medium density model. and the extracted particle vibration velocity One half of the product of the squares of the modulus.

[0073] Regarding the potential energy component, i.e., the instantaneous potential energy density The value depends on the elastic properties and stress state of the medium. The waveguide effect analysis module 304 is based on the optimal three-dimensional heterogeneous velocity model. Medium density model and the extracted pressure wave field Calculations show that the value equals the pressure wave field. The square of the amplitude divided by twice the bulk modulus of the medium, where the bulk modulus is obtained through the medium density model. With the optimal three-dimensional heterogeneous velocity model It is obtained by multiplying the squares of the products.

[0074] Based on the two components mentioned above, the waveguide effect analysis module 304 superimposes the instantaneous kinetic energy density and the instantaneous potential energy density at the same point to obtain the instantaneous total energy density. The calculation formula is as follows: By performing the above calculations on all simulation times and all spatial grid points within the simulation period, the waveguide effect analysis module 304 constructs a complete dataset containing three spatial dimensions and one temporal dimension, i.e., a four-dimensional energy density field. The dataset stores... Value quantification describes the entire process of seismic wave energy in a heterogeneous attenuating medium, from excitation and diffusion to capture by geological structures.

[0075] Having obtained a four-dimensional energy density field describing energy flow (i.e., the instantaneous total energy density at all locations and times) After generating the dataset, the waveguide effect analysis module 304 needs to automatically and objectively identify waveguide structures with energy trapping characteristics from this massive dataset. Traditional visual interpretation methods are highly subjective and inefficient. To address this, this invention proposes a novel quantitative index, the waveguide locking index, to measure the ability of any point in space to capture and maintain seismic wave energy.

[0076] Step S601: The waveguide effect analysis module 304 performs time integration based on the four-dimensional energy density field to obtain the time-integrated energy of each spatial grid point. The core characteristic of waveguide structures is their ability to confine energy within a specific region for an extended period, resulting in accumulated energy in these regions far exceeding that of the surrounding medium throughout the entire earthquake event. To quantify this characteristic, the waveguide effect analysis module 304 analyzes the energy of each spatial grid point within the study area. Along the time axis, the instantaneous total energy density Integrating is performed to obtain the spatial grid points. Time integral energy : ,in, Spatial grid points Time integral energy at that point This is the total duration of the forward simulation, which needs to be long enough to cover the complete process from the source excitation to the energy decay in the main region. After the above steps, the four-dimensional energy density field is reduced to a three-dimensional time-integrated energy field. This visually highlights the spatial regions through which earthquake energy mainly flows and resides during its propagation.

[0077] In step S602, the waveguide effect analysis module 304 acquires the peak energy at each spatial grid point. Besides the accumulated energy, the intensity of energy capture is also a crucial manifestation of the waveguide effect. The waveguide effect analysis module 304 traverses the stored four-dimensional energy density field to identify each spatial grid point. Total instantaneous energy density experienced throughout the entire simulation period The maximum value is denoted as peak energy. .

[0078] Step S603, waveguide effect analysis module 304 constructs waveguide locking index ( The waveguide locking index is designed as a comprehensive indicator that simultaneously reflects both energy accumulation effects (duration) and energy trapping intensity (peak value). The waveguide effect analysis module 304 analyzes each spatial point within the study area. By integrating energy over time With peak energy Multiply to calculate spatial grid points waveguide locking index : ,in, Spatial grid points Waveguide locking index at the location.

[0079] Based on the waveguide locking index field, the waveguide effect analysis module 304 can effectively distinguish between real waveguide structures with energy trapping characteristics and non-critical regions that only serve as transient paths. A high waveguide locking index value characterizes the spatial grid points. It possesses the dual characteristics of high energy peak and long duration.

[0080] In step S604, the waveguide effect analysis module 304 filters waveguide structures according to a preset threshold. The waveguide effect analysis module 304 sets the filtering threshold. , will satisfy The spatial grid points were identified as waveguide structure points. The spatial collection of these waveguide structure points constitutes the waveguide channel revealed by this microseismic event.

[0081] Step S701, the waveguide effect analysis module 304 constructs a binary spatial mask. The waveguide effect analysis module 304 uses a screening threshold... Generate a key region mask with the same size as the 3D mesh. For each spatial grid point within the study area. ,like Then assign a value Otherwise, assign a value. Meanwhile, the waveguide effect analysis module 304 generates a background region mask through a logical NOT operation. .

[0082] Step S702: The waveguide effect analysis module 304 counts the number of grid cells in the region. To perform subsequent energy normalization, the waveguide effect analysis module 304 performs the following steps: and Summing across the entire space yields the total number of grid points in the key region. and the total number of grid points in the background area .

[0083] Step S703: Waveguide effect analysis module 304 extracts the average energy time series of the key region. Waveguide effect analysis module 304 traverses the four-dimensional energy density field, for each simulation time step... Using a key region mask to measure the instantaneous total energy density of the entire space Spatial filtering and averaging were performed to obtain the energy evolution curves of the key regions. The calculation formula is expressed as follows: ; in, Indicates the critical region at simulation time. Average instantaneous energy density, This represents the total number of grid points within the critical area. Represents the entire three-dimensional simulation space domain. For critical area masking, when spatial grid points Waveguide locking index at the location The value is 1 if it is true, and 0 otherwise.

[0084] In step S704, the waveguide effect analysis module 304 extracts the average energy time series of the background region as a reference. Using the same processing logic as in step S703, the waveguide effect analysis module 304 utilizes a background region mask. and the total number of grid points in the background area The background region is calculated at each simulation time. The average instantaneous energy density was used to construct the energy evolution curve of the background region. The calculation formula is as follows: ,in, Indicates the background area at simulation time Average instantaneous energy density, This indicates the total number of grid points within the background area. As a mask for the background area, when the spatial grid points The value is 1 if it belongs to the background area, and 0 otherwise.

[0085] By comparing the energy evolution curves in key areas Energy evolution curves of the background region ,like The decay rate is significantly lower than Furthermore, maintaining a high energy level confirms that the selected key regions possess energy trapping characteristics.

[0086] To identify whether a rock mass system is in a critical instability state, it is necessary to separate the micro-fluctuation characteristics under the macro-energy decay trend based on the critical slowing theory. The precursor identification module 306 performs detrending processing on the non-stationary energy decay sequence to extract the stationary residual sequence.

[0087] Step S801: The precursor identification module 306 extracts the analysis time window. The precursor identification module 306 analyzes the energy evolution curves in key regions. The tailwave attenuation segment, located after the peak time and meeting the signal-to-noise ratio requirements, is selected as the analysis object. The selected time interval is set as follows. ,in This is the starting point after the energy peak. The cutoff time is specified. The precursor identification module 306 extracts data within the time interval to construct and analyze the state sequence. This reflects the transient response of the rock mass medium affected by the waveguide effect during the energy dissipation process.

[0088] In step S802, the precursor identification module 306 calculates the deterministic trend term. To separate long-term energy dissipation behavior from short-term random fluctuations, the precursor identification module 306 uses a Gaussian kernel function to analyze the state sequence. Perform convolutional smoothing. Convolutional smoothing filters out high-frequency noise and preserves the overall shape and trend of the sequence. The calculation is achieved through discrete convolution, as shown in the following formula: , in, This indicates the current point in time when the trend is being calculated. This represents the summation traversal time points within the valid time window. Indicates the effective time window The total number of time steps contained within. This is a Gaussian kernel function, used to define the weights of the contribution to the trend at different time points. Its expression is: ,in, Indicates time difference, The smoothing bandwidth is a key parameter that controls the smoothness of the trend fit.

[0089] In this embodiment, smooth bandwidth The value needs to be set according to the characteristic time scale of the sequence. The precursor identification module 306 selects... The value should be large enough to ensure that only high-frequency random fluctuations are filtered out without overfitting the short-term oscillations of the system, thus guaranteeing the extracted trend term. It accurately reflects the physical baseline of energy dissipation in the rock mass.

[0090] In step S803, the precursor identification module 306 performs a difference operation to extract the random residual sequence. This yields a deterministic trend term. Then, the precursor identification module 306 will analyze the state sequence. A point-by-point subtraction operation is performed with the trend term to remove the background trend and retain the fluctuation component reflecting the microdynamic behavior of the system. The resulting residual sequence... The calculation is as follows: ,in, Indicates at time The energy residual sequence.

[0091] Residual sequence The macroscopic energy decline trend caused by the geometric diffusion of seismic waves and the inherent attenuation of the medium is removed, transforming the data into a quasi-stationary sequence fluctuating around zero. This step forms the mathematical basis for subsequent calculations of critical slowing indicators such as variance and lag first-order autocorrelation coefficient, ensuring that the calculation of statistical characteristics reflects only changes in the system's internal stability, rather than interference from external macroscopic trends. By extracting the residual sequence, this invention achieves a dimensional shift from macroscopic energy propagation monitoring to microscopic system stability analysis.

[0092] After obtaining the residual series that has eliminated the interference of macroeconomic trends. Subsequently, the precursor identification module 306 further quantifies the memory effect of the rock mass system at the critical point. This invention uses the sliding window technique to calculate the hysteresis first-order autocorrelation coefficient of the residual sequence, thereby capturing the dynamic characteristics of the system's reduced recovery rate to disturbances.

[0093] Step S901: The precursor identification module 306 constructs a time-sliding window. Given the non-stationarity of earthquake coma characteristics, the precursor identification module 306 uses the residual sequence... Set the window length to A sliding window. The module sets the sliding step size. The window is gradually moved backward from the start of the sequence until it covers the entire effective analysis period, and a series of local residual subsequences are extracted.

[0094] Step S902, the precursor identification module 306 calculates local statistical features. For the center time is The sliding window module extracts local residual subsequences within the window. To correct for local mean drift and normalize fluctuation amplitude, the precursor identification module 306 first calculates the local mean of the data within the sliding window. and time domain variance Local mean Calculated as Time-domain variance Calculated as .

[0095] In step S903, the precursor identification module 306 calculates the first-order autocorrelation coefficient. After completing the centering process of the local data, the precursor identification module 306 calculates the Pearson correlation coefficient between the residual sequence within the window and its own sequence lagged by one time step, which is the first-order autocorrelation coefficient. The calculation formula is as follows: ; in, Indicates time The lagged first-order autocorrelation coefficient, with a range of values. , For the first in the window Each residual data point corresponds to a simulation time. , For the first in the window Each residual data point corresponds to a simulation time. , This is the window length.

[0096] When the system is far from the critical point, the residuals exhibit memoryless random noise. Approaching zero; as the system nears instability, the critical slowing effect causes a high positive correlation between the residuals at adjacent time points. It rises and approaches 1.

[0097] In step S904, the precursor identification module 306 generates the evolution curve of the first-order autocorrelation coefficient. The precursor identification module 306 then uses the curves calculated by each sliding window... By connecting them in chronological order, the evolution curve of the first-order autocorrelation coefficient with lag is generated. The upward trend and amplitude level of the evolution curve of the first-order autocorrelation coefficient with lag directly quantify the evolution process of the internal dynamic stability of the rock mass, serving as a key indicator for determining waveguide risk.

[0098] Besides the first-order autocorrelation coefficient, another precursory feature revealed by the critical slowdown theory is the abnormally increased amplitude of fluctuations in the system's state variables. When a rock waveguide system approaches its instability critical point, the curvature of the system's potential energy function tends to flatten, and the restoring force weakens, leading to background noise of equal intensity generating state fluctuations with larger amplitudes. This invention quantifies the changes in the intensity of micro-fluctuations within the rock mass medium by calculating the evolution of the time-domain variance of the residual sequence, using this as an auxiliary indicator for identifying waveguide lock-in risk.

[0099] In step S905, the precursor identification module 306 reuses a time sliding window to segment the residual sequence. To ensure the correspondence of statistical indicators on the time axis, the sliding window mechanism in step S901 is used. The precursor identification module 306 uses the same window length. and sliding step size For residual sequences Perform a traversal. For each central time point... The module uses a sliding window to obtain local residual subsequences within the window. This ensures the consistency of the variance index and the lagged first-order autocorrelation coefficient index in terms of spatiotemporal resolution.

[0100] In step S906, the precursor identification module 306 calculates the temporal variance value within each window. For each sliding window, the precursor identification module 306 uses the local mean... The data is centered, and the average of the sum of squares of the deviations of the residual data from the mean within the window is calculated. This average is the time-domain variance. The calculation formula is as follows: ; in, Indicates at time The time-domain variance, For the first in the sliding window Each residual data point This represents the local mean of the residual data within the current sliding window. This represents the total number of data points within the sliding window.

[0101] In step S907, the precursor identification module 306 generates a variance evolution curve and performs trend analysis. The precursor identification module 306 calculates all the variance curves obtained from the sliding window... Arranged according to time series, a time-domain variance evolution curve is constructed. The time-domain variance evolution curve reflects the increasing instability of the rock mass system during energy dissipation. If an abnormal upward trend is detected in the time-domain variance evolution curve during the tail phase, i.e. This indicates a decrease in system recovery stiffness, an increase in stress concentration within the waveguide structure, and a tendency for the system to approach critical instability. The precursor identification module 306 synchronously outputs the time-domain variance evolution curve and the lag first-order autocorrelation coefficient evolution curve, providing multi-dimensional support for early warning decision-making.

[0102] Refer to Appendix 2, Figure 2 This is a schematic diagram of a time sliding window division and data interception mechanism according to an embodiment of the present invention.

[0103] To capture the dynamic characteristics of the rock mass system evolving from a stable state to a critical unstable state, the precursor identification module 306 constructs a parameter-adjustable sliding window mechanism to detect residual sequences. The refined temporal domain segmentation ensures the spatiotemporal resolution of subsequent statistical indicator calculations.

[0104] In step S908, the precursor detection module 306 determines the window length and overlap step size. Window length The selection needs to balance statistical stability and temporal resolution. The precursor identification module 306 sets the window length based on the dominant frequency characteristics of the residual sequence. This ensures that the time span it covers includes at least one complete main period of the residual fluctuation. Simultaneously, a sliding step size is set. The step size setting controls the speed at which the window moves along the time axis, and determines the smoothness and time resolution of the indicator curve.

[0105] Step S909, the precursor identification module 306 calculates the total number of windows. This is based on the total number of data points in the residual sequence. Given the set window parameters, the precursor identification module 306 calculates the total number of sliding windows generated within the effective analysis period. The boundary conditions for subsequent iterative processing are established, as shown in the following formula: ; in, This indicates the total number of sliding windows. Representing the residual sequence Total number of sampling points This indicates the length of the window within a single sliding window. Indicates the sliding step size. This indicates the floor function.

[0106] In step S910, the precursor identification module 306 performs subsequence extraction based on index mapping. The precursor identification module 306 establishes a subsequence extraction method based on window sequence numbers. The mapping relationship to the original sequence time index. For the first... A sliding window ( ), which in the original residual sequence Starting index in and terminating index Defined as: ; ; Based on the above index range, the module extracts the corresponding local residual subsequences. : ; Local residual subsequence This refers to the step described above used to calculate the local mean. Time domain variance and the first-order lag autocorrelation coefficient The input dataset.

[0107] In step S911, the precursor identification module 306 performs high overlap rate data multiplexing. To improve the ability to detect sudden precursors, this invention sets... This causes data overlap between adjacent windows. Window overlap rate. The definition is as follows: ; Precursor identification module 306 will It remains within a preset range. This high overlap mechanism ensures that the residual data at each time step is sampled multiple times and contributed to the statistical calculations of different time windows, thereby ensuring the continuity of the generated evolution curve and avoiding information loss due to window boundary segmentation.

[0108] Step S1001, the precursor identification module 306 normalizes the time-domain variance index. This is because the first-order autocorrelation coefficient is lagging behind. The range of values ​​is while time domain variance The magnitudes of the two variables differ depending on the absolute intensity of the input energy field. To achieve linear superposition, the precursor identification module 306 uses the variance level of the background region or stationary phase as a benchmark to perform dimensionless processing on the time-domain variance and calculate the normalized variance exponent. The calculation formula is as follows: ; in, Indicates time The normalized variance index, with a range of values. , The time-domain variance of the current window. and These represent the statistical minimum and maximum values ​​of the time-domain variance within the effective analysis period, respectively.

[0109] In step S1002, the precursor identification module 306 constructs a comprehensive instability probability index. Based on the critical slowdown characteristic, the precursor to system instability is characterized by a first-order autocorrelation coefficient approaching 1 and an increasing variance. The precursor identification module 306 introduces a weighted coupling algorithm to fuse two independent indicators into a comprehensive instability probability index. The formula is as follows: ; in, For a moment The comprehensive instability probability index, and For the weighting coefficients, satisfying This is used to adjust the weight of memory effect and wave intensity in the identification process (for hard rock waveguides, it is usually taken as...). ), This is a step function; it takes the value 1 when the input value is greater than 0, and 0 otherwise. (Comprehensive instability probability index) Ensure that risk is only factored in when the system exhibits a positive correlation, excluding negative correlation fluctuations.

[0110] Step S1003: The precursor detection module 306 sets dual threshold criteria and persistence checks. To improve the robustness of the warning, this invention employs dual constraints of amplitude and trend. The precursor detection module 306 presets a yellow warning threshold. and red alert threshold .

[0111] For calculation time To trigger an alert, the following conditions must be met simultaneously: Condition 1: Amplitude exceeds the limit.

[0112] , in, Corresponding to yellow or red thresholds.

[0113] Condition 2: The trend is upward.

[0114] The precursor identification module 306 calculates the index in recent times. Average rate of change within a window ,Require: .

[0115] In step S1004, the precursor identification module 306 outputs a graded early warning signal.

[0116] when Furthermore, when the trend of growth is satisfied, the waveguide structure is determined to have entered a subcritical state, and a yellow warning is output to indicate that stress is beginning to concentrate. when Furthermore, when the trend of increasing intensity is met, the waveguide structure is determined to have entered a critical instability state, and a red alert is issued, indicating the formation of a waveguide locking effect, obstructed energy dissipation, and a risk of rockburst, requiring depressurization measures. By combining criteria, wavefield energy data is converted into risk level signals, enabling early identification of waveguide locking disasters.

[0117] In this embodiment, a deep hard rock tunnel excavation face with obvious layered joint structure is selected as the monitoring object. The energy accumulation process induced by far-field disturbance is simulated to verify the effectiveness of the multi-dimensional coupling early warning mechanism.

[0118] Step S1101: The precursor identification module 306 initializes the monitoring parameters. Based on the background noise level and rock acoustic characteristics of the monitoring area, the module sets the window length. (Sampling points), covering approximately twice the main fluctuation period; set the sliding step size. Sampling points are used to achieve window overlap. In the joint early warning criterion, the weight of the first-order autocorrelation coefficient with lag is set. Normalized variance weights Set a yellow alert threshold. (Subcritical state), red alert threshold =0.85 (critical instability state).

[0119] In step S1102, the precursor identification module 306 performs dynamic feature extraction of the real-time data stream. The precursor identification module 306 extracts the residual sequence during the tailwave decay phase after the elastic disturbance occurs. In the early stages of monitoring ( The rock mass is in the elastic recovery stage, with a first-order autocorrelation coefficient lagging behind. Oscillates around 0, normalized variance index At this point, the overall instability probability index If the value remains below 0.3 and the threshold is not triggered, the system determines that the region is in a safe dissipation state.

[0120] Step S1103, the precursor identification module 306 captures critical slowing precursors and triggers a graded early warning. At that moment, the waveguide locking effect became apparent. Monitoring showed... Monotonous rise, to When the value exceeds 0.6, it indicates that the system's recovery rate to micro-perturbations decreases; at the same time It rises above 0.5. At this point... Exceed And the rate of change The precursor detection module 306 generates a yellow warning signal, marking the area as a stress accumulation zone.

[0121] In step S1104, the precursor detection module 306 determines that waveguide locking has occurred and issues an emergency alarm. to Within the range, standing wave resonance forms inside the waveguide structure. Approaching the saturation value of 0.92, It reached a peak of 0.95. (Comprehensive Instability Probability Index) Breakthrough 1.0. The precursor detection module 306 determines that the region has entered a waveguide-locked critical state, triggering a red alarm and sending a control command. This method provides early warning when the energy has not reached the destructive limit but the system stability has undergone a qualitative change. The specific threshold can be adaptively adjusted according to geological conditions and historical data, which falls within the technical scope of this invention.

Claims

1. A method for predicting the dynamic risk of rockburst by integrating real-time microseismic data, applied in a system comprising a microseismic monitoring sensor array (10), a data transmission unit (20), and a central processing server (30), characterized in that, Includes the following steps: S1, the data splitting and preprocessing module (301) separates the acquired continuous vibration signal into a background noise data stream and a micro-vibration event data stream; S2, the medium velocity model construction module (302) uses the background noise data stream to construct a real-time three-dimensional heterogeneous velocity model that reflects the stress state of the rock mass; S3, The attenuation parameter iterative inversion module (303) uses the real-time three-dimensional heterogeneous velocity model as the initial condition and uses the microseismic event data stream to invert the optimal quality factor model. S4, Energy Field Reconstruction and Screening Module (304) generates a four-dimensional energy density field based on the real-time three-dimensional heterogeneous velocity model and the optimal quality factor model, calculates the waveguide locking index to screen key regions and extracts the energy evolution curve of key regions; S5. The dynamic characteristic analysis module (305) performs detrending processing on the energy evolution curve of the key region to obtain the residual sequence, and calculates the lag first-order autocorrelation coefficient and time-domain variance of the residual sequence. S6. The comprehensive early warning decision module (306) determines the critical slowdown characteristics and triggers an early warning based on the evolution trend of the first-order autocorrelation coefficient and the time-domain variance.

2. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 1, characterized in that, In S1, the data splitting preprocessing module (301) calculates the ratio of the short-time to the long-time average of the continuous vibration signal; The background noise data stream and the microseismic event data stream are extracted based on the comparison result between the short-time average ratio and the preset threshold. The background noise data stream is subjected to detrending processing, zero-phase bandpass filtering, time-domain normalization processing, and spectral whitening processing to serve as the input data for S2.

3. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 1, characterized in that, The S2 step specifically involves, The medium velocity model construction module (302) calculates and superimposes the cross-correlation function of the background noise data streams between sensor pairs based on the principle of seismic interferometry to recover the empirical Green's function; Envelope analysis is performed on the empirical Green's function to pick up the observed travel time, and the theoretical travel time is calculated using the equation of process function; A system of linear equations reflecting the difference between the observed travel time and the theoretical travel time is constructed, and regularization constraints are introduced. The real-time three-dimensional heterogeneous velocity model is obtained by iteratively solving the equations using a tomographic imaging algorithm.

4. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 1, characterized in that, The S3 step is specifically as follows: The attenuation parameter iterative inversion module (303) constructs the viscous acoustic wave equation based on the real-time three-dimensional heterogeneous velocity model; An initial attenuation parameter model is set, and based on the viscosonic wave equation, the attenuation parameter model is used to perform reverse time wavefield extension on the microseismic event data stream to generate normalized imaging results and calculate wavefield focusing entropy. The gradient of the wavefield focusing entropy relative to the attenuation parameter model is calculated using the adjoint state method. The attenuation parameter model is then iteratively updated using the gradient until convergence. The converged attenuation parameter model is then converted into the optimal quality factor model.

5. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 1, characterized in that, In step S4, the step of generating the four-dimensional energy density field includes: The energy field reconstruction and screening module (304) performs full-waveform positive wave field simulation on the microseismic event data stream under the optimal quality factor model and the real-time three-dimensional heterogeneous velocity model; For each spatial grid point at each simulation moment, the instantaneous total energy density is calculated by superimposing the kinetic energy density and potential energy density from the full-waveform forward wave field simulation results. Based on the distribution of the instantaneous total energy density in the spatiotemporal dimension, the four-dimensional energy density field is constructed.

6. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 5, characterized in that, In step S4, the steps of calculating the waveguide locking index to screen key regions and extracting the energy evolution curves of the key regions include: The energy field reconstruction and filtering module (304) integrates the instantaneous total energy density along the time axis for each of the spatial grid points in the four-dimensional energy density field to obtain the time integral energy, and extracts the peak energy of the spatial grid points at all the simulation times; Multiplying the time-integrated energy by the peak energy yields the waveguide locking index of the spatial grid point; The key region is formed by selecting spatial grid points whose waveguide locking index is greater than or equal to a preset screening threshold, and the four-dimensional energy density field is spatially filtered using the spatial location information of the key region to extract the energy evolution curve of the key region.

7. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 1, characterized in that, In step S5, obtaining the residual sequence includes: The dynamic characteristic analysis module (305) extracts analysis data from the tail wave attenuation segment of the energy evolution curve in the key region; The analysis data is smoothed by convolution using a Gaussian kernel function, and a deterministic trend term is calculated. The analytical data and the deterministic trend term are subjected to point-by-point differencing to remove the deterministic trend term, thereby obtaining the residual sequence.

8. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 1, characterized in that, In step S5, calculating the first-order lag autocorrelation coefficient and the time-domain variance includes: The dynamic feature analysis module (305) moves a sliding window on the residual sequence to extract local residual subsequences; For each sliding window, the local mean of the local residual subsequence is calculated, and the average of the sum of squares of the deviations of the local residual subsequence from the local mean is calculated to obtain the time-domain variance corresponding to the sliding window; The variance of the time domain is used to perform variance normalization on the local residual subsequence within the sliding window, and the Pearson correlation coefficient between the local residual subsequence and the sequence lagged by one time step is calculated to obtain the first-order autocorrelation coefficient corresponding to the sliding window.

9. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 1, characterized in that, In step S6, determining the critical slowing characteristic includes: The integrated early warning decision module (306) performs amplitude normalization processing on the time domain variance to obtain the normalized variance index; A weighted coupling algorithm is introduced to perform a weighted summation of the first-order autocorrelation coefficient and the normalized variance index to generate a comprehensive instability probability index. In the weighted summation process, a positive value constraint is set, and the first-order autocorrelation coefficient is included in the calculation of the comprehensive instability probability index only when the first-order autocorrelation coefficient is positive. Using the comprehensive instability probability index as a quantitative indicator to characterize system stability, the growth trend of the comprehensive instability probability index over time is determined as the critical slowdown characteristic.

10. The intelligent rockburst dynamic risk prediction method integrating real-time microseismic data according to claim 9, characterized in that, In S6, the logic for triggering the warning is as follows: The integrated early warning decision module (306) monitors the integrated instability probability index in real time; When the comprehensive instability probability index exceeds the preset yellow warning threshold and shows an increasing trend, a yellow warning signal is output. When the comprehensive instability probability index exceeds the preset red warning threshold and shows an increasing trend, a red alarm signal is output.

Images