3D LOTEM inversion based on initial-model optimization—a case study of deep karst detection

Abstract

The long-offset transient electromagnetic method (LOTEM) offers a large depth of investigation and high sensitivity to subsurface resistivity variations, making it valuable for deep resource exploration, oil and gas prospecting, and engineering investigations. As exploration targets move to greater depths, increasingly undulating surface and subsurface interfaces and more complex structural settings make LOTEM responses more difficult to relate accurately to the true subsurface resistivity distribution. Under such conditions, 3D LOTEM inversion faces greater challenges in background-structure representation, target-anomaly recovery, and result reliability. To address the inadequate representation of background resistivity structures by conventional initial models under complex backgrounds, and the resulting impact on inversion reliability, this study proposes an improved 3D LOTEM inversion strategy based on initial-model optimization. Prior information capable of reflecting the characteristics of background resistivity distribution and structural relationships is used to guide the construction of the initial model, thereby providing a more suitable model basis for the early stage of inversion. In synthetic examples, a background initial model is first constructed from station-wise 1D inversion results, and a horizon-constrained initial model is then further established by incorporating horizon information. The results show that the proposed strategy can significantly improve the recovery of both background structures and target anomalies. Finally, the strategy is applied to LOTEM field data for deep karst investigation in a shale gas exploration area in southwestern China, where a deep low-resistivity anomaly is identified near the target interval. Combined with the regional geological setting, karst development patterns, and well-seismic constraints, the anomaly is inferred to represent a deep karst-developed zone and provides geophysical support for drilling planning. Overall, the proposed strategy improves the representation of background resistivity structures, reduces the influence of background structures on target recovery, and enhances the recovery performance and reliability of inversion results, thereby providing technical support for methodological improvement of 3D LOTEM inversion and for deep-target detection in complex resistivity settings.

1 Introduction

The long-offset transient electromagnetic method (LOTEM) uses a grounded-wire source to inject pulsed current into the subsurface and, after current shut-off, records the horizontal electric field and vertical magnetic field generated by induced eddy currents, thereby providing a more accurate characterization of subsurface resistivity structure (Figure 1). As a secondary-field method, LOTEM offers a large depth of investigation and high sensitivity to low-resistivity anomalies, and has therefore been widely applied in deep resource exploration, oil and gas prospecting, and engineering investigations (Keller et al., 1984; Strack et al., 1990; Strack, 1992; Khan et al., 2018; Commer et al., 2006; Cheng et al., 2023; Haroon et al., 2015; Xie et al., 2016; Yan et al., 2018; Wang et al., 2023; Rezaei Mirghahed et al., 2024; Carter et al., 2021).

FIGURE 1

As exploration targets extend into deeper and more complex subsurface electrical settings, responses from deep targets are commonly superimposed on background responses. Accordingly, the focus of LOTEM research has gradually shifted from response-pattern analysis to three-dimensional inversion and quantitative interpretation of complex subsurface resistivity structures. In this context, the limitations of conventional inversion and interpretation methods based on one-dimensional assumptions have become increasingly evident (Stephan et al., 1991; Hördt et al., 1992; Müller et al., 2002). On the one hand, 1D inversion generally assumes horizontally layered media and neglects the influence of lateral resistivity variations on electromagnetic responses. In areas with complex resistivity structures and pronounced lateral heterogeneity, this assumption is often difficult to satisfy (Hördt et al., 1992; Scholl, 2005; Oldenburg et al., 2013). When strong lateral resistivity contrasts or localized anomalous bodies are present, fitting the observed data may introduce non-physical vertical layering unrelated to the true geology, leading to spurious anomalies in the inversion results. Such bias can result in erroneous estimates of anomaly position, geometry, and resistivity, thereby substantially reducing interpretation reliability (Hördt et al., 1992; Scholl, 2005; Oldenburg et al., 2013). On the other hand, in practical exploration settings, sedimentary cover, mudstone caprocks with large thickness variations, and multiple stacked strata are commonly developed, and their electrical structures typically exhibit strong heterogeneity and multiscale characteristics. These overburden units not only produce pronounced resistivity layering in space, but also modify electromagnetic responses through their resistivity contrasts with deep host rocks and anomalous bodies. As a result, shielding and distortion effects can suppress the effective responses of deep targets (Haroon et al., 2015; Tang et al., 2011; Hou et al., 2017; Xue et al., 2021; Xue et al., 2022). Under such conditions, 1D inversion cannot effectively separate background responses associated with the overburden from localized responses related to deep anomalies, making it difficult to achieve both reliable recovery of stratigraphic structure and detailed delineation of deep targets (Xie et al., 2016; Yan et al., 2018; Oldenburg et al., 2013). Therefore, under complex background conditions, 1D inversion results alone are often insufficient for refined characterization of deep targets.

To improve interpretation under complex geological conditions, extensive research has been carried out on three-dimensional forward and inverse modeling of transient electromagnetic data. Three-dimensional algorithms based on finite-difference, finite-volume, and finite-element methods have developed rapidly (Wang and Hohmann, 1993; Commer and Newman, 2004; Haber et al., 2007; Newman and Commer, 2005), providing theoretical and technical support for detailed characterization of complex resistivity structures (Oldenburg et al., 2013; Ren et al., 2018; Qi et al., 2020; Sun et al., 2021). Among them, vector finite-element methods based on unstructured tetrahedral meshes are particularly well suited to rugged topography and irregular geological bodies because of their flexible local refinement capability, and have shown good applicability in practical data processing (Li et al., 2020; Liu et al., 2019; Wang et al., 2023). Even so, under complex exploration conditions, substantial challenges remain in converting efficient 3D forward modeling and high-accuracy numerical discretization into stable and reliable inversion results. On the one hand, because 3D inversion is strongly nonlinear and non-convex, its performance depends heavily on the initial model (Oldenburg et al., 2013; Haber et al., 2007; Yang et al., 2014). In complex resistivity settings, improving forward accuracy usually requires large-scale mesh discretization, especially local refinement near the source and receivers, which greatly increases the number of model parameters and makes the solution space of the inversion objective function more complicated. When prior information is lacking, the inversion can easily become trapped in local minima, resulting in slow convergence or even failure to converge, which limits the practical effectiveness of 3D inversion (Newman and Commer, 2005; Yang et al., 2014). On the other hand, in real applications, widely developed sedimentary cover and multiple stacked strata impose more complicated background constraints on the inversion. Strong resistivity contrasts between the overburden, deep host rocks, and anomalous bodies tend to dominate the electromagnetic response, so the inversion preferentially fits the background structure and becomes less sensitive to deep anomalies (Haroon et al., 2015; Xie et al., 2016; Yan et al., 2018). Under these conditions, when a homogeneous half-space is used as the initial model, the early stage of inversion is often dominated by background-structure adjustment. This increases the proportion of large-scale non-target updates, and non-physical resistivity adjustments may then be introduced to compensate for data residuals, ultimately degrading the recovery of anomaly position, boundaries, and stratigraphic relationships (Hördt et al., 1992; Scholl, 2005; Oldenburg et al., 2013). Previous studies have shown that a reasonable initial model can significantly improve the convergence path of the inversion objective function (Ye et al., 2013; Newman and Commer, 2005; Yang et al., 2014). Therefore, for 3D inversion in complex exploration settings, constructing an initial model that better reflects realistic geological conditions has become a key issue. Although complete, high-resolution 3D prior information is rarely available in practice, the layered-response characteristics contained in transient electromagnetic data provide an important source of information for initial-model optimization. Making full use of the resistivity trends indicated by station-wise 1D inversion results to guide the construction of a reasonable 3D initial model is therefore a practical way to improve inversion performance under complex geological conditions.

Based on these considerations, this study investigates initial-model optimization for 3D LOTEM inversion under complex background conditions. To address the inadequate representation of background structure by a homogeneous half-space initial model and the resulting impact on inversion reliability, we propose an improved 3D LOTEM inversion strategy based on initial-model optimization. Station-wise 1D inversion results are used to assign reasonable background resistivity contrasts to the initial-model space, and horizon information is further incorporated to constrain the spatial distribution of resistivity. Unlike methods that impose prior constraints directly during inversion, the present study uses prior information primarily for initial-model construction, so that it does not continue to participate in model updating during iteration. This places lower requirements on the stability and constraint strength of the prior information, and is better suited to practical datasets in which prior information varies in source, completeness, and reliability. Following this idea, a progressive set of synthetic experiments is used to evaluate the recovery of background structures and target anomalies under different initial-model conditions. Quantitative metrics are then used to assess the degree of improvement in inversion results and how it changes with background complexity, while also discussing the way prior information is used and the conditions under which the strategy is applicable. Finally, the proposed strategy is applied to field LOTEM data for deep karst investigation to verify its applicability under complex background conditions.

2 Three-dimensional LOTEM forward modeling

2.1 Time-domain vector finite-element formulation on unstructured tetrahedral meshes

Starting from Maxwell’s equations combined with the constitutive relations, the electric field satisfies the diffusion-type governing equation given in Equation 1 (Hu, 2021; Um, 2011; Li et al., 2020; Cai et al., 2022):where denotes the impressed source current density at time t and source location , is the electrical conductivity, and is the magnetic permeability of free space (taken as for non-magnetic media).

To accurately represent complex topographic variations and subsurface heterogeneous electrical structures, the computational domain is discretized using vector finite-element formulation on unstructured tetrahedral meshes (Figure 2a). Specifically, first-order edge-based vector basis functions are adopted to represent while enforcing tangential continuity across element faces (Wang et al., 2023; Um, 2011).

FIGURE 2

The numbering of nodes and edges of a tetrahedral element is shown in Figure 2b. For the tetrahedral element ,the electric field is expanded using edge-based vector interpolation functions as follows:where is the degree of freedom associated with the edge of the element, and is the corresponding edge-based vector basis function, which is expressed in Equation 3:where and denote the vertex indices associated with the edge, represents the edge length, and denotes the scalar nodal shape function used to construct the edge-based vector basis. By applying a standard Galerkin weak formulation to Equation 2 and assembling the element contributions, the following semi-discrete matrix system is obtained:where is the stiffness matrix, is the conductivity-weighted mass matrix, and represents the source term. For the element, the corresponding element matrices and vectors are defined in Equations 5–7:

For the time derivatives of the electric field and the source current in Equation 4, an unconditionally stable second-order backward Euler scheme is adopted for temporal discretization, yielding:where denotes the time steps, is the current time step, and is the time-step ratio. Substituting Equations 8, 9 into Equation 4 and collecting terms at time level , we obtain Equation 10, with the corresponding source term given in Equation 11:where is the source current amplitude, is the unit vector indicating the current direction, and denotes the length of each discretized electric-dipole segment along the source wire.

For computational efficiency, Equation 10 can be rearranged into the three-term recurrence form given in Equations 12–14:

After global assembly over all elements and time discretization, the resulting block-banded system can be written as follows (Wang et al., 2023; Liu et al., 2019):

To advance Equation 15 in time, an initial condition for the electric field, is required. For an LOTEM survey with a step-off current waveform, is conveniently initialized using the corresponding quasi-static (DC) solution, which is obtained from the electric potential (Wang et al., 2023; Liu et al., 2019). Accordingly, the initial electric field can be expressed as Equation 16:where denotes the quasi-static electric field, is the spatial sampling/interpolation operator used to evaluate the field at the required locations, and is the electric potential. The electric potential satisfies (Equation 17):

To solve the associated boundary-value problem on a truncated domain, we assume the artificial boundary is placed sufficiently far from the transmitter and impose a homogeneous Dirichlet boundary condition on the outer boundary , yielding Equations 18, 19 (Wang et al., 2023; Liu et al., 2019):where is the source current amplitude and denotes the transmitter (source) location. Using as the test (weighting) function and applying integration by parts, the weak form of Equation 18 can be written as:

The electric potential is discretized using a standard nodal finite-element formulation on the tetrahedral mesh. Within an element, is approximated by the nodal basis expansion:

where denotes the nodal potential at the node of element, and is the corresponding scalar nodal shape function. Substituting Equation 21 into Equation 20, and assembling the element contributions yields the element-level system in Equation 22:where and are given by Equations 23, 24:

After global assembly of all tetrahedral elements with the prescribed degree-of-freedom numbering, the discretized system can be written as Equation 25:where is the global stiffness matrix for the potential equation, is the nodal potential vector, and is the corresponding right-hand-side load vector.

2.2 Verification of forward-modeling accuracy

To verify the accuracy of the three-dimensional time-domain finite-element forward solver used in this study, we consider a representative layered geoelectrical model (Figure 3a) and compare the numerical responses with the corresponding analytical solutions. The model consists of three layers: a resistive surface layer, a conductive intermediate layer, and a resistive basement. The thicknesses and resistivities of these layers are 500 m and 100 Ω m, 100 m and 20 Ω m, and 100 Ω m for the basement, respectively.

FIGURE 3

A grounded-wire source with a length of 100 m is adopted. The transient electric-field responses are evaluated at an observation point located 400 m from the source. The transmitter current is set to 1 A, and the Ex component at the observation point is computed, which is a standard benchmark for layered-earth LOTEM responses.

Figure 3b compares the numerical results with the analytical solutions. The two curves show close agreement over the considered time range, and the relative error remains within 2%, confirming the numerical accuracy of the forward-modeling implementation.

3 Three-dimensional LOTEM inversion methodology and evaluation

3.1 3D LOTEM inversion using the L-BFGS algorithm

In geophysical inversion problems, the objective function is typically composed of two components:where measures the weighted data misfit and penalizes departures from a reference model ; is the regularization parameter; accounts for data uncertainties; and is the roughness (smoothing) operator.

The objective function is minimized using the L-BFGS optimization method (Wang et al., 2023; Liu et al., 2019). Differentiating Equation 26 with respect to yields the gradient:where is the Jacobian (sensitivity) matrix; . In practice, the action of on a vector is computed efficiently by an adjoint-state procedure, significantly reducing memory usage and computational cost.

To derive , we start from the discrete forward system written in operator form as Equation 28:where collects the discrete electric-field unknowns, is the system matrix, and represents the source term.

Taking the derivative of Equation 28 with respect to the model parameters gives:where compactly denotes the parameter-dependent forcing term induced by model perturbations.

The predicted data are obtained by sampling/interpolating the field solution, so the sensitivity of the data can be expressed as:

The interpolation operator can be expressed as , where Q interpolates the electric-field solution to receiver locations and .

Substituting Equations 29, 30 into Equation 27 gives:

As a result, the adjoint forward equation is obtained as Equation 32:

After obtaining the vector u through backward time stepping, it can be substituted into Equation 31 to compute , and the gradient of the objective function can then be evaluated.

For unstructured tetrahedral meshes, the roughness matrix is constructed as Equation 33:where is the number of nearest neighboring cells; is the volume of the inversion cell; represents the neighboring cell associated with the inversion cell; is the maximum volume among all inversion cells; denotes the total volume of the neighboring cells associated with the inversion cell; and is the Euclidean distance between the neighboring cell and the inversion cell.

To ensure physically plausible model parameters, we first determine a feasible conductivity range from available geological and geophysical information, and then impose lower/upper bounds on each inversion cell in logarithmic space. For the inversion cell, we introduce a bounded parameter mapping and write the constraint in the form of Equation 34:

where is the model parameter of the cell, and and are the prescribed lower and upper bounds, respectively. Using the chain rule, the gradient with respect to the transformed variable can be rewritten as Equation 35:

By introducing a new model parameter vector , the constrained optimization problem is transformed into an unconstrained optimization problem. Using the L-BFGS optimization algorithm, the updated model parameters can be expressed as Equation 36:

3.2 Verification of inversion accuracy

To examine inversion behavior in a simple embedded-anomaly setting, we design a single anomalous-body model (Model 1; Figure 4a). The resistivity of air is set to 10⁸ Ω m, the background half-space resistivity is 100 Ω m, and the anomalous body resistivity is 10 Ω m. The anomalous body has dimensions of 1,000 m × 1,000 m × 400 m, with its top buried at a depth of 600 m.

FIGURE 4

A grounded-wire source with a length of 2000 m is used, and the distance between the source and the anomaly center is 2,300 m. Nine survey lines are deployed above the target with a line spacing of 200 m. Each line contains 13 receivers with a station spacing of 200 m, yielding 117 observation points in total. The observation time window spans 0.01 m to 0.1 s with 101 time channels. The Ex responses computed from forward modeling are used as input data for the 3D inversion. To account for measurement uncertainty, 5% Gaussian noise is added to all synthetic data prior to inversion.

Figure 4b shows the RMS misfit curve during the inversion of Model 1. The RMS decreases continuously from an initial value of 4.83 and falls below 1 after 18 iterations. The inversion results are shown in Figures 4c–f. A distinct low-resistivity anomaly is recovered at the target location, and the anomaly is well represented in terms of its central position, lateral extent, and burial depth, with a relatively clear geometry. Taken together, the RMS curve and the recovered model indicate that the proposed 3D LOTEM inversion method achieves stable data fitting and accurately recovers the spatial position and resistivity characteristics of the target anomaly, demonstrating good inversion performance.

3.3 Analysis of synthetic examples

3.3.1 Anomalous body within a single layer

As shown in Figure 5a, two overburden layers are added to the half-space model to form Model 2. The first layer is a high-resistivity overburden with a resistivity of 2000 Ω m and a thickness of 200 m. The second layer has a resistivity of 400 Ω m and a thickness of 200 m. The third layer is the host rock, with a resistivity of 100 Ω m, within which the anomalous body is fully embedded with a resistivity of 10 Ω m. All other parameters are identical to those of Model 1. The observation time ranges from 0.01 m to 0.1 s, and electric-field responses are calculated at 101 time channels.

FIGURE 5

Three-dimensional inversion is carried out using forward responses contaminated with 5% Gaussian noise. A homogeneous half-space is adopted as the initial model, with an initial resistivity of 400 Ω m, corresponding to the average background resistivity. As shown by the RMS convergence curve in Figure 5b, the RMS decreases continuously from an initial value of 37.5 and reaches 0.97 after 48 iterations. The inversion results are presented in Figures 5c–f. Overall, the recovered model reproduces the layered electrical structure and the position of the low-resistivity anomalous body. More specifically, the recovered low-resistivity anomaly broadly corresponds to the true model in spatial position, but its geometry is noticeably distorted, with expanded lateral extent and diffuse boundaries. At the same time, the layering of the overburden is not clearly recovered. The intermediate overburden layer is only weakly expressed in the inversion result, fails to form a distinct boundary, and its resistivity information is not effectively recovered.

3.3.2 Anomalous body spanning multiple layers

In practical exploration settings, deep targets such as concealed ore bodies, water-bearing karst features, and hydrocarbon reservoirs commonly occur in an interlayered form. To further examine the influence of cross-layer structural relationships on inversion recovery, Model 3 is constructed as shown in Figure 6a. Based on Model 2, the thickness of the two overburden layers is increased to 400 m, while the position of the anomalous body remains unchanged, so that the anomaly is located between the two overburden layers. All other model parameters are kept the same.

FIGURE 6

The inversion is carried out using a homogeneous half-space initial model with a resistivity of 800 Ω m. As shown by the RMS convergence curve in Figure 6b, the RMS decreases continuously from an initial value of 33.9 and reaches 0.94 after 64 iterations. Compared with Model 2, Model 3 requires more iterations, and the RMS reduction shows a more pronounced step-like pattern. The inversion results are presented in Figures 6c–f. A distinct low-resistivity anomaly is recovered within the target depth interval, but its geometry deviates noticeably from the true model. Relative to the true model, the anomaly boundaries are more diffuse and the cross-layer relationship is difficult to distinguish clearly. At the same time, the layered background becomes further weakened, and the relative positional relationship between the background structure and the target anomaly is more difficult to resolve than in Model 2.

3.4 Initial-model optimization and analysis of synthetic examples

3.4.1 Initial-model optimization strategy

Synthetic inversion results show that, under different background conditions, a homogeneous half-space initial model can still achieve reasonably good data fitting, but the improvement in model-space recovery remains limited. Noticeable deviations are still present in anomaly boundaries, background structures, and the relative positional relationship between the background and the target anomaly. On this basis, prior-guided initial-model construction is further investigated to evaluate how initial-model optimization improves background-structure representation and inversion recovery.

To improve the representation of background structure in the initial model, prior information capable of reflecting background resistivity characteristics is first introduced. In this study, station-wise 1D inversion results are used as background prior information and mapped onto the 3D mesh to construct a background initial model based on 1D results. Because 1D inversion results may contain local oscillations caused by near-surface interference or inversion instability, they are smoothed before being introduced into the initial model. For the station, the smoothed 1D result at depth is expressed by Equation 37, with the Gaussian weight defined in Equation 38:where denotes the set of neighboring stations associated with the station, and is the Gaussian weight, defined as:where is the horizontal distance between stations and , and is the smoothing width.

The smoothed 1D inversion results are then mapped onto the 3D inversion mesh. For the cell in the mesh, the model parameter obtained by inverse-distance weighting is given in Equation 39:where is the number of neighboring stations used in the interpolation. In this study, the four nearest stations are used, i.e., . The interpolation weight is defined in inverse-distance-squared normalized form as Equation 40:where is the horizontal distance between the center of the cell and the station. The resulting initial model constructed from the 1D results is therefore written as Equation 41:

Under realistic geological conditions, background structure is reflected not only in the overall resistivity trend with depth, but also in the spatial geometry of stratigraphic interfaces, the lateral continuity within layers, and the relatively stable resistivity distribution within certain intervals. On the basis of , horizon information is further introduced to construct a horizon-constrained initial model. Using well data, seismic interpretation results, and other geological information, the 3D inversion domain can be partitioned into horizon units as Equation 42:where denotes the entire inversion domain and is the spatial subdomain corresponding to the horizon unit.

Based on well-log resistivity, layer statistics, and background resistivity information, a reference trend function varying with the depth trend is specified for each horizon unit as Equation 43:where and are the coefficients for the horizon unit and can be determined jointly from well data, layer statistics, and station-wise 1D inversion results. For the cell in the horizon unit, the background model is then combined with the horizon reference trend as Equation 44:where is a weighting coefficient. To further enhance continuity within the same horizon unit and suppress local anomalous variations introduced by interpolation, in-layer smoothing is applied using Equation 45:where is the smoothing parameter, and the normalized weight is defined in Equation 46:

The final horizon-constrained initial model is written as Equation 47:

Compared with the initial model constructed solely from the 1D inversion results, the horizon-constrained initial model further incorporates interface geometry into the model-building stage. This strengthens the representation of background structure and reduces mismatch-driven background adjustments during the early stage of inversion. The resulting inversion workflow is shown in Figure 7.

FIGURE 7

3.4.2 Analysis of synthetic examples

3.4.2.1 Anomalous body within a single layer

As discussed in the preceding section, in Model 2 the anomalous body is located within a single layer, and the main limitations of the recovered result are reflected in two aspects: diffuse anomaly boundaries and insufficient recovery of background structure. Specifically, the anomaly boundaries become spread out, and stratigraphic features are incorporated into the low-resistivity recovery around the anomaly. On this basis, an initial model constructed from 1D inversion results is first adopted to examine how initial-model optimization improves the recovery of background structure and target anomaly.

Figure 8a shows the background initial model constructed from the 1D inversion results. Here, the 1D inversion results of Model 2 without a local anomaly are mapped onto the 3D inversion mesh to form the background initial model. It can be seen that the 1D inversion results clearly reproduce the overall layered resistivity trend.

FIGURE 8

Figure 8b compares the RMS convergence curves for different initial-model settings. After adopting the background initial model constructed from the 1D results, the initial RMS is greatly reduced and decreases rapidly in the early stage, reaching an acceptable data misfit level after 44 iterations. Compared with the inversion result of Model 2, the iteration count is significantly reduced. As shown by the corresponding 3D inversion result, after using the initial model constructed from the 1D results, the stratigraphic relationship between the shallow high-resistivity layer and the underlying background layer becomes clearer, and the resistivity layering of the background structure near the target zone maintains better continuity. The low-resistivity anomaly remains concentrated near the target interval, with its lateral extent clearly reduced and boundary diffusion suppressed. At the same time, the relative positional relationship between the background structure and the target anomaly becomes easier to distinguish.

3.4.2.2 Anomalous body spanning multiple layers

As indicated by the inversion results obtained with the conventional 3D inversion scheme, in Model 3 the anomalous body is embedded between layers, and the main limitations of the recovered result are reflected in two aspects: the coupling between background responses near the interfaces and anomaly responses, and the blurred recovery of cross-layer structural relationships. The layered background becomes further weakened, and the relative positional relationship between the background structure and the target anomaly is more difficult to distinguish. On this basis, the effect of the 1D-result-based initial model on inversion recovery is first examined, and the improvement brought by further incorporating horizon information is then analyzed, with emphasis on the recovery of interface-adjacent zones and cross-layer relationships.

The analysis first considers the initial model constructed from the 1D inversion results. As shown by the RMS convergence curve in Figure 9a, after adopting the 1D-result-based initial model, the initial RMS is substantially lower than that obtained with the homogeneous half-space initial model in the conventional 3D inversion, the RMS decreases more rapidly in the early stage, and then enters a low-value plateau, reaching an acceptable data misfit level after 47 iterations.

FIGURE 9

Figures 9b–e show the corresponding inversion results obtained using the initial model constructed from the 1D results. Compared with the conventional 3D inversion using a homogeneous half-space initial model, the shallow high-resistivity layer is recovered more completely, and the interface undulations and thickness variations are better expressed over the whole section. The relative positional relationship between the anomalous body and the background structure becomes easier to distinguish, and the background low-resistivity values introduced by stratigraphic features are suppressed. A closer comparison of the anomaly itself and the recovery near the interfaces shows that the diffuse low-resistivity recovery around the upper anomaly is reduced and the main body of the anomaly becomes more concentrated, whereas the lower anomaly begins to form a layered recovery pattern that is closer to the true model.

At the same time, this stage of improvement still has clear limitations. Diffuse recovery remains within a certain range around the anomaly. Although the lower anomaly is less enhanced than before, its lateral extent, boundary sharpness, and correspondence with the upper anomaly still differ from the true model. The resistivity recovery of the middle layer remains weak, and no clear resistivity contrast is formed in the interface-adjacent zone. The overall structural relationship between the anomaly and the background in this zone is therefore still not adequately recovered. This indicates that, although the 1D-result-based initial model improves the background resistivity trend, its effect on the recovery of interface geometry near cross-layer anomalies remains limited. It is therefore necessary to further introduce horizon geometry to strengthen the recovery of such relationships.

Figure 10a shows the constructed horizon-constrained initial model, in which the blue curves represent the 1D inversion results and the red curves represent the synthetic model stratigraphy. For the buried geometry of the overburden interface, the 1D inversion results capture the vertical resistivity trend, but the interface position, thickness, and transition relationship with the underlying background are blurred. After horizon information is introduced, the initial model preserves the resistivity trend while strengthening the interface geometry, so that the spatial position and thickness of the middle layer are more clearly represented and the overall background structural framework becomes more reasonable.

FIGURE 10

As shown by the RMS convergence curve in Figure 10b, compared with the initial model constructed solely from the 1D results, the horizon-constrained initial model not only yields a lower initial RMS at the start of inversion, but also maintains a faster rate of decrease during the first few iterations, finally reaching an acceptable data misfit level after 38 iterations. Figures 10c–f present the final recovered results obtained with the horizon-constrained initial model. Compared with the inversion result based on the 1D-result-based initial model, the improvement brought by the horizon-constrained strategy is first reflected in stronger recovery near the interfaces. The boundaries of the upper anomaly become sharper and the outward diffusion of low resistivity is weakened; the lower anomaly is further better constrained, and the lateral extent and geometric spread of the anomaly recovery become closer to the true model. At the same time, the recovery in interface-adjacent zones becomes clearer, and the mixed transition zone between the background structure and the target anomaly is reduced.

3.4.2.3 Undulating layered model

Based on the regular layered-background model, an undulating layered model (Model 4) is further introduced to examine more directly the influence of topographic relief and interface undulation on inversion results and to evaluate the effect of the proposed strategy under more complex background conditions. As shown in Figure 11a, all physical parameters, source layout, and observation system remain the same as those of Model 3, except that the overburden and its underlying interface are given an undulating geometry.

FIGURE 11

FIGURE 12

Using the homogeneous half-space initial model and the horizon-constrained initial model, two inversions are carried out for Model 4. Figure 11b compares the RMS convergence curves of the two inversion schemes. Under the homogeneous half-space initial model, the initial RMS is higher than that of Model 3 and decreases to 0.96 after 63 iterations, but the overall reduction is slower and the inversion converges less efficiently, indicating lower inversion efficiency. Under the horizon-constrained initial model, the initial RMS is markedly reduced, and an acceptable data misfit level is reached after 45 iterations.

Figures 11c–e compare the inversion results obtained with the homogeneous half-space initial model. A relatively distinct low-resistivity anomaly is recovered at the target location, and its spatial position and overall extent broadly correspond to the true model. The main discrepancies in the recovered result are concentrated in the overall spatial relationship among the topographic background, stratigraphic interfaces, and the target anomaly. More specifically, the recovery of stratigraphic layering remains weak, the transition zone associated with the interface is not clearly delineated, and the topographic undulation and interface geometry are not adequately recovered. In contrast, the inversion results obtained with the horizon-constrained initial model in Figures 11d,e show that the improvement is first reflected in the representation of background structure and stratigraphic relationships. The shallow high-resistivity layer, the middle transition zone, and the deeper relatively low-resistivity background become more clearly expressed, background disturbances near the target anomaly are weakened, and the mixed relationship among the background structure, stratigraphic interfaces, and the target anomaly is reduced.

3.4.3 Analysis of optimization performance and method evaluation

To compare more clearly the improvement in inversion performance achieved by the initial-model optimization strategy, a quantitative evaluation is carried out, in addition to the qualitative analysis presented above, from two aspects: target-anomaly recovery and background-structure recovery.

Let denote the logarithmic resistivity of the recovered cell, and let denote the logarithmic resistivity of the anomaly in the true model. Using the true anomaly extent expanded outward by 100 m as the search region, the recovered cells satisfying are identified. To distinguish between recovery of the anomaly core and that of its margins, a Gaussian weight based on resistivity similarity is introduced as Equation 48:where is the weighting parameter, taken as 0.15 in this study. With this choice, recovered cells with logarithmic resistivity within of the true value retain relatively high weights, whereas the weights decrease markedly for cells deviating by or more. In this way, the evaluation of anomaly position and volume places greater emphasis on the recovered portions that are closer to the true value.

Based on the above definition, the weighted center-position error of the recovered anomaly is defined as Equation 49:

The geometrically weighted center of the recovered anomaly is calculated using Equation 50:where is the geometrically weighted center of the recovered anomaly, is the geometric center of the true anomaly body, and and are the volume and center coordinates of the cell, respectively.

To compare anomaly-volume recovery, the weighted concentration is further defined as Equation 51:

For background-structure recovery, the root-mean-square error of logarithmic resistivity within the background region is used to evaluate the overall recovery effect, as defined in Equation 52:where and denote the recovered and true values of the background cells, respectively, and is the number of cells included in the statistics.

A quantitative comparison of the inversion results based on these metrics is summarized in Table 1. Taken together, the quantitative metrics in Table 1 show that initial-model optimization produces a clear improvement in 3D LOTEM inversion. From the perspective of computational efficiency, compared with conventional inversion, the optimized strategy significantly reduces the number of iterations and still provides substantial savings in overall computation time after accounting for the cost of initial-model construction. From the perspective of recovery quality, for Model 2 the background resistivity trend is represented more reasonably, the target anomaly is recovered in a more concentrated manner, the anomaly position error is markedly reduced, and peripheral low-resistivity spreading is suppressed. For Model 3, the initial model constructed from 1D inversion results enhances concentrated recovery of the anomaly near the target zone, while the introduction of horizon constraints further improves the spatial recovery of the background structure and yields clearer recovery of the cross-layer anomaly and the disturbed interval between layers. For Model 4, the horizon-constrained initial model produces a more pronounced improvement in background structure and stratigraphic relationships, and also makes anomaly recovery more concentrated around the target zone.

These results indicate that the proposed initial-model optimization strategy improves the model basis before inversion by making use of prior information, and that its effectiveness depends closely on the quality of the prior information and the way it is used. In the synthetic examples, the background resistivity trend can be obtained from station-wise 1D inversion results, and the horizon information can be extracted directly from the model parameters; therefore, the initial-model construction is based on complete and ideal prior conditions. For field data, however, these two types of information usually need to be supplemented by integrating station-wise 1D inversion results, well-log data, seismic horizon interpretation, and existing geological or electromagnetic interpretations. Horizon relationships can then be established jointly from control wells, seismic horizons, and regional resistivity models. Compared with methods that introduce prior constraints directly into the inversion, the present strategy uses prior information for initial-model construction, and is therefore better suited to practical datasets in which prior information varies in source, completeness, and reliability. On the one hand, the reasonableness of the prior information directly affects how well the initial model represents the background structure. On the other hand, because prior information does not continue to participate in model updating during iteration, the requirements on its stability and constraint strength are relatively lower, making the strategy more suitable for practical applications.

Overall, the proposed initial-model optimization strategy is well suited to 3D LOTEM inversion under complex background conditions, because it improves the model basis before inversion and provides a more reasonable starting point for subsequent recovery.

4 Inversion analysis of field LOTEM data for deep karst detection

To investigate deep concealed karst in the study area, LOTEM field surveys were carried out in a shale gas block in southwestern China. In this region, carbonate strata are continuously developed, and under humid climatic conditions karstification is active. Structural activity and fluvial incision have produced a complex karst system dominated by fractures, dissolution pores, and karst cavities. In the target interval, karst development is closely related to the distribution of carbonate strata, fracture activity, and dissolution intensity, and its spatial distribution is highly irregular, which may lead to drilling problems such as sticking, bit dropping, and deviation. Therefore, targeted geophysical investigation of deep concealed karst can provide a basis for identifying its spatial distribution, assessing karst-related engineering risk, and optimizing well deployment, and thus has clear geological and engineering significance.

For this target, a LOTEM survey was conducted in the study area. The transmitter was arranged approximately east-west, with a total length of 2,500 m. The survey layout and geoelectrical background of the study area are shown in Figure 12. A total of 12 transmitting stations, 11 survey lines, and 231 measurement points were completed. The station spacing was 40 m, the line spacing was 60 m, and the component parallel to the transmitter direction was recorded.

Following the initial-model optimization strategy proposed above, 1D inversion was first performed on the field data. Representative results for selected stations are shown in Figure 13a. The 1D inversion results clearly reflect the overall variation trend of background resistivity with depth in the study area, but pronounced local fluctuations are observed among different stations, as well as anomalous jumps at zero depth and discontinuities between survey lines. These features are mainly related to topographic relief and near-surface interference. If the 1D inversion results were mapped directly to the 3D mesh to construct the initial model, the correspondence between layers would become unclear, and local resistivity jumps and background-structure fluctuations near the target interval would be amplified. To address this issue, data from a parameter well at the western edge of the study area and seismic horizon interpretations constrained by several wells passing through the target interval and the central part of the area were incorporated (Figure 13b). Using the parameter-well resistivity log and the seismic interpretation, a reference depth trend constrained by horizon geometry was constructed for each stratigraphic interval, and then compared with the 1D inversion results station by station. In this way, the resistivity layering between stations was kept consistent with the stratigraphic interfaces, ensuring that the same layer maintained a corresponding spatial relationship among different stations. Subsequently, in-layer smoothing was applied to the non-targeted conductive anomalies, fluctuations, and discontinuities superimposed on the 1D inversion results at the corresponding horizons, thereby weakening local fragments and discontinuous features and forming a resistivity-depth relationship consistent with the stratigraphic framework. The final initial model is shown in Figure 13e.

FIGURE 13

The inversion converged after 182 iterations. As shown in Figure 14a, the RMS decreased gradually from an initial value of 199 and became stable near 7.9, where the inversion was terminated. Compared with the synthetic examples above, the field-data inversion starts from a much higher RMS and converges more slowly, indicating higher noise levels and a more complex subsurface background. In the study area, the background resistivity structure shows a distinct layered pattern: the shallow section is dominated by high-resistivity background, and the local high-resistivity layer has relatively good continuity; with increasing depth, the background resistivity gradually transitions into a middle-resistivity interval. A deep low-resistivity anomaly is identified in the northern part of the study area at depths of 1,400–1900 m, mainly near the target interval, and is recovered with relatively good continuity. The anomaly is shallower and broader and gradually narrows with depth. The overlying background disturbance is weak, whereas the surrounding resistivity background at depth becomes more convergent. To further examine its spatial distribution, horizontal slices at different depths are shown in Figure 14c. The anomaly begins to appear at about 1,400 m depth, mainly between Lines L07 and L09, with an areal extent of about 300 m × 200 m. As depth increases to about 1,650 m, the anomaly becomes smaller and shifts northwestward. At about 1900 m depth, no distinct low-resistivity anomaly is observed on the slice.

FIGURE 14

Combined with the structural setting of the study area, the pattern of karst development, and well-seismic constraints, the deep low-resistivity anomaly has a relatively clear geological basis. The anomaly is mainly located between Lines L04 and L09, close to the target interval, and its depth range corresponds well to the development interval of the carbonate sequence as well as to local fracture and fault zones. According to regional structural evolution and karst development characteristics, the occurrence of water-bearing fractures and karst cavities is closely related to dissolution along fractures and fracture-related pathways. Carbonate strata provide the material basis for karst development, whereas fracture damage zones and their surrounding fractured intervals are favorable for the formation of water-conducting channels and thus for the development of karst. Therefore, under the combined influence of stratigraphy, structure, and hydrogeological conditions, this anomaly is expressed in the resistivity data as a deep low-resistivity feature. The anomaly maintains relatively good continuity on slices at different depths and is characterized overall by a broad upper part and a gradually convergent lower part. Considering its spatial distribution, resistivity characteristics, and its neighboring relationship with the fracture damage zone, this deep low-resistivity anomaly is inferred to represent a deep karst-developed zone in the study area and can provide a geophysical basis for drilling deployment.

5 Conclusion

Overall, the proposed initial-model optimization strategy improves inversion recovery under complex background conditions by using prior information in initial-model construction, thereby making the representation of background structure more reasonable. The strategy is applicable to practical datasets with different levels of prior information. When prior information is more complete, the background structure in the initial model is represented more reasonably, and the improvement in inversion results becomes more pronounced.

References

• 1

  CaiH.LiuM.ZhouJ.LiJ.HuX. (2022). Effective 3D transient electromagnetic inversion using finite-element method with a parallel direct solver. Geophysics87, E377–E392. 10.1190/geo2021-0630.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  CarterS.XieX.ZhouL.YanL. (2021). Hybrid monte carlo 1-D joint inversion of LOTEM and MT. J. Appl. Geophys.194, 104424. 10.1016/j.jappgeo.2021.104424

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  ChengM.YangD. K.LuoQ. (2023). Interpreting surface large-loop time-domain electromagnetic data for deep mineral exploration using 3D forward modeling and inversion. Minerals13, 34. 10.3390/min13010034

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  CommerM.NewmanG. A. (2004). A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources. Geophysics69, 1192–1202. 10.1190/1.1801936

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  CommerM.HelwigS. L.HördtA.SchollC.TezkanB. (2006). New results on the resistivity structure of merapi volcano (indonesia) derived from three-dimensional restricted inversion of long-offset transient electromagnetic data. Geophys. J. Int.167, 1172–1187. 10.1111/j.1365-246x.2006.03182.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  HaberE.OldenburgD. W.ShekhtmanR. (2007). Inversion of time-domain three-dimensional electromagnetic data. Geophys. J. Int.171, 550–564. 10.1111/j.1365-246x.2007.03365.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  HaroonA.AdrianJ.BergersR.GurkM.TezkanB.MammadovA. L.et al (2015). Joint inversion of long-offset and central-loop transient electromagnetic data: application to a mud volcano exploration in perekishkul, Azerbaijan. Geophys. Prospect.63, 478–494. 10.1111/1365-2478.12157

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  HördtA.DruskinV. L.KnizhnermanL. A.StrackK. M. (1992). Interpretation of 3-D effects in long-offset transient electromagnetic (LOTEM) soundings in the münsterland area, Germany. Geophysics57, 1127–1137. 10.1190/1.1443327

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  HouD.XueG.ZhouN.YanS. (2017). “The shielding effect of low resistivity layer in TEM,” in Technology and Application of Environmental and Engineering Geophysics (Singapore: Springer), 145–150.

  • Google Scholar
  • View reference in article

• 10

  HuZ. (2021). Three-Dimensional Forward and Inverse Modeling of Marine Time-Domain Electromagnetics Based on Unstructured Finite Elements. Changchun, China: Jilin University. 10.27162/d.cnki.gjlin.2021.000478

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  KellerG. V.PritchardJ. I.JacobsonJ. J.HarthillN. (1984). Megasource time-domain electromagnetic sounding methods. Geophysics49, 993–1009. 10.1190/1.1441743

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  KhanM. Y.XueG. Q.ChenW. Y.ZhongH. S. (2018). Analysis of long-offset transient electromagnetic (LOTEM) data in time, frequency, and pseudo-seismic domain. J. Environ. Eng. Geophys.23, 15–32. 10.2113/jeeg23.1.15

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  LiJ.HuX.CaiH.LiuY. (2020). A finite-element time-domain forward-modeling algorithm for transient electromagnetics excited by grounded-wire sources. Geophys. Prospect.68, 1379–1398. 10.1111/1365-2478.12917

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  LiuY.YinC.QiuC.HuiZ.ZhangB.RenX.et al (2019). 3-D inversion of transient EM data with topography using unstructured tetrahedral grids. Geophys. J. Int.217, 301–318. 10.1093/gji/ggz014

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  MüllerM.HördtA.NeubauerF. M. (2002). Internal structure of mount merapi, Indonesia, derived from long-offset transient electromagnetic data. J. Geophys. Res. Solid Earth107, ECV 2-1–ECV 2-14. 10.1029/2001jb000148

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  NewmanG. A.CommerM. (2005). New advances in three-dimensional transient electromagnetic inversion. Geophys. J. Int.160, 5–32. 10.1111/j.1365-246x.2004.02468.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  OldenburgD. W.HaberE.ShekhtmanR. (2013). Three-dimensional inversion of multisource time-domain electromagnetic data. Geophysics78, E47–E57. 10.1190/geo2012-0131.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  QiY.LiX.YinC.LiH.QiZ.ZhouJ.et al (2020). 3-D time-domain airborne EM inversion for a topographic Earth. IEEE Trans. Geoscience Remote Sens.60, 2000113. 10.1109/tgrs.2020.3036084

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  RenX.YinC.MacnaeJ. Z. B.LiuY.ZhangB. (2018). 3D time-domain airborne electromagnetic inversion based on secondary field finite-volume method. Geophysics83, E219–E228. 10.1190/geo2017-0585.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  Rezaei MirghahedB.Dehghan MonfaredA.RanjbarA. (2024). Enhanced petrophysical evaluation through machine learning and well-logging data in an Iranian oil field. Sci. Rep.14, 28941. 10.1038/s41598-024-80362-w

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  SchollC. (2005). The Influence of Multidimensional Structures on the Interpretation of LOTEM Data with One-Dimensional Models and the Application to Data from Israel. Cologne, Germany: Universität zu Köln.

  • Google Scholar
  • View reference in article

• 22

  StephanA.SchniggenfittigH.StrackK. M. (1991). Long-offset transient EM sounding north of the rhine–ruhr coal district, Germany. Geophys. Prospect.39, 505–525. 10.1111/j.1365-2478.1991.tb00325.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  StrackK. M. (1992). Exploration with Deep Transient Electromagnetics. Amsterdam, Netherlands: Elsevier.

  • Google Scholar
  • View reference in article

• 24

  StrackK. M.LüschenE.KötzA. W. (1990). Long-offset transient electromagnetic (LOTEM) depth soundings applied to crustal studies in the black forest and swabian alb, Federal Republic of Germany. Geophysics55, 834–842. 10.1190/1.1442897

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  SunX.WangY.YangX.WangY. (2021). Three-dimensional transient electromagnetic inversion with optimal transport. J. Inverse Ill-Posed Problems30, 549–565. 10.1515/jiip-2020-0159

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  TangX. G.HuW. B.YanL. J. (2011). Topographic effects on long-offset transient electromagnetic response. Appl. Geophys.8, 277–284. 10.1007/s11770-011-0297-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  UmE. S. (2011). Three-Dimensional Finite-Element Time-Domain Modeling of the Marine Controlled-Source Electromagnetic Method. Stanford, CA, USA: Stanford University.

  • Google Scholar
  • View reference in article

• 28

  WangT.HohmannG. W. (1993). A finite-difference time-domain solution for three-dimensional electromagnetic modeling. Geophysics58, 797–809. 10.1190/1.1443465

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  WangX.CaiH.LiuL.RevilA.HuX. (2023). Three-dimensional inversion of long-offset transient electromagnetic method over topography. Minerals13, 908. 10.3390/min13070908

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  XieX. B.ZhouL.YanL. J.HuW. B. (2016). Remaining oil detection with time-lapse long offset & window transient electromagnetic sounding. Oil Geophys. Prospect.51, 605–612. 10.13810/j.cnki.issn.1000-7210.2016.03.024

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  XueG.ZhouN.WangR.LiuH.GuoW. (2021). Exploration of lead–zinc deposits using electromagnetic method: a case study in fengtai ore deposits in western China. Geol. J.56, 3314–3321. 10.1002/gj.4103

  • CrossRef
  • Google Scholar

• 32

  XueG.ChenW.WuX.YanS.GuoW. (2022). A near-source electromagnetic method for deep ore explorations. Minerals12, 1208. 10.3390/min12101208

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  YanL. J.ChenX. X.TangH.XieX. B.ZhouL.HuW. B.et al (2018). Continuous TDEM for monitoring shale hydraulic fracturing. Appl. Geophys.15, 26–34. 10.1007/s11770-018-0661-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  YangD.OldenburgD. W.HaberE. (2014). 3D inversion of airborne electromagnetic data parallelized and accelerated by local mesh and adaptive soundings. Geophys. J. Int.196, 1492–1507. 10.1093/gji/ggt465

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  YeT.ChenX. B.YanL. J. (2013). Refined techniques for data processing and two-dimensional inversion in magnetotelluric (III): using the impressing method to construct starting model of 2D magnetotelluric inversion. Chin. J. Geophys.56, 3596–3606. 10.6038/cjg20131034

  • CrossRef
  • Google Scholar
  • View reference in article