Multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR and AWMF

Abstract

The application discloses a multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR and AWMF, relates to the technical field of microseismic signal processing, and comprises the following steps: acquiring an original signal; decomposing the original signal to obtain a plurality of intrinsic mode function (IMF) components; dividing the IMF components into a first component, a second component and a third component; performing denoising on the first component, the second component and the third component, and reconstructing to obtain a preliminary denoised signal; reconstructing the preliminary denoised signal in phase space and denoising; filtering the IMF components; analyzing an event window, and smoothing a non-event window to obtain a final denoised signal; the VMD decomposition process is automatic and reproducible, and the signal-to-noise ratio and reconstruction accuracy are significantly improved on a plurality of types of synthetic signals and field microseismic signals; by combining phase space reconstruction and local SVD, residual coherent noise and random disturbance are effectively suppressed, and the main lobe structure and chaotic characteristics of the microseismic signal are maintained; and by wavelet-morphological joint filtering, fine suppression of high-frequency noise and ringing is realized.

Description

Technical Field

[0001] This invention relates to the field of microseismic signal processing technology, and in particular to a multi-stage event-holding microseismic signal noise reduction method based on WOA-VMD, PSR and AWMF. Background Technology

[0002] With the advancement of major infrastructure projects in Southwest China in recent years, the geological environment and rock mass occurrence at these project sites have become increasingly extreme and complex. Under extreme conditions such as complex geological structures, high ground stress, and high temperature and pressure, the combined effects of excavation disturbances make the energy accumulation and release process in rock masses more intense and unpredictable, leading to frequent dynamic disasters such as rockbursts. In recent years, microseismic monitoring technology, with its high sensitivity and real-time performance, has become an important technical means for monitoring rock mass fracturing processes and providing early warning of rockburst disasters. However, in practical applications in underground engineering such as tunnels, the raw data is severely interfered with by noise from complex environments and instrument current, resulting in a low signal-to-noise ratio in the acquired microseismic signals. This directly affects the reliability of subsequent feature extraction, spatial positioning, and source mechanism inversion, posing a greater challenge to the accuracy of disaster early warning at engineering sites. While single noise reduction methods can suppress noise components in the signal to some extent, they are prone to signal detail loss, transient information weakening, and noise residue in complex noise structures. Therefore, research on adaptive noise reduction methods with high fidelity and robustness for microseismic signals is crucial for improving the accuracy of monitoring and early warning of rockburst disasters in underground engineering. Summary of the Invention

[0003] The purpose of this invention is to design a multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF to solve the above problems.

[0004] The present invention achieves the above objectives through the following technical solutions:

[0005] A multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF is characterized by comprising:

[0006] S1. Obtain the original signal s raw (t), and perform preprocessing;

[0007] S2. Optimize the key parameters of variational mode decomposition (VMD), and use the optimized key parameters to perform variational mode decomposition (VMD) on the preprocessed original signal to obtain multiple intrinsic mode function (IMF) components.

[0008] S3. Divide all IMF components into the first component, the second component, and the third component. The first component is the component dominated by the effective signal, the second component is the component that is a mixture of signal and noise, and the third component is the component dominated by the noise signal.

[0009] S4. The first component, the second component, and the third component are denoised using the first denoising strategy, the second denoising strategy, and the third denoising strategy, respectively. The denoised first component, the second component, and the third component are then reconstructed to obtain the initial denoised signal v(t).

[0010] S5. Reconstruct the phase space of the initial denoised signal v(t); and within the local neighborhood of the reconstructed phase space, use singular value decomposition (SVD) to separate the principal component and noise component of the reconstructed initial denoised signal, remove the noise component, and retain the principal component as the IMF component w after secondary denoising. k (t);

[0011] S6, regarding IMF components w k (t) Perform adaptive wavelet-morphological filtering to obtain the filtered signal;

[0012] S7. Analyze the event window in the filtered signal, smooth the non-event window, and obtain and output the final noise-reduced signal.

[0013] The beneficial effects of this invention are as follows: By introducing a unified, dimensionless joint fitness function and employing the whale optimization algorithm to adaptively optimize the core parameters of VMD, it replaces the traditional empirical parameter setting method, automating and reproducible the VMD decomposition process, and significantly improving the signal-to-noise ratio and reconstruction accuracy for multiple types of synthetic signals and field microseismic signals; through an IMF classification strategy that integrates multiple features such as energy, correlation, Gaussian p-value, and sample entropy, it achieves transparent and interpretable partitioning of signal / mixed / noise modes, which is more stable and applicable to different noise scenarios compared to traditional IMF selection rules based on experience or manual judgment; by combining phase space reconstruction with local SVD and applying... The constraints of accumulated energy, main frequency consistency, and spectral entropy effectively suppress residual coherent noise and random disturbances, while maintaining the main lobe structure and chaotic characteristics of the microseismic signal (the maximum Lyapunov exponent change is less than about 7%, and the recurrence rate and other indicators remain within the original range). Through a wavelet-morphological joint filtering framework, the wavelet basis, decomposition level, threshold intensity, and morphological structural elements are adaptively selected under the triple driving forces of IMF category, event interval, and SNR level. This achieves fine suppression of high-frequency noise and ringing without excessively smoothing event edges and the main lobe. Experiments show that the spectral entropy can be reduced by about 9–21%, the zero-crossing rate by about 33–91%, and the residual high-frequency energy by about 69–85%. Attached Figure Description

[0014] Figure 1 The processing flow of a multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF;

[0015] Figure 2 Comparison of time-domain characteristics of the first group of microseismic signals before and after noise reduction;

[0016] Figure 3 Comparison of time-domain characteristics of the second group of microseismic signals before and after noise reduction;

[0017] Figure 4 Comparison of time-domain characteristics of the first group of microseismic signals before and after noise reduction;

[0018] Figure 5 Comparison of time-domain characteristics of the second group of microseismic signals before and after noise reduction;

[0019] Figure 6 Comparison of phase space probability density distribution before and after denoising of microseismic signals;

[0020] Figure 7 Comparison of phase space dynamic vector fields before and after denoising of microseismic signals;

[0021] in, Figure 1 , Figure 2 , Figure 3 , Figure 4 and Figure 5 In the diagram, (a) represents the original signal and (b) represents the denoised signal. Figure 6 and Figure 7 (a) is group 1 and (b) is group 2. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of 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, not all, of the embodiments of the present invention. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0023] Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0024] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0025] In the description of this invention, it should be understood that the terms "upper," "lower," "inner," "outer," "left," "right," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship commonly used when the product of this invention is in use, or the orientation or positional relationship commonly understood by those skilled in the art. They are only used to facilitate the description of this invention and to simplify the description, and are not intended to indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0026] Furthermore, the terms "first," "second," etc., are used only to distinguish descriptions and should not be interpreted as indicating or implying relative importance.

[0027] In the description of this invention, it should also be noted that, unless otherwise explicitly specified and limited, terms such as "set" and "connection" should be interpreted broadly. For example, "connection" can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium; it can be a connection within two components. Those skilled in the art can understand the specific meaning of the above terms in this invention according to the specific circumstances.

[0028] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0029] like Figure 1 As shown, a multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF includes:

[0030] S1. Obtain the original signal s raw (t) is then preprocessed; the preprocessing specifically includes DC removal and normalization. Alternatively, the original signal s can be modified based on the sensor's frequency response and engineering experience. raw (t) Apply a broadband bandpass filter (e.g., 20–2000Hz) to avoid extremely low frequency drift and high frequency noise interference that far exceeds the sensor bandwidth. The preprocessed signal is denoted as s(t).

[0031] S2. Optimize the key parameters of Variational Mode Decomposition (VMD), and use the optimized key parameters to perform VMD on the preprocessed original signal to obtain multiple Intrinsic Mode Function (IMF) components; specifically including:

[0032] S21. Initialize the population of the Whale Optimization Algorithm (WOA), set the position vector of each whale, and the position vector corresponds to the combination of parameters to be optimized in the Submodal Decomposition (VMD) [K, α], where K is the number of modes and α is the penalty factor;

[0033] S22. Simultaneously considering reconstruction error and signal-to-noise ratio, define a fitness function f(θ), using the reconstruction coefficients and signal-to-noise ratio of the decomposed IMF components as evaluation indicators, expressed as: ,in, , θ=[K,α,τ], The reconstruction error is represented as: The signal-to-noise ratio is calculated in dB and then mapped to [0,1], expressed as: , , and These are the weighting coefficients for the reconstruction error term and the signal-to-noise ratio term, respectively. For signal-to-noise ratio loss, Here is the parameter vector to be optimized; τ is the update step size, used to control the parameter update magnitude during the WOA search process; ε is a positive number to prevent the denominator from approaching 0, clip(x,0,1)=min(max(x,0),1), corresponding to the linear calibration of SNR∈[−20,80]dB; in this embodiment, we take... =0.6, =0.4, primarily considering reconstruction error while also taking into account high SNR. Signal-to-noise ratio in decibels. To be After linear calibration, then through the function The dimensionless signal-to-noise ratio index obtained by restricting it to the interval [0,1].

[0034] S23. Perform WOA (Whale Optimization Algorithm) iterative optimization. By simulating the contraction and spiral predation behavior of humpback whales, search for the optimal θ and finally output the optimal parameter θ* to obtain the optimal parameter combination [K,α] of the fitness function.

[0035] S24. Using the optimal parameter combination [K, α], perform variational mode decomposition (VMD) on the preprocessed original signal, decomposing it into K eigenmode functions u. k (t), each u k (t) represents the narrowband AM-FM component, with a center angular frequency of ω. k VMD is achieved by minimizing the sum of the bandwidths of all modes, expressed as: Furthermore, by constructing the augmented Lagrangian function and using the alternating direction multiplier method (ADMM) iterative solution, the eigenmode functions u0 were obtained. k (t) and its center frequency ω k The intrinsic mode functions u k (t) represents the IMF components distributed from high frequency to low frequency as: IMF1, IMF2, ..., IMFK.

[0036] S3. By analyzing the characteristic indicators of each IMF component, all IMF components are divided into three components: the first component, the second component, and the third component. The first component is dominated by effective signals, the second component is a mixture of signal and noise, and the third component is dominated by noise signals. Specifically, it includes:

[0037] S31. Calculate the characteristic indices of the IMFi components, including the relative energy E. i The absolute value of the Pearson correlation coefficient r with the original signal, and the sample entropy S. i (m=2, r=0.2·σ) i ) and Gaussian property value p i For N≤5000, the Gaussianity value of each IMF component is calculated using the Shapiro-Wilk test. For other cases, the Gaussianity value of each IMF component is calculated using the D'Agostino K² test. Gaussian properties Characterizing whether it is closer to Gaussian white noise, Gaussian property value This is used to characterize the degree to which each IMF component deviates from Gaussian white noise, and is used in subsequent comprehensive score calculations. They participate in IMF classification in the form of [these methods].

[0038] S32. Normalize the feature indicators; specifically, normalize each of the above feature indicators in the IMF dimension, that is, perform 0-1 linear normalization on the same feature value sequence corresponding to all IMF components.

[0039] S33. Analyze the comprehensive score based on the normalized feature indicators. i , represented as: Score i =Ê i +r̂ i +(1−p̂ i )−0.5·Ŝ i Ê i 、r̂ i 、p̂ i and Ŝ i The normalized relative energy, absolute value of Pearson correlation coefficient, sample entropy, and Gaussian p-value are respectively assigned a weight of 1 for relative energy, absolute value of Pearson correlation coefficient, and non-Gaussian p-value, and a weight of -0.5 for sample entropy to enhance the structured components.

[0040] S34. Based on the overall score i K-means clustering (K) of IMFi c =3), obtain the clustering results, sort the clustering results according to the average score, and map the clustering results to the first component, second component or third component.

[0041] S4. The first component, the second component, and the third component are denoised using the first denoising strategy, the second denoising strategy, and the third denoising strategy, respectively. The denoised first component, the second component, and the third component are then reconstructed to obtain the initial denoised signal v(t).

[0042] The first noise reduction strategy is as follows: when the IMF belongs to the first component (signal-dominant), it is generally retained entirely; if it is necessary to suppress a small amount of residual noise, a mild threshold is set for its envelope energy or amplitude. For example, take the average absolute amplitude of the component. Times: Only the local amplitude is lower than Weakening only occurs when the event window is not active;

[0043] The second noise reduction strategy is to apply adaptive wavelet thresholding to the second component, which is a mixture of signal and noise. Specifically, based on the local signal-to-noise ratio estimation result of this component, the number of wavelet decomposition layers is adaptively selected (not exceeding 4 layers), and the noise standard deviation of each layer is estimated using the MAD method. Press again Calculate the soft threshold, where λ(⋅) is a coefficient function that monotonically decreases with the local signal-to-noise ratio (e.g., linear interpolation within [−5,20] dB), thereby using a larger threshold in high-noise regions and a smaller threshold in event regions.

[0044] The third noise reduction strategy is to employ strong threshold suppression for the third component, which is predominantly noise; that is, a strategy close to "hard zeroing". Specifically, a threshold can be set. , ,like If noise dominates and there are no obvious signs of an event, the entire third component can be set to zero.

[0045] S5. Reconstruct the phase space of the initial denoised signal v(t); and within the local neighborhood of the reconstructed phase space, use singular value decomposition (SVD) to separate the principal component and noise component of the reconstructed initial denoised signal, remove the noise component, and retain the principal component as the IMF component w after secondary denoising. k (t); specifically including:

[0046] S51. Calculate the auto-mutual information of the initial denoised signal v(t), and find the first significant local minimum lag as a candidate time delay ζ; if there is no significant minimum within the preset lag range [1, N / 10], then use the lag point where the autocorrelation function R(ζ) decays to R(0) / e as the candidate time delay ζ; if the dominant frequency f̂ is known, the constraint can be added: ζ ≥ round(F s / f̂) / 4, to avoid aliasing; "significant local minimum" refers to... Within the search range, mutual information curve The first local minimum point appears, and the function value at this point decreases by at least 5%–10% compared to the previous local maximum.

[0047] S52. The pseudo nearest neighbor (FNN) method is used to select the embedding dimension m. The dimension is increased from 2 to an upper limit (e.g., 12). Specifically, when the proportion of FNN (m) is lower than a certain threshold (e.g., 1–3%), and the difference between FNN (m) and FNN (m+1) is not large (e.g., ≤0.5%), this dimension is determined as the current embedding dimension m. Furthermore, m can be appropriately adjusted for different frequency bands; for example, for low-frequency components, m can be set to a lower value. low =min(m+4,20), for high-frequency components, take m high =max(m−4,6);

[0048] S53. Using a sliding window, the initial noise-reduced signal v(t) of length N is divided into multiple segments of length W, where W can be 0.1N and is limited to [20, N / 8]. This step can effectively improve local stability and resolution.

[0049] S54. Construct the trajectory matrix X∈R for all sub-segments. L×m , represented as: Where L = N−(m−1)ζ;

[0050] S55. Perform singular value decomposition (SVD) on each trajectory matrix X, as follows: , where Σ is a diagonal singular value matrix, and the singular values ​​are sorted in descending order;

[0051] S56. Determine the truncation quantity r, retain the first r singular components, represented as: Specifically, this involves calculating the cumulative energy ratio of the first r singular values, and determining the ratio when it exceeds a preset threshold η. region Stop at time, η region Based on the time window region category: Event main window region: η E ≈0.98η; P-wave and other leader regions: η P ≈η; Background noise region: η B ≈0.90η; where η can be taken as approximately 0.95; to ensure the consistency constraint of the dominant frequency, the dominant frequency f is calculated for both the reconstructed local signal and the original local signal. domIf the deviation between the two exceeds a predetermined ratio (e.g., 10–15%), then r is appropriately increased until the main frequency consistency is met or r reaches the preset maximum value; compare the spectral entropy of the local signal before and after reconstruction. If the spectral entropy increases significantly after reconstruction, it indicates that the energy distribution is more diffuse. Then r is appropriately reduced to avoid the main lobe energy leaking out of the band.

[0052] S57, retain the singular component X denoise The time series is projected back to one dimension by diagonal mean, and the singular components X of all segments are retained. denoise The IMF components are spliced ​​or overlapped on the time axis and then normalized to obtain the IMF components after secondary denoising. k (t).

[0053] S6, regarding IMF components w k (t) Perform adaptive wavelet-morphological filtering to obtain the filtered signal; specifically including:

[0054] S61, for the IMF components after secondary denoising w k (t) Perform wavelet decomposition and obtain wavelet detail coefficients at different scales. The total number of wavelet decomposition levels is L and the number of decomposition levels is ℓ. Specifically, for high-frequency IMFs (center frequency close to half the sampling frequency): preferentially select wavelet bases such as db2, db4, and sym4 that have compact support and good time-domain resolution. The upper limit of the number of decomposition levels is L. cap =2; For intermediate frequency IMF: Near-symmetric wavelet bases such as sym5 and sym8 are preferred, L cap =3; For the entire microseismic record: L cap ≤5, the actual number of layers is min(dwt_max_level(N,ψ),L cap );

[0055] S62. For each detail coefficient dℓ, calculate the MAD noise estimate σ. ℓ , represented as: σ ℓ =median(|dℓ|) / 0.6745, is the median of the absolute values ​​of the detail coefficients for this layer, and 0.6745 is the scaling factor for converting the median absolute deviation to the standard deviation. Estimate the standard deviation of the noise in the details of this layer;

[0056] S63. Calculate the filter threshold T ℓ , represented as: T ℓ =α ℓ ·σ ℓ ; where threshold coefficient The value is related to the IMF category and decomposition level. The total number of layers L and whether the time position is within the event window are related factors. For example, for low-level details of signal-type IMFs, a smaller value is used. (e.g., 0.25), a soft threshold is used to preserve edges; for mid-frequency IMFs, a soft threshold strategy with layer depth increasing is adopted: T ℓ =0.9· •(0.7+0.2·ℓ / L); A larger hard threshold is used for high-frequency IMFs (e.g. ≈1.4), enhance high-frequency noise suppression; for global branches, appropriately reduce the threshold for low-frequency details within the event window (e.g., ∈{0.4,0.8}), while increasing the threshold outside the event and at higher-level details (e.g. ≈1.5). The soft / hard threshold rules described above are applied to the detail coefficients, while the approximation coefficients are generally left unprocessed or only slightly smoothed.

[0057] S64. Based on the filtering threshold T ℓ For the IMF components after secondary denoising w k (t) is used for filtering and noise reduction;

[0058] S65, The IMF components after secondary denoising... k (t) Perform one-dimensional grayscale mathematical morphology operations, design dilation and erosion operators, and base them on the IMF components w k (t) Feature selection: single opening / closing or cascaded opening / closing operations yield the morphological filtering result m. k (t); the opening operation performs erosion followed by expansion, denoted as , represented as: The closing operation involves performing expansion followed by erosion, denoted as... , represented as: One-dimensional grayscale mathematical morphological operations further suppress spike noise and fill local troughs. Opening operations can effectively eliminate positive pulses smaller than the structuring element while preserving the main structural features of the signal. Closing operations can effectively eliminate negative pulses smaller than the structuring element while filling small troughs in the signal. Specifically, based on the frequency band of the IMF and the sampling frequency f... s Choose a flat set of structure element lengths: High-frequency IMF: primarily open operation, structure element length set {3,5}; Mid-frequency IMF: open-close combination (open then closed or closed then open), length set {3,5,7}; Low-frequency IMF: primarily closed operation, length set {5,9,13}; When the sampling frequency is greater than 1kHz, it can be proportional to f. s / 1000 is a linear amplification of the above length, rounded to an odd number to maintain symmetry;

[0059] Let the structuring element be g, and the morphological verification signal be x(t). Erosion calculates the minimum value within a local window, thereby "shrinking" the signal and eliminating positive impulses. The dilation operation calculates the maximum value within a local window, thereby "expanding" the signal and eliminating negative impulses. ;

[0060] S66, For each IMF component w k (t), the combined IMF component after secondary denoising w k (t), Filtering and noise reduction results y k (t) and morphological branch output m k (t), perform multi-branch fusion to obtain the fused IMF component z k (t), denoted as: z k (t)=β1w k (t)+β2y k (t)+β3m k (t), where β1+β2+β3=1, and adaptively adjusted according to IMF category, event window position and local SNR: β1 and β2 have larger weights near the main lobe of the event to ensure waveform fidelity; β2 and β3 have larger weights at times dominated by high-frequency noise to enhance noise reduction.

[0061] S67. All IMF components z k Summing (t) yields the processed, denoised micro-vibration signal ŝ(t).

[0062] S7. Analyze the filtered signal to identify the event window and smooth the non-event window to obtain and output the final denoised signal. Specifically, based on the envelope energy, short-time energy, and amplitude threshold of the denoised microseismic signal ŝ(t), identify the start and end times of the microseismic events to form an event window, and record the position and duration of the main lobe of the event. Smoothing is performed in the non-event window, while over-smoothing is avoided in the event window, thereby ensuring that the output signal retains the key signal events without distortion while maintaining the uniformity of the energy criterion, and thus obtaining and outputting the final denoised signal.

[0063] In engineering scenarios where no real "clean signal" is available, to objectively evaluate the effectiveness of the noise reduction method of this invention, this embodiment uses no-reference statistical indicators to compare and analyze the microseismic signals before and after noise reduction. Specifically, within the event window and the entire recording range, the following four types of statistical characteristics of the original signal and the noise-reduced signal are calculated and their changes are compared: spectral entropy SE: reflects spectral complexity; kurtosis K: reflects peak sharpness and peak prominence; zero-crossing rate ZCR: reflects high-frequency fluctuations and noise content; standard deviation SD: reflects the amplitude of amplitude variation and whether it is overly smoothed to a certain extent.

[0064] Building upon the aforementioned evaluation based on non-reference statistical indicators such as spectral entropy, kurtosis, zero-crossing rate, and standard deviation, this embodiment further verifies the ability of the present invention to preserve the intrinsic dynamic characteristics of microseismic signals during noise reduction. This is achieved by jointly evaluating dynamic indicators related to phase space and recurrence analysis. Specifically, these include: maximum Lyapunov exponent, attractor complexity, recurrence rate, determinism, average length of the diagonal of the recurrence plot, trajectory smoothness, and trajectory regularity.

[0065] Among them, the maximum Lyapunov exponent and attractor complexity mainly reflect the degree of chaos and the phase space attractor structure of the system; recurrence rate and determinism characterize the recurrence behavior and orderliness of the trajectory in phase space; and average diagonal length, trajectory smoothness, and trajectory regularity measure the continuity and regularity of the trajectory from the perspective of temporal geometry. By comparing the changes in the above dynamic indicators between the original signal and the denoised signal, it can be determined whether the denoising algorithm, while suppressing noise, destroys the inherent nonlinear dynamic characteristics of the microseismic signal.

[0066] In the embodiments of this application, to verify the applicability and superiority of the proposed WOA-VMD-PSR-AWMF multi-stage noise reduction method under actual engineering conditions, two sets of typical microseismic acceleration records collected by a microseismic monitoring system deployed at the tunnel engineering site were selected as test signals. The above signals were sampled at a frequency of 4kHz and were acquired under complex underground environments and construction disturbances, containing significant environmental and instrument noise components. The overall signal-to-noise ratio was low, and these signals can represent typical working conditions of actual underground engineering microseismic monitoring data.

[0067] As shown in Table 1, after processing with the WOA-VMD-PSR-AWMF multi-stage denoising method proposed in this invention, the maximum Lyapunov exponent of each group of microseismic signals only fluctuated slightly, with the typical variation range controlled within ±7%. This indicates that the denoising process has little impact on the overall nonlinear chaotic characteristics of the signal, and the dynamic essence of the system is preserved. Meanwhile, the attractor complexity, recurrence rate, and other indicators remain basically within the original range, indicating that the phase space structure has not been significantly simplified or destroyed, and the denoised signal can still reflect the intrinsic dynamic behavior of the rock mass fracture process well.

[0068] Correspondingly, indicators such as determinism, average diagonal length, trajectory smoothness, and trajectory regularity all showed a significant improvement trend after denoising: the determinism of most records increased by about 7% to 45%, the average diagonal length and trajectory smoothness increased by more than 40%, and some samples saw increases exceeding 70%. This indicates that denoising effectively suppressed random disturbances in the time series, improved the coherence and predictability of the trajectory in the phase space evolution process, and made the main lobe structure and energy evolution process of microseismic events clearer, which is beneficial for subsequent source mechanism analysis and disaster gestation process identification.

[0069] Table 2 compares the performance of the method of this invention with other denoising algorithms such as EMD-PSR-AWMF, VMD-PSR-AWMF, I-CEEMDAN-PSR-AWMF, EMD-WT-SVD, and VMD-WT-SVD in terms of statistical characteristics such as spectral entropy (SE), kurtosis (K), zero-crossing rate (ZCR), and standard deviation (SD). It can be seen that the spectral entropy obtained by the method of this invention is generally lower than that of other comparative methods, with some records showing a reduction of approximately 9% to 21%, indicating that the signal spectrum energy distribution is more concentrated after denoising, and high-frequency noise energy is effectively suppressed. The zero-crossing rate index decreases significantly, with a typical reduction of approximately 33% to 91%, reflecting a substantial reduction in high-frequency oscillation components and random spike noise, resulting in a smoother and more continuous signal waveform.

[0070] In terms of kurtosis, the method of this invention achieves higher kurtosis values ​​in most samples, showing a significant improvement compared to traditional methods. This indicates that the processed signal has an advantage in preserving or enhancing the characteristics of pulse-type events, and the microseismic main lobe and its transient impact components are not excessively smoothed or weakened. The standard deviation index generally decreases slightly or remains within a reasonable range, indicating that while effectively suppressing noise, the overall signal energy level is not significantly compressed, which is beneficial to the stability of subsequent event identification and energy estimation.

[0071] Figure 2 and Figure 3 The time-domain waveforms of typical microseismic records before and after noise reduction are shown. It can be seen that after processing by the method of this invention, the background noise and high-frequency fine oscillations that are widespread in the original records are significantly reduced, the baseline is more stable, and the waveform contours are smooth and clear. The arrival time of the main lobe, the amplitude of the main peak, and the subsequent decay process of the microseismic events in each record are completely preserved without obvious distortion or over-smoothing, which reflects the event preservation characteristics of this method.

[0072] Figure 4 and Figure 5 The comparison results of the corresponding signals in the frequency domain and time-frequency domain are presented. Compared with the original signal, the energy of the denoised signal within the main frequency band remains basically unchanged or is slightly concentrated, while the broadband high-frequency noise energy outside the main frequency band is significantly attenuated, and the spectral lines are more concentrated and compact. The high-frequency tails of some records are almost completely suppressed, further confirming the quantitative conclusions of the spectral entropy and zero-crossing rate indices. The method of this invention improves the signal-to-noise ratio and highlights the main frequency components without introducing significant spectral distortion or frequency shift.

[0073] Figure 6The comparison results of the phase space probability density distribution of the microseismic signal before and after denoising are shown. It can be observed that after denoising, the distribution range of the phase space point cloud converges, the scattered points caused by noise on the periphery are significantly reduced, and the probability density concentrates from a dispersed state towards the vicinity of the main attractor. However, the overall shape of the attractor and the main distribution area remain basically consistent. This is consistent with the results in Table 1, where the maximum Lyapunov exponent and attractor complexity only change slightly, indicating that the method of this invention effectively preserves the intrinsic phase space structure of the system while reducing random noise.

[0074] Figure 7 This corresponds to a comparison of the phase space dynamic vector fields before and after denoising. It can be seen that the original signal's vector field contains many local regions with directional disorder and abrupt amplitude changes, reflecting strong noise interference. After denoising, the vector field streamlines are smoother and more continuous, the vector directions are more consistent with the main rotation / contraction trend, local extreme vectors are significantly reduced, and the system evolution trajectory exhibits a clearer, more ordered structure. This phenomenon is highly consistent with the significant improvement in determinism, trajectory smoothness, and trajectory regularity, indicating that the method of this invention can significantly improve the interpretability of system evolution without destroying the original dynamic attractor.

[0075] In summary, through comparisons of dynamic characteristic parameters (maximum Lyapunov exponent, recurrence rate, determinism, etc.), statistical characteristic indices (spectral entropy, kurtosis, zero-crossing rate, standard deviation, etc.), and multiple perspectives in the time domain, frequency domain, and phase space, it can be seen that the proposed WOA-VMD-PSR-AWMF multi-stage denoising method outperforms existing typical denoising algorithms in terms of high-frequency noise suppression, event waveform preservation, and dynamic characteristic maintenance. This method can significantly improve the signal-to-noise ratio and energy concentration of microseismic signals while preserving the nonlinear dynamic properties of the original signal as much as possible, providing a more reliable input data foundation for subsequent microseismic event identification, source mechanism inversion, and rockburst disaster early warning.

[0076] Table 1 compares the dynamic and statistical characteristic parameters of the microseismic signal before and after noise reduction.

[0077]

[0078] Table 2 compares the denoising results of microseismic signals under different denoising algorithms.

[0079]

[0080] The technical solutions of the present invention are not limited to the specific embodiments described above. Any technical modifications made in accordance with the technical solutions of the present invention fall within the protection scope of the present invention.

Claims

1. A multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF, characterized in that, include: S1. Obtain the original signal s raw (t), and perform preprocessing; S2. Optimize the key parameters of variational mode decomposition (VMD), and use the optimized key parameters to perform variational mode decomposition (VMD) on the preprocessed original signal to obtain multiple intrinsic mode function (IMF) components. S3. Divide all IMF components into the first component, the second component, and the third component. The first component is the component dominated by the effective signal, the second component is the component that is a mixture of signal and noise, and the third component is the component dominated by the noise signal. S4. The first component, the second component, and the third component are denoised using the first denoising strategy, the second denoising strategy, and the third denoising strategy, respectively. The denoised first component, the second component, and the third component are then reconstructed to obtain the initial denoised signal v(t). S5. Reconstruct the phase space of the initial noise-reduced signal v(t); Within the local neighborhood of the reconstructed phase space, singular value decomposition (SVD) is used to separate the principal components and noise components of the initially denoised signal after reconstruction. The noise component is removed, and the principal components are retained as the IMF components after secondary denoising. k (t); S6, regarding IMF components w k (t) Perform adaptive wavelet-morphological filtering to obtain the filtered signal; S7. Analyze the event window in the filtered signal, smooth the non-event window, and obtain and output the final noise-reduced signal.

2. The multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF according to claim 1, characterized in that, S2 includes: S21. Initialize the population of the Whale Optimization Algorithm (WOA), set the position vector of each whale, and the position vector corresponds to the combination of parameters to be optimized in the Submodal Decomposition (VMD) [K, α], where K is the number of modes and α is the penalty factor; S22. Define the fitness function f(θ), using the reconstruction coefficients and signal-to-noise ratio of the decomposed IMF components as evaluation metrics, expressed as: ,in, , θ=[K,α,τ], For reconstruction error, and These are the weighting coefficients for the reconstruction error term and the signal-to-noise ratio term, respectively. For signal-to-noise ratio loss, Let τ be the vector of parameters to be optimized, and τ be the update step size. S23. Perform WOA (Whale Optimization Algorithm) iterative optimization. By simulating the contraction and spiral predation behavior of humpback whales, search for the optimal θ and finally output the optimal parameter θ* to obtain the optimal parameter combination [K,α] of the fitness function. S24. Using the optimal parameter combination [K,α], perform variational mode decomposition (VMD) on the preprocessed original signal to obtain IMF components distributed from high frequency to low frequency, denoted as: IMF1, IMF2, ..., IMFK.

3. The multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF according to claim 1, characterized in that, In S3, by analyzing the characteristic indicators of each IMF component, all IMF components are divided into the first component, the second component, and the third component.

4. The multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF according to claim 3, is characterized in that, S3 includes: S31. Calculate the characteristic indices of the IMFi components, including the relative energy E. i The absolute value of the Pearson correlation coefficient r with the original signal i Sample entropy S i Gaussian property value p i ; S32. Normalize the feature indicators; S33. Analyze the comprehensive score based on the normalized feature indicators. i , represented as: Score i =Ê i +r̂ i +(1−p̂ i )−0.5·Ŝ i Ê i 、r̂ i 、p̂ i and Ŝ i These are the normalized relative energy, the absolute value of the Pearson correlation coefficient, the sample entropy, and the Gaussianity value, respectively. S34. Based on the overall score i Perform K-means clustering on the IMF components IMFi to obtain the clustering results, and map the clustering results to the first, second, or third component.

5. The multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF according to claim 1, characterized in that, The first noise reduction strategy is to retain all noise or remove noise values ​​with amplitudes less than a first preset threshold. The first component is set to zero; the second noise reduction strategy is to process the second component using an adaptive noise reduction algorithm; the third noise reduction strategy is to set the amplitude to less than a second preset threshold. The third component is set to zero.

6. The multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF according to claim 1, characterized in that, S5 includes: S51. Calculate the self-mutual information of the initial denoised signal v(t) and find the first significant local minimum lag as a candidate delay ζ; if there is no significant minimum within the preset lag range [1, N / 10], then use the lag point where the autocorrelation function R(ζ) decays to R(0) / e as the candidate delay ζ; a significant local minimum refers to... Within the search range, mutual information curve The first local minimum point appears, and the function value at this point decreases by at least 5%–10% compared to the previous local maximum; S52. The embedding dimension m is selected using the pseudo nearest neighbor method (FNN). S53. Using a sliding window, the initial noise-reduced signal v(t) of length N is divided into multiple sub-segments of length W; S54. Construct the trajectory matrix X∈R for all sub-segments. L×m , represented as: Where L = N−(m−1)ζ; S55. Perform singular value decomposition (SVD) on each trajectory matrix X, as follows: , where Σ is a diagonal singular value matrix, and the singular values ​​are sorted in descending order; S56. Determine the truncation quantity r, retain the first r singular components, represented as: ; S57, retain the singular component X denoise The time series is projected back to one dimension by diagonal mean, and the singular components X of all segments are retained. denoise The IMF components are spliced ​​or overlapped on the time axis and then normalized to obtain the IMF components after secondary denoising. k (t).

7. The multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF according to claim 1, characterized in that, S6 includes: S61, for the IMF components after secondary denoising w k (t) Perform wavelet decomposition and obtain wavelet detail coefficients at different scales. The total number of wavelet decomposition levels is L and the decomposition level is ℓ. S62. For each detail coefficient dℓ, calculate the MAD noise estimate σ. ℓ , represented as: σ ℓ =median(|dℓ|) / 0.6745; S63. Calculate the filter threshold T ℓ , represented as: T ℓ =α ℓ ·σ ℓ , This refers to the value of the threshold coefficient; S64. Based on the filtering threshold T ℓ For the IMF components after secondary denoising w k (t) is used for filtering and noise reduction; S65, The IMF components after secondary denoising... k (t) Perform one-dimensional grayscale mathematical morphology operations, design dilation and erosion operators, and base them on the IMF components w k (t) Feature selection: single opening / closing or cascaded opening / closing operations yield the morphological filtering result m. k (t); Opening operation This means that corrosion occurs first, followed by expansion, which is represented as... Closing operation The process involves expansion followed by corrosion, represented as: ; S66, For each IMF component w k (t), the combined IMF component after secondary denoising w k (t), Filtering and noise reduction results y k (t) and morphological branch output m k (t), perform multi-branch fusion to obtain the fused IMF component z k (t), denoted as: z k (t)=β1w k (t)+β2y k (t)+β3m k (t), where β1+β2+β3=1; S67. All IMF components z k Summing (t) yields the processed, denoised micro-vibration signal ŝ(t).

8. The multi-stage event-preserving microseismic signal denoising method based on WOA-VMD, PSR, and AWMF according to claim 1, characterized in that, In S7, based on the envelope energy, short-time energy, and amplitude threshold of the denoised microseismic signal ŝ(t), the start and end times of microseismic events are identified, an event window is formed, and the position and duration of the main lobe of the event are recorded. Smoothing is performed in the non-event window to obtain and output the final noise-reduced signal.

Images