Algorithm for inversion of resistivity logging-while-drilling data in 2D pixel-based model

Multi-frequency and multi-component extra-deep azimuthal resistivity measurements with depth of investigation of a few tens of meters provide advanced possibilities for mapping of complex reservoir structures. Inversion of the induction measurements set becomes an important technical problem. We present a regularized Levenberg–Marquardt algorithm for inversion of resistivity measurements in a 2D environment model with pixel-based resistivity distribution. The cornerstone of the approach is an efficient parallel algorithm for computation of resistivity tool signals and its derivatives with respect to the pixel conductivities using volume integral equation method. Numerical tests of the suggested approach demonstrate its feasibility for near real time inversion.


Introduction
Drilling high-angle and horizontal wells sections over the last 30 years has led to significant improvement of logging-while-drilling (LWD) resistivity tools. Advanced multi-frequency and multicomponent extra-deep azimuthal (or directional) resistivity tools transmit in real time to the surface a few tens of measurements with depth of investigation of up to a few tens of meters. Such rich dataset is used for imaging of formation resistivity distribution in the vicinity of well path and provides detecting and mapping reservoir boundaries for geosteering adjustments and updating reservoir model. The inversion-based interpretation is common method for mapping resistivity distribution, because visual analysis of several tens of measurements with different depth of investigation and azimuthal sensitivity is not always possible.
Real-time geosteering decisions are usually based on 1D layer-cake resistivity model which is obtained as a result of data inversion. Azimuthal (or directional) resistivity measurements are provided by transmitter-receiver configuration with the coils that is transverse or tilted to the tool axis. Responses are collected at several equally spaced sectors and these data makes it possible to resolve the azimuth of a bed boundary in simple cases. Provided set of multi-frequency and multi-component azimuthal measurements can be sensitive to several deviated bed boundaries along the well path and cannot be fully described by 1D layer-cake model. This behavior is usually observed in complex geological scenarios, such as faults, pinching-outs, unconformities, and scours. Joint interpretation of bulk and azimuthal resistivity measurements in complex scenarios based on a 2D model is preferable for interpretation.
This paper describes an efficient inversion approach based on the parallel algorithm for computation of resistivity tool responses and its derivatives with respect to the pixel conductivities using volume  [1]. The solution of the inverse problem is performed by the regularized Levenberg-Marquardt algorithm to fit resistivity measurements in a 2D environment model.

2D pixel models
The model is two-dimensional and conductivity distribution does not change along one axis. This axis is the strike axis represented by Y. Conductivity in the plane of the model (XZ) is a piecewise constant distribution within rectilinear grid (figure 1). Tool trajectory is an arbitrary 3D curve. For this model we can perform Fourier transform along the Y axis and turns our problem into a set of twodimensional problems [2]. Each of resulting 2D problems corresponds to its spatial harmonics. This technique allows for calculating faster and using less memory than direct 3D simulation; moreover the resulting set of 2D problems can be effectively parallelized.

Inversion algorithm
In the pixel model, the inverse problem is reduced to finding the optimal distribution of conductivity over pixels ( ⃗). The solution to the problem is based on minimizing the sum of the regularizing term ( ) and the root-mean-square residual ( ) of the simulated ( ) and measured ( ) tool responses, normalized to the tool measurement error Here N is the total number of signals at all measurement points used to solve the inverse problem; is the normalized residual of the i-th signal.
It is advisable to make the conductivity dimensionless before optimization. An arbitrary monotonic function ς (σ) transforming the interval ( , ) into (0, 1) can be used for this purpose. The total variation of the "dimensionless conductivity" was used as a regularizing term (4).
Here r is the regularization coefficient, , is the "dimensionless conductivity" of the pixel located in the i-th column and j-th row of the regular grid (figure 1).
Expression (1) is minimized by the Levenberg-Marquardt method. The n-th iteration of optimization is reduced to solving a linear system of equations to determine the displacement vector of the dimensionless conductivity ∆ ⃗ ( (6) Here ( ⃗ ) and ( ⃗ ) are the gradient and the Hessian matrix of the sum of the normalized residual and the regularizing term; is a non-negative constant that is different for each step. The formulas for calculating the gradient and the Hessian matrix are given below.
Here ⃗ is the vector of residuals , is the Jacobi matrix of ⃗ , is the Hessian matrix for . The contribution to the Hessian matrix from the signal residual is calculated approximately under the assumption of the dominant role of the term 2 in comparison with the term 2 ∑ . The contribution from the regularizing term is calculated explicitly. The expression for calculating the components of the Jacobi matrix is shown below.
Here i is the signal index, j is the pixel index. Thus, solving the optimization problem using the Levenberg-Marquardt method is technically reduced to the sequential solution of a set of direct problems for calculating the values of signals and their derivatives with respect to the conductivity. Since these direct problems cannot be solved in parallel, the only way to achieve high computational performance of the inverse problem is to use an efficient parallelized solver.
(12) Background magnetic fields are calculated exactly using rigorous analytical solution for 1D-layered background model. With this approach, the 2D pixel-based model is treated as an extension of the 1Dlayered model. The anomalous electromagnetic field can be presented as an integral over the excess currents [3][4][5].
where and are the electric and magnetic Green's tensors for current as a source, calculated for 1D-layered background model. Integral equations approach was adapted for two-dimensional modeling by Fourier transform along the strike axis Y. Fourier transformation along Y axis is represented by a wave above the notation of a variable (̃). Here we used symmetric form of Fourier transform: In the spectral domain ( , , ) electromagnetic fields (̃, ̃) are also a superposition of background (̃, ̃) and anomalous fields (̃, ̃) . ̃=̃+̃, (16) ̃=̃+̃.
(17) Applying the Fourier transforms to (4) and (5), we obtain the following integral equations: (21) To find electric field distribution a direct method was used. We discretize and solve integral equation above numerically using LU factorization. Receiving coils of EM tools measure magnetic fields only along their axis thus to calculate magnetic field measured by receivers we multiply equation (3) with vector of receiver's directions (⃗⃗).
⃗ ⃗⃗ • ⃗⃗ = ⃗ ⃗⃗ • ⃗⃗ + ⃗ ⃗⃗ • ⃗⃗ (22) To calculate full magnetic field in receivers we perform inverse Fourier transform of anomalous magnetic field and add the background magnetic field to the result: where ̃( ⃗ ) • ⃗⃗ is the magnetic field measured by the receivers; ⃗ is the coordinate of the receiver in the model plane. Magnetic Green's tensor multiply with vector of the receiver's directions (̃( ⃗ , ⃗ ′ ) • ⃗⃗), can be found using the reciprocity principle. In the spectral domain ( , , ) this leads to the following relation: where ̃( ⃗ ′ , ⃗ , − ) is the electric field induced by the receivers in the background model with unity magnetic moment. Tool signals ( ) are calculated from the magnetic field values in receivers ̃( ⃗ ) • ⃗⃗. To calculate the derivatives of signals with respect to conductivity ( ), it is also necessary to calculate the derivatives of the fields in receivers with respect to conductivity (̃( ). To do this, use the expression below [5]: where ̃( ⃗ ′ ) is the electric field induced by the receivers with unity magnetic moment; ̃( ⃗ , ⃗ ′ ) is the magnetic Green's tensors in the spectral domain ( , , ) corresponds to our two dimensional model.
where ̃( ⃗ ′ , ⃗ , − ) is the electric field induced by the receivers with unity magnetic moment.

Method implementation and performance
Base of the described method an efficient parallel solver (PACK2D) was developed [1]. The solver was optimized for resistivity LWD tool with multiple frequencies and multiple transmitter-receiver configurations in 2D pixel-based model crossed by an arbitrary well trajectory. Since simulation of tool signals in trajectory interval is reduced to a set of independent two-dimensional problems corresponding to different frequency, subintervals and spatial harmonics, the problem set is parallelized. The solver was developed on C++ using Intel®MKL [6], OpenMP [7] and CUDA [8] technology and can perform computation on CPU and GPU.
The inversion grid differs from the computational grid one. In the forward modeling each twodimensional problem solved on its own non-regular grid in the model plane. Forward modeling grid size is optimized to achieve necessary accuracy depending on frequency and transmitter to receiver spacing. The continuous equations (10-26) is discretized by the Galerkin method on this non-regular grids. Inversion grid is rectilinear and non-regular grids used in the forward modeling is nested in it.
As an example consider a synthetic dataset of extra-deep azimuthal and bulk resistivity measurements by DAR and EDAR tools [9] in the anticline model (figures 2, 3). The model parameters are shown in table 1. Azimuth and dip in the table 1 are the angles describing the tool trajectory. Azimuth is the angle between the trajectory and the axis Y and dip is the angle between the trajectory projection on the plane XZ and axis Z.  Synthetic data is generated for the trajectory interval of 20 m with 21 measure points using in-house code based on surface integral equations method Pie2d [10,11]. Then inversion using the described above algorithm on this data is performed. Initial model for the inversionis two-layered with the same up and down resistivity ( 1 , 2 ) as in the anticline model. Results of inversion of the DAR and EDAR tools signals on 4900 pixels grid are shown below.     The inversion time for 20 meters interval with DAR and EDAR signals at four frequencies is six hours using workstation with Xeon E5-1620 V3 processor.

Conclusions
An efficient algorithm for inversion of resistivity tool responses in 2D pixel model crossed by an arbitrary trajectory is presented. The approach base is signal misfit minimization with Levenberg-Marquardt algorithm using volume integral equation method for forward modeling. The achieved inversion speed is a few hours per 20 meter trajectory interval using standard mid-level desktop hardware. To achieve real-time computation performance, the algorithm should be ported on HPC cluster.