Seismic signal denoising method based on PSO-VMD combined improved wavelet transform

Abstract

The invention discloses a seismic signal denoising method based on PSO-VMD combined improved wavelet transform. The seismic signal denoising method comprises the following steps: S1, obtaining noisy and non-noisy seismic signals; s2, searching an optimal decomposition layer number K and a penalty factor alpha of VMD decomposition by using a particle swarm algorithm; s3, performing VMD decomposition on the seismic signal to obtain K modal functions; s4, kurtosis values of the K modal functions are calculated, and the modal functions containing noise and not containing noise are screened out; s5, retaining the original modal component of the modal function which is not smaller than the set kurtosis threshold value, and performing improved wavelet transform denoising on the modal function which is smaller than the set kurtosis threshold value to obtain an approximate function of the seismic signal and the wavelet function; and S6, reconstructing the modal function after noise reduction and the reserved modal function to obtain a seismic signal after noise reduction. According to the method, more effective signals can be reserved while similar calculation efficiency is kept, and the method has good stability.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic signals, and in particular relates to a seismic signal denoising method using a PSO-VMD combined with improved wavelet transform. Background Art

[0002] During the acquisition process, seismic signals are often interfered with by various noises, which come not only from the instrument itself but also from environmental factors or complex factors in the propagation of seismic waves.

[0003] The presence of noise greatly reduces the quality of seismic signals, affecting subsequent data processing and seismic wave analysis. Improving the signal-to-noise ratio (SNR) to effectively remove noise and retain valid signals has always been one of the core tasks in seismic signal processing.

[0004] There are several main types of seismic signal denoising techniques, including time domain methods, frequency domain methods, wavelet transform methods, matrix decomposition methods, etc. Each of these methods has its own advantages and disadvantages.

[0005] Variational mode decomposition (VMD) automatically extracts modal functions from signals through an optimization algorithm and can effectively avoid the modal aliasing problem in the EMD method, with good time-frequency resolution. However, the denoising effect of VMD is affected by the number of decomposition layers (K) and the penalty factor (α). The selection of these parameters is crucial to the denoising effect. Particle swarm optimization (PSO) is an efficient global optimization algorithm that can quickly find the optimal solution in a complex parameter space. Wavelet transform can analyze signals and noise at multiple scales and is suitable for local noise removal in seismic data. However, it requires the selection of an appropriate wavelet basis and threshold. Excessive denoising may cause signal loss. Therefore, studying VMD, PSO, and wavelet transform denoising algorithms is of great significance for improving seismic data resolution and subsequent data interpretation.

[0006] At present, Chinese invention patent CN119537772A discloses a seismic data reconstruction method based on VMD filtering and U-net image segmentation. This method performs spectral analysis on post-stack seismic data to clarify the amplitude and frequency characteristics of bad tracks, thereby effectively removing bad tracks. Variational mode decomposition (VMD) is then used to filter the seismic data after bad track removal, significantly improving the signal-to-noise ratio of the seismic data. A U-net network structure model is used for data reconstruction, and the model is converged through multiple training and iterations to achieve accurate reconstruction of the seismic data. Numerical experiments have shown that this method has good results in processing bad track removal and data reconstruction, and can effectively improve the quality and accuracy of seismic data.

[0007] For example, Chinese invention patent CN115436469A discloses a liner ultrasonic echo signal denoising method based on PSO-VMD. This method first constructs an ultrasonic echo signal with double echoes through simulation signals, and uniformly adds Gaussian white noise to simulate the actual noise environment; uses a variational mode decomposition algorithm to analyze the influence of different decomposition scales and penalty factors on the signal decomposition effect, and further uses a particle swarm optimization algorithm (PSO) to automatically select the optimal decomposition scale and penalty factor, thereby optimizing the VMD algorithm; uses the optimized VMD algorithm to decompose the ultrasonic echo signal collected on site, screens the effective IMF components, and reconstructs the signal through correlation coefficient calculation and threshold selection, ultimately achieving ultrasonic echo signal denoising;

[0008] For example, Chinese invention patent CN116975527A discloses a fault diagnosis method based on soft thresholding and wavelet transform denoising. This method collects vibration signals and divides the data into training, test, and validation sets. The signals are preprocessed by wavelet transform to obtain high-frequency and low-frequency components, which are expressed through a two-dimensional time-frequency graph. A nonlinear layer of soft thresholding and attention mechanism is constructed to adaptively obtain the threshold of each channel feature. Fault diagnosis is performed in combination with a ResNet network model, ultimately determining the fault type and location. Experiments have shown that the fault diagnosis model based on this method has significant advantages in removing noise and improving signal quality, and can efficiently and accurately classify different types of faults.

[0009] Variational mode decomposition (VMD) and wavelet transform are commonly used methods for seismic signal denoising. However, existing VMD and wavelet transform algorithms still have certain limitations. VMD is easily affected by the number of decomposition levels (K) and the penalty factor (α) during the decomposition process, while wavelet transform requires the selection of an appropriate wavelet basis and threshold. The parameter selection has a significant impact on the denoising effect. Therefore, it is of great significance to study an efficient and high-precision seismic signal denoising method. Summary of the Invention

[0010] The purpose of the present invention is to overcome the defects of the existing technology and provide a seismic signal denoising method based on a PSO-VMD combined improved wavelet transform. The particle swarm optimization (PSO) algorithm is used to optimize the decomposition level (K) and the penalty value (α). The seismic signal is decomposed by VMD to obtain multiple modal functions. On this basis, the threshold function in the wavelet transform is improved and the wavelet coefficients are smoothly adjusted. While maintaining similar computational efficiency, this method can retain more effective signals and has good stability. It is particularly suitable for denoising seismic signals containing complex noise and has significant advantages.

[0011] The object of the present invention is achieved through the following technical solutions:

[0012] A seismic signal denoising method based on PSO-VMD combined with improved wavelet transform includes the following steps:

[0013] Step S1, obtaining seismic signals with and without noise;

[0014] Step S2, using particle swarm optimization to search for the optimal decomposition level K and penalty factor α of VMD decomposition;

[0015] Step S3, performing VMD decomposition on the seismic signal to obtain K mode functions;

[0016] Step S4: Calculate the kurtosis values ​​of the K modal functions, compare the kurtosis values ​​of the K modal functions with a set kurtosis threshold, and screen out the modal functions containing noise and those not containing noise;

[0017] Step S5: retaining the original modal components of the modal functions that are not less than the set kurtosis value threshold, performing improved wavelet transform denoising on the modal functions that are less than the set kurtosis value threshold, and obtaining an approximate function of the seismic signal and the wavelet function;

[0018] Step S6: reconstruct the denoised modal function and the retained modal function to obtain a denoised seismic signal.

[0019] In one embodiment, in step S1, the original noise-free seismic signal x(t) and the noisy seismic signal f(t)=x(t)+c(t) are read in, where c(t) is the noise and t is the time parameter, to obtain noisy and noise-free seismic signals.

[0020] In one embodiment, in step S2, the calculation formula is as follows:

[0021] v i,j (t+1)=ωv i,j (t)+c1rand1(pBest i,j -x ij (t))+c2rand2(gBest i,j -x i,j (t));

[0022] x i,j (t+1)=x i,j (t)+v i,j (t+1);

[0023] Among them, v ij (t) is the velocity of particle i in dimension j at time t, ω is the inertia weight, which controls the influence of the particle's last velocity, and c1 is the particle's optimal position toward the individual pBest i,j Adjustment strength, x ij(t) is the position of the particle, that is, the decomposition level K and penalty value α in VMD, and c2 is the particle's adjustment gBest towards the global optimal position. i,j strength, rand1 and rand2 are random numbers, pBest i,j is the individual optimal position of particle i in dimension j, gBest i,j is the global optimal position of particle i in dimension j, that is, the optimal solution.

[0024] In one embodiment, step S3 includes:

[0025] Step S301: Perform VMD decomposition on the seismic signal f(t) to obtain K IMF functions u with different center frequency and limited bandwidth. k (t), the constraint equations of VMD decomposition are as follows:

[0026]

[0027] st∑ k u k (t) = f(t);

[0028] Among them, u k (t) represents the K IMF components obtained by decomposition, ω k Represents the frequency center of each component, * is the convolution symbol, It means to find the derivative of the function with respect to time t, δ(t) is the unit pulse function, represents the square of the gradient 2 norm;

[0029] Step S302: introducing a penalty factor α and a Lagrangian multiplication operator λ to obtain an augmented Lagrangian function L;

[0030]

[0031] Where L is the Lagrangian function, α is the penalty factor, Indicates the derivative of the function with respect to time t, ω k represents the frequency center of each component, δ(t) is the unit pulse function, Indicates the derivative of the function with respect to time t, * is the convolution symbol, and λ(t) is the Lagrange operator. represents the square of the gradient 2 norm;

[0032] Step S303: Solve the Lagrangian function using the alternating direction multiplier algorithm and iterate to obtain u k 、ω k ;

[0033]

[0034] Where ^ represents Fourier transform, λ is the Lagrange multiplication operator, and u k is the modal function, ω k is the center frequency, n represents the number of iterations;

[0035] Step S304: Determine whether to stop iterating and updating the modal function u k , the judgment formula is as follows:

[0036]

[0037] Among them, n represents the number of iterations, u k is the modal function, θ is the given discrimination accuracy, and ω is the center frequency.

[0038] In one embodiment, in step S4, calculating the kurtosis values ​​of the K modal functions includes:

[0039]

[0040] Among them, m i is the i-th sample value of the seismic signal, N is the number of sample points of the seismic signal, μ is the mean of the signal, and kurtosis is the kurtosis value.

[0041] In one embodiment, in step S4, filtering out the modal functions containing noise and those not containing noise includes:

[0042] If the modal function u k (t) is less than the set kurtosis value threshold, it is judged to be a modal function containing noise. If the modal function u k If (t) is not less than the set kurtosis value threshold, it is judged to be a modal function that does not contain noise.

[0043] In one embodiment, in step S5, the modal functions with a kurtosis value less than a preset threshold are subjected to improved wavelet transform denoising to obtain an approximate function of the seismic signal and the wavelet function, including:

[0044]

[0045] Among them, the scale variable a controls the expansion and contraction of the wavelet function, and the translation variable b controls the translation of the wavelet function. is the wavelet basis function, u k (t) is the mode function, W u (a, b) is the approximate function obtained by doing the inner product of the seismic signal and the wavelet basis.

[0046] In one embodiment, step S5 further includes performing multi-scale analysis on the selected modal functions and calculating the scale coefficient c of the modal function using a suitable wavelet basis function. j,i and wavelet coefficient dj,i ,

[0047]

[0048] Among them, c j,i is the low-frequency approximate component at scale j and position i, d j,i is the high-frequency detail component at scale j and position i, φ j,i (t) is the scaling function, ψ j,i (t) is the wavelet function, and u(t) is the modal function.

[0049] In one embodiment, step S5 further includes:

[0050] Set the number of wavelet decomposition layers and calculate the threshold for each layer, including:

[0051]

[0052] σ=media(|d j,N |) / 0.6745;

[0053] Among them, σ is the noise variance of the signal to be processed, N is the number of wavelet coefficients with higher frequencies under the scale change, j represents the decomposition scale, media(|d j,N |) is the median of the absolute values ​​of all wavelet coefficients at scale j, and ρ is the threshold value;

[0054] Quantize the detail coefficients that are smaller than the threshold, including:

[0055]

[0056] Among them, d * j,i is the quantized detail coefficient, d j,i is the original detail coefficient, ρ n is the threshold value of each layer, sign is the sign function;

[0057] For the original detail coefficient d j,i Perform quantization processing to obtain the quantized detail coefficient d * j,i , and then reconstruct to get the denoised modal function u k (t) * ,include:

[0058]

[0059] Among them, u k (t) * is the mode function after denoising, c iis the approximation coefficient, φ(tk) is the scaling function, ψ(2 j tk) is an approximate function, d * j,i is the quantized detail coefficient, c i is the approximate coefficient.

[0060] In one embodiment, step S6 includes:

[0061] The denoised mode function and the unprocessed mode function are reconstructed to obtain the denoised seismic signal:

[0062] ∑ k u k (t)+u k (t) * =f(t);

[0063] In the formula, u k (t) is the retained mode function, u k (t) * is the mode function after denoising.

[0064] The beneficial effects of the present invention are:

[0065] This method targets random noise in seismic signals. The invention uses a particle swarm optimization algorithm to optimize the decomposition layer number (K) and penalty value (α), decomposes the seismic signal through VMD, and obtains multiple modal functions. On this basis, the threshold function in the wavelet transform is improved and the wavelet coefficients are smoothly adjusted. While maintaining similar computational efficiency, this method can retain more effective signals and remove more noise. It has high stability and robustness, and has good stability. It is particularly suitable for denoising seismic signals containing complex noise and has significant advantages. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] The present invention will be described in more detail below based on embodiments and with reference to the accompanying drawings, wherein:

[0067] Figure 1 Shows a flow chart of the seismic signal denoising method of the present invention;

[0068] Figure 2 Schematic diagram showing the original seismic signal and the noisy seismic signal;

[0069] Figure 3 Schematic diagram showing the decomposition of VMD into K modal functions (K=5);

[0070] Figure 4 A schematic diagram of the denoising results of the VMD decomposition method of the present invention is shown;

[0071] Figure 5A schematic diagram of the wavelet basis (db3) of the present invention is shown;

[0072] Figure 6 shows a schematic diagram of the denoising result of the present invention;

[0073] In the drawings, like reference numerals are used for like parts, but the drawings are not necessarily true to scale. DETAILED DESCRIPTION

[0074] The present invention will be further described below with reference to the accompanying drawings.

[0075] The present invention provides a seismic signal denoising method based on PSO-VMD combined with improved wavelet transform, such as Figure 1 As shown, the following steps are included:

[0076] Step S1: Obtain seismic signals with and without noise to read in seismic signals f(t), such as Figure 2 As shown;

[0077] f(t)=x(t)+c(t);

[0078] Among them, f(t) is the seismic signal with noise added, x(t) is the seismic signal without noise, c(t) is the added random noise, and t is the time parameter;

[0079] Step S2, using particle swarm optimization to search for the optimal decomposition level K and penalty factor α of VMD decomposition;

[0080] v i,j (t+1)=ωv i,j (t)+c1rand1(pBest i,j -x ij (t))

[0081] +c2rand2(gBest i,j -x i,j (t))

[0082] x i,j (t+1)=x i,j (t)+v i,j (t+1)

[0083] In the formula, v ij (t) is the velocity of particle i in dimension j at time t, ω is the inertia weight, which controls the influence of the particle's last velocity, and c1 is the particle's optimal position toward the individual pBest i,j Adjustment strength, x ij (t) is the position of the particle, that is, the decomposition level K and penalty value α in VMD, and c2 is the particle's adjustment gBest towards the global optimal position. i,jstrength, rand1 and rand2 are random numbers, pBest i,j is the individual optimal position of particle i in dimension j, gBest i,j is the global optimal position of particle j in dimension i, that is, the optimal solution;

[0084] Step S3: Perform VMD decomposition on the seismic signal f(t) to obtain K IMF functions u with different center frequency and limited bandwidth. k (t), such as Figure 3 As shown, the constraint equations of VMD decomposition are as follows:

[0085]

[0086] Among them, u k (t) represents the K IMF components obtained by decomposition; ω k Represents the frequency center of each component; * is the convolution symbol; represents the derivative of the function with respect to time t; δ(t) is the unit pulse function; represents the square of the gradient 2 norm;

[0087] Further, the method includes step S301:

[0088] In order to obtain the optimal solution of the above constraint equation, the penalty factor α and the Lagrangian multiplication operator λ are introduced to obtain the augmented Lagrangian function L:

[0089]

[0090] Where L is the Lagrangian function, α is the penalty factor, Indicates the derivative of the function with respect to time t, ω k represents the frequency center of each component, δ(t) is the unit pulse function, Indicates the derivative of the function with respect to time t, * is the convolution symbol, λ(t) is the Lagrangian operator, represents the square of the gradient 2 norm;

[0091] The Lagrangian function is solved by the alternating direction multiplier algorithm, and u is obtained by iteration. k 、ω k ;

[0092]

[0093] Where ^ represents Fourier transform, λ is the Lagrange multiplication operator, and u k is the modal function, ω k is the center frequency, n represents the number of iterations;

[0094] Determine whether to stop iteratively updating the modal function uk :

[0095]

[0096] In the formula, n represents the number of iterations, u k is the modal function, θ is the given discrimination accuracy, and ω is the center frequency;

[0097] Step S4: Calculate the kurtosis values ​​of the K modal functions, compare the kurtosis values ​​of the K modal functions with a set kurtosis threshold, and screen out the modal functions containing noise and those not containing noise;

[0098]

[0099] In the formula, m i is the i-th sample value of the seismic signal, N is the number of sample points of the seismic signal, μ is the mean value of the signal, and kurtosis is the kurtosis value;

[0100] Among them, if the modal function u k (t) is less than the set kurtosis value threshold, it is judged to be a modal function containing noise. If the modal function u k If (t) is not less than the set kurtosis value threshold, it is judged to be a modal function that does not contain noise;

[0101] Step S5: Figure 5 As shown in the figure, the original modal components of the modal functions that are not less than the set kurtosis value threshold are retained, and the modal functions that are less than the set kurtosis value threshold are subjected to improved wavelet transform denoising.

[0102] Get the approximate function of the seismic signal and the wavelet function;

[0103]

[0104] Among them, the scale variable a controls the expansion and contraction of the wavelet function, and the translation variable b controls the translation of the wavelet function. is the wavelet basis function, u k (t) is the mode function, W u (a, b) are approximate functions obtained by inner product of seismic signal and wavelet basis;

[0105] Perform multi-scale analysis on the selected modal functions and use appropriate wavelet basis functions to calculate the scale coefficient c of the modal function j,i and wavelet coefficient d j,i ;

[0106]

[0107] In the formula, c j,i is the low-frequency approximate component at scale j and position i, d j,iis the high-frequency detail component at scale j and position i, φ j,i (t) is the scaling function, ψ j,i (t) is the wavelet function, u(t) is the modal function;

[0108] Set the number of layers of wavelet decomposition and calculate the threshold of each layer. The threshold function calculation formula is:

[0109]

[0110] σ=media(|d j,N |) / 0.6745;

[0111] Among them, σ is the noise variance of the signal to be processed, N is the number of wavelet coefficients with higher frequencies under the scale change, j represents the decomposition scale, media(|d j,N |) is the median of the absolute values ​​of all wavelet coefficients at scale j, and ρ is the threshold value;

[0112] Compare the detail coefficients obtained at each layer with the threshold value, and quantize the detail coefficients that are less than the threshold value. The threshold function is expressed as:

[0113]

[0114] Among them, d * j,i is the quantized detail coefficient, d j,i is the original detail coefficient, ρ n is the threshold value of each layer, sign is the sign function;

[0115] According to the threshold function, the original detail coefficient d j,i Perform quantization processing to obtain the quantized detail coefficient d * j,i , and then reconstruct to get the denoised modal function u k (t) * :

[0116]

[0117] Among them, u k (t) * is the mode function after denoising, c i is the approximation coefficient, φ(tk) is the scaling function, ψ(2 j tk) is an approximate function, d * j,i is the quantized detail coefficient, c i is the approximate coefficient;

[0118] Step S6: reconstruct the denoised modal function and the retained modal function to obtain the denoised seismic signal, such as Figure 6 As shown:

[0119] ∑ k u k (t)+u k (t) * =f(t);

[0120] Among them, u k (t) is the unprocessed mode function, u k (t) * is the processed modal function.

[0121] It should be noted that this paper proposes a seismic signal denoising method based on the PSO-VMD joint improved wavelet transform. For random noise, the analysis and calculation results show that the VMD denoising effect is affected by the decomposition layer number K and the penalty factor α. When the K value is too large, the decomposition result will have a mixing phenomenon, and when the K value is too small, it will affect the focus of the time-frequency mapping, which is not conducive to capturing the center frequency of the mode. The penalty factor α controls the data fidelity. An inappropriate α value will lead to the loss of modal components. The signal-to-noise ratio (SNR) of the denoised signal obtained using VMD is 10.97dB. Figure 4 As shown;

[0122] This method automatically selects the optimal K and α values ​​through the particle swarm optimization (PSO) algorithm, then uses VMD to decompose the signal to obtain multiple modal functions, and then performs denoising through an improved wavelet transform. The final signal-to-noise ratio (SNR) is 13.47dB. This method not only improves the denoising effect, but also effectively retains more valid signals and removes more noise. It has high stability and robustness. Due to the intelligence and stability of the algorithm, this method has significant denoising advantages, especially in complex noise environments.

[0123] In the description of the present invention, it should be understood that the terms "upper", "lower", "bottom", "top", "front", "back", "inside", "outside", "left", "right", etc., indicating directions or positional relationships, are based on the directions or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and therefore cannot be understood as limiting the present invention.

[0124] Although the present invention is described herein with reference to specific embodiments, it should be understood that these embodiments are merely illustrative of the principles and applications of the invention. It should be understood that many modifications may be made to the illustrative embodiments, and that other arrangements may be devised, without departing from the spirit and scope of the invention as defined by the appended claims. It should be understood that the various dependent claims and features described herein may be combined in ways other than those described in the original claims. It should also be understood that features described in conjunction with individual embodiments may be employed in conjunction with other described embodiments.

Claims

1. A seismic signal denoising method based on PSO-VMD combined with improved wavelet transform, characterized in that: The steps include: Step S1, obtaining seismic signals with and without noise; Step S2, using particle swarm optimization to search for the optimal decomposition level K and penalty factor α of VMD decomposition; Step S3, performing VMD decomposition on the seismic signal to obtain K mode functions; Step S4: Calculate the kurtosis values ​​of the K modal functions, compare the kurtosis values ​​of the K modal functions with a set kurtosis threshold, and screen out the modal functions containing noise and those not containing noise; Step S5: retaining the original modal components of the modal functions that are not less than the set kurtosis value threshold, performing improved wavelet transform denoising on the modal functions that are less than the set kurtosis value threshold, and obtaining an approximate function of the seismic signal and the wavelet function; Step S6: reconstruct the denoised modal function and the retained modal function to obtain a denoised seismic signal.

2. The seismic signal denoising method of PSO-VMD combined with improved wavelet transform according to claim 1 is characterized in that: In step S1, the original noise-free seismic signal x(t) and the noisy seismic signal f(t)=x(t)+c(t) are read in, where c(t) is the noise and t is the time parameter, to obtain the noisy and noise-free seismic signals.

3. The seismic signal denoising method of PSO-VMD combined with improved wavelet transform according to claim 2 is characterized in that: In step S2, the calculation formula is as follows: v i,j (t+1)= ωv i,j (t)+c1rand1(pBest i,j -x ij (t))+c2rand2(gBest i,j -x i,j (t)); x i,j (t+1)=x i,j (t)+v i,j (t+1); Among them, v ij (t) is the velocity of particle i in dimension j at time t, ω is the inertia weight, which controls the influence of the particle's last velocity, and c1 is the particle's optimal position toward the individual pBest i,j Adjustment strength, x ij (t) is the position of the particle, that is, the decomposition level K and penalty value α in VMD, and c2 is the particle's adjustment gBest towards the global optimal position. i,j strength, rand1 and rand2 are random numbers, pBest i,j is the individual optimal position of particle i in dimension j, gBest i,j is the global optimal position of particle i in dimension j, that is, the optimal solution.

4. The seismic signal denoising method of PSO-VMD combined with improved wavelet transform according to claim 3 is characterized in that: Step S3 includes: Step S301: Perform VMD decomposition on the seismic signal f(t) to obtain K IMF functions u with different center frequency and limited bandwidth. k (t), the constraint equations of VMD decomposition are as follows: s.t.∑ k u k (t)=f(t); Among them, u k (t) represents the K IMF components obtained by decomposition, ω k Represents the frequency center of each component, * is the convolution symbol, It means to find the derivative of the function with respect to time t, δ(t) is the unit pulse function, represents the square of the gradient 2 norm; Step S302: introducing a penalty factor α and a Lagrangian multiplication operator λ to obtain an augmented Lagrangian function L; Where L is the Lagrangian function, α is the penalty factor, Indicates the derivative of the function with respect to time t, ω k represents the frequency center of each component, δ(t) is the unit pulse function, Indicates the derivative of the function with respect to time t, * is the convolution symbol, and λ(t) is the Lagrange operator. represents the square of the gradient 2 norm; Step S303: Solve the Lagrangian function using the alternating direction multiplier algorithm and iterate to obtain u k 、ω k ; Where ^ represents Fourier transform, λ is the Lagrange multiplication operator, and u k is the modal function, ω k is the center frequency, n represents the number of iterations; Step S304: Determine whether to stop iterating and updating the modal function u k , the judgment formula is as follows: Among them, n represents the number of iterations, u k is the modal function, θ is the given discrimination accuracy, and ω is the center frequency.

5. The seismic signal denoising method of PSO-VMD combined with improved wavelet transform according to claim 4 is characterized in that: In step S4, the kurtosis values ​​of the K modal functions are calculated, including: Among them, m i is the i-th sample value of the seismic signal, N is the number of sample points of the seismic signal, μ is the mean of the signal, and kurtosis is the kurtosis value.

6. The seismic signal denoising method of PSO-VMD combined with improved wavelet transform according to claim 5 is characterized in that: In step S4, the modal functions containing noise and excluding noise are screened out, including: If the modal function u k (t) is less than the set kurtosis value threshold, it is judged to be a modal function containing noise. If the modal function u k If (t) is not less than the set kurtosis value threshold, it is judged to be a modal function that does not contain noise.

7. The seismic signal denoising method of claim 6 using PSO-VMD combined with improved wavelet transform, characterized in that: In step S5, the improved wavelet transform is performed on the modal functions with a kurtosis value less than a preset threshold value to remove noise, and an approximate function of the seismic signal and the wavelet function is obtained, including: Among them, the scale variable a controls the expansion and contraction of the wavelet function, and the translation variable b controls the translation of the wavelet function. is the wavelet basis function, u k (t) is the mode function, W u (a, b) is the approximate function obtained by doing the inner product of the seismic signal and the wavelet basis.

8. The seismic signal denoising method of PSO-VMD combined with improved wavelet transform according to claim 7 is characterized in that: Step S5 also includes performing multi-scale analysis on the selected modal functions and calculating the scale coefficient c of the modal function using a suitable wavelet basis function. j,i and wavelet coefficient d j,i , Among them, c j,i is the low-frequency approximate component at scale j and position i, d j,i is the high-frequency detail component at scale j and position i, φ j,i (t) is the scaling function, ψ j,i (t) is the wavelet function, and u(t) is the modal function.

9. The seismic signal denoising method of PSO-VMD combined with improved wavelet transform according to claim 8, characterized in that: Step S5 further includes: Set the number of wavelet decomposition layers and calculate the threshold for each layer, including: σ=media(|d j,N |) / 0.6745; Among them, σ is the noise variance of the signal to be processed, N is the number of wavelet coefficients with higher frequencies under the scale change, j represents the decomposition scale, media(|d j,N |) is the median of the absolute values ​​of all wavelet coefficients at scale j, and ρ is the threshold value; Quantize the detail coefficients that are smaller than the threshold, including: Among them, d * j,i is the quantized detail coefficient, d j,i is the original detail coefficient, ρ n is the threshold value of each layer, sign is the sign function; For the original detail coefficient d j,i Perform quantization processing to obtain the quantized detail coefficient d * j,i , and then reconstruct to get the denoised modal function u k (t) * ,include: Among them, u k (t) * is the mode function after denoising, c i is the approximation coefficient, φ(tk) is the scaling function, ψ(2 j tk) is an approximate function, d * j,i is the quantized detail coefficient, c i is the approximate coefficient.

10. The seismic signal denoising method of PSO-VMD combined with improved wavelet transform according to claim 9, characterized in that: Step S6 includes: The denoised mode function and the unprocessed mode function are reconstructed to obtain the denoised seismic signal: ∑ k u k (t)+u k (t) * =f(t); In the formula, u k (t) is the retained mode function, u k (t) * is the mode function after denoising.

Images