A deep rock mass disaster prediction method and system

Abstract

This invention relates to the field of rock mass disaster monitoring and early warning technology, and particularly to a method and system for predicting deep rock mass disasters. The method includes: acquiring deep rock mass microseismic monitoring signals and performing variational mode decomposition; obtaining a denoised microseismic target signal through weighted reconstruction based on the product of the inverse of kurtosis and the energy contribution rate; generating a parameter set through grid search within the delay time and embedded dimension space, and reconstructing the phase space of the denoised microseismic target signal; calculating the correlation integral to generate a double logarithmic coordinate curve; identifying linear scaling regions and determining the fitting confidence weights based on the second derivative of the data points; obtaining a comprehensive fractal scaling index through weighted averaging of the microseismic fractal scaling index using linear fitting; and determining that a disaster is imminent when the rate and magnitude of decrease of the comprehensive fractal scaling index both exceed a set threshold within a continuous time window. This invention can reduce subjective errors in parameter selection and improve the accuracy of disaster early warning.

Description

Technical Field

[0001] This invention relates to the field of rock mass disaster monitoring and early warning technology, and more specifically, to a method and system for predicting deep rock mass disasters. Background Technology

[0002] Deep rock masses exist in complex geological environments characterized by high ground stress, high ground temperature, and high karst water pressure. Disturbance from engineering excavation can easily trigger sudden dynamic disasters such as rock bursts and rockbursts. The operation of rock drills, drilling and blasting, and the electromagnetic environment at deep construction sites generate strong background noise. This mechanical and environmental noise masks the weak precursor signals released by micro-fractures in the rock mass, resulting in a low signal-to-noise ratio for monitoring signals. Under these conditions, conventional time-domain statistical or frequency-domain analysis methods are insufficient to accurately reflect the internal damage evolution of the rock mass, leading to difficulties in disaster prediction due to lag and a high rate of misjudgment.

[0003] In the prior art, Chinese patent document CN106814396B discloses a noise reduction and filtering method for mine microseismic signals based on variational mode decomposition. This method performs variational mode decomposition on noisy microseismic signals and calculates the cross-correlation coefficients between the original signal and each variational mode component. Components with a center frequency higher than a set value and a cross-correlation coefficient lower than a set threshold are filtered out as noise. The remaining components are then reconstructed by summing them with equal weights. While this method uses cross-correlation coefficients and center frequency as the selection criteria, it employs equal weighting for the retained components during the reconstruction stage, failing to establish a differentiated weighting mechanism based on the inherent statistical characteristics of each component, such as kurtosis and energy contribution rate. Under strong noise conditions, this may result in the retention of some noise-dominant components or the suppression of low-energy components containing precursor information, leading to insufficient signal-to-noise ratio after noise reduction.

[0004] Phase space reconstruction technology and fractal scaling indices can be used to identify the nonlinear evolution of signals in multidimensional space. When microcracks penetrate the rock mass and the system tends towards macroscopic instability, the fractal scaling index of the signal usually shows a decreasing trend. However, the above methods have the following limitations: the delay time and embedding dimension required for phase space reconstruction usually depend on empirical selection or determination by a single criterion. Different parameter selections have a significant impact on the calculation results of the fractal scaling index, and a single parameter combination cannot cover the multi-scale nonlinear characteristics of complex systems. When solving for the fractal scaling index, the delineation of the linear scaling region relies on manual visual judgment and lacks quantitative identification criteria based on the local curvature of the curve. In addition, there is a lack of confidence assessment and weighted fusion mechanism based on the fitting quality of the scaling region among multiple fractal scaling indices obtained under different parameter conditions, resulting in large fluctuations in the calculation results and limiting the accuracy of disaster early warning. Summary of the Invention

[0005] To address the technical problems of existing deep rock mass monitoring signal noise reduction lacking modal weighting and fractal scaling index calculation relying on empirical parameters, this invention provides solutions in the following aspects.

[0006] In a first aspect, the present invention provides a method for predicting deep rock mass disasters, comprising: S1, acquiring deep rock mass microseismic monitoring signals and performing variational mode decomposition to obtain intrinsic mode functions reflecting the characteristics of different frequency bands of microseismic waves; performing weighted reconstruction based on the product of the reciprocal of the kurtosis value representing the microseismic impact characteristics of each intrinsic mode function and the energy contribution rate representing the microseismic radiation energy to obtain a denoised microseismic target signal; performing grid search within a preset delay time and embedded dimension space to generate a parameter set; and using the parameter set to perform phase space reconstruction on the denoised microseismic target signal to obtain a reconstructed phase space system representing the rock mass damage evolution; S2, calculating the correlation integral of any reconstructed phase space in the reconstructed phase space system at different scale radii to generate... The double logarithmic coordinate curve is used to calculate the second derivative of the data points. Continuous intervals where the absolute value of the second derivative of the data points is lower than a preset threshold are identified as linear scaling regions reflecting the scale-free evolution of microseismic rupture. The fitting confidence weight is determined based on the length of the linear scaling region and the statistics of the absolute value of the second derivative of the data points within the linear scaling region. S3. Linear fitting is performed within the linear scaling region to obtain the curve slope, and the curve slope is used as the microseismic fractal scaling index. The microseismic fractal scaling index of the reconstructed phase space system is weighted and averaged using the fitting confidence weight to obtain the comprehensive fractal scaling index. When the rate of decrease and the magnitude of decrease of the comprehensive fractal scaling index within a continuous time window both exceed the set threshold, it is determined that a deep rock mass disaster is about to occur.

[0007] This invention constructs a modal saliency index by multiplying the inverse of the kurtosis value of each intrinsic mode function by its energy contribution rate. After normalization, this index is used as a weighting coefficient to reconstruct each intrinsic mode function. Components with low kurtosis values ​​and high energy contribution rates receive greater weights, while the weights of noise-dominant components are reduced. This improves the signal-to-noise ratio of the target signal while preserving micro-fracture precursor information. A grid search is performed within a two-dimensional parameter space defined by the delay time and embedding dimension. Phase space reconstruction is performed using each set of parameters in the parameter set. The fitting confidence weights for each set of parameters are obtained based on the length of the linear scaling region and the average absolute value of the second derivative of the data points. A weighted average is then used to obtain the comprehensive fractal scaling index. This reduces the nonlinear feature extraction bias caused by a single empirical parameter, and improves the computational stability of the comprehensive fractal scaling index.

[0008] Preferably, the step of acquiring deep rock mass microseismic monitoring signals and performing variational mode decomposition to obtain eigenmode functions reflecting the characteristics of different frequency bands of microseismic waves includes: Set the number of decomposition modes, the quadratic penalty factor, and the convergence tolerance for variational mode decomposition. Initialize the center frequency and eigenmode functions. Iteratively update each eigenmode function and its corresponding center frequency using the alternating direction multiplier method until the update error of each eigenmode function is lower than the convergence tolerance. Stop the iteration and output all eigenmode functions obtained by decomposition.

[0009] Preferably, obtaining the noise-reduced micro-vibration target signal includes: The fourth central moment of each intrinsic mode function is calculated and divided by the fourth power of the standard deviation of the intrinsic mode function to obtain the kurtosis value representing the transient impact intensity of the microseismic event. The reciprocal of the kurtosis value is obtained. The sum of squares of the amplitudes of each intrinsic mode function is calculated and divided by the sum of the sums of squares of the amplitudes of all intrinsic mode functions to obtain the microseismic radiation energy contribution rate. The reciprocal of the kurtosis value of each intrinsic mode function is multiplied by the microseismic radiation energy contribution rate to obtain the modal prominence index reflecting the effective rupture information of the microseismic event. The modal prominence index of each intrinsic mode function is divided by the sum of the modal prominence indices of all modal prominence indices to obtain the weighting coefficient. Each intrinsic mode function is multiplied by the corresponding weighting coefficient and then summed to obtain the noise-reduced microseismic target signal. The weight of the noise-dominant component is reduced in the weighted reconstruction.

[0010] Preferably, the step of generating a parameter set by performing a grid search within a preset delay time and embedding dimension space includes: A one-dimensional delay time series is generated by determining the lower and upper limits of the delay time and a fixed search step size. A one-dimensional embedding dimension series is generated by determining the lower and upper limits of the embedding dimension and a fixed search step size. A discrete grid matrix is ​​constructed by orthogonally permuting each value in the one-dimensional delay time series and each value in the one-dimensional embedding dimension series. All coordinate nodes of the discrete grid matrix are extracted into a parameter set composed of multiple combinations of delay time and embedding dimension.

[0011] Preferably, identifying the continuous intervals where the absolute value of the second derivative of the data points is lower than a preset threshold as linear scaling regions reflecting the scale-free evolution of microseismic rupture includes: resampling the abscissa of the double logarithmic coordinate curve at equal intervals, obtaining the second difference result for each discrete data point on the curve using the central difference algorithm, dividing the second difference result by the square of the abscissa step size to obtain the second derivative of the data point, taking the absolute value of the second derivative of each data point, comparing the absolute value of the second derivative of each data point with a set constant scalar linearity threshold in turn, extracting all data points whose absolute value is less than the constant scalar linearity threshold to form an initial set, searching for each connected component segment that is continuous from beginning to end without interruption within the initial set, calculating the number of data points contained in each connected component segment, extracting the single connected component segment containing the most data points, and labeling the lateral interval range between the data points at both ends of the connected component segment as a linear scaling region.

[0012] Preferably, determining the fitting confidence weights includes: calculating the absolute value of the difference between the x-coordinates of the data points at both ends of the linear scaling region as the length of the linear scaling region; counting the total number of discrete data points contained within the linear scaling region; summing the absolute values ​​of the second derivatives of all discrete data points within the linear scaling region and dividing by the total number of discrete data points to obtain the average value of the absolute values ​​of the second derivatives of the data points; calculating the reciprocal of the sum of the average value of the absolute values ​​of the second derivatives of the data points and a preset positive constant; multiplying the reciprocal by the length of the linear scaling region to obtain the initial value of the fitting confidence weights; dividing the initial value of the fitting confidence weights of each reconstructed phase space by the sum of the initial values ​​of the fitting confidence weights of all reconstructed phase spaces; and setting the output result as the fitting confidence weights.

[0013] Preferably, the determination that a deep rock mass disaster is about to occur includes: The maximum value and the current endpoint value of the comprehensive fractal scaling index within a continuous time window reflecting the microseismic disaster incubation period are extracted. The difference between the maximum value and the current endpoint value is used to obtain a numerical scalar as the decrease magnitude reflecting the tendency of microseismic sources to spatially cluster. Based on the least squares fitting algorithm, the comprehensive fractal scaling index sequence within the continuous time window is fitted with the time coordinates using a first-order linear fitting operation. The negative value of the slope parameter of the fitted equation is extracted as the decrease rate representing the speed of dimensional reduction and reconstruction of the rock mass fracture system. When the decrease magnitude is greater than the magnitude alarm threshold and the decrease rate is greater than the rate alarm threshold, an alarm output command is triggered and a deep rock mass disaster prediction report including rockburst or large deformation is generated to determine that a deep rock mass disaster is about to occur.

[0014] The magnitude of the decline and the rate of decline represent the decay state of the comprehensive fractal scaling index from the perspectives of cumulative amount and rate of change, respectively. An alarm is triggered when both indicators exceed the set threshold, which can reduce the probability of false alarms compared to judging by a single indicator.

[0015] Preferably, after acquiring the microseismic monitoring signal of the deep rock mass and before performing variational mode decomposition to obtain the intrinsic mode function, the method further includes mean removal and normalization preprocessing of the microseismic monitoring signal of the deep rock mass; the number of decomposition modes of the variational mode decomposition is determined in the following way: The number of decomposed modes is gradually increased starting from a small number. When the difference between the center frequency of the newly added intrinsic mode function and the center frequency of the existing intrinsic mode function is less than the preset frequency interval threshold, it is determined that frequency overlap has occurred, the increase is stopped, and the number of decomposed modes in the previous step is taken as the final value.

[0016] Preferably, a fixed-length sliding time window is used to process the deep rock mass monitoring signal data stream to calculate the comprehensive fractal scaling index in each time window. If the rate of decrease in multiple consecutive time windows is greater than the rate alarm threshold and the cumulative decrease of the comprehensive fractal scaling index in multiple consecutive time windows is greater than the amplitude alarm threshold, it is determined that a deep rock mass disaster is about to occur.

[0017] Secondly, the present invention provides a deep rock mass disaster prediction system, including a processor and a memory, wherein the memory stores computer program instructions, and when the computer program instructions are executed by the processor, the above-mentioned deep rock mass disaster prediction method is implemented.

[0018] By adopting the above technical solution, a computer program for predicting deep rock mass disasters is generated and stored in a memory so that it can be loaded and executed by a processor. A terminal device can then be made based on the memory and processor for convenient use.

[0019] The beneficial effects of this invention are as follows: This invention addresses the problems of low signal-to-noise ratio in deep rock mass monitoring signals, reliance on empirically selected phase space reconstruction parameters, and interference from subjective parameters in the fractal scaling index. It constructs a modal prominence index based on the product of the inverse of kurtosis and the energy contribution rate to perform weighted reconstruction of intrinsic mode functions, suppressing noise components while preserving micro-fracture precursor information, thus improving the signal-to-noise ratio of the target signal. Furthermore, it performs grid search within a two-dimensional parameter space composed of delay time and embedding dimension, and performs phase space reconstruction for each set of parameters. Combined with a fitting confidence weight constructed based on the length of the linear scaling region and the average absolute value of the second derivative, it performs a weighted average of the fractal scaling index, reducing the interference of subjective parameter selection on the calculation results and improving the stability of the comprehensive fractal scaling index.

[0020] Furthermore, by simultaneously monitoring the magnitude and rate of decrease of the comprehensive fractal scaling index within a continuous time window, an alarm is triggered and a disaster prediction report is generated when both dimensions exceed the set threshold. This establishes disaster determination based on the combined conditions of magnitude and rate, reducing the false alarm rate caused by single-indicator determination and improving the accuracy and advance warning time of deep rock mass disaster warnings. Attached Figure Description

[0021] Figure 1 A flowchart of a method for predicting disasters in deep rock masses; Figure 2 A schematic diagram comparing the reciprocal of the kurtosis value of each intrinsic mode function with its energy contribution rate; Figure 3 This is a diagram showing the performance indicators of different experimental schemes. Detailed Implementation

[0022] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0023] This invention discloses a method for predicting deep rock mass disasters, referring to... Figure 1 This includes steps S1-S3: S1. Microseismic signal decomposition and phase space reconstruction.

[0024] Time-series signals from deep rock mass microseismic monitoring were acquired using a microseismic sensor network, and the acquired signals underwent mean removal and normalization preprocessing. Variational mode decomposition (VMD) was then performed on the preprocessed signals. VMD is a signal processing method that decomposes an input signal into several narrowband eigenmode functions; its details will not be elaborated here.

[0025] Set the number of decomposition modes, the quadratic penalty factor, and the convergence tolerance for variational mode decomposition. Initialize the center frequency and eigenmode functions. Iteratively update each eigenmode function and its corresponding center frequency using the alternating direction multiplier method until the update error of each eigenmode function is lower than the convergence tolerance. Stop the iteration and output all eigenmode functions obtained by decomposition.

[0026] The number of decomposed modes is determined using the center frequency method. The criterion for the center frequency method is as follows: starting with a small number of decomposed modes, the number is gradually increased. When the difference between the center frequency of a newly added intrinsic mode function and the center frequency of an existing intrinsic mode function is less than a preset frequency interval threshold, frequency overlap is determined, the increase is stopped, and the previous number of decomposed modes is used as the final value. In this embodiment, the number of decomposed modes is set within the integer range of 4 to 8, and 6 is chosen in this embodiment. Considering that too small a number of decomposed modes will cause signal components from different frequency bands to be mixed in the same intrinsic mode function and cannot be separated, while too large a number of decomposed modes will introduce spurious components and increase the computational burden, choosing 6 can balance decomposition accuracy and computational efficiency. The secondary penalty factor is set to 2000 in this embodiment, and in other embodiments, it can be adjusted within the range of 1000 to 3000 according to the signal characteristics. When the second-order penalty factor is too small, the bandwidth constraints of each intrinsic mode function weaken, and spectral overlap between different modes is more likely to occur; when the second-order penalty factor is too large, the bandwidth constraints are too strong, which may truncate wideband signal components into narrowband fragments and lose effective information. The convergence tolerance is set to 1×10. -5 This is used to control the iteration termination condition. After completing the parameter configuration, each eigenmode function and center frequency is initialized to a zero matrix.

[0027] After entering the iterative solution phase, the alternating direction multiplier method updates each eigenmode function and its corresponding center frequency in the frequency domain. In each iteration, the frequency band center and amplitude-frequency distribution of the modes are adjusted according to the current residual signal and the spectral characteristics of each mode. The sum of the squares of the Fourier spectrum changes of each mode in the two iterations is calculated as the update error, and it is determined whether the update error is lower than the convergence tolerance of 1×10⁻⁶. -5 The loop stops when the update errors of all modes meet the termination condition, and all intrinsic mode functions obtained from the decomposition are output.

[0028] Existing variational mode decomposition (VMD) denoising methods use equal-weighted superposition of the retained intrinsic mode functions (IMFs), failing to distinguish the differences in contribution of each component to the target signal. This embodiment employs a differentiated weighting method. After obtaining all IMFs, the target signal is obtained by weighted reconstruction based on the product of the reciprocal of the kurtosis value and the energy contribution rate of each IMF. The fourth-order central moment of each IMF is calculated and divided by the fourth power of the IMF standard deviation to obtain the kurtosis value, and the reciprocal of the kurtosis value is calculated. The kurtosis value is a statistic used to measure the thickness and sharpness of the tail of the signal probability distribution; taking its reciprocal can reduce the contribution of high-kurtosis noise components in subsequent reconstruction. Taking a certain IMF as an example, if the calculated kurtosis value is 4.5, then the reciprocal of the kurtosis value is 0.222.

[0029] The energy contribution rate is obtained by calculating the sum of squared amplitudes of each intrinsic mode function (IMF) and dividing it by the sum of the squared amplitudes of all IMFs. The energy contribution rate reflects the proportion of a single mode in the total energy of the original signal and ranges from 0 to 1. The modal saliency index is obtained by multiplying the reciprocal of the kurtosis value of each IMF by the energy contribution rate. (Refer to...) Figure 2 Intrinsic mode functions 2 and 3 both have high reciprocals of kurtosis and energy contribution rates, corresponding to signal modes containing precursor information of micro-fractures in rock mass; intrinsic mode functions 5 and 6 have low values ​​for both indicators, mainly background noise components.

[0030] The weighting coefficient is obtained by dividing the modal saliency index of each intrinsic mode function by the sum of all modal saliency indices. The modal saliency index combines two statistical properties: the reciprocal of the kurtosis value and the energy contribution rate. The target signal is obtained by multiplying each intrinsic mode function by its corresponding weighting coefficient and then summing them.

[0031] A parameter set is generated by performing a grid search within a preset delay time and embedding dimension space. A one-dimensional delay time sequence is generated by determining the lower and upper limits of the delay time and a fixed search step size. Similarly, a one-dimensional embedding dimension sequence is generated by determining the lower and upper limits of the embedding dimension and a fixed search step size. In this embodiment, the lower limit of the delay time is set to 1, the upper limit to 15, and the fixed search step size to 2, generating an arithmetic sequence of delay times from 1, 3, 5 to 15. The lower limit of the embedding dimension is set to 3, the upper limit to 10, and the fixed search step size to 1, generating a continuous integer sequence from 3 to 10. The delay time ranges from 1 to 15, and the embedding dimension ranges from 3 to 10.

[0032] A discrete grid matrix is ​​constructed by orthogonally combining each value in the one-dimensional delay time series with each value in the one-dimensional embedding dimension series. All coordinate nodes of the discrete grid matrix are then extracted into a parameter set consisting of multiple combinations of delay time and embedding dimension. In this embodiment, the orthogonal permutation combination generates a parameter set containing 64 different parameter combinations.

[0033] The target signal is reconstructed using each set of parameters in the parameter set to obtain a reconstructed phase space system. Phase space reconstruction is based on Tukens' theorem; by setting specific delay times and embedding dimensions, a time delay matrix is ​​constructed to map a one-dimensional time series to a multi-dimensional state space. In this embodiment, 64 sets of parameters generate 64 reconstructed phase spaces, forming a reconstructed phase space system.

[0034] S2, Calculation of correlation integral and confidence weighting.

[0035] For any reconstructed phase space in the reconstructed phase space system, the Grassberger-Procaccia algorithm is used to calculate the correlation integral. The lower limit of the scale radius sequence is set to 5% of the minimum Euclidean distance between point pairs in the reconstructed phase space, and the upper limit is set to 20% of the maximum distance. Within this range, 100 scale radius values ​​are generated at equal logarithmic intervals. The Grassberger-Procaccia algorithm is a well-known method for estimating the correlation dimension of chaotic attractors, and will not be elaborated upon here. The proportion of all point pairs in the phase space whose distance is less than each scale radius is calculated as the correlation integral value.

[0036] By taking the natural logarithm of the scale radius sequence and the correlation integral values, a discrete data point array of logarithmic scale radius and logarithmic correlation integral is constructed, generating a double logarithmic coordinate curve.

[0037] Assuming that, under a certain set of delay time and embedding dimension parameters, the spacing of the original data points on the logarithmic scale radius axis of the reconstructed phase space double logarithmic coordinate curve varies from 0.02 to 0.15, it needs to be converted to an evenly spaced distribution. Setting the resampling step size to 0.05, the original data points are mapped to evenly spaced nodes using spline interpolation. After evenly resampling the abscissa of the double logarithmic coordinate curve, the central difference algorithm is used to obtain the second-order difference result for each discrete data point on the curve. Since the data points on the horizontal logarithmic scale of the double logarithmic coordinate curve may not be evenly distributed, spline interpolation is used to resample the data points into a uniformly distributed point sequence with a constant abscissa step size. After completing the evenly spaced resampling, the central difference scheme is used to obtain the second-order difference result for each evaluation point and its adjacent points, and the second-order difference result is divided by the square of the abscissa step size to obtain the second derivative of the data point. The second derivative of the data point represents the local curvature of the double logarithmic coordinate curve at each node.

[0038] The absolute value of the second derivative of each data point is taken, and then compared sequentially with a set constant scalar linearity threshold. It is worth noting that the setting of the constant scalar linearity threshold directly affects the recognition accuracy of the linear scaling region: a threshold that is too large may introduce data points from non-linear regions, leading to increased fitting deviation; a threshold that is too small may result in an overly narrow identified linear scaling region, causing the loss of effective fitted data. In this embodiment, the constant scalar linearity threshold is set to 0.1; in other embodiments, it can be adjusted within the range of 0.05 to 0.2 according to the signal characteristics.

[0039] Extract all data points whose absolute value is less than a constant scalar linearity threshold to form an initial set. Within this initial set, search for connected component segments that are consecutive and unbroken. A connected component segment is a subset of data points in the initial set whose adjacent indices are arranged consecutively without any breaks. Calculate the number of data points contained in each connected component segment, extract the single connected component segment containing the most data points, and define the horizontal interval between the data points at both ends of the connected component segment as the linear scaling region.

[0040] The fitting confidence weight is determined by multiplying the length of the linear scaling region by the average absolute value of the second derivative of the data points within the linear scaling region plus the reciprocal of a preset positive constant. The absolute value of the difference between the x-coordinates of the data points at both ends of the linear scaling region is calculated as the length of the linear scaling region. A typical dataset from a reconstructed phase space is used for verification: the logarithmic x-coordinate spans from -3.2 to -1.5. The length of the linear scaling region is... This length indicates that the target rock mass system maintains self-similarity within a scale range from -3.2 to -1.5, covering the main scale interval for microcrack propagation.

[0041] The total number of discrete data points contained within the linear scaling region is counted. The absolute values ​​of the second derivatives of all discrete data points within the linear scaling region are summed and divided by the total number of discrete data points to obtain the average value of the absolute values ​​of the second derivatives of the data points. In this embodiment, the linear scaling region of a certain reconstructed phase space contains 15 data points, and the average value of the absolute values ​​of the second derivatives of the data points is 0.02.

[0042] This invention constructs fitting confidence weights by combining the length of the linear scaling region and the curve smoothness to evaluate the reliability of fractal features under different reconstruction parameters. The reciprocal of the sum of the average absolute values ​​of the second derivatives of the data points and a preset positive constant is calculated. In this embodiment, the preset positive constant is set to 1 × 10⁻⁶. -6 This is used to ensure that the denominator is not zero. The smaller the average absolute value of the second derivative of the data points, the larger the reciprocal, indicating that the segment of the double logarithmic coordinate curve is closer to the ideal straight line. The reciprocal is multiplied by the length of the linear scaling region to obtain the initial value of the fitting confidence weight.

[0043] Divide the initial value of the fitting confidence weight for each reconstructed phase space by the sum of the initial values ​​of the fitting confidence weights for all reconstructed phase spaces, and set the output as the fitting confidence weight. After normalization, the sum of the fitting confidence weights for all reconstructed phase spaces is 1.

[0044] S3. Calculation of microseismic fractal scaling index and disaster assessment.

[0045] Log-scale radius and log-correlation integral data points within the linear scaling region are extracted, and least-squares linear regression is performed for fitting. The slope of the fitted curve is extracted and used as the fractal scaling exponent of the reconstructed phase space at a specific delay time and embedding dimension. Least-squares linear regression is a well-known method for fitting data points by minimizing the sum of squared residuals, which will not be elaborated here. The above operation is repeated for the reconstructed phase space system generated by the parameter set to obtain all fractal scaling exponents and their corresponding fit confidence weights.

[0046] The weighted average is calculated by multiplying each fractal scaling index by its corresponding fit confidence weight and summing the results, then dividing by the sum of all fit confidence weights to obtain the comprehensive fractal scaling index. It is worth noting that if the fractal scaling indices of each parameter group are directly averaged without weighting by fit confidence, parameter combinations with poor fit quality will receive the same weight as those with good fit quality, leading to increased volatility in the comprehensive fractal scaling index. Since the sum of the fit confidence weights after normalization is 1, the weighted average calculation is equivalent to directly multiplying each fractal scaling index by its fit confidence weight and then summing the results.

[0047] A fixed-length sliding time window is used to process the monitoring signal data stream of deep rock mass to calculate the comprehensive fractal scaling index within each time window. In this embodiment, the time window step size is set to 1 hour, and the window length is set to 24 hours. When a large number of microcracks initiate and connect within the rock mass, and the system tends towards macroscopic instability and failure, the comprehensive fractal scaling index will show a downward trend.

[0048] Extract the maximum and current endpoint values ​​of the composite fractal scaling index within a continuous time window. In the first time window, the first calculated value of the composite fractal scaling index sequence within that window is used as the initial value for both the maximum and current endpoint values. The difference between the maximum and current endpoint values ​​is taken as the decrease magnitude. Verification is performed using typical parameters: the composite fractal scaling index decreases from 2.15 in the stationary period to 1.85, a decrease of 0.3.

[0049] The least squares fitting algorithm is used to perform a first-order linear fitting operation between the comprehensive fractal scaling index sequence within a continuous time window and the time coordinates. The negative value of the slope parameter of the fitted equation is extracted as the descent rate. The reason for taking a negative value for the slope parameter is that the comprehensive fractal scaling index decreases as the disaster approaches, and the slope of the fitted line is negative. Taking a negative value turns it into a positive value, which is then used for comparison with the threshold.

[0050] It should be noted that the amplitude alarm threshold and rate alarm threshold are calibrated based on the critical failure mechanics experimental laws of regional geological strata. The calibration process is as follows: A time series of the comprehensive fractal scaling index prior to historical disaster events in the target mining area is collected. The decrease amplitude and rate of decrease of the comprehensive fractal scaling index within a continuous time window before each disaster event are statistically analyzed. The lower quartile of the historical samples is used as the initial threshold value, and then fine-tuned according to the balance requirements of the missed detection rate and false alarm rate. In this embodiment, the amplitude alarm threshold is set to 0.25, and in other embodiments, it can be adjusted within the range of 0.15 to 0.35 according to the geological characteristics of the mining area; the rate alarm threshold is set to 0.01 / h, and in other embodiments, it can be adjusted within the range of 0.005 / h to 0.02 / h. A smaller amplitude alarm threshold increases sensitivity but raises the false alarm rate; a larger threshold may miss true precursors.

[0051] An alarm is triggered and a disaster prediction report is generated when both the descent rate and magnitude alarm thresholds are greater than the magnitude alarm threshold, indicating an impending deep rock mass disaster. An alarm trigger signal is sent to the mine safety monitoring center system when the descent rate exceeds the rate alarm threshold and the cumulative descent magnitude exceeds the magnitude alarm threshold across multiple consecutive time windows. In this embodiment, an alarm is triggered when the descent rate exceeds 0.01 / h across three consecutive time windows and the cumulative descent magnitude exceeds 0.25.

[0052] Verification was conducted using actual operational data: When the decrease in the comprehensive fractal scaling index is 0.3, which is greater than the magnitude alarm threshold of 0.25, and the decrease rate is 0.015 / h, which is greater than the rate alarm threshold of 0.01 / h, the system activates a hardware relay to trigger an audible and visual alarm output command, assembles the precursor indicator data and the dangerous time interval information, and exports a structured disaster prediction report.

[0053] The experiment used a real continuous dataset of rock mass dynamics collected by a microseismic monitoring system in a deep mining area, containing 500 sets of precursor signals of rock mass microfractures and conventional background noise samples. Three control schemes were set up for comparison: the first scheme was the complete scheme, including weighted reconstruction based on the inverse of kurtosis and energy contribution rate according to variational mode decomposition, phase space grid search, and a confidence-weighted fitting mechanism based on the second derivative of data points; the second scheme was a control scheme without the weighted reconstruction module, where all intrinsic mode functions obtained from decomposition were added with equal weights to restore the signal; the third scheme was a control scheme without the phase space grid search and confidence-weighted modules, using only a fixed delay time of 5 and an embedding dimension of 6 to calculate a single fractal scaling exponent and perform a double threshold judgment.

[0054] The experiment extracted the target signal signal-to-noise ratio, disaster early warning accuracy, and average early warning time as evaluation indicators. (Refer to...) Figure 3 The performance index comparison results of the three schemes are as follows: Figure 3 As shown in the figure, regarding the signal-to-noise ratio (SNR) of the target signal, the first complete scheme achieved an SNR of 28.5 dB, the second comparative scheme reduced the reconstructed signal SNR to 19.2 dB, and the third comparative scheme achieved an SNR of 20.1 dB. Regarding disaster early warning accuracy, the first complete scheme achieved a disaster early warning accuracy of 96.2%, the second comparative scheme reduced it to 83.4%, and the third comparative scheme achieved 78.5%. Regarding average early warning time, the first complete scheme achieved an average early warning time of 14.5 hours, the second comparative scheme shortened it to 9.2 hours, and the third comparative scheme achieved 7.8 hours. The three comparisons demonstrate that modal weighted reconstruction and multi-parameter optimization mechanisms improve signal quality, prediction accuracy, and early warning timeliness.

[0055] This invention also discloses a deep rock mass disaster prediction system, including a processor and a memory. The memory stores computer program instructions, which, when executed by the processor, implement a deep rock mass disaster prediction method according to the present invention.

[0056] The system also includes other components well known to those skilled in the art, such as communication buses and communication interfaces, the settings and functions of which are known in the art and will not be described in detail here.

[0057] In the description of this specification, "multiple" or "several" means at least two, such as two, three or more, unless otherwise expressly and specifically defined.

Claims

1. A method for predicting disasters in deep rock masses, characterized in that, include: S1. Acquire deep rock mass microseismic monitoring signals and perform variational mode decomposition to obtain intrinsic mode functions reflecting the characteristics of different frequency bands of microseismic waves. Based on the product of the reciprocal of the kurtosis value representing the microseismic impact characteristics of each intrinsic mode function and the energy contribution rate representing the microseismic radiation energy, perform weighted reconstruction to obtain the denoised microseismic target signal. Perform grid search within the preset delay time and embedded dimension space to generate a parameter set. Use the parameter set to reconstruct the phase space of the denoised microseismic target signal to obtain the reconstructed phase space system representing the rock mass damage evolution. S2. Calculate the correlation integral of any reconstructed phase space in the reconstructed phase space system at different scale radii to generate a double logarithmic coordinate curve. Calculate the second derivative of the data points on the double logarithmic coordinate curve and identify continuous intervals where the absolute value of the second derivative of the data points is lower than a preset threshold as linear scaling regions reflecting the scale-free evolution of microseismic rupture. Determine the fitting confidence weight based on the length of the linear scaling region and the statistics of the absolute value of the second derivative of the data points within the linear scaling region. S3. Perform linear fitting within the linear scaling region to obtain the curve slope and use the curve slope as the microseismic fractal scaling index. Calculate the comprehensive fractal scaling index by weighted averaging of the microseismic fractal scaling index of the reconstructed phase space system using the fitting confidence weight. When the rate of decrease and the magnitude of decrease of the comprehensive fractal scaling index within a continuous time window both exceed the set threshold, it is determined that a deep rock mass disaster is about to occur.

2. The method for predicting deep rock mass disasters according to claim 1, characterized in that, The process of acquiring deep rock mass microseismic monitoring signals and performing variational mode decomposition to obtain intrinsic mode functions reflecting the characteristics of different frequency bands of microseismic waves includes: Set the number of decomposition modes, the quadratic penalty factor, and the convergence tolerance for variational mode decomposition. Initialize the center frequency and eigenmode functions. Iteratively update each eigenmode function and its corresponding center frequency using the alternating direction multiplier method until the update error of each eigenmode function is lower than the convergence tolerance. Stop the iteration and output all eigenmode functions obtained by decomposition.

3. The method for predicting deep rock mass disasters according to claim 1, characterized in that, The obtained noise-reduced micro-vibration target signal includes: The fourth-order central moment of each intrinsic mode function is calculated and divided by the fourth power of the standard deviation of the intrinsic mode function to obtain the kurtosis value representing the transient impact intensity of the microseismic event. The reciprocal of the kurtosis value is obtained. The sum of squares of the amplitudes of each intrinsic mode function is calculated and divided by the sum of the sums of the squares of the amplitudes of all intrinsic mode functions to obtain the microseismic radiation energy contribution rate. The reciprocal of the kurtosis value of each intrinsic mode function is multiplied by the microseismic radiation energy contribution rate to obtain the modal prominence index reflecting the effective rupture information of the microseismic event. The modal prominence index of each intrinsic mode function is divided by the sum of the modal prominence indices of all modal prominence indices to obtain the weighting coefficient. Each intrinsic mode function is multiplied by the corresponding weighting coefficient and then summed to obtain the noise-reduced microseismic target signal.

4. The method for predicting deep rock mass hazards according to claim 1, characterized in that, The step of generating a parameter set by performing a grid search within a preset delay time and embedding dimension space includes: A one-dimensional delay time series is generated by determining the lower and upper limits of the delay time and a fixed search step size. A one-dimensional embedding dimension series is generated by determining the lower and upper limits of the embedding dimension and a fixed search step size. A discrete grid matrix is ​​constructed by orthogonally permuting each value in the one-dimensional delay time series and each value in the one-dimensional embedding dimension series. All coordinate nodes of the discrete grid matrix are extracted into a parameter set composed of multiple combinations of delay time and embedding dimension.

5. The method for predicting deep rock mass disasters according to claim 1, characterized in that, The step of identifying continuous intervals where the absolute value of the second derivative of a data point is lower than a preset threshold as linear scaling regions reflecting the scale-free evolution of microseismic rupture includes: resampling the abscissa of the double logarithmic coordinate curve at equal intervals, then using the central difference algorithm to obtain the second difference result for each discrete data point on the curve, dividing the second difference result by the square of the abscissa step size to obtain the second derivative of the data point, taking the absolute value of the second derivative of each data point, comparing the absolute value of the second derivative of each data point with a set constant scalar linearity threshold, extracting all data points whose absolute value is less than the constant scalar linearity threshold to form an initial set, searching for each connected component segment that is continuous from beginning to end without interruption within the initial set, calculating the number of data points contained in each connected component segment, extracting the single connected component segment containing the most data points, and labeling the lateral interval range between the data points at both ends of the connected component segment as the linear scaling region.

6. The method for predicting deep rock mass disasters according to claim 1, characterized in that, The determination of the fitting confidence weights includes: calculating the absolute value of the difference between the x-coordinates of the data points at both ends of the linear scaling region as the length of the linear scaling region; counting the total number of discrete data points contained within the linear scaling region; summing the absolute values ​​of the second derivatives of all discrete data points within the linear scaling region and dividing by the total number of discrete data points to obtain the average value of the absolute values ​​of the second derivatives of the data points; calculating the reciprocal of the sum of the average value of the absolute values ​​of the second derivatives of the data points and a preset positive constant; multiplying the reciprocal by the length of the linear scaling region to obtain the initial value of the fitting confidence weights; dividing the initial value of the fitting confidence weights of each reconstructed phase space by the sum of the initial values ​​of the fitting confidence weights of all reconstructed phase spaces; and setting the output result as the fitting confidence weights.

7. The method for predicting deep rock mass disasters according to claim 1, characterized in that, The determination that a deep rock mass disaster is about to occur includes: The maximum value and the current endpoint value of the comprehensive fractal scaling index within a continuous time window reflecting the microseismic disaster incubation period are extracted. The difference between the maximum value and the current endpoint value is used to obtain a numerical scalar as the decrease magnitude reflecting the tendency of microseismic sources to spatially cluster. Based on the least squares fitting algorithm, the comprehensive fractal scaling index sequence within the continuous time window is fitted with the time coordinates using a first-order linear fitting operation. The negative value of the slope parameter of the fitted equation is extracted as the decrease rate representing the speed of dimensional reduction and reconstruction of the rock mass fracture system. When the decrease magnitude is greater than the magnitude alarm threshold and the decrease rate is greater than the rate alarm threshold, an alarm output command is triggered and a deep rock mass disaster prediction report including rockburst or large deformation is generated to determine that a deep rock mass disaster is about to occur.

8. The method for predicting deep rock mass disasters according to claim 1, characterized in that, After acquiring the microseismic monitoring signals of deep rock masses and before performing variational mode decomposition to obtain the intrinsic mode functions, the process also includes mean removal and normalization preprocessing of the deep rock mass microseismic monitoring signals; the number of decomposition modes in the variational mode decomposition is determined in the following way: The number of decomposed modes is gradually increased starting from a small number. When the difference between the center frequency of the newly added intrinsic mode function and the center frequency of the existing intrinsic mode function is less than the preset frequency interval threshold, it is determined that frequency overlap has occurred, the increase is stopped, and the number of decomposed modes in the previous step is taken as the final value.

9. A method for predicting deep rock mass hazards according to claim 7, characterized in that, A fixed-length sliding time window is used to process the deep rock mass monitoring signal data stream to calculate the comprehensive fractal scaling index within each time window. If the rate of decrease in multiple consecutive time windows is greater than the rate alarm threshold and the cumulative decrease in the comprehensive fractal scaling index in multiple consecutive time windows is greater than the amplitude alarm threshold, it is determined that a deep rock mass disaster is about to occur.

10. A deep rock mass disaster prediction system, characterized in that, include: The processor and memory, the memory storing computer program instructions, implement a deep rock mass disaster prediction method according to any one of claims 1 to 9 when the computer program instructions are executed by the processor.

Images