Fast hierarchical inversion for borehole resistivity measurements in high-angle and horizontal wells using ADNN-AMLM
1 Fast hierarchical inversion for borehole resistivity measurements in high-angle and horizontal wells using ADNN-AMLM
Yizhi Wu a,b,c , Yiren Fan a,b,c, * a School of Geosciences, China University of Petroleum, Qingdao 266580, China b Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China c CNPC Key Laboratory for Well Logging, China University of Petroleum, Qingdao 266580, China
Abstract
With the rapid development of deep learning, intelligent schemes have been gradually introduced to solve various inverse nonlinear problems. In this paper, we combine an efficient adaptive deep neural network (ADNN) framework with an adaptive modified Levenberg-Marquardt (AMLM) algorithm based on a three-layer inversion model to exact the formation resistivity and invasion depth from the measurements of array laterolog. The ADNN presented in this paper can achieve the 2D/3D fast forward modeling of array laterolog. The AMLM algorithm and a hierarchical inversion scheme are adopted to improve the anti-noise ability and convergence in complex logging environments, as well as achieve the fast and accurate reconstruction of longitudinal resistivity profiles in high-angle (HA)/horizontal (HZ) wells. Numerical simulations show that ADNN-based forward modeling only takes 0.021 s for each logging point, and the maximum relative error is less than 2%. The three-layer inversion model can eliminate the effect of the surrounding bed and improve the inversion accuracy in thinly layered formations. The AMLM inversion algorithm can effectively suppress the influence of noise, and takes only 10 steps to achieve convergence.
Key words
High-angle/horizontal (HA/HZ) wells; Array laterolog; Inversion; Adaptive deep neural network(ADNN); Adaptive modified Levenberg-Marquardt(AMLM) algorithm
In the drilling process of high-angle (HA)/horizontal (HZ) wells, the pressure difference between the drilling fluid and permeable layer will lead to the radial invasion of mud filtrate. Due to the differences in mud filtrate and the formation resistivities, the information obtained by dual lateral ∗ Corresponding author. School of Geosciences, China University of Petroleum, Qingdao 266580, China. E-mail address: [email protected] (Y. Fan). logging may not match the true resistivity of the formation (Liu et al.,1993; Mezzatesta et al.,1995; Zhao et al.,2019; Zhang et al.,2000), leading to inaccurate reservoir judgments. Therefore, Schlumberger developed a new high-resolution array lateral tool, HRLA (Itskovich et al.,1998), which provides six apparent resistivity measurements with different depths of investigation (DOIs) and is widely used in the evaluation of oil and gas reservoirs with HA/HZ wells. However, in thinly formations, the measured laterolog arrays are still seriously affected by the relative dipping angle between the well trajectory and the formation, drilling fluid invasion, the surrounding bed and other factors, thus, the true formation resistivity cannot be accurately estimated (Frenkel et al.,1995; Galli et al.,2002; Ni et al., 2017; Hu et al., 2019). Therefore, accurately extracting the real formation resistivity from the apparent resistivity information is important in evaluations of oil and gas reservoirs. Initially, to solve the resistivity distortion problem, a correction method was widely used by interpreters to eliminate the effects of boreholes, surrounding beds and invasion (Griffiths et al.,2000; Phelps et al.,2007; Maurer et al.,2009). However, this technique, which is based on single-factor superposition correction, is only applicable to thick layers in vertical wells (Abdel et al.,2011; Wu et al.,2019). Therefore, geophysical inversion technology that can extract the true resistivity of formation was introduced. Inversion techniques are an excellent way to process HRLA data since they allow for the accurate and simultaneous consideration of both the formation resistivity distribution (radially and vertically) and the dipping angle. Due to the strong nonlinear relationship between the logging response and formation model, it is difficult to extract true values of formation parameters (invasion depth, invasion zone resistivity and undisturbed formation resistivity). Therefore, based on least squares theory, researchers applied the deterministic algorithm in laterolog array data processing and successfully reconstructed the formation resistivity profile (Hakvoort et al.,1998; Wang et al.,2009;Yao et al.,2010; Olabode et al.,2014; Hu et al.,2018; Wang et al.,2019). The objective of this algorithm is to obtain the global minimum through continuous iterations in the direction of gradient descent based on the objective function (Anderson et al.,2001; Yamashita et al.,2001). However, the deterministic algorithm has some disadvantages that cannot be ignored in 3D formation environments: (1) the 3D forward is time-consuming; (2) the calculation of the Jacobian matrix is memory-consuming;(3) the convergence and global optimization ability of the algorithm are strongly dependent on the initial values (Sergio et al.,2016; Li et al., 2019). In addition, with the rapid development of intelligent algorithms, a new perspective for processing logging data has emerged. Compared with deterministic algorithms, intelligent optimization methods can not only eliminate the Jacobian matrix but also achieve global optimization based on an arbitrary initial value, which greatly improves the data processing efficiency (Ewa et al.,2015; Zhu et al.,2019; Feng et al.,2019; Wang et al.,2019; Li et al.,2020). Such algorithms include simulated annealing algorithms, differential evolution algorithms, particle swarm optimization algorithms and Bayesian algorithms. However, in all of these algorithms, forward modeling is still the key restricting factor for the inversion efficiency. Therefore, some researchers have introduced deep learning (DL) into geophysical inversion (Jin et al.,2019; Raissi et al.,2019; Bai et al.,2020). As the foundation of DL, deep neural networks (DNNs) technology has a strong learning ability and can be used as to establish surrogate models for forward simulation (Wu et al.,2019; Li et al.,2019; Kinjal et al.,2019). However, various inversion problems must still be solved, including how to obtain accurate logging databases, how to build efficient neural network frameworks, and how to ensure sufficient inversion precision. In this paper, a new hierarchical joint inversion strategy is proposed for array laterolog measurements. Based on the 3D finite element method (3D-FEM), an efficient adaptive deep neural network (ADNN) framework is constructed, and the adaptive modified Levenberg-Marquardt (AMLM) algorithm with cubic convergence is introduced, which ensures the speed and accuracy of inversion for array lateral measurements. To simplify the workflow, a domain decomposition method is introduced to decompose the large logging data set into several subregions, which improves the inversion efficiency. The paper is organized as follows. In section 2, the principles of array laterolog, the 3D-FEM modeling method and the ADNN fast modeling method are introduced. In section 3, the AMLM algorithm, regional decomposition technology, uncertainty evaluation method and specific hierarchical inversion workflow are proposed. In section 4, synthetic examples are presented to verify the stability and efficiency of the inversion scheme. In the final section, we extend the inversion scheme to anisotropic formations.
In this section, we briefly review the operating mode of the array laterolog tool. Then the signal definition and forward modeling algorithm are presented. The array laterolog tool provides measurements with different DOIs by using multifocusing configurations. Fig. 1 shows the basic framework of the HRLA introduced by Schlumberger in 1998, and it consists of a main electrode A0 , six pairs of shielding electrodes A1 - A6 ( A1' - A6' ), and two pairs of monitoring electrodes M1 and M2 ( M1' and
M2' ). These electrodes are symmetric with respect to the main current electrode A0 . The current injected from the main electrode is focused by the adjacent shielding electrodes and then flows deeply into the formation. By changing the transceiver combinations, we can obtain 6 apparent resistivity curves with different DOIs, namely, RLA0 - RLA6.
Fig. 1. Configuration of the HRLA
The apparent resistivities can be calculated by measuring the potential of the monitoring electrode. Taking the shallow detection mode
RLA0 as an example, the potential difference between two monitoring electrodes on the same side need to be measured, as shown in Eq. (1). In other detection modes, only the potential at a given monitoring electrode need to be measured, as shown in Eq. (1). The apparent resistivities are calculated as follows: ( )( )
M M RLARLA RLA RLA
UR K I = ( ) ( )( ) ( )
1, 2 , 5
M RLAiRLAi i RLAi
UR K iI = = ( ) where, i represents the detection mode; R RLA0 and R
RLAi are the apparent resistivities under different detection modes; K RLA0 and K i are the electrode system coefficients of the HRLA under different detection modes; U M1M2 is the potential difference between monitoring electrodes ( M1 and M2 ); U M is the potential of M
1; and I and I are the emission current intensities of the main electrode A . The key to determining the apparent resistivity of the HRLA is to obtain the potential of the monitoring electrode M Since the HRLA is powered by a low-frequency alternating current, the logging response can be reduced to a direct-current field calculation. In the cylindrical coordinate system ( r , φ , z ), the differential equation can be expressed as (Zhang, 1996): r U U Ur r Rr z z + + = ( ) where U is the measured potential and is the resistivity of the formation. It should be noted that the current and potential are continuous at the interface of the subspace with different resistivities, which satisfies the following boundary conditions: U UU Un n − +− + = = ( ) where, ‘-’ and ‘+’ represents the possible directions at the interface and n is the unit normal vector at the interface. Furthermore, Dirichlet and Neumann boundary conditions must be satisfied at the outer boundary and the surface of the insulated electrode: UUn = = ( ) where, Г and Г represent the outer boundary and the surface of the insulated electrode, respectively. On the surfaces of the main electrode and shielded electrodes, the equipotential boundary conditions should be satisfied: i A i
U dS In − = ( ) where, U Ai is the potential of electrode A i , ρ is the resistivity of the region where the electrode A i is located, Г is the surface of the main or a shielded electrode, n is the unit normal vector at the interface of A i and I i is the emission current intensity. To obtain Eq. (3), the partial differential equation-based problem is transformed into an extreme function-based problem:
62 2 22 =0 i ii
U U UU d I Ur r z = + + − ( ) where, is the solution area, I i and U i respectively are the current and potential of each electrode. Fig. 2. FEM mesh: (a) Mesh of a 3D cylinder model; (b) Mesh structure at x = 0 in the 3D cylinder model To obtain the minimum value of the above function, the finite element method (FEM) is adopted in this paper to discretize the solution area. In the direction of r , the equal-spacing division method is used inside the borehole and nonequal-spacing is used in the formation. In the r- plane, equal-spacing is used for division and the z direction is divided with nonuniform spacing. In this way, the solution region is divided into a certain number of hexahedrons. However, in the three-dimensional finite element method (3D-FEM), tetrahedrons are widely used because of their high calculation accuracy and adaptability, as shown in Fig. 2. Therefore, we further divide these hexahedra into tetrahedra. Thus, the continuous differential equations and boundary conditions with infinite degrees of freedom are transformed into equations containing finite node variables. Finally, the front-line solution method is adopted to solve the linear equations with large sparse matrix, and the potential values of each node can be obtained: K F = ( ) where, i m ei K K = = , m is the number of elements ; , , TN U U = , , , TN F I I = and N is the number of nodes. Since the number of shielded and returned electrodes can vary when the instrument is in different detection modes, the principle of electric field superposition is used to calculate the real electric field. The objective of the method is to generate 7 subfields by assuming that every electrode in the HRLA independently emits current. It should be noted that the current injected by main electrode A0 is assumed to be 1A, and the current intensities injected by the other 6 electrodes are I A1 ~I A6 , respectively. Finally, seven subfields are superposed into real electric fields: A A A A A A
U U I U I U I U I U I U I U = + + + + + + ( ) where, U ~ U are the distribution functions of the seven subfields. To obtain the six weighted coefficients in Eq. (9), six equations generated based on the six detection modes of the HRLA should be satisfied. Taking Mode 6 ( RLA5 ) as an example, the two monitoring electrodes on the same side should have the same electric potential, the shielded or backflow electrodes should have the same electric potential, and the instrument current should be constant. The corresponding formula set is expressed as follows: M MA A A A AA A A A A A A
U UU U U U UI I I I I I I = = = = = + + + + + + = ( ) where, U M and U M are the potential of monitoring electrodes M M U A ~ U A5 are the potential of shielding electrodes A1 ~ A5 . In this paper, to verify the accuracy of the 3D-FEM algorithm, we compare the results simulated by 3D-FEM with those of the physical simulation software COMSOL. The random formation model is given as follows. The borehole diameter is 8 inches, and the borehole is filled with drilling mud with a resistivity of 0.1 Ω.m. The horizontal and vertical resistivities of the middle layer are 20 Ω . m and 80 Ω . m, respectively. In addition, the dipping angle is 80°. The resistivity of the surrounding layer on both sides is 2 Ω.m. The simulation results are shown in Fig. 3a. The scatter points illustrate the COMSOL simulation results and the solid line is the 3D-FEM simulation results. Fig. 3b shows the mean relative error (MRE) at each test point. The maximum error is less than 0.2%, which reflects the high accuracy of 3D-FEM. Fig. 3. (a)Simulated results of 3D-FEM and COMSOL and (b)MRE between them
Traditionally, inversion is time-consuming because hundreds or thousands of 3D-FEM modeling function invocations are involved. In this paper, we propose a rapid forward-modeling approach based on the ADNN framework. Due to their strong generalization characteristics and learning capabilities, deep neural networks (DNNs) are widely used in DL. Specifically, multilayer neural networks, mature DL frameworks, have been widely applied in natural language processing, image target detection, and data prediction (LeCun et al.,2015). A classic multilayer neural network is composed of three layers, as shown in Fig. 4, including an input layer, a hidden layer and an output layer. The neurons in each layer are connected by weights, biases and activation functions. Generally, the three-layer neural network (TL-NN) framework can approximate a simple continuous function. However, for the forward modeling of array laterolog, which contains multiple inputs and outputs, the learning capability of a TL-NN is too weak, and a TL-NN cannot be used as a surrogate model. Therefore, in this paper, an ADNN is constructed to predict the responses of the HRLA in HA/HZ wells. The specific process is as follows: (1) label the formation model, and create a database; (2) train the database, and verify the accuracy of the ADNN; (3) predict the responses of the HRLA. Fig. 4. Three-layer neural network framework
Fig. 5. Three-layer training or inverted model
When drilling through a thinly layered bed in HA/HZ wells, the effect of the surrounding bed on the array laterolog responses cannot be ignored. Therefore, a three-layer model, as shown in Fig. 5, is considered in DL. There are six parameters in the model, including the dipping angle ( ) , the resistivities of the surrounding and target layers ( Rs and Rt ), and the invasion depth ( Di ), resistivity ( Rxo ), and thickness of the target layer ( H ). To decrease the number of neurons in the input layer, the training model is divided into 10 submodels with different dipping angles. The range of dipping angles is from 0° to 90°, and the interval is 10°. Furthermore, a sand-shale sequence is mainly considered in this paper, and it is assumed that the resistivity of the surrounding bed is low. Therefore, the ranges of parameters used in every submodel are expressed in Table 1. In this manner, the database for each submodel is established based on 3D-FEM. It should be noted that we randomly selected 95% of the data for training and the other data for validation of the network. Table 1 The ranges of parameters used in the trained model Parameter Di (m) H (m) Rt (Ω.m) Rs (Ω.m) Rxo (Ω.m) Range 0 ~ 2 0~5 0.1~200 0.1~50 0.1 ~30 In the ADNN, the loss function
L(x) is defined based on maximum likelihood estimation, as shown in Eq. (11), and the sigmoid function, as shown in Eq. (12), is set as the activation function. ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
1, , , ( ( , , )) 1 1 ( , , ) ni L w b x p y log g w p b p x y log g w p b p xn = = − + − − ( ) ( ) ( ) z = exp -z + ( ) where, g ( w,b,x ) = σ ( w*x + b ) , w and b are the weights and biases, respectively; p is the number of iterations; n is the number of samples; x is the vector of the input layer (here, x= ( Rs,Di,Rxo,Rt,H ) T ) ; and y is the vector of the output layer (here, y = log ( RLA1,RLA2,RLA3,RLA4,RLA5 ) T ) . There are two learning mechanisms in the ADNN, including a feedforward mechanism (FF) and a back-propagation mechanism (BP). The loss function can be calculated through the FF mechanism. Then, the weights and biases of the ADNN will be adjusted through BP based on the chain rule until the loss function reaches a manually selected accuracy threshold, as shown in Eq. (13). ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ,1 ,1 k kk k kk kk k k
L w p b pw p w p w pL w p b pb p b p b p + = − + = − ( ) where, k is the number of hidden layers and α is the learning rate; generally, α< In addition, to accelerate the training of the ADNN, the momentum term is introduced in each iteration, as follows: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , 1 , 11 1, 1 , 11 1 k k k kk k k kk k k kk k k k
L w p b p L w p b pw p w p w p w pL w p b p L w p b pb p b p b p b p − − + = − − − − − + = − − − ( ) where η is the momentum coefficient in the interval (0, 1). Then, an adaptive learning rate is introduced to improve the training accuracy and speed: ( ) ( ) p p + = ( ) ( ) ( ) L w b x p L w b x potherwise + ( ) where, ζ is the adaptive coefficient and α ( p ) is the learning rate for the p th iteration. Finally, a genetic algorithm (GA) is adopted to search the suitable weights and biases through BP, which ensures the accuracy of the ADNN. A detail introduction to GA was provided by Byungwhan (2005). Now, it is convenient to present the complete algorithm as follows. Algorithm of the ADNN Input: Rs , Di, Rxo , Rt and H . Maximum iterations Ite = α = η = Terminates iteration threshold
Goal =
Step 1.
Normalize the input data
Step 2.
Set p = 1, and initialize the weights and biases using the GA Step 3.
Compute
L(w, b, x, p). If L Step 4. Calculate w k (p+1) and b k (p+1) based on (14) Step 5. Calculate L(w, b, x, p+1) based on (11) Step 6. Choose ζ as follows: ( ) ( ) L w b x p L w b x potherwise + Step7. Set α(p+1) = α(p) Step8. Set p = p+1 . If k<=Ite , go to Step 3; otherwise, stop. Step9. Output and save the ADNN architecture In this section, we test the precision of the trained ADNN. The borehole diameter of the testing model is 8 inches filled with drilling mud with a resistivity of 0.1 Ω.m, and the other parameters are expressed in Tabel 2. The measurement point is located at the middle of model. Table 2 Parameters of the testing model Di (m) (degree (°)) H (m) Rt (Ω.m) Rs (Ω.m) Rxo (Ω.m) Model1 0 ~ 1.5 0 3 50 2 3 Model2 0.5 0 ~ 90 3 20 2 3 Model3 0.5 0 0~5 50 2 3 Model4 0.5 0 3 3 ~ 150 2 3 The effects of Di , , H and Rt on the prediction accuracy of the ADNN are analyzed in Fig. 6(a)-(d). The scatter points are the predicted results of the ADNN, and the solid line is the 3D-FEM simulation results. The predicted results of the ADNN agree well with the simulation results of 3D-FEM, and the maximum relative errors in Fig.6 are 1.93%, 0.02%, 0.88% and 0.39%, respectively. In addition, it takes approximately 0.021 s to compute values for one logging point with the ADNN; this computational time is nearly 50 times faster than 3D-FEM algorithm. Fig. 6. Comparison of the simulated results of the ADNN and 3D-FEM in various formations with different (a) Di , (b) , (c) H and (d) Rt/Rxo values The traditional Levenberg-Marquardt (L-M) algorithm can improve the global convergence of inversion because it combines the Gauss-Newton approach with the steepest descent method. However, the observed measurements of array laterolog are complicated by the effects of noise, which can limit the global convergence and the efficiency of inversion. Therefore, the adaptive modified Levenberg-Marquardt method (AMLM) with cubic convergence is proposed to optimize anti-noise performance and convergence. The cost function in inversion is as follows: ( ) ( ) ( ) ( ) m 22 21 i pi C x f x e x x x = = = + − ( ) where, x is the inverted vector of the N th layer model, ( ) 1, 1, 1, 1, , , , , Tn n n n x Di Rxoh Rth Di Rxoh Rth = , T is the transpose of matrix; ‖𝑒(𝑥)‖ is the L2 norm of the difference between the measured and forward modeling values; x p is a reference vector containing the parameters for step p -1; and λ is a regularization parameter. The final term in the cost function is used to suppress measurement noise and eliminate pathological effects in inversion. The purpose of inversion is to obtain an x * value that minimizes the cost function; i.e., x* = argmin{ C(x) }. At every iteration, the trial step d k is executed as in the L-M method: ( ) ( ) ( ) ( ) T k k T k kk k d J x C x J x J x I − = − + ( ) where, k is the iteration number, I is the identity matrix and J ( x k ) is the Jacobian at x k . If J(x) is Lipschitz continuous and nonsingular at the solution, then the convergence of the L-M method is quadratic. However, it is difficult to satisfy the inherent nonsingular characteristics in the inversion of array laterolog with HA/HZ wells. Therefore, Fan (2012) introduced the modified L-M method(MLM) with cubic convergence. The specific modified approach involves executing a pseudotrial step k d at each iteration: ( ) ( ) ( ) ( ) T k k T k kk k d J x C y J x J x I − = − + ( ) where, y k is the solution of ( k +1) th step in LM algorithm, y k =x k +d k and λ k is the regularization parameter. In addition, an adaptive regularization parameter selection method was proposed by Keyvan Amini (2015) to optimize global convergence: ( ) ( ) kkk C xC xk = + ( ) ( ) k k k k C x = ( ) where, ω k is the adaptive coefficient of the regularization parameter, which is calculated as follows : k kk k kk r pr p pmax m −−− = ( ) where, m is a positive constant and r k is the gain coefficient, which is defined as follows : ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) k kk k k k k k k k kk k k k C x C x d dr C x C x J x d C y C y J x d − + += − + + − + ( ) Therefore, the final trial step at the k th iteration s k is : ( ) ( ) k k k kTk k kk T Tk k k k s d dd dmin d J x J x d = + = + ( ) where, is the constraint factor, generally, >1. To verify the robustness of the AMLM method, an invaded formation model with Di = 0.5 m is given in this paper. The resistivities of the invaded and uninvaded zone are 2 Ω.m and 20 Ω.m, respectively. The inversion parameter vector is x = ( Di , Rxo , Rth ) T , and two initial value vectors are x = (0.8, 50, 100) T and x = (0.3, 10, 50) T . T he cost function varies with the number of iterations and is shown in Fig. 7. Notably, Fig. 7a shows that the AMLM algorithm will converge faster than the L-M method when the initial value is close to the true value. Additionally, Fig. 7b shows that the AMLM algorithm improves the global convergence when the initial value is 25 times larger than the true value. Figure. 7. Evolution of the data fitting (misfit) norm using the (a) AMLM and (b) LM methods In this paper, a hierarchical joint inversion scheme is proposed to simplify and accelerate the inversion of measured laterolog arrays in HA/HZ wells. In the first stage of inversion, 3D-FEM modeling is replaced with ADNN modeling to obtain the first inverted results and reduce time consumption based on the AMLM method. Then, in the second stage of inversion, as a strict forward algorithm, 3D-FEM is used to guarantee the precision of the inverted results. Because the initial values in the second stage are close to the true initial values, this step is also not time consuming. In this work, we adopt the domain decomposition method to divide the entire logging data set corresponding to long depth intervals that include a significant number of formation layers into many windows, as shown in Fig. 8. It should be noted that the three-layer model is included in each window. The thickness of the middle layer can extend to the boundaries, which can be determined by the inflection points of measurements, gamma rays or spontaneous potentials. In addition, the thicknesses of the upper and lower layers are set to semi-infinite in the inversion process. To improve the inversion precision, the inverted results for the middle and lower layers in the i th window are set as the initial values of the upper and middle models in the ( i +1) th window. The initial values of the lower model in the ( N +1) th window are artificially assigned. In particular, the inversion process can be repeated two or more times until the stopping iteration conditions are satisfied. After the j th processing step for a whole data set, the inverted results will be regarded as the initial values of the ( j +1) th processing step to improve the accuracy of the results. The inversion workflow is shown in Fig. 9. Figure. 8. Domain decomposition scheme Figure. 9. Workflow of the hierarchical joint inversion scheme In this paper, the covariance matrix is adopted to evaluate the uncertainty of the inversion results, and the formula is as follows: ( ) ( ) * * * * T E x E x x E x = − − ( ) where, x * is the final solution of inversion and E denotes the expectation value. When the noise satisfies the Gaussian distribution, Θ can be approximated by the Fisher information matrix: ( ) ( ) 11 * 1 * Ts d J x J x −− − + ( ) where, ξ = θ d -2 λ * ; J ( x * ) is the Jacobi matrix corresponding to the solution x * ; s and d are the numbers of inverted parameters and data points, respectively; Θ s = I s ; Θ d = θ d2 I d ; θ d is the standard deviation of the data points; I d and I s are identity matrices I s ∈ R s×s and I d ∈ R d×d , respectively. The Jacobi matrix can be decomposed by the singular value decomposition method as J = USV T . Additionally, Eq. (26) can be simplified as follows: Tm V S V + ( ) where, I m is an identity matrix and * is the regularization parameter. Finally, the uncertainty of parameters in the i th layer can be evaluated based on the square root of the diagonal terms of the estimator's covariance matrix: i ii = ( ) In this section, 2D synthetic data generated by HRLA are processed to verify the effectiveness of the new proposed inversion scheme, and the detailed implementations are given. Pseudo random noise of different intensities is added to the synthetic data. The well-known Oklahoma (OK) benchmark models with different relative dipping angles ( ), illustrated in Fig. 10a-b, are represented to show the effectiveness of the inversion scheme. The model has 17 layers in a depth segment of 30 m. The formation includes multilayer invasion and consists of alternating conductive and resistive layers with thicknesses ( H ) from 0.3 m to 4.0 m. The borehole diameter is 8 inches filled with drilling mud with a resistivity of 0.1 Ω.m. The invasion depth Di varies from 0 m to 1.3 m. The resistivity of invaded zones Rxo varies from 0.2 Ω.m to 10 Ω.m. The resistivity of uninvaded zones Rt varies from 0.2 Ω.m to 100 Ω.m. Fig. 10c-d show the corresponding measurements simulated with the 3D-FEM code, and the dipping angles are 30 o and 80 o . From the responses shown in Fig. 10c-d, the apparent resistivities are distorted because of the effects of the surrounding thinly layered beds, e.g., the 2 nd bed. In addition, due to the effect of mud invasion, the apparent resistivities are more than five times lower than the true formation resistivity, even in thick beds, e.g., the 6 th , 8 th and 15 th beds. This characteristic becomes increasingly obvious as the dipping angle increases. Therefore, it is necessary to apply inversion to measurements to reconstruct the true resistivity profile of the formation. In this case, we assume that the locations of formation boundaries are determined from other logging information and only invert the invaded-zone resistivity Rxo , uninvaded -zone resistivity Rt , and invasion radius Di . The inversion process is executed two times for a 30 m segment with 0.1 m sampling interval. The time consumption is 87 s on a single i7-9700 processor with 8 GB of RAM. Fig. 11 shows the inverted results based on the inversion scheme proposed in this paper. Fig. 10. OK model with different (a)Di, (b) Rxo and Rt and array laterolog responses for dipping angles of (c) 30 and (d) 80 Fig. 11 shows that the inverted Rt values are more accurate than the inverted Rxo values . The mean relative errors for Rt and Rxo are 4.11% and 6.68%, respectively. With increasing dipping angle, the logging response is also affected by the surrounding beds, which increases the relative error. In general, the mean relative errors of the inverted results at 30° and 80° are 3.46% and 4.66%, respectively. Overall, the inverted results agree well with the original formation values. It is worth noting that the inversion model adopted in this paper includes invasion in every layer. Due to the existence of multiple solutions, the inverted Di in an uninvaded bed is not zero, but the inverted Rxo and Rt values are very close to zero, e.g., as observed for the 9 th bed. Figure. 11. Model and inverted (a) Di and (b) Rxo and Rt at dipping angles of 30 o and 80 o To restore the complex downhole environment and investigate the effects of noise on the inversion scheme, we assume that the measurements are contaminated by pseudorandom noise and that the inversion parameters include Di , Rxo and Rt . Both the formation boundaries and synthetic data are contaminated with noise as follows: ( ) , 1, ,161 , 1, , ; 1, , j jobs obsmn mn z z jRa Ra m M n N = + = = + = = ( ) where, z j is the boundary position and specifies the intensity of pseudorandom noise. is random noise realized with a uniform distribution in the interval [1, 1]. Based on an OK benchmark model with a 30° dipping angle and a 3% pseudorandom noise added to the boundary positions, Fig. 12 shows the measurements for different intensities. From the left panel to the right panel in Fig. 12, the intensities are 5%, 10% and 25%. The inverted results are shown in Fig. 13. It is obvious that the stronger the noise is, the larger the relative error. Specifically, we observe that Rxo_inv , the inverted resistivity of the invaded zones, is generally more accurate than Rt_inv , the inverted resistivity of the uninvaded zones. When the noisy intensity is 25%, the mean relative errors from top to bottom in the OK model for Rxo and Rt are 7.72% and 9.11%, respectively. Furthermore, the influence of noise on Rt in thinly layered beds cannot be completely eliminated, e.g., as observed for the 15 th bed. The relative error for Rt is 33.64%. Generally, in this case, although the boundary positions and measurements are contaminated with noise, the inverted results based on the hierarchical joint inversion scheme still agree well with the values in the original formation model, which verifies the robustness of the scheme proposed in this paper. Figure. 12. The array laterolog responses for noise intensities of (a) 5%,(b) 10%,(c) 25% Figure. 13. Inverted results for noise-containing data:(a) Inverted Di ;(b) Inverted Rxo and Rt To evaluate the uncertainty of the inverted results after obtaining the global solution, we introduce 25% pseudorandom noise as an example. Fig. 14 displays the uncertainty of the inverted Di , Rxo and Rt . Notably, the uncertainty of Di approaches 1 in some uninvaded formation, e.g., the 1 st , 9 th and 13 th beds. In some beds with large invasion depths, such as the fifteenth layer, the uncertainty of Rt increases, and the precision of Rt_inv decreases. In particular, when the invasion depth approximately 0.4 m, the uncertainties of Rxo and Rt are relatively close. Overall, the uncertainty of Rt is greater than that of Rxo . Figure. 14. Uncertainly evaluation of (a) inverted Di ;(b) inverted Rxo and Rt Because the size of sedimentary particles is irregular and thin interbeds develops, the anisotropic characteristics are obvious in reservoirs. Therefore, in this section, we extend the inversion scheme to 3D anisotropic formations to investigate its efficiency. Fig. 15 shows the responses of the HRLA in an anisotropic formation with different dipping angles, and the dipping angle ranges from 0 o to 90 o . The formation is penetrated by a borehole with a diameter of 8 inches, and the borehole is filled with mud with a resistivity of 0.1 Ω.m. The horizontal resistivity Rh is 10 Ω.m. The anisotropic coefficients in Fig. 15a and Fig. 15b are 2 and 5, respectively. The curves with different DOIs become dispersed with increasing anisotropic coefficient. The differences among the five curves reverse and become positive ( RLA1 > RLA2 > RLA3 > RLA4 > RLA5 ). The reversal angle is located in the range of 40 o -65 o . Fig. 15. Influence of anisotropy on the responses of array laterolog:(a) =2;(b) =5 A three-layer anisotropic model with a dipping angle of 80 is used to verify the suitability of the ADNN approach. In the model, the two surrounding beds are homogeneously isotropic, and the middle bed is homogeneously anisotropic. The thicknesses of the surrounding and middle layers are semi-infinite and 2 m, respectively. The resistivities of two surrounding beds are 2 Ω·m. The horizonal resistivity of the middle layer is 2 Ω·m. Fig. 16 shows a comparison between the data predicted by the ADNN ( RLA _ ADNN ) and the data simulated with by 3D-FEM ( RLA _ ). The measurement point is located in middle of the model. The effect of the anisotropic coefficient on the prediction accuracy is shown in Fig. 16a. The horizontal resistivity is 5 Ω·m, and ranges from 1 to 2. The maximum relative error between RLA _ ADNN and RLA _ is 0.13%. The effect of horizontal resistivity Rh on the prediction accuracy is shown in Fig. 16b. When is set to 2, the maximum relative error is 0.07%. Therefore, it is feasible to regard ADNN as a surrogate model for 3D-FEM forward simulation. Fig. 16. Comparison of the simulated results of the ADNN and 3D-FEM in anisotropic formations with different (a) and (b) Rh In this experiment, we apply the inversion scheme to a simple three-layer deviated anisotropic model, and the details of the implementation are the same as those listed in Sec. 4 . Fig. 17a shows the responses and resistivity distribution of the three-layer model. In the model, the two surrounding beds are homogeneously isotropic, and the middle bed is homogeneously anisotropic. The thicknesses of the surrounding and middle layers are semi-infinite and 3 m, respectively. The resistivities of the two surrounding beds are 2 Ω·m, and the horizontal and vertical resistivities are 20 Ω·m and 100 Ω·m, respectively. The dipping angle is 80°. The inverted results are presented in Fig. 17b. In this case, accurate formation parameters can be obtained after 6 s. From the results, the resistivity profile can be accurately reconstructed, which indicates the robustness of the inversion scheme for 3D formation. Fig. 17. (a) Model and corresponding array laterolog responses and (b) inverted Rh and Rt In this paper, a fast hierarchical joint inversion scheme is proposed to extract true formation resistivities from array laterolog measurements in 2D/3D formation with HA/HZ wells. Based on the simulated results and inversion experiments, the following conclusions can be obtained: (1) A database established with 3D forward simulator is used to develop an ADNN model, which can be used as a surrogate model for 3D-FEM. An adaptive learning rate is introduced into the ADNN architecture to improve the predicted accuracy to 98%, and the model only requires 0.021 s to perform calculations at each logging point. (2) Compared to the LM algorithm, the AMLM algorithm can optimize anti-noise performance and global convergence in inversion; notably, the proposed method will converge within 15 steps, even if the initial value is 25 times larger than the true value. (3) A domain decomposition method considering the three-layer model is combined with the AMLM algorithm to improve the inversion accuracy and avoid time consumption in inversion. This approach is nearly 50 times faster than the traditional inversion scheme. (4) Even at a noise level of 25%, the formation resistivity profile can be accurately reconstructed, and the mean relative error of the inverted Rt values is less than 10%. Moreover, the inversion scheme is robust when applied to anisotropic formations. Although only the HRLA is discussed and processed in this paper, the proposed inversion scheme is also applicable to other borehole measurements. Acknowledgements This work is supported by Natural Science Foundation of China (41674131, 41574118, 41974146, 41904109). References Anderson BI. Modeling and inversion methods for the interpretation of resistivity logging tool response. Delft University Press. 2001. Bai YN, Chen W, Chen J. Deep learning methods for solving linear inverse problems: Research directions and paradigms. Signal Processing. 2020; 177:1–23. Byungwhan K, Sungmo K. GA-optimized backpropagation neural network with multi-parameterized gradients and applications to predicting plasma etch data. Chemometrics and Intelligent Laboratory Systems. 2005; 79(1-2):123–128. Ewa GZ, Robert S, Macidj S, et al. A hybrid method for inversion of 3D DC resistivity logging measurements. Nat Comput. 2015; 14:355–374. Frenkel MA, Mezzatesta A.G. Rapid 2-D inversion of resistivity logging data. In: SEG annual meeting, 8–13 October 1995, Houston, Texas; 1995. Feng J, Ni XW, Yang Q, et al. Research on Array Lateral Logging Real-Time Inversions Based on Hybrid Simulated Annealing Algorithms. Petroleum Drilling Techniques. 2019; 47(5):121–126. Fan JY. The Modified Levenberg Marquardt Method for Nonlinear Equations with Cubic Convergence. Mathematics of Computation. 2012; 81(277):447–466. Galli MT, Gonfalini M, Mele M, et al. Resistivity modeling of array laterolog tools: an application in an offshore Norway clastic reservoir. In: SPE annual technical conference and exhibition, 29 September–02 October 2002, San Antonio, Texas; 2002. Griffiths R, Barber T, Faivre O, et al. Optimal evaluation of formation resistivities using array induction and array laterolog tools. In: SPWLA 41th annual logging symposium, 4–7 June 2000, Dallas, Texas; 2000. Hakvoort RG, Fabris A, Frenkel MA, et al. Field measurements and inversion results of the high-definition lateral log. In: SPWLA 39th annual logging symposium, 26–29 May 1998, Keystone, Colorado; 1998. Hu S, Chen L and Wang J. Fast inversion of array laterolog measurements in an axisymmetric medium. Applied Geophysics. 2019; 16:39–548. Hu XF, Fan YR. Huber inversion for logging-while-drilling resistivity measurements in high angle and horizontal wells. Geophysics. 2018; 83(4):D113–D125. Itskovich GB, Mezzatesta AG, Strack KM, et al. High-definition lateral log - resistivity device: basic physics and resolution. In: SPWLA 39th annual logging symposium, 26–28 May 1998, Keystone, Colorado;1998. Jin YC, Wu XQ, Chen JF, et al. Using a physics-driven deep neural network to solve inverse problems for LWD azimuthal resistivity masurements. In: SPWLA 60th annual logging symposium, 15–19 June 2019, The Woodlands, Texas; 2019. Kinjal DG, Vallega V, Maniar H, et al. A Deep-Learning Approach for Borehole Image Interpretation. In: SPWLA 60th annual logging symposium, 15–19 June 2019, The Woodlands, Texas; 2019. Keyvan A, Faramarz R. A modified two steps Levenberg–Marquardt method for nonlinear equations. Journal of Computational and Applied Mathematics. 2015; 288:341–350. Liu QH. Reconstruction of two-dimensional axisymmetric inhomogeneous media. IEEE Transactions on geoscience and remote sensing. 1993;31(3); 587–594. Li H, Wang L. Fast Modeling and Practical Inversion of Laterolog-Type Downhole Resistivity Measurements. IEEE Trans Geosci Remote Sens. 2019; 57(1): 120–127. Li H, Wang H, Wang L, et al. A modifed Boltzmann Annealing Differential Evolution algorithm for inversion of directional resistivity logging-while-drilling measurements. J Petrol Sci Eng. 2020; 188:1–10. Li H, Liu G, Yang SS, et al. Automated resistivity inversion and formation geometry determination in high-angle and horizontal wells using deep learning techniques. In: SPWLA 60th annual logging symposium, 15–19 June 2019, The Woodlands, Texas; 2019. LeCun Y, Bengio Y, Hinton G. Deep learning. Nature. 2015; 521:436–444. Mezzatesta A.G, Eckard M.H, Strack K.M. Integrated 2-D Interpretation of Resistivity Logging Measurements by Inversion Methods. In: SPWLA 36th annual logging symposium, 26–29 June 1995, Paris, France; 1995. Maurer H, Antonov Y, Corley Y, et al. Advanced Processing For A New Array Laterolog Tool. In: SPWLA 50th annual logging symposium, 21–24 June 2009, The Woodlands, Texas; 2009. Ni XW, Xu GY, Feng JM, et al. Forward response of array lateral logging in anisotropic reservoir of inclined shaft. Fault-Block Oil & Gas Field. 2017;24(5); 637–641. Olabode L, Carlos TV, William EP. Inversion-based petrophysical interpretation of logging-while-drilling nuclear and resistivity measurements. Geophysics. 2013; 78(6):D473–D489. Philip DR, Zhou ZQ, et al. Array processing—a new method to detect and correct errors on array resistivity logging tool measurements. Chinese Journal of Oceanology and Limnology. 2007;25(1); 51–58. Raissi M, Perdikaris P, Karniadakis GE. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics. 2019; 378:686–707. Shafy TA, Fattah A, Corley B, et al. Application of a new multi laterolog logging tool in the Phiops field of Egypt, a case study. In: OMC 10th annual symposium, 23–25 March 2011, Ravenna, Italy; 2011. Sergio R, Ignacio M, David P. A quadrature-free method for simulation and inversion of 1.5D direct current (DC) borehole measurements. Comput Geosci. 2016; 20:1301–1318. Wu ZG, Fan YR, Wang JW, et al. Application of 2.5-D Finite Difference Method in Logging-While-Drilling Electromagnetic Measurements for Complex Scenarios. IEEE Geosci Remote Sens Lett. 2019; 17(4):577–581. Wang GL, Carlos TV, Jesús MS, et al. Fast 2D inversion of large borehole EM induction data sets28