A method for estimating the instantaneous frequency spectrum of magnetotelluric signals

Abstract

This invention provides a method for instantaneous spectrum estimation of magnetotelluric signals, relating to the field of geophysical electromagnetic exploration signal processing technology. The invention acquires four magnetotelluric signals and selects one as the signal to be processed; it constructs a local mean operator and a sieving operator using a double-average filtering function, repeatedly sieving to obtain a first-order temporary IMF component; it obtains the phase using Hilbert transform and calculates the instantaneous frequency using forward differential; it constructs a resampling function based on the instantaneous frequency and interpolates to obtain a resampling sequence; it performs discrete Fourier transforms on the resampling sequence and the filtered loop sequence respectively to obtain two spectra; it iterates by multiplying corresponding terms of the two spectra, and after convergence, performs an inverse transform to obtain the IMF component, which is then discarded for further decomposition; it shapes the instantaneous phase waveform to obtain a non-negative instantaneous frequency and constructs the time spectrum, repeating the processing on the remaining three signals. This method improves convergence, avoids negative frequencies, and enhances the reliability of the time spectrum.

Description

Technical Field

[0001] This invention relates to the field of geophysical electromagnetic exploration signal processing technology, and in particular to a method for instantaneous spectrum estimation of magnetotelluric signals. Background Technology

[0002] Magnetotelluric sounding uses natural electromagnetic fields as the field source. The observed signals are characterized by non-stationarity, complexity, multiple sources, weak effective signals, and susceptibility to human noise interference. Regarding reliable impedance estimation and denoising, existing technologies have formed several approaches, including time-domain denoising and time-domain impedance estimation, frequency-domain denoising and impedance estimation, and time-spectrum processing-based denoising and impedance estimation. Among them, the time-spectrum approach can dynamically reflect the changes in the spectrum of non-stationary signals over time and suppress false anomalies.

[0003] Current methods for obtaining the instantaneous spectrum are mostly based on the sliding window short-time spectrum or the instantaneous spectrum based on the Hilbert transform. The short-time Fourier transform is limited by a fixed time window, making it difficult to balance time and frequency resolution. The Hilbert transform is an integral transform, and the obtained instantaneous frequency is easily an average result within the integration range. Furthermore, the phase derivative may have a negative frequency without physical meaning, making it difficult to reflect the superposition of multiple frequencies at the same moment. At the same time, the components that rely on recursive decomposition may also have problems such as mode mixing or insufficient physical meaning. Therefore, a technical solution that can stably decompose non-stationary signals and reliably construct the instantaneous spectrum is still needed. Summary of the Invention

[0004] To overcome the shortcomings of existing technologies, the purpose of this invention is to provide a method for estimating the instantaneous spectrum of magnetotelluric signals. By combining improved fast resampling iterative filtering with phase waveform shaping, the decomposition convergence is enhanced and negative frequencies are avoided, thereby improving the reliability and applicability of the instantaneous frequency and time-spectrum results.

[0005] To achieve the above objectives, the present invention provides the following solution:

[0006] A method for instantaneous spectrum estimation of magnetotelluric signals includes:

[0007] S1. Acquire four magnetotelluric signals and select one of them as the signal to be processed;

[0008] S2. Define a local mean operator and a sieving operator using a double average filter function, and repeatedly sieve the signal to be processed to obtain the first-order temporary IMF component.

[0009] S3. The phase is obtained by Hilbert transformation of the first-order temporary IMF component, and the instantaneous frequency is obtained by forward differential.

[0010] S4. Calculate the extreme point spacing function based on the instantaneous frequency and construct the resampling function. Obtain the resampling sequence by cubic spline interpolation. Obtain the spectrum of the resampling sequence by discrete Fourier transform of the resampling sequence. Convolve the double average filter function to construct the filter loop sequence. Obtain the filter spectrum by discrete Fourier transform of the filter loop sequence.

[0011] S5. Multiply the corresponding terms of the resampled sequence spectrum and the filter spectrum one by one and iterate until the convergence condition is met. Then, perform an inverse Fourier transform to obtain the first-order IMF component. Remove the first-order IMF component from the signal to be processed and repeat S2 to S5 to extract each order of IMF component until the remaining signal has no more than two extreme points or reaches the maximum decomposition order.

[0012] S6. Regularize the IMF components and derivative sequences of each order to construct a complex sequence and obtain the instantaneous phase;

[0013] S7. After shaping the instantaneous phase into a monotonically increasing form, the forward differential is used to obtain the non-negative instantaneous frequency. The time spectrum is constructed based on the time and the non-negative instantaneous frequency.

[0014] S8, and the remaining three magnetotelluric signals are processed from S2 to S7 to obtain the time spectrum.

[0015] Preferably, four magnetotelluric signal data are acquired, including:

[0016] Acquire electric field data in the north-south direction and electric field data in the east-west direction;

[0017] Acquire magnetic track data in the north-south direction and the east-west direction of the magnetic field;

[0018] The electrical track data and the magnetic track data are used together as the magnetotelluric signal data.

[0019] Preferably, step S2 includes:

[0020] S21. Set the signal to be processed as the screening input sequence for this round of screening;

[0021] S22. The local mean operator is used to perform low-pass filtering on the sieving input sequence to obtain a local mean sequence;

[0022] S23. Subtract the local mean sequence from the screening input sequence to obtain the screening result for this round;

[0023] S24. Update the current round of screening results to the screening input sequence for the next round of screening, and repeat steps S22 to S24 until the maximum number of iterations is reached or the convergence parameter meets the threshold, to obtain the first-order temporary IMF component; wherein, the maximum number of iterations is 500; the convergence parameter is the ratio of the L2 norm of the difference between two adjacent rounds of screening results to the L2 norm of the next screening result; the threshold is 0.05.

[0024] Preferably, step S3 includes:

[0025] Perform a Hilbert transform on the first-order temporary IMF component to obtain a Hilbert transform sequence;

[0026] Construct a complex sequence using the first-order temporary IMF component as the real part and the Hilbert transform sequence as the imaginary part;

[0027] The phase is determined based on the complex sequence.

[0028] The instantaneous frequency is obtained by differentiating the phase using forward differential, and the last term of the instantaneous frequency is set to the same value as the previous term of the instantaneous frequency.

[0029] Preferably, step S4 includes:

[0030] The extreme point spacing function is calculated based on the instantaneous frequency;

[0031] A resampling function is constructed based on the extreme point spacing function, and the resampling point set is determined by the resampling function;

[0032] The signal to be processed is subjected to cubic spline interpolation on the set of resampled points to obtain a resampled sequence;

[0033] Perform a discrete Fourier transform on the resampled sequence to obtain the spectrum of the resampled sequence;

[0034] The filter function is taken as the reconvolution of the double average filter function and a filter loop sequence is constructed. The filter loop sequence is then subjected to a discrete Fourier transform to obtain the filter spectrum.

[0035] Preferably, the corresponding terms of the resampled sequence spectrum and the filter spectrum are multiplied sequentially and iterated until the convergence condition is met. Then, the inverse Fourier transform is performed to obtain the first-order IMF component, including:

[0036] S51. Use the spectrum of the resampled sequence as the spectrum of the current iteration;

[0037] S52. Multiply the current iteration spectrum with the filter spectrum by corresponding terms to obtain the updated iteration spectrum;

[0038] S53. Determine whether the difference between the updated iteration spectrum and the current iteration spectrum meets the allowable error. If not, repeat step S52 using the updated iteration spectrum as the new current iteration spectrum. The convergence condition includes that the difference between the updated iteration spectrum and the current iteration spectrum meets the allowable error.

[0039] S54. When the tolerance error is met or the iteration limit is reached, perform an inverse Fourier transform on the iterative spectrum that meets the conditions to obtain the first-order IMF component; wherein, the iteration limit is 1000; and the tolerance error is 0.05.

[0040] Preferably, the first-order IMF component is removed from the signal to be processed, and S2 to S5 are repeated to extract IMF components of each order until the remaining signal has no more than two extreme points or the maximum decomposition order is reached, including:

[0041] When the number of remaining signal extreme points does not exceed two, stop extracting the IMF components of each order one by one;

[0042] When the extraction order of each IMF component exceeds the maximum decomposition order, the extraction of each IMF component is stopped; wherein, the maximum decomposition order is an integer from ten to twenty.

[0043] Preferably, step S6 includes:

[0044] For each IMF component, a smooth curve is used to connect the maximum or minimum points of each IMF component to form the first envelope.

[0045] Based on the first envelope, each order of IMF component is regularized to obtain the regularized IMF component.

[0046] For the derivative sequence corresponding to each IMF component, a smooth curve is used to connect the maximum or minimum points of the derivative sequence to form a second envelope;

[0047] The derivative sequence is regularized based on the second envelope to obtain a regularized derivative sequence.

[0048] A complex sequence is constructed using the regularized IMF components and the regularized derivative sequence, and the instantaneous phase corresponding to each order of IMF components is obtained from the complex sequence.

[0049] Preferably, step S7 includes:

[0050] The extreme points of the instantaneous phases corresponding to each order of IMF components are numbered in chronological order.

[0051] When the extreme point numbered one is a maximum value, the instantaneous phase is shaped according to the first phase shaping rule;

[0052] When the extreme point numbered one is a minimum value, the instantaneous phase is shaped according to the second phase shaping rule so that the shaped instantaneous phase is monotonically increasing.

[0053] The non-negative instantaneous frequency is obtained by using forward differential differentiation on the shaped monotonically increasing instantaneous phase.

[0054] A time spectrum is constructed with time as the horizontal axis, the non-negative instantaneous frequency as the vertical axis, and the instantaneous amplitude or instantaneous phase as the spectral value; wherein the instantaneous amplitude is the absolute value of the IMF component of the corresponding order.

[0055] Preferably, the filter function is a reconvolution of the double-average filter function, and a filter loop sequence is constructed. A discrete Fourier transform is then performed on the filter loop sequence to obtain the filter spectrum, including:

[0056] The filter function is re-regularized to obtain a re-regularized filter function, wherein the re-regularization constant used for the re-regularization process is approximated to one.

[0057] A filter loop sequence is constructed based on the re-regularized filter function, and the length of the filter loop sequence is set to be the same as the length of the resampling sequence.

[0058] The filter spectrum is obtained by performing a discrete Fourier transform on the filtered cyclic sequence.

[0059] The present invention discloses the following technical effects:

[0060] This invention adaptively decomposes the signal to be processed using a local mean operator and a sieving operator. It can extract temporary components and components of various orders based on the local changes of the signal itself, even when the magnetotelluric signal is non-stationary, the effective signal is weak, and it is subject to complex noise interference. This avoids the problem of difficulty in balancing time resolution and frequency resolution caused by fixed time windows or fixed basis functions, thereby more accurately depicting the details of the spectrum changing over time and improving the adaptability of instantaneous spectrum estimation to complex scenarios.

[0061] This invention sets a maximum number of iterations or a stopping condition where the convergence parameter meets a threshold during the repeated screening process. After removing the first-order component, the remaining signal is used as a new signal to be processed to extract subsequent components one by one. This can suppress component distortion caused by over-screening or under-screening, reduce the risk of mode mixing, and make the obtained components have more stable scale separation and clearer frequency direction, providing a more reliable component basis for subsequent time-spectrum construction.

[0062] This invention employs forward difference to obtain instantaneous frequency by phase differentiation and performs end-padding on the instantaneous frequency sequence. Simultaneously, it performs regularization on each order component and its derivative sequence before constructing the instantaneous phase. This reduces the impact of endpoint effects and local abrupt changes on phase and instantaneous frequency calculations, reduces abnormal jumps in instantaneous frequency, and improves the continuity and interpretability of the instantaneous frequency trajectory.

[0063] This invention calculates the extreme point spacing function based on instantaneous frequency and constructs a resampling function. It uses cubic spline interpolation to obtain a resampling sequence and forms a resampling sequence spectrum. At the same time, it constructs a filter loop sequence by reconvolution of the double average filter function and forms a filter spectrum. Furthermore, it obtains the components through iterative operation of multiplying the corresponding terms of the spectrum and inverse transformation, thereby achieving the focusing of effective frequency components and the suppression of noise components. Compared with the instantaneous spectrum acquisition method that only relies on integral transformation, the instantaneous spectrum estimation result is more stable and has a higher signal-to-noise ratio.

[0064] This invention performs waveform shaping on the instantaneous phases corresponding to each order component to make them monotonically increasing, and then differentiates to obtain non-negative instantaneous frequencies. The time spectrum is constructed using time, non-negative instantaneous frequencies, and instantaneous amplitudes or instantaneous phases. This avoids interference from frequency components such as negative frequencies that lack physical meaning. At the same time, the same process is performed on the other three magnetotelluric signal data to achieve consistent time-frequency characterization of the four signals, thereby providing a more reliable time-frequency basis for subsequent denoising and impedance estimation. Attached Figure Description

[0065] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0066] Figure 1 A flowchart of the method provided in an embodiment of the present invention;

[0067] Figure 2 This is a schematic diagram of the technical route provided for an embodiment of the present invention. Detailed Implementation

[0068] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0069] The purpose of this invention is to provide a method for instantaneous spectrum estimation of magnetotelluric signals, which can achieve high time-frequency resolution instantaneous spectrum estimation for complex and non-stationary magnetotelluric signals. The obtained time-frequency spectrum is stable, reliable, and easier to interpret.

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

[0071] Figure 1 The method flowchart provided in the embodiments of the present invention is as follows: Figure 1 As shown, this invention provides a method for instantaneous spectrum estimation of magnetotelluric signals, comprising:

[0072] S1. Acquire four magnetotelluric signals and select one of them as the signal to be processed;

[0073] S2. Define a local mean operator and a sieving operator using a double average filter function, and repeatedly sieve the signal to be processed to obtain the first-order temporary IMF component.

[0074] S3. The phase is obtained by Hilbert transformation of the first-order temporary IMF component, and the instantaneous frequency is obtained by forward differential.

[0075] S4. Calculate the extreme point spacing function based on the instantaneous frequency and construct the resampling function. Obtain the resampling sequence by cubic spline interpolation. Obtain the spectrum of the resampling sequence by discrete Fourier transform of the resampling sequence. Convolve the double average filter function to construct the filter loop sequence. Obtain the filter spectrum by discrete Fourier transform of the filter loop sequence.

[0076] S5. Multiply the corresponding terms of the resampled sequence spectrum and the filter spectrum one by one and iterate until the convergence condition is met. Then, perform an inverse Fourier transform to obtain the first-order IMF component. Remove the first-order IMF component from the signal to be processed and repeat S2 to S5 to extract each order of IMF component until the remaining signal has no more than two extreme points or reaches the maximum decomposition order.

[0077] S6. Regularize the IMF components and derivative sequences of each order to construct a complex sequence and obtain the instantaneous phase;

[0078] S7. After shaping the instantaneous phase into a monotonically increasing form, the forward differential is used to obtain the non-negative instantaneous frequency. The time spectrum is constructed based on the time and the non-negative instantaneous frequency.

[0079] S8, and the remaining three magnetotelluric signals are processed from S2 to S7 to obtain the time spectrum.

[0080] Figure 2 This is a schematic diagram of the technical route provided in the embodiments of the present invention, such as... Figure 2 As shown, the implementation steps of the technical approach in this embodiment are as follows:

[0081] Step 1: Import the magnetotelluric signal from the instrument data card into the computer. Based on the difference in the direction of the four recorded electromagnetic signals, they are denoted as Ex (north-south direction), Ey (east-west direction), Hx (north-south direction), and Hy (east-west direction).

[0082] Step 2: Take one of Ex, Ey, Hx, and Hy and denote it as x(t). Based on x(t), define a local mean operator L. The functional expression of the local mean operator L is L(x)(t):

[0083] (1)

[0084] In the above formula, N is the length of the signal x(t), τ is an intermediate variable, and a is a low-pass filter function. Different values ​​of a low-pass filter function can affect the decomposition result and the decomposition rate. In this embodiment, considering both the reliability of the decomposition result and the ease of implementation, a is chosen as a double-average filter function.

[0085] (2)

[0086] The functional expression of the sieving operator S can be defined based on the functional expression of the local mean operator L:

[0087] (3)

[0088] Repeat the screening process described above, that is:

[0089] (4)

[0090] In the above formula, x m The result obtained by filtering signal x through m times (x m (t) is the functional expression of the sieving result. The functional expressions of subsequent parameters are similar to the parameter explanations here, and will not be repeated here. When m reaches the maximum number of iterations m max (Set according to signal complexity; typically set to 500 in magnetotelluric signal processing), or stop screening when the convergence parameter SD is less than or equal to a given threshold. SD is defined as:

[0091] (5)

[0092] In the above formula, ||·||2 represents the 2-norm, and in this embodiment, the threshold is set to 0.05. The value of x after stopping the screening is taken. m As the first-order temporary IMF component of signal x, it is denoted as I1(t).

[0093] Step 3: Perform Hilbert transform on component I1(t) (In this embodiment, Hilbert transform is used to process the temporary IMF component because the instantaneous frequency obtained in this step is only for solving the spacing function, not the final result of the instantaneous spectrum, and the requirements for time resolution and physical meaning are not high). The formula is as follows:

[0094] (6)

[0095] Take the sequences before and after the transformation, and construct the complex sequence I1. * (t), the formula is as follows:

[0096] (7)

[0097] In the above formula, i is the imaginary unit, then the phase of the complex sequence is:

[0098] (8)

[0099] Using forward difference pairs of complex sequences I1 * The phase derivative of (t) yields:

[0100] (9)

[0101] In the above formula, ω1(t) is the instantaneous frequency. Since forward differencing reduces the sequence length, this embodiment takes the same value for the instantaneous frequency at t=N-1 as the instantaneous frequency of its preceding term. Exemplarily, steps 2 and 3 of this embodiment are processes of obtaining the time-frequency representation of the signal using iterative filtering. Their function is to obtain the change of the highest frequency component of the signal over time, and to calculate the resampling interval function accordingly.

[0102] Step 4: Based on the results obtained above, perform resampling processing to obtain the resampling sequence h and the filtering function f. First, based on ω1(t), the spacing function l(t) of the extreme points can be obtained (in this embodiment, the spacing function is obtained from the instantaneous frequency of the temporary IMF component, which is a value that changes continuously with time, rather than a traditional fixed length, and its change is consistent with the highest frequency component of the current signal. This ensures that the sampling spacing can closely follow the signal change in each iteration of resampling, avoiding unnecessary additional computation due to oversampling, avoiding loss of high-frequency signals, and having good adaptability). The formula is:

[0103] (10)

[0104] In the above formula, ξ0 is the spacing factor. In this embodiment, ξ0 = 1 is taken to effectively reflect high-frequency changes. Based on l(t), the resampling function G can be obtained. -1 (t):

[0105] (11)

[0106] Take the filter parameter M=G -1 (1) The resampling function h(x) can be obtained by cubic spline interpolation of x(t) on the point set G(y), where y=Mt / N is the normalized grid and G(y) is the G -1 The reciprocal of (y), that is, for any y>0, we have:

[0107] (12)

[0108] Define the value of h1 in the first iteration of the resampled sequence as h(G(y)), and perform a discrete Fourier transform on it to obtain its spectrum:

[0109] (13)

[0110] In the above formula, k is the frequency grid, with values ​​of 0, 1, 2, ..., N-1; Y is the maximum value M(N-1) / N that the standard grid y can take; and H1 is the discrete Fourier spectrum of h1. The filter function f is taken as a reconvolution of the double-average filter function a. (The reconvolution filter used in this embodiment has advantages such as convergence, positive definiteness, and easy solution of eigenvalues ​​during the calculation process.) That is:

[0111] (14)

[0112] In the above formula, s is the support length of the filter function f, usually taken as s = [N / M], where [·] is the floor operator. Based on this, the filter cyclic sequence f is defined. s The formula is:

[0113] (15)

[0114] In the above formula, f s The length is N, P=M / N. ∑f(x(t)-x(τ)) is the re-regularization constant, which is approximated as P=1 in this embodiment. Based on this, the filter spectrum is defined as F. s The formula is:

[0115] (16)

[0116] Step 5: Based on the results obtained in Step 4, the second iteration value of the resampled sequence spectrum is:

[0117] (17)

[0118] In the above formula, "·" indicates multiplication of corresponding terms. Repeating the above operation of multiplying corresponding terms...

[0119] (18)

[0120] Obtain H2, H3, ..., H in sequence m Until:

[0121] or (19)

[0122] In the above formula, m is the current iteration number, m max The maximum number of iterations is given, and is usually taken to be greater than 1000 to ensure sufficient convergence. δ is a given tolerance error, usually 0.05. The spectrum H of the resampling function is taken when the above conditions are met. m The inverse Fourier transform of α1, which is the first-order IMF component, is given by the following formula:

[0123] (20)

[0124] Step 6: Remove the first-order IMF component from the signal x(t) to obtain the remaining signal r, as shown in the formula:

[0125] (twenty one)

[0126] Using r instead of x(t) in step 2, the calculation process from step 2 to step 5 is performed. That is, the obtained IMF components are removed from the signal x(t), and the samples are substituted into the above steps for repeated sampling and iterative filtering to extract the IMF components α2, α3, ..., α1 one by one. n When the extreme point of the residual r does not exceed 2 or the number of iterations n>n max The iteration stops when n is reached. max Given the maximum number of iterations, typically an integer between 10 and 20, the expression for r is:

[0127] (twenty two)

[0128] In the above formula, R represents the total number of IMF components currently obtained.

[0129] Step 7: Take the IMF component obtained in the previous step, and differentiate it using forward difference. The formula is:

[0130] (twenty three)

[0131] Then, the IMF component α was analyzed. n and its derivative sequence α n Perform regularization processing separately:

[0132] (twenty four)

[0133] (25)

[0134] In the above formula, F 1n F 2n α n With α n The regularized sequence has q(t) and p(t) as α. n With α n An envelope of ' (connected by a smooth curve α) n With α n (The maximum or minimum value of ′ is obtained).

[0135] Step 8: Based on the sequence F obtained in Step 5 1n F 2n Construct complex sequence F n

[0136] (26)

[0137] In the above formula, i is the imaginary unit, and it is based on the obtained sequence F. n The instantaneous phase of signal x at time t can be obtained as follows (in this embodiment, the complex sequence is constructed using a sequence and its derivative sequence regularization method, avoiding the use of the integral form of Hilbert transform, thus resulting in higher phase and frequency time resolution):

[0138] (27)

[0139] The instantaneous phase extrema points obtained above are numbered in chronological order. When the first extremum point is a maximum value, the phase of the j-th point is:

[0140] (28)

[0141] In the above formula, j is the number of the extreme point, and Φ(t) is the transformed phase. When the first extreme point is a minimum, we have:

[0142] (29)

[0143] The instantaneous frequency ω can be obtained by performing forward differential, as shown in the formula:

[0144] (30)

[0145] The instantaneous amplitude is α n The absolute value of (t) (Although the instantaneous frequency is obtained by differentiation in this embodiment, the result obtained after phase waveform shaping transformation is a monotonically increasing function, and there will be no negative frequency without physical meaning. At the same time, the phase is transformed in units of π, and the result will not change the actual value of the phase).

[0146] Step 9: Plot time t on the horizontal axis and instantaneous frequency ω on the vertical axis, using color depth to represent the magnitude of instantaneous amplitude or instantaneous phase, thereby obtaining the amplitude-frequency and phase-frequency information of the IMF. Map the amplitude-frequency and phase-frequency information of IMFs of different orders onto the same coordinate axis using the same representation method described above, thus obtaining the instantaneous amplitude spectrum and instantaneous phase spectrum (i.e., the time spectrum of x) of signal x, respectively.

[0147] Step 10: Process the other magnetotelluric signals according to the procedures in steps 2 to 9, and then organize them to obtain the time spectrum of the magnetotelluric signal.

[0148] In a specific implementation, this method decomposes magnetotelluric signals using an improved fast resampling iterative filter. Unlike recursive decomposition relying on envelope fitting, this method obtains the moving average of the signal through a low-pass filter and extracts intrinsic mode function components within the iterative filtering framework. This provides a clearer mathematical basis for the decomposition process, improving the physical meaning and convergence stability of the decomposed components. Furthermore, this method introduces an adaptive mechanism for the resampling interval, allowing the filtering and resampling processes to adjust to local signal variations. This enhances the adaptability to complex, non-stationary signals and reduces errors caused by the ambiguity and instability of the decomposition results.

[0149] In its specific implementation, this method employs a strategy combining phase waveform shaping and phase differentiation for instantaneous spectrum estimation. First, a complex sequence is constructed based on the components, and the phase is calculated. Then, without altering the rate of phase change, the phase waveform is shaped to be monotonically increasing, thus avoiding the appearance of physically meaningless negative frequencies. Subsequently, forward differential is used to differentiate the phase to obtain the instantaneous frequency, and the reduction in sequence length caused by the final differential is padded to ensure the continuous usability of the instantaneous frequency sequence. Compared to methods that rely solely on integrals to obtain average results, this strategy obtains the instantaneous frequency change at the corresponding moment, which is more conducive to reflecting the transient characteristics of the signal and improving time resolution and interpretability.

[0150] In its specific implementation, this method replaces the envelope fitting step in traditional recursive decomposition with iterative filtering and employs a reconvolution filter for spectral iteration updates. This effectively suppresses mode aliasing and reduces the dispersion of the same frequency component across components or the presence of multiple frequency components in a single component, thereby improving the frequency resolution and spectral reliability of the instantaneous spectrum. Furthermore, performing the same process on four magnetotelluric signal data channels yields consistent time-spectral characterization, providing a more stable and reliable time-frequency basis for subsequent denoising, anomaly identification, and impedance estimation. This approach is particularly suitable for real-world observation scenarios with high noise levels and significant non-stationary characteristics.

[0151] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0152] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method of magnetotelluric signal instantaneous spectral estimation, characterized in that, include: S1. Acquire four magnetotelluric signals and select one of them as the signal to be processed; S2. Define a local mean operator and a sieving operator using a double average filter function, and repeatedly sieve the signal to be processed to obtain the first-order temporary IMF component. S3. The phase is obtained by Hilbert transformation of the first-order temporary IMF component, and the instantaneous frequency is obtained by forward differential. S4. Calculate the extreme point spacing function based on the instantaneous frequency and construct the resampling function. Obtain the resampling sequence by cubic spline interpolation. Obtain the spectrum of the resampling sequence by discrete Fourier transform of the resampling sequence. Convolve the double average filter function to construct the filter loop sequence. Obtain the filter spectrum by discrete Fourier transform of the filter loop sequence. S5. Multiply the corresponding terms of the resampled sequence spectrum and the filter spectrum one by one and iterate until the convergence condition is met. Then, perform an inverse Fourier transform to obtain the first-order IMF component. Remove the first-order IMF component from the signal to be processed and repeat S2 to S5 to extract each order of IMF component until the remaining signal has no more than two extreme points or reaches the maximum decomposition order. S6. Regularize the IMF components and derivative sequences of each order to construct a complex sequence and obtain the instantaneous phase; S7. After shaping the instantaneous phase into a monotonically increasing form, the forward differential is used to obtain the non-negative instantaneous frequency. The time spectrum is constructed based on the time and the non-negative instantaneous frequency. S8, the remaining three magnetotelluric signals are processed from S2 to S7 to obtain the time spectrum; The double-average filter function satisfies: where N is the length of the signal to be processed.

2. The magnetotelluric signal time -frequency spectrum estimation method according to claim 1, characterized in that, Four magnetotelluric signal data were acquired, including: Acquire electric field data in the north-south direction and electric field data in the east-west direction; Acquire magnetic track data in the north-south direction and the east-west direction of the magnetic field; The electrical track data and the magnetic track data are used together as the magnetotelluric signal data.

3. The MT signal time spectrum estimation method according to claim 1, characterized in that, Step S2 includes: S21. Set the signal to be processed as the screening input sequence for this round of screening; S22. The local mean operator is used to perform low-pass filtering on the sieving input sequence to obtain a local mean sequence; S23. Subtract the local mean sequence from the screening input sequence to obtain the screening result for this round; S24. Update the current round of screening results to the screening input sequence for the next round of screening, and repeat steps S22 to S24 until the maximum number of iterations is reached or the convergence parameter meets the threshold, to obtain the first-order temporary IMF component; wherein, the maximum number of iterations is 500; the convergence parameter is the ratio of the L2 norm of the difference between two adjacent rounds of screening results to the L2 norm of the next screening result; the threshold is 0.

05.

4. The magnetotelluric signal time -frequency spectrum estimation method according to claim 1, characterized in that, Step S3 includes: Perform a Hilbert transform on the first-order temporary IMF component to obtain a Hilbert transform sequence; Construct a complex sequence using the first-order temporary IMF component as the real part and the Hilbert transform sequence as the imaginary part; The phase is determined based on the complex sequence. The instantaneous frequency is obtained by differentiating the phase using forward differential, and the last term of the instantaneous frequency is set to the same value as the previous term of the instantaneous frequency.

5. The instantaneous spectrum estimation method for magnetotelluric signals according to claim 1, characterized in that, Step S4 includes: The extreme point spacing function is calculated based on the instantaneous frequency; A resampling function is constructed based on the extreme point spacing function, and the resampling point set is determined by the resampling function; The signal to be processed is subjected to cubic spline interpolation on the set of resampled points to obtain a resampled sequence; Perform a discrete Fourier transform on the resampled sequence to obtain the spectrum of the resampled sequence; The filter function is taken as the reconvolution of the double average filter function and a filter loop sequence is constructed. The filter loop sequence is then subjected to a discrete Fourier transform to obtain the filter spectrum.

6. The instantaneous spectrum estimation method for magnetotelluric signals according to claim 1, characterized in that, The resampled sequence spectrum and the corresponding terms of the filter spectrum are multiplied term by term and iterated. After the convergence condition is met, the inverse Fourier transform is used to obtain the first-order IMF components, including: S51. Use the spectrum of the resampled sequence as the spectrum of the current iteration; S52. Multiply the current iteration spectrum with the filter spectrum by corresponding terms to obtain the updated iteration spectrum; S53. Determine whether the difference between the updated iteration spectrum and the current iteration spectrum meets the allowable error. If not, repeat step S52 using the updated iteration spectrum as the new current iteration spectrum. The convergence condition includes that the difference between the updated iteration spectrum and the current iteration spectrum meets the allowable error. S54. When the tolerance error is met or the iteration limit is reached, perform an inverse Fourier transform on the iterative spectrum that meets the conditions to obtain the first-order IMF component; wherein, the iteration limit is 1000; and the tolerance error is 0.

05.

7. The magnetotelluric signal time -frequency spectrum estimation method according to claim 1, characterized in that, The first-order IMF component is removed from the signal to be processed, and S2 to S5 are repeated to extract IMF components of each order until the remaining signal has no more than two extreme points or reaches the maximum decomposition order, including: When the number of remaining signal extreme points does not exceed two, stop extracting the IMF components of each order one by one; When the extraction order of each IMF component exceeds the maximum decomposition order, the extraction of each IMF component is stopped; wherein, the maximum decomposition order is an integer from ten to twenty.

8. The instantaneous spectrum estimation method for magnetotelluric signals according to claim 1, characterized in that, Step S6 includes: For each IMF component, a smooth curve is used to connect the maximum or minimum points of each IMF component to form the first envelope. Based on the first envelope, each order of IMF component is regularized to obtain the regularized IMF component. For the derivative sequence corresponding to each IMF component, a smooth curve is used to connect the maximum or minimum points of the derivative sequence to form a second envelope; The derivative sequence is regularized based on the second envelope to obtain a regularized derivative sequence. A complex sequence is constructed using the regularized IMF components and the regularized derivative sequence, and the instantaneous phase corresponding to each order of IMF components is obtained from the complex sequence.

9. The instantaneous spectrum estimation method for magnetotelluric signals according to claim 1, characterized in that, Step S7 includes: The extreme points of the instantaneous phases corresponding to each order of IMF components are numbered in chronological order. When the extreme point numbered one is a maximum value, the instantaneous phase is shaped according to the first phase shaping rule; When the extreme point numbered one is a minimum value, the instantaneous phase is shaped according to the second phase shaping rule so that the shaped instantaneous phase is monotonically increasing. The non-negative instantaneous frequency is obtained by using forward differential differentiation on the shaped monotonically increasing instantaneous phase. A time spectrum is constructed with time as the horizontal axis, the non-negative instantaneous frequency as the vertical axis, and the instantaneous amplitude or instantaneous phase as the spectral value; wherein the instantaneous amplitude is the absolute value of the IMF component of the corresponding order.

10. The instantaneous spectrum estimation method for magnetotelluric signals according to claim 5, characterized in that, The filter function is a reconvolution of the double-average filter function, and a filter loop sequence is constructed. The filter loop sequence is then subjected to a discrete Fourier transform to obtain the filter spectrum, including: The filter function is re-regularized to obtain a re-regularized filter function, wherein the re-regularization constant used for the re-regularization process is approximated to one. A filter loop sequence is constructed based on the re-regularized filter function, and the length of the filter loop sequence is set to be the same as the length of the resampling sequence. The filter spectrum is obtained by performing a discrete Fourier transform on the filtered cyclic sequence.

Images