A method for estimating magnetotelluric impedance

Abstract

The disclosed embodiments provide a method for estimating magnetotelluric impedance, which belongs to the field of surveying and specifically includes: Step 1: Performing a Fourier transform on electromagnetic time series data of multiple electromagnetic field components to obtain a spectrum array corresponding to each electromagnetic field component; Step 2: Identifying the interference frequency of each electromagnetic field component based on the spectrum array; Step 3: Combining the interference frequency and deviation value; Step 4: Attenuating the spectrum array using the combined interference frequency and deviation value; Step 5: Restoring the electromagnetic time series data using the attenuated spectrum array to obtain a transformation result; Step 6: Estimating magnetotelluric impedance based on the transformation result. The disclosed solution improves the adaptability and accuracy of magnetotelluric impedance estimation.

Description

Technical Field

[0001] The embodiments of the present disclosure relate to the field of surveying, and in particular to a method for estimating magnetotelluric impedance. Background Art

[0002] Magnetotellurics (MT) is an exploration method that uses natural electromagnetic signals to detect underground targets. It is widely used in mineral exploration, oil and gas exploration, engineering exploration, and disaster monitoring. MT offers advantages such as deep exploration, convenient operation, and high exploration efficiency. However, it is susceptible to electromagnetic interference, particularly in areas with frequent human activities. To address the increasing severity of electromagnetic interference, researchers at home and abroad have adopted various noise reduction methods. Gamble et al. first proposed and utilized a remote reference denoising method, which can suppress irrelevant noise in acquired signals. Egbert et al. used an improved weighted least squares algorithm for impedance estimation and compared the performance of related algorithms. Tang Jingtian et al. used mathematical morphology filtering to denoise time series, removing strong interference noise with a well-defined pattern. Garcia et al. used the continuous wavelet transform to select data from AMT dead bands, addressing the problem of weak signals in these dead bands. Tang Jingtian et al. introduced the Hilbert-Huang transform into MT data processing, using empirical mode decomposition to remove some low-frequency interference. These noise removal methods can achieve some results for certain data types, but the processing effect is not ideal for strong interference data with long duration, strong interference amplitude and obvious periodicity.

[0003] It can be seen that there is an urgent need for a method for estimating the magnetotelluric impedance with high detection adaptability and accuracy. Summary of the Invention

[0004] In view of this, an embodiment of the present disclosure provides a method for estimating magnetotelluric impedance, which at least partially solves the problems of poor detection adaptability and accuracy in the prior art.

[0005] The present disclosure provides a method for estimating magnetotelluric impedance, including:

[0006] Step 1: Perform Fourier transform on electromagnetic time series data of multiple electromagnetic field components to obtain a spectrum array corresponding to each electromagnetic field component;

[0007] Step 2, identifying the interference frequency of each electromagnetic field component based on the spectrum array;

[0008] Step 3, combine the interference frequency and deviation value;

[0009] Step 4, attenuate the spectrum array using the combined interference frequency and deviation value;

[0010] Step 5, using the attenuated spectrum array to restore the electromagnetic time series data to obtain the transformation result;

[0011] Step 6: Estimating the magnetotelluric impedance based on the transformation results.

[0012] According to a specific implementation of the embodiment of the present disclosure, the electromagnetic field components are at least two of Ex, Ey, Hx, Hy, and Hz, where Ex is the electric field signal component in the X direction, Ey is the electric field signal component in the Y direction, Hx is the magnetic field signal component in the X direction, Hy is the magnetic field signal component in the Y direction, and Hz is the magnetic field signal component in the Z direction. The X direction and the Y direction are horizontal directions and orthogonal to each other, and the Z direction is a vertical direction.

[0013] The electromagnetic time series data is sample point data of each electromagnetic field component collected at a preset sampling rate.

[0014] According to a specific implementation of the embodiment of the present disclosure, step 2 specifically includes:

[0015] Step 2.1, calculate the power spectrum array Pi(w) of each electromagnetic field component based on the spectrum array;

[0016] Step 2.2, filter Pi(w) using the median filter method to obtain filtered data Pi1(w);

[0017] Step 2.3, calculate the difference ΔPi(w) between Pi(w) and Pi1(w);

[0018] Step 2.4: Select all interference frequency arrays Wr that are greater than the threshold Y0 from ΔPi(w), and record the ΔPi(w) corresponding to Wr into the deviation value Dr.

[0019] According to a specific implementation of the embodiment of the present disclosure, step 3 specifically includes:

[0020] Step 3.1, merge the interference frequency array Wr of each electromagnetic field component into an array Ws by taking the union method, and only retain one identical frequency value in Ws;

[0021] In step 3.2, the composite deviation value of each frequency is calculated according to the order of the frequencies in the Ws array, and stored in the array Ds in order.

[0022] According to a specific implementation of the embodiment of the present disclosure, the composite deviation value corresponds to the maximum value of the deviation values ​​of all components of the frequency.

[0023] According to a specific implementation of the embodiment of the present disclosure, step 4 specifically includes:

[0024] Step 4.1, obtain the number of elements N of the interference frequency Ws and initialize i=1;

[0025] Step 4.2, obtain the i-th frequency Wi in Ws and the i-th deviation value Di in Ds;

[0026] Step 4.3, calculate the attenuation coefficient Ri using the deviation value Di;

[0027] Step 4.4, attenuate the spectrum value Ai corresponding to each electromagnetic field component frequency Wi of the electromagnetic field, that is, Ai = Ai*Ri;

[0028] Step 4.5: Determine the value of i. If i is equal to N, complete the attenuation of each spectrum array. If i is not equal to N, set i = i + 1 and jump to step 4.2 to continue processing.

[0029] According to a specific implementation of the embodiment of the present disclosure, the attenuation frequency value of each electromagnetic field component is the same, and the attenuation coefficient of each electromagnetic field component with the same frequency is the same.

[0030] According to a specific implementation of the embodiment of the present disclosure, Ri and Di are inversely correlated, and the calculation formula of Ri is Ri=10 -Di .

[0031] The magnetotelluric impedance estimation scheme in the embodiment of the present disclosure includes: step 1, performing Fourier transform on electromagnetic time series data of multiple electromagnetic field components to obtain a spectrum array corresponding to each electromagnetic field component; step 2, identifying the interference frequency of each electromagnetic field component based on the spectrum array; step 3, combining the interference frequency and deviation value; step 4, attenuating the spectrum array using the combined interference frequency and deviation value; step 5, restoring the electromagnetic time series data using the attenuated spectrum array to obtain a transformation result; and step 6, estimating the magnetotelluric impedance based on the transformation result.

[0032] The beneficial effects of the embodiments of the present disclosure are:

[0033] 1. It can automatically identify strong interference signals in magnetotelluric signals, eliminating the need for manual analysis and processing of massive data and improving data processing efficiency;

[0034] 2. It can suppress the periodic strong interference signal in the magnetotelluric signal and improve the accuracy of magnetotelluric impedance estimation;

[0035] 3. Using the anti-correlation algorithm between the frequency attenuation coefficient and the deviation value, the greater the interference, the stronger the suppression, achieving the effect of adaptive denoising;

[0036] 4. The multi-component same-frequency and same-coefficient attenuation method ensures that the calculation results will not be distorted. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] In order to more clearly illustrate the technical solutions of the embodiments of the present disclosure, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present disclosure. For ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0038] Figure 1 A schematic diagram of a flow chart of a method for estimating magnetotelluric impedance provided by an embodiment of the present disclosure;

[0039] Figure 2 A schematic diagram of a process for combining four component interference frequencies and deviation values ​​provided in an embodiment of the present disclosure;

[0040] Figure 3 A schematic diagram of a process for attenuating a single component spectrum provided in an embodiment of the present disclosure. DETAILED DESCRIPTION

[0041] The embodiments of the present disclosure are described in detail below with reference to the accompanying drawings.

[0042] The following describes the embodiments of the present disclosure through specific examples, and those skilled in the art can easily understand other advantages and effects of the present disclosure from the contents disclosed in this specification. Obviously, the described embodiments are only a part of the embodiments of the present disclosure, rather than all of the embodiments. The present disclosure can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present disclosure. It should be noted that, in the absence of conflict, the following embodiments and features in the embodiments can be combined with each other. Based on the embodiments in the present disclosure, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present disclosure.

[0043] It should be noted that various aspects of the embodiments within the scope of the appended claims are described below. It should be apparent that the aspects described herein can be embodied in a wide variety of forms, and any specific structure and / or function described herein is merely illustrative. Based on this disclosure, it should be understood by those skilled in the art that an aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects described herein can be used to implement an apparatus and / or practice a method. In addition, other structures and / or functionalities other than one or more of the aspects described herein can be used to implement this apparatus and / or practice this method.

[0044] It should also be noted that the illustrations provided in the following embodiments are only schematic illustrations of the basic concept of the present disclosure. The illustrations only show components related to the present disclosure and are not drawn according to the number, shape and size of components in actual implementation. In actual implementation, the type, quantity and proportion of each component can be changed at will, and the component layout type may also be more complicated.

[0045] Additionally, in the following description, specific details are provided to provide a thorough understanding of the examples. However, one skilled in the art will appreciate that the aspects described can be practiced without these specific details.

[0046] The embodiments of the present disclosure provide a method for estimating magnetotelluric impedance, which can be applied to the process of detecting underground targets using natural electromagnetic signals in scenarios such as mineral exploration, oil and gas exploration, engineering exploration, and disaster monitoring.

[0047] See also Figure 1 , is a flow chart of a method for estimating magnetotelluric impedance provided by an embodiment of the present disclosure. Figure 1 As shown, the method mainly includes the following steps:

[0048] Step 1: Perform Fourier transform on electromagnetic time series data of multiple electromagnetic field components to obtain a spectrum array corresponding to each electromagnetic field component;

[0049] Optionally, the electromagnetic field components are at least two of Ex, Ey, Hx, Hy, and Hz, where Ex is the electric field signal component in the X direction, Ey is the electric field signal component in the Y direction, Hx is the magnetic field signal component in the X direction, Hy is the magnetic field signal component in the Y direction, and Hz is the magnetic field signal component in the Z direction, the X direction and the Y direction are horizontal directions and orthogonal to each other, and the Z direction is a vertical direction;

[0050] The electromagnetic time series data is sample point data of each electromagnetic field component collected at a preset sampling rate.

[0051] For example, the electromagnetic field components of the magnetotelluric method that need to be processed are Ex, Ey, Hx, and Hy, a total of 4 components, and their time series data are X1, X2, X3, and X4, the time series length is N, and the sampling rate is S. Then, X1, X2, X3, and X4 can be subjected to Fourier transform of length N to obtain the spectrum arrays Y1, Y2, Y3, and Y4 of the 4 electromagnetic field components.

[0052] Step 2, identifying the interference frequency of each electromagnetic field component based on the spectrum array;

[0053] Furthermore, the step 2 specifically includes:

[0054] Step 2.1, calculate the power spectrum array Pi(w) of each electromagnetic field component based on the spectrum array;

[0055] Step 2.2, filter Pi(w) using the median filter method to obtain filtered data Pi1(w);

[0056] Step 2.3, calculate the difference ΔPi(w) between Pi(w) and Pi1(w);

[0057] Step 2.4: Select all interference frequency arrays Wr that are greater than the threshold Y0 from ΔPi(w), and record the ΔPi(w) corresponding to Wr into the deviation value Dr.

[0058] In specific implementation, the obtained spectrum array can be used to calculate the power spectrum arrays P1, P2, P3, and P4 of the four electromagnetic field components. The formula for converting the spectrum value Yi of the i-th frequency into the power spectrum value Pi is P i =20×log 10 Y i , copy the power spectrum arrays P1, P2, P3, and P4 to the four newly applied power spectrum arrays P11, P21, P31, and P41, and perform median filtering on P11, P21, P31, and P41 in a distributed manner, with a filter length of FL = N / 30. Calculate the difference between the unfiltered power spectrum and the filtered power spectrum, and store them in the difference arrays ΔP1, ΔP2, ΔP3, and ΔP4 respectively. Then, select the element values ​​​​greater than the threshold Y0 from the ΔP1, ΔP2, ΔP3, and ΔP4 arrays and store them in the Dr[i] array (i = 0 to 3), and store the index values ​​of the elements in the difference array in the Wr[i] array (i = 0 to 3). Preferably, Y0 = 10 is generally taken.

[0059] Step 3, combine the interference frequency and deviation value;

[0060] Furthermore, the step 3 specifically includes:

[0061] Step 3.1, merge the interference frequency array Wr of each electromagnetic field component into an array Ws by taking the union method, and only retain one identical frequency value in Ws;

[0062] In step 3.2, the composite deviation value of each frequency is calculated according to the order of the frequencies in the Ws array, and stored in the array Ds in order.

[0063] Optionally, the composite deviation value corresponds to the maximum value of the deviation values ​​of all components of the frequency.

[0064] When implementing it specifically, Figure 2As shown, the four Wr[i] arrays can be merged into the array Ws, and the four Dr[i] arrays can be merged into the array Ds to complete the merging of the interference frequencies and deviation values. The final merged frequency array and deviation value array are Ws and Ds, respectively. Wr and Dr in the figure are two-dimensional arrays, storing the interference frequency values ​​and deviation values ​​of the four components, respectively. "Get the index k of Wr[i][j] in Ws" means searching for the value of Wr[i][j] in the Ws array. If the value is found, the index of the value in Ws is returned; otherwise, -1 is returned.

[0065] Step 4, attenuate the spectrum array using the combined interference frequency and deviation value;

[0066] Furthermore, the step 4 specifically includes:

[0067] Step 4.1, obtain the number of elements N of the interference frequency Ws and initialize i=1;

[0068] Step 4.2, obtain the i-th frequency Wi in Ws and the i-th deviation value Di in Ds;

[0069] Step 4.3, calculate the attenuation coefficient Ri using the deviation value Di;

[0070] Step 4.4, attenuate the spectrum value Ai corresponding to each electromagnetic field component frequency Wi of the electromagnetic field, that is, Ai = Ai*Ri;

[0071] Step 4.5: Determine the value of i. If i is equal to N, complete the attenuation of each spectrum array. If i is not equal to N, set i = i + 1 and jump to step 4.2 to continue processing.

[0072] Furthermore, the attenuation frequency value of each electromagnetic field component is the same, and the attenuation coefficient of each electromagnetic field component with the same frequency is the same.

[0073] Furthermore, Ri and Di are inversely correlated, and the calculation formula for Ri is Ri=10 -Di .

[0074] In specific implementation, the values ​​in Ws and Ds and the specific attenuation process can be used to process the spectrum arrays Y1, Y2, Y3, and Y4 until the attenuation is completed. The specific process of attenuating a single component spectrum is as follows: Figure 3 As shown, Y in the figure represents the spectrum array pointer of a certain component. For example, if the second component is processed, Y=Y2.

[0075] Step 5, using the attenuated spectrum array to restore the electromagnetic time series data to obtain the transformation result;

[0076] During specific implementation, an inverse Fourier transform of length N may be performed on the spectrum arrays Y1, Y2, Y3, and Y4, and the transform results may be stored in X1, X2, X3, and X4, respectively.

[0077] Step 6: Estimating the magnetotelluric impedance based on the transformation results.

[0078] In specific implementation, conventional magnetotelluric impedance estimation can be performed based on the transformation results to obtain parameters such as magnetotelluric impedance and apparent resistivity.

[0079] The magnetotelluric impedance estimation method provided in this embodiment automatically identifies all interference frequencies by performing power spectrum analysis on each component of the electromagnetic signal, calculates the deviation value of the interference frequency, and then synchronously attenuates the spectrum of each component at the interference frequency according to the magnitude of the deviation value to obtain a signal after filtering out interference. Finally, impedance estimation processing is performed, thereby improving the adaptability and accuracy of magnetotelluric impedance estimation.

[0080] The units involved in the embodiments described in this disclosure may be implemented by software or hardware.

[0081] It should be understood that various parts of the present disclosure can be implemented in hardware, software, firmware, or a combination thereof.

[0082] The above description is merely a specific embodiment of the present disclosure, but the scope of protection of the present disclosure is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in this disclosure should be included in the scope of protection of the present disclosure. Therefore, the scope of protection of the present disclosure should be based on the scope of protection of the claims.

Claims

1. A method for estimating magnetotelluric impedance, characterized in that: include: Step 1: Perform Fourier transform on electromagnetic time series data of multiple electromagnetic field components to obtain a spectrum array corresponding to each electromagnetic field component; Step 2, identifying the interference frequency of each electromagnetic field component based on the spectrum array; Step 3, combine the interference frequency and deviation value; Step 4, attenuate the spectrum array using the combined interference frequency and deviation value; Step 5, using the attenuated spectrum array to restore the electromagnetic time series data to obtain the transformation result; Step 6: Estimating the magnetotelluric impedance based on the transformation results.

2. The method according to claim 1, characterized in that , the electromagnetic field components are at least two of Ex, Ey, Hx, Hy, and Hz, where Ex is the electric field signal component in the X direction, Ey is the electric field signal component in the Y direction, Hx is the magnetic field signal component in the X direction, Hy is the magnetic field signal component in the Y direction, and Hz is the magnetic field signal component in the Z direction, the X direction and the Y direction are horizontal directions and orthogonal to each other, and the Z direction is the vertical direction; The electromagnetic time series data is sample point data of each electromagnetic field component collected at a preset sampling rate.

3. The method according to claim 1, characterized in that , the step 2 specifically includes: Step 2.1, calculate the power spectrum array Pi(w) of each electromagnetic field component based on the spectrum array; Step 2.2, filter Pi(w) using the median filter method to obtain filtered data Pi1(w); Step 2.3, calculate the difference ΔPi(w) between Pi(w) and Pi1(w); Step 2.4: Select all interference frequency arrays Wr that are greater than the threshold Y0 from ΔPi(w), and record the ΔPi(w) corresponding to Wr into the deviation value Dr.

4. The method according to claim 1, characterized in that , the step 3 specifically includes: Step 3.1, merge the interference frequency array Wr of each electromagnetic field component into an array Ws by taking the union method, and only retain one identical frequency value in Ws; In step 3.2, the composite deviation value of each frequency is calculated according to the order of the frequencies in the Ws array, and stored in the array Ds in order.

5. The method according to claim 4, characterized in that , the composite deviation value corresponds to the maximum value of the deviation values ​​of all components of the frequency.

6. The method according to claim 1, characterized in that , the step 4 specifically includes: Step 4.1, obtain the number of elements N of the interference frequency Ws and initialize i=1; Step 4.2, obtain the i-th frequency Wi in Ws and the i-th deviation value Di in Ds; Step 4.3, calculate the attenuation coefficient Ri using the deviation value Di; Step 4.4, attenuate the spectrum value Ai corresponding to each electromagnetic field component frequency Wi of the electromagnetic field, that is, Ai = Ai*Ri; Step 4.5: Determine the value of i. If i is equal to N, complete the attenuation of each spectrum array. If i is not equal to N, set i = i + 1 and jump to step 4.2 to continue processing.

7. The method according to claim 6, characterized in that ,The attenuation frequency value of each electromagnetic field component is the same, and the attenuation coefficient of the same frequency of each electromagnetic field component is the same.

8. The method according to claim 6, characterized in that ,Ri and Di are inversely correlated, and the calculation formula of Ri is Ri=10-Di.

Images