A microseismic signal denoising method based on WPD-EMD-WPD

Abstract

The present invention discloses a microseismic signal denoising method based on WPD-EMD-WPD, comprising the following steps: step S1: obtaining a noisy microseismic signal X={s(1), s(2), ..., s(N)}, and pre-denoising the noisy microseismic signal using a denoising method using a wavelet packet transform to obtain a pre-denoised microseismic signal X′={s(1)′, s(2)′, ..., s(N)′}, wherein N is the signal length; step S2: performing empirical mode decomposition on the pre-denoised microseismic signal X′={s(1)′, s(2)′, ..., s(N)′} to obtain a set of IMF components with frequencies from high to low: imf1, imf2, ..., imf M ; Find the boundary IMF component between noise and effective signal through the correlation coefficient between each IMF component and the pre-denoising microseismic signal X′; Step S3: Perform wavelet packet transform denoising on the IMF component before the boundary, and reconstruct the denoised IMF component and the remaining IMF component to obtain the denoised signal X″={s(1)″,s(2)″,...,s(N)″}. The present invention overcomes the lack of adaptability of wavelet packet transform and solves the problem of low EMD decomposition quality under low signal-to-noise ratio, providing accurate conditions for subsequent analysis of microseismic signals.

Description

Technical Field

[0001] The present invention relates to the intersection of hydropower station safety and signal processing, and in particular to a noise suppression method for microseismic signals in hydropower station caverns with low signal-to-noise ratio and random non-stationary, and is a microseismic signal denoising method based on WPD-EMD-WPD. Background Art

[0002] Microseismic monitoring technology originated from mining engineering. As a three-dimensional "body" monitoring method, microseismic uses sensors to receive rock microfracture signals and then captures the precursor information of rock instability and destruction. Microseismic monitoring technology can capture the internal fracture information of rock mass in real time, directly revealing the evolution process of crack initiation, development, expansion, penetration and macroscopic deformation and destruction, and is an important tool for studying the evolution process of rock damage. It is now mainly used in large-scale projects such as coal mining, large hydropower stations, and underground cavern excavation. Due to the use of microseismic monitoring systems, the number of deaths caused by coal mining in my country has dropped dramatically each year. However, the working conditions in actual projects are complex, and the signals picked up by each channel will be interfered by noise, affecting the first arrival time picking, source location and parameter inversion. Therefore, microseismic signal noise reduction is the basis for carrying out microseismic monitoring and rock stability evaluation, and has very important engineering significance.

[0003] Microseismic signals have the characteristics of non-stationarity, randomness, low signal-to-noise ratio, and complexity. Although Fourier transform can suppress noise to a certain extent, it is mainly applicable to stationary random signals and has poor noise reduction effect on sudden microseismic signals. At present, the noise reduction methods used for non-stationary signals include wavelet transform and empirical mode decomposition (EMD decomposition); however, wavelet packet transform lacks adaptability, and EMD decomposition has low quality under low signal-to-noise ratio, which brings many inconveniences to the noise reduction of microseismic signals. Summary of the invention

[0004] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a microseismic signal denoising method based on WPD-EMD-WPD, which overcomes the lack of adaptability of wavelet packet transform and solves the problem of low EMD decomposition quality under low signal-to-noise ratio, thereby providing accurate conditions for subsequent analysis of microseismic signals.

[0005] The object of the present invention is achieved through the following technical solution: A microseismic signal denoising method based on WPD-EMD-WPD, comprising the following steps:

[0006] Step S1: Obtain a noisy microseismic signal X = {s(1), s(2), ..., s(N)}, and use a wavelet packet transform denoising method to pre-denoise the noisy microseismic signal to obtain a pre-denoised microseismic signal X′ = {s(1)′, s(2)′, ..., s(N)′}, where N is the number of sampling points;

[0007] Step S2: Perform empirical mode decomposition on the pre-denoised microseismic signal X′={s(1)′, s(2)′, ..., s(N)′} to obtain a set of IMF components with high to low frequencies: imf1, imf2, ..., imf M ; Find the boundary IMF component between noise and effective signal through the correlation coefficient between each IMF component and the pre-denoised microseismic signal X′;

[0008] Step S3: Perform wavelet packet transform to denoise the IMF components before the boundary IMF components, reconstruct the denoised IMF components and the remaining IMF components to obtain the denoised signal X″, assuming that the obtained signal X″={s(1)″, s(2)″, ..., s(N)″}.

[0009] Among them, WPD refers to the wavelet packet denoising method. First, the wavelet packet denoising (WPD) method is used for pre-denoising, and then the empirical mode decomposition (EMD) method is used to decompose the IMF. The boundary IMF component between the noise and the effective signal is found through the correlation coefficient between each IMF and the microseismic signal X′ after pre-denoising. Finally, the IMF component before the boundary is subjected to wavelet packet denoising (WPD), and the denoised IMF component and the remaining IMF component are reconstructed to obtain the denoised signal; therefore, the entire process in this application is referred to as WPD-EMD-WPD denoising.

[0010] Furthermore, in the step S1, when the denoising method of wavelet packet transform is used to pre-denoise the noisy microseismic signal, the db9 wavelet and the decomposition layer number of 3 layers that can reflect the characteristics of the original signal are selected;

[0011] In step S3, when performing wavelet packet transform denoising on the IMF component before the boundary IMF component, the db9 wavelet and the number of decomposition layers of 3 are selected.

[0012] Furthermore, the step S2 includes the following sub-steps:

[0013] S201. Perform empirical mode decomposition on the pre-denoised microseismic signal X′={s(1)′,s(2)′,...,s(N)′} to obtain a set of IMF components with high to low frequencies: imf1,imf2,...,imf M ;

[0014] S202. Calculate the correlation coefficient between the i-th IMF component and the pre-denoised microseismic signal X′:

[0015]

[0016] Among them, Cov is the covariance function, Var is the variance function;

[0017] The i-th IMF component is equivalent to being composed of many points. Substituting the value corresponding to each time point into the formula, the correlation coefficient between the IMF component and the signal X′ is calculated;

[0018] S203. When i=1, 2, ..., M, repeat step S202 to obtain each IMF component and the signal X′, wherein the IMF component with the largest correlation coefficient with the signal X′ is the boundary IMF component.

[0019] Furthermore, the step S3 comprises:

[0020] S301. Let imf1, imf2, ..., imf M The jth component is the boundary IMF component; j is an integer between 1 and M; for imf1, imf2..., imf j Each IMF component in is subjected to wavelet packet transform to reduce noise, and imf1′, imf2′..., imf j ′;

[0021] S302. The denoised IMF components imf1′, imf2′, ..., imf j ′ and the residual IMF component imf j+1 ,imf j+2 ,...,imf M Add together to reconstruct the denoised signal and obtain the reconstructed signal X″:

[0022] X″=imf1′+imf2′+...+imf j ′+imf j+1 +imf j+2 ,...,imf M ;

[0023] Assume that the obtained signal is represented by X″={s(1)″,s(2)″,...,s(N)″}.

[0024] Furthermore, the microseismic signal denoising method further includes: calculating a signal-to-deviation ratio (SDR) between the noisy microseismic signal X and the denoised signal X″ as a basis for judging the denoising effect. The signal-to-deviation ratio (SDR) is calculated as follows:

[0025]

[0026] Wherein, s(n) is the value of the nth sampling point in the noisy microseismic signal X, s(n)″ is the value of the nth sampling point in the signal X″, and n=1, 2,…N.

[0027] The beneficial effects of the present invention are as follows: the noise reduction method using EMD decomposition combined with wavelet packet transform is adaptive, and there is no need to repeatedly set the wavelet basis function and the number of wavelet decomposition layers to achieve the optimal noise reduction effect. The uncertainty problems of the selection of basis function type and the number of decomposition layers in wavelet decomposition are solved; the noise reduction effect of the wavelet packet combined with EMD decomposition is better than that of the traditional single noise reduction algorithm, which is mainly manifested in that the high-frequency noise of the signal is effectively suppressed, the signal details are rich, and more useful signals are retained. The proposed signal-to-deviation ratio can better reflect the noise reduction effect, and has a good correspondence with the signal-to-noise ratio and the root mean square error. The smaller the signal-to-deviation ratio, the more the noise is suppressed, so the deviation from the signal without noise is smaller, and the noise reduction effect is more obvious. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 is a flow chart of the method of the present invention;

[0029] Figure 2 The waveform diagram of the original signal and the noisy signal in the embodiment;

[0030] Figure 3 Spectrum diagram of the original signal and the noisy signal in the embodiment;

[0031] Figure 4 Wavelet packet denoising waveform and spectrum diagram in the embodiment;

[0032] Figure 5 The EMD noise reduction waveform and spectrum diagram in the embodiment;

[0033] Figure 6 The noise reduction waveform and spectrum diagram of the present invention in the embodiment;

[0034] Figure 7 Schematic diagram of comparison between the original microseismic signal in the embodiment and the waveform and spectrum after noise reduction of the present invention;

[0035] Figure 8 The waveform diagram and spectrum diagram of EMD noise reduction of microseismic signal in the embodiment;

[0036] Fig. 9 The waveform diagram and spectrum diagram of the microseismic signal wavelet packet denoising in the embodiment;

[0037] Fig.10 It is a time-frequency spectrum diagram of the original microseismic signal in the embodiment and the microseismic signal with noise reduction of the present invention;

[0038] Fig.11 It is a schematic diagram of the principle of wavelet decomposition and reconstruction;

[0039] Fig.12 Schematic diagram of correlation coefficients of IMF components in the embodiment. DETAILED DESCRIPTION

[0040] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings, but the protection scope of the present invention is not limited to the following.

[0041] The present invention first uses wavelet packets to pre-denoise noisy microseismic signals. This step suppresses the influence of noise on EMD decomposition and reduces the phenomenon of modal aliasing. Then, the EMD of the pre-denoised microseismic signal is decomposed into a group of IMFs with high to low frequencies. The boundary IMF components between noise and effective signals are found through the correlation coefficients between each IMF and the pre-denoised microseismic signal. The IMF components before the boundary are subjected to wavelet packet denoising. Finally, the denoised IMF components and the remaining IMF components are reconstructed to obtain the denoised signal. Specifically:

[0042] like Figure 1 As shown, a microseismic signal denoising method based on WPD-EMD-WPD includes the following steps:

[0043] Step S1: Obtain a noisy microseismic signal X = {s(1), s(2), ..., s(N)}, and use a wavelet packet denoising method to pre-denoise the noisy microseismic signal to obtain a pre-denoised microseismic signal X′ = {s(1)′, s(2)′, ..., s(N)′}, where N is the number of sampling points;

[0044] Step S2: Perform empirical mode decomposition on the pre-denoised microseismic signal X′={s(1)′, s(2)′, ..., s(N)′} to obtain a set of IMF components with high to low frequencies: imf1, imf2, ..., imf M ; Find the boundary IMF component between noise and effective signal through the correlation coefficient between each IMF component and the original signal;

[0045] Step S3: Perform wavelet packet denoising on the IMF components before the boundary, reconstruct the denoised IMF components and the remaining IMF components to obtain the denoised signal X″, assuming that the obtained signal X″={s(1)″, s(2)″, ..., s(N)″}.

[0046] Among them, WPD refers to the wavelet packet denoising method. First, the wavelet packet denoising (WPD) method is used for pre-denoising, and then the empirical mode decomposition (EMD) method is used to decompose the IMF. The boundary IMF component between the noise and the effective signal is found through the correlation coefficient between each IMF and the microseismic signal X′ after pre-denoising. Finally, the IMF component before the boundary is subjected to wavelet packet denoising (WPD), and the denoised IMF component and the remaining IMF component are reconstructed to obtain the denoised signal; therefore, the entire process in this application is referred to as WPD-EMD-WPD denoising.

[0047] Furthermore, in step S1, when the wavelet packet denoising method is used to pre-denoise the noisy microseismic signal, the db9 wavelet and 3 decomposition layers that can reflect the characteristics of the original signal are selected.

[0048] like Fig.11 As shown, S0 is equivalent to the original signal (X). S0 can be decomposed by selecting the appropriate wavelet basis (db9 wavelet) and the number of decomposition layers (3 layers). The figure shows a three-layer decomposition.

[0049] S0 is decomposed into L3, H3, L4, H4, L5, H5, L6, H6. (S0 = L3 + H3 + L4 + H4 + L5 + H5 + L6 + H6). The threshold quantization of the decomposition coefficients is to process L3, H3, L4, H4, L5, H5, L6, H6 to obtain L'3, H'3, L'4, H'4, L'5, H'5, L'6, ​​H'6. The wavelet coefficients of the noise are small, and the wavelet coefficients of the noise are smaller than the wavelet coefficients of the signal. By selecting a suitable threshold in advance, the wavelet coefficients greater than the threshold are considered to be generated by the signal and should be retained. The wavelet coefficients less than the threshold are considered to be generated by the noise and set to zero to achieve the purpose of denoising. The reconstruction process is to add L'3, H'3, L'4, H'4, L'5, H'5, L'6, ​​H'6 to obtain the denoised signal S'0.

[0050] In step S3, when performing wavelet packet transform denoising on the IMF component before the boundary IMF component, the db9 wavelet and the number of decomposition layers of 3 are selected.

[0051] Furthermore, the step S2 includes the following sub-steps:

[0052] S201. Perform empirical mode decomposition on the pre-denoised microseismic signal X′={s(1)′,s(2)′,...,s(N)′} to obtain a set of IMF components with high to low frequencies: imf1,imf2,...,imf M ;

[0053] In the embodiment of the present application, the specific process of EMD decomposition includes:

[0054] Step A1: Find all extreme points of the original signal a(t), and use the cubic spline difference function to connect the maximum points and the minimum points respectively to form two upper and lower envelopes;

[0055] Step A2: Calculate the mean of the upper and lower envelopes b(t);

[0056] Step A3: Calculate the difference c1(t)=a(t)-b1(t);

[0057] Step A4: Check whether c1(t) meets the IMF restriction conditions: If it meets the IMF restriction conditions, r1(t) = a(t) - c1(t); If it does not meet the conditions, c1(t) is used as a new signal and steps A1 to A3 are repeated until c1(t) meets the IMF conditions;

[0058] In step A4, the IMF's limiting conditions are:

[0059] Condition 1: imfi(t) must satisfy that the sum of the number of local extreme values ​​of the signal is equal to the number of zero crossings or differs by at most 1;

[0060] Condition 2: imfi(t) must satisfy the condition that the average value of the lower envelope defined by the local extremum of the signal is zero;

[0061] Step A5: If r1 is a monotonic function, the decomposition is completed. If r1 is not a monotonic function, r1 is used as the original function to repeat steps A1 to A5. After the decomposition is completed, the signal is decomposed into n IMF components and 1 trend term r n , expressed as:

[0062]

[0063] S202. Calculate the correlation coefficient between the i-th IMF component and the pre-denoised microseismic signal X′:

[0064]

[0065] Among them, Cov is the covariance function, Var is the variance function;

[0066] The i-th IMF component is equivalent to being composed of many points. Substituting the value corresponding to each time point into the formula, the correlation coefficient between the IMF component and the signal X′ is calculated;

[0067] S203. When i=1, 2, ..., M, repeat step S202 to obtain each IMF component and the signal X′, wherein the IMF component with the largest correlation coefficient with the signal X′ is the boundary IMF component.

[0068] Furthermore, the step S3 comprises:

[0069] S301. Let imf1, imf2, ..., imf M The jth component is the boundary IMF component; j is an integer between 1 and M; for imf1, imf2..., imf j Each IMF component in is subjected to wavelet packet transform to reduce noise, and imf1′, imf2′..., imf j ′;

[0070] In the process of wavelet packet denoising, the decomposition level is 3; the wavelet basis type is db9 wavelet

[0071] S302. The denoised IMF components imf1′, imf2′, ..., imf j ′ and the residual IMF component imf j+1 ,imf j+2 ,...,imf M Add together to reconstruct the denoised signal and obtain the reconstructed signal X″:

[0072] X″=imf1′+imf2′+...+imf j ′+imf j+1 +imf j+2 ,...,imf M ;

[0073] Assume that the obtained signal is represented by X″={s(1)″,s(2)″,...,s(N)″}.

[0074] like Fig.12 As shown, multiple imfs can be obtained through EMD decomposition, and in this embodiment, there are 8. As shown in the figure, one of the decomposed imf1 is imf1. imf1 is equivalent to being composed of many points. The amplitude corresponding to each time point is substituted into the following formula to calculate the correlation coefficient between imf1 and the original signal. The same principle is used to obtain the correlation coefficients of other imfs and the microseismic signal X′ after pre-noise reduction. Finally, the correlation coefficients of each imf and the microseismic signal X′ after pre-noise reduction can be obtained. It can be seen that IMF3 is the decomposition component (the correlation coefficient reaches the maximum value), and then this article continues to perform wavelet packet denoising on the first two imfs, and then reconstructs them with the last six imf components;

[0075] Furthermore, the microseismic signal denoising method further includes: calculating a signal-to-deviation ratio (SDR) between the noisy microseismic signal X and the denoised signal X″ as a basis for judging the denoising effect. The signal-to-deviation ratio (SDR) is calculated as follows:

[0076]

[0077] Wherein, s(n) is the value of the nth sampling point in the noisy microseismic signal X, s(n)″ is the value of the nth sampling point in the signal X″, and n=1, 2,…N.

[0078] In the first embodiment of the present application, a simulated signal is used for denoising verification.

[0079] The Ricker wavelet simulating earthquake is used for simulation analysis, and the expression is:

[0080] f(t)=(1-2π 2 fp 2 t 2 )exp-(π 2 f p 2 t 2 )

[0081] Where f(t) is the amplitude, and fp is the spectrum peak frequency. Take fp = 30 Hz, the sampling frequency is 1 kHz, and the number of sampling points is 1000. Add Gaussian white noise with a signal-to-noise ratio of 4 dB to the Ricker wavelet. The waveform diagram of the original signal and the noisy signal in the embodiment is as follows: Figure 2 As shown, Figure 2 (a) is the original waveform. Figure 2 (b) is the noise waveform diagram;

[0082] Figure 3 is a spectrum diagram of the original signal and the noisy signal in the embodiment, wherein Figure 3 (a) is the original signal spectrum, it can be seen that the main frequency is 30HZ, Figure 3 (b) is the spectrum of the noisy signal, and it can be seen that there is a lot of high-frequency noise. Figure 4 As shown in Figure 2, wavelet packet transform denoising effectively suppresses high-frequency noise, where Figure 4 (a) is the waveform after wavelet packet transform denoising, as shown in Figure 4 As shown in (a), although this method can make the waveform smoother, there are still burrs; Figure 4 (b) is the spectrum after wavelet packet denoising, as shown in Figure 4 As shown in (b), although this method suppresses high-frequency noise to a certain extent, some high-frequency components still exist. Figure 5 (a) is the EMD noise reduction waveform. Figure 5 (b) is the EMD denoised spectrum;

[0083] The waveform and spectrum diagram after denoising of the present invention and the evaluation index values ​​after denoising of the three methods are shown as follows: Figure 6 As shown in Table 1, Figure 6 (a) is a noise reduction waveform diagram of the present invention, Figure 6 (b) is a spectrum diagram of noise reduction of the present invention;

[0084] Table 1 Comparative analysis of noise reduction effects

[0085]

[0086] from Figure 4 , 6 The denoising effect diagram shows that the signal after wavelet packet transform denoising has burrs and residual high-frequency noise, while the signal after denoising of the present invention is closer to the original signal and has less residual high-frequency noise.

[0087] Table 1 shows more detailed results: Compared with the other two methods, the signal after noise reduction of the present invention has the lowest signal-to-bias ratio, the highest signal-to-noise ratio, and the smallest root mean square error, which reflects the practicality and more excellent denoising performance of the present invention, and also shows that the signal-to-bias ratio has a good correspondence with other noise reduction indicators, providing a new idea for evaluating the noise reduction effect.

[0088] In the second embodiment of the present application, microseismic signals are used for denoising verification.

[0089] The actual data comes from the microseismic monitoring of an underground cavern of a large hydropower project. Noisy microseismic signals are randomly extracted from the monitoring data for analysis. The invention of this paper is used to perform noise reduction on the waveform of the microseismic signal. The waveform and spectrum before and after noise reduction are shown in the figure below. Figure 7 (a) to (d) are shown, and the signal-to-deviation ratio (SDR) of the microseismic signal before and after noise reduction using three different methods is calculated. The calculation results are shown in Table 2. Figure 7 (a) It can be seen that the original microseismic signal contains a large amount of random noise, which greatly interferes with subsequent research work such as first arrival time picking, frequency characteristic analysis, and focal mechanism interpretation. Figure 8 The waveform diagram and spectrum diagram of EMD noise reduction of microseismic signal in the embodiment; wherein, Figure 8 (a) is the EMD noise reduction waveform of the microseismic signal, where: Figure 8 (b) is the EMD denoised spectrum of the microseismic signal; Fig. 9 The waveform diagram and spectrum diagram of the microseismic signal wavelet packet denoising in the embodiment are shown in FIG. Fig. 9 (a) is the waveform of microseismic signal denoised by wavelet packet. Fig. 9 (b) is the wavelet packet denoised spectrum of the microseismic signal;

[0090] like Fig.10 The time-frequency spectrum of the microseismic signal before and after noise reduction is shown. Fig.10 (a) is the time spectrum of the original microseismic signal. Fig.10 (b) is the time-frequency spectrum of the microseismic signal after noise reduction by the method of the present invention. Fig. 9 (b) After that, the waveform and spectrum of the microseismic signal are clearer, and the dominant spectrum is more obvious. As can be seen from Table 2, the signal after denoising by the WPD-EMD-WPD method has a lower signal-to-bias ratio, indicating that the signal after denoising is less deviated from the microseismic signal without noise.

[0091] Table 2 Analysis of noise reduction effect of microseismic signals

[0092]

[0093] In summary, the present invention improves the denoising ability of microseismic signals, overcomes the lack of adaptability of wavelet packet transform, and avoids the modal aliasing phenomenon caused by traditional EMD decomposition. Two experiments, noisy Ricker wavelet denoising and actual microseismic signal denoising, show that the method of the present invention can further improve the signal-to-noise ratio of the denoised signal and reduce the root mean square error. It is suitable for denoising of random non-stationary microseismic signals with low signal-to-noise ratio, and the proposed signal-to-noise ratio has a good correspondence with other noise reduction effect indicators, providing a new idea for evaluating the noise reduction effect.

[0094] The above description shows and describes a preferred embodiment of the present invention, but as mentioned above, it should be understood that the present invention is not limited to the form disclosed herein, and should not be regarded as excluding other embodiments, but can be used in various other combinations, modifications and environments, and can be modified within the scope of the invention concept described herein through the above teachings or the technology or knowledge of the relevant field. And the changes and modifications made by those skilled in the art do not depart from the spirit and scope of the present invention, and should be within the scope of protection of the claims attached to the present invention.

Claims

1. A microseismic signal denoising method based on WPD-EMD-WPD, characterized by: The following steps are involved: Step S1: Obtain a noisy microseismic signal X = {s(1), s(2), ..., s(N)}, and use a wavelet packet transform denoising method to pre-denoise the noisy microseismic signal to obtain a pre-denoised microseismic signal X′ = {s(1)′, s(2)′, ..., s(N)′}, where N is the number of sampling points; Step S2: Perform empirical mode decomposition on the pre-denoised microseismic signal X′={s(1)′, s(2)′, ..., s(N)′} to obtain a set of IMF components with high to low frequencies: imf1, imf2, ..., imf M ; Find the boundary IMF component between noise and effective signal through the correlation coefficient between each IMF component and the pre-denoised microseismic signal X′; The step S2 comprises the following sub-steps: S201. Perform empirical mode decomposition on the pre-denoised microseismic signal X′={s(1)′,s(2)′,...,s(N)′} to obtain a set of IMF components with high to low frequencies: imf1,imf2,...,imf M ; S202. Calculate the correlation coefficient between the i-th IMF component and the pre-denoised microseismic signal X′: Among them, Cov is the covariance function, Var is the variance function; The i-th IMF component is equivalent to being composed of many points. Substituting the value corresponding to each time point into the formula, the correlation coefficient between the IMF component and the signal X′ is calculated; S203. When i=1, 2, ..., M, repeat step S202 to obtain each IMF component and signal X', wherein the IMF component with the largest correlation coefficient with signal X' is the boundary IMF component; Step S3: Perform wavelet packet transform to denoise the IMF components before the boundary IMF components, reconstruct the denoised IMF components and the remaining IMF components to obtain the denoised signal X″, assuming that the obtained signal X″={s(1)″, s(2)″, ..., s(N)″}.

2. The microseismic signal denoising method based on WPD-EMD-WPD according to claim 1, characterized in that: In the step S1, when the denoising method of wavelet packet transform is used to pre-denoise the noisy microseismic signal, the db9 wavelet and the decomposition layer number of 3 layers that can reflect the characteristics of the original signal are selected; In the step S3, when performing wavelet packet transform denoising on the IMF component before the boundary IMF component, the db9 wavelet and the number of decomposition layers of 3 are selected.

3. The microseismic signal denoising method based on WPD-EMD-WPD according to claim 1, characterized in that: The step S3 comprises: S301. Let imf1, imf2, ..., imf M The jth component is the boundary IMF component; j is an integer between 1 and M; for imf1, imf2..., imf j Each IMF component in is subjected to wavelet packet transform to reduce noise, and imf1′, imf2′..., imf j ′; S302. The denoised IMF components imf1′, imf2′, ..., imf j ′ and the residual IMF component imf j+1 ,imf j+2 ,...,imf M Add together to reconstruct the denoised signal and obtain the reconstructed signal X″: X″=imf1′+imf2′+...+imf j ′+imf j+1 +imf j+2 ,...,imf M ; Assume that the obtained signal is represented by X″={s(1)″,s(2)″,...,s(N)″}.

4. The microseismic signal denoising method based on WPD-EMD-WPD according to claim 1, characterized in that: The microseismic signal denoising method further includes: calculating a signal-to-deviation ratio (SDR) between the noisy microseismic signal X and the denoised signal X″ as a basis for judging the denoising effect. The SDR is calculated as follows: Wherein, s(n) is the value of the nth sampling point in the noisy microseismic signal X, s(n)″ is the value of the nth sampling point in the signal X″, and n=1, 2,…N.

Images