A denoising method for non-stationary and non-linear magnetotelluric sounding signals

Abstract

The present invention provides a denoising method for non-stationary and non-linear magnetotelluric sounding signals, which relates to the technical field of geophysical exploration. The method includes obtaining magnetotelluric sounding signals from a magnetotelluric sounding instrument, selecting one channel from the obtained magnetotelluric sounding signals, denoising it, and repeating the above operations to perform denoising processing on the remaining obtained magnetotelluric sounding signals. By combining variational mode decomposition (VMD) and independent component analysis (ICA), the present invention realizes the denoising of magnetotelluric sounding signals from multiple perspectives of time-frequency domain and statistical characteristics. When the effective signal is interfered by non-stationary noise with high energy for a long time, this method can ensure the high quality of the reconstructed signal, effectively reduce the influence of artificial parameter selection, be robust to noise, and at the same time provide a clear physical explanation, overcoming the limitations of traditional denoising methods and providing a more reliable and accurate data processing means for the field of geophysical exploration.

Description

Technical Field

[0001] The present invention relates to the technical field of geophysical exploration, and particularly to a denoising method for non-stationary and non-linear magnetotelluric sounding signals. Background Art

[0002] The magnetotelluric (MT) sounding method, as an effective geophysical exploration technology, has always been highly regarded by scientific researchers and has been widely used in the research of the crust and upper mantle structures, the exploration of oil and gas mineral resources, and earthquake monitoring. The MT method records the time variation of the natural electromagnetic field on the earth's surface and deduces the conductivity distribution of the underground medium using the principle of electromagnetic induction. Although the MT method has a relatively deep penetration depth and high resolution, its signal is extremely vulnerable to natural and man-made noises during the acquisition process, seriously reducing the signal-to-noise ratio of MT data and possibly causing distortion of the apparent resistivity and phase curves, affecting the inversion and interpretation of magnetotelluric data. Common man-made noises such as high-voltage power lines and radio transmission towers will produce long-term strong interference on the observation of magnetotelluric sounding signals. Therefore, researching effective denoising technologies is of great significance for improving the quality of MT data and enhancing the interpretation accuracy.

[0003] Since the magnetotelluric method was proposed by Tikhonov and Cagniard in the 1950s, significant progress has been made in research on noise suppression, and a variety of methods have been derived. According to the different processing objects, magnetotelluric denoising techniques can be divided into three categories: time domain, frequency domain, and time-frequency domain: 1. Traditional methods such as the cross-power spectrum method with a distant reference channel and the Robust statistical method, which are both implemented in the impedance estimation process in the frequency domain. They can respectively achieve the removal of uncorrelated noise and suppress the influence brought by the extreme values of "fly points"; 2. Time domain analysis methods such as time domain impedance estimation based on Bayesian estimation, which updates the estimation of parameters by combining prior knowledge and observed data, is a framework for statistical inference under uncertainty. Compared with the frequency domain method based on Fourier transform, it has more advantages in dealing with non-stationary data; The mathematical morphology filtering technology has the ability to accurately identify specific noise types. When the useful signal and noise share the same frequency band, this method shows better robustness compared with the traditional band-pass filter; In addition, machine learning methods such as artificial neural networks, through their powerful pattern recognition and feature extraction capabilities, can capture the local features of signals from massive data, automatically learn and adjust model parameters to accurately denoise signals in the time domain; 3. Time-frequency domain analysis methods include wavelet analysis and Hilbert-Huang Transform (HHT). Wavelet analysis decomposes signals into components of different time scales through the stretching and translation of wavelet basis functions. Hilbert-Huang Transform adaptively obtains intrinsic mode functions through empirical mode decomposition and obtains the time-frequency spectrum through Hilbert transform. Both methods can obtain time-frequency characteristics. When estimating impedance, combining distant reference data and robust estimation can effectively suppress coherent noise while maintaining a high time-frequency resolution, so they have received extensive attention and research from relevant researchers.

[0004] Although the above methods are all feasible in theory, when dealing with the field measured data of magnetotelluric sounding, due to the non-stationary, non-linear and non-minimum phase characteristics of the signals, the traditional Fourier transform-based spectral analysis has certain limitations. In the time-domain denoising methods, Bayesian estimation depends on the setting of prior knowledge, which may be based on experience or subjective judgment and may affect the objectivity of the final results. In addition, in practical applications, due to the complexity of the magnetotelluric system and incomplete data information, it may be difficult to select an appropriate Bayesian model; the effect of mathematical morphology filtering strongly depends on the selection of the structuring element, which requires professional knowledge and experience, and different structuring elements may be needed for different signals; machine learning methods for training deep learning models usually require a large amount of computing resources, including high-performance GPUs and a large amount of memory. In addition, its internal working mechanism and decision-making process are difficult to explain, which may be a problem in the denoising of magnetotelluric data that requires clear physical explanations. During the time-frequency analysis process, the underground structure is complex and variable. The effect of wavelet decomposition depends on the selection of the wavelet basis in the current environment, and the threshold setting has a significant impact on the denoising results. Since time-frequency domain analysis methods are particularly suitable for processing non-stationary signals and can obtain detailed information at different time scales and frequency levels, and have unique advantages in identifying local features occurring in signals such as instantaneous events or transient phenomena, the EMD and its improved methods have become the mainstream methods in current time-frequency analysis and are widely used. However, empirical mode decomposition lacks a wide range of standardization and theoretical basis, which to a certain extent limits its accuracy and reliability in the decomposition process. Especially when dealing with the problem of mode mixing, EMD may encounter difficulties, and different signal modes may interfere with each other, resulting in the inability to clearly distinguish the mode components, and the decomposition results may not always be directly related to physical phenomena.

[0005] Therefore, it is very necessary to design a denoising method for non-stationary and non-linear magnetotelluric sounding signals. Summary of the Invention

[0006] In order to overcome the deficiencies of the prior art, the purpose of the present invention is to provide a denoising method for non-stationary and non-linear magnetotelluric sounding signals.

[0007] To achieve the above purpose, the present invention provides the following solutions:

[0008] The present invention provides a denoising method for non-stationary and non-linear magnetotelluric sounding signals, including:

[0009] Obtain magnetotelluric sounding signals from a magnetotelluric sounding instrument;

[0010] Select one channel from the obtained magnetotelluric sounding signals and denoise it;

[0011] Repeat the above operations to denoise the remaining magnetotelluric sounding signals obtained.

[0012] Preferably, obtain magnetotelluric sounding signals from a magnetotelluric sounding instrument, specifically:

[0013] Obtain 5-channel magnetotelluric sounding signals from a magnetotelluric sounding instrument. Among them, the first channel and the second channel are the measured values of the electric channels in the north-south direction and the east-west direction respectively, the third channel and the fourth channel are the measured values of the magnetic channels in the north-south direction and the east-west direction respectively, and the fifth channel is the measured value of the vertical magnetic channel.

[0014] Preferably, select one channel from the obtained magnetotelluric sounding signals and denoise it, specifically:

[0015] Perform VMD decomposition on the magnetotelluric signal, define the intrinsic mode function as an amplitude-frequency modulation and frequency modulation signal as the IMF to be optimized k component, that is, the k-th order intrinsic mode component;

[0016] Construct a constrained variational optimization problem for the intrinsic mode function and solve it;

[0017] Combined with Parseval's theorem and Fourier isometric transformation, transform it into the frequency domain, use the alternating direction multiplier method in matlab for iterative optimization, calculate the frequency domain expressions of each mode in the frequency domain, and obtain the time domain expressions of the IMFs with different center frequencies based on the inverse Fourier transform k component;

[0018] Based on the obtained intrinsic mode components of each order, construct the observation matrix of ICA;

[0019] Perform pre-whitening and centering operations on the observation matrix;

[0020] Construct an optimization function of the weight vector based on negative entropy and kurtosis, use the gradient ascent method to iterate the weight vector, repeat the iterative operation, obtain n weight vectors, form a demixing matrix, and calculate the estimated value of the signal source;

[0021] Evaluate the IMFs components after ICA processing from the aspects of morphology and time-frequency characteristics respectively, identify and suppress the noise components, retain the useful signal components, recombine the useful signal components, construct the denoised magnetotelluric sounding signal, and evaluate its performance.

[0022] Preferably, perform VMD decomposition on the magnetotelluric signal, define the intrinsic mode function as an amplitude-frequency modulation and frequency modulation signal as the IMF to be optimized k component, that is, the k-th order intrinsic mode component, specifically:

[0023] Define the intrinsic mode function as an amplitude-frequency modulation and frequency modulation signal as the IMF to be optimizedk The component, namely the k-th order eigenmode component, where and are respectively the instantaneous amplitude and instantaneous phase of the IMF k component, and the instantaneous frequency is obtained by differentiating the phase, which is:

[0024] ;

[0025] Use to fit the input signal.

[0026] Preferably, construct a constrained variational optimization problem to solve for the intrinsic mode function, specifically:

[0027] For construct a constrained variational optimization problem to solve, ensuring that each mode satisfies certain constraint conditions, which are:

[0028] ;

[0029] In the formula, is the k-th order mode component, is the instantaneous frequency, is an input magnetotelluric sounding signal, is the summation operation for k = 1, 2,..., K, is the operation of taking the partial derivative with respect to time, is the convolution operation, is the square of the second norm, is to find the minimum value, is the constraint condition, and δ(t) is the impulse function.

[0030] Preferably, combined with Parseval's theorem and Fourier isometric transformation, transform it into the frequency domain, and use the alternating direction multiplier method in matlab for iterative optimization to calculate the frequency domain expressions of each mode in the frequency domain, and obtain the time domain expressions of the IMF k components with different center frequencies, specifically:

[0031] Combined with Parseval's theorem and Fourier isometric transformation, transform it into the frequency domain, and use the alternating direction multiplier method in matlab for iterative optimization to make their energy in the frequency domain mainly concentrated in the predetermined frequency band, and obtain the and of each mode obtained in the frequency domain, where the expression after alternating optimization iteration is:

[0032] ;

[0033] ;

[0034] After Fourier transform, they are all functions of frequency ω at this time. The original x(t) and y(t) are transformed into x(ω) and y(ω). The subscript k represents the k-th IMF component, i represents the current iteration number as the i-th, the superscript n is the total iteration number, λ is the Lagrange multiplier used in the iterative calculation, α is the penalty factor, and ω k is the central frequency of the k-th IMF;

[0035] The IMF with different central frequencies is obtained by taking the inverse Fourier transform and taking the real part. Among them, the expression of the inverse Fourier transform is: k The time-domain expression of, where the expression of the inverse Fourier transform is:

[0036] ;

[0037] Select and update the modes according to the set number of modes K. When the threshold of the iteration number or the residual energy is satisfied, the iterative process ends, that is:

[0038] ;

[0039] In the formula, is the threshold, and finally K IMF components are obtained.

[0040] Preferably, based on the obtained eigenmode components of each order, construct the observation matrix of ICA, specifically:

[0041] Obtain the eigenmode components of each order in the time-frequency domain obtained by decomposition, and construct the observation matrix of ICA X = [IMF1, IMF2,..., IMF n , where n independent source signals are represented as S = [s1, s2,..., s n , and under the action of the mixing matrix A = [a1, a2,..., a n , the observed values X = [x1, x2,..., x n are obtained, which is represented by the following matrix equation:

[0042] .

[0043] Preferably, construct an optimization function of the weight vector based on negentropy and kurtosis, use the gradient ascent method to iterate the weight vector, repeat the iterative operation, obtain n weight vectors, form the demixing matrix, and calculate the estimated value of the signal source, specifically:

[0044] Construct an optimization function of the weight vector based on negentropy and kurtosis, which is:

[0045] ;

[0046] In the formula, represents kurtosis, which is defined as , represents negative entropy and is approximately , where is a non - quadratic function, is the standard normal distribution, is the operation of finding the mathematical expectation;

[0047] The gradient ascent method is used to iteratively update the weight vector;

[0048] Repeat the iterative operation to obtain n weight vectors, which are combined into the demixing matrix W = [w1, w2, …, w n , and according to the observed value, the amplitude recovery is performed by the least - squares method to obtain the estimated value of the signal source , where n is the number of estimated source signals, and each column vector of

[0049] Preferably, the IMFs components after ICA processing are evaluated respectively from the morphological and time - frequency characteristics, the noise components are identified and suppressed, the useful signal components are retained, the useful signal components are recombined to construct the denoised magnetotelluric sounding signal, and its performance is evaluated, specifically:

[0050] The IMFs after ICA processing are evaluated respectively from the morphological and time - frequency characteristics, the noise components are identified and suppressed, the useful signal components are retained, the useful signal components are recombined to construct the denoised magnetotelluric sounding signal, and the normalized cross - correlation coefficient, signal - to - noise ratio and reconstruction error are used as evaluation indexes to evaluate its performance.

[0051] According to the specific embodiments provided by the present invention, the following technical effects are disclosed:

[0052] The present invention provides a denoising method for non - stationary and non - linear magnetotelluric sounding signals. The method includes obtaining magnetotelluric sounding signals from a magnetotelluric sounding instrument, selecting one channel from the obtained magnetotelluric sounding signals for denoising, and repeating the above operation to perform denoising processing on the remaining obtained magnetotelluric sounding signals. By combining variational mode decomposition (VMD) and independent component analysis (ICA), the present invention realizes the denoising of magnetotelluric sounding signals from multiple perspectives of time - frequency domain and statistical characteristics. When the effective signal is interfered by non - stationary noise with high energy for a long time, this method can ensure the high quality of the reconstructed signal, effectively reduce the influence of artificial parameter selection, be robust to noise, and at the same time provide a clear physical interpretation, overcoming the limitations of traditional denoising methods and providing a more reliable and accurate data processing means for the field of geophysical exploration. The advantages of the present invention include:

[0053] 1. Through the variational mode decomposition technique, the present invention can adaptively decompose signals without prior knowledge or setting specific filters, which provides a flexible method for processing complex signals.

[0054] 2. By independent component analysis, the modal aliasing problem is further solved, showing significant effectiveness in dealing with the signal-to-noise mixture within the same frequency band. For long-existing and powerful interference noises, it can not only cope with the instantaneous changes of signals, but also stably identify and suppress interference on a long time scale, ensuring the quality and reliability of magnetotelluric sounding data denoising.

[0055] 3. The introduction of ICA reduces the dependence on the selection of VMD parameters, reduces the influence of human factors on the denoising results, enhances the robustness of the method, and enables it to work stably under complex conditions.

[0056] 4. The denoising process and results have more explicit physical meanings. Each independent signal component may be associated with a specific physical process. For example, high-frequency components may reveal shallow aquifers, low-frequency components may indicate deep rock structures, and the interference of power frequency noise is often at 50 Hz and its harmonics, making the denoised signal easier to be interpreted as specific geological structures or processes. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required in the embodiments. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0058] Figure 1 It is a schematic flow chart of the denoising method for non-stationary and non-linear magnetotelluric sounding signals provided by the embodiments of the present invention.

[0059] Figure 2 It is a schematic flow chart of the VMD decomposition process.

[0060] Figure 3 It is a schematic flow chart of the ICA separation process.

[0061] Figure 4 It is a schematic flow chart of the signal screening and reconstruction process. DETAILED DESCRIPTION OF THE INVENTION

[0062] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0063] By combining variational mode decomposition (VMD) and independent component analysis (ICA), the present invention realizes denoising of magnetotelluric sounding signals from multiple perspectives of time-frequency domain and statistical characteristics. When the effective signal is interfered by long-term high-energy non-stationary noise, this method can ensure the high quality of the reconstructed signal, effectively reduce the influence of artificial parameter selection, be robust to noise, and at the same time provide a clear physical explanation, overcoming the limitations of traditional denoising methods and providing a more reliable and accurate data processing means for the field of geophysical exploration.

[0064] To make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0065] Figure 1 The flowchart of the method provided for the embodiments of the present invention is as Figure 1 shown. The present invention provides a denoising method for non-stationary and non-linear magnetotelluric sounding signals, including:

[0066] Step 100: Obtain magnetotelluric sounding signals from a magnetotelluric sounding instrument;

[0067] Step 200: Select one channel from the obtained magnetotelluric sounding signals and perform denoising on it;

[0068] Step 300: Repeat the above operation to perform denoising on the remaining obtained magnetotelluric sounding signals.

[0069] In step 100, obtaining magnetotelluric sounding signals from a magnetotelluric sounding instrument specifically means:

[0070] Obtain 5 magnetotelluric sounding signals from a magnetotelluric sounding instrument. Among them, the first and second channels are the measured values of the electric channels in the north-south and east-west directions respectively, the third and fourth channels are the measured values of the magnetic channels in the north-south and east-west directions respectively, and the fifth channel is the measured value of the vertical magnetic channel.

[0071] In step 200, selecting one channel from the obtained magnetotelluric sounding signals and performing denoising on it specifically means:

[0072] Step 201: Perform VMD decomposition on the magnetotelluric sounding signals;

[0073] Step 202: Perform ICA separation on the signal after VMD decomposition;

[0074] Step 203: Perform signal screening and reconstruction on the signal after ICA separation;

[0075] As Figure 2 shown, in Step 201, perform VMD decomposition on the magnetotelluric sounding signal, specifically:

[0076] Step 2011: Define the intrinsic mode function as an amplitude - frequency - modulated signal as the IMF k component to be optimized, that is, the k - th order intrinsic mode component, where and are the instantaneous amplitude and instantaneous phase of the IMF k respectively, and the instantaneous frequency can be obtained by taking the derivative of the phase ), and fit the input signal with ;

[0077] Step 2012: Solve the constructed constrained variational optimization problem in Step 2011 to ensure that each mode satisfies certain constraint conditions, such as having higher energy near the center frequency of the frequency band:

[0078] ;

[0079] In the formula, is the k - th order mode component, is the instantaneous frequency, is an input magnetotelluric sounding signal of one channel, is the summation operation for k = 1, 2, …, K, is the operation of taking the partial derivative with respect to time, is the convolution operation, is the square of the two - norm, is to find the minimum value, is the constraint condition, and δ(t) is the impulse function (Dirac function);

[0080] Step 2013: Combine the Parseval theorem and Fourier isometric transformation, transform it into the frequency domain, and use the alternating direction multiplier method in matlab for iterative optimization, so that their energy in the frequency domain is mainly concentrated in the predetermined frequency band, and the and of each mode obtained in the frequency domain are obtained. Through the inverse Fourier transform and taking the real part, the time - domain expressions of the IMF k with different center frequencies are obtained to solve the optimization problem in Step 4;

[0081] Among them, the expression after alternating optimization iteration is:

[0082] ;

[0083] ;

[0084] After Fourier transform, they are all functions of frequency ω at this time. The original x(t) and y(t) are transformed into x(ω) and y(ω). The subscript k represents the k-th IMF component; i represents the current iteration number as the i-th; the superscript n is the total iteration number; λ is the Lagrange multiplier, which is used in the iterative calculation; α is the penalty factor; ω k is the central frequency of the k-th IMF;

[0085] The expression of the inverse Fourier transform is:

[0086] ;

[0087] Select and update the modes according to the set number of modes K. When the threshold of the iteration number or the residual energy is satisfied, the iterative process ends, that is:

[0088] ;

[0089] Finally, K IMF components are obtained.

[0090] As Figure 3 shown, in step 202, the signal after VMD decomposition is separated by ICA. Specifically:

[0091] Step 2021: Take each order of intrinsic mode components obtained from the signal decomposition in step 2013, that is, K IMF components, to construct the observation matrix X = [IMF1, IMF2,..., IMF n , assume that n independent source signals are represented as S = [s1, s2,..., s n , and the observed values X = [x1, x2,..., x n are obtained under the action of the mixing matrix A = [a1, a2,..., a n . It can be represented by the following matrix equation:

[0092] .

[0093] Step 2022: Perform pre-whitening and centering operations on the observed data matrix . Pre-whitening makes the algorithm easier to converge and can improve stability. Centering eliminates the DC component in the data, enabling the ICA algorithm to better identify and separate statistically independent signal sources;

[0094] Step 2023: Construct an optimization function for the weight vector based on negative entropy and kurtosis:

[0095] ;

[0096] In the formula, represents kurtosis, defined as , represents negative entropy, which can be approximated as , where is a non - quadratic function, is the standard normal distribution, is the operation of taking the mathematical expectation;

[0097] Use the gradient ascent method to iteratively update the weight vector;

[0098] Step 2024: Repeat the iterative operation in Step 2023 to obtain n weight vectors, where n is the number of estimated source signals, and combine them into a demixing matrix W = [w1, w2, …, w n , and perform amplitude recovery using the least - squares method according to the observed value to obtain the estimated value of the signal source , each column vector of

[0099] is the estimated value of the signal source; Figure 4 As shown in

[0100] Step 203: For the signals separated by ICA, perform signal screening and reconstruction, specifically: k Step 2031: Evaluate the IMFs (IMF is the abbreviation of Intrinsic Mode Function, and n is the k - th IMF (one of IMF1, IMF2, …, to IMF n in n ), and IMFs refers to all IMFs, that is, Intrinsic Mode Functions) from the aspects of morphology and time - frequency characteristics respectively, identify and suppress the noise components, and retain the useful signal components;

[0101] Step 2032: Re - combine the useful signal components retained in Step 2031 to construct the denoised magnetotelluric sounding signal;

[0102] Step 2033: Evaluate the performance of the denoised signal obtained in Step 2032, and use the normalized cross - correlation coefficient (NCC), signal - to - noise ratio (SNR), and reconstruction error (E) as evaluation indicators to ensure that the denoising effect meets the expectations;

[0103] Among them, the evaluation indicators are respectively:

[0104] ;

[0105] ;

[0106] ;

[0107] In the formula, is the original data, is the data after denoising processing. The higher the signal-to-noise ratio, the closer the NCC is to 1, and the smaller the reconstruction error E is, indicating that the denoising effect is better.

[0108] In step 300, repeat the above operations to perform denoising processing on the remaining magnetotelluric sounding signals obtained. Specifically:

[0109] Repeat the above operations to perform denoising processing on the remaining 4 channels of magnetotelluric sounding signals.

[0110] For Figure 2 , Figure 3 , Figure 4 in the content is explained. Among them, is the original data of each channel of magnetotelluric sounding, i is the number of channels, d is the current iteration number, D is the preset maximum iteration number, K is the set number of modes, k is the IMF order of the current decomposition, is the k-th order IMF component of the current signal, is the center frequency of this component, is the regularization parameter, SVD is the singular value decomposition operation, is the separation matrix, and each column vector in it is the extraction vector, is the extraction vector in the current iteration, is the updated extraction vector, is the set tolerance error parameter, is the estimated value of the signal matrix of the i-th channel, and each column is respectively each signal source estimated by this channel, is the signal estimated value of the current channel.

[0111] In this specification, each embodiment is described in a progressive manner. The key point of each embodiment is to illustrate the differences from other embodiments. For the same and similar parts between each embodiment, reference can be made to each other.

[0112] In this article, specific examples are used to elaborate on the principle and implementation manner of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea; at the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manner and application scope. In summary, the content of this specification should not be construed as a limitation to the present invention.

Claims

1. A denoising method for non-stationary and non-linear magnetotelluric sounding signals, characterized in that, Including: Obtaining magnetotelluric sounding signals from a magnetotelluric sounding instrument; Selecting one channel from the obtained magnetotelluric sounding signals and denoising it, specifically: Perform VMD decomposition on magnetotelluric signals, and define the intrinsic mode function as an amplitude-frequency modulated signal as the IMF to be optimized k component, that is, the k-th order intrinsic mode component; Solving a constrained variational optimization problem for intrinsic mode functions; Combined with Parseval's theorem and Fourier isometric transformation, it is transformed into the frequency domain. The alternating direction multiplier method is used for iterative optimization in Matlab to calculate the frequency domain expressions of each mode in the frequency domain, and the time domain expressions of the IMFs with different center frequencies are obtained based on the inverse Fourier transform k The time domain expressions of the components; Constructing an observation matrix for ICA based on the obtained intrinsic mode components of each order; Performing pre-whitening and centering operations on the observation matrix; Constructing an optimization function for the weight vector based on negentropy and kurtosis, iterating the weight vector using the gradient ascent method, repeating the iterative operation, obtaining n weight vectors, forming a demixing matrix, and calculating the estimated value of the signal source; Evaluating the IMFs components after ICA processing respectively from the aspects of morphology and time-frequency characteristics, identifying and suppressing noise components, retaining useful signal components, recombining the useful signal components, constructing a denoised magnetotelluric sounding signal, and evaluating its performance; Repeating the above operations to perform denoising processing on the remaining magnetotelluric sounding signals obtained; 2. The method according to claim 1, wherein Obtaining magnetotelluric sounding signals from a magnetotelluric sounding instrument, specifically: Obtaining 5 magnetotelluric sounding signals from a magnetotelluric sounding instrument, where the first and second channels are the electrical channel measurement values in the north-south and east-west directions respectively, the third and fourth channels are the magnetic channel measurement values in the north-south and east-west directions respectively, and the fifth channel is the vertical magnetic channel measurement value.

3. The method according to claim 2, wherein Perform VMD decomposition on magnetotelluric signals, and define the intrinsic mode function as an amplitude-frequency modulated signal as the IMF to be optimized k component, that is, the k-th order intrinsic mode component, specifically: Define the Intrinsic Mode Function The amplitude-modulated and frequency-modulated signal is used as the IMF to be optimized k component, that is, the k-th order intrinsic mode component, where and are the instantaneous amplitude and instantaneous phase of the IMF k component respectively, and the instantaneous frequency is obtained by differentiating the phase, which is: ； Use to fit the input signal.

4. The method according to claim 3, wherein Solving a constrained variational optimization problem for intrinsic mode functions, specifically: For Solving the constrained variational optimization problem to ensure that each mode satisfies certain constraint conditions, which are: ； In the formula, is the k-th order modal component, is the instantaneous frequency, is an input magnetotelluric sounding signal, is the summation operation for k = 1, 2, …, K, is the operation of taking the partial derivative with respect to time, is the convolution operation, is the square of the second norm, is the operation of finding the minimum value, is the constraint condition, and δ(t) is the impulse function.

5. The method according to claim 4, wherein Combined with Parseval's theorem and Fourier isometric transformation, it is transformed into the frequency domain, and the alternating direction multiplier method is used for iterative optimization in Matlab to calculate the frequency domain expressions of each mode in the frequency domain, and the time domain expressions of the IMF components with different center frequencies are obtained based on the inverse Fourier transform. Specifically: k The time domain expression of the component is as follows: Combined with Parseval's theorem and Fourier isometric transformation, it is transformed into the frequency domain, and the alternating direction multiplier method is used for iterative optimization in Matlab, so that their energy in the frequency domain is mainly concentrated in the predetermined frequency band, and the sum of each mode obtained in the frequency domain is as follows. The expression after alternating optimization iteration is: ； ； After Fourier transform, they are all functions of frequency ω at this time. The original x(t) and y(t) are transformed into x(ω) and y(ω). The subscript k represents the k-th IMF component, i represents that the current iteration number is the i-th, the superscript n is the total number of iterations, λ is the Lagrange multiplier used in the iterative calculation, α is the penalty factor, ω k is the central frequency of the k-th IMF; The time-domain expressions of IMFs with different center frequencies are obtained by taking the inverse Fourier transform and then taking the real part. Among them, the expression of the inverse Fourier transform is as follows: k ​ ； Selecting and updating modes according to the set number of modes K, and ending the iterative process when the threshold of the number of iterations or residual energy is met, that is: ； wherein, is the threshold value, and finally K IMF components are obtained.

6. The method according to claim 3, characterized in that, Constructing an observation matrix for ICA based on the obtained intrinsic mode components of each order, specifically: Obtain the eigenmode components of each order in the time-frequency domain after decomposition, and construct the observation matrix X of ICA as X = [IMF1, IMF2, …, IMF n , where n independent source signals are denoted as S = [s1, s2, …, s n , and under the action of the mixing matrix A = [a1, a2, …, a n , the observations X = [x1, x2, …, x n are obtained, which is represented by the following matrix equation: 。 7. The method according to claim 6, wherein Constructing an optimization function for the weight vector based on negentropy and kurtosis, iterating the weight vector using the gradient ascent method, repeating the iterative operation, obtaining n weight vectors, forming a demixing matrix, and calculating the estimated value of the signal source, specifically: The optimization function for constructing the weight vector based on negentropy and kurtosis is: ； In the formula, represents kurtosis, defined as , represents negative entropy, approximated as , where is a non - quadratic function, is the standard normal distribution, is the operation of taking the mathematical expectation; Iterating the weight vector using the gradient ascent method; Perform iterative operations to obtain n weight vectors and combine them into a demixing matrix W = [w1, w2, …, w n , and perform amplitude recovery using the least squares method based on the observed value to obtain the estimated value of the signal source , where n is the number of estimated source signals, and each column vector of is the estimated value of the signal source.

8. The method according to claim 7, wherein Evaluating the IMFs components after ICA processing respectively from the aspects of morphology and time-frequency characteristics, identifying and suppressing noise components, retaining useful signal components, recombining the useful signal components, constructing a denoised magnetotelluric sounding signal, and evaluating its performance, specifically: Evaluating the IMFs after ICA processing respectively from the aspects of morphology and time-frequency characteristics, identifying and suppressing noise components, retaining useful signal components, recombining the useful signal components, constructing a denoised magnetotelluric sounding signal, and using the normalized cross-correlation coefficient, signal-to-noise ratio, and reconstruction error as evaluation indicators to evaluate its performance.

Images