Improved EKF (Extended Kalman Filter) and wavelet packet collaborative ultrasonic echo signal joint noise reduction method

Abstract

The invention provides an ultrasonic echo signal joint noise reduction method based on improved EKF and wavelet packet cooperation, belongs to the technical field of ultrasonic detection, and solves the problems that the noise reduction effect is poor in a traditional ultrasonic signal noise reduction method and noise reduction and feature retention are difficult to balance in the wavelet packet noise reduction process. Modeling the ultrasonic echo signal and the noise signal to obtain an ultrasonic echo signal with noise; extracting a target ultrasonic echo signal by using an improved extended Kalman filter to obtain an observation equation; predicting and updating on the basis of the observation equation to obtain ultrasonic echo signal estimation; optimizing the noise covariance matrix and the observation noise covariance matrix by adopting a differential evolution algorithm to obtain an optimal parameter, executing extended Kalman filter filtering based on the optimal parameter, and generating an ultrasonic echo signal estimated value after preliminary denoising; and carrying out wavelet packet threshold denoising to obtain a final denoised ultrasonic echo signal.

Description

Technical Field

[0001] The present invention relates to an improved ultrasonic echo signal joint denoising method combining EKF and wavelet packet collaboration, belonging to the technical field of ultrasonic testing. Background Art

[0002] As an important means in non-destructive testing (NDT), ultrasonic testing is widely used in the identification and evaluation of internal defects in materials. After exciting ultrasonic waves, by collecting and analyzing the echo signals, tiny defects such as cracks, pores, and inclusions can be effectively detected, thus ensuring the structural integrity and service life of the materials. However, in practical applications, the collected ultrasonic signals are often accompanied by various noise interferences, severely restricting the detection and evaluation of defects. For example, when detecting coarse-grained materials by ultrasonic waves, due to the presence of a large number of reflecting grain boundaries in the materials, a large amount of noise generated by grain scattering, i.e., particle noise, often appears in the echo signals, resulting in the defect echo being possibly submerged in the noise signals. And the frequency band presented by this kind of noise is very similar to that of the target echo signal, making it impossible to be eliminated by classical denoising methods (such as time averaging or linear filtering). Therefore, effectively removing the noise in the ultrasonic echo signals becomes a key step in improving the detection accuracy.

[0003] Shi Qian, Li Qiufeng, Zhou Ruiqi et al. disclosed in the invention patent CN103901115A "An ultrasonic testing method for coarse-grained materials based on EMD combined with wavelet threshold denoising". This method combines empirical mode decomposition (EMD) and wavelet threshold denoising techniques. First, the ultrasonic signal is adaptively decomposed into multiple IMFs by EMD, and each IMF reflects the different frequency band characteristics of the signal. Then, wavelet threshold denoising is applied to the low-frequency IMF components to effectively reduce the noise while retaining the useful signals. This method is particularly suitable for environments with low signal-to-noise ratios, can significantly improve the signal quality and signal-to-noise ratio, and makes the signal characteristics originally masked by noise prominent. However, when screening the effective signal IMF components, when the high-frequency part overlaps with the target frequency band, a part of valuable content will be lost because the high-frequency IMF components are completely removed.

[0004] Du Qiaoling, Zhang Bin and others disclosed "An Ultrasonic Echo Denoising Method Based on CEEMDAN-Wavelet" in Patent CN118035647A. This method combines Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) and wavelet threshold denoising technology. First, the noisy ultrasonic signal is decomposed into multiple Intrinsic Mode Functions (IMFs) by CEEMDAN, and the sample entropy, correlation coefficient and root mean square error of each IMF are calculated. The entropy weight method is used to rank the comprehensive scores of these indicators, and the IMF components with more effective signals are selected. Then, the selected IMF components are further processed by wavelet threshold denoising. Finally, the processed IMF components are superimposed and reconstructed to complete denoising. This method improves the adaptability of ultrasonic signal processing and is applicable to ultrasonic detection in complex environments. However, since the noise content is high and the signal components are complex in relatively high-frequency bands (such as IMF0 and IMF1), simply using CEEMDAN and wavelet threshold denoising technology may not be able to effectively retain the content of high-frequency effective signals. In addition, the wavelet threshold denoising method has lower efficiency and accuracy compared to the wavelet packet threshold denoising when processing multi-scale or multi-band signals such as ultrasonic echo signals. In terms of the threshold selection method, the soft threshold will cause the signal to be over-smoothed, and the hard threshold will cause signal energy loss. Summary of the Invention

[0005] The present invention aims to solve the problems of poor noise reduction effect in traditional ultrasonic signal noise reduction and the difficulty in balancing denoising and feature retention in the wavelet packet denoising process, and further proposes a combined ultrasonic echo signal joint denoising method that combines an improved EKF and wavelet packet.

[0006] The technical solution adopted by the present invention to solve the above problems is as follows: The present invention includes the following steps:

[0007] Step 1: Collect ultrasonic echo signals, model the ultrasonic echo signals and noise signals to obtain the noisy ultrasonic echo signal s(t);

[0008] Step 2: Use the improved extended Kalman filter to extract the target ultrasonic echo signal to obtain the observation equation;

[0009] Step 3: Based on the observation equation, perform prediction and update to obtain the estimated value of the ultrasonic echo signal

[0010] Step 4: Use the differential evolution algorithm to optimize the noise covariance matrix Q k and the observation noise covariance matrix R k to obtain the optimal parameters ((Q * , R * ), and based on the optimal parameters (Q * , R *The preliminary estimated value of the ultrasonic echo signal Perform extended Kalman filter filtering to generate an estimated value of the ultrasonic echo signal after preliminary denoising

[0011] Step 5: The estimated value of the ultrasonic echo signal after preliminary denoising Perform wavelet packet threshold denoising to obtain the finally denoised ultrasonic echo signal

[0012] Preferably, step 1 specifically includes:

[0013] Step 1.1: Obtain the functional form of the analog signal and use it as the original signal f(t);

[0014] Step 1.2: Obtain the frequency-domain representation of the received signal, where the frequency-domain received signal includes defect reflection signals, scattering noise, and white noise;

[0015] Step 1.3: Transform the original signal f(t) to the frequency through fast Fourier transform, and calculate the frequency vector f according to the number of sampling points n and time step Δt of the original signal f(t);

[0016] Step 1.4: Generate scattering noise N1(f) that conforms to a Gaussian distribution and the frequency-domain representation of the scattering noise according to the single-scattering model;

[0017] Step 1.5: Perform conjugate symmetry filling on the negative frequency part of the scattering noise and white noise f < 0;

[0018] Step 1.6: Generate additional white noise N2(f) and perform conjugate symmetry filling, and superimpose the conjugate symmetry-filled scattering noise and white noise in the frequency domain to obtain the frequency-domain representation of the total noise;

[0019] Step 1.7: Transform the frequency-domain representation of the total noise to the time domain through inverse Fourier transform to obtain the time-domain noise Noise(t);

[0020] Step 1.8: Adjust the noise amplitude of the time-domain noise Noise(t) until the target signal-to-noise ratio is reached to obtain the time-domain noise Noise scaled (t);

[0021] Step 1.9: Scale the time-domain noise Noise scaled (t) and add it to the original signal f(t) to obtain the ultrasonic echo signal s(t) with noise;

[0022] The expression of the functional form of the analog signal is:

[0023]

[0024] In Equation (1), A1 and A2 represent amplitude factors, a1 and a2 are bandwidth factors, τ1 and τ2 are waveform arrival times, f1 and f2 are the center frequencies of the signals, and φ1 and φ2 are phases;

[0025] The expression for the frequency-domain representation of the received signal is:

[0026]

[0027] In Equation (2), is the reflected signal from the defect, N1(f) is the scattering noise obeying a Gaussian distribution, which reflects the scattering effect caused by internal particles of the material, H(f) is the frequency response of the transducer, α0 is the attenuation coefficient of the material, and N2(f) is the additionally introduced white noise;

[0028] The calculation formula for the frequency response H(f) of the transducer is:

[0029]

[0030] In Equation (3), f0 is the center frequency and σ is the bandwidth;

[0031] The expression for the frequency-domain representation of the scattering noise is:

[0032]

[0033] The expression for conjugate filling of the negative frequency part of the scattering noise and the white noise is:

[0034] ScatteringNoiseFull(f) = ScatteringNoise*(-f) (5);

[0035] The expression for conjugate filling of the additional white noise N2(f) is:

[0036]

[0037] The expression for the frequency-domain representation of the total noise:

[0038] Noise f (f) = ScatteringNoiseFull(f) + WhiteNoiseFull(f) (7);

[0039] The calculation formula for the ultrasonic echo signal s(t) with noise is:

[0040] s(t) = f(t) + Noise scaled (t) (8).

[0041] Preferably, Step 2 specifically includes:

[0042] Step 2.1: Define the state vector x of the extended Kalman filter k , where the state vector includes the amplitude, frequency, phase, frequency change rate, frequency acceleration, and frequency attenuation degree that can all reflect signal information;

[0043] Step 2.2: Based on the state vector x k Establish the state transition equation x of the system at the time step Δt k+1 , and expand the state transition equation x k+1 to obtain the state space model;

[0044] Step 2.3: Establish the observation equation to map the state vector to the representation of the observed signal;

[0045] The expression of the state vector x k is:

[0046]

[0047] In formula (9), A 1,k and A 2,k are the amplitudes of two frequency components, f 1,k and f 2,k are the frequencies of two frequency components, φ 1,k and φ 2,k are the corresponding phases respectively, and are the frequency change rates of two signal components, and are the frequency accelerations of two signal components, α 1,k and α 2,k are the amplitude attenuation rates of two signal components;

[0048] The expression of the state transition equation x k+1 is:

[0049] x k+1 = f(x k , Δt) + w k (10);

[0050] The expression of the state space model x k+1 is:

[0051]

[0052] In formulas (10) and (11), f(·) is the nonlinear state transition function used to describe the change of the state vector over time, w k is the process noise, representing the random fluctuation in the change of the system state, and satisfies w k ~ N(0, Q k ), among the variables, the amplitude Ai Exponentially decay over time at a damping rate α i with a frequency f i vary over time at a rate of change of frequency and a frequency acceleration Update, with a phase φ i Update based on the current frequency and frequency acceleration, with the rate of change of frequency vary over time at a frequency acceleration Update, with the frequency acceleration and an amplitude decay rate α i Set to a constant;

[0053] The expression of the observation equation is:

[0054] z k = h(x k ) + v k = A 1,k cos(φ 1,k ) + A 2,k cos(φ 2,k ) + v k (12);

[0055] In formula (12), z k is the observation signal at time step k, h(·) is the non - linear observation function used to describe how the state vector generates the observation signal, and v k is the observation noise, which follows a normal distribution N(0, R k ).

[0056] Preferably, step 3 specifically includes:

[0057] Step 3.1: After establishing the observation equation, enter the prediction phase. For each time step k, use the state transition equation to predict the state vector and covariance matrix of the observation function at the next moment. Apply the prediction phase to the system model to obtain the prior estimate of the current state and complete the prediction process;

[0058] Step 3.2: After the prediction is completed, enter the update phase. Take the newly predicted observation value z k as the signal observation value at the current time step, calculate the difference between the signal observation value and the predicted value at the current time step to obtain the innovation y k ;

[0059] Step 3.3: Calculate the difference distribution between the observation value and the predicted value at the current time step to obtain the innovation covariance S k , where the innovation covariance S k is the variance of y k ;

[0060] Step 3.4: Based on the covariance matrix, observation function, and innovation covariance S at the current moment k Calculate the Kalman gain K k , and based on the Kalman gain K k Update the new state vector and covariance matrix to complete the update process;

[0061] Step 3.5: Obtain all updated state vectors and construct a signal estimate value based on all updated state vectors

[0062] The update expression for the prior estimate is:

[0063]

[0064] The update expression for the covariance matrix is:

[0065]

[0066] In formulas (13) and (14), is the prior estimate at the k-th time step, is the prior covariance matrix at the k-th time step, is the Jacobian matrix of the state transition function, is the process noise covariance matrix;

[0067] The innovation y k is calculated as:

[0068]

[0069] In formula (15), is the observation prediction value based on the predicted state vector;

[0070] The innovation covariance S k is calculated as:

[0071]

[0072] In formula (16), is the Jacobian matrix of the observation function, which can linearize the nonlinear system, is the observation noise covariance matrix;

[0073] The Kalman gain K k is calculated as:

[0074]

[0075] The expressions for updating the new state vector and covariance matrix are:

[0076]

[0077] In formulas (18) and (19), is the posterior estimate at the k-th time step, is the posterior covariance matrix at the k-th time step, and I is the identity matrix;

[0078] Signal estimated value The calculation formula of is:

[0079]

[0080] In formula (20), A 1,k and A 2,k are respectively the final estimated values of the amplitudes of two frequency components, and φ 1,k and φ 2,k are respectively the final estimated values of the corresponding phases.

[0081] Preferably, step 4 specifically includes:

[0082] Step 4.1: During the optimization process, with the goal of maximizing the signal-to-noise ratio SNR out of the filtered output, the negative value of the output signal-to-noise ratio, Objective(Q, R) = -SNR out is used as the objective function, and the combination of the parameter vectors Q k and R k is used as the population, and each individual combination (Q i , R i ) in the population is used as the output value;

[0083] Step 4.2: Define the search ranges of the parameter vectors Q k and R k , and initialize the vector as the starting point of the population for differences, where N is the population size;

[0084] Step 4.3: Perform a mutation operation to generate a mutation vector v i for each set of parameters (Q i , R i,G );

[0085] Step 4.4: If the random number rand(0, 1) < C r , then generate an experimental vector u i,G through a crossover operation; if rand(0, 1) > C r , then directly copy the current individual as the experimental vector, i.e., u i,G = x i,G ;

[0086] Step 4.5: Based on the preset crossover probability Cr Perform a selection operation to compare the fitness values of the experimental vector and the current vector. If Objective(u i,G ) < Objective(x i,G ), then x i,G+1 = u i,G . If Objective(u i,G ) > Objective(x i,G ), then x i,G+1 = x i,G ;

[0087] Step 4.6: Repeat Steps 4.3 - 4.5 until the maximum number of iterations is reached or the fitness value converges to obtain the optimal parameters (Q * , R * ) that minimize the objective function. Use the optimal parameters (Q * , R * ) to perform extended Kalman filter filtering to generate an estimated value of the pre - denoised ultrasonic echo signal

[0088] Output the calculation formula for the signal - to - noise ratio SNR out is:

[0089]

[0090] In formula (21), f is the original signal, is the estimated signal after filtering;

[0091] The mutation vector v i,G 's calculation formula is:

[0092] v i,G = x r1,G + F·(x r2,G - x r3, G) (22);

[0093] In formula (22), x r1,G , x r2,G , x r3,G are different individuals randomly selected from the population, and F is the scaling factor.

[0094] Preferably, Step 5 specifically includes:

[0095] Step 5.1: Perform wavelet packet transform on the estimated value of the pre - denoised ultrasonic echo signal to obtain the detailed coefficients of the pre - denoised ultrasonic echo signal in all frequency bands;

[0096] Step 5.2: Obtain the wavelet packet coefficient c in the detailed coefficients through the Gaussian mixture modelj,k Distribution;

[0097] Step 5.3: Estimate the parameters of the Gaussian mixture model according to the expectation-maximization algorithm, and obtain the optimal parameters of the Gaussian mixture model by combining the maximized likelihood function. Calculate the posterior probability that each wavelet packet coefficient c belongs to the noise component based on the Gaussian mixture model and the optimal parameters; j,k

[0098] Step 5.4: Define a non-linear threshold function λ(c j,k ) according to the posterior probability. The non-linear threshold function λ(c j,k ) adaptively adjusts the threshold size according to the local statistical characteristics of the signal;

[0099] Step 5.5: Apply the non-linear threshold function λ(c j,k ) to each wavelet packet coefficient c j,k ) for processing to obtain the coefficient c'; j,k ;

[0100] Step 5.6: After completing the non-linear threshold processing, use the inverse wavelet packet transform to reconstruct the processed coefficient c′ j,k into the finally denoised signal

[0101] The estimated value of the preliminarily denoised ultrasonic echo signal The expression of the wavelet packet transform representation is:

[0102]

[0103] In formula (23), ψ j,k (t) is the basis function of the wavelet packet transform, c j,k are the coefficients of each frequency band, j is the decomposition level, J is the maximum decomposition level, and k is the frequency band index;

[0104] The expression of the distribution of the wavelet packet coefficient c j,k is:

[0105]

[0106] In formula (24), M is the number of Gaussian components, π m is the mixing weight of the m-th Gaussian component, satisfying is a Gaussian distribution with a mean of μ m and a variance of ;

[0107] The calculation formula for the posterior probability is:

[0108] ​

[0109] In formula (25), N is the set of Gaussian components corresponding to the noise;

[0110] The expression of the non-linear threshold function λ(c j,k ) is:

[0111] (c j,k ) = α·σ m ·Φ -1 {1 - P(noise|c j,k ) (26);

[0112] In formula (26), α is the adjustment factor, σ m is the standard deviation of the noise component, and Φ -1 is the inverse cumulative distribution function of the standard normal distribution;

[0113] The calculation formula for the coefficient c' j,k is:

[0114]

[0115] The finally denoised signal The calculation formula for it is:

[0116]

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

[0118] (1) Compared with other conventional noise reduction methods, this method uses the EKF, based on the recursive process of prediction and update, and accurately estimates the state of the signal through the dynamic model. Compared with the modal decomposition method, it can better retain the high-frequency signal components and effectively reduce these non-linear granular noises.

[0119] (2) This method globally optimizes the process noise covariance matrix Q and the observation noise covariance matrix R of the EKF through the differential evolution algorithm to ensure that the EKF can still maintain a high signal-to-noise ratio in a high-frequency noise environment.

[0120] (3) Combining the dynamic state estimation of the EKF and the multi-band denoising technology of the wavelet packet, while effectively reducing systematic and non-Gaussian noises, it can accurately distinguish and retain high-frequency effective information through the multi-scale decomposition of the wavelet packet, thereby significantly improving the signal-to-noise ratio.

[0121] (4) In the selection of the threshold method for wavelet packet threshold denoising, this method uses the non-linear threshold function method based on the Gaussian mixture model. Compared with the commonly used soft threshold or hard threshold, this method can achieve a better balance between noise reduction and feature retention. Description of the Drawings

[0122] Figure 1 Schematic diagram of the process of an improved EKF combined with wavelet packet collaborative ultrasonic echo signal joint denoising method provided by the present invention;

[0123] Figure 2 Effect comparison diagram of the signal denoising of the present invention and other denoising methods;

[0124] Figure 3 Different signal processing stages of the present invention under the condition of 10 dB noise. Specific implementation manner

[0125] Combined with Figure 1 and Figure 3 Describe this implementation manner. As Figure 1 shown, the steps of an improved EKF combined with wavelet packet collaborative ultrasonic echo signal joint denoising method described in this implementation manner include:

[0126] S1: Collect ultrasonic echo signals and model the ultrasonic echo signals and noise signals;

[0127] S101: During the ultrasonic detection process, the signals received by the receiving transducer are ultrasonic echo signals. These echo signals not only contain the reflected waves caused by defects, but may also include interface waves returned from the material boundary. Therefore, in order to be closer to the actual collected echo signal form, this implementation manner uses a signal composed of two oscillation components to model the actual collected echo signals, and the echo signals can be modeled as time-shifted Gaussian signals; therefore, the functional form f(t) of the analog signal can be expressed as:

[0128]

[0129] In formula (1), A1 and A2 represent amplitude factors, a1 and a2 are bandwidth factors, τ1 and τ2 are waveform arrival times, f1 and f2 are the center frequencies of the signals, and φ1 and φ2 are phases;

[0130] S102: Particle noise has a significant impact on the signals in actual ultrasonic detection, and its characteristics are composed of two parts: scattering noise N1(f) and white noise N2(f) in the frequency domain; the expression of the received signal in the frequency domain is:

[0131]

[0132] In formula (2), is the reflected signal from the defect, N1(f) is the scattering noise obeying the Gaussian distribution, reflecting the scattering effect caused by particles inside the material, H(f) is the frequency response of the transducer, α0 is the attenuation coefficient of the material, and N2(f) is the additional white noise introduced;

[0133] S103: Add the modeled frequency-domain noise to the original signal f(t). The original signal needs to be transformed to the frequency domain through the Fast Fourier Transform (FFT), and then the frequency vector f is calculated according to the number of sampling points n and the time step Δt of the original signal. The frequency response H(f) of the transducer is a Gaussian filter centered at the center frequency f0 with a bandwidth of σ, and its expression is:

[0134]

[0135] In formula (3), f0 is the center frequency and σ is the bandwidth;

[0136] S104: Generate the scattering noise N1(f) that conforms to the Gaussian distribution according to the single-scattering model. Its real and imaginary parts are independent and identically distributed Gaussian random variables. The scattering noise is calculated in the frequency domain through the following formula:

[0137]

[0138] S105: Since the ultrasonic echo signal is real and its Fourier transform has conjugate symmetry, it is necessary to perform conjugate symmetric filling on the negative frequency parts of the scattering noise and white noise; for f < 0, the expression is as follows:

[0139] ScatteringNoiseFull(f) = ScatteringNoise * (-f) (5);

[0140] S106: Generate additional white noise N2(f), whose frequency-domain representation is independent of frequency, and also perform conjugate symmetric filling. The formula is as follows:

[0141]

[0142] S107: Superimpose the scattering noise and white noise in the frequency domain to obtain the frequency-domain representation of the total noise as shown below:

[0143] Noise f (f) = ScatteringNoiseFull(f) + WhiteNoiseFull(f) (7);

[0144] S108: Transform the frequency-domain representation of the total noise back to the time domain through the inverse Fourier transform to obtain the time-domain noise Noise(t). In the further verification and evaluation process, in order to ensure that the added noise reaches the target signal-to-noise ratio (SNR), it is necessary to adjust the amplitude of the noise. The adjusted noise signal is represented as Noise scaled (t), and the specific adjustment steps are provided in the experimental section. Scale the noise and add it to the original signal to obtain the ultrasonic echo signal s(t) with noise. The formula is as follows:

[0145] s(t) = f(t) + Noise scaled (t) (8).

[0146] S2: Define the state - space model and the observation equation;

[0147] This embodiment uses the Extended Kalman Filter (EKF) as the core algorithm to initially extract the target echo signal and further improve it. Compared with the modal decomposition noise reduction method and the wavelet transform noise reduction method, EKF can achieve stronger robustness in the face of complex and non - Gaussian noise while retaining a more complete frequency band.

[0148] As a recursive algorithm, EKF updates the state vector and covariance matrix continuously through linearization according to the observed quantity and the current state estimate at each time step, and calculates the Kalman gain, so as to estimate the details and trends of the system state, specifically including:

[0149] S201: Define the state vector x of EKF k , which contains all relevant variables that can reflect signal information, such as amplitude, frequency, phase, frequency change rate, frequency acceleration, and frequency attenuation degree, as follows:

[0150]

[0151] In formula (9), a 1,k and A 2,k are the amplitudes of two frequency components, f 1,k and f 2,k are the frequencies of two frequency components, φ 1,k and φ 2,k are the corresponding phases respectively, and are the frequency change rates of two signal components, and are the frequency accelerations of two signal components, α 1,k and α 2,k are the amplitude attenuation rates of two signal components;

[0152] S202: The state vector reflects the variable state at the current time point, while the state - transition equation describes the changes of these variables over time. The specific state - transition equation x k+1 As shown in formula (10), the further extended state - space model can be expressed as formula (11), which describes the state evolution of the system after the time step Δt;

[0153] x k+1 = f(x k , Δt) + w k(10);

[0154]

[0155] In equations (10) and (11), f(·) is a non-linear state transition function used to describe the change of the state vector over time, and w k is the process noise, which characterizes the random fluctuations in the system state change and satisfies w k ~N(0, Q k ). Among the variables, the amplitude A i exponentially decays over time at a damping rate α i , the frequency f i is updated over time at a frequency change rate and a frequency acceleration , the phase φ i is updated based on the current frequency and frequency acceleration. The frequency change rate is updated over time at a frequency acceleration , and the frequency acceleration and the amplitude decay rate α i are set to constants;

[0156] S203: After defining the state transition equation, an observation equation needs to be established next to describe the observation process of the system. The observation equation maps the state vector to the observation signal and is expressed as equation (12); since each echo in the signal can be modeled as a time-shifted Gaussian echo wavelet, the observation equation is constructed as a model composed of the superposition of two frequency components, which is used to reflect the actual measured signal;

[0157] z k = h(x k ) + v k = A 1,k cos(φ 1,k ) + A 2,k cos(φ 2,k ) + v k (12);

[0158] In equation (12), z k is the observation signal at time step k, h(·) is the non-linear observation function used to describe how the state vector generates the observation signal, and v k is the observation noise, which follows a normal distribution N(0, R k ).

[0159] S3: Predict the state vector and covariance matrix at the next time step;

[0160] In this embodiment, the process of EKF recursive estimation is achieved by combining the prediction step and the update step, continuously optimizing the estimation of the state vector, thereby improving the accuracy of the estimation of the original signal model, including the following steps:

[0161] S301: Enter the prediction stage. For each time step k, use the state transition equation to predict the state vector and covariance matrix at the next moment, as shown in formula (13). The prediction step provides a prior estimate of the current state through the application of the system model, providing a basis for the subsequent update step:

[0162]

[0163] In formulas (13) and (14), is the prior estimate at the k-th time step, is the prior covariance matrix at the k-th time step, is the Jacobian matrix of the state transition function, is the process noise covariance matrix;

[0164] S302: After the prediction step is completed, enter the update stage. In this stage, once a new observation value z k , that is, the signal observation value at the current time step, is obtained, EKF will adjust the state vector and covariance matrix. First, calculate the innovation y k , that is, the difference between the observation value and the predicted value, as shown in formula (15):

[0165]

[0166] In formula (15), is the observation prediction value based on the predicted state vector;

[0167] Next, calculate another important variable in the prediction stage, that is, the innovation covariance S k (InnovationCovariance), which is the variance of the innovation y k , used to measure the distribution of the difference between the observation value and the predicted value, as shown in formula (16):

[0168]

[0169] In formula (16), is the Jacobian matrix of the observation function, which can linearize the nonlinear system, is the observation noise covariance matrix;

[0170] S303: In the EKF, the Kalman gain is the weight balance of the prediction and the observed value in the state estimation, which is crucial for optimizing the accuracy of the state estimation. With a reasonable Kalman gain, the EKF can maintain the prediction accuracy and effectively utilize the new observed information. After obtaining the covariance matrix, the observation function, and the innovation covariance, the Kalman gain Kk is further calculated, and the formula is shown in (17):

[0171]

[0172] S304: After calculating the Kalman gain K k the new state vector and covariance matrix are updated:

[0173]

[0174] In formulas (18) and (19), is the posterior estimate at the k-th time step, is the posterior covariance matrix at the k-th time step, and I is the identity matrix;

[0175] S305: After the EKF process ends, all updated state vectors can be obtained, that is, the relevant variables (such as amplitude, frequency, phase, etc.); through these variables, the final signal estimate value is constructed and obtained, and the formula is shown in (20):

[0176]

[0177] In formula (20), A 1,k and A 2,k are the final estimated values of the amplitudes of the two frequency components respectively, φ 1,k and φ 2,k are the final estimated values of the corresponding phases respectively. However, the EKF depends on an accurate noise model, including the estimation of the process noise Q k and the observation noise R k . If the estimation of these noise covariances is inaccurate, it may lead to over-suppression or under-suppression of the noise by the filter, resulting in a large amount of noise remaining in the output signal.

[0178] S4: The differential evolution algorithm is used to optimize the parameters Q k and R k to obtain the optimal parameters (Q * , R * ). Based on the optimal parameters (Q * , R * ), the extended Kalman filter is executed for filtering to generate the estimated value of the preliminarily denoised ultrasonic echo signal

[0179] To further improve the noise reduction performance of the EKF, this embodiment uses the differential evolution algorithm to optimize the process noise covariance matrix Q involved in the EKF k and the observation noise covariance matrix R k ; The differential evolution algorithm (DE) is a population-based heuristic search method for global optimization, especially suitable for high-dimensional, non-linear and non-convex optimization problems, and can effectively explore the parameter space to find the global optimal solution.

[0180] S401: In this optimization, the goal is to maximize the signal-to-noise ratio SNR of the filtered output out , which is defined as follows:

[0181]

[0182] In formula (21), f is the original signal, is the estimated signal after filtering;

[0183] S402: To achieve maximization, this embodiment sets the objective function to the negative value of the output signal-to-noise ratio, that is, Objective(Q,R) = -SNR out ; By maximizing the signal-to-noise ratio, it can be ensured that the filtered signal is as close as possible to the original signal and the noise is suppressed;

[0184] S403: To solve the problem of inaccurate noise model estimation, this embodiment uses the differential evolution (DE) algorithm to optimize the process noise covariance matrix Q in the prediction and update steps of the EKF k and the observation noise covariance matrix R k ; The DE algorithm gradually approaches the optimal solution by maintaining a population composed of multiple candidate solutions and applying mutation and crossover operations to generate new solutions in each generation. Each individual in the population represents a possible solution, and the algorithm improves the overall optimization effect through the cooperation and competition among individuals, including the following steps:

[0185] S40301: Construct a population composed of multiple individuals. In the DE algorithm, the population is a set composed of several parameter vectors, and each individual represents a possible solution;

[0186] S40302: Define the search ranges of Q k and R k , and initialize a set of parameter vectors as the starting point of the differential population, where N is the population size;

[0187] S40303: Perform mutation operations. For each set of parameters (Q i , R i ), generate the mutated vector v i,G :

[0188] v i,G = x r1,G + F·(x r2,G - x r3,G ) (22);

[0189] In formula (22), x r1,G , x r2,G , x r3,G are different individuals randomly selected from the population, and F is the scaling factor;

[0190] S40304: If the random number rand(0,1) < C r , then generate the experimental vector u i,G through the crossover operation; if rand(0,1) > C r , then directly copy the current individual as the experimental vector, that is, u i,G = x i,G ;

[0191] S40305: To ensure the diversity and effectiveness of the trial vectors, in this embodiment, the crossover condition is based on a preset crossover probability C r . In the selection operation, by comparing the fitness values of the trial vector and the current vector, determine the individuals in the next generation population. If Objective(u i,G ) < Objective(x i,G ), then x i,G+1 = u i,G ; otherwise x i,G+1 = x i,G . The selection operation ensures that the fitness of individuals in each generation does not decrease;

[0192] S40306: Repeat S40303 - S40305 for repeated evaluation until the maximum number of iterations is reached or the fitness value converges, to obtain the optimal parameters (Q * , R * ) that minimize the objective function. Use the optimal parameters (Q * , R * ) to perform extended Kalman filter filtering to generate an estimated value of the pre - denoised ultrasonic echo signal Thus, the best noise reduction effect is achieved. Through parameter optimization, the filter can adaptively adjust in different noise environments, ensuring the robustness and efficiency of the filtering effect. Especially in a high - frequency noise environment, it can still retain effective signal information.

[0193] The EKF combines predicted and actual observation data, gradually corrects the state estimate, and achieves preliminary restoration of the original signal, which is also equivalent to the denoising process. However, since errors are introduced during the linearization of the non-linear system by the EKF, certain noise remains in the frequency domain.

[0194] S5: Perform wavelet packet threshold denoising and output the finally denoised ultrasonic echo signal

[0195] Wavelet packet denoising (WPD) is used to further eliminate the residual noise of the EKF. WPD decomposes the signal into sub-signals of different frequency bands through multi-scale decomposition, effectively distinguishing the signal characteristics of each scale. At each decomposition level, the coefficients of each frequency band are processed by thresholding to further eliminate the noise not completely removed by the EKF, including the following steps:

[0196] S501: For the signal to be processed after the EKF Perform wavelet packet transform to obtain the detailed coefficients of the signal in each frequency band. The signal representation of the wavelet packet transform is given in formula (23):

[0197]

[0198] In formula (23), ψ j,k (t) is the basis function of the wavelet packet transform, c j,k is the coefficient of each frequency band, j is the decomposition level, J is the maximum decomposition level, and k is the frequency band index;

[0199] S502: According to the non-linear threshold function method of the Gaussian mixture model (GMM), the threshold method of the wavelet packet denoising method can be further optimized. Through the flexible modeling ability of the GMM for the multi-modal distribution of wavelet packet coefficients, the wavelet packet decomposition can more effectively distinguish between noise and signal components. The distribution of the wavelet packet coefficient c j,k can be described by the Gaussian mixture model, and the formula is as shown in (24):

[0200]

[0201] In formula (24), M is the number of Gaussian components, π m is the mixing weight of the m-th Gaussian component, satisfying is a Gaussian distribution with a mean of μ m and a variance of ;

[0202] S503: The parameters of the Gaussian mixture model Estimation is performed by the Expectation-Maximization (EM) algorithm. The EM algorithm can obtain the optimal model parameters by maximizing the likelihood function to describe the statistical characteristics of wavelet packet coefficients; subsequently, based on the GMM model, calculate the posterior probability that each coefficient c j,k belongs to the noise component:

[0203]

[0204] In formula (25), N is the set of Gaussian components corresponding to the noise;

[0205] S504: In this embodiment, a non-linear threshold function λ(c j,k ) is defined according to the posterior probability, and the formula is as shown in (26):

[0206] (c j,k ) = α·σ m ·Φ -1 {1 - P(noise|c j,k ) (26);

[0207] In formula (26), α is a regulation factor, σ m is the standard deviation of the noise component, Φ -1 is the inverse cumulative distribution function of the standard normal distribution. This threshold function can adaptively adjust the threshold size according to the local statistical characteristics of the signal, so as to achieve the best noise reduction effect under different frequency bands and different signal intensities;

[0208] S505: After obtaining the threshold function, apply the non-linear threshold function to each wavelet packet coefficient c j,k . This non-linear threshold processing can not only effectively suppress the noise component, but also retain the main features of the signal, avoiding the signal distortion problem that may be caused by the traditional hard threshold method. The specific formula is as shown in (27)

[0209]

[0210] S506: After completing the non-linear threshold processing, use the inverse wavelet packet transform to reconstruct the processed coefficient c' j,k into the denoised signal The specific formula is as shown in (28):

[0211]

[0212] In S5, the multi-layer decomposition ability of wavelet packet transform enables it to meticulously analyze the frequency domain characteristics of signals. By introducing the non-linear threshold function method based on Gaussian mixture model into wavelet packet threshold denoising, a better balance between denoising and feature preservation is achieved compared to hard and soft thresholds, ensuring that the edges and details of the signal after EKF filtering are retained, thereby significantly enhancing the overall signal-to-noise ratio and clarity of the signal. Figure 3 The four images from top to bottom respectively show different signal processing stages under 10 dB noise conditions: the first image is the original signal, the second image is the signal with added noise, the third image is the signal after applying EKF, and the last image shows the signal effect after combining EKF and wavelet packet threshold denoising technology. These images intuitively demonstrate the improvement of signal quality at each step of processing.

[0213] Embodiment

[0214] Combination Figure 1-2 This embodiment will be described as follows. Figure 1 As shown, an improved EKF combined with wavelet packet collaborative ultrasonic echo signal joint denoising method described in this embodiment includes:

[0215] First, this embodiment models an ultrasonic echo simulation signal for testing the denoising effect. This signal is composed of two exponentially decaying cosine waves superimposed, and the specific expression is as follows:

[0216] s(t) = 2·e -3|t - τ| ·cos(2π·10·(t - τ)) + e -6|t-2 |·cos(2π·10·(t - 2)) (29);

[0217] In formula (29), τ = 7 seconds marks the arrival time of the first oscillating component, and t is the time variable; this embodiment sets the time axis range to [-1, 15] seconds and divides it into 1000 sampling points, where the time step is dt = t[1] - t[0]; the component with a larger amplitude and relatively lower frequency is used to represent the end face wave reflected from the material boundary, and the component with a smaller amplitude and relatively higher frequency is used to represent the defect wave.

[0218] Next, to simulate the actually received signal, this embodiment adds granular noise to the simulation signal; in the program, the granular noise is achieved by combining scattering noise and white noise in the frequency domain. The scattering noise simulates the signal attenuation at a central frequency of 10 Hz, with a standard deviation of 2 Hz and a very small attenuation coefficient of 1×10 -6To control the frequency distribution and attenuation rate, white noise provides randomness with a uniform frequency distribution; after synthesizing these two types of noise in the frequency domain, the noise amplitude is adjusted to achieve a set signal-to-noise ratio, ranging from 0 dB to 10 dB, and then inverse Fourier transformed back to the time domain to form the final noisy signal.

[0219] After receiving the noisy signal, in this embodiment, the EKF is used for preliminary noise reduction processing; first, the initial state vector and state covariance matrix of the EKF are set. In this embodiment, more general initial values are used and the covariance is increased to reflect uncertainty. Therefore, the initial state vector x0 is set as:

[0220]

[0221] Let this state vector be x k The system state is updated and predicted recursively through time; through the state transition equation and the observation equation, the EKF continuously optimizes the state estimation, effectively processes non-Gaussian noise, and retains the original form of the signal; finally, the updated state vector is used to reconstruct the signal, and parameters such as the amplitude and phase of each frequency component are estimated to construct the estimated original form of the signal.

[0222] During the iteration process, the parameter optimization of the EKF is achieved through the differential evolution algorithm. A cost function needs to be defined, which receives the process noise covariance matrix Q and the measurement noise covariance matrix R of the EKF as input parameters. Specifically, the search ranges of Q and R are set. In this example, they are defined as Q ∈ [10 -6 , 10 -2 and R ∈ [10 -4 , 10 -1 . The SNR between the denoised signal and the original signal is calculated during each iteration of the EKF process. The differential evolution algorithm searches for the optimal Q and R within the preset parameter range to achieve the highest signal-to-noise ratio.

[0223] Next, in order to further filter out the residual noise, this embodiment adopts a non-linear threshold wavelet packet denoising technique based on the Gaussian mixture model: first, the signal is decomposed by wavelet packet using the'sym4' wavelet with a maximum decomposition level of 6; then, the signal coefficients obtained by decomposition are standardized, and a Gaussian mixture model containing 6 components is fitted to these normalized coefficients. By comparing the mean of each component with the overall mean, the components below the average are identified as noise; next, the threshold is dynamically adjusted by estimating the noise probability of each coefficient, and then the non-linear threshold function λ(c j,k ) is used to selectively shrink or eliminate the coefficients, and the reconstructed signal is obtained through inverse wavelet packet transform. This method not only effectively suppresses the noise components but also better retains the main features of the signal, significantly improving the signal-to-noise ratio.

[0224] In addition, to verify the technical effects of this embodiment, the following simulation experiment was designed:

[0225] Additive white Gaussian noise with different intensities was added to the analog signal. The addition of granular noise was carried out based on the signal-to-noise ratio (SNR), which gradually increased from 0 dB to 10 dB. The implementation method was to first calculate the average power P of the generated initial granular noise n(t) noiseinitial :

[0226]

[0227] To adjust the noise power to the required P noise , calculate the scaling factor α and scale the noise n(t), so as to add noises with different energies (signal-to-noise ratios) to the original signal:

[0228]

[0229] n scaled (t) = α·n(t) (33);

[0230] s(t) = f(t) + n scaled (t) (34);

[0231] After adding noises of 0 to 10 decibels to the signal, in order to comprehensively evaluate the performance of different noise reduction methods, this embodiment selected a variety of mainstream noise reduction technologies as control objects; specifically including the empirical mode decomposition (EMD) noise reduction method, the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) noise reduction method, the combined CEEMDAN and wavelet packet threshold (WPD) noise reduction method, and only the EKF noise reduction method, to compare with the method proposed in this embodiment. Using the above five methods respectively, the same original signal f1(t) was denoised, and the signal-to-noise ratio SNR of the denoised signal and the original signal was calculated respectively out to compare the effects.

[0232] After the experiment was completed, in order to visually display the noise reduction effect of this method under different input noise levels, this embodiment plotted the relationship curve between the input signal-to-noise ratio (SNR_in) and the average output signal-to-noise ratio (SNR_out). In addition, for a comprehensive performance evaluation, this embodiment compared and analyzed the proposed method with several of the above ultrasonic echo signal noise reduction methods, and the relevant results are as Figure 2 shown, clearly showing the noise reduction performance of each method under different input noise levels, Figure 2Clearly demonstrates the superior performance of the method of the present invention compared to other methods under different SNR input conditions. The blue line represents the noise reduction method of the present invention, showing a significant improvement in the signal-to-noise ratio under all test conditions; Table 1 shows the signal-to-noise ratio (SNR_out) of the signals generated by various noise reduction methods under different noise conditions.

[0233] Table 1

[0234]

[0235] As can be seen from Table 1, the EKF + wavelet packet threshold noise reduction method of the present invention achieves the highest signal-to-noise ratio under all noise conditions, significantly superior to other methods. By calculating the average value, after adding 0 to 10 dB of granular noise, the average signal-to-noise ratio of the method proposed in this study is increased by about 12.2% compared to the mainstream CEEMDAN combined with wavelet packet threshold noise reduction method.

[0236] In summary, the improved joint noise reduction method combining EKF and wavelet packet proposed by the present invention shows significant advantages in ultrasonic echo signal noise reduction; through the signal form restoration technology of EKF combined with the multi-scale signal processing technology of wavelet packet, this method realizes efficient noise removal and effective retention of signal features. The introduction of the DE algorithm improves the accuracy of EKF prediction, and the introduction of the Gaussian mixture model further improves the accuracy of threshold processing, making the noise reduction effect more superior.

[0237] The above is only a preferred embodiment of the present invention and does not impose any form of limitation on the present invention. Although the present invention has been disclosed above with a preferred embodiment, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to the above-disclosed technical content to form equivalent embodiments within the scope of the technical solution of the present invention. However, as long as it does not depart from the content of the technical solution of the present invention and is based on the technical essence of the present invention, any simple modification, equivalent replacement, and improvement of the above embodiments are still within the protection scope of the technical solution of the present invention.

Claims

1. An improved EKF combined with wavelet packet synergy ultrasonic echo signal joint denoising method, characterized in that: The steps of the improved EKF combined with wavelet packet collaborative ultrasonic echo signal joint denoising method include: Step 1: Collect ultrasonic echo signals, model the ultrasonic echo signals and noise signals, and obtain ultrasonic echo signals with noise s(t); Step 2: Use the improved extended Kalman filter to extract the target ultrasonic echo signal and obtain the observation equation; Step 3: Predict and update based on the observation equation to obtain a preliminary estimate of the ultrasonic echo signal Step 4: Use differential evolution algorithm to calculate the noise covariance matrix Q in the extended Kalman filter process k and the observation noise covariance matrix R k Optimize and obtain the optimal parameters (Q * ,R * ), based on the optimal parameter (Q * ,R * ) is a preliminary estimate of the ultrasonic echo signal Perform extended Kalman filter filtering to generate an estimate of the ultrasonic echo signal after preliminary denoising Step 5: Estimation of the ultrasonic echo signal after preliminary denoising Perform wavelet packet threshold denoising to obtain the final denoised ultrasonic echo signal 2. The method for joint denoising of ultrasonic echo signals by combining improved EKF with wavelet packet synergy according to claim 1, characterized in that: Step 1 specifically includes: Step 1.1: Get the functional form of the analog signal and use it as the original signal f(t); Step 1.2: Obtain a frequency domain representation of the received signal, where the frequency domain received signal includes a defect reflection signal, scattered noise, and white noise; Step 1.3: Convert the original signal f(t) to frequency through fast Fourier transform, and calculate the frequency vector f according to the number of sampling points n and the time step Δt of the original signal f(t); Step 1.4: Based on the single scattering model, generate the scattering noise N1(f) that conforms to the Gaussian distribution and the frequency domain representation of the scattering noise; Step 1.5: Perform conjugate symmetric filling on the negative frequency parts of the scattered noise and white noise f<0; Step 1.6: Generate additional white noise N2(f) and perform conjugate symmetric filling, superimpose the conjugate symmetric filled scattered noise and white noise in the frequency domain, and obtain the frequency domain representation of the total noise; Step 1.7: Convert the frequency domain representation of the total noise into the time domain representation by inverse Fourier transform to obtain the time domain noise Noise(t); Step 1.8: Adjust the noise amplitude of the time domain noise Noise(t) until the target signal-to-noise ratio is reached to obtain the time domain noise Noise scaled (t); Step 1.9: Convert the time domain noise scaled (t) is scaled and added to the original signal f(t) to obtain an ultrasonic echo signal s(t) with noise; The expression of the function form of the analog signal is: In formula (1), A1 and A2 represent amplitude factors, a1 and a2 are bandwidth factors, τ1 and τ2 are waveform arrival times, f1 and f2 are the center frequencies of the signals, and φ1 and φ2 are phases; The frequency domain expression of the received signal is: In formula (2), is the reflection signal from the defect, N1(f) is the scattered noise obeying Gaussian distribution, reflecting the scattering effect caused by particles inside the material, H(f) is the frequency response of the transducer, α0 is the attenuation coefficient of the material, and N2(f) is the additional white noise introduced; The frequency response H(f) of the transducer is calculated as: In formula (3), f0 is the center frequency and σ is the bandwidth; The frequency domain expression of scattered noise is: The expression for conjugate filling of the negative frequency part of scattered noise and white noise is: ScatteringNoiseFull(f)=ScatteringNoise * (-f) (5); The expression of the additional white noise N2(f) conjugate filling is: The expression of total noise in frequency domain is: Noise f (f)=ScatteringNoiseFull(f)+WhiteNoiseFull(f) (7); The calculation formula of ultrasonic echo signal s(t) with noise is: s(t)=f(t)+Noise scaled (t) (8)。 3. The method for joint denoising of ultrasonic echo signals by combining improved EKF with wavelet packet synergy according to claim 1, characterized in that: Step 2 specifically includes: Step 2.1: Define the state vector x of the extended Kalman filter k , where the state vector includes all the amplitude, frequency, phase, frequency change rate, frequency acceleration and frequency attenuation that can reflect the signal information; Step 2.2: Based on the state vector x k Establish the state transfer equation x of the system at time step Δt k+1 , for the state transfer equation x k+1 Expand to obtain the state space model; Step 2.3: Establish the observation equation and map the state vector to the observation signal representation; The state vector x k The expression is: In formula (9), A 1,k and A 2,k are the amplitudes of the two frequency components, f 1,k and f 2,k are the frequencies of the two frequency components, φ 1,k and φ 2,k are the corresponding phases, and is the frequency change rate of the two signal components, and is the frequency acceleration of the two signal components, α 1,k and α 2,k is the amplitude attenuation rate of the two signal components; State transfer equation x k+1 The expression is: x k+1 =f(x k ,Δt)+w k (10); State space model x k+1 The expression is: In formulas (10) and (11), f(·) is the nonlinear state transfer function, which is used to describe the change of the state vector over time, and w k is process noise, which characterizes the random fluctuations in the system state change and satisfies w k ~N(0,Q k ), in the variables, the amplitude A i With time, the damping rate α i Exponential decay, frequency f i Frequency change rate over time and frequency acceleration Update, phase φ i Update based on current frequency and frequency acceleration, frequency change rate Acceleration with frequency over time Update, frequency acceleration and the amplitude attenuation rate α i Set to a constant; The expression of the observation equation is: z k =h(x k )+v k =A 1,k cos(φ 1,k )+A 2,k cos(φ 2,k )+v k (12); In formula (12), z k is the observation signal at time step k, h(·) is the nonlinear observation function, which is used to describe how the state vector generates the observation signal, v k is the observation noise, which conforms to the normal distribution N(0,R k ).

4. The method for joint denoising of ultrasonic echo signals by combining improved EKF with wavelet packet synergy according to claim 1, characterized in that: Step 3 specifically includes: Step 3.1: After establishing the observation equation, enter the prediction phase. For each time step k, use the state transfer equation to predict the state vector and covariance matrix of the observation function at the next moment. Apply the prediction phase to the system model to obtain the prior estimate of the current state and complete the prediction process. Step 3.2: After the prediction is completed, enter the update phase and update the predicted new observation value z k As the signal observation value at the current time step, calculate the difference between the signal observation value and the predicted value at the current time step to obtain the innovation y k ; Step 3.3: Calculate the difference distribution between the observed value and the predicted value at the current time step to obtain the innovation covariance S k , where the innovation covariance S k for y k The variance of Step 3.4: Based on the covariance matrix, observation function and innovation covariance S at the current moment k Calculate the Kalman gain K k , based on the Kalman gain K k Update the new state vector and covariance matrix to complete the update process; Step 3.5: Get all updated state vectors and construct signal estimates based on all updated state vectors The update expression of the prior estimate is: The update expression of the covariance matrix is: In formulas (13) and (14), is the prior estimate of the kth time step, is the prior covariance matrix of the kth time step, is the Jacobian matrix of the state transfer function, is the process noise covariance matrix; Innovation k The calculation formula is: In formula (15), is the observed predicted value based on the predicted state vector; Innovation Covariance S k The calculation formula is: In formula (16), is the Jacobian matrix of the observation function, which can linearize the nonlinear system. is the observation noise covariance matrix; Kalman gain K k The calculation formula is: The expression for updating the new state vector and covariance matrix is: In formulas (18) and (19), is the posterior estimate at the kth time step, is the posterior covariance matrix of the kth time step, and I is the identity matrix; Signal Estimation The calculation formula is: In formula (20), A 1,k and A 2,k are the final estimates of the amplitudes of the two frequency components, φ 1,k and φ 2,k are the final estimated values ​​of the corresponding phases, respectively.

5. The method for joint denoising of ultrasonic echo signals by combining improved EKF with wavelet packet according to claim 1, characterized in that: Step 4 specifically includes: Step 4.1: During the optimization process, the output signal-to-noise ratio (SNR) after filtering is maximized. out As the target, the negative value of the output signal-to-noise ratio is Objective(Q,R)=-SNR out As the objective function, the parameter vector Q k and R k As a population, each individual combination (Q i ,R i ) as the output value; Step 4.2: Define the parameter vector Q k and R k The search range, initialization vector As the starting point of the population for difference, where N is the population size; Step 4.3: Perform mutation operation for each set of parameters (Q i ,R i ) Generate the mutation vector v of the corresponding parameters i,G ; Step 4.4: If the random number rand(0,1) <C r , then the experimental vector u is generated by crossover operation i,G ; If rand(0,1)>C r , then directly copy the current individual as the experimental vector, that is, u i,G =x i,G ; Step 4.5: Based on the preset crossover probability G r Perform a selection operation and compare the fitness value of the experimental vector with the current vector. If Objective(u i,G )<Objective(x i,G ), then x i,G+1 =u i,G , if objective(u i,G )>Objective(x i,G ), then x i,G+1 =x i,G ; Step 4.6: Repeat steps 4.3 to 4.5 until the maximum number of iterations is reached or the fitness value converges, and the optimal parameter (Q * ,R * ), using the optimal parameter (Q * ,R * ) performs extended Kalman filter filtering to generate an estimate of the ultrasonic echo signal after preliminary denoising Output signal-to-noise ratio SNR out The calculation formula is: In formula (21), f is the original signal, is the estimated signal after filtering; Mutation vector v i,G The calculation formula is: v i,G =x r1,G +F·(x r2,G -x r3,G ) (22); In formula (22), x r1,G 、x r2,G 、x r3,G are different individuals randomly selected from the population, and F is the scaling factor.

6. The method for joint denoising of ultrasonic echo signals by combining improved EKF with wavelet packet synergy according to claim 1, characterized in that: Step 5 specifically includes: Step 5.1: Estimation of the ultrasonic echo signal after preliminary denoising Perform wavelet packet transform to obtain the ultrasonic echo signal after preliminary denoising Detailed coefficients in all frequency bands; Step 5.2: Obtain the wavelet packet coefficients c in the detailed coefficients through the Gaussian mixture model j,k Distribution of Step 5.3: Parameters of the Gaussian mixture model according to the expectation maximization algorithm Estimate and combine the maximum likelihood function to obtain the optimal parameters of the Gaussian mixture model, and calculate each wavelet packet coefficient c based on the Gaussian mixture model and the optimal parameters j,k The posterior probability of belonging to the noise component; Step 5.4: Define a nonlinear threshold function λ(c j,k ), nonlinear threshold function λ(c j,k ) Adaptively adjust the threshold value according to the local statistical characteristics of the signal; Step 5.5: For each wavelet packet coefficient c j,k Applying a nonlinear threshold function λ(c j,k ) is processed to obtain the coefficient c' j,k ; Step 5.6: After completing the nonlinear threshold processing, use the inverse wavelet packet transform to transform the processed coefficients c' j,k Reconstruct the final denoised signal Estimated value of ultrasonic echo signal after preliminary denoising The expression of wavelet packet transform is: In formula (23), ψ j,k (t) is the basis function c of wavelet packet transform j,k is the coefficient of each frequency band, j is the number of decomposition levels, J is the maximum number of decomposition levels, and k is the frequency band index; Wavelet packet coefficient c j,k The distribution expression is: In formula (24), M is the number of Gaussian components, π m is the mixing weight of the mth Gaussian component, satisfying The mean is μ m , the variance is Gaussian distribution of The formula for calculating the posterior probability is: In formula (25), N is the set of Gaussian components corresponding to the noise; Nonlinear threshold function λ(c j,k ) is: (c j , k )=a·s m ·F -1 {1-P(noise|c j,k ) (26); In formula (26), α is the adjustment factor, σ m is the standard deviation of the noise component, Φ -1 is the inverse cumulative distribution function of the standard normal distribution; Coefficient c' j,k The calculation formula is: The final denoised signal The calculation formula is:

Images