A New Approximation Modeling Method for the Triaxial Induction Logging in Planar-Stratified Biaxial Anisotropic Formations

Abstract

A novel and efficient modeling approach has been developed for simulating the responses of triaxial induction logging (TIL) in layered biaxial anisotropic (BA) formations. The core of this innovative technique lies in analytically calculating the primary fields within a homogeneous medium and approximating the scattered fields within layered formations. The former involves employing a two-level subtraction technique. Initially, the first-level subtraction entails altering the direction of the Fourier transform to mitigate the integral singularity of the spectral fields, particularly in high-angle and horizontal wells. Conversely, the second-level subtraction aims to further optimize integral convergence by creating an equivalent unbounded transverse isotropic (TI) formation and eliminating the corresponding spectral fields. With the two-level subtractions, the convergence of the spectral field has been enhanced by more than six orders of magnitude. Additionally, a strict recursive algorithm and approximation method are developed to compute the scattered fields in layered biaxial anisotropic media. The rigorous algorithm is based on a modified amplitude propagator matrix (MAPM) approach and serves as the benchmark for the approximation method. In contrast, the approximation method exploits the similarity between the spectral scattered field of the TI medium and the BA medium, establishing corresponding equivalent layered TI models for each magnetic component. Since the scattered field in TI models only involves a one-dimensional semi-infinite integral, the computational complexity is significantly reduced. Numerical simulation examples demonstrate that the new simulation method is at least two orders of magnitude faster than the current modeling approach while maintaining computational precision error within 0.5%. This significantly improved simulation efficiency provides a solid foundation for expediting the logging data processing.

1. Introduction

Formation anisotropy is prevalent in deep-water reservoirs and can be classified into micro-anisotropy and macro-anisotropy. The former is typically induced by directional grain deposition or micro-fracture development, whereas the latter occurs when the formation is thinner than the vertical resolution of the logging tool [1,2,3]. Accurately measuring formation resistivity around the borehole has always been a significant challenge in geophysical logging due to the influence of formation anisotropy [4,5]. Due to the special characteristics of the deep-water reservoir structure, high angle (HA) and horizontal (HZ) wells have become the main drilling method [6]. Existing electrical logging tools are primarily designed for isotropic formations, resulting in biased responses and inaccurate readings in HA and HZ wells with anisotropy. The triaxial induction logging (TIL) technique was developed to address the need to detect resistivities in different directions within the formation [7,8]. TIL measures the tensor magnetic field and provides rich formation information [9]. However, the responses of TIL in deviated wells are often non-linear, necessitating the use of fast inversion techniques for qualitative interpretation and quantitative evaluation of these data [10,11]. Existing inversion methods for TIL typically utilize parametric regularization gradient optimization algorithms, which require extensive forward modelling to compute the Jacobian matrix. Therefore, a fast and accurate forward simulator is essential for the rapid reconstruction of anisotropic resistivity from complex TIL data, considering its rich information content [12,13,14].

The fractured shale or tight sand formations can be depicted by a biaxial anisotropic (BA) model, where the conductivity varies in x, y and z directions. Considering the complexity of the forward model, the simulation of TIL in BA formations can be categorized into three types: numerical methods, hybrid methods and analytical methods [15,16]. The numerical and hybrid methods offer better accommodation with the complex borehole and formation structure. However, due to their low computational efficiency, these algorithms are seldom utilized in the practical logging data processing. Neglecting the effects of borehole, invasion and coil size, the analytical methods are commonly used to obtain TIL responses in one-dimensional (1D) media. The foundation of analytical methods lies in describing the propagation of electromagnetic (EM) waves in planar-stratified BA media. To this end, the generalized reflection matrix method and the coefficient propagation matrix method are developed [17,18,19]. The former involves complex formula derivation, while the latter often suffers from the issue of numerical overflow. Consequently, there still lacks a concise and stable algorithm for calculating the magnetic fields in layered BA media.

For the forward modeling of TIL in 1D formation, efficiently and accurately converting EM fields from the spectral domain to the spatial domain is a challenge. This conversion, also known as the inverse Fourier transform (IFT) or Sommerfeld integral, is significantly impacted by the integrand convergence [20,21,22]. IFT methods can be classified into the direct integral method and digital filter integral method. The former, e.g., Gauss–Laguerre integral and cubic spline interpolation, can offer controllable computation precision, albeit requiring excessive nodes for HA wells with slowly converging integrals. In contrast, the latter, such as the fast Hankel transform algorithm, can effectively alleviate slow integrand convergence, and has been widely used in the geophysical fields [23,24]. Unfortunately, both methods are inadequate for handling the case of HZ well case where the integral is non-convergent [25,26]. Hu et al. (2018) addressed this singularity by separating the primary fields from the total fields [27]. Similarly, Li et al. (2020) subtracted part of the primary spectral fields by constructing an equivalent transversely isotropy (TI) model [28]. The periodicity of kernel function is another factor dominating the modeling speed of TIL in BA medium. In homogeneous BA formation, Deng et al. (2020) utilized the field symmetry at four integral quadrants, giving rise to the computation speed improved by three times [29]. Later, Wang et al. (2023b) extended above property to a layered BA media [30]. However, the existing forward modeling techniques of TIL in layered BA medium are still limited in the efficient processing of the double IFT, and unable to meet the needs of real time data processing.

First, this paper introduces the principle of TIL, which explains the relationship between apparent conductivities and measured magnetic fields. Then, a two-level singularity subtraction scheme is established to significantly improve the convergence of the spectral primary field in homogeneous media. In Section 3, both the strict derivation and the approximation method for calculating the spectral scattered field in layered BA media are presented, and the accuracy and efficiency of the approximation method are verified. Finally, the accuracy and stability of the new forward modeling method are validated, and the characteristics of TIL responses under different anisotropy are discussed.

2. Physics of Triaxial Induction Logging

The TIL tool usually consists of a transmitter (T), bucking (B) coils and a receiver (R), all of which are composed of three mutually orthogonal coils oriented at the x, y and z directions, respectively. The specific structure and antenna configuration of the TIL tool are shown in Figure 1a. These tools typically operate within a frequency of 10 kHz to 200 kHz, with negligible impact from the coil size. TIL tools are capable of measuring the formation tensor voltage $\widehat{\mathit{V}}$,

$$\widehat{V}=\left[\begin{array}{ccc}{V}_{xx}& V_{xy}& V_{xz}\\ V_{yx}& V_{yy}& V_{yz}\\ V_{zx}& V_{zy}& V_{zz}\end{array}\right],$$

(1)

where V_ij signifies the voltage component of j-oriented receiver emitted from an i-oriented transmitter. The voltage tensor $\widehat{V}$ is further converted into formation apparent conductivity σ_a by utilizing the tool coefficient K,

$${\sigma }_a=\left[\begin{array}{ccc}{\sigma }_a^{xx}& {\sigma }_a^{xy}& {\sigma }_a^{xz}\\ {\sigma }_a^{yx}& {\sigma }_a^{yy}& {\sigma }_a^{yz}\\ {\sigma }_a^{zx}& {\sigma }_a^{zy}& {\sigma }_a^{zz}\end{array}\right]=-{\displaystyle \frac{\widehat{V}}{K}},$$

(2)

$$K={\displaystyle \frac{{\left(\omega {\mu }_0\right)}^2}{16\pi L}}\left[\begin{array}{ccc}2A_{T_x}{A}_{R_x}{N}_{T_x}{N}_{R_x}& 2A_{T_x}{A}_{R_y}{N}_{T_x}{N}_{R_y}& 4A_{T_x}{A}_{R_z}{N}_{T_x}{N}_{R_z}\\ 2A_{T_y}{A}_{R_x}{N}_{T_y}{N}_{R_x}& 2A_{T_y}{A}_{R_y}{N}_{T_y}{N}_{R_y}& 4A_{T_y}{A}_{R_z}{N}_{T_y}{N}_{R_z}\\ 4A_{T_z}{A}_{R_x}{N}_{T_z}{N}_{R_x}& 4A_{T_z}{A}_{R_y}{N}_{T_z}{N}_{R_y}& A_{T_z}{A}_{R_z}{N}_{T_z}{N}_{R_z}\end{array}\right],$$

(3)

where μ₀ is the free-space magnetic permittivity, and ω is the angular frequency; L represents the spacing between the T and R; A_T/A_R and N_T/N_R represent the area and turns of the T/R, respectively.

The essence of TIL lies in measuring the magnetic field tensor $\overline{H}$, $\widehat{V}=i\omega \mu A_T{A}_R{N}_T{N}_R\overline{H}$. The measurement of $\overline{H}$ in the tool coordinate system (X, Y, Z) not only depends on the spacing and operating frequency, but also relies on the geometric relationship between the TIL tool and the formation. As depicted in Figure 1b, the direction of the TIL tool is arbitrary put with respect to the principal axis of magnetic tensor H in the formation coordinate system (x, y, z). The magnetic field tensor at two coordinate systems is connected by a rotation matrix R,

$$\overline{H}=R^{-1}HR,R=\left[\begin{array}{ccc}cos\alpha cos\beta & -cos\alpha cos\beta & sin\alpha cos\beta \\ cos\alpha sin\beta & cos\beta & sin\alpha sin\beta \\ 0& 0& cos\alpha \end{array}\right],$$

(4)

where α denotes the dipping angle between the tool and the normal of formation. β represents the azimuthal angle of the tool relative to formation.

3. Computation of Primary TIL Response in Homogeneous Formation

3.1. Spectral Primary Field in Homogeneous BA Medium

Assuming a harmonic source with e^-iωt dependence, the Maxwell’s equations satisfy,

$$\nabla \times E\left(r\right)=i\omega {\mu }_0\left[H\left(r\right)+M\cdot \delta \left(r\right)\right],\nabla \times H\left(r\right)={\sigma }^{*}\cdot E\left(r\right),$$

(5)

where M is the magnetic moment, M = M_xe_x + M_ye_y + M_ze_z. E and H are electric and magnetic fields in the spatial domain, respectively. r (x, y, z) is the source position, and δ(·) is the Dirac function. μ₀ is the free-space magnetic permittivity. σ* is the complex conductivity tensor, which can be expressed as follows

$${\sigma }^{*}=diag\left({\sigma }_x^{*},{\sigma }_y^{*},{\sigma }_z^{*}\right),{\sigma }_i^{*}={\sigma }_i-i\omega {\epsilon }_i,$$

(6)

where σ is the formation conductivity, and Ɛ represents the dielectric constant.

Since the heterogeneity of formations only exists along the z-direction, the spatial EM fields can be converted to the spectral domain via a twofold Fourier transform. Hence, the tangential EM component satisfy the following expression,

$${\displaystyle \frac{\partial }{\partial z}}\tilde{S}\left(k,z\right)+U\tilde{S}\left(k,z\right)=\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right]\delta \left(z-z_a\right),\tilde{S}={\left[\begin{array}{cccc}{\tilde{E}}_x& {\tilde{E}}_y& -{\tilde{H}}_y& {\tilde{H}}_x\end{array}\right]}^T,$$

(7)

where the matrix $\tilde{S}$ is composed of spectral EM components to be solved, and k_x and k_y correspond to the wave numbers along the x- and y-axes, respectively. z_a signifies the vertical position of the dipole source. U is a 4 × 4 system matrix $U=\left[\begin{array}{cc}0& u_1\\ u_2& 0\end{array}\right]$, and d₁ and d₂ are the 2 × 1 source vector. The system matrix has four eigenvalues (λ_I, λ_II, −λ_I, −λ_II), and the corresponding eigenvectors form the eigenmatrix C, $C=\left[\begin{array}{cc}{c}_{11}& c_{11}\\ c_{21}& -c_{21}\end{array}\right]$, where

$$\begin{array}{l}{u}_1=\left[\begin{array}{cc}i\omega \mu -{\displaystyle \frac{k_x^2}{{\sigma }_z^{*}}}& -{\displaystyle \frac{k_x{k}_y}{{\sigma }_z^{*}}}\\ -{\displaystyle \frac{k_x{k}_y}{{\sigma }_z^{*}}}& i\omega \mu -{\displaystyle \frac{k_y^2}{{\sigma }_z^{*}}}\end{array}\right],d_1=\left[\begin{array}{c}i\omega \mu M_y\\ -i\omega \mu M_x\end{array}\right],\\ u_2=\left[\begin{array}{cc}-{\sigma }_x^{*}-{\displaystyle \frac{i{k}_y^2}{\omega \mu }}& {\displaystyle \frac{i{k}_x{k}_y}{\omega \mu }}\\ {\displaystyle \frac{i{k}_x{k}_y}{\omega \mu }}& -{\sigma }_y^{*}-{\displaystyle \frac{i{k}_x^2}{\omega \mu }}\end{array}\right],d_2=\left[\begin{array}{c}i{k}_y{M}_z\\ -i{k}_x{M}_z\end{array}\right].\end{array}$$

(8)

Mathematically, the solution of Equation (7) can be written as $\tilde{S}=\beta {\tilde{S}}^1+{\tilde{S}}^2,\beta =\{\begin{array}{c}0,source\text{-}free\hfill \\ 1,source\hfill \end{array}$. ${\tilde{S}}^1$ represents the spectral primary field, and ${\tilde{S}}^2$ denotes the spectral scattered fields. In a homogeneous BA medium, applying the Fourier transform along the z-direction, the spectral primary field can be further written as,

$${\tilde{S}}^1\left(k_x,k_y,k_z\right)={\tilde{S}}^p\left(k_x,k_y,k_z\right)\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right],{\tilde{S}}^p\left(k_x,k_y,k_z\right)={\displaystyle \frac{{\left(U+i{k}_{z}I\right)}^{*}}{\left|U+i{k}_{z}I\right|}},$$

(9)

where k_z corresponds to the wave numbers along the z-axe, and the superscript * demotes the adjoint matrix. Performing the IFT on Equation (9) and then combining with the residue theorem, the analytical expression of the spectral primary fields in homogeneous BA medium can be obtained (Yuan et al., 2010),

$${\tilde{S}}^p(k_x,k_y,z)=\{\begin{array}{c}{\displaystyle \frac{{\left(U-{\lambda }_{I}I\right)}^{*}{e}^{-{\lambda }_I\left(z-z_a\right)}}{2{\lambda }_I\left({\lambda }_{II}^2-{\lambda }_I^2\right)}}-{\displaystyle \frac{{\left(U-{\lambda }_{II}I\right)}^{*}{e}^{-{\lambda }_{II}\left(z-z_a\right)}}{2{\lambda }_{II}\left({\lambda }_{II}^2-{\lambda }_I^2\right)}},z>z_a\\ {\displaystyle \frac{{\left(U+{\lambda }_{I}I\right)}^{*}{e}^{{\lambda }_I\left(z-z_a\right)}}{2{\lambda }_I\left({\lambda }_{II}^2-{\lambda }_I^2\right)}}-{\displaystyle \frac{{\left(U+{\lambda }_{II}I\right)}^{*}{e}^{{\lambda }_{II}\left(z-z_a\right)}}{2{\lambda }_{II}\left({\lambda }_{II}^2-{\lambda }_I^2\right)}},z<z_a\end{array}.$$

(10)

Generally, the contribution of primary fields exceeds 95% in the total field (Løseth and Ursin B, 2007). Hence, achieving a swift and precise primary field calculation is imperative. In order to derive the spatial primary field, it is mandatory to implement the following double IFT,

$$S^1\left(r\right)={\displaystyle \frac{1}{4{\pi }^2}}{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }{\tilde{S}}^1\left(k_x,k_y,z\right)}}{e}^{i{k}_{x}x+i{k}_{y}y}d{k}_{x}d{k}_y.$$

(11)

To mitigate the computational complexity, the spectral EM fields at four quadrants are observed to exhibit symmetric or antisymmetric expressions in the homogeneous medium. Recently, Wang et al. (2023b) extends this symmetry analysis to the layered BA media. Hence, the original infinite integral in Equation (11) is simplified to a semi-infinite one, resulting in a three-quarters reduction in the computational burden for the IFT.

Table 1. Symmetry of spectral magnetic components at four quadrants.

	Components	${\mathit{H}}_{\mathit{x}\mathit{x}}$	${\mathit{H}}_{\mathit{y}\mathit{y}}$	${\mathit{H}}_{\mathit{z}\mathit{z}}$	${\mathit{H}}_{\mathit{x}\mathit{y}},{\mathit{H}}_{\mathit{y}\mathit{x}}$	${\mathit{H}}_{\mathit{x}\mathit{z}},{\mathit{H}}_{\mathit{z}\mathit{x}}$	${\mathit{H}}_{\mathit{y}\mathit{z}},{\mathit{H}}_{\mathit{z}\mathit{y}}$
Quadrants
1st		+1	+1	+1	+1	+1	+1
2nd		−1	−1	−1	+1	+1	−1
3rd		+1	+1	+1	+1	−1	−1
4th		−1	−1	−1	+1	−1	+1

further illustrates the relationships of spectral magnetic fields at different quadrants. The spectral field obtained in the first quadrant is used as the reference (denoted by a plus sign) and compared to the spectral fields in the other quadrants with it, indicating the opposite direction with a minus sign. For instance, when calculating H_xx, the simplified integral takes the following form,

$$H_{xx}(r)={\displaystyle \frac{1}{4{\pi }^2}}{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\infty }{\tilde{H}}_{xx}(k_x,k_y,z)(e_1-e_2+e_3-e_4)d{k}_{x}d{k}_y}},$$

(12)

where ${\tilde{H}}_{xx}(k_x,k_y,z)$ represents the spectral field of xx component in the first quadrant (k_x > 0, k_y > 0), $e_1=e^{i{k}_{x}x+i{k}_{y}y},e_2=e^{-i{k}_{x}x+i{k}_{y}y},e_3=e^{-i{k}_{x}x-i{k}_{y}y},e_4=e^{i{k}_{x}x-i{k}_{y}y}$.

3.2. The Two-Level Subtraction of Spectral Primary Fields

Although the computation of the double IFT can be greatly reduced by exploiting the symmetry of the spectral primary field, the high oscillation and low convergence of the integral kernel usually result in poor accuracy and high computational complexity of the double semi-infinite integral. In order to improve the convergence of the integral kernel function, a two-level subtraction scheme based on a homogenous BA model is introduced in this section. The first-level subtraction is based on the principle of coordinate transform, which determines the direction of the optimal Fourier transform by the position of the receiver. Specifically, the projection lengths of the tool are acquired along the x-, y- and z-axes, and then Fourier transform is performed along the two directions with the smaller length. The second-level subtraction involves constructing an equivalent TI medium, subtracting the contribution of the spectral primary fields of TI medium and evaluating the residual spectral primary fields. The TI medium can be regarded as a special case of the BA medium, i.e., σ_x = σ_y, and the spectral primary fields of the two model have similar convergence properties. Since the spectral primary fields of the TI medium are analytical, the spectral primary fields of homogeneous BA medium can be written as follows,

$$S^1\left(r\right)={\displaystyle \frac{1}{4{\pi }^2}}{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }{\tilde{S}}_{res}^1}{e}^{i\left(l{k}_l+q{k}_q\right)}d{k}_{l}d}{k}_q+S_{TI}^1\left(r\right),{\tilde{S}}_{res}^1={\tilde{S}}^1-{\tilde{S}}_{TI}^1,$$

(13)

where ${\tilde{S}}_{res}^1$ is the residual field in a BA medium after subtracting the primary field of a TI medium. l and q represent the directions of the Fourier transform, while k_l and k_q denote the wave numbers corresponding to l and q.

In order to evaluate the effectiveness of the two-level subtraction scheme, a homogeneous BA model is established with conductivity σ = diag{1, 0.5, 0.25} S/m. Figure 2 shows the TIL responses of the coaxial magnetic component for α = 0°, 30°, 45°, 60°, 75° and 90° with the operation frequency of 10 kHz and a spacing of 0.5 m. The upper floating layer in Figure 2 represents the spectral primary fields without subtraction, while the lower floating layer depicts the spectral primary fields after the first-level subtraction. Evidently, the convergence of the spectral primary fields decreases with increasing α and fails in HZ wells. However, the first-level subtraction scheme notably enhances the convergence of spectral primary fields in HA wells and resolves the non-convergence issue in HZ wells, albeit proving impractical for low-dipping wells (α ≤ 45°).

In contrast, Figure 3 compares the results obtained using the first-level subtraction with those using the second-level subtraction, where the equivalent TI medium of the second-level subtraction has the horizontal conductivity σ_h = sqrt (σ_x* σ_y) and the vertical conductivity σ_v = σ_z. It can be found that the second-level subtraction further enhances the convergence of the spectral primary fields by over four orders of magnitude, making it suitable for any dipping wells. Therefore, the two-level subtraction scheme effectively mitigates the integral singularity and computation complexity, increasing the calculation speed of the primary fields in the homogeneous BA medium by over ten times.

3.3. Optimal Construction of Equivalent TI Medium

To ensure the best convergence enhancement, this subsection explores the optimal construction of equivalent TI models for both coaxial and coplanar magnetic components. Figure 4 shows the convergence of the spectral magnetic components using different equivalent TI models in the vertical well, where k_ρ = sqrt ($k_x^2$ + $k_y^2$). The black solid line signifies the convergence of the spectral primary field after the first-level subtraction, while the other lines represent the convergence of the residual fields after the second-level subtraction using different equivalent models. Here, four equivalent TI media are considered, all models sharing identical vertical conductivity, σ_v = σ_z. In contrast, the horizontal conductivities of four models are set as σ_h = sqrt (σ_x* σ_y), σ_h = σ_y, σ_h = σ_x and σ_h = sqrt (${\sigma }_x^2$ + ${\sigma }_y^2$). Additionally, two scenarios of BA cases (σ_x > σ_y and σ_x= < σ_y) are considered. Figure 4a–c uses a homogeneous BA medium with conductivity of [σ_x, σ_y, σ_z] = [1, 0.5, 0.25] S/m, while Figure 4d–f is computed with [σ_x, σ_y, σ_z] = [0.5, 1, 0.25] S/m. Taking the ${\tilde{H}}_{xx}$ component as an example, it is observed that employing four distinct equivalent TI model constructions yields comparable results for the scenario where σ_x > σ_y, with the approach of σ_h = σ_y slightly outperforming the other three methods. This trend presents in the scenario where σ_x < σ_y. Likewise, the optimal TI model construction for the ${\tilde{H}}_{yy}$ and ${\tilde{H}}_{zz}$ components are σ_h = σ_x and σ_h = sqrt (σ_x* σ_y), respectively.

To further clarify the effectiveness of the aforementioned equivalent model, Figure 5 illustrates the convergence of magnetic components after the second-level subtraction in HZ wells. In contrast to the vertical wells, notable disparities emerge in the effect of second-level subtraction when constructing different equivalent TI models for the HZ wells. The subtracted spectral fields with the optimal TI models in the vertical well converge much faster. Therefore, the optimal equivalent TI model for coaxial and coplanar components in HZ wells remains consistent with that in vertical wells, regardless of variations in the dip angle.

4. Rapid Computation of Scattered Fields in Layered BA Formations

4.1. The Modified Amplitude Propagator Matrix (MAPM) Algorithm

In contrast to the primary fields, analytical expression for the scattered fields in layered BA formations are unavailable, posing significant challenges for the fast computations of TIL responses. To drive the scattered fields, a mandatory approach is a recursive algorithm, such as the amplitude propagator matrix algorithm. In layered BA formation, the scattered fields at the j-th layer can be formulated as follows,

$${\tilde{S}}_j^2=C_j\cdot \left[\begin{array}{cc}{e}^{-{\Lambda }_j\left(z-z_{j-1}\right)}& 0\\ 0& e^{{\Lambda }_j\left(z-z_j\right)}\end{array}\right]\cdot \left[\begin{array}{c}{v}_j\\ w_j\end{array}\right],$$

(14)

where

$$v_j={\left[v_{I,j},v_{II,j}\right]}^T,w_j={\left[w_{I,j},w_{II,j}\right]}^T,{\Lambda }_j=diag\left\{{\lambda }_{I,j},{\lambda }_{II,j}\right\},C_j=\left[\begin{array}{cc}{c}_{11,j}& c_{11,j}\\ c_{21,j}& -c_{21,j}\end{array}\right].$$

(15)

assuming that the number of layers in formation is n, and the propagation of EM waves in a layered medium is shown in Figure 6. In theory, v_j and w_j are the amplitude of upgoing and downgoing waves at the lower and upper boundaries of the j-th layer, respectively. Since there is no upgoing wave in the first layer, v_I,1 = v_II,1 = 0. Similarly, there is no downgoing wave in the last layer, w_n = 0.

Considering the continuity condition satisfied by the amplitude of upgoing/downgoing waves at each boundary, the propagation of EM waves in a layered medium can be obtained. Take the boundary at z = z_j as an example, the continuous tangential EM fields of adjacent layers are formulated as follows

$$F_j\cdot \left[\begin{array}{c}{v}_j\\ w_j\end{array}\right]-G_{j+1}\cdot \left[\begin{array}{c}{v}_{j+1}\\ w_{j+1}\end{array}\right]=\left[{\zeta }_{j+1}{\tilde{S}}_{z_j\ge z_a}^p\left(z_j\right)-{\zeta }_j{\tilde{S}}_{z_j<z_a}^p\left(z_j\right)\right]\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right],$$

(16)

where if the source embedded in the i-th layer, ${\zeta }_i=1$ otherwise ${\zeta }_i=0$. F_j and G_j+1 are written as follows,

$$F_j=C_j\cdot \left[\begin{array}{cc}{e}^{-{\Lambda }_j\left(z_j-z_{j-1}\right)}& \\ & 1\end{array}\right],G_{j+1}=C_{j+1}\cdot \left[\begin{array}{cc}1& \\ & e^{{\Lambda }_{j+1}\left(z_j-z_{j+1}\right)}\end{array}\right].$$

(17)

In the traditional coefficient propagator matrix algorithm, the amplitudes of upgoing and downgoing waves in each layer are typically treated as a propagator pair, [v_j, w_j]^T. When solving for the amplitude matrix, we need to obtain the inverse of matrices F and G. Unfortunately, problems arise when λ takes on a large value, such as in the case of thick layers or high wave numbers. In such situations, F⁻¹ and G⁻¹ tend to be ill-conditioned, leading to a numerical overflow issue.

To accurately stably and stably determine the unknown wave amplitudes, a modified amplitude propagator matrix (MAPM) algorithm is proposed. This algorithm establishes the relationship between three constructive amplitude terms of upgoing and downgoing waves at three adjacent boundaries. Hence, Equation (16) is written as follows,

$$\begin{array}{l}{R}_{2j-1}{x}_{2j-2}+Q_{2j-1}{x}_{2j-1}+T_{2j-1}{x}_{2j}=o_{2j-1}\\ R_{2j}{x}_{2j-1}+Q_{2j}{x}_{2j}+T_{2j}{x}_{2j+1}=o_{2j}\end{array}$$

(18)

where x_2j−1 = w_j, x_2j = v_j+1, and

$$\{\begin{array}{c}{R}_{2j-1}=\left(c_{11,j}{c}_{21,j+1}+c_{11,j+1}{c}_{21,j}\right)e^{-{\Lambda }_j\left(z_j-z_{j-1}\right)}\\ Q_{2j-1}=c_{11,j}{c}_{21,j+1}-c_{11,j+1}{c}_{21,j}\\ T_{2j-1}=-c_{11,j+1}{c}_{21,j+1}-c_{11,j+1}{c}_{21,j+1}\\ \begin{array}{l}{o}_{2j-1}=c_{21,j+1}\left[{\zeta }_{j+1}{\tilde{S}}_{z_j\ge z_a}^p\left(z_j\right)-{\zeta }_j{\tilde{S}}_{z_j<z_a}^p\left(z_j\right)\right]d_1\\ +c_{11,j+1}\left[{\zeta }_{j+1}{\tilde{S}}_{z_j\ge z_a}^p\left(z_j\right)-{\zeta }_j{\tilde{S}}_{z_j<z_a}^p\left(z_j\right)\right]d_2\end{array}\end{array},\{\begin{array}{c}{R}_{2j}=c_{11,j}{c}_{21,j}+c_{11,j}{c}_{21,j}\\ Q_{2j}=c_{11,j}{c}_{21,j+1}-c_{11,j+1}{c}_{21,j}\\ T_{2j}=\left(-c_{11,j+1}{c}_{21,j}-c_{11,j}{c}_{21,j+1}\right)e^{{\Lambda }_{j+1}\left(z_j-z_{j+1}\right)}\\ \begin{array}{l}{o}_{2j-1}=c_{21,j}\left[{\zeta }_{j+1}{\tilde{S}}_{z_j\ge z_a}^p\left(z_j\right)-{\zeta }_j{\tilde{S}}_{z_j<z_a}^p\left(z_j\right)\right]d_1\\ -c_{11,j}\left[{\zeta }_{j+1}{\tilde{S}}_{z_j\ge z_a}^p\left(z_j\right)-{\zeta }_j{\tilde{S}}_{z_j<z_a}^p\left(z_j\right)\right]d_2\end{array}\end{array}.$$

(19)

assuming the MD source is located at the m-th layer, and its vertical position is z_a, as illustrated in Figure 6. Applying Equation (18) to all boundaries, the amplitude of upgoing/downgoing waves at any boundary can be obtained recursively. For example, if we start from z₁ and recursively backward to z_n-1, the x_2n-1 can be obtained directly,

$$x_{2n-1}={\left(-1\right)}^{2n}{\displaystyle \prod _{i=2m}^{2n-1}{R}_i{\left(Q_i-R_i{T^{\prime }}_{i-1}\right)}^{-1}}\left[\left(Q_{2a}+R_{2a}{T^{\prime }}_{2a-1}\right){\psi }_1-R_{2m-1}{\psi }_2\right]{{T^{\prime }}_{2a}}^{-1}{\left(1-{T^{\prime }}_{2n-1}{Q}_{2n}^{-1}{R}_{2n}\right)}^{-1}$$

(20)

where

$$\begin{array}{l}{\psi }_1=c_{21,m-1}\left[{\zeta }_a{\tilde{S}}_{z_{m-1}\ge z_a}^p\left(z_{m-1}\right)-{\zeta }_{m-1}{\tilde{S}}_{z_{m-1}<z_a}^p\left(z_{m-1}\right)\right]s_1-\\ c_{11,m-1}\left[{\zeta }_a{\tilde{S}}_{z_{m-1}\ge z_a}^p\left(z_{m-1}\right)-{\zeta }_{m-1}{\tilde{S}}_{z_{m-1}<z_a}^p\left(z_{m-1}\right)\right]s_2\end{array}$$

(21)

$${\psi }_2=c_{21,m}{\zeta }_m{\tilde{S}}^p\left(z_a\right)s_1+c_{11,m}{\zeta }_m{\tilde{S}}^p\left(z_a\right)s_2$$

(22)

$${T^{\prime }}_j=\{\begin{array}{c}{Q}_j^{-1}{T}_j,j=1\\ {\left(Q_j-R_j{T^{\prime }}_{j-1}\right)}^{-1}{T}_j,j=2,\dots ,2n-3\end{array}$$

(23)

where the superscript −1 represents the inverse operator of the matrix. It should be highlight that the computational stability of Equations (20)–(23) primarily depends on $Q_j^{-1}$. Since $Q_j^{-1}$ is composed of the eigenvectors of the system matrix U, this term always exists, ensuring that there is no numerical overflow problem during the calculation of the MAPM algorithm.

In order to verify the accuracy of the MAPM algorithm, a three-layer model is depicted in Figure 7, where the upper and lower surrounding rocks contain a TI medium, the middle layer features a BA medium, and the specific conductivity of each layer is shown in the figure. Figure 8b–f represents the components’ responses of TIL for different dip angles in the three-layer model. The solid lines are the computation results of the MAPM algorithm, and the scattered points represent the simulation results of the finite element method (FEM) with the operation frequency of 10 kHz and a spacing of 0.5 m. Obviously, the calculation results of the two methods coincide in any dipping wells, which verifies the accuracy of the MAPM method.

Although the recursive nature of the MAPM method is straightforward, it necessitates a significant number of operations on 2 × 2 matrices and does not circumvent the need to solve the double IFT, which leads to the large calculation time of the MAPM method in the case of a large number of layers.

4.2. Approximating the Scattered Fields with Layered TI Model

An approximation method with the layered TI model is built to further accelerate the solving of the scattered fields in the layered BA formation. Hence, the calculation for the scattered field of the layered TI model is well developed, enabling the simplification of the dual infinite integral into a one-dimensional infinite integral, commonly referred to as the Sommerfeld integral, which effectively reduces the computational complexity. In addition, as mentioned in Section 3.3, the spectral scattered fields in layered media also satisfy the symmetry and anti-symmetry. Therefore, the infinite integral can be further simplified to a semi-infinite integral,

$$S_{TI}^2\left(r\right)={\displaystyle \frac{1}{4{\pi }^2}}{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }{\tilde{S}}_{TI}^2\cdot }{e}^{i\left(l{k}_l+q{k}_q\right)}d{k}_{l}d}{k}_q={\displaystyle \frac{1}{4{\pi }^2}}{\displaystyle {\int }_0^{\infty }{\tilde{S}}_{TI}^2{J}_n(k_{\rho }\rho )}d{k}_{\rho },$$

(24)

where J_n(·) is the n-th Bessel function, and n = 0, 1. $S_{TI}^2$ and ${\tilde{S}}_{TI}^2$ represent the spatial and spectral scattered field in the TI model, respectively. After the above simplification, the computational speed of the scattered field can be increased by one order of magnitude.

In order to verify the feasibility of the approximation, the magnetic dipole source and the three receiver points P₁, P₂ and P₃ are placed at the positions shown in Figure 7. The convergence of the spectral primary field and the spectral scattered field of coaxial and coplanar components at point P₁ are represented by the black and blue solid lines in Figure 9a–c, respectively. It can be noticed that the spectral scattered field converges rapidly and contributes a smaller percentage, which allows the spectral scattered fields to be approximated without compromising the calculation accuracy of the total field. The blue scatter plots in Figure 9a–c show the spectral scattered field convergence at P₁ in the equivalent TI model for the coaxial and coplanar components by the equivalent formulae proposed in Section 3.3. It can be observed that the convergence process of the spectral scattered fields in the BA model overlaps with the equivalent TI model, confirming the feasibility of the approximation method. In addition, the approximate effect of the spectral scattered fields at P₂ and P₃ is also given in Figure 9a–c, which further verifies the applicability of the approximation method.

4.3. Best Approximation Model Analysis

Generally, azimuthal angle effects are neglected in layered media, resulting in only five magnetic components (H_xx, H_yy, H_zz, H_xz and H_zx) exhibiting non-zero responses. The optimal construction of the equivalent TI medium of coaxial and coplanar components is given by analyzing the convergence of spectral primary fields in the previous sections. In order to verify whether these TI media can still obtain the best results when approximating the scattered fields, an error analysis is performed using the model shown in Figure 7. Figure 10 shows the average errors of the scattered fields of coaxial and coplanar components calculated by approximation and the MAPM algorithm in the model shown in Figure 7 under different dipping wells. Two BA cases, σ_x > σ_y ([σ_x, σ_y, σ_z] = [0.25, 0.1111, 0.0667] S/m) and σ_x < σ_y ([σ_x, σ_y, σ_z] = [0.1111, 0.25, 0.0667] S/m) are used in this example. It can be observed that the best approximation effect of H_xx can be obtained by σ_h = σ_y, which is consistent with the analysis results in Section 3.3. The same conclusion can be reached through the analysis of H_yy and H_zz. This illustrates that the optimal construction of the spectral primary field is also applicable to the approximation of the spectral scattered fields.

Analogous to the coaxial and coplanar components, approximate formulae are available to calculate the scattered fields for the cross magnetic components (H_xz and H_zx). Figure 11 shows the average errors of the cross components scattered fields under the same simulation conditions as in Figure 10. For the case of σ_x > σ_y, the most accurate approximation of the scattered field for the H_xz and H_zx components can be obtained by using the σ_h = σ_y for any dipping wells. In contrast, the σ_h = sqrt (σ_x* σ_y) method can obtain the optimal equivalent TI medium model for the H_xz and H_zx components for the case of σ_x < σ_y. In summary,

Table 2. Optimal construction of the equivalent TI model for different magnetic components.

	Components	${\mathit{H}}_{\mathit{x}\mathit{x}}$	${\mathit{H}}_{\mathit{y}\mathit{y}}$	${\mathit{H}}_{\mathit{z}\mathit{z}}$	${\mathit{H}}_{\mathit{x}\mathit{z}},{\mathit{H}}_{\mathit{z}\mathit{x}}$
Cases
σx > σy		σh = σy	σh = σx	${\sigma }_h=\sqrt{{\sigma }_x*{\sigma }_y}$	σh = σy
σx < σy					${\sigma }_h=\sqrt{{\sigma }_x^2+{\sigma }_y^2}$

shows the optimal construction of the equivalent TI model used to calculate the scattered fields for coaxial, coplanar and cross magnetic components.

5. Numerical Examples

5.1. Accuracy Analysis

In order to verify the accuracy of the approximation method, the TIL response of the model in Figure 7 was calculated. Assuming the operating frequency of the TIL tool is 100 kHz and the source spacing is 0.8 m, Figure 12a–e shows the imaginary part of the magnetic components calculated by utilizing the equivalent TI medium and BA medium in different dipping angles when β = 0°. In Figure 12, the solid lines represent the calculated results of the magnetic components of BA medium, and the dotted lines represent the magnetic components results obtained after calculating the scattered field using the equivalent TI medium. Obviously, the results of the two models coincide with each other, which proves the accuracy of the approximation method. Figure 13a–e further illustrates the absolute errors calculated in Figure 12. It can be observed that the maximum absolute error of each magnetic component approximation is less than 0.5%, and the error mainly concentrated at the interface.

To further verify the applicability of the approximation method, Figure 14 shows the simulation results under different azimuth angles when the α = 80°. It can be observed that the approximation is still effective.

Table 3. Calculation speed comparison and average errors of each component response between the MAPM and approximation algorithms in different dipping and azimuth angles.

α, β (°)	Average Errors (%)					Relative Time Cost (MAPM (s))	Relative Time Cost (Approximation(s))	Speed-Up Factor
	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{x}\mathit{x}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{y}\mathit{y}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{z}\mathit{z}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{x}\mathit{z}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{z}\mathit{x}}$
0, 0	0.203	0.212	0.198	0.221	0.224	0.0321	0.0008	41.25
45, 0	0.234	0.231	0.222	0.262	0.260	0.0382	0.0010	37.34
90, 0	0.257	0.248	0.242	0.277	0.269	0.0739	0.0022	33.56
80, 30	0.264	0.252	0.245	0.279	0.271	0.0858	0.0025	33.89
80, 45	0.268	0.257	0.247	0.280	0.272	0.0937	0.0029	32.51
80, 75	0.276	0.267	0.253	0.284	0.275	0.1309	0.0044	30.05

and

Table 4. The minimum errors and maximum errors of each component response between the MAPM and approximation algorithms in different dipping and azimuth angles.

α, β (°)	Minimum Errors (%)					Maximum Errors (%)
	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{x}\mathit{x}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{y}\mathit{y}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{z}\mathit{z}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{x}\mathit{z}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{z}\mathit{x}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{x}\mathit{x}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{y}\mathit{y}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{z}\mathit{z}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{x}\mathit{z}}$	${\mathit{\sigma }}_{\mathit{a}}^{\mathit{z}\mathit{x}}$
0, 0	0.022	0.023	0.022	0.024	0.024	0.448	0.468	0.437	0.488	0.495
45, 0	0.025	0.025	0.024	0.028	0.028	0.468	0.462	0.444	0.524	0.52
90, 0	0.028	0.027	0.026	0.030	0.029	0.514	0.496	0.484	0.554	0.538
80, 30	0.028	0.028	0.027	0.030	0.029	0.528	0.504	0.49	0.558	0.542
80, 45	0.029	0.028	0.027	0.030	0.029	0.536	0.514	0.494	0.56	0.544
80, 75	0.030	0.029	0.0275	0.030	0.029	0.552	0.534	0.506	0.568	0.551

show the comparison of simulation speed and calculation accuracy of each component response between the MAPM and approximation algorithms in different dipping and azimuth angles. In this example, the specific parameters of the computing platform used are as follows, the operating system is Windows 10 64 bit Professional 21H2, the CPU is an Intel (R) Core (TM) i9-9900, and the memory capacity is 16 GB. The simulation speed of scattered fields can be increased by more than 30 times utilizing the approximation.

5.2. Applicability Analysis

To verify the applicability and stability of the new forward modeling method, the Oklahoma model with BA was established, as shown in Figure 15. The thickness of the model ranges from 0.6 m to 5.5 m, including high and low resistance layers, and has different cases of BA. Figure 16 shows the TIL responses of coaxial, coplanar and cross magnetic components in the Oklahoma model with different frequencies using the MAPM method and approximation method. In this example, the dipping angle is 85° and the spacing is 1.0 m. It can be found that the approximation is applicable to arbitrary cases of formation and tool parameter conditions. The following conclusions can be drawn by analyzing the following: (I) when the frequency increases, each component signal becomes more sensitive to conductivity changes; (II) when σ_x > σ_y, each component curve exhibits the greatest separation, indicating a stronger influence of formation conductivity in the x-direction; (III) within the HZ well, the responses of cross components are zero in thick layers, and the influence of anisotropy is mainly reflected at the boundaries, while the responses of ${\sigma }_a^{xz}$ and ${\sigma }_a^{zx}$ have opposite patterns; and (IV) variations in horizontal anisotropy result in significant differences in coplanar component responses, while coaxial components are minimally affected. Overall, the different component responses show different sensitivities to different cases of BA, providing an opportunity to accurately characterize formation anisotropy.

6. Discussion

This paper introduces a fast and novel forward modeling algorithm for TIL in layered BA formations. The proposed method employs a two-level singularity subtraction technique for the rapid computation of spectral primary fields. Additionally, a concise recursive method (MAPM) and a fast approximation method are developed for the calculation of spectral scattered fields. This method significantly improves the forward modeling speed of TIL responses in layered BA media. However, it also has some limitations, as outlined below:

6.1. Advantages

• The new approximation method is at least two orders of magnitude faster than the current modeling approach while maintaining a computational precision error within 0.5% after processing the primary and scattered fields, making real-time data processing possible.
• The new approximation method is highly applicable to a wide range of formation conditions (including dip angle, thickness, and anisotropy) and tool parameters (such as frequency and spacing).
• In addition to the approximation method, this paper also presents a concise and stable MAPM algorithm for the rigorous solution of scattered fields in layered BA media. This algorithm can be extended to more complex anisotropic models, such as cross anisotropy, full anisotropy.

6.2. Limits

• The approximation method is currently applicable only to layered BA formations and is not yet suitable for more complex anisotropic cases, such as cross anisotropy and full anisotropy. Additionally, the method is designed for actual formations where the coefficient of horizontal anisotropy is small due to formation development and depositional features. As a result, it may not be applicable to anisotropic models in other fields.
• The coil system of commercial TIL tools is typically designed with a single transmitter and two receivers. Due to the varying distances between the transmitting coil and the two receiving coils, it is necessary to separately approximate the formations at the locations of the two coils during the measurement process, making the procedure quite complex.

7. Conclusions

A fast and novel forward modeling algorithm of TIL in layered BA formation is presented in this paper. The two-level singularity subtraction on the primary fields has greatly narrowed the integral range and reduced the computational complexity. The first-level subtraction scheme works well in HA/HZ wells, whereas it is inapplicable to low deviated wells. By contrast, the second-level subtraction scheme can not only accommodate wells with arbitrary dipping angles, but also improves the convergence of spectral fields by four orders of magnitude.

The MAPM algorithm is concise and stable and can accurately calculate the spectral scattered fields for a layered BA medium with any thickness and arbitrary dipping wells. The optimal construction of the equivalent TI model is defined to approximate the scattered fields of coaxial, coplanar and cross components. The optimal construction formulae corresponding to H_xx, H_yy and H_zz are σ_h = σ_y, σ_h = σ_x and σ_h = sqrt (σ_x* σ_y) for the any BA cases, respectively. For the case of σ_x > σ_y, H_xz and H_zx can be most accurately approximated by σ_h = σ_y. Otherwise, their optimal construction formula is σ_h = sqrt (${\sigma }_x^2$ + ${\sigma }_y^2$).

Numerical examples demonstrate that the new simulation method is at least two orders of magnitude faster than the current modeling approach while maintaining computational precision error within 0.5% after processing the primary and scattered fields. Meanwhile, the new method has good performance in any formation conditions and tool parameters, which can meet the requirements of real time data processing. In the future, we will try to extend the fast computation method to the full anisotropic formations to realize the fast forward and inversion of TIL in the full anisotropic formations.