A method, device and medium for reducing noise of multi-level magnetotelluric data

Abstract

The embodiments of the present specification disclose a data denoising method, device and medium for multi-level magnetotelluric data, which relate to the field of noise reduction processing technology. The method includes: acquiring an original magnetotelluric data sequence, determining first noise reduction data, performing variational mode decomposition on the first noise reduction data, and determining second noise reduction data; determining a first signal-to-noise ratio corresponding to the second noise reduction data, and when the first signal-to-noise ratio is less than a preset signal-to-noise ratio reference value, dynamically fusing residual signals to determine a fused residual signal; performing noise reduction processing on the fused residual signal to determine fused residual noise reduction data, performing residual recovery based on the fused residual noise reduction data, and determining current noise reduction data after residual recovery; and outputting the current noise reduction data when the second signal-to-noise ratio corresponding to the current noise reduction data is not less than the signal-to-noise ratio reference value. This method can meet the signal recovery requirements of weak effective information in magnetotelluric data, enhance the fidelity of the signal, and improve the capabilities of deep noise processing and weak signal protection.

Description

Technical Field

[0001] The present invention relates to the field of noise reduction processing technology, and in particular to a method, device and medium for data noise reduction of multi-level magnetotelluric data. Background Art

[0002] Magnetotelluric exploration, an important means of detecting deep geological structures, collects natural electromagnetic signals with wide bandwidth, weak energy, and nonlinear and nonstationary characteristics. Frequency-domain noise reduction technology relies on the spectral separation characteristics of signal and noise. However, under low signal-to-noise ratio conditions, the effective spectrum is often completely submerged by noise, resulting in distortion of many frequency points and high phase errors in the measured apparent resistivity-phase curve. While time-domain noise reduction methods can circumvent the problem of spectral aliasing, they require a pre-defined empirical noise model (such as the assumption of a fixed power frequency period), have poor adaptability to non-periodic pulse interference (square wave, triangle wave, and sudden noise). Furthermore, uniformly filtering out abnormal waveforms results in effective signal loss, severely weakening the deep geological response.

[0003] To overcome traditional limitations, an existing technique (patent CN 119046617A) proposes a combined denoising method for industrial pumps based on SVD and GWO-VMD. This method uses multiresolution singular value decomposition (MSVD) for initial denoising, then optimizes variational mode decomposition (VMD) parameters using the Grey Wolf algorithm (GWO). Finally, modal components are selected based on the kurtosis-correlation coefficient. However, this existing technique still has significant limitations in denoising multi-level magnetotelluric (MT) data. Industrial pump signals differ significantly from MT signals: industrial pump signals are periodic, while MT data are impulsive and harmonic. This requires protecting weak signals and addressing ultra-wideband dynamic noise in the 0.01-1000 Hz range. The residual signal after denoising may contain valid signals, especially for weak MT signals. Directly discarding the residual signal can result in loss of valid signals. Moreover, in the existing multi-level magnetotelluric data noise reduction process, the noise reduction method is fixed, the purpose of MT signals is ignored, and the different impacts on MT signals under different geological analysis requirements are not considered.

[0004] Therefore, the signal characteristics of the MT signal and the geological analysis requirements are ignored during the MT signal denoising process. The residual signal is discarded after denoising, resulting in effective signal loss, and the denoising effect cannot meet the geological analysis requirements. Summary of the Invention

[0005] One or more embodiments of this specification provide a data denoising method, device, and medium for multi-layer magnetotelluric data, which are used to solve the following technical problems: ignoring the signal characteristics and geological analysis requirements of the MT signal during the denoising process, discarding the residual signal after denoising, resulting in loss of effective signal, and the denoising effect cannot meet the geological analysis requirements.

[0006] One or more embodiments of this specification adopt the following technical solutions:

[0007] One or more embodiments of the present specification provide a data denoising method for multi-level magnetotelluric data, the method comprising: acquiring an original magnetotelluric data sequence, performing resolution singular value decomposition on the original magnetotelluric data sequence to determine first denoised data, performing variational mode decomposition on the first denoised data using a predetermined set of optimization parameters to determine second denoised data; determining a first signal-to-noise ratio corresponding to the second denoised data, and when the first signal-to-noise ratio is less than a preset signal-to-noise ratio reference value, dynamically fusing a residual signal based on the first denoised data and the second denoised data to determine a fused residual signal; performing denoising on the fused residual signal to determine corresponding fused residual denoised data, performing residual recovery based on the fused residual denoised data, and determining current denoised data after residual recovery; and outputting the current denoised data when the second signal-to-noise ratio corresponding to the current denoised data is not less than the signal-to-noise ratio reference value.

[0008] One or more embodiments of this specification provide a data denoising device for multi-level magnetotelluric data, including:

[0009] at least one processor; and,

[0010] a memory communicatively connected to the at least one processor; wherein,

[0011] The memory stores instructions that can be executed by the at least one processor. The instructions are executed by the at least one processor to enable the at least one processor to perform the above method.

[0012] One or more embodiments of this specification provide a non-volatile computer storage medium storing computer-executable instructions, wherein the computer-executable instructions are configured to execute the above method.

[0013] At least one of the above technical solutions adopted in the embodiments of this specification can achieve the following beneficial effects: through the technical solutions of the embodiments of this specification, the original magnetotelluric data sequence corresponding to the target geological analysis area in the target geological analysis requirement information is collected, ensuring the regional matching of the signal source and the analysis requirement; the MSVD algorithm is introduced to perform preliminary noise reduction on the non-stationary signal MT signal using multi-level decomposition. Compared with the SVD of single-level linear decomposition, the low-frequency noise in the MT signal can be processed at multiple scales and high resolution, and the target signal can be effectively retained while suppressing the noise; in response to problems such as modal aliasing generated in the decomposition process of the residual pulse signal, the Thomson function is introduced on the basis of the VMD algorithm, and the components exceeding the threshold are adaptively soft-thresholded to avoid hard truncation distortion, which is more effective than the traditional VMD algorithm. The synthesis of pulse signals is effectively suppressed; the high convergence characteristic of the GWO algorithm is used to optimize the key parameters in VMD. During the optimization process, the characteristic vector library of electromagnetic interference in the target geological analysis area is introduced as noise prior information. By constructing a composite fitness function that integrates envelope entropy, decomposition results and noise feature matching, target-guided adjustment of VMD parameters is achieved, ensuring that the noise prior knowledge matches the target geological analysis area. This not only greatly reduces the possibility of modal aliasing, but also improves computational efficiency and significantly shortens calculation time by optimizing the number of modes; the dynamic residual feedback method and weighted fusion mechanism are used to further process the effective signals in the residuals, which can meet the signal recovery requirements of weak effective information in magnetotelluric data, enhance signal fidelity, and improve the ability of deep noise processing and weak signal protection. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] In order to more clearly illustrate the embodiments of this specification or the technical solutions in the prior art, the following briefly introduces the drawings required for the embodiments or the description of the prior art. Obviously, the drawings described below are only some of the embodiments described in this specification. For those skilled in the art, other drawings can be obtained based on these drawings without inventive work. In the drawings:

[0015] Figure 1 A flowchart of a method for denoising multi-level magnetotelluric data provided in an embodiment of this specification;

[0016] Figure 2 A schematic diagram of evaluation parameters after processing noise signals using different methods provided in the embodiments of this specification;

[0017] Figure 3 A schematic diagram of the denoising effect of measured magnetotelluric data provided in an embodiment of this specification;

[0018] Figure 4This is a schematic diagram of the structure of a multi-level magnetotelluric data noise reduction device provided in an embodiment of this specification. DETAILED DESCRIPTION

[0019] To help those skilled in the art better understand the technical solutions in this specification, the following will provide a clear and complete description of the technical solutions in the embodiments of this specification, in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of this specification, not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of this specification without creative work should fall within the scope of protection of this specification.

[0020] The embodiments of this specification provide a method for reducing noise in multi-level magnetotelluric data. It should be noted that the execution subject in the embodiments of this specification can be a server or any device with data processing capabilities. Figure 1 A flow chart of a method for denoising multi-level magnetotelluric data provided in an embodiment of this specification is shown in FIG. Figure 1 As shown, it mainly includes the following steps:

[0021] Step S101: determine target geological analysis requirement information, collect an original magnetotelluric data sequence corresponding to a target geological analysis area in the target geological analysis requirement information, perform resolution singular value decomposition on the original magnetotelluric data sequence to determine first denoised data, perform variational mode decomposition on the first denoised data using a predetermined set of optimization parameters, and determine second denoised data.

[0022] The optimization parameter set is determined based on a pre-built electromagnetic interference feature library corresponding to the target geological analysis area.

[0023] In one embodiment of this specification, target geological analysis requirements are determined, including the geological analysis type and target geological analysis area. Based on the target geological analysis area, a corresponding raw magnetotelluric data sequence is collected. The raw MT data sequence can be obtained and preprocessed using a V5-2000 MT depth sounder. It should be noted that the magnetotelluric data, multi-level magnetotelluric data, and MT data in this embodiment represent the same type of data, and the raw MT data sequence herein is the target noise reduction data requiring noise reduction. Based on the characteristics of MT data, to improve noise reduction, this embodiment of this specification preprocesses the raw MT data sequence after acquisition. In one embodiment of this specification, preprocessing includes de-averaging and segmentation. De-averaging aims to eliminate the DC component of the signal, uniformly distributing the data around zero for easier processing. Segmentation can reduce the impact of noise and enhance the stability of subsequent spectral analysis and impedance tensor estimation. Both processes play an important role in improving signal quality, reducing the impact of outliers, and enhancing the effectiveness of subsequent noise reduction methods.

[0024] In one embodiment of the present specification, the pre-processed data is subjected to MSVD decomposition, and the detailed components D k The MSVD reconstruction data is obtained by superposition. The steps of MSVD reconstruction are as follows:

[0025] Step (a): Assume there is a 1-D data X(t) = [ ], where N is the length of the 1-D data. In the embodiment of this specification, N represents the total number of sampling points of the original magnetotelluric data sequence, that is, the length of the original magnetotelluric data sequence. Construct a multidimensional Hankel matrix of X(t) :

[0026] (1)

[0027] Step (b): Perform singular value decomposition (SVD) to obtain the singular value matrix .

[0028] Step (c): Select the main singular value according to the size of the singular value, reconstruct the matrix using the SVD inverse operation, and restore the reconstructed matrix to the approximate component and fine weight .

[0029] Step (d): Return to step (a) and construct the matrix using the approximate components , the loop stops when the number of decomposition levels reaches the optimal level.

[0030] According to the above steps, the MT data sequence that needs to be denoised is subjected to MSVD decomposition, and the detailed components of each layer are The data are superimposed to obtain MSVD reconstructed data, namely the first noise reduction data.

[0031] The first denoised data is subjected to variational modal decomposition using a predetermined set of optimization parameters to determine the second denoised data, specifically comprising: performing parameter optimization on the variational modal decomposition using a Grey Wolf optimization algorithm and a pre-constructed electromagnetic interference feature library corresponding to the target geological analysis area to determine an optimization parameter set, wherein the optimization parameter set includes the number of modal decompositions and a quadratic penalty factor; performing variational modal decomposition on the first denoised data using the optimization parameter set to obtain modal components of the modal decomposition number; calculating an instantaneous frequency threshold for each of the modal components, and determining the second denoised data using a Thomson function and the instantaneous frequency threshold.

[0032] In one embodiment of this specification, the Grey Wolf Optimizer (GWO) algorithm is used to jointly optimize the modal decomposition number k and the quadratic penalty factor α of the VMD algorithm. During the optimization process, a feature vector library of electromagnetic interference (including power frequency harmonics, square waves, triangle waves, etc.) in the target geological analysis area is added as noise prior information. The first denoised data after MSVD denoising is decomposed into the optimal modal components using the selected k and α through VMD, and the instantaneous frequency threshold of each modal component is calculated. The Thomson function is used to calculate the optimal modal component. The portion of each modal component exceeding the threshold is weighted, and the portion below the threshold is retained. The optimized modal components are summed to reconstruct the MT data main sequence, which is the second noise reduction data.

[0033] The variational modal decomposition is parameterized and an optimized parameter set is determined using the Grey Wolf optimization algorithm and a pre-constructed electromagnetic interference feature library corresponding to the target geological analysis area. Specifically, the method includes: collecting magnetotelluric noise samples of various noise types within the target geological analysis area to extract electromagnetic interference feature vectors corresponding to each noise type, and constructing an electromagnetic interference feature library based on the electromagnetic interference feature vectors corresponding to each noise type; determining a noise matching evaluation index corresponding to the decomposed modal component based on the decomposed modal component after the resolution singular value decomposition and a plurality of the electromagnetic interference feature vectors; and determining a fitness function using the noise matching evaluation index and a preset envelope entropy function, so as to optimize the variational modal decomposition parameters using the Grey Wolf optimization algorithm and determine the optimized number of modal decompositions and the quadratic penalty factor.

[0034] In one embodiment of the present specification, the steps of the VMD algorithm are as follows:

[0035] First, the overall constrained variational model of VMD is formulated as follows:

[0036] (2)

[0037] in, is the signal to be analyzed, is the modal function, is the center frequency of each mode, The unit impulse function is used for Hilbert transform, Represents the partial derivative with respect to time. Introducing the quadratic penalty factor and Lagrange multipliers , converting the constrained variational problem into an unconstrained variational problem, the model is expressed as follows:

[0038] (3)

[0039] Where, To ensure the accuracy of signal reconstruction, Make the constraints more stringent. Finally, VMD will analyze the signal Adaptively decompose into Narrowband modal components . In the process of studying VMD decomposition, it was found that the accuracy of VMD decomposition is mainly affected by the number of modal decompositions k and the quadratic penalty factor α. These two parameters must be determined in advance in practical applications. Therefore, the embodiment of this specification uses the GWO algorithm to optimize the key parameters of VMD. And in the optimization process, a library of characteristic vectors of typical electromagnetic interference (including power frequency harmonics, square waves, triangle waves, etc.) is introduced as noise prior information. In the fitness function based on envelope entropy, the matching evaluation index of the current decomposition result and the predefined noise characteristic vector is superimposed to achieve target-guided adjustment of VMD parameters.

[0040] The GWO algorithm is used to optimize the key parameters of VMD. The mathematical model of the distance D between the gray wolf and its prey and the gray wolf's position update X(t+1) is as follows:

[0041]

[0042]

[0043]

[0044]

[0045]

[0046] Where t is the current iteration number, , X are the position vectors of the prey and the wolf respectively; A and C are coefficients; is the convergence factor, which linearly increases from 2 to 0 with each iteration; r1 and r2 are random vectors generated in [0,1], and M is the maximum number of iterations. The wolf pack, led by wolf α, hunts its prey. Its hunting model is as follows:

[0047]

[0048]

[0049]

[0050] In the formula, α, β, and δ represent the three alpha wolves, and ω is regarded as the eliminated wolf. 、 、 are the direction vectors between the prey and the three wolf leaders. Z represents the distance between the prey and the wolf pack. 、 、 is the distance between the alpha, beta, delta and the prey. 、 、 Represent the positions of the three alpha wolves. Vectors A and C are coefficient vectors. Z (t+1) represents the distance between the prey and the wolf pack after t iterations.

[0051] Magnetotelluric signals are often mixed with various noise types, including power frequency harmonics, transient pulses, and square wave oscillations. This is particularly true for different geological analysis areas, where diverse noise types may exist. Traditional fixed-threshold filtering methods are inadequate for the noise characteristics of magnetotelluric signals. Different noise types exhibit significant time-frequency variations, such as narrowband power frequency noise and broadband burst pulse noise. A single filtering model cannot adaptively distinguish and process these noise types. Directly discarding specific singular value components can also accidentally damage weak, useful signals, such as low-frequency components of geological tectonic responses. By collecting typical noise samples within the target geological analysis area (such as substation power frequency interference, lightning pulses, and switching power supply oscillations), their electromagnetic interference feature vectors are extracted to form a priori knowledge base, i.e., a regional electromagnetic interference feature library corresponding to the target geological analysis area. Each modal component after SVD decomposition is compared with the feature library, and a noise matching index is calculated to accurately identify the noise properties of the component. This avoids the simplistic component discarding or retention that is often the case with traditional methods, and instead achieves intelligent weighted fusion based on the perception of noise types within the analysis area.

[0052] Specifically, based on the geographic analysis information of the target geological analysis area, collection points are deployed within high-interference areas within the region. For example, high-precision MT collection stations are deployed in strong-interference areas such as substations and industrial areas. These stations continuously collect pure noise fragments, potentially exceeding 10,000 groups, covering core noise types such as 50 / 60 Hz power frequency harmonics and their multiples, square wave pulses generated by switchgear operation, and triangular wave pulses corresponding to inverter disturbances. Each sample is labeled with the noise type, intensity level (dBμV / m), and time segment. Anti-aliasing sampling technology is used with a sampling rate of ≥1kHz to ensure complete capture of pulse details. Noise events are then isolated using a sliding time window, which can be 0.5s long with a 30% overlap.

[0053] Feature extraction is performed on the collected geomagnetic noise samples to obtain the electromagnetic interference feature vector corresponding to each noise type. In one example, a 4096-point FFT transform is performed on each noise segment. The energy contribution of the 0-100Hz frequency band is calculated, the total energy of the 0-100Hz frequency band is calculated, and the low-frequency energy of the 0-50Hz frequency band is calculated. The band energy contribution is the ratio of the total energy to the low-frequency energy. The band energy contribution is used to quantify the degree of high-frequency noise pollution. For example, the power frequency harmonic band energy contribution is ≈1.2, and the impulse noise band energy contribution is >2.5. Furthermore, time-domain kurtosis analysis is performed on the noise samples to extract kurtosis characteristics. Generally, square wave pulses have sharp transitions with a kurtosis of >15, while triangular wave pulses have smooth transitions with a kurtosis of ≈3. Furthermore, zero-crossing rate detection is performed on each noise sample to determine the number of times the signal crosses zero level per unit time. Based on these three features, feature extraction is performed on multiple samples of each noise type to construct a feature triplet corresponding to each noise type and determine the electromagnetic interference feature vector. It should be noted that the median of each feature can be used as the template value to construct an electromagnetic interference feature library in the above manner. Based on the decomposed modal components after resolution singular value decomposition and the multiple electromagnetic interference feature vectors in the electromagnetic interference feature library, the noise matching evaluation index corresponding to the decomposed modal components is determined. In this process, the same processing as the feature library construction is performed on each decomposed modal component, with segmented and synchronous triple features calculated using the same window length. The Euclidean distance between the decomposed modal component and each electromagnetic interference feature vector in the electromagnetic interference feature library is calculated, and the minimum Euclidean distance is taken as the noise matching evaluation index corresponding to this decomposed modal component.

[0054] Through the above technical solution, traditional methods rely on manually set thresholds or fixed-band filtering, which are difficult to deal with variable-frequency harmonics and composite noise. By matching the characteristic vectors of the noise in the target geological analysis area, automatic identification and separation of noise types are achieved, effectively improving the detection rate of power-frequency noise; the feature library supports online updates, such as adding high-speed rail interference samples, which can dynamically expand the noise type identification capability.

[0055] The fitness function is determined using the noise matching evaluation index and a preset envelope entropy function, specifically including: generating a dynamic weight coefficient based on the proportional relationship between the current envelope entropy value and the preset maximum reference entropy value, adjusting the noise matching evaluation index based on the dynamic weight coefficient, and determining an environmental noise reference item; determining a signal purity item based on the envelope entropy function, and determining the fitness function through the signal purity item and the environmental noise reference item.

[0056] In one embodiment of the present specification, the improved envelope entropy function and the noise matching evaluation index are used as the fitness function, and the mathematical model of the standard envelope entropy is:

[0057] (4)

[0058] is the envelope entropy; is the normalized envelope energy of the modal component (IMF) obtained by the j-th variational mode decomposition; is the envelope energy of the IMF, obtained through Hilbert transform. Formula (4) extracts the time domain envelope signal of the modal component through Hilbert transform, calculates the normalized probability of its energy distribution, and then obtains the Shannon entropy value that characterizes the randomness of the signal.

[0059] A two-dimensional function including signal purity evaluation and noise separation evaluation is constructed to dynamically balance the guiding role of the two evaluation indicators in parameter optimization. The overall form of the fitness function is as follows:

[0060] (5)

[0061] in, is the signal purity term determined based on the envelope entropy function, To adjust the noise matching evaluation index based on the dynamic weight coefficient, the environmental noise reference item is determined. Specifically, { u i ( t )} is the IMF component set; F IMF is the characteristic vector of all current IMFs; F noise,j Represents the feature vector of the j-th noise template. Among them, Dist(F IMF ,F noise,j ) = ∥F IMF -F noise,j ∥2 is the noise matching evaluation index, which is used to measure the similarity between the characteristics of the IMF in the decomposition result and the noise template, that is, the minimum distance between the feature vector of the decomposition result and the feature vector of the template in the noise library.

[0062] According to the current envelope entropy value and the preset maximum reference entropy value The proportional relationship is used to generate dynamic weight coefficients, and a negative correlation adjustment mechanism between the envelope entropy value and the weight coefficient is established, so that the fitness function automatically switches the optimization focus under different noise environments. The calculation method of the dynamic weight coefficient is as follows:

[0063] (6)

[0064] Among them, λ is the dynamic weight coefficient, λ0 is the maximum weight factor, and the initial value is set to 0.5. Indicates the maximum reference value of envelope entropy.

[0065] Calculating the instantaneous frequency threshold for each modal component specifically includes: determining the median value of all sampling points of the modal component, calculating the absolute deviation of the value of each sampling point in the modal component relative to the median value, and determining the median dispersion index of the absolute deviation; calculating the signal kurtosis distribution morphological parameter of the modal component, and determining the morphological adjustment coefficient based on the ratio of a preset adjustment constant to the signal kurtosis distribution morphological parameter; determining a basic threshold based on the ratio of the median dispersion index to the morphological adjustment coefficient, and performing environmental compensation on the basic threshold using the total number of sampling points obtained in advance to determine the instantaneous frequency threshold.

[0066] Through the Thomson function and instantaneous frequency threshold, the Thomson function is a function used in the field of electromagnetic fields to process abnormal values ​​exceeding the threshold in the modal component. The portion of each modal component exceeding the threshold is weighted, and the portion of the modal component that does not exceed the threshold is retained. The optimized modal components are added together to reconstruct the MT data main sequence, which is the second noise reduction data. The specific implementation process of weighting the portion of each modal component exceeding the threshold is as follows:

[0067] Frequency threshold and Thomson function for:

[0068] (7)

[0069] (8)

[0070] Where N represents the total number of sampling points of the original magnetotelluric data sequence, is the value of the jth sampling point in the i-th IMF taken in the entire time domain, For The median of the sequence is the median value of the jth sampling point in the i-th IMF. is the absolute median deviation of the ith IMF, i.e., the absolute deviation, which can detect the protruding impulse noise. is the adaptive coefficient of the i-th IMF component, that is, the morphological adjustment coefficient, , is the kurtosis of each IMF component, i.e., the signal kurtosis distribution morphology parameter, and c is the adjustment coefficient set to 0.6745.

[0071] by The process of removing impulse noise from magnetotelluric data for a threshold is as follows:

[0072] (9)

[0073] Formula (9) represents the processing method for each sampling point j of each IMF component (i.e., the i-th modal component). For each sampling point j of each IMF component, first determine whether the absolute value of the sampling point is greater than the threshold ,Right now If it is greater than the threshold, the first rule is executed ; If it is not greater than the threshold (i.e. less than or equal to), then execute the second rule . is the instantaneous frequency threshold calculated for the i-th modal component. This threshold is used to distinguish normal signal points from abnormal noise points, usually the critical value of impulse noise. When a sampling point is considered to be a noise point (i.e., exceeds the threshold), the Thomson function is used. Apply a weighting (falloff) to the point. is a value between 0 and 1 that exponentially attenuates the value exceeding the threshold, thereby suppressing noise. Points that do not exceed the threshold are considered valid signals, so their original values ​​are retained. It should be noted that this formula is actually a soft thresholding process. Compared to hard thresholding, which directly sets points exceeding the threshold to zero, soft thresholding, by multiplying by a factor less than 1, can more smoothly attenuate noise and avoid signal distortion.

[0074] Step S102: Determine a first signal-to-noise ratio corresponding to the second denoised data and a signal-to-noise ratio reference value corresponding to the target geological analysis requirement. When the first signal-to-noise ratio is less than a preset signal-to-noise ratio reference value, dynamically fuse the residual signal based on the first denoised data and the second denoised data to determine a fused residual signal, so as to analyze weak signals in the residual signal and avoid discarding weak signals.

[0075] In one embodiment of the present specification, a first signal-to-noise ratio (SNR) corresponding to the second noise reduction data is determined, and the signal-to-noise ratio (SNR) formula is:

[0076] (10)

[0077] Where, and are the original effective signal and the denoised synthetic time series, i.e., the amplitude of the data series corresponding to the second denoised data, and n is the sampling point. When the first signal-to-noise ratio is less than the preset signal-to-noise ratio reference value, it means that the noise reduction effect at this time does not meet the requirements.

[0078] In the field of magnetotelluric exploration, a critical threshold for ensuring the basic usability of the apparent resistivity-phase curve is that the useful signal power is at least three times the noise power (4.7dB). However, obtaining reliable data that can be used for detailed inversion interpretation generally requires a higher signal-to-noise ratio. The preset signal-to-noise ratio reference value needs to be determined based on different working conditions and industry standards. It should be noted that the signal-to-noise ratio reference value in the embodiments of this specification can be directly set to 8dB based on experience. Alternatively, a signal-to-noise ratio reference value corresponding to the target geological analysis requirements can be set. When the first signal-to-noise ratio is not less than the preset signal-to-noise ratio reference value, the second de-noised data is used as the final de-noised output signal.

[0079] To determine the signal-to-noise ratio reference value corresponding to the target geological analysis requirements, the following methods can be used: Combined with the target geological body parameters within the target geological analysis area, which include the resistivity and burial depth of the target geological body, the geological signal-to-noise ratio reference value corresponding to the geological conditions is set based on empirical data. For example, through testing, statistics can be collected to determine the multiple signal-to-noise ratios that meet the analysis requirements when performing geological analysis under different geological conditions. Through statistical analysis, the geological signal-to-noise ratio reference values ​​corresponding to different geological conditions are obtained, and a mapping relationship table is constructed. The geological conditions in this table can be differentiated through geological classification.

[0080] After determining the geological signal-to-noise ratio reference value corresponding to the geological conditions, the signal-to-noise ratio is corrected based on the task type of the target geological analysis task. It should be noted that the task types here include structural boundary positioning, lithologic identification, and fluid detection. A corresponding signal-to-noise ratio correction value is pre-set for each task type, such as 0 dB for structural boundary positioning, 1.5 ± 0.3 dB for lithologic identification, and 2.5 ± 0.5 dB for fluid detection. This ensures that the obtained signal-to-noise ratio reference value matches the target geological analysis requirements.

[0081] Because effective magnetotelluric signals (such as responses from deep fault zones) are typically weak at the microvolt level and easily masked by noise, the initial noise reduction process may oversuppress useful components that overlap with the noise band, such as signals in the 0.1-1 Hz low-frequency band. When the initial signal-to-noise ratio (SNR) falls below the preset SNR reference value, it indicates that the residual noise energy in the current noise reduction result exceeds the effective signal, and the subsequent residual recovery process is performed.

[0082] Residual recovery is used to solve the technical problem of the interweaving of stubborn noise and deep weak signals in magnetotelluric exploration. When the signal-to-noise ratio of the noise reduction channel output does not reach the preset threshold, it means that the conventional noise reduction process is difficult to balance the contradiction between noise suppression and signal fidelity. Excessive noise suppression will damage the microvolt-level geological response implicit in the residual, while retaining the residual will cause the noise to continue to pollute the effective signal. Through an innovative dynamic fusion mechanism, the two types of residual signals left over from the initial noise reduction are converted into a resource library that can be purified for the second time, achieving accurate separation of noise and weak signals. The traditional method of directly discarding the residuals actually permanently discards some effective signals, especially leading to the irreversible loss of key information such as the fault zone characteristics and thin coal seam interfaces of deep strata responses.

[0083] According to the first noise reduction data and the second noise reduction data, the residual signal is dynamically fused to determine the fused residual signal, specifically including: determining the first residual data corresponding to the first noise reduction data and the second residual data corresponding to the second noise reduction data; respectively calculating the first time domain energy data of the first residual data and the second time domain energy data of the second residual data to generate a dynamic weight parameter combination based on the first time domain energy data and the second time domain energy data; and performing weighted fusion on the first residual data and the second residual data through the dynamic weight parameter combination to determine the fused residual signal.

[0084] In one embodiment of the present specification, the first residual data corresponding to the first denoised data is the mixed component remaining after the MSVD preliminary denoising, and the first denoised data is the MSVD reconstructed signal. , through the formula , calculate the first residual data, where is the first residual data, is the original sampling sequence, i.e. the original magnetotelluric data sequence. The second denoised data is the reconstructed signal after VMD secondary denoising. , the second residual data corresponding to the second noise reduction data is the high-frequency pulse and oscillation noise remaining after the second noise reduction. , calculate the second residual data .

[0085] Based on the first time domain energy data and the second time domain energy data, a dynamic weight parameter combination is generated, specifically including: determining the total energy of the first time domain energy data and the second time domain energy data, and generating a first weight coefficient corresponding to the first residual data based on the ratio of the first time domain energy data to the total energy; and determining a second weight coefficient corresponding to the second residual data based on the first weight coefficient and a preset normalization constraint.

[0086] In one embodiment of the present specification, the fusion residual signal The formula is:

[0087] (11)

[0088] in, 、 To set weights based on residual signal energy, the default normalization constraint is to satisfy .

[0089] (12)

[0090] in, is the residual between the original data and the data after MSVD denoising (i.e. the residual of the first denoising), The residual between the reconstructed signal data after MSVD-VMD joint denoising and the MSVD reconstructed signal (i.e., the residual of the second denoising) is weighted fused by the two residuals using formula (11) to obtain the fused residual signal . Calculate the first time domain energy data of the first residual data and the second time domain energy data of the second residual data respectively. The signal energy here is defined as the sum of the squares of the amplitudes of each sampling point. Represents the squared delta-normal (i.e., energy) of the signal. When determining the first weight coefficient in the dynamic fusion weight combination, the squared delta-normal value of the first residual data is divided by the sum of the squared delta-normal values ​​of the residual signal and the squared delta-normal value of the second residual data. Using the dynamic weight parameter combination, the first and second residual data are weightedly fused to determine the fused residual signal. When the MSVD residual energy ratio is high, low-frequency noise processing is enhanced. When the MSVD-VMD residual energy ratio is high, high-frequency noise suppression is enhanced.

[0091] Through the above-mentioned residual recovery method, the cyclic regeneration of signal components is achieved. The low-frequency weak signals implicit in the MSVD residuals and the high-frequency pulses remaining in the VMD residuals are dynamically weighted and fused. The low-frequency weak signals implicit in the SVD residuals are weighted and enhanced because they carry deep geological responses. The high-frequency pulses remaining in the VMD residuals are gradually attenuated due to their pure noise properties, allowing the deep weak signals to be recovered and significantly repairing the distortion of the apparent resistivity curve.

[0092] Step S103 , performing noise reduction processing on the fused residual signal to determine corresponding fused residual noise reduction data, performing residual recovery based on the fused residual noise reduction data, and determining current noise reduction data after residual recovery.

[0093] In one embodiment of the present specification, the fused residual signal after fusion is denoised using the MSVD-VMD main denoising channel. After obtaining the denoised fused residual denoised data, the fused residual denoised data is subjected to residual recovery and merged with the main signal to form a reconstructed MT data sequence. The main signal here is the reconstructed signal after secondary denoising by MSVD-VMD, i.e., the second denoised data, to determine the current denoised data after residual recovery. It should be noted that in subsequent iterative cycles, the second denoised data is used as the original MT data sequence for denoising and residual recovery. The recovered residual signal is merged with the reconstructed signal after secondary denoising by MSVD-VMD in the current cycle to determine the denoised data after residual recovery.

[0094] Step S104 : when the second signal-to-noise ratio corresponding to the current noise reduction data is not less than the signal-to-noise ratio reference value, output the current noise reduction data to perform geological analysis on the target geological analysis area based on the current noise reduction data.

[0095] In one embodiment of the present specification, a signal-to-noise ratio is calculated for the current denoised data after residual recovery. If a second signal-to-noise ratio corresponding to the current denoised data is not less than a reference signal-to-noise ratio value, the current denoised data is output, and the denoising process ends. If the second signal-to-noise ratio corresponding to the current denoised data is less than the reference signal-to-noise ratio value, the current denoised data is output, and the denoising process ends.

[0096] The method further includes: after determining the current denoised data after residual recovery, performing a denoising count operation to determine the current denoising iteration number; and outputting the current denoising data when the current denoising iteration number reaches a preset iteration number threshold.

[0097] In one embodiment of the present specification, after determining the current denoised data after residual recovery, a denoising count operation is performed to determine the current denoising iteration count. In other words, the denoising process after residual recovery is considered the first denoising step. When the current denoising iteration count reaches a preset threshold, the current denoised data is output to avoid excessive iterations. Subsequent geological analysis of the target geological analysis area is performed based on the current denoised data to ensure signal availability.

[0098] In order to further analyze the noise reduction performance of the embodiments of this specification, different evaluation parameters are used to conduct an in-depth analysis of its stability and effectiveness. Figure 2 This is a schematic diagram of the evaluation parameters after processing noise signals using different methods provided in the embodiments of this specification. Figure 3 To further verify the denoising capability of the embodiment of this specification, in order to simulate the effect of noise on the effective signal at different signal-to-noise ratio levels, noise of different intensities is introduced into the simulated signal. Figure 3 This is a schematic diagram of the denoising effect of measured magnetotelluric data provided in the embodiment of this specification, as shown in FIG. Figure 3 As shown in Figure 1, similar and continuous noise interference can be observed in the Ex and Ey channels, while complex noise mixed with pulse signals exists in the Hx and Hy magnetic channels. The amplitude of these interference signals far exceeds the effective MT signal, resulting in a low signal-to-noise ratio for the MT data. The MSVD-VMD algorithm with dynamic residual feedback is used to reconstruct the denoising results, as shown in Figure 1. Figure 3 Middle (e)- Figure 3 As shown in (h), the noise amplitude is reduced, the square wave noise in the electrical channel is effectively suppressed, and the effective signal is preserved. The waveform of the magnetic channel is smoother and more continuous, without abnormal spikes, and closely matches the time domain characteristics of natural magnetotelluric signals. These experimental results demonstrate the effectiveness and superiority of the embodiments of this specification in noise data processing.

[0099] Through the technical solutions of the embodiments of this specification, the original magnetotelluric data sequence corresponding to the target geological analysis area in the target geological analysis requirement information is collected, ensuring the regional matching between the signal source and the analysis requirement; the MSVD algorithm is introduced to perform preliminary noise reduction on the non-stationary signal MT signal by multi-level decomposition. Compared with the SVD of single-level linear decomposition, the low-frequency noise in the MT signal can be processed at multiple scales and with high resolution, which can effectively retain the target signal while suppressing the noise; in order to solve the problems such as modal aliasing generated in the decomposition process of the residual pulse signal, the Thomson function is introduced on the basis of the VMD algorithm, and the components exceeding the threshold are adaptively soft-thresholded to avoid hard truncation distortion, which can more effectively suppress the synthesis of the pulse signal than the traditional VMD algorithm; using GW The high convergence characteristics of the O algorithm are used to optimize key parameters in VMD. During the optimization process, a library of characteristic vectors of electromagnetic interference within the target geological analysis area is introduced as noise prior information. By constructing a composite fitness function that integrates envelope entropy, decomposition results, and noise feature matching, target-guided adjustment of VMD parameters is achieved. This ensures that the noise prior knowledge matches the target geological analysis area. This not only significantly reduces the possibility of modal aliasing, but also improves computational efficiency and significantly shortens calculation time by optimizing the number of modes. A dynamic residual feedback method and a weighted fusion mechanism are used to further process the effective signals in the residuals, meeting the signal recovery requirements of weak effective information in magnetotelluric data, enhancing signal fidelity, and improving the capabilities of deep noise processing and weak signal protection.

[0100] The embodiment of this specification also provides a multi-level magnetotelluric data noise reduction device, such as Figure 4 As shown, the device includes: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions that can be executed by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform the above method.

[0101] The embodiments of this specification also provide a non-volatile computer storage medium storing computer executable instructions, wherein the computer executable instructions are configured to execute the above method.

[0102] The various embodiments in this specification are described in a progressive manner. Similar portions between the various embodiments can be referenced to each other, and each embodiment focuses on the differences from the other embodiments. In particular, the device, apparatus, and non-volatile computer storage medium embodiments are generally similar to the method embodiments, so their descriptions are relatively simplified. For relevant details, refer to the descriptions of the method embodiments.

[0103] The foregoing description is merely one or more embodiments of this specification and is not intended to limit this specification. It will be apparent to those skilled in the art that various modifications and variations may be made to one or more embodiments of this specification. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of one or more embodiments of this specification are intended to be within the scope of the claims of this specification.

Claims

1. A method for denoising multi-level magnetotelluric data, characterized in that: The method comprises: Determining target geological analysis requirement information, collecting an original magnetotelluric data sequence corresponding to a target geological analysis area in the target geological analysis requirement information, performing resolution singular value decomposition on the original magnetotelluric data sequence to determine first de-noised data, and performing variational mode decomposition on the first de-noised data using a predetermined set of optimization parameters to determine second de-noised data, wherein the set of optimization parameters is determined based on a pre-constructed electromagnetic interference feature library corresponding to the target geological analysis area; determining a first signal-to-noise ratio corresponding to the second de-noised data and a signal-to-noise ratio reference value corresponding to the target geological analysis requirement; and when the first signal-to-noise ratio is less than a preset signal-to-noise ratio reference value, dynamically fusing a residual signal based on the first de-noised data and the second de-noised data to determine a fused residual signal, so as to facilitate analysis of weak signals in the residual signal and avoid discarding the weak signals; performing noise reduction processing on the fused residual signal to determine corresponding fused residual noise reduction data, performing residual recovery based on the fused residual noise reduction data, and determining current noise reduction data after residual recovery; When the second signal-to-noise ratio corresponding to the current noise reduction data is not less than the signal-to-noise ratio reference value, the current noise reduction data is output to perform geological analysis on the target geological analysis area based on the current noise reduction data.

2. The method for denoising multi-level magnetotelluric data according to claim 1, characterized in that: Performing variational mode decomposition on the first denoised data using a predetermined set of optimization parameters to determine second denoised data specifically includes: Optimizing the variational modal decomposition parameters using a Grey Wolf optimization algorithm and a pre-built electromagnetic interference signature library corresponding to the target geological analysis area to determine an optimized parameter set, wherein the optimized parameter set includes the modal decomposition quantity and a quadratic penalty factor; Performing variational modal decomposition on the first denoised data using the optimized parameter set to obtain the modal decomposition number of modal components; An instantaneous frequency threshold is calculated for each of the modal components, and the second noise reduction data is determined by using a Thomson function and the instantaneous frequency threshold.

3. The method for denoising multi-level magnetotelluric data according to claim 2, characterized in that: The variational modal decomposition is optimized using the Grey Wolf optimization algorithm and a pre-built electromagnetic interference signature library corresponding to the target geological analysis area to determine an optimized parameter set, specifically including: Collecting magnetotelluric noise samples of various noise types in the target geological analysis area to extract electromagnetic interference feature vectors corresponding to each noise type, and constructing an electromagnetic interference feature library based on the electromagnetic interference feature vectors corresponding to each noise type; Determining, based on the decomposed modal components after the resolution singular value decomposition and the plurality of electromagnetic interference eigenvectors, a noise matching evaluation index corresponding to the decomposed modal components, so as to quantify the environmental noise interference in the target geological analysis area; The fitness function is determined by using the noise matching evaluation index and a preset envelope entropy function, so as to optimize the parameters of the variational modal decomposition through the gray wolf optimization algorithm and determine the optimized modal decomposition quantity and quadratic penalty factor.

4. The method for denoising multi-level magnetotelluric data according to claim 3, characterized in that: The fitness function is determined based on the noise matching evaluation index and the preset envelope entropy function, specifically including: Generating a dynamic weight coefficient according to a proportional relationship between a current envelope entropy value and a preset maximum reference entropy value, and adjusting the noise matching evaluation index based on the dynamic weight coefficient to determine an environmental noise reference item; Based on the envelope entropy function, a signal purity term is determined, and the fitness function is determined by the signal purity term and the environmental noise reference term.

5. The method for denoising multi-level magnetotelluric data according to claim 2, characterized in that: Calculating an instantaneous frequency threshold for each of the modal components specifically includes: Determine the median value of all sampling points of the modal component, calculate the absolute deviation of each sampling point value in the modal component relative to the median value, and determine the median dispersion index of the absolute deviation; Calculating a signal kurtosis distribution morphological parameter of the modal component to determine a morphological adjustment coefficient based on a ratio of a preset adjustment constant to the signal kurtosis distribution morphological parameter; A basic threshold is determined according to the ratio of the median dispersion index to the morphological adjustment coefficient, and the instantaneous frequency threshold is determined by performing environmental compensation on the basic threshold using the total number of sampling points acquired in advance.

6. The method for denoising multi-level magnetotelluric data according to claim 1, characterized in that: Dynamically fusing the residual signal according to the first denoised data and the second denoised data to determine a fused residual signal specifically includes: Determining first residual data corresponding to the first denoised data and second residual data corresponding to the second denoised data; respectively calculating first time-domain energy data of the first residual data and second time-domain energy data of the second residual data, so as to generate a dynamic weight parameter combination based on the first time-domain energy data and the second time-domain energy data; The first residual data and the second residual data are weightedly fused by combining the dynamic weight parameters to determine the fused residual signal.

7. The method for denoising multi-level magnetotelluric data according to claim 6, characterized in that: Generating a dynamic weight parameter combination based on the first time-domain energy data and the second time-domain energy data specifically includes: Determine the total energy of the first time-domain energy data and the second time-domain energy data, and generate a first weight coefficient corresponding to the first residual data based on a ratio of the first time-domain energy data to the total energy; A second weight coefficient corresponding to the second residual data is determined according to the first weight coefficient and a preset normalization constraint condition.

8. The method for denoising multi-level magnetotelluric data according to claim 1, characterized in that: The method further comprises: After determining the current denoising data after residual recovery, a denoising count operation is performed to determine the current denoising iteration number; When the current denoising iteration number reaches a preset iteration number threshold, the current denoising data is output.

9. A data denoising device for multi-level magnetotelluric data, characterized in that: The device comprises: at least one processor; and, a memory communicatively connected to the at least one processor; wherein, The memory stores instructions that can be executed by the at least one processor. The instructions are executed by the at least one processor to enable the at least one processor to perform the method according to any one of claims 1 to 8.

10. A non-volatile computer storage medium storing computer executable instructions, characterized in that: The computer executable instructions are configured to execute the method according to any one of claims 1 to 8.

Images