Static Shift Correction and Fractal Characteristic Analysis of Time-Frequency Electromagnetic

1. Introduction

In electromagnetic exploration, the static shift effect is a phenomenon in which the presence of electrically heterogeneous bodies near the surface causes the accumulation of static charges on their surfaces due to the horizontal electric field, distorting the amplitude of the measured electric field and leading to a shift in apparent resistivity with logarithmic frequency, severely affecting the effectiveness of exploration [8,10,13].

To mitigate this effect, several methods have been developed, primarily falling into three categories: First, correction methods that incorporate prior information. For example, statistical variations in profile resistivity can be obtained through outcrop geology and resistivity measurements in the study area, or through statistical analysis of drilling data [8] and high-quality resistivity sounding curves [4,10,12,15]. However, such prior information is difficult to match perfectly with electromagnetic exploration areas and can be costly. Second, low-pass filtering methods. Since static shifts exhibit high-frequency characteristics, low-pass filtering can be applied to suppress such effects. The electromagnetic array profiling method, originally proposed by Bostick [4], uses spatial filtering for correction. Torres-Verdín and Bostick [17] analyzed the secondary magnetic field and electric field anomalies of shallow subsurface anomalies using the Born approximation method. However, spatial filtering methods are highly sensitive to the input parameters selected for the filter and heavily rely on the experience of data processors, often resulting in over- or under-correction, affecting exploration results. Third, correction using other data unaffected by static shifts. This involves using simultaneously collected electromagnetic data unaffected by static shifts, such as independently collected transient electromagnetic data [1,11,14,22], for correction, which has shown some effectiveness. Although this method avoids human interference, it requires independent acquisition systems, increasing workload and potentially being affected by survey point locations and acquisition parameters. Tournerie et al. [16] proposed a method to estimate apparent resistivity using the spatial distribution of phases, but this method requires at least 30 survey points and points not affected by static shift effects, making it unsuitable for independent or small-scale applications. To address these issues, this paper proposes a static shift correction method based on time-frequency electromagnetic (TFEM) data.

Time-frequency electromagnetic (TFEM) technology is a novel technique that combines time and frequency information for electromagnetic exploration [24,25]. Developed over more than two decades to meet the needs of oil and gas exploration transitioning from structural trap exploration to lithological and other hidden trap exploration, it has become a mature technology for oil and gas exploration and oil-bearing trap detection, widely applied with good exploration results [9,19,21,25,26]. TFEM technology simultaneously measures horizontal electric and vertical magnetic fields, leveraging the complementary advantages of electric and magnetic fields to enhance exploration effectiveness. However, current static shift correction methods in TFEM still use traditional filtering techniques, which may result in under- or overcorrection. To fully utilize the advantage of TFEM’s simultaneous acquisition of vertical magnetic and horizontal electric fields for efficient static-shift correction, this paper proposes a new static shift correction method. Additionally, traditional static shift correction methods lack quantitative standards to judge the suitability of correction. To solve this problem, we introduce the fractal characteristic analysis.

Fractal characteristic analysis has been gradually applied to geophysical data analysis in recent years [2,5,18,20,23]. Multifractal analysis, which studies the global characteristics of a system from its local properties, mainly discusses the distribution laws of probabilities using statistical physics methods. Halsey et al. [6]proposed the mathematical framework for the multifractal spectrum, discussing the singularities of fractal measures and the mathematical foundation of multifractal analysis. We introduce this method to perform multifractal spectrum analysis on electric field data, analyzing the fractal dimension characteristics of static shift effects caused by shallow subsurface electrical heterogeneities, and further exploring the possibility of identifying geological anomaly features.

2. Materials and Methods

2.1. Principle of Static Shift Effects

The secondary field response of the shallow surface anomaly $(x_0,y_0,z_0)$ in the far field case is used as an example to analyze the mechanism of the static shift effect. In electromagnetic exploration, the observed electromagnetic field can be written as the sum of the primary and secondary fields:

$$E=E_1+E_2$$

(1)

$$H=H_1+H_2$$

(2)

where $E_0$ and $H_0$ represent the primary field of the electric and magnetic fields, respectively, and $E_2$ represents the secondary field. Conductivity anomalies in the subsurface can be represented in a similar form:

$$\sigma ={\sigma }_0+\Delta \sigma$$

(3)

where the uniform earth conductivity is ${\sigma }_0$, and the difference between the anomaly and the background conductivity is $\Delta \sigma$.

Under far-field conditions, the relationship between the electromagnetic fields can be described by a system of Maxwell’s equations:

$$\nabla \times E=-i\omega \mu H$$

(4)

$$\nabla \times H=\sigma E$$

(5)

where i is the imaginary unit, $\omega$ is the angular frequency of an electromagnetic wave, and $\mu$ is the magnetic permeability in a vacuum.

The Helmholtz equation for the secondary electric field can be obtained by substituting Equations (1)–(3) into (4) and (5):

$${\nabla }^2{E}_2+i\omega \mu E_2=i\omega \mu \Delta \sigma E$$

(6)

If $\Delta \sigma <<\sigma$, the Helmholtz equations can be written in integral form with the help of Green’s function, using the first-order Born approximation and replacing the total field with a primary field:

$$E_2\left(\overrightarrow{r}\right)={\int }_v\Delta \sigma \left(r_0\right)G(r,r_0)E_1{e}^{-ik{z}_0}dv$$

(7)

The above equation can be written as:

$${\displaystyle \frac{{E_2}_j\left(\overrightarrow{r}\right)}{E_1}}={\int }_v{\rho }_n\left({\overrightarrow{r}}_0\right)K_{ej}(\overrightarrow{r},{\overrightarrow{r}}_0)dv$$

(8)

The kernel function of the MT electric field orthogonal problem is:

$$K_{ej}(\overrightarrow{r},{\overrightarrow{r}}_0)=-\Delta {\sigma }_0{G}_j(\overrightarrow{r},{\overrightarrow{r}}_0)E_1{e}^{-ik{z}_0}$$

(9)

If the coordinates of the surface receiving point are $(x_0,y_0,0)$, then it can be set:

$$K_{ej}(\overrightarrow{r},{\overrightarrow{r}}_0)={\lambda }_{ej}(x_0-x,y_0-y,z_0)$$

(10)

In the Cartesian coordinate system, equation(1.8) can be written as:

$${\displaystyle \frac{{E_2}_j\left(\overrightarrow{r}\right)}{E_1}}={\int }_0^{\infty }[{\lambda }_{ej}(x,y,z_0)\ast {\rho }_n(x,y,z_0)]d{z}_0$$

(11)

where * denotes a two-dimensional convolution about x and y. To obtain a more intuitive algebraic form, we perform a Fourier transform of the above equation to obtain the form in the wave number domain and set $\xi ,\eta ,\zeta$ to represent the wave number in the $x,y,z$ direction, respectively.

$${\displaystyle \frac{E_j(\xi ,\eta )}{E_1}}={\int }_0^{\infty }[{\Lambda }_{ej}(\xi ,\eta ,z_0)\ast P(\xi ,\eta ,z_0)]d{z}_0$$

(12)

where $E_j(\xi ,\eta )$ denotes the Fourier transform of ${E_2}_j\left(\overrightarrow{r}\right)$, $\Lambda$ denotes the Fourier transform of $\lambda$, and $P(\xi ,\eta ,z_0)$ denotes the Fourier transform of ${\rho }_n(x,y,z_0)$. If $P(\xi ,\eta ,z_0)$ is regarded as the input and ${\displaystyle \frac{E_j(\xi ,\eta )}{E_0}}$ as the output, $\Lambda$ can be regarded as the transfer function of the linear system, and similarly, the linear integral equation of the quadratic magnetic field can be obtained in the wave number domain:

$${\displaystyle \frac{H_j(\xi ,\eta )}{H_0}}={\int }_0^{\infty }[{\Lambda }_{hj}(\xi ,\eta ,z_0)\ast P(\xi ,\eta ,z_0)]d{z}_0$$

(13)

The derivation process is omitted, and the electric-field Ex component transfer function is given directly:

$${\Lambda }_{ex}(\xi ,\eta )={\displaystyle \frac{i\omega \mu {\sigma }_0}{i\zeta }}{r}_1(\xi ,\eta )e^{-i(k+\zeta )z_0}+{\displaystyle \frac{{\xi }^2}{i\zeta }}{e}^{-i+\left(k\zeta \right)z_0}$$

(14)

The magnetic field $H_z$ component kernel function is:

$${\Lambda }_{hz}(\xi ,\eta )=-\sqrt{i\omega \mu {\sigma }_0}{r}_2(\xi ,\eta )e^{-i(k+\zeta )z_0}$$

(15)

where $r_1$, $r_2$ are parameters related to the wave number. Observation of the above equation shows that the electric field transfer function has two terms. The first term, which is frequency dependent, is referred to as the dynamic term and the second term, which is frequency independent, is the static term. For the anomaly with fixed depth $z_0$, the dynamic term for the data in the frequency domain varies with frequency $\omega$, while the static term is the same for different frequency data, and the static term is the reason for the static shift effect. And Equation (15) shows that there is no static term in the vertical component of the magnetic field, only a dynamic term. Thus, the magnetic field does not produce static shift effects.

2.2. Time-Frequency Electromagnetic Static Shift Correction Method

Using the fact that the magnetic field or magnetic induction components are not affected by the static shift effect, the electric field is corrected using apparent resistivity as an intermediary. The correction factor is initially derived by averaging the first three (or five) apparent resistivity values of the early magnetic field and the high-frequency electric field, respectively. This correction factor is then applied to the full frequency-domain data. Assuming a total of n measurement points, the specific method is outlined in Figure 1.

2.3. Multifractal Analysis of Electromagnetic Anomalies

Based on the theoretical analysis of static shift generation, it is evident that near-surface inhomogeneous conductivity bodies have different effects on the electric and magnetic fields. Therefore, we introduce the multifractal analysis method to distinguish the electric and magnetic field anomaly characteristics caused by near-surface inhomogeneous conductivity bodies.

Multifractals, also known as multiscale fractals, represent an advanced development of fractal theory in physics. The multifractal approach studies the overall characteristics of a system by examining its local features. Statistical physics methods are utilized primarily to discuss the distribution patterns of probability. The idea behind multifractal singularity spectrum analysis is to characterize the singularity properties of the normalized measure (distribution) defined on a fractal set, to express its scaling behavior. The singularity of the distribution is represented by the singularity strength $\alpha$ and the corresponding distribution density $f\left(\alpha \right)$. The multifractal spectrum is then determined by the range of values of $\alpha$ and the function $f\left(\alpha \right)$. In order to calculate the multifractal singularity spectrum of electromagnetic field data, the box-counting method [6] is used in this paper:

For any given set of 1D electromagnetic field data, denoted as X, we divide it into several boxes with a scale $\epsilon$, covering the data. The probability measure of the normalized number of boxes $P_i\left(\epsilon \right)$ is then spatially defined.

$$P_i\left(\epsilon \right)={\displaystyle \frac{N_i\left(\epsilon \right)}{N}}$$

(16)

where $N_i$ is the electromagnetic field data (mass, number of points) in the i-th box , and N is the total sum of the 1D electromagnetic field data (mass, number of points). As in $\epsilon \to 0$, it can be assumed that the data are uniformly distributed across the boxes. For multifractal measures, for the i-th box on scale $\epsilon$, $P_i$ follows the following proportional relationship:

$$p_i\left(\epsilon \right)\sim {\epsilon }^{{\alpha }_i}$$

(17)

$$N_{\alpha }\left(\epsilon \right)\sim {\epsilon }^{f\left(\alpha \right)}$$

(18)

where ${\alpha }_i$ is the singularity scaling exponent and a partition function is defined as follows:

$$Z_q\left(\epsilon \right)=\sum _i^{n_i\left(\epsilon \right)}{P}_i^q$$

(19)

where $n_i\left(\epsilon \right)$ is the number of boxes on the scale $\epsilon$, and q is referred to as the weight factor. Its value can be arbitrary, but when $\left|q\right|$ reaches a certain level, further increases will no longer affect the multifractal spectrum.

Mass index:

$$\tau q=\underset{\epsilon \to 0}{lim}{\displaystyle \frac{log{Z}_q\left(\epsilon \right)}{log\epsilon }}$$

(20)

According to the Legendre transform, we can derive the multifractal singularity spectrum relationship for the electromagnetic field data:

$$f\left({\alpha }_q\right)=q{\alpha }_q-{\tau }_q$$

(21)

$${\alpha }_q={\displaystyle \frac{d{\tau }_q}{dq}}$$

(22)

Taking ${\alpha }_q$ as the horizontal coordinate and $f\left({\alpha }_q\right)$ as the vertical coordinate, the multifractal spectrum of the data from the electromagnetic field in 1D X can be obtained. The width of the multifractal spectrum, denoted as $\Delta \alpha$, is referred to as the inhomogeneity index:

$$\Delta \alpha ={\alpha }_{\max }-{\alpha }_{min}$$

(23)

For continuous and differentiable electromagnetic field data, $\Delta \alpha$ is a quantitative measure of the total rate of change (gradient) of all data points in the entire data set. Data that are not differentiable, such as apparent resistivity variations in electromagnetic methods, serve as the overall quantitative indicator of the data’s inhomogeneity.

The width of the distribution density, denoted $f\left(\alpha \right)$, is referred to as the roughness index:

$$\Delta f\left(\alpha \right)={f\left(\alpha \right)}_{\max }-{f\left(\alpha \right)}_{\min }$$

(24)

where $\Delta f\left(\alpha \right)$ represents the difference in the distribution density extrema of all data points within the entire electromagnetic field data. The smoother and more uniform the data, the smaller the roughness index.

The inhomogeneity index $\Delta \alpha$ and the roughness index $\Delta f\left(\alpha \right)$ are two key parameters that describe the multifractal characteristics of the data. When the electromagnetic field data do not exhibit multifractal characteristics, both $\Delta \alpha$ and $\Delta f\left(\alpha \right)$ are zero. The more pronounced the multifractal characteristics of the data, the larger the values of $\Delta \alpha$ and $\Delta f\left(\alpha \right)$.

3. Results

3.1. Example of Static shift Correction in a Theoretical Model

We design the 2.5D subsurface anomaly model Figure 2(a) and calculate its electromagnetic response Figure 2(c), and as a control, we first calculate the electromagnetic response of anomalies without a shallow surface Figure 2(b) and compare it with the results of the processing using this algorithm Figure 2(d).

It can be observed that the current algorithm effectively suppresses static shift, and the results closely match the theoretical results without static shift effects.

3.2. Example of Static Shift Correction on Measured Data

The time-frequency electromagnetic field measurement data Figure 3 from a specific region shows an apparent resistivity pseudo-section with a significant amount of longitudinal strip-like anomalies.

After processing with the algorithm proposed in this paper, the results shown in Figure 4 reveal significant changes in the electrical properties near the measurement points 10 to 17, which may correspond to regions with localized electrical property discontinuities in the subsurface. Near measurement point 65, there appears to be a noticeable low-resistivity layer. In the figure, it is evident that the proposed algorithm effectively eliminates most of the static shift, allowing the profile to clearly reflect the resistivity characteristics of the subsurface.

3.3. Multifractal Spectrum Analysis of Electric Field in a Simple Electrical Model

Using the two parameters in the multifractal spectrum, we can distinguish deep anomalies that are obscured by static shift effects caused by shallow inhomogeneous bodies in the electric field data and also identify magnetic field anomalies induced by shallow anomalous bodies.

A 2.5D electrical profile model is designed as follows: The plate-like anomaly has the same depth as the upper surface of the shallow anomaly, and the resistivity is the same. The locations of the measurement points are shown in Figure 5, where the TM configuration is used to measure the electric field component $E_x$ and the vertical magnetic field Hz. Figure 6 (a) shows the curve of the components of the electric field $E_x$ component curve, and Figure 6(b) presents the multifractal spectrum.

In the data of Figure 6(a), the electric field curves obtained at measurement points over the shallow subsurface anomaly and the vertical plate-shaped anomaly in the simple electrical model are quite similar in shape. In the presence of numerous shallow anomalies, the differentiation of them requires experience. In Figure 6(b), it can be seen that the multifractal spectrum of the shallow subsurface anomaly differs very little from the case without anomalies, indicating that the resistivity anomaly of the shallow subsurface anomaly manifests primarily as a static component effect. In contrast, the multifractal spectrum of the vertical plate-like anomaly shows significant differences, reflecting the presence of a true anomaly.

3.4. Multifractal Feature Analysis of Field Data

As previously discussed, the differences in multifractal spectra between shallow subsurface anomalies and plate-like anomalies substantiate that variations in the underlying anomalous bodies lead to different inhomogeneity and roughness characteristics in field data. We conducted a multifractal spectrum analysis on the field data before and after correction (Figure 3 and Figure 4). After correction, significant changes were observed in apparent resistivity values; however, the two multifractal scaling exponents exhibited minimal variation (Figure 7). This suggests that our static effect correction process effectively preserved the gradient variations between apparent resistivity values at different frequencies, along with information related to roughness and other inherent characteristics.

By performing a multifractal spectrum analysis on the survey line (Figure 8 and Figure 9), it can be observed that the values at both ends of the line are significantly higher than those in the middle section, indicating that the complexity of the data is higher at the ends, suggesting more complex geological conditions. Although there are numerous longitudinal strip-like anomalies in the middle section, the broadly distributed and stable low values indicate that the anomalies in this segment of the profile are primarily caused by static shift effects.

We also performed Gaussian mixture clustering in the multifractal spectrum, resulting in three categories of measurement points (Figure 10). The yellow regions are the most numerous and represent the most common data points in the profile. The cluster center of this region reflects the average data complexity of most of the measurement points. The points in the red region exhibit higher data complexity, indicating that the subsurface electrical structure beneath these points is more complex. In contrast, the points in the green region have very low complexity, suggesting a simpler subsurface electrical structure beneath these locations.

In Figure 11, we marked the high-complexity and low-complexity measurement points, as identified by clustering, on the corrected apparent resistivity profile. Red indicates high-complexity data points, and green indicates low-complexity data points. It is clearly observed that the red dots correspond to areas with more significant electrical variations, while the green dots exhibit relatively small changes in both electrical properties and frequency. This shows that the multifractal spectrum effectively distinguishes the complexity of the subsurface electrical structure in the field data.

4. Discussion

From a mathematical perspective, the multifractal spectrum is highly sensitive to the overall gradient (inhomogeneity) and smoothness (roughness) of the data. The author envisions that this characteristic could be applied to geophysical inversion algorithms in the future. During the iterative solution process, not only would the measured data values be fitted, but also the multifractal spectrum of the forward response would be matched with that of the measured data. This approach could potentially improve both the efficiency and the accuracy of the inversion process. As inversion not only considers numerical restoration, but also requires the restoration of the relative relationships between the data. Multifractal spectrum analysis provides valuable insights for further data analysis.

5. Conclusions

We proposed a new time-frequency electromagnetic static shift correction method, taking advantage of the advantages of simultaneous measurement of electromagnetic fields in the time-frequency electromagnetic method. The apparent resistivity defined by the high-frequency magnetic field is used to correct the electric field. The method was validated using field data, demonstrating that our algorithm effectively eliminates the impact of static shift effects. In addition, we discussed the applicability of the multifractal algorithm in detecting the complexity of electromagnetic field data. We showed that, after static shift correction using this method, the changes in multifractal parameters are minimal, indicating that the method preserves the complete electrical information of the subsurface while correcting for static shift. Multifractal spectrum analysis can intuitively extract subsurface electrical complexity information from electrical field or apparent resistivity data, which not only evaluates the effectiveness of static shift correction from multiple dimensions but also provides important guidance for subsequent inversion processes.