Anomaly Field Extraction Based on Layered-Earth Model and Equivalent Eddy Current Inversion: A Case Study of Borehole zk506 in Baishiquan, Xinjiang

TL;DR

This paper introduces a new method to extract and interpret geophysical data from boreholes to locate hidden ore deposits using a layered-earth model and eddy current inversion.

Contribution

The novel approach combines a layered-earth model with equivalent eddy current inversion to extract pure anomalies from borehole electromagnetic data.

Findings

The method successfully extracted pure anomalies from BHTEM data in the zk506 borehole.

The nickel ore body was identified as dipping northwestward southeast of the zk506 borehole.

Validation with nearby boreholes confirmed the accuracy of the inversion results.

Abstract

Based on the concept of an equivalent eddy current, anomaly field inversion provides an efficient and rapid inversion method for borehole transient electromagnetic (BHTEM) measurements. It enables the utilization of the equivalent eddy current to rapidly process and interpret BHTEM data. This method allows for the accurate determination of the central position and spatial distribution of anomalies. However, the equivalent eddy current method is solely applicable to the inversion of the anomaly field. Given that the measured data frequently contain strong background field information, it is challenging to directly apply the equivalent eddy current approach to the inversion and interpretation of the measured data. To address the aforementioned issues, in this study, we innovatively put forward a method. Specifically, we utilize the response of the layered earth to simulate the background field and subtract the background field from the measured data through the “difference method” to extract the anomaly field. Subsequently, by integrating the equivalent eddy current method, the inversion of the anomaly field was accomplished. Eventually, a rapid quantitative inversion and interpretation of BHTEM data were achieved. We applied this approach to extract the pure anomaly from the measured data of the zk506 borehole in the Baishiquan mining area, Xinjiang, and then conducted equivalent eddy current inversion. The spatial position and distribution characteristics of the concealed ore bodies near the zk506 borehole were successfully pinpointed. Validation by the zk507 and zk508 boreholes confirmed that the main anomaly of the nickel ore body is positioned in the southeast of the boreholes, dipping northwestward. This outcome validates the accuracy of the BHTEM inversion interpretation and rectifies the geological understanding obtained from the zk506 single borehole. It demonstrates the effectiveness and significance of the pure anomaly extraction based on the layered-earth model and equivalent eddy current inversion in the exploration of high-conductivity sulfide ores.

Figures

Full text

1. Introduction

Mineral resources constitute the bedrock for a nation’s economic development. Exploration of these resources is of paramount importance in safeguarding national resource security. At present, with the rapid depletion of shallow mineral deposits, prospecting efforts are gradually shifting to greater depths. When confronted with complex terrains and the need to detect deep-seated ore bodies, traditional surface exploration methods encounter challenges. For instance, issues such as limited accuracy and significant interference arise [1,2,3]. These problems demand higher standards for both the depth and precision of detection [4,5]. Addressing the limitations of traditional surface exploration methods in complex terrains and deep-seated ore body detection, the BHTEM method emerges as a viable solution. Operating on a ground-emission and borehole-reception principle, the BHTEM method [6] offers distinct advantages. Firstly, the earth acts as a natural filter, effectively reducing ground-based interference [7]. Secondly, with the borehole probe in close proximity to anomalies, it captures stronger signals [8]. This approach overcomes the drawbacks of surface-based techniques, providing a precise and in-depth tool for mineral exploration [9]. Over the years, the BHTEM method has proven highly successful in the exploration of sulfide-type copper/nickel deposits [10,11,12,13]. It has become an essential technique for identifying deep-seated and concealed ore bodies, as demonstrated by numerous studies [14,15,16,17]. With the development of artificial intelligence technology, machine learning methods have been increasingly applied to geophysical data processing, significantly enhancing data interpretation capabilities [18,19]. This highlights the importance of new information processing technologies in geophysical exploration [20,21,22]. The primary objective of BHTEM detection is to pinpoint the location of blind ore bodies at the bottom or adjacent to the borehole, thereby providing crucial support for ore body tracing. Currently, there are three main positioning methods. The first method involves determining the anomaly center through vector intersection, leveraging the magnetic field response generated by induced eddy currents [23]. This approach is relatively fast; however, it falls short in discerning the spatial distribution of anomalies. The second method is the equivalent eddy current inversion [24,25,26]. By equating the anomaly response to that of an eddy current, this method can rapidly ascertain the spatial location and distribution of anomalies. Nevertheless, it necessitates obtaining the pure anomaly response of the underground abnormal body. The third method is 3D inversion [27,28,29]. This technique can yield accurate information regarding anomaly location, spatial distribution, and resistivity. However, it involves generating a global mesh, which typically demands substantial computing resources. Additionally, the data from a single borehole is rather sparse in the inversion space, failing to impose effective constraints on the global grid. As a result, the accuracy and uniqueness of the inversion are compromised, making it difficult to achieve satisfactory results in practical applications.

In addressing these challenges, this study presents an innovative technology that integrates anomaly field extraction with equivalent eddy current inversion. Initially, this approach utilizes a layered-earth model to fit the background, effectively isolating the pure response of the underground ore body. Subsequently, the equivalent eddy current method is employed to model this pure response, enabling the precise determination of the location of deep-seated anomalies. This not only represents an advancement in the equivalent eddy current inversion technique but also substantially enhances its practical applicability. Finally, through a comprehensive field case study, this method is applied to the BHTEM-measured data. It successfully predicted the spatial location and distribution characteristics of the ore body, and the results were validated through drilling.

2. Theoretical Methods and Case Studies Verification

2.1. Equivalent Eddy Current Theory

The BHTEM method generally adopts a ground-transmitting and borehole-receiving configuration, as depicted in Figure 1. When a step current is applied to the transmitting loop, it generates a primary magnetic field. This primary magnetic field induces a secondary current within the underground conductive rocks or ore bodies. Once the transmitting current is abruptly terminated, the receiving coil captures the time-varying signal of the secondary magnetic field. By meticulously analyzing these signals and applying the principle of electromagnetic induction, which relates the induced electromotive force to the rate of change of the secondary magnetic field, we can infer the resistivity information of geological bodies at different depths underground.

Figure 2 depicts the physical process of induced eddy current generation and dissipation within a conductive thin plate under a primary magnetic field. When a conductive thin plate is placed in a primary magnetic field and the transmitting process suddenly stops, according to Faraday’s law, to maintain the influence of the primary field before shutdown, induced eddy currents will arise within the conductive thin plate. These induced eddy currents form a series of eddy currents inside the plate. Resembling the shape of the plate, these loops vary in intensity and size. Initially, these loops are mainly distributed at the edges of the plate and are relatively large. As time passes, they gradually spread toward the center of the conductor and eventually dissipate due to heat loss. At any moment during this process, the induced current distribution within the conductive thin plate can be approximated by the field generated by an equivalent eddy current [24]. Therefore, the characteristics of a specific conductive thin plate can be described through an equivalent eddy current (Figure 3). This loop, located in free space, can be defined by seven parameters, and its magnetic field in free space can be calculated using Formulas (1) and (2).

$$[eqn]$$ $$[eqn]$$

Among them, $[eqn]$ , $[eqn]$ is the radial component of the magnetic field, $[eqn]$ is the axial component of the magnetic field, $[eqn]$ and $[eqn]$ are the first and second types of elliptic integrals, $[eqn]$ is the radius of the circular eddy current, $[eqn]$ is the free-space magnetic permeability, and $[eqn]$ is the current intensity of the eddy current.

2.2. Background Field Calculation

The layered-earth model is employed in the simulation research of the background field. The rectangular-loop transmission source of transient electromagnetics can be approximated by multiple electric dipoles of a certain length. The response of the electric dipoles is derived using Kerry Key’s [30] method. The vector potential generated by a horizontal electric dipole oriented in the x-direction at any point on a layered medium is as follows:

$$[eqn]$$ $$[eqn]$$

The corresponding magnetic field components are expressed as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Among them, $[eqn]$ represents the integration kernel in the x-direction, $[eqn]$ is the vector potential in the x-direction, $[eqn]$ is the integration kernel in the z-direction, and $[eqn]$ denotes the vector potential in the z-direction. $[eqn]$ stands for the Bessel function of the first kind of order zero. $[eqn]$ , $[eqn]$ , and $[eqn]$ denote the magnetic field strengths in the x, y, and z directions, respectively.

The magnetic field response obtained through the above steps is a complex integral containing Bessel functions, which is solved using Hankel digital filtering calculation [31]. Its general expression can be written as follows:

$$[eqn]$$

Among them, $[eqn]$ is the sampling interval, $[eqn]$ is the number of filter points, and $[eqn]$ is the filter coefficient. After completing the calculation of the frequency domain field, the cosine transform is used to convert the field values into the time domain. Taking $[eqn]$ as an example, the numerical discrete approximation expression is written as follows:

$$[eqn]$$

Among them, $[eqn]$ is the vertical magnetic field, $[eqn]$ is the angular frequency, $[eqn]$ is the sampling delay, and $[eqn]$ is the $[eqn]$ .

Similarly, the electromagnetic field response of a directional dipole can be obtained.

To validate the accuracy of the algorithm, Maxwell(5.4.53.17914), a widely recognized and industry-proven commercial software for transient electromagnetics, was selected for comparison (for detailed information, refer to https://www.electromag.com.au/maxwell.php (accessed on 5 April 2025)). The designed model, as presented in Table 1, consists of five layered-earth models. The transmitting and receiving parameters adopted are as follows: the transmission frame size is 200 m $[eqn]$ 200 m, the transmission current is 1A, the transmission time base is 6.25 Hz, the sampling time is Crone 50 ms (there are 36 sampling channels), the rising ramp is 0.00001 ms, the falling ramp is 0.00001 ms (there are 100 measuring points), and the interval is 5 m. Unless otherwise specified, the delay chosen for the numerical simulations below is Crone 50 ms.

As can be observed from Figure 4a–e, the response computed by the self-compiled program and those from the Maxwell lerio module under different layered models are generally in good agreement. With the exception that the amplitude of the first channel shows a slight difference, the other channels are fitted well. Table 2 indicates that the maximum average relative error for each measurement point within the same channel is 3.96%, and the maximum relative error of a single measurement point in a single channel is 4.87%.

2.3. Example of Anomaly Field Extraction

In the previous section, the validity of the background field algorithm has been verified. Next, we will integrate anomaly field extraction (background field elimination) with the equivalent eddy current technique to perform an inversion test on the theoretical model, so as to evaluate its applicability.

In this study, we employ the BHTEM response of a conductive thin plate within a uniform half-space as the test data. The resistivity of the uniform half-space is set at 300 Ω·m, the dimensions of the transmitting frame are 200 m × 200 m, and the sampling delay is 50 ms. To be consistent with actual working parameters, the rising ramp is set to 0.25 ms and the falling ramp to 0.5 ms. The conductive thin plate has dimensions of 50 m × 50 m, with its center coordinates at (−50, −50, −100). Its longitudinal conductivity is 500 S, the inclination angle is 25°, and the inclination rotation angle is 30°. Under the same background resistivity, the response obtained through the self-compiled program is larger than that obtained using Maxwell (as shown in Figure 5a). Besides considering that the background half-space and thin plate coupling diminish the overall response, the primary influencing parameter is the falling ramp. During the program algorithm verification phase, the falling ramp was set to 0.00001 ms, indicating an impulse response for the emission waveform. However, when the falling ramp was adjusted to 0.5 ms, significant differences were observed in the early stage between the responses calculated by the self-compiled program and Maxwell (Figure 4, ch1–ch9). This discrepancy arises from the finer segmentation of the falling ramp in our program, where a 0.5 ms falling ramp is divided into 30 calculations, leading to a prolonged computation time. In contrast, Maxwell’s software uses a relatively fast but crude method for calculating falling ramp to increase efficiency. Consequently, the calculation results show discrepancies in the early stage, while the fitting is more precise in the middle stage (as depicted in Figure 6, channels 10–26). It is clearly evident from Figure 5a that the response generated by our self-compiled program effectively depicts the characteristics of the background field of the test model. The abnormal field, obtained via the “difference method”, closely approximates the equivalent eddy current field (Figure 5b–d). The spatial positions of the eddy currents corresponding to different channels after inversion are presented in Figure 5e–f. Significantly, the size of these loops is slightly smaller than that of the conductive thin plate in the Maxwell model. Nevertheless, the central coordinates, inclination, and azimuth, which represent the spatial location, distribution, and shape of the anomaly, are essentially the same. These parameters are regarded as typical, and their average values are tabulated in Table 3.

As can be seen from Table 3, the average values of typical parameters for the eddy current indicate that the anomaly location and spatial distribution align with the thin plate parameters in the designed model. This calculation example demonstrates that the anomaly extraction based on the layered-earth and equivalent eddy current inversion is feasible.

2.4. Robustness Evaluation

When simulating the background field for measured data, the selection of parameters is of critical importance. Therefore, it is necessary to evaluate the perturbation caused by minor changes in formation resistivity on the response. We designed three homogeneous half-space models (Table 4) to assess the impact of parameter variations on abnormal field extraction and inversion results by altering the resistivity values of the homogeneous half-space. In the calculation of the model containing a conductive thin plate, the background model is selected as Model 1 in Table 4, and the thin plate model in Section 2.3 is adopted for the thin plate.

Under disturbance conditions, the anomaly extraction results of the “difference method” show that the overall shape of the abnormal values is consistent with the response of the thin plate in free space. The three types of disturbances caused by changing the background resistivity values have limited effects on anomaly extraction, and their impacts only occur in the early stage (Figure 6, channels 1–9), without contributing to the responses in the middle and late stages, and ultimately have almost no impact on inversion. This also demonstrates that the anomaly extraction based on the “difference method” has good robustness. Of course, since transient electromagnetic (TEM) is sensitive to low-resistivity bodies, this method needs to be used with caution in the presence of low-resistivity backgrounds.

3. Baishiquan BHTEM Survey

3.1. Geological Background

The strata exposed in the Baishiquan copper/nickel mining area consist of an ancient metamorphic rock series. These ancient basement rocks have undergone numerous tectonic movements, which provide essential conditions for the formation of ore deposits in this area. The Baishiquan copper/nickel mine is located 3 km south of the Shaquanzi fault, where neutral, basic, and ultrabasic rocks occupy the fault zones, serving as primary conduits for magmatic activity in the area [32,33]. The ore-bearing rock mass consists of a basic–ultramafic complex, characterized by an elliptical lithological unit at its core as the predominant component (Figure 5a). The eastern contact zone displays an arcuate distribution, dipping from northwest to southwest and forming an inward-steepened basinal structure. Copper/nickel mineralization is predominantly hosted within ultramafic lithologies of the complex [34], which show peripheral concentration with minor intra-complex occurrences. Principal lithological units comprise pyroxenite, olivine pyroxenite, plagioclase peridotite, and lherzolite, with lherzolite demonstrating optimal ore potential. To date, 28 copper/nickel ore bodies have been identified in the Baishiquan deposit, the most significant being designated as Ni1, Ni12, and Ni13 [35]. These ore bodies display lenticular, stratoid, and vein-type morphologies, characterized by distinct branching/recombination structures (Figure 5b). Most ore bodies exhibit surface exposure or subcropping features attributed to shallow hypogene positions. Dimensionally constrained, the ore bodies range from 10 to 200 m in strike length (locally > 300 m), with thicknesses averaging 1–15 m (maximum > 15 m).

To evaluate the resource potential systematically, mineral exploration entities have deployed exploration boreholes in the central sector of the mining district, including borehole ZK506 (Figure 7a). Drilling intersected three mineralized intervals of copper/nickel sulfide mineralization. Lithostratigraphic logging delineates a vertical sequence comprising gabbroic norite, mineralized peridotite, and plagiogranite, with the gabbroic norite recurrently intercalated within the sequence. Significant stratiform mineralization occurs specifically at the petrologic contact between gabbroic norite and plagiogranite. Preliminary structural analysis indicates a southeastward dip of the ore zone (Figure 8a, uncorrected geological profile). In order to delineate the spatial configuration and axial geometry of the geophysical anomaly associated with borehole ZK506, BHTEM surveys were conducted within the drill hole.

3.2. BHTEM Data Acquisition

Data acquisition was conducted utilizing a 2.4 kW Digital Pulsed Electromagnetic (PEM) system configured with the following parameters: emission current maintained at 15 A, a time base duration of 20 ms, and a turn-off time of 500 μs. The transmitter system utilizes a central loop configuration of 100 m × 100 m, with the borehole collocated at the geometric center of the loop.

It is evident that there are three distinct abnormal manifestations in the z-component (as shown in Figure 8). These are located at the positions marked by the blue horizontal line at 340 m, the red horizontal line at 365 m, and the green line at 390 m. By analyzing the morphology and amplitude of these anomalies in the early, middle, and late periods, it becomes clear that each of these three anomalies has its own unique characteristics.

The anomaly marked by the blue horizontal line is a low-amplitude negative anomaly, likely indicating a small-scale or low-grade ore body near the borehole. The red-marked anomaly initially shows a small positive amplitude before transitioning to extremely strong negative values in mid-to-late stages, interpreted as a borehole anomaly intersecting the ore head. Its intense mid-to-late-stage response, especially the extremely strong negative anomaly in the late stage, suggests a large-scale or high-grade ore body near the borehole. The green horizontal line marks a positive anomaly with strong early-to-medium-stage signals (weaker than the red anomaly) but no late-stage response, identified as an in-hole anomaly with the borehole position essentially centered within it, indicating a medium-scale or low-grade ore body.

At the depth of 370 m, the negative anomaly response of the z-component indicates the presence of a large-scale or high-grade blind ore body near the borehole. Combined with the characteristics of the x- and y-component curves, this anomaly is determined to be located in the third quadrant. Additionally, the positive response values at one end of the zero point of the horizontal component along the reverse “S”-shaped curve are all greater than the negative response values, suggesting a tendency for this anomaly to tilt toward the west and north directions.

3.3. Pure Anomaly Extraction and Equivalent Eddy Current Inversion

Based on the geological logging of borehole zk506, the lithology and thickness of the formation can be ascertained. Regarding the formation resistivity value, it draws reference from the copper/nickel mining area in the east of Huangshan, which has similar metallogenic conditions (Table 5).

Based on the statistical values of the background resistivity of rocks and ores, the initial formation background thickness and resistivity values were provided (Figure 9b). Through man–machine interactive trial calculations, the formation resistivity was appropriately adjusted until the fitting of the background field met the requirements (Figure 10a). Then, the background field was subtracted from the measured data to obtain the pure anomaly field (Figure 10b). Based on the morphological characteristics of BHTEM-measured curves and borehole logging, three main anomalies can be identified in this borehole. Among them, the primary anomaly occurs at 370 m depth, located in the southeast direction of the borehole. The anomaly at 340 m exhibits the smallest amplitude and is positioned northeast of the borehole. The 390 m anomaly represents an in-borehole anomaly with a depth range between 340 m and 370 m. Generally, the eddy current radius should be smaller than the transmitting frame radius, and the current intensity should be set to relatively low values. The inclination and rotation angles are determined according to curve morphology [26]. The initial parameters of the given eddy current are detailed in Table 6.

4. Discussion

From the fitting curves of each component (Figure 10a), it is evident that the early-stage curves following anomaly extraction exhibit varying degrees of distortion. In contrast, the middle- and late-stage curves are minimally influenced by background factors and can effectively reflect the true anomaly. Consequently, the central coordinates and two angle parameters associated with the middle- and late-stage anomalies carry greater practical significance. These parameters have been systematically tabulated in Table 7 for clarity and reference. Figure 11a–d show the spatial distribution characteristics of the eddy current at different angles.

It can be observed from Figure 9b that, after the extraction of the anomaly field, an anomaly with a narrow opening and large amplitude appears at 370 m on the z-component curve. This suggests that the anomaly exhibits significant extension along the inclined direction. The circular eddy current model fails to adequately simulate this anomaly, resulting in the inability to fit the z-component at its maximum value, although the other components are well fitted. When combined with the shape analysis of the x- and y-component curves, it was determined that the main anomaly (marked by the red horizontal line in Figure 8) was positioned in the southeast direction of the borehole. However, the interpretation of these anomalies based on the survey results of BHTEM does not conform to the geological logging of the zk506 borehole at that time. According to the drilling outcomes of the zk506 borehole, geologists at that time believed that the anomalies had a gentle inclination toward the southeast and were nearly horizontal (Figure 9a). Through the logging of the zk507 and zk508 boreholes constructed in the subsequent stage, geologists recognized that the main ore body was located in the southeast direction of the zk506 borehole and inclined toward the northwest (Figure 7b). This finding was consistent with the ore body tendency and orientation inferred by the BHTEM, thus demonstrating the significance and necessity of developing the BHTEM method in the exploration of sulfide copper/nickel ore deposits.

The late response of the z-component of the measured data is strong and sharp. Considering that it is caused by the uneven conductivity and long extension of the local anomaly, the thin plate can basically reflect the spatial distribution of the anomaly. Figure 8a presents the original logging details of borehole zk506. At that time, there were three ore layers. The first layer was lean ore, while the second and third layers were high-grade ore. Prior to conducting BHTEM measurements, geological personnel considered Layer 3 to have considerable thickness and predicted it to be the primary ore layer, determining that the ore body exhibited a gently inclined southeastward attitude. After the BHTEM survey, based on the characteristics of the z-component response curve, it is inferred that the extremely strong Anomaly-2 is the main anomaly (Figure 8) which suggests that there is a large-scale or high-grade blind ore in the vicinity of the borehole. Comprehensive inference based on combined x- and y-component responses indicates that the second-layer ore body intersecting the ore head is positioned in the borehole’s southeast direction and displays a gentle northwestward dip. Subsequent drilling operations on zk507 and zk508 confirmed this conclusion (Figure 9c), and ultimately, geological experts revised the geological profile of zk506. It is worth noting that, after the rock and mineral testing, the grade of the third abnormal body was downgraded in the subsequent compilation. This is consistent with the interpretation of the measured curve, meaning that the scale is large but the grade is relatively low. As a result, the z-component response curve only shows reactions in the early and middle stages and exhibits no abnormality in the late stages.

Based on the distribution of the inverted eddy currents (Figure 10), the anomaly at 370 m is identified as the main anomaly. It is situated in the southeast of the borehole, which is in accordance with the description in the geological logging of borehole zk506. Notably, following the actual inversion of the high-grade Ni17 ore body anomaly, the center of the anomaly was slightly uplifted, falling within the depth range of 360 m to 365 m. This indicates that the anomaly in the southeast direction exhibited an upward-trending pattern. Notably, due to the spatial separation between the firing position of the launch frame and the Ni17-1 ore body (Figure 7b), full excitation coupling with the ore body could not be achieved. Additionally, borehole zk506 was located at a significant distance from the Ni17-1 ore body, resulting in the Ni17-1-induced anomalies in the measured data being both small in magnitude and attenuated, thus rendering them difficult to identify.

5. Conclusions

This paper proposes an anomaly field extraction and equivalent eddy current inversion method based on layered-earth responses, enabling rapid quantitative inversion of BHTEM data. The correctness and applicability of this method are verified through theoretical calculation examples and borehole field measurements, providing a simple, rapid, and feasible inversion-interpretation approach for BHTEM data processing and interpretation and further improving the technical level of inversion interpretation. However, this method also has shortcomings and limitations. First, it has strong applicability to simple geological models, but there may be deviations in the inversion and interpretation of complex geological models. Therefore, in future research, advanced technologies such as machine learning can be used to implement 3D inversion to more accurately invert the electrical structure of subsurface media. Second, if high-conductivity layers or saturated aquifers exist in the shallow surface, they may act as interference bodies and mask ore body anomalies. Therefore, multi-parameter geophysical surveys may be an effective way to improve the credibility of interpretations.