Denoising method for non-stationary nonlinear magnetotelluric sounding signal

Abstract

The invention provides a denoising method for non-stationary nonlinear magnetotelluric sounding signals, and relates to the technical field of geophysical exploration, and the method comprises the steps: obtaining magnetotelluric sounding signals from a magnetotelluric depth finder, selecting one from the obtained magnetotelluric sounding signals, carrying out the denoising of the selected magnetotelluric sounding signals, repeating the above operation, and carrying out the denoising of the selected magnetotelluric sounding signals. And carrying out de-noising processing on the other acquired magnetotelluric sounding signals. VMD (variational mode decomposition) and ICA (independent component analysis) are combined, denoising of the magnetotelluric sounding signals is realized from multiple angles of time-frequency domain and statistical property, and when effective signals are interfered by long-time high-energy non-stationary noise, the method can ensure high quality of reconstructed signals, effectively reduce the influence of manual parameter selection, and improve the denoising performance of the magnetotelluric sounding signals. Compared with the prior art, the method has the advantages that the method is high in robustness to noise, clear physical interpretation is provided, limitation of a traditional denoising method is overcome, and a more reliable and accurate data processing means is provided for the field of geophysical exploration.

Description

Technical Field

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

[0002] As an effective geophysical exploration technology, magnetotelluric (MT) sounding has always been highly valued by scientific researchers and has been widely used in the fields of crust and upper mantle structure research, oil and gas mineral resource exploration, and earthquake monitoring. The MT method records the temporal changes of the natural electromagnetic field on the earth's surface and uses the principle of electromagnetic induction to deduce the conductivity distribution of underground media. Although the MT method has deep penetration and high resolution, its signal is easily affected by natural and artificial noise during the acquisition process, which seriously reduces the signal-to-noise ratio of MT data and may cause distortion of apparent resistivity and phase curves, affecting the inversion and interpretation of magnetotelluric data. Common human noise such as high-voltage power lines and radio transmission towers will cause long-term strong interference to the observation of magnetotelluric sounding signals. Therefore, studying effective denoising technology is of great significance to improving the quality of MT data and enhancing the accuracy of interpretation.

[0003] Since Tikhonov and Cagniard proposed the magnetotelluric method in the 1950s, research on noise suppression has made significant progress and derived a variety of methods. Depending on the processing object, magnetotelluric denoising technology 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. These methods are all implemented in the impedance estimation process in the frequency domain, which can respectively remove non-correlated noise and suppress the influence of extreme values ​​of "flying points"; 2. Time domain analysis methods such as time domain impedance estimation based on Bayesian estimation, which updates the estimate of parameters by combining prior knowledge and observation data, is a framework for statistical inference under uncertainty. Compared with the frequency domain based on Fourier transform, Method, which has more advantages in processing non-stationary data; 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 than the traditional bandpass filter; in addition, machine learning methods such as artificial neural networks, through their powerful pattern recognition and feature extraction capabilities, can capture the local characteristics of signals from massive data, automatically learn and adjust model parameters to accurately denoise the signal from the time domain; 3. Time-frequency domain analysis methods include wavelet analysis and Hilbert-Huang Transform (HHT). Wavelet analysis decomposes the signal into components of different time scales by scaling and translating the wavelet basis function. The Hilbert-Huang transform adaptively obtains the intrinsic mode function through empirical mode decomposition and obtains the time-frequency spectrum by Hilbert transform. Both methods can obtain time-frequency characteristics. When estimating impedance, combined with remote reference data and robust estimation, it can effectively suppress coherent noise while maintaining a high time-frequency resolution. Therefore, it has received extensive attention and research from relevant researchers.

[0004] Although the above methods are feasible in theory, when processing field data of magnetotelluric sounding, the non-stationary nonlinearity and non-minimum phase characteristics of the signal make the traditional spectrum analysis based on Fourier transform have certain limitations. In the time domain denoising method, 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 result. In addition, in practical applications, it may be difficult to select a suitable Bayesian model due to the complexity of the magnetotelluric system and incomplete data information; the effect of mathematical morphology filtering strongly depends on the selection of structural elements, which requires professional knowledge and experience, and different structural elements may be required for different signals; machine learning methods usually require a lot of computing resources to train deep learning models, including high-performance GPUs and a lot of memory. In addition, its internal working mechanism and decision-making process are difficult to explain, which may be a problem in magnetotelluric data denoising that requires a clear physical explanation. In the process of time-frequency analysis, the underground structure is complex and changeable, the effect of wavelet decomposition depends on the selection of wavelet basis in the current environment, and the threshold setting has a significant impact on the denoising result. Since the time-frequency domain analysis method is particularly suitable for processing non-stationary signals, it can obtain detailed information at different time scales and frequency levels, and has unique advantages in identifying local features in the signal, such as instantaneous events or transient phenomena, making EMD and its improved methods the mainstream methods in time-frequency analysis and widely used. However, empirical mode decomposition lacks extensive standardization and theoretical basis, which to some extent limits its accuracy and reliability in the decomposition process. In particular, when dealing with modal aliasing problems, EMD may encounter difficulties. Different signal modes may interfere with each other, resulting in the inability to clearly distinguish modal components, and the decomposition results may not always be directly related to physical phenomena.

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

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

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

[0008] The present invention provides a denoising method for non-stationary nonlinear magnetotelluric sounding signals, comprising:

[0009] Acquire magnetotelluric sounding signals from magnetotelluric sounding instruments;

[0010] Select one of the acquired magnetotelluric sounding signals and denoise it;

[0011] Repeat the above operation to perform denoising on the remaining magnetotelluric sounding signals.

[0012] Preferably, the magnetotelluric sounding signal is obtained from the magnetotelluric sounding instrument, specifically:

[0013] Five magnetotelluric sounding signals were obtained from the magnetotelluric sounding instrument, among which the first and second tracks were the north-south and east-west electric track measurement values, the third and fourth tracks were the north-south and east-west magnetic track measurement values, and the fifth track was the vertical magnetic track measurement value.

[0014] Preferably, one of the acquired magnetotelluric sounding signals is selected and denoised, specifically:

[0015] Perform VMD decomposition on the magnetotelluric signal and define the intrinsic mode function as the amplitude-frequency modulation signal as the IMF to be optimized. k Component, i.e., the k-th order eigenmode component;

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

[0017] Combining Parseval's theorem and Fourier isometric transform, we convert the frequency domain into the frequency domain, use the alternating direction multiplier method in MATLAB for iterative optimization, calculate the frequency domain expression of each mode in the frequency domain, and obtain IMF with different center frequencies based on the inverse Fourier transform. k Time domain expressions of components;

[0018] Based on the obtained eigenmode components of each order, the ICA measurement matrix is ​​constructed;

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

[0020] An optimization function of a weight vector is constructed based on negative entropy and kurtosis, and the weight vector is iterated using a gradient ascent method. The iterative operation is repeated to obtain n weight vectors, form an unmixing matrix, and calculate the estimated value of the signal source;

[0021] The IMFs components after ICA processing are evaluated from the perspective of morphology and time-frequency characteristics, the noise components are identified and suppressed, the useful signal components are retained, the useful signal components are recombined, the denoised magnetotelluric sounding signal is constructed, and its performance is evaluated.

[0022] Preferably, the magnetotelluric signal is subjected to VMD decomposition processing, and the intrinsic mode function is defined as an amplitude-frequency modulated signal as the IMF to be optimized. k Component, that is, the k-th order eigenmode component, is specifically:

[0023] Define the eigenmode function For the AM / FM signal as the IMF to be optimizedk Component, that is, the k-th order eigenmode component, where and IMF k The instantaneous amplitude and instantaneous phase of the component, and the instantaneous frequency are obtained by differentiating the phase, which is:

[0024] ;

[0025] use Fit the input signal.

[0026] Preferably, a constrained variational optimization problem is constructed for the intrinsic mode function to be solved, specifically:

[0027] right Construct a constrained variational optimization problem to solve, ensuring that each mode satisfies certain constraints, which is:

[0028] ;

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

[0030] Preferably, Parseval theorem and Fourier isometric transform are combined to convert to the frequency domain, and the alternating direction multiplier method is used in MATLAB for iterative optimization to calculate the frequency domain expression of each mode in the frequency domain, and the IMF with different center frequencies is obtained based on the inverse Fourier transform. k The time domain expression of the component is:

[0031] Combining Parseval's theorem and Fourier isometric transform, they are converted to the frequency domain, and the alternating direction multiplier method is used in MATLAB for iterative optimization, so that their energy in the frequency domain is mainly concentrated in the predetermined frequency band. and , where the expression after alternating optimization iteration is:

[0032] ;

[0033] ;

[0034] After Fourier transform, they are now functions of frequency ω. The original x(t) and y(t) are transformed into x(ω) and y(ω). The subscript k represents the kth IMF component, i represents the current iteration number is the ith, the superscript n is the total number of iterations, λ is the Lagrange multiplier, used in iterative calculations, α is the penalty factor, and ω k is the center frequency of the kth IMF;

[0035] By taking the inverse Fourier transform and taking the real part, we can get IMFs with different center frequencies. k The time domain expression of , where the expression of the inverse Fourier transform is:

[0036] ;

[0037] The mode is selected and updated according to the set mode number K. When the threshold of the number of iterations or residual energy is met, the iteration 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, an ICA measurement matrix is ​​constructed, specifically:

[0041] Obtain the decomposed intrinsic mode components of each order in the time-frequency domain and construct the ICA measurement matrix X = [IMF 1 ,IMF 2 ,…,IMF n ], where n independent source signals are represented as S=[s 1 ,s 2 ,…,s n ], in the mixing matrix A=[a 1 ,a 2 ,…,a n ], we get the observed value X=[x 1 ,x 2 ,…,x n ], which can be expressed as the following matrix equation:

[0042] .

[0043] Preferably, an optimization function of a weight vector is constructed based on negative entropy and kurtosis, and the weight vector is iterated using a gradient ascent method. The iterative operation is repeated to obtain n weight vectors, form an unmixing matrix, and calculate the estimated value of the signal source, specifically:

[0044] The optimization function for constructing the weight vector based on negative entropy and kurtosis is:

[0045] ;

[0046] In the formula, represents the kurtosis, defined as , represents negative entropy, which is approximately ,in, is a non-quadratic function, is the standard normal distribution, Operations to find mathematical expectations;

[0047] The weight vector is iterated using the gradient ascent method;

[0048] Repeat the iterative operation to obtain n weight vectors, which are combined into the unmixing matrix W=[w 1 ,w 2 ,…,w n ], according to the observed The value of is restored by the least square method to obtain the signal source Estimated value of , where n is the number of estimated source signals, Each column vector of is the estimated value of the signal source.

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

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

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

[0052] The present invention provides a denoising method for non-stationary nonlinear magnetotelluric sounding signals, the method comprising acquiring magnetotelluric sounding signals from a magnetotelluric sounding instrument, selecting one of the acquired magnetotelluric sounding signals, denoising it, repeating the above operation, and denoising the remaining acquired magnetotelluric sounding signals. The present invention combines variational mode decomposition (VMD) and independent component analysis (ICA) to achieve denoising of magnetotelluric sounding signals from multiple angles of time-frequency domain and statistical characteristics. When the effective signal is interfered by long-term high-energy non-stationary noise, the method can ensure the high quality of the reconstructed signal, effectively reduce the influence of artificial parameter selection, be robust to noise, and provide a clear physical explanation at the same time, overcome the limitations of traditional denoising methods, and provide a more reliable and accurate data processing method for the field of geophysical exploration. The advantages of the present invention include:

[0053] 1. Through variational mode decomposition technology, 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. The modal aliasing problem is further solved through independent component analysis, which shows significant effectiveness in processing signal-noise mixing in the same frequency band. For long-standing and powerful interference noise, it can not only cope with the instantaneous changes of the signal, but also stably identify and suppress interference over a long time scale, ensuring the quality and reliability of magnetotelluric sounding data denoising;

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

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

[0057] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0058] Figure 1A schematic flow chart of a denoising method for non-stationary nonlinear magnetotelluric sounding signals provided in an embodiment of the present invention;

[0059] Figure 2 It is a schematic diagram of VMD decomposition process;

[0060] Figure 3 This is a schematic diagram of the ICA separation process;

[0061] Figure 4 Schematic diagram of signal screening and reconstruction process. DETAILED DESCRIPTION

[0062] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0063] The present invention combines variational mode decomposition (VMD) and independent component analysis (ICA) to achieve 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, the method can ensure the high quality of the reconstructed signal, effectively reduce the influence of artificial parameter selection, be robust to noise, and provide a clear physical explanation at the same time, thus overcoming the limitations of traditional denoising methods and providing a more reliable and accurate data processing method for the field of geophysical exploration.

[0064] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0065] Figure 1 A flow chart of a method provided by an embodiment of the present invention, such as Figure 1 As shown, the present invention provides a denoising method for non-stationary nonlinear magnetotelluric sounding signals, comprising:

[0066] Step 100: Acquire a magnetotelluric sounding signal from a magnetotelluric sounding instrument;

[0067] Step 200: Select one of the acquired magnetotelluric sounding signals and perform denoising on it;

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

[0069] In step 100, a magnetotelluric sounding signal is obtained from a magnetotelluric sounding instrument, specifically:

[0070] Five magnetotelluric sounding signals were obtained from the magnetotelluric sounding instrument, among which the first and second tracks were the north-south and east-west electric track measurement values, the third and fourth tracks were the north-south and east-west magnetic track measurement values, and the fifth track was the vertical magnetic track measurement value.

[0071] In step 200, one of the acquired magnetotelluric sounding signals is selected and denoised, specifically:

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

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

[0074] Step 203: performing signal screening and reconstruction on the signal separated by ICA;

[0075] like Figure 2 As shown, in step 201, the magnetotelluric sounding signal is subjected to VMD decomposition processing, specifically:

[0076] Step 2011: Define the intrinsic mode function For the AM / FM signal as the IMF to be optimized k Component, that is, the k-th order eigenmode component, where and IMF k The instantaneous amplitude and instantaneous phase of the wave, and the instantaneous frequency can be obtained by taking the derivative of the phase ),use Fitting the input signal;

[0077] Step 2012: For step 2011 Construct a constrained variational optimization problem to solve, ensuring that each mode meets certain constraints, such as having higher energy near the center frequency of the frequency band:

[0078] ;

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

[0080] Step 2013: Combine Parseval's theorem and Fourier isometric transform to convert to the frequency domain, and use the alternating direction multiplier method in MATLAB to perform iterative optimization so that their energy in the frequency domain is mainly concentrated in the predetermined frequency band. and , by taking the inverse Fourier transform and taking the real part, we can get IMFs with different center frequencies k The time domain expression of , solves the optimization problem in step 4;

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

[0082] ;

[0083] ;

[0084] After Fourier transform, they are now functions of frequency ω. The original x(t) and y(t) are transformed into x(ω) and y(ω). The subscript k represents the kth IMF component; i represents the current iteration number is the ith; the superscript n is the total number of iterations; λ is the Lagrange multiplier, which is used in iterative calculations; α is the penalty factor; ω k is the center frequency of the kth IMF;

[0085] The expression of inverse Fourier transform is:

[0086] ;

[0087] The mode is selected and updated according to the set mode number K. When the threshold of the number of iterations or residual energy is met, the iteration process ends, that is:

[0088] ;

[0089] Finally, K IMF components are obtained.

[0090] like Figure 3 As shown, in step 202, ICA separation is performed on the signal after VMD decomposition processing, specifically:

[0091] Step 221: Take the intrinsic mode components of each order obtained by the signal decomposition in step 2013, that is, K IMF components, and construct the ICA observation matrix X = [IMF 1 ,IMF 2 ,…,IMF n ], let n independent source signals be expressed as S=[s 1 ,s 2 ,…,s n ], in the mixing matrix A=[a 1 ,a2 ,…,a n ], we get the observed value X=[x 1 ,x 2 ,…,x n ] can be expressed by the following matrix equation:

[0092] .

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

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

[0095] ;

[0096] In the formula, represents the kurtosis, defined as , represents negative entropy, which can be approximated as ,in, is a non-quadratic function, is the standard normal distribution, Operations to find mathematical expectations;

[0097] The weight vector is iterated using the gradient ascent method;

[0098] Step 2024: Repeat the iterative operation of step 2023 to obtain n weight vectors, where n is the number of estimated source signals, and combine them into an unmixing matrix W = [w 1 ,w 2 ,…,w n ], according to the observed The value of is restored by the least square method to obtain the signal source Estimated value of , Each column vector of is the estimated value of the signal source;

[0099] like Figure 4 As shown, in step 203, the signal after ICA separation is screened and reconstructed, specifically:

[0100] Step 2031: IMFs (IMF is the abbreviation of Intrinsic Mode Function, IMF) processed by ICA are analyzed from the perspective of morphology and time-frequency characteristics. k is the kth IMF (from IMF 1, IMF 2 , …, to the IMF n One of them), IMFs refers to all IMFs, namely Intrinsic Mode Functions) to evaluate, identify and suppress noise components, and retain useful signal components;

[0101] Step 2032: Recombining the useful signal components retained in step 2031 to construct a denoised magnetotelluric sounding signal;

[0102] Step 2033: Perform performance evaluation on the denoised signal obtained in step 2032, using the normalized cross correlation coefficient (NCC), signal-to-noise ratio (SNR) and reconstruction error (E) as evaluation indicators to ensure that the denoising effect reaches the expected result;

[0103] The evaluation indicators are:

[0104] ;

[0105] ;

[0106] ;

[0107] In the formula, is the original data, is the data after denoising. 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, the above operation is repeated to perform denoising on the remaining magnetotelluric sounding signals obtained, specifically:

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

[0110] right Figure 2 , Figure 3 , Figure 4 The contents in are described, among which, is the original data of each channel of magnetotelluric sounding, i is the channel number, d is the current iteration number, D is the preset maximum iteration number, K is the set mode number, k is the current decomposed IMF order, 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 a separation matrix, in which each column vector To extract the 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, where each column are each signal source estimated by the channel, It is the signal estimation value of the current channel.

[0111] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0112] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A denoising method for non-stationary nonlinear magnetotelluric sounding signals, characterized in that: include: Acquire magnetotelluric sounding signals from magnetotelluric sounding instruments; Select one of the acquired magnetotelluric sounding signals and perform denoising on it, specifically: Perform VMD decomposition on the magnetotelluric signal and define the intrinsic mode function as the amplitude-frequency modulation signal as the IMF to be optimized. k Component, i.e., the k-th order eigenmode component; Construct a constrained variational optimization problem for the intrinsic mode function and solve it; Combining Parseval's theorem and Fourier isometric transform, we convert the frequency domain into the frequency domain, use the alternating direction multiplier method in MATLAB for iterative optimization, calculate the frequency domain expression of each mode in the frequency domain, and obtain IMF with different center frequencies based on the inverse Fourier transform. k Time domain expressions of components; Based on the obtained eigenmode components of each order, the ICA measurement matrix is ​​constructed; Perform pre-whitening and centering operations on the observation matrix; An optimization function of a weight vector is constructed based on negative entropy and kurtosis, and the weight vector is iterated using a gradient ascent method. The iterative operation is repeated to obtain n weight vectors, form an unmixing matrix, and calculate the estimated value of the signal source; The IMFs components after ICA processing are evaluated from the perspective of morphology and time-frequency characteristics, the noise components are identified and suppressed, the useful signal components are retained, the useful signal components are recombined, the denoised magnetotelluric sounding signal is constructed, and its performance is evaluated; Repeat the above operation to perform denoising on the remaining magnetotelluric sounding signals.

2. The method according to claim 1, characterized in that The magnetotelluric sounding signal is obtained from the magnetotelluric sounding instrument, specifically: Five magnetotelluric sounding signals were obtained from the magnetotelluric sounding instrument, among which the first and second tracks were the north-south and east-west electric track measurement values, the third and fourth tracks were the north-south and east-west magnetic track measurement values, and the fifth track was the vertical magnetic track measurement value.

3. The method according to claim 2, characterized in that Perform VMD decomposition on the magnetotelluric signal and define the intrinsic mode function as the amplitude-frequency modulation signal as the IMF to be optimized. k Component, that is, the k-th order eigenmode component, is specifically: Define the eigenmode function For the AM / FM signal as the IMF to be optimized k Component, that is, the k-th order eigenmode component, where and IMF k The instantaneous amplitude and instantaneous phase of the component, and the instantaneous frequency are obtained by differentiating the phase, which is: ； use Fit the input signal.

4. The method according to claim 3, characterized in that The constrained variational optimization problem is constructed for the intrinsic mode function and solved as follows: right Construct a constrained variational optimization problem to solve, ensuring that each mode satisfies certain constraints, which is: ； In the formula, is the kth order modal component, is the instantaneous frequency, is an input magnetotelluric sounding signal, To perform the sum operation for k=1,2,…,K, To find the partial derivative with respect to time, is the convolution operation, To find the square of the two-norm, To find the minimum value, is the constraint condition, and δ(t) is the impulse function.

5. The method according to claim 4, characterized in that Combining Parseval's theorem and Fourier isometric transform, we convert the frequency domain into the frequency domain, use the alternating direction multiplier method in MATLAB for iterative optimization, calculate the frequency domain expression of each mode in the frequency domain, and obtain IMF with different center frequencies based on the inverse Fourier transform. k The time domain expression of the component is: Combining Parseval's theorem and Fourier isometric transform, they are converted to the frequency domain, and the alternating direction multiplier method is used in MATLAB for iterative optimization, so that their energy in the frequency domain is mainly concentrated in the predetermined frequency band. and , where the expression after alternating optimization iteration is: ； ； After Fourier transform, they are now functions of frequency ω. The original x(t) and y(t) are transformed into x(ω) and y(ω). The subscript k represents the kth IMF component, i represents the current iteration number is the ith, the superscript n is the total number of iterations, λ is the Lagrange multiplier, used in iterative calculations, α is the penalty factor, and ω k is the center frequency of the kth IMF; By taking the inverse Fourier transform and taking the real part, we can get IMFs with different center frequencies. k The time domain expression of , where the expression of the inverse Fourier transform is: ； The mode is selected and updated according to the set mode number K. When the threshold of the number of iterations or residual energy is met, the iteration process ends, that is: ； In the formula, is the threshold, and finally K IMF components are obtained.

6. The method according to claim 3, characterized in that: Based on the obtained eigenmode components of each order, the ICA measurement matrix is ​​constructed, specifically: Obtain the decomposed intrinsic mode components of each order in the time-frequency domain and construct the ICA measurement matrix X=[IMF1,IMF2,…,IMF n ], where n independent source signals are represented by S=[s1,s2,…,s n ], in the mixing matrix A=[a1,a2,…,a n ], we get the observed value X=[x1,x2,…,x n ], which can be expressed as the following matrix equation: 。 7. The method according to claim 6, characterized in that The optimization function of the weight vector is constructed based on negative entropy and kurtosis. The weight vector is iterated using the gradient ascent method. The iterative operation is repeated to obtain n weight vectors to form an unmixing matrix and calculate the estimated value of the signal source. Specifically: The optimization function for constructing the weight vector based on negative entropy and kurtosis is: ； In the formula, represents the kurtosis, defined as , represents negative entropy, which is approximately ,in, is a non-quadratic function, is the standard normal distribution, Operations to find mathematical expectations; The weight vector is iterated using the gradient ascent method; Repeat the iterative operation to obtain n weight vectors and combine them into an unmixing matrix W=[w1,w2,…,w n ], according to the observed The value of is restored by the least square method to obtain the signal source Estimated value of , where n is the number of estimated source signals, Each column vector of is the estimated value of the signal source.

8. The method according to claim 7, characterized in that The IMFs components after ICA processing are evaluated from the perspective of morphology and time-frequency characteristics, the noise components are identified and suppressed, the useful signal components are retained, the useful signal components are recombined, the denoised magnetotelluric sounding signal is constructed, and its performance is evaluated, specifically: The IMFs after ICA processing are evaluated from the perspective of morphology and time-frequency characteristics, the noise components are identified and suppressed, the useful signal components are retained, the useful signal components are recombined, and the denoised magnetotelluric sounding signals are constructed. The normalized cross-correlation coefficient, signal-to-noise ratio and reconstruction error are used as evaluation indicators to evaluate its performance.

Images