Logging while drilling electromagnetic wave resistivity data can be used for fluid identification and fluid saturation calculation in reservoir evaluation. In order to realize fast inversion to obtain true formation parameters in real-time, it is usually necessary to simplify complex borehole and formation environment (three-dimensional problems) into a series of one-dimensional (1D) formation models^[1,2,3,4,5]. This method is feasible in logging while drilling (LWD) because of two main reasons. First, in the LWD process, the influence of borehole and drilling fluid invasion is small and can be neglected. Second, the complex wellbore trajectory can be constrained by progressive window opening and simplified to a series of simple-shaped wellbore trajectories (i.e. a series of 1D formation models). The traditional 1D fast LWD electromagnetic wave resistivity numerical simulation is mainly based on the transverse isotropic formation model with vertical symmetry axis (VTI). Under low-energy sedimentary environment, such as deep-sea and deep-lake deposits, the formation can present standard VTI medium characteristics (Fig. 1a). VTI medium is a special case of uniaxial anisotropic medium, which has the same resistivity R_h in any horizontal direction, different resistivity R_v in vertical direction, and in general, R_h≤R_v. The calculation of electromagnetic field in this kind of medium is more mature. By introducing the Hertz potential function, the wave equation of electromagnetic field can be simplified into two uncoupled scalar wave equations. The final expression of electromagnetic field is a single integral containing 0-order and 1-order Bessel functions^[5,6,7].

With the deepening of petroleum exploration, complicated oil and gas reservoirs have become the focus of research, which limits the application of VTI stratigraphic model. Not all strata can be simplified to VTI model. The strata formed in high energy sedimentary environment (river and desert environment) usually contain inclined beddings (Fig. 1b). In this circumstance, the strata are transversely isotropic (TTI) strata with tilted symmetry axis, which have the same resistivity along the bedding direction but different resistivity in the direction perpendicular to the bedding. If it is still simplified as VTI stratigraphic model, the parameters such as resistivity obtained by inversion is no longer reliable. If the resistivity is used for qualitative evaluation and quantitative calculation, the interpretation accuracy will drop. Clear understanding of electromagnetic wave logging response characteristics in TTI medium can make the qualitative and quantitative logging interpretation more comprehensive and effective.

The theory and calculation of electromagnetic field based on TTI medium are very complex, and there is no relevant literature published in China. Anderson et al.^[8-9] from Schlumberger briefly discussed the response characteristics of induction logging to formation anisotropy caused by inclined bedding. Wang et al.^[10,11,12] gave the forward and inverse numerical simulation algorithm of multi-component induction logging to formation anisotropy caused by inclined bedding and analyzed response characteristics of multi-component induction logging under different inclined bedding conditions. Anderson and Wang et al.^{[8,9,10,11,12]} all pointed out that the influence of inclined bedding on log response characteristics should not be neglected. Especially in the application of actual data, suitable formation model should be established by considering the actual formation conditions. In this study, the authors used anisotropic angle and azimuth to characterize the anisotropic properties of TTI medium with inclined bedding, and deduced the solution of electromagnetic field in TTI medium based on dyadic Green's function. Then, taking half-space medium as an example, the reflection and transmission characteristics of electromagnetic wave in VTI medium and TTI medium were analyzed by numerical simulation method. The limitation of traditional Hertz potential and generalized reflection coefficient methods were pointed out. The reliability of the algorithm was verified, and the LWD electromagnetic wave responses under different borehole deviations, anisotropic angles and anisotropic azimuths were figured out.

In an infinitely thick anisotropic medium, the wave equation of the electric field can be expressed as:

In VTI medium, the conductivity tensor of anisotropic formation can be expressed as:

In TTI medium, anisotropic angle Ψ and anisotropic azimuth χ can be defined to characterize its anisotropic property (Fig. 2), then the conductivity tensor ${{\mathsf{\vec{\sigma }}}_{\mathrm{b}}}$ in TTI medium can be expressed as:

${{\mathrm{\vec{R}}}_{\mathrm{ }\!\!\psi\!\!\text{ }}}=\left( \begin{matrix} \cos \psi & 0 & \sin \psi \\ 0 & 1 & 0 \\ -\sin \psi & 0 & \cos \psi \\\end{matrix} \right)$, ${{\mathrm{\vec{R}}}_{\mathrm{ }\!\!\chi\!\!\text{ }}}=\left( \begin{matrix} \cos \chi & \sin \chi & 0 \\ -\sin \chi & \cos \chi & 0 \\ 0 & 0 & 1 \\\end{matrix} \right)$

Using Fourier transform, the expression of electric field in space domain is:

Equation (4) can be further rewritten as:

where

$\mathrm{W}=\mathrm{k}\cdot \mathrm{k}+i\omega \mu {{\mathsf{\vec{\sigma }}}_{\mathrm{c}}}$

$\mathrm{k}=\left( \begin{matrix} 0 & -{{k}_{z}} & {{k}_{y}} \\ {{k}_{z}} & 0 & -{{k}_{x}} \\ -{{k}_{y}} & {{k}_{x}} & 0 \\ \end{matrix} \right)$

${{\mathrm{W}}^{-1}}=\frac{adj\left( \mathrm{W} \right)}{\left| \mathrm{W} \right|}$

The determinant |W| is a fourth-order polynomial of k_z as follows:

where a, b, c, d, e are expressions of k_x and k_y. It can be further expressed as follows:

When |W|=0, there exist four solutions for k_z, representing two types of waves. $k_{z,\text{I}}^{+}$ and $k_{z,\text{II}}^{+}$ are the upward z-direction wavenumber of type I wave and type II wave. $k_{z,\text{I}}^{-}$ and $k_{z,\text{II}}^{-}$ are the downward z-direction wavenumber of type I wave and type II wave^[13,14].

Substituting equation (7) into equation (5):

where R=r-r°, dk=dk_xdk_ydk_z.

By using residue theorem, when the position z of the receiving coil is above the source position z', the electric field E(R) is:

where

${{\mathrm{r}}_{\mathrm{s}}}=\left( x-{x}' \right)\hat{x}+\left( y-{y}' \right)\hat{y}$

${{\mathrm{k}}_{\mathrm{s}}}={{k}_{x}}\hat{x}+{{k}_{y}}\hat{y}$

$\text{d}{{\mathrm{k}}_{\mathrm{s}}}=\text{d}{{k}_{x}}\text{d}{{k}_{y}}$

$\mathrm{W}\left( k_{z,\text{I}}^{+} \right)={{\left. \frac{adj\left[ \mathrm{W}\left( k_{z,\text{I}}^{+} \right) \right]}{a\left( {{k}_{z}}-k_{z,\text{I}}^{-} \right)\left( {{k}_{z}}-k_{z,\text{II}}^{+} \right)\left( {{k}_{z}}-k_{z,\text{II}}^{-} \right)} \right|}_{{{k}_{z}}=k_{z,\text{I}}^{+}}}$

$\mathrm{W}\left( k_{z,\text{II}}^{+} \right)={{\left. \frac{adj\left[ \mathrm{W}\left( k_{z,\text{II}}^{+} \right) \right]}{a\left( {{k}_{z}}-k_{z,\text{I}}^{-} \right)\left( {{k}_{z}}-k_{z,\text{I}}^{+} \right)\left( {{k}_{z}}-k_{z,\text{II}}^{-} \right)} \right|}_{{{k}_{z}}=k_{z,\text{II}}^{+}}}$

When the receiving coil position z is below the source position z', the electric field E(R) is:

where

$\mathrm{W}\left( k_{z,\text{I}}^{-} \right)={{\left. \frac{adj\left[ \mathrm{W}\left( k_{z,\text{I}}^{-} \right) \right]}{a\left( {{k}_{z}}-k_{z,\text{I}}^{+} \right)\left( {{k}_{z}}-k_{z,\text{II}}^{+} \right)\left( {{k}_{z}}-k_{z,\text{II}}^{-} \right)} \right|}_{{{k}_{z}}=k_{z,\text{I}}^{-}}}$

$\mathrm{W}\left( k_{z,\text{II}}^{-} \right)={{\left. \frac{adj\left[ \mathrm{W}\left( k_{z,II}^{-} \right) \right]}{a\left( {{k}_{z}}-k_{z,\text{I}}^{-} \right)\left( {{k}_{z}}-k_{z,\text{I}}^{+} \right)\left( {{k}_{z}}-k_{z,\text{II}}^{+} \right)} \right|}_{{{k}_{z}}=k_{z,\text{II}}^{-}}}$

Suppose the medium has N layers, the source is in S_th layer, the receiving coil is in the L_th layer. The electric field in the position of the receiving coil is:

δ_LS is Kronecker delta function, when L=S, δ_LS=1; when L≠S, δ_LS=0; E_d(R) is the direct electrical filed; E_s(R) is the scatter electrical field (including reflected and transmitted field). Hence, as long as the direct field and scattering field are derived, the electric field in each layer can be obtained. The corresponding magnetic field component can be obtained from the electric field component by the following expression:

The deduction of electromagnetic field in multi-layer TTI medium can be referred according to Wang et al.^[10,11,12]. Due to the limited space, we will not discuss in detail here.

Based on the deduction in the previous section, the electromagnetic field at the receiving coil in TTI medium can be obtained, and then the voltage of the receiving coil can be calculated. Taking the structure of traditional LWD electromagnetic wave with single-emission and double-receiving coils as example, the voltages at the two receiving coils are V₁ and V₂, the amplitude ratio EATT and phase difference ∆ϕ are defined as follows:

The phase difference and amplitude ratio can be converted into amplitude ratio resistivity R_ad and phase difference resistivity R_ps by using the calibration plate of phase difference and amplitude ratio.

Firstly, the calibration plate of phase difference and amplitude ratio resistivity is used to verify the validity of this algorithm. Fig. 3 shows the calibration plate at traditional operation frequency of 2 MHz and spacing from transmitting coil to recording point of 40.64 cm (16 inches). It can be seen that the computation results based on dyadic Green's function perfectly match with the results using Hertz's potential function method. Secondly, two five-layer formation models are used to validate the proposed algorithm (with well deviation of 0° and 30° respectively with each layer thickness of 2.0 m, corresponding resistivity parameters are in Table 1). It can be seen from Fig. 4a and Fig. 4b that the proposed algorithm is also applicable to multi-layer media.

In order to analyze the reflection and transmission characteristics of electromagnetic wave at interfaces of different medium, the reflection and transmission coefficients in half-space are given. Fig. 5 shows the half-space models under different conditions, and Table 2 shows the corresponding conductivity parameters. σ_1u and σ_2u represent the upper layer conductivities, σ_1d and σ_2d represent the lower layer conductivities. If the upper layer is isotropic, then σ_1u=σ_2u; if it is anisotropic, then σ_1u≠σ_2u. It is the same for the lower layer.

In isotropic and VTI formations, Ψ=0, χ=0. In TTI formation, Ψ≠0. Due to limited space, only the case when χ=0 is presented here (similar with the cases of reflection and transmission when χ≠0, no more details). Fig. 6 shows the variation of Fresnel reflection and transmission coefficients with incident angle at the interface of different models. R_i,j and T_i,j represent Fresnel reflection and transmission coefficients, respectively. i and j denote the types of electromagnetic waves (I or II type).

From Fig. 6a to Fig. 6c, it can be seen that when the media on both sides of the interface are isotropic or VTI media, i.e. Ψ =0 and χ =0, R_Ⅰ,Ⅱ=R_Ⅱ,Ⅰ=T_Ⅰ,Ⅱ=T_Ⅱ,Ⅰ=0. This means only R_Ⅰ,Ⅰ, R_Ⅱ,Ⅱ, T_Ⅰ,Ⅰ and T_Ⅱ,Ⅱ exist and there are no cross-reflection and transmission waves. That is to say, when the incident wave is type I, the reflected and transmitted I waves will be produced; when the incident wave is type II, the situation is similar. Therefore, when both sides of the interface are isotropic or VTI media, the reflection and transmission of type I and type II waves wouldn’t couple at the interface. The electromagnetic fields of type I and type II waves (also known as TE and TM waves, which are transverse electric wave and transverse magnetic wave) can be calculated separately. This is the theoretical basis of traditional calculation of TE and TM electromagnetic fields based on VTI media.

From Fig. 6d to 6f, it can be seen that when TTI stratum exists on one side of the interface, the cross components (R_Ⅰ,Ⅱ, R_Ⅱ,Ⅰ, T_Ⅰ,Ⅱ, T_Ⅱ,Ⅰ) are no longer zero. This means that after reflection and transmission, the type I wave and type II wave couple, and the electromagnetic field is the coupling field of two types of waves. In this case, the traditional electromagnetic field calculation method based on VTI medium is no longer applicable. The method of dyadic Green's function proposed in this paper can be used.

Assuming that the anisotropic formation is infinitely thick, the resistivity in one direction is R₁=2.0 Ω•m, and the resistivity in the other direction is R₂=20.0 Ω•m. The operation frequency is 2.0 MHz and the spacing from the transmitting coil to recording point is 40.64 cm (16 inches). We use phase difference to analyze the influence of the anisotropic angle and anisotropic azimuth. Fig. 7 shows the variation of traditional LWD electromagnetic wave phase difference with well deviation angle θ. When ψ=0, the phase difference decreases monotonously with the increase of θ, while the anisotropic azimuth has no effect on phase difference with different θ. When ψ ≠0, if χ=0, there exists a critical well deviation angle θ_c=90-ψ, that is, when ψ is 30°, 45° and 60°, θ_c is 60°, 45°and 30° respectively (Fig. 7a). When θ≤θ_c, the phase difference decreases with the increase of the deviation angle; when θ>θ_c, the phase difference increases with the increase of the deviation angle. If χ≠0, there still has a critical well deviation angle θ_c with the same change rule as that shown in Fig. 7a (Fig. 7b), but θ_c≠90-ψ.

Fig. 8 shows the variation of phase difference of traditional LWD electromagnetic wave with anisotropy angle ψ. Firstly, when θ=0, the phase difference decreases monotonously with the increase of ψ, and χ has no effect on the variation of phase difference. When θ≠0, if χ=0, there exists a critical anisotropic angle ψ_c=90-θ, that is, when θ is 30°, 45° and 60°, ψ_c is 60°, 45° and 30° respectively (Fig. 8a). When ψ ≤ ψ_c, the phase difference decreases with the increase of anisotropy angle. When ψ>ψ_c, the phase difference increases with the increase of anisotropic angle. If χ≠0, there still has a critical anisotropic angle ψ_c (Fig. 8b), the variation rule is the same as that shown in Fig. 8a, but ψ_c≠90-θ.

It is worth to note that the ‘horn’ produced by LWD electromagnetic responses at the interface is mainly caused by factors such as well deviation and resistivity contrast (the electric field on both sides of the interface is discontinuous^[15]), but in fact the anisotropic angle and azimuth are also important factors affecting the ‘horn’. Fig. 9a is a two-layer formation model and its corresponding resistivity parameters. The upper layer is isotropic (χ₀ =0, ψ₀ =0), and the lower layer is uniaxial anisotropic.

Firstly, the anisotropic azimuth χ₁ of the lower layer is set at 0 to investigate the traditional LWD electromagnetic responses at the interface at different ψ₁ and θ. When the lower layer is VTI medium, i.e. ψ₁=0 and χ₁=0 (Fig. 9b and 9c), it can be seen that with the increase of well deviation, both the amplitude ratio resistivity R_ad and phase difference resistivity R_ps of the lower formation increase. When the deviation angle is less than 60°, there will be no horns at the interface. If the lower formation is TTI medium, and ψ₁=30° (Fig. 10a and Fig.10b), when θ is 30°, 45° and 60° respectively, compared with Fig. 9b and Fig. 9c, the amplitude of R_ad and R_ps in the lower layer change obviously. In addition, R_ad produces obvious horn at the interface. When θ=60°, R_ps also produces obvious horn at the interface. Therefore, compared with VTI formation, under the same deviation angle, the existence of anisotropic angle may cause horns at the interface.

Based on the two-layer formation model and parameters shown in Fig. 9a, the LWD electromagnetic responses at the interface are investigated when the lower formation has different anisotropic azimuths and χ₁≠0°. It can be seen from Fig. 11, the horns of R_ad and R_ps at the interface vary with the anisotropic azimuth, so the anisotropic azimuth is one of the factors causing horns.

Fig. 12 shows two five-layer models with well deviation of 45°. The specific formation parameters are shown in Table 3. There is no inclined bedding in each layer of Fig. 12a, while layer 2 and layer 4 in Fig. 12b have inclined beddings. The anisotropic angles and azimuths of layer 2 and layer 4 in Fig. 12b are different. Fig. 13 shows the traditional LWD electromagnetic responses at the transmission frequency of 2.0 MHz and the spacing from transmission coil to the recording point of 40.64 cm (16 inches). It can be seen from Fig. 13, although the resistivity parameters of the two formation models are the same, their anisotropic characteristics are different, so their amplitude of LWD electromagnetic resistivity responses and response characteristics at the interface are different. Clearly, if the formation is always simplified to VTI formation without considering the change of its internal bedding, neither the LWD response nor the formation real resistivity from inverse modeling can characterize the true formation electrical characteristics.

In order to satisfy the demand of fast resistivity inversion, traditional LWD electromagnetic wave forward modeling is mainly based on VTI formation model, but the actual formation is more complicated. VTI formation is only a special case of anisotropic formation. Based on dyadic Green's function, a set of forward calculation methods for LWD electromagnetic responses in multi-layer TTI formation are derived in detail. The algorithm is validated by numerical examples. The results show that the proposed algorithm is not only applicable to VTI formations, but also to TTI formations, and has strong universality.

Through numerical simulation of reflection and transmission characteristics of electromagnetic wave in half space, it is concluded that when both sides of the interface are isotropic or VTI formation, the electromagnetic fields generated by two types of electromagnetic waves (type I and type II) wouldn’t couple, and the forward simulation of LWD can adopt the traditional Hertzian potential function method. When one side of the interface is TTI formation, the electromagnetic fields generated by type-I wave and type-II wave will couple, and the forward numerical simulation can adopt the algorithm presented in this paper.

The effects of anisotropic angle and anisotropic azimuth on traditional LWD electromagnetic responses can’t be ignored. In TTI formation, there is a critical well deviation angle θ_c and critical anisotropy angle ψ_c. When θ≤θ_c, the phase difference decreases with the increase of well deviation angle; when θ>θ_c, the phase difference increases with the increase of well deviation angle; when ψ≤ψ_c, the phase difference decreases with the increase of anisotropy angle; when ψ>ψ_c, the phase difference increases with the increase of anisotropy angle. In addition, when χ=0, the θ_c = 90-ψ and ψ_c =90-θ.

The horn produced by LWD electromagnetic responses at the interface is not only related to well deviation and resistivity contrast etc., but also affected by anisotropic angle and azimuth. For example, in strata with inclined beddings, horn could occur even at low well deviation.

E—electrical field, V/m;

E_d(R)—direct field, V;

E_s(R)—scatter field, V;

$\mathrm{\tilde{E}}\left( {{k}_{x}},{{k}_{y}},{{k}_{z}} \right)$—wave number domain of electrical field;

EATT—amplitude ratio, dB;

f—frequency, Hz;

H—magnetic field, A/m;

i, j—electromagnetic wave types (type I and type Ⅱ);

k—wave number, dimensionless;

k_x—wave number in x direction, rad/m;

k_y—wave number in y direction, rad/m;

k_z—wave number in z direction, rad/m;

$k_{z,\text{I}}^{+}$, $k_{z,\text{II}}^{+}$ — wave number of type I and typeⅡ in z upward direction, rad/m;

$k_{z,\text{I}}^{-}$, $k_{z,\text{II}}^{-}$—wave number of type I and typeⅡ in z downward direction, rad/m;

M—magnetic source, Am²;

M₀—magnetic source strength, Am²;

r—receiver position of (x, y, z), m;

r°—source position of (x°, y°, z°), m;

R_ad—amplitude ratio resistivity, Ω·m;

R_h—horizontal resistivity, Ω·m;

R_i,j—Fresnel reflection coefficient;

R_ps—phase difference resistivity, Ω·m;

R_v—vertical resistivity, Ω·m;

${{\mathrm{\vec{R}}}_{\mathrm{ }\!\!\psi\!\!\text{ }}}$—rotation matrix of anisotropic dip, (°);

${{\mathrm{\vec{R}}}_{\mathrm{ }\!\!\chi\!\!\text{ }}}$—rotation matrix of anisotropic azimuth, (°);

T_i,j—Fresnel transmission coefficient;

V₁ and V₂—induced voltages of two receiving coils, V;

|V₁| and |V₂|—amplitudes of V₁ and V₂, V;

μ—permeability in free space, H/m;

ω—angular frequency, rad/s;

σ_h—horizontal conductivity, S/m;

σ_ν—vertical conductivity, S/m;

Ψ—anisotropic dip, (°);

χ—anisotropic azimuth, (°);

${{\mathrm{\vec{ }\!\!\sigma\!\!\text{ }}}_{\mathsf{c}}}$—conductivity tensor in VTI medium, S/m;

${{\mathrm{\vec{ }\!\!\sigma\!\!\text{ }}}_{\mathsf{b}}}$—conductivity tensor in TTI medium, S/m;

δ_LS—Kronecker delta function, if L=S, δ_LS=1, if L≠S, δ_LS=0;

∆ϕ—phase difference, (°);

Φ₁ and Φ₂—phase angles of V₁ and V₂, (°);

θ—well deviation angle, (°);

θ_c—critical well deviation angle, (°);

ψ_c—critical anisotropic dip, (°);

σ_1u and σ_2u—resistivities in two directions of upper layer in the two-layer stratigraphic model, S/m;

σ_1d and σ_2d—resistivities in two directions of lower layer in the two-layer stratigraphic model, S/m;

—Hamiltonian operator.