Micro-seismic wave wavelet denoising method based on VMD decomposition information guidance

Abstract

The invention relates to the technical field of signal processing, in particular to a micro-seismic wave wavelet denoising method based on VMD decomposition information guidance, which comprises the following steps: determining the optimal number of layers of VMD decomposition, and performing VMD decomposition on a micro-seismic waveform to obtain an intrinsic mode function; calculating the corresponding kurtosis based on the intrinsic mode function, judging to obtain a signal component and a noise component, and reconstructing the signal component; performing multi-window long and short window ratio calculation on the signal component, and picking up a P-wave position point of the signal component configured at the long and short window ratio; clustering the P-wave position set by using an improved DBSCAN algorithm to obtain an optimal P-wave position point; dividing the reconstructed signal component into a pure noise segment and an effective feature signal segment; calculating a noise standard deviation and a wavelet threshold based on the pure noise segment; determining the number of decomposition layers of the wavelet basis function; and performing wavelet decomposition and denoising on the reconstructed signal according to the decomposition layer number of the wavelet basis function, the noise standard deviation and the wavelet threshold. According to the method, the noise in the noisy micro-seismic signal is effectively removed.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of signal processing, and particularly relates to a microseismic wavelet denoising method based on VMD decomposition information guidance. BACKGROUND

[0002] As an important monitoring means for early warning of rock burst in the fields of mine, water conservancy and hydropower, traffic tunnel (road) and other engineering, the core of microseismic monitoring technology is to capture the microseismic signals generated by rock mass rupture through high-sensitivity sensors. These signals usually exhibit low-frequency, non-stationary transient waveforms, and the effective frequency band and noise are significantly overlapped. In the actual engineering environment, the microseismic signals are easily affected by multi-source interference such as mechanical vibration, electromagnetic interference and transmission link noise, resulting in a sharp decrease in signal-to-noise ratio (SNR). Especially when the rock mass releases low energy in the early stage of rupture, the effective signal may be completely submerged in the background noise, which directly restricts the accuracy of subsequent key analyses such as source location and energy level calculation.

[0003] In the prior art, the hybrid denoising method based on signal decomposition and threshold optimization has become the mainstream research direction. The Chinese patent "Power harmonic denoising method based on improved sparrow algorithm optimized variational mode decomposition", publication number CN118378027A, improves the effect of sparrow algorithm on variational mode decomposition parameter optimization by introducing Logistic-Tent mapping, balance factor and other strategies. The sparrow algorithm itself has a complex calculation process, and the introduction of new strategies greatly increases the calculation amount. The spectral correlation coefficient is used to screen the components. For the case where the noise amplitude is large (masking signal), the probability of failure of this method is greater. The Chinese patent "Signal denoising method based on improved empirical mode decomposition and wavelet threshold function", publication number CN114970602B, discloses a signal denoising method combining the improved CEEMDAN method with the improved wavelet threshold function. However, the method of obtaining the target wavelet threshold function by using the preset adjustment parameters through the hard threshold function and the soft threshold function has low universality. In the face of signals superimposed with complex noise, the decomposed components are more complex, and the effective signal may be mistakenly deleted. The Chinese invention patent "Ultrasonic echo signal joint denoising method improved EKF combined with wavelet packet cooperation", publication number CN120179990A, solves the problems of poor denoising effect in traditional ultrasonic signal denoising and the difficulty in balancing denoising and feature preservation in the wavelet packet denoising process. This method mainly relies on differential evolution algorithm to optimize the noise covariance matrix, which has high computational complexity; and the problem of high-frequency signal processing is not solved. The Chinese patent "Heart sound signal denoising method with improved wavelet threshold", publication number CN119517062A, discloses a method for improving the wavelet denoising threshold. The standard deviation of each component noise is determined according to the average frequency of the CEEMDAN decomposition components, and then the selection method of the wavelet threshold is improved. In the formula for constructing the wavelet threshold, the parameter d has a great influence on the final denoising result, but it needs to be empirically set, and the universality is insufficient. The Chinese patent "GPR signal denoising method based on variational mode decomposition and singular spectrum analysis", publication number CN113887398B, discloses a method for selecting the K value of variational mode decomposition based on the minimum energy loss ratio. The K value corresponding to the minimum energy loss ratio is determined as the optimal modal number. The process of finding the K value with the minimum energy loss ratio has a large amount of calculation, and there may be local optimum. For complex noise conditions, there may be a situation where the K value is too large.

[0004] It is worth noting that although the existing methods have proposed innovative methods in the field of signal denoising, they generally have insufficient utilization of the characteristics of the noise itself, use uniform parameters to process different types of noisy microseismic signals, and have strong subjectivity due to the pre-set parameters, resulting in low universality of the methods. Some methods innovatively propose parameter optimization methods, but have high computational complexity and complexity, which restricts the universality and efficiency of engineering applications. SUMMARY

[0005] According to the technical problems proposed above, a microseismic wave denoising method based on VMD decomposition information guidance is provided. The present application mainly uses VMD decomposition, multi-window STA / LTA and P-wave position guided noise segment segmentation technology, combined with VMD decomposition information to adaptively determine the wavelet basis function and the number of decomposition layers, so as to use the noise characteristics of a single microseismic signal waveform itself to generate appropriate denoising parameters, effectively remove various types of noise in the noisy microseismic signal, and achieve noise control.

[0006] The technical means adopted by the present application are as follows: The microseismic wave denoising method based on VMD decomposition information guidance comprises the following steps: Microseismic waveforms are collected by microseismic detection sensors during tunnel excavation; The optimal number of VMD decompositions is determined, and the microseismic waveforms are decomposed based on the optimal number of VMD decompositions to obtain intrinsic mode functions; The corresponding kurtosis of the intrinsic mode functions is calculated based on the intrinsic mode functions, and the intrinsic mode functions with kurtosis greater than the kurtosis threshold are determined as signal components, and the intrinsic mode functions with kurtosis less than or equal to the kurtosis threshold are determined as noise components, and the signal components are reconstructed; The signal components are calculated by a multi-window long-short window ratio, the P-wave position points of the signal components in the long-short window ratio configuration are picked up, and the P-wave positions are constructed as a P-wave position set; The P-wave position set is clustered by an improved DBSCAN algorithm to obtain the optimal P-wave position points; The reconstructed signal components are divided into pure noise segments and effective feature signal segments by using the optimal P-wave position points; Based on the pure noise segments, the noise standard deviation and the wavelet threshold are calculated; The decomposition layer number of the wavelet basis function is determined according to the number of VMD decompositions; According to the decomposition layer number of the wavelet basis function, the noise standard deviation and the wavelet threshold, the reconstructed signal is wavelet decomposed and denoised.

[0007] Further, the determination process of the number of VMD decompositions comprises: The microseismic waveforms are Fourier transformed to obtain a frequency domain representation; The peak maximum points of the frequency domain representation are detected, and the effective peak points are obtained by screening according to the preset amplitude threshold and distance threshold, the number of peak points is counted, and the initial mode decomposition layer number is set as the number of peak points; Within a preset search radius, a layer number search range is constructed centered on an initial modal decomposition layer number, for each candidate layer number within the search range, the microseismic time-domain waveform signal is decomposed into K modal components using a variational modal decomposition algorithm, a layer number score function corresponding to the candidate layer number K is calculated according to the K modal components, the layer number score function is used to comprehensively evaluate the reconstruction error under the decomposition and the correlation degree between the modal components, the calculation formula of the layer number score function is:

[0008] wherein, is a weight coefficient, is a layer number score function, is a total sampling point number, is a microseismic waveform, is a candidate variational modal decomposition layer number, is a kth eigenmodal function, is a ith eigenmodal function component, is a jth eigenmodal function component, represents an inner product operation, represents a correlation coefficient between modes; Within the search range, a candidate layer number that makes the layer number score function minimum is selected as an optimal layer number of the VMD decomposition, the optimal layer number of the VMD decomposition is:

[0009] wherein, is an optimal layer number of the VMD decomposition, is an initial modal decomposition layer number, is a search radius.

[0010] Further, the calculation formula of the kurtosis is:

[0011] wherein, is a kurtosis of a kth eigenmodal function component, is a standard deviation of the kth eigenmodal function component, is a mean of k eigenmodal function components, is a total sampling point number, is an amplitude of the kth eigenmodal function component at a discrete time point n, and the calculation formula of the standard deviation of the modal component is: .

[0012] Further, the calculation formula of the reconstructed signal component is:

[0013] wherein, is the kth intrinsic modal function, is the optimal number of VMD decomposition, is the reconstructed signal component, is the selection coefficient of the kth modal component, is the summation variable.

[0014] Further, the multi-window long-short window ratio calculation of the signal component, picking up the P-wave position point of the signal component in the long-short window ratio configuration, comprises: calculating the multi-window long-short window ratio according to the window configuration, the window configuration comprising a first parameter combination, a second parameter combination and a third parameter combination, calculating the energy feature ratio of the signal component at each sampling point under each parameter combination based on the window configuration, and the calculation formula of the energy feature ratio being:

[0015] wherein, is the energy feature ratio value calculated at sampling point n under the mth window configuration of the kth signal component, is the length of the short-time window of the mth window configuration, is the length of the long-time window of the mth window configuration, is the signal amplitude of the kth signal component at sampling point n i, is the signal amplitude of the kth signal component at sampling point n j, is the index of the current time, is the point position of the forward backtracking; calculating a dynamic threshold value based on the energy feature ratio, and the calculation formula of the dynamic threshold value being:

[0016] wherein, is, is a basic threshold value, is an adjustment factor, is the standard deviation of the energy feature ratio of the kth signal component in the mth window, is the arithmetic mean of the energy feature ratio of the kth signal component in the mth window; the sampling point at which the energy feature ratio first exceeds the corresponding dynamic threshold value is recorded as the P-wave position point.

[0017] Further, the P-wave position set is clustered using an improved DBSCAN algorithm to obtain the optimal P-wave position point, comprising: The neighborhood radius, the core point minimum neighbor number and the core point determination condition are defined, the neighborhood radius is equal to the sampling frequency of the microseismic waveform divided by 20, and the calculation formula of the core point minimum neighbor number is:

[0018] Wherein, The core point minimum neighbor number, The total number of signal components, and the formula of the core point determination condition is:

[0019] Wherein, Any other point in the P-wave position set, The P-wave position point currently waiting for determination, The domain radius; Based on the neighborhood radius, the core point minimum neighbor number and the core point determination condition, the core point search, the density connection expansion and the cluster division are performed on the P-wave position set, and the cluster result is output; From the cluster result, the cluster with the most points is selected and recorded as the maximum cluster; The candidate points in the maximum cluster are weighted and collected, and the optimal P-wave position point is calculated, and the calculation formula of the optimal P-wave position point is:

[0020] Wherein, The optimal P-wave position point, The maximum cluster, The weight.

[0021] Further, the calculation formula of the weight is:

[0022] Wherein, The kurtosis of the kth eigenmode function component, , And The first window weight coefficient, the second window weight coefficient and the third window weight coefficient are respectively.

[0023] Further, the noise standard deviation and the wavelet threshold value are calculated based on the pure noise segment, including The noise standard deviation is calculated by the MAD method; The wavelet threshold value is calculated based on the standard deviation and the general threshold value of the information entropy correction:

[0024] Wherein, The wavelet threshold value, The noise standard deviation, is the number of the jth layer coefficient, is the entropy weight of each layer coefficient.

[0025] Further, the decomposition layer number of the wavelet base function is determined according to the number of the VMD decomposition layer, comprising: extract the center frequency of the intrinsic mode function to form a center frequency set; based on the center frequency set, calculate the minimum center frequency and the maximum center frequency; based on the minimum center frequency and the maximum center frequency, calculate the decomposition layer number of the wavelet base function:

[0026] wherein, is the decomposition layer number of the wavelet base function, is the minimum center frequency, is the maximum center frequency, is the maximum allowed layer number, is the sampling frequency.

[0027] Compared with the prior art, the present application has the following advantages: The microseismic wavelet denoising method based on VMD decomposition information guidance provided by the present application realizes the "noise treatment noise" adaptive denoising without presetting fixed parameters and with signal noise characteristics driving parameter generation, effectively improves the universality and engineering application efficiency of the method for different noises.

[0028] Based on the above reasons, the present application can be widely popularized in the field of signal processing. It is suitable for denoising of microseismic signals in the fields of microseismic monitoring and engineering blasting monitoring to improve the accuracy of P-wave first arrival detection and the effect of subsequent microseismic information analysis. BRIEF DESCRIPTION OF DRAWINGS

[0029] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced. Obviously, the drawings in the following description are some embodiments of the present application, and those skilled in the art can obtain other drawings according to these drawings without creative labor.

[0030] Figure 1 is the flow chart of the microseismic wavelet denoising method based on VMD decomposition information guidance of the present application.

[0031] Fig. 2(a) is an example of microseismic wave with noise.

[0032] Figure 2(b) is an example of denoised microseismic waveform of the present application.

[0033] Figure 3 Flow chart for determining the number of decomposition layers of the present application.

[0034] Figure 4 Schematic diagram of VMD decomposition and classification of the present application.

[0035] Figure 5 Schematic diagram of multi-window adaptive threshold STA / LTA detection of the present application.

[0036] Figure 6 Schematic diagram of P-wave position set and clustering result of the present application.

[0037] Figure 7 Example of electrical signal denoising of the present application. DETAILED DESCRIPTION

[0038] In order to enable those skilled in the art to better understand the technical solutions of the present application, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor should fall within the scope of protection of the present application.

[0039] It should be noted that the terms "first", "second", etc. in the specification and claims of the present application and the above-mentioned drawings are used to distinguish similar objects, and do not necessarily describe a specific order or sequence. It should be understood that the data thus used can be interchanged under appropriate circumstances, so that the embodiments of the present application described herein can be implemented in an order other than that illustrated or described herein. In addition, the terms "include" and "have" and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device that includes a series of steps or units does not necessarily have to be limited to those steps or units clearly listed, but can include other steps or units that are not clearly listed or inherent to these processes, methods, products or devices.

[0040] As shown in Figure 1 The present application provides a microseismic wavelet denoising method based on VMD decomposition information guidance, and the steps are specifically as follows: S1. Collecting microseismic waveform by using microseismic detection sensor during tunnel excavation process.

[0041] Due to the complex construction process and environment of the drill and blast tunnel, the microseismic waveform occurring in the surrounding rock collected by the microseismic monitoring sensor is usually superimposed with different types of noise such as electrical signals, power frequency interference noise, mechanical noise and the like due to various working conditions.

[0042] The sampling frequency of the microseismic monitoring sensor is 4000 sampling points per second. The collected microseismic waveform is cut to have a uniform length of 4000, and a single waveform with a duration of 1 second is obtained. The original noisy microseismic waveform and the waveform after denoising are shown in FIG. 2(a) and (b).

[0043] S2. Determine the optimal number of layers of VMD decomposition, and decompose the microseismic waveform based on the optimal number of layers of VMD decomposition to obtain the intrinsic mode function.

[0044] As shown in FIG. 1, the optimal number of layers is determined by the initial number of layers and the optimal number of layers decomposition. Figure 3 As shown in FIG. 1, the optimal number of layers is determined by the initial number of layers and the optimal number of layers decomposition.

[0045] Specifically, the optimal number of layers of VMD decomposition is determined by the following method: First step, Fourier transform the microseismic waveform to obtain the frequency domain representation , and perform frequency domain analysis. Wherein, is the microseismic waveform, is the imaginary unit, is the frequency, is the time. The formula is a known formula for converting a signal from the time domain to the frequency domain.

[0046] Second step, detect the peak maximum point of the frequency domain representation, and screen the effective peak point according to the preset amplitude threshold and distance threshold, count the number of peak points, and set the initial modal decomposition layer number to the number of peak points .

[0047] Specifically, the peak maximum point of the frequency domain waveform is detected, the local maximum point is searched, and the peak value of the local maximum point needs to be greater than 0.1 times the peak maximum value, and the distance between adjacent local maximum points needs to be greater than 0.05 times the signal length, so as to avoid too many peaks being detected. The number of peaks meeting the requirements is counted , and the initial modal decomposition layer number is set to . Through the above requirements, the initial layer number meeting the requirements is 6.

[0048] Third step, within the preset search radius , the layer number search range is constructed with the initial modal decomposition layer number as the center. is the search radius, which is set to 3, and if If less than 3, set the minimum search range as 3, for each candidate layer number within the search range, use the variational mode decomposition algorithm to decompose the microseismic time-domain waveform signal into K modal components According to the K modal components, calculate the layer number score function corresponding to the candidate layer number K, the layer number score function is used to comprehensively evaluate the reconstruction error under the decomposition and the correlation degree between the modal components, and the calculation formula of the layer number score function is:

[0049] Wherein, is a weight coefficient, is 0.1 times the overall variance of the signal, is the layer number score function, is the total number of sampling points, is the candidate variational mode decomposition layer number, is the kth intrinsic modal function, is the ith intrinsic modal function component, is the jth intrinsic modal function component, represents the inner product operation, represents the correlation coefficient between the modes.

[0050] In this embodiment, after setting the initial layer number as 6, the preset layer number range is set as 3, 9, all available integers 3, 4, 5, 6, 7, 8, 9 between 3 and 9 are traversed, and then decomposition is performed. The time-domain example waveform is decomposed from 3 to 9, and each time the modal component corresponding to the modal number is decomposed.

[0051] Fourth step, in the search range, select the candidate layer number that makes the layer number score function minimum as the optimal layer number of VMD decomposition, and the optimal layer number of VMD decomposition is:

[0052] Wherein, is the optimal layer number of VMD decomposition, is the initial modal decomposition layer number, is the search radius.

[0053] In this embodiment, the calculation result is that when the decomposition layer number is 9, the score function value is minimum, and the decomposition effect is best.

[0054] S3. Calculate the corresponding kurtosis based on the intrinsic modal function, and determine the intrinsic modal function with kurtosis greater than 3 as a signal component, and determine the intrinsic modal function with kurtosis less than or equal to 3 as a noise component, and reconstruct the signal component.

[0055] Based on the kurtosis index, the 9 modal components decomposed from the example waveform are classified, and the decomposition result is as followsFigure 4 The kurtosis is shown in FIG. 6 (where Figure 4 The middle gray part is the noise component). The calculation formula of the kurtosis is:

[0056] wherein, is the kurtosis of the kth eigenmode function component, is the standard deviation of the kth eigenmode function component, is the mean of the k eigenmode function components, is the total number of sampling points, is the amplitude of the kth eigenmode function component at the discrete time point n, and the calculation formula of the standard deviation of the mode component is: .

[0057] As a preferred mode of the present embodiment, the kurtosis threshold is set to 3.

[0058] The calculation formula of the reconstructed signal component is:

[0059] wherein, is the kth eigenmode function, is the optimal number of layers of VMD decomposition, is the reconstructed signal component, is the selection coefficient of the kth mode component, is the summation variable.

[0060] The calculation formula of the selection coefficient of the kth mode component is:

[0061] wherein, is the kurtosis threshold.

[0062] S4. Perform multi-window STA / LTA calculation on the signal component, pick up the P-wave position point of the signal component configured by the long-short time window ratio, and construct the P-wave position as a P-wave position set.

[0063] Specifically, S4 includes: First step, calculate the signal component Calculate the multi-window STA / LTA (unit: sampling points) according to the window configuration, and the window configuration includes a first parameter combination, a second parameter combination and a third parameter combination. Based on the window configuration, calculate the energy feature ratio of the signal component at each sampling point under each parameter combination. The calculation formula of the energy feature ratio is:

[0064] wherein, Energy feature ratio calculated at sampling point n for the mth window configuration of the kth signal component, Short time window length configured for the mth window configuration, Long time window length configured for the mth window configuration, Signal amplitude at sampling point n-i for the kth signal component, Signal amplitude at sampling point n-j for the kth signal component, Signal amplitude at sampling point n-i for the kth signal component, Signal amplitude at sampling point n-j for the kth signal component, Index of current time, Point of backtracking. n-i and n-j are used to calculate the energy ratio of the current time using historical points. m is the number of windows, which is 1, 2, 3.

[0065] The short windows of the three sets of windows are (20, 50, 50) in turn; the long windows of the three sets of windows are (200, 200, 500) in turn; that is, Window1 is (20, 200), which is used to capture the initial peak mid-frequency component: Window2 is (50, 200), which is used to balance the noise and suppress the low-frequency component: Window3 is (50, 500) which is used to capture the long-period waveform. Second step, based on the energy feature ratio, calculate the dynamic threshold, the calculation formula of the dynamic threshold is:

[0066] Wherein, is, is the basic threshold, is the adjustment factor, is the standard deviation of the energy feature ratio of the kth signal component in the mth window, is the arithmetic mean of the energy feature ratio of the kth signal component in the mth window.

[0067] As a preferred mode of the present embodiment, the basic threshold is taken as 2.2, and the adjustment factor is taken as 0.7.

[0068] Third step, record the sampling point where the energy feature ratio first exceeds the corresponding dynamic threshold as the P wave position point. Through kurtosis classification, a total of 5 effective signal components are obtained, and each component picks up 3 effective P wave positions through multi-window STA / LTA detection, and the picking result is as shown in Figure 5 .

[0069] Fourth step, record all triggered Wave position points as a set, denoted as ; wherein represents the P wave position point detected by the signal component in the mth window configuration.

[0070] S5. The P-wave position set is clustered using an improved DBSCAN algorithm to obtain the optimal P-wave position point.

[0071] The first step is to define the neighborhood radius, the minimum number of neighbors of the core point, and the core point determination condition. The neighborhood radius is equal to the sampling frequency of the microseismic waveform divided by 20. The calculation formula of the minimum number of neighbors of the core point is:

[0072] wherein, is the minimum number of neighbors of the core point, is the total number of signal components, and the formula of the core point determination condition is:

[0073] wherein, is any other point (sample index point) in the P-wave position set, used to determine whether it is within the domain of, is the P-wave position point to be determined (also a sample index point), representing a candidate P-wave arrival time point in the microseismic waveform, is the domain radius, used to define the maximum distance threshold between points. Here, the value is the sampling frequency of the microseismic waveform divided by 20, and the result is in units of sample points, representing the allowed deviation range in time index.

[0074] The second step is to perform core point search, density connection expansion, and cluster division on the P-wave position set based on the neighborhood radius, the minimum number of neighbors of the core point, and the core point determination condition, and output the clustering result.

[0075] The third step is to select the cluster with the most points from the clustering result, denoted as the maximum cluster.

[0076] The fourth step is to perform a weighted set on the candidate points in the maximum cluster, and calculate the optimal P-wave position point. The calculation formula of the optimal P-wave position point is:

[0077] wherein, is the optimal P-wave position point, the P-wave position candidate points are divided into several clusters (clusters), is the maximum cluster, i.e., the cluster containing the most candidate points, is the weight. The calculation formula of the weight is:

[0078] wherein, is the k-th eigenmode function component kurtosis, , and The first window weight coefficient, the second window weight coefficient, and the third window weight coefficient are respectively.

[0079] The clustering result is shown as Figure 6 .

[0080] S6. Using the optimal P-wave position point, the reconstructed signal component is divided into a pure noise segment and an effective feature signal segment.

[0081] Specifically, according to the position of the obtained optimal P-wave, the reconstructed signal is divided into two segments, i.e., a pure noise segment before the P-wave position point and an effective feature signal segment after the P-wave position point. The effective feature signal segment contains the main part of the microseismic signal and the calculation information.

[0082] The signal length before the P-wave is regarded as pure noise, denoted as ; and the signal length after the P-wave is regarded as an effective feature signal segment, denoted as .

[0083] S7. Based on the pure noise segment, the noise standard deviation and the wavelet threshold are calculated.

[0084] Specifically, the noise standard deviation is calculated by the MAD method.

[0085] The wavelet threshold is calculated based on the standard deviation and the general threshold value corrected by the information entropy:

[0086] wherein, is the wavelet threshold, is the noise standard deviation, is the number of the jth layer coefficient, is the entropy weight of each layer coefficient. The wavelet threshold is selected as , and the basis function is selected as db4.

[0087] S8. The decomposition layer number of the wavelet basis function is determined according to the number of VMD decomposition layers.

[0088] First step, extract the center frequency of the intrinsic mode function to form a center frequency set.

[0089] Let the center frequency set of the modal components obtained by VMD decomposition be: .

[0090] Second step, based on the center frequency set, calculate the minimum center frequency and the maximum center frequency . The minimum center frequency of the waveform is 120.00 Hz, and the maximum center frequency is 594 Hz.

[0091] According to the similarity between the wavelet base and the microseismic signal, the db wavelet is selected as the wavelet base. The wavelet base is selected according to the frequency band range requirement, and the screening criteria are as follows:

[0092] Third step, based on the minimum center frequency and the maximum center frequency, the decomposition level of the wavelet base function is calculated:

[0093] Wherein, is the decomposition level of the wavelet base function, is the minimum center frequency, is the maximum center frequency, is the maximum allowed number of layers, which is 10 layers here. is the sampling frequency. The decomposition level only needs to make the maximum and minimum center frequencies covered after decomposition.

[0094] Fourth step, parameter regularization, . Through calculation, it is found that when the decomposition level is 6, the minimum frequency and the maximum frequency requirements can be completely covered.

[0095] S9. According to the decomposition level of the wavelet base function, the noise standard deviation and the wavelet threshold, the reconstructed signal is wavelet decomposed and denoised.

[0096] Using the adaptively selected wavelet base function and the decomposition level , the discrete wavelet transform (DWT) is performed on the reconstructed signal to obtain the wavelet coefficients of each layer: approximation coefficient (low frequency part of the layer) detail coefficient (high frequency part of the layer) Wavelet decomposition formula:

[0097] Wherein, and are the translation and scaling forms of the scale function and the wavelet function respectively, is the approximation coefficient of the jth layer, is the detail coefficient of the jth layer, is the decomposition level index, is the translation position index, is the input discrete signal, is the time sampling point index.

[0098] Using the threshold processed coefficient and reserved approximation coefficients , inverse discrete wavelet transform (IDWT) to obtain the de-noised signal :

[0099] Figure 7 is a schematic diagram of the de-noised result of the present embodiment.

[0100] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, but not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A wavelet denoising method for microseismic waves guided by VMD decomposition information, characterized in that, Includes the following steps: Microseismic waveforms are collected using microseismic sensors during tunnel excavation. Determine the optimal number of layers for VMD decomposition, and perform VMD decomposition on the micro-vibration waveform based on the optimal number of layers to obtain the intrinsic mode functions; The kurtosis of the intrinsic mode functions is calculated based on the intrinsic mode functions. Intrinsic mode functions with kurtosis greater than the kurtosis threshold are determined to be signal components, and intrinsic mode functions with kurtosis less than or equal to the kurtosis threshold are determined to be noise components. The signal components are then reconstructed. The signal components are subjected to multi-window long-short window ratio calculation, and the P-wave position points of the signal components in the long-short window ratio configuration are picked out. The P-wave positions are then constructed into a P-wave position set. The P-wave location set is clustered using the improved DBSCAN algorithm to obtain the optimal P-wave location points. Using the optimal P-wave location, the reconstructed signal components are divided into pure noise segments and effective feature signal segments; Based on the pure noise segment, calculate the noise standard deviation and wavelet threshold; The number of decomposition layers of the wavelet basis functions is determined based on the number of layers of the VMD decomposition. Based on the decomposition level of the wavelet basis function, the noise standard deviation, and the wavelet threshold, wavelet decomposition and denoising are performed on the reconstructed signal.

2. The microseismic wavelet denoising method based on VMD decomposition information guided by claim 1, characterized in that, The process of determining the number of layers in the VMD decomposition includes: Perform a Fourier transform on the micro-vibration waveform to obtain its frequency domain representation; The peak maxima points represented in the frequency domain are detected, and effective peak points are selected based on preset amplitude thresholds and distance thresholds. The number of peak points is counted, and the initial modal decomposition layer number is set as the number of peak points. Within a preset search radius, a search range is constructed centered on the initial modal decomposition layer number. For each candidate layer within the search range, the microseismic time-domain waveform signal is decomposed into K modal components using a variational modal decomposition algorithm. A layer score function corresponding to the candidate layer K is calculated based on these K modal components. This layer score function is used to comprehensively evaluate the reconstruction error under this decomposition and the correlation between the modal components. The formula for calculating the layer score function is as follows: in, These are the weighting coefficients. For the layer-based scoring function, This represents the total number of sampling points. It is a micro-vibration waveform. represents the candidate variational mode decomposition layer number. Let k be the intrinsic mode function. For the i-th intrinsic mode function component, For the j-th intrinsic mode function component, This represents the inner product operation. This represents the correlation coefficient between modes; Within the search range, candidate layers that minimize the layer score function are selected as the optimal layer number for the VMD decomposition. The optimal layer number for the VMD decomposition is: in, The optimal number of layers for VMD decomposition. This represents the initial number of mode decomposition layers. The search radius is [value].

3. The microseismic wavelet denoising method based on VMD decomposition information guided by claim 1, characterized in that, The formula for calculating the kurtosis is: in, Let be the kurtosis of the k-th intrinsic mode function component. Let be the standard deviation of the k-th intrinsic mode function component. Let the mean of the k intrinsic mode function components be denoted as . This represents the total number of sampling points. Let n be the amplitude of the k-th intrinsic mode function component at discrete time point n. The formula for calculating the standard deviation of the mode component is: 。 4. The microseismic wavelet denoising method based on VMD decomposition information guided by claim 1, characterized in that, The formula for calculating the reconstructed signal components is as follows: in, Let k be the intrinsic mode function. The optimal number of layers for VMD decomposition. For the reconstructed signal components, The selection coefficient for the k-th modal component is... For the summation variable.

5. The microseismic wavelet denoising method based on VMD decomposition information guided by claim 1, characterized in that, The step of calculating the long-to-short window ratio of the signal components and picking the P-wave position points of the signal components at the long-to-short window ratio configuration includes: The ratio of long to short windows in a multi-window configuration is calculated based on the window configuration, which includes a first parameter combination, a second parameter combination, and a third parameter combination. For each parameter combination based on the window configuration, the energy characteristic ratio of the signal component at each sampling point is calculated. The formula for calculating the energy characteristic ratio is: in, Configure the energy characteristic ratio calculated at the nth downsampling point for the mth group window of the kth signal component. The short time window length configured for the m-th group of windows. The long window length configured for the m-th group of windows. For the k-th signal component at sampling point n The signal amplitude at point i, For the k-th signal component at sampling point n The signal amplitude at point j, For the current time index, The points for backward tracing; Based on the energy characteristic ratio, a dynamic threshold is calculated, and the formula for calculating the dynamic threshold is as follows: in, for, Based on the threshold, As a regulating factor, Let be the standard deviation of the energy eigenvalue ratio of the k-th signal component in the m-th window. The arithmetic mean of the energy eigenvalues ​​of the k-th signal component in the m-th window; The sampling point where the energy characteristic ratio first exceeds the corresponding dynamic threshold is recorded as the P-wave location point.

6. The microseismic wavelet denoising method based on VMD decomposition information guided by claim 1, characterized in that, The step of clustering the P-wave location set using the improved DBSCAN algorithm to obtain the optimal P-wave location points includes: Define the neighborhood radius, the minimum number of neighbors for the core point, and the core point determination criteria. The neighborhood radius is equal to the sampling frequency of the microseismic waveform divided by 20. The formula for calculating the minimum number of neighbors for the core point is as follows: in, The minimum number of neighbors for the core point. Given the total number of signal components, the formula for determining the core point is: in, For any other point in the set of P-wave locations, This is the current P-wave location point awaiting determination. The radius of the domain; Based on the neighborhood radius, the minimum number of neighbors of the core point, and the core point determination criteria, the core point search, density connectivity expansion, and clustering are performed on the P-wave location set, and the clustering results are output. From the clustering results, the cluster with the most points is selected and denoted as the largest cluster; The candidate points in the maximum cluster are weighted and the optimal P-wave location is calculated. The formula for calculating the optimal P-wave location is as follows: in, This is the optimal P-wave location point. For the largest cluster, For weights.

7. The microseismic wavelet denoising method based on VMD decomposition information guided by claim 6, characterized in that, The formula for calculating the weight is: in, Let be the kurtosis of the k-th intrinsic mode function component. , and These are the weighting coefficients for the first window, the second window, and the third window, respectively.

8. The microseismic wavelet denoising method based on VMD decomposition information guided by claim 1, characterized in that, The calculation of noise standard deviation and wavelet threshold based on pure noise segments includes... Calculate the noise standard deviation using the MAD method; Calculate the wavelet threshold based on a universal threshold corrected for standard deviation and information entropy: in, For wavelet threshold, The standard deviation of noise. Let j be the number of coefficients in the j-th layer. The entropy values ​​of each layer are weighted.

9. The microseismic wavelet denoising method based on VMD decomposition information guided by claim 1, characterized in that, The determination of the decomposition level of the wavelet basis function based on the number of VMD decomposition levels includes: Extract the center frequencies of the intrinsic mode functions to form a set of center frequencies; Based on the set of center frequencies, calculate the minimum center frequency and the maximum center frequency; Based on the minimum and maximum center frequencies, calculate the number of decomposition levels of the wavelet basis functions: in, The decomposition level of the wavelet basis functions is denoted as . For the minimum center frequency, For the maximum center frequency, The maximum allowed number of floors, The sampling frequency.

Images