Logging-while-drilling (LWD) electromagnetic (EM) measurements, which can provide both amplitude ratio and phase difference resistivities, have been routinely used for geosteering horizontal wells^[1,2]. Unfortunately, complicated responses such as curve separation and electromagnetic “horns” are usually encountered, which makes the quantitative or even qualitative evaluation of these measurements very difficult. Moreover, the limited curves and lack of azimuth identification of LWD EM data further add difficulty to the prediction of nearby formation boundary^[3].

To date, the processing of LWD EM data from horizontal wells mainly depends on the one-dimensional (1D) planar stratified formation model and gradient algorithm^[4,5]. Researchers abroad have tried to correct well deviation effect on EM logging data and make inversion of anisotropic formation resistivities^[5,6]. However, this data processing method requires the boundary information from other logging data, and therefore it only works for post-drilling interpretation. Meanwhile, the gradient algorithm relies heavily on the initial model, which may result in non-unique inversion results. In order to determine the boundaries around the borehole, the interactive interpretation method is commonly used in horizontal well geosteering and formation evaluation. In this method, by manually adjusting the boundary positions and running the forward simulation accordingly, the formation models with optimal fit between measured and simulated data are obtained^[7]. Another important feature of the interactive method is that the logging data of pilot or offset wells and known geological information are introduced to effectively reduce the uncertainty of the interpretation results. Theoretically, the essence of this method is manual optimization by trial and error, which is time consuming and limited in precision. Hence, it is urgent to develop a new method to significantly improve the efficiency of current interpretation method and the prediction accuracy of bed boundaries.

Rapid 1D forward modeling is the kernel of the processing of LWD EM data from horizontal well. In 1D modeling, the main difficulty lies in the integration of the kernel function^{[8,9,10,11,12]}, which is highly oscillating and converging slowly in horizontal wells. Hence, the direct integral method is no longer applicable. Fast Hankel transform (FHT) method is widely adopted to simulate the response of LWD EM logging currently^[13]. However, FHT method still has many problems: (1) the integral singularity is inevitable when the transmitter and receiver are at the same vertical position^[14]; (2) the integral precision is difficult to control; and (3) the fields at different recording points require repeated derivations. Therefore, a new and efficient integral method needs to be developed to avoid repeated computation of array coil responses and to improve the speed and accuracy of forward modeling of LWD EM measurements.

In this study, first, the forward model was established for the simulation of LWD EM in horizontal well. An optimal boundary match technique in combination with an efficient and direct Sommerfeld integral method were then proposed to accelerate the simulation of LWD EM measurements. Based on the 1D simulator, an interactive inversion method was developed by combining interactive model adjusting scheme and gradient optimization algorithm to process the LWD EM data in horizontal wells to predict formation boundaries. Finally, the forward and inversion methods were applied to process the field data.

LWD EM logging can provide both amplitude ratio and phase difference apparent resistivities by using the coaxial antenna structure with one transmitter and two receivers. These measurements usually have tool spacing of less than 1.5 m and can detect the formations several meters around the borehole. Generally, the formations around the borehole are regarded as parallel with each other^[15,16,17]. Thus, the 3D forward model should consider the borehole, tool, formations and relative dip (θ) between tool axis and formation, as shown in Fig. 1a. This model can be simplified to 1D when simulating the LWD EM responses. This is out of two considerations. One, the effect of borehole on the logging response is usually small or can be neglected^[4]; the other, the borehole effect can be readily corrected.

In order to validate the feasibility of aforementioned dimensionality reduction of the forward model, we established a two-layer formation model in which the tool crossed the formation at an inclination of 80°. Fig. 2 shows the computed apparent resistivity curves from 3D/1D algorithms of the forward model with and without borehole. In the simulation, the borehole radius was set at 0.101 6 m and three types of drilling fluid, oil-based, fresh water-based and salty water-based with resistivity of 1000 Ω•m, 1.0 Ω•m and 0.1 Ω•m, respectively, were considered. It can be seen from the figure that the LWD EM responses of the well using oil-based drilling fluid are identical to those of the 1D model without borehole, validating the accuracy and feasibility of the dimension reduction method. By contrast, the apparent resistivities of the well using water-based drilling fluid differ widely from the 1D result, and the lower the drilling fluid resistivity, the bigger the difference is. Then the calibration chart for the model with borehole was made to correct the LWD EM response. The corrected resistivities are shown in Fig. 2c and Fig. 2d. It can be seen that the logging responses with borehole after correction almost overlap with the 1D results, which indicates the model dimension reduction method also works for the well with water-based drilling fluid.

In the 1D planar-stratified medium shown in Fig. 1b, the EM field radiated by a point source can be solved by recursive algorithm. For the specific case of LWD EM, the magnetic field H at the receiver can be regarded as the superposition of the fields from the horizontal magnetic dipole (HMD) and the vertical magnetic dipole (VMD), which is:

Theoretically, VMD only excites transverse electric (TE) wave in 1D medium, whereas HMD generates both TE and transverse magnetic (TM) waves. Meanwhile, the two waves from magnetic dipole decouple. Here, we assume the transmitter is located at the origin of the formation coordinate system. When the transmitter and receiver are at m-th and n-th layers, Eq. (1) can be written as follows,

where, $k_{n, \mathrm{~h}, z}=\sqrt{i \omega \mu \sigma_{n, \mathrm{~h}}-k_{\rho}^{2}}, k_{n, \mathrm{v}, Z}=\sqrt{i \omega \mu \sigma_{n, \mathrm{v}}-k_{\rho}^{2}}$.

For simplicity, this paper only gives the expressions of the propagation coefficient of TM waves. Detailed derivation of the propagation coefficient of TE waves is presented in Reference [9]. For the TM wave from HMD, F^TM,h can be expressed as follows,

The general reflection and transmission coefficients at the bed boundary can be written as:

Replacing TM,h of the above equation with TE,h or TE,v and substituting k_n,h,z for λ_nk_n,v,z, the general reflection coefficients of TE waves $\tilde{R}_{n, n+1}^{\mathrm{TE}, \mathrm{h}}$ and $\tilde{R}_{n, n+1}^{\mathrm{TE,v}}$ can be readily obtained, the local reflection coefficients $\tilde{R}_{n, n+1}^{\mathrm{TE}, \mathrm{h}}$ and $\tilde{R}_{n, n+1}^{\mathrm{TE,v}}$ could be expressed as:

To compute the magnetic field at the receiver, the Sommerfeld integral of Eq. (2) should be solved. When the relative dips are 0°, 30°, 60°, 80° and 90°, Fig. 3 shows the integrand as a function of k_ρ in the two-layer model. Here, the tool spacing and operation frequency are set at 0.5 m and 2 MHz, and the transmitter coil is fixed at the bed boundary. Obviously, the integral kernel function oscillates more significantly and converges more slowly as the relative dip increases. When the tool is parallel to the bed boundary, the kernel function doesn’t converge, so accurate computation of the total field is still challenging.

The oscillating property of the total field is caused by the scattered and primary fields, and the primary field has analytical expression. Fig. 3b presents the convergence of the scattered spectral field (h^s). Compared with the total field, the scattered spectral field has weaker oscillation. For the same formation model, the relationship between the integral spectrum of the scattered field and the distance of the tool to formation boundary (D) is given in Fig. 4. Apparently, the kernel function drops quickly and converges fast as the D increases. When the relative dip is less than 60° or the D is bigger than 0.1 m, setting the upper integral bound at 100 can ensure the computational error of the magnetic field at less than 0.1%.

Theoretically, the scattered field can be treated as the superposition of multiple reflection/transmission waves. The contribution of multiple waves decreases quickly as the times of reflection/transmission increase. By using this property, we have come up with a new integral method. The new method includes three steps: (1) subtract the contribution of the first reflection/transmission wave from the spectral scattered field; (2) analytically calculate the integral of the subtracted term; and (3) compute the integral of the kernel function of residual field. The scattered field H^s comes from two parts i.e. the integral of 0-th and 1-st order Bessel function,

Separating the first reflection/transmission wave from the scattered field, $H_{0}^{s}$ becomes

After the field separation, more integral terms arise and the integral of first reflection/transmission field has no analytical solution. The forward problem then becomes how to efficiently compute the first reflection/transmission wave.

When k_ρ approaches to infinity, k_n,z approximates to k_ρ and the general reflection coefficient almost equals the local reflection coefficient. Thus, the local reflection coefficients of TE and TM waves can be expanded as,

Neglecting the high-order term of the above equation, $xi_{1}^{1}$ and $xi_{2}^{1}$ can be approximated as:

In Eq. (13), $\tilde{\xi}_{1}$ and $\tilde{\xi}_{2}$ are independent of the integral invariant. The Sommerfled integral of Eq. (13) can be analytically calculated using the following identity:

where $k_{n}=\sqrt{i \omega \mu \sigma_{n, \mathrm{~h}}}$

Based on Eq. (13) and Eq. (14), Eq. (11) can be transformed as:

In the above equation, the second term is the approximated first reflection/transmission wave. Thus, the issue is converted to integrating the first term of Eq. (15). Similarly, the integral over 1-st order Bessel function can be calculated.

In order to validate the effectiveness of the new integral method, we built a formation model shown in Fig. 5a. Fig. 5b shows the convergence of the scattered field before and after treated with the new integral method at the angles between tool axis and formation of 90° and 89°. Obviously, when the first reflection/transmission wave are approximated and subtracted, the kernel function of remaining scattered field drops quickly with the increase of k_ρ. Meanwhile, the new Sommerfeld integral method can make the convergence of the intergral kernel function increase by more than three orders of magnitude. Hence, the integration interval set to [0, 50] is adequate to guarantee the engineering precision. Table 1 further compares the computational efficiency and accuracy of the new direct integral method and traditional FHT method. When the computational error is 1×10^-4, the FHT method needs at least 260 filtering points, whereas 40 GLQ points are enough for the new method^[18]. Compared with traditional FHT method, the new method has much higher integral accuracy and is six times faster. In addition, the calculation speed of the new method is hardly affected by the number of receivers, and therefore is more applicable to the simulation of array LWD EM responses.

The total layer number of the forward model is an important factor affecting the simulation speed of LWD EM. The layer number is usually adversely proportional to the modeling speed^[15]. Generally, the geosteering model contains a large number of formation layers, and proper cut-off must be applied in the forward modeling of LWD EM. Fixed cut-off distance is usually used due to its easy implementation. For instance, we only consider the boundaries which are within 10 m from the logging tool. Since the detection scope of the LWD EM tool varies with the formation resistivity and operation frequency, this traditional boundary cut-off method may lead to low computational accuracy or long time consumption. Therefore, we have worked out a new cut-off method to adaptively match the bed boundaries.

In homogeneous formation, the attenuation of EM filed is often indicated by skin depth (δ),

In 1D medium, the decaying ratio (s) of EM filed at the bed boundary is a function of skin depth, formation thickness and transmission coefficient. Assuming the source is located at m-th layer, s at the boundary of l-th layer above the tool can be approximately expressed as:

When s is less than 1%, the EM waves almost vanish and the formations further away only have little effect on the tool responses. Hence, s=1% is taken as the threshold for cutting off the bed boundary. To verify the feasibility of the above cut-off method, an anisotropic Oklahoma model shown in Fig. 6a was built^[9]. Fig. 6 shows the computed amplitude ratio resistivity curves and corresponding error with different boundary cut-off methods as the tool penetrates the beds from top to bottom. It can be seen when the cut-off distance is fixed to 8 m, the simulation result of tradition method is consistent with the result of fine interpretation. When the cut-off distance is reduced to 4 m, the tradition method drops quickly in computation precision. Especially at lateral depth of 25-170 m, the relative error of traditional boundary cut-off method is usually above 5%. In contrast, the relative computation error of the self-adaptive method is less than 0.5%, validating the accuracy and stability of the new cut-off method.

To further demonstrate the advantages of the new method, Fig. 7 compares the total layer numbers needed for the forward model when using the two boundary cut-off methods. Compared with traditional method, the new method can adaptively cut off the encountered formations and therefore the total layers of the forward model significantly reduce. By using the new method, the modeling speed of LWD EM is accelerated more than 1 time. Using the Intel (R) i7-8700 CPU, the new modeling method can calculate the responses of more than 16 000 logging points per second. Therefore, the modeling algo- rithm can satisfy the requirement of real-time interactive inversion of the LWD EM data of horizontal wells.

To determine the bed boundary positions, this section presents an interactive inversion method of LWD EM, and the detailed workflow is shown in Fig. 8. The inversion consists of three key parts, i.e. establishment of formation model according to the data of offset well, manual adjustment of the inversion model, and gradient optimization of bed boundary. The first part involves extracting the resistivity and lithology sequence from the offset/pilot wells to build geosteering/interpretation model. Combining the reference model with the priori constraints from nearby logging points (known inversion results), the initial inversion model is established. Then, using gradient algorithm and 1D forward modeling, the bed boundary is optimized. If the inversion results don’t meet the computational precision, the bed boundary positions are manually adjusted and then the gradient optimization is used. The inversion stops until the simulated responses agree well with the measured data.

The inversion of LWD EM data can be converted into solving the least square between the measured data R_a and simulated response S. The objective function of the inversion meets the following relation:

The updating of the regularization parameter satisfies:

Regularized Gauss-Newton algorithm is adopted to get the minimum of the objective function. At t-th iteration, the bed boundaries near the borehole are determined by the following relation:

For the regularized gradient method, the inversion precision usually depends on the precision of initial inversion values. To accurately predict the bed boundary, the interactive inversion with multiple initial values was taken in this study. The strategy is realized by using the following mode: (1) the bed boundary information of nearby logging point serves as the initial model value; (2) based on the previous inversion result, by using Ciflog software, the formation structure is manually adjusted to build a new initial inversion model; (3) repeat step (2) until the inversion results reach the precision required. When adjusting the initial value, we adhere to the following principles: (1) if the apparent resistivity value is much higher than the resistivity of the formation sequence, the distance to upper or lower formation boundary should be smaller; (2) when the tool is at the high resistivity layer, if the simulated response is smaller than field data, the distance d of the log tool to boundary should be shortened, otherwise d is increased. In general, around three to five interactive adjustments and gradient iterations are adequate to get accurate boundary positions. Taking advantages of the 1D fast modeling algorithm, this interactive inversion is capable of processing tens or hundreds of logging points per second, enough to meet the demand of real-time processing of LWD EM data.

The interactively forward and inverse modeling was used in a horizontal well of Jilin Oil field. The targeted layer is siltstone with thin shale interbeds. The formation structure varies slowly along the horizontal direction, and the relative angles between the borehole and formation boundary are between 60° and 133°. Based on the logging data from an offset well, the formation sequences were extracted and then used to establish the geological model shown in Fig. 9. At lateral distance of 214.6 m, the GR value drops quickly and apparent resistivity increases to 20 Ω•m. Based on these evidences, it was predicted this was the landing position. However, it is hard to determine the tool positions from the resistivity curves after landing. The seismic profile was used to constrain the formation strike, the boundary positions where the tool moved into or out of the targeted zone were determined from the GR logs. Track 1 presents the predicted and observed GR data of the interpreted formation model. It can be seen that the two curves almost overlap together. But, when using this interpretation model, the simulated LWD EM responses are far from the observed data, unable to meet the demand of fine interpretation.

The LWD EM inversion method considering detailed LWD EM tool configuration (i.e. MCR from Schlumberger company) was used to accurately derive the bed boundaries around the borehole. After point-by-point interactive inversion, the final interpreted formation model of the horizontal well is shown in Fig. 10. Compared with Fig. 9, the latter model is smoother and the simulated phase resistivity agrees better with the field data, further verifying the accuracy and reliability of the inversion method. The interpretation results demonstrate that the majority of the horizontal well section lies in the targeted zone. Hence, the whole horizontal well section was fractured. This horizontal well had a daily oil production of 86.5 m³ initially, proving the targeted bed is a pure oil layer.

A real-time forward and inverse modeling method of LWD EM has been developed and used to process the field data to solve the problem of bed boundary prediction in horizontal wells. The following findings have been reached:

By approximating and subtracting the first reflection/transmission waves from the scattered field, the convergence of the Sommerfeld integral method can be increased by more than three orders of magnitude. Compared with traditional integral method, the new method can be used for wells with arbitrary dipping angles and has much higher integral accuracy and efficiency. It is more applicable to the simulation of array electric logging responses.

The self-adaptive boundary cut-off method can significantly reduce the layers of the forward model, and balance the tradeoffs between the computation speed and precision. By combining the new integral method with the boundary cut-off technique, the forward modeling method of LWD EM presented in this paper can compute more than 16 000 logging points per second, satisfying the need of real-time forward and inverse modeling.

By introducing priori information of offset/pilot wells and using the classic responses of LWD EM, the interactive inversion of LWD EM logging in horizontal wells is capable of providing real-time and accurate bed boundaries around the borehole and reducing the uncertainty and non-uniqueness of interpretation results.

A_n—amplitude of spectral field in n-th formation layer, dimensionlss;

C(m)—objective function of the inversion;

d—formation thickness, m;

D—vertical distance of the tool to formation boundary, m;

F^TE,h, F^TE,v—propagation coefficients of TE wave from horizontal and vertical sources within the same layer of receiver;

F^TM,h—propagation coefficient of TM wave from a horizontal source within the same layer of receiver;

h—integral kernel function of magnetic field, A/m;

h^s—integral kernel function of scattered magnetic field;

H—magnetic field intensity, A/m;

H_pq—p-direction magnetic field component radiated by q-direction source, A/m;

H^s—scattered magnetic field intensity, A/m;

$H_{0}^{s}$—scattered magnetic field intensity from 0-th order Bessel function, A/m;

$H_{1}^{s}$—scattered magnetic field intensity from 1-st order Bessel function, A/m;

I—identity matrix;

J—Jacobian matrix;

J₀—0-th order Bessel function;

J₁—1-st order Bessel function;

k—wave number, dimensionless;

k_ρ—integral variable;

k_n,h,z, k_n,v,z—horizontal and vertical wave numbers in n-th layer;

m, n, l—number of layers;

m—parameter vector needs to be inverted;

m_ref—reference model;

N—Total number of layers in the model;

Ο—higher-order infinitesimal function;

R_a—measured apparent resistivity, Ω•m;

R_P33H—phase difference resistivity at operating frequency of 2 MHz and spacing of 0.838 2 m (33 in), Ω·m;

$R_{n, n+1}^{\mathrm{TE,h}}$,${R}_{n, n+1}^{\mathrm{TE}, \mathrm{v}}$—local reflection coefficients of upgoing TE wave by horizontal and vertical sources at z_n;

${R}_{n, n+1}^{\mathrm{TE}, \mathrm{h}}$,$ {R}_{n+1, n}^{\mathrm{TE}, \mathrm{h}}$—local reflection coefficients of upgoing and downgoing TM waves at z_n;

s—attenuation ratio of electromagnetic wave radiated by magnetic source at the bed boundary;

S(m)—simulated response, Ω·m;

t—iteration times;

$T_{n,n+1}^{TM,h}$—local transmission coefficient of upgoing TM wave at z_n;

z—vertical coordinate, m;

z_n—vertical coordinate at the n-th bed boundary, m;

$\tilde{R}_{n, n+1}^{\mathrm{TM}, \mathrm{h}},\tilde{R}_{n+1,n}^{\mathrm{TM}, \mathrm{h}}$—general reflection coefficients of upgoing and downgoing TM waves at z_n;

$\tilde{R}_{n, n+1}^{\mathrm{TE}, \mathrm{h}}$,$\tilde{R}_{n, n+1}^{\mathrm{TE}, \mathrm{v}}$—general reflection coefficients of upgoing TE wave by horizontal and vertical sources at z_n;

σ_n,h, σ_n,v—horizontal and vertical conductivity in the n-th layer, S/m;

χ—regularization parameter;

χ₀—initial regularization parameter;

ρ—radial distance, m;

δ—skin depth, m;

θ—relative dipping angle between tool axis and normal direction of formation, (°);

ω—angular frequency, rad/s;

μ—magnetic conductivity of free space, H/m;

λ—anisotropic coefficient of formation resistivity; ξ₁,

ζ₁—amplitudes of upgoing spectral waves in the integrand of $H_{0}^{s}$ and $H_{1}^{s}$;

ξ₂, ζ₂—amplitudes of downgoing spectral waves in the integrand of $H_{0}^{s}$ and $H_{1}^{s}$;

$\xi_{1}^{1}$,$\xi_{2}^{1}$—amplitude of first reflection or transmission wave;

$\tilde{\xi}_{1}$—amplitude of subtracted upgoing wave;

$\tilde{\xi}_{2}$—amplitude of subtracted downgoing wave;

$\|\cdot\|_{2}^{2}$—find the 2-norm of the vector