- Email: [email protected]
Two dimensional combined inversion of short- and long-normal dc resistivity well log data
Journal of Applied Geophysics 73 (2011) 130–138
Contents lists available at ScienceDirect
Journal of Applied Geophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j a p p g e o
Two dimensional combined inversion of short- and long-normal dc resistivity well log data Emin U. Ulugergerli ⁎ Çanakkale Onsekiz Mart University, Faculty of Engineering and Architecture, Department of Geophysical Engineering, Çanakkale, Turkey
a r t i c l e
i n f o
Article history: Received 19 October 2009 Accepted 8 December 2010 Available online 23 December 2010 Keywords: Well-logging Direct current resistivity 2D inversion Finite element method
a b s t r a c t A method is presented for the two-dimensional combined inversion of short- and long-normal tool direct current resistivity data with symmetry. The forward problem is solved using the finite element method in the cylindrical coordinates system. The inverse problem is solved using a conjugate gradient technique with the partial derivatives obtained using reciprocity. The parameters were obtained by means of both conjugate gradient relaxation and conventional conjugate gradient method. The solution of this highly underdetermined inverse problem is stabilized using Tikhonov regularization and the scheme yields a blurred image of the subsurface. The scheme is tested using synthetic data and field data. Tests using synthetic data suggest that traces of the horizontal boundaries are delineated in the range of the exploration distance while the resolution of vertical boundaries depends upon the solution regularization. Application to field data shows that additional information is necessary for resolving the resistivity structure when there are low resistivity contrasts between formation units. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Well logging applications produce useful information about geological formations where the thicknesses of formations are larger than the tool's size. Yang and Ward (1984) generalized the investigation results of Daniels (1978) to show the relation between electrodes spacing and layer thickness; if the layer is thicker than ten times the electrode spacing, the apparent resistivity will be closer to the true layer resistivity. If the thickness is less than six times the electrode spacing, estimating true resistivity will be cumbersome and departure curves (e.g. Dutta, 1994) may be employed to estimate the true resistivity values. Departure curves are similar to normalized apparent resistivity curves in surface application of the direct current resistivity (DCR) method. In reality, features such as borehole fillings, invaded zone, layering, and shoulder beds create a complex environment. The apparent resistivity curves represent filtered and simplified images of geo-electrical structures. Multi-dimensional modeling is a convenient way to image complex subsurface structures. As well known for surface measurements, two-dimensional (2D) modeling of the DCR method requires the solution of a boundary value problem in Fourier transformed space (e.g. Dey and Morrison, 1979). However, this requirement is not necessary when radial symmetry is assumed. The axisymmetric subsurface model indicates that the well
⁎ Tel./fax: +90 286 2180018x2207, +90 286 2180541. E-mail address: [email protected]. 0926-9851/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2010.12.004
itself together with the invasion zones and the resistive structures have radial symmetry in the range of exploration and the resistivity variations may occur in both vertical and radial directions. This condition assumes that the structures have concentric ring-like shapes, centered on the well axis. Axisymmetric layout of a well allows one to employ a cylindrical coordinate system. In such a case, either analytic solutions or numerical techniques can be used in the evaluation of electric tool responses. If the concentric ring-like resistivity structures have infinite thickness, the problem is reduced to one-dimension (1D) and an analytical solution will provide the required solution (e.g. Daniels, 1978; Yang and Ward, 1984). If the thicknesses are finite then, either the finite element method (FEM, e.g. Cozzolino and da Conceição da Silva, 2007), the finite differences method (FDM, e.g. Roy and Dutta, 1997; Towle et al., 1988), the integral equation method (IEM, e.g. Zhang and Shen, 1984) or a combination of all of these methods can be employed. In this paper, the FEM with bilinear rectangular element has been used for problem discretization and to solve the axisymmetric earth model. The modeling scheme takes into account only one half of an infinite number of planes in the domain. The use of non-uniform grid spacing allows for modeling large scale structures. The inversion of DCR logging data for subsurface resistivity distribution is a highly under-determined problem. While only the surveyed zone is considered in the inversion process, equations for hundreds of model parameters are solved per DCR reading. Hence, the problem requires regularization. The resulting regularized problem is then solved using the conjugate gradient (CG) method. Two conventional approaches are considered in this paper. First one is
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
CG relaxation (e.g. Mackie and Madden, 1993; Zhang et al., 1995) and will be called as “relaxation” approach in the text. Second solution is yielded via conventional CG (e.g. Zhdanov, 2002; Zhdanov and Tolstaya, 2004) and will be called as “conventional” approach in the text. The definitions are used solely for the purpose of abbreviations. For computational efficiency, the principle of reciprocity is used for the calculation of partial derivatives. Calculating the partial derivatives of potentials with respect to the model parameters for each reading point requires thousands of forward solutions for each inversion iteration. The reciprocity method given by deLugao and Wannamaker (1996) reduces this necessity to performing two forward solutions per reading; one for current sources and another for a unit source at the receiver site. As mentioned previously, simulation of DC logs in axiallysymmetric media was studied extensively. Many of the publications usually focus on departure curves. Some publications related to inversion may be found either in proceedings (e.g. LaBrecque et al., 2006) or special journals (e.g. Pardo et al., 2007). Unlike the induction logs (e.g. Liu and Lin, 2002), interpretation stage of DC logging usually consists of plotting and evaluating the apparent resistivity curves with using departure curves or performing single-log inversion only. The innovative part of the paper is the combined inversion of short and long normal DC data. This approach, then, tested on low-contrast case which is the bane of inversion. The resulting inversion schemes will be evaluated using synthetic and field data. Observable data for short normal (16 in.) tools and long normal (64 in.) tools are calculated for simple but representative borehole environments. The field data considered are from a low-contrast geoelectrical situation. Both data sets are used in the inversion routine. Tests on synthetic data indicated that the single log inversion of short normal data recovers structural pattern in the range of the exploration radius. On the other hand, the inversion of long normal data delineates only the anomalous zone. Besides the single log inversion, it is also demonstrated that joint inversion of short- and long-normal data produces a blurred picture of the subsurface, and that while horizontal boundaries can be delineated in the vicinity of the well, the recovery of structures with low conductivity contrast is a hard task requiring some a priori information. 2. Solution of the forward problem The governing equation for the flow of a steady state electric current in a non-uniform medium (Dey and Morrison, 1979; Mufti, 1976) is given by ∇ ðσ ðx; y; zÞ∇νðx; y; zÞÞ =
∂Q ðx; y; zÞ ; ∂t
ð1Þ
where σ is the conductivity, υ is the electrical potential and Q is the charge density. For the axisymmetric case, Eq. (1) can be rewritten as ∂ ∂νðr; zÞ ∂ ∂νðr; zÞ 1 ∂νðr; zÞ σ ðr; zÞ + σ ðr; zÞ + σ ðr; zÞ r ∂r ∂r ∂z ∂z ∂r =
131
measurement zone of the mesh has an equal spacing in the vertical (z) direction up to the end of the surveyed zone followed by a gradually increase of cell thicknesses. Similarly, the width of the cells gradually increases in the radial (r) direction. The presented test models used 2 in. (≈5 cm) spacing in the vertical direction which mimic borehole logging surveys with 2 in. intervals. The radial intervals vary between 0.125 and 100 in. Application of the FEM to Eq. (2) produces a set of equations, ð3Þ
Gν = Si
where G is the global or stiffness matrix, v is the vector of potentials and Si is the current source vector for the ith measurement. Eq. (3) is solved for each source location. These computations were implemented in MATLAB. Since the program was developed in MATLAB, a builtin-function of back slash (\), based on Gaussian elimination with partial pivoting is used to solve Eq. (3). The code uses four electrodes in the computation to obtain a complete simulation. Therefore, the possible effect of the shoulder bed or upper layers may also be considered. The distance between N and M electrodes is set at 15.24 m (50 ft) and B electrode placed at the surface. The geometric factor is set alike to the conventional four electrodes system to handle the deviation in the near surface measurements, k=
4π 1 1 1 AM − AN − BN
+
1 BM
ð4Þ
Due to the fact that all the distances are very large compared to the AM distance (Fig. 1), the operative distance is AM and is set to either 0.4 m (≈16 in.) or 1.6 m (≈64 in.), representing the short and long normal tools, respectively. Once the potential and geometric factors are calculated, the apparent resistivities are obtained as ρa = k
∇VMN : I
ð5Þ
Fig. 2a compares the forward solution from the present code with the analytical result of Dutta (1994). The 1D model considers mud and two zones; the mud, flushed and uncontaminated zones which have resistivities of 1, 10 and 0.5 Ω m, respectively. The first two zones are 1 and 1.5 m thick, respectively. A 6.5 m thick invaded zone has linear variation with radial distance. Details of the analytical calculation can be seen in Dutta (1994). Additionally, Fig. 2b presents the apparent resistivity curves of short-normal logging for the case where there is no borehole. The target layer is 80 in. (≈2 m) thick and has a resistivity of 10 or 100 Ω m in models 1 and 2 respectively. Models 1 and 2 have background resistivities of 100 and 10 Ω m, respectively. In both cases the borehole radius is 8 in. (≈0.2 m). Model 2 also has a
ð2Þ
Io δðr Þδðz−zc Þ 2π r
where r is the radial distance from well and z defines the depth. r = 0 is the centre of the well while z = 0 is the surface. Current, Io is injected at depth zc in the well (r = 0). v is the potential created by the current. FEM with bilinear rectangular elements with four nodes (e.g. Kwon and Bang, 2000) is used to solve Eq. (2). Various types of boundary conditions can be applied to edges of the mesh. Since the finite element method uses all nodes in calculations, I did not impose any external values to boundaries but set all borders away from the region of interest. Hence, the
Fig. 1. Basic model for borehole (not scaled). Right side is the calculation domain. Rings on toll indicate the electrodes.
132
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
optimization problem. The problem is underdetermined and needs to be regularized. The objective function to minimize therefore incorporates stabilizing functional (Tikhonov and Arsenin, 1977), viz: F = ϕðpÞ + α sðpÞ = min
Fig. 2. (a) Comparison of analytical (Dutta, 1994) and numerical (FEM) solution of Eq. (2) for 1D model. Resistivity of invasion zone obeys linear variation with radial distance. Axes are resistivity (Ω m) and transmitter–receiver spacing (AM). (b) Sample curves for the case with a 80″ (~ 2 m) thick layer. Both cases the borehole radius is 8″ (~0.2 m). Model 2 also has a transition zone of 50 Ω m, radius of which is 40″ (~1 m).
transition zone of 50 Ω m, the radius of which is 40 in. (≈1 m). Since the ratio of layer thickness to tool length (80/16) is smaller than 10 in Fig. 2b, the apparent resistivity curve deviates from the true value. An inversion technique will be employed to obtain the formation resistivity value in such a case instead of using the traditional correction charts.
3. Methodology for 2D joint inversion of short- and long-normal tool responses The inverse problem is to recover the conductivity distribution around the well from the typically few measured dc resistivity welllogging data. Since the mesh defines the electrical resistivity structure permeated by the well, the logarithms of conductivities of each cell in the 2D mesh are taken as the model parameters (p) while the logarithms of apparent resistivities serve as data (d) in our nonlinear
ð6Þ
where ϕ(p) is a measure of fit between observed and calculated data. α is a regularization parameter which provides a balance between ϕ(p) and stabilizing functional s(p). The solution to Eq. (6) is subject to various publications. The outlines of the researches focuses on that the selection of regularization functional, s(p) and that the investigation of algorithms to compute regularized solution. Starting from the first one, a wide range of stabilizing functional is available in the literature (e.g. Zhdanov, 2002). Depending on the nature of the inverse problem taken into consideration, different stabilizer can be employed in Eq. (6). Main expectations from a stabilizer functional can be summarized as follows (e.g. Zhdanov, 2002; Zhdanov and Tolstaya, 2004); it should (a) stabilize the illposed inverse problem, (b) reduce the number of possible models responds of which fit the data, (c) generate a compact image of the geo-electrical model and (d) not smear out the small features. Equivalent expressions in practice are all fallows, in all multidimensional modeling studies, the lack of enough information content in data, which leads to solution to an underdetermined problem (statement a), urges a priory model to be used. Then, the Tikhonov regularization theory seeks a model close to that priory (or reference) model (statement b). The most cases a reference model cannot be enough to get an appropriate solution. Thus some additional constraints can be provided by, for instance, the minimum norm of the second order finite difference operator, also called maximum smoothness stabilizing, through which neighbor parameters are forced to be correlated with each other. This constraint results in a model presents smooth resistivity distribution. As an alternative to smooth variation, focusing or minimum support (MS) or minimum gradient support stabilizing functionals could be considered where detailed structural information (statements c and d) is required (e.g. Portniaguine and Zhdanov, 1999). In our case, the sharp boundaries are needed to describe the complex resistivity distribution. Among a large family stabilizers, MS were proved to be effective to recover a model consists of discontinuous images (Zhdanov and Tolstaya, 2004). The stabilizing functional, s(p) is given as 2 sðpÞ = ‖Wm Ws p−pref ‖
ð7Þ
where ||.|| represent the two-norm, Wm is a model weighting matrix, diagonal matrix Ws is the stabilizing operator. p and pref are the logarithms of the current model and the reference model, respectively. The weighting for MS stabilizer, 0
1
B C 1 B C Ws = diagB 1 = 2 C; 2 @ A 2 +ε pi−1 −pref
ð8Þ
where ε is the focusing parameter. It is calculated, at each iteration, by dividing the minimum conductivity value in the mesh by a hundred. Note that Ws depends on the parameters from previous iteration. Once decided the stabilizer, the second step requires the selection of inversion algorithm. Rodi and Mackie (2001) compared nonlinear conjugate gradient (NLCG), Gauss-Newton (GN) and linear conjugate gradient (LCG) algorithms to compute regularized solution of the 2D magnetotellurics inverse problem. They concluded that the conjugate gradients (CG) based algorithms are superior to the GN algorithm. Later, Candansayar (2008a) included the comparison of different stabilizing functionals to the same inverse problem. He stated
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
weakness and superiority of the algorithms in terms of stabilizing functionals. Both studies declare that the CG based algorithms is favorable for the multi-dimensional inverse problems. Before going further, first term in Eq. (6) is also misfit functional is given as 2
ϕðpÞ = ‖Wd ðd−F ðpÞÞ‖
ð9Þ
where F(p) is the forward solution given by Eqs. (3) and (5). Diagonal matrix, Wd consists of positive weightings which are a function of the measure of the uncertainties of the observed data. Since we do not have any information on data errors, Wd set identity and is dropped from Eq. (9) for simplicity. Linearization of F(p) yields (e.g. Haber and Oldenburg, 1998), 2
ϕðpÞ = ‖ J∇p−Δd‖ ;
ð10Þ
where J is the sensitivity matrix, Δp are the model updates, the vector Δd is the discrepancy between the logarithm of observed and logarithm of calculated apparent resistivities. For the joint inversion problem, the partitioned matrix, J = ½ JSN JLN
T
ð11aÞ
and the partitioned vector, Δd = ½ΔdSN ΔdLN
T
ð11bÞ
represent the short-normal and long-normal contributions to the system of equations. There are many procedures to obtain the solution to Eq. (6) but the CG approach (e.g. Pidlisecky et al., 2007; Zhang et al., 1995; Zhdanov, 2002) is adopted in this study. Basic idea is that one needs to find Δp which, in turn, minimize the Eq. (6). The model updates are sought either in one step (e.g. Zhdanov, 2002) or by means of relaxation at each iteration of the inversion (e.g. Mackie and Madden, 1993). These procedures, in this study, will be named as the conventional and the relaxation solution, respectively and details of each will be given later in the text. In any case the solution to Eq. (6) yields perturbation (or updates) for parameter pi. pi +1 = pi + b Δp
ð12Þ
where b is a line-search parameter (e.g. Pidlisecky et al., 2007) which controls the magnitude of the update vector. Besides the scaling by b, the Δp values are truncated as:
These bounds are selected arbitrarily and prevent large steps in the model update process. Finally, the model weighting matrix Wm is defined as (e.g. Portniaguine and Zhdanov, 1999), k T Wm = diag J J
inversion iteration and maximum number of relaxation iteration, respectively. The empty lines at right panel indicate that the corresponding lines at left panel are excluded. Relaxation solution
Conventional solution
For k = 1 to ninv J calculation and Update W Weight J r = d − F(pk − 1) Io = JTr − αWTW(pk − 1 − pref) r0 = Io; dp0 = 0; b = 1 For j = 1 to line search rk = JTr − αWTW(pk − 1 − pref) For i = 1 to nrex beta = rTi ri/rTi − 1ri − 1 Ii = r i + beta Ii − 1 step = rTi ri/ITi [JTJ + αW] Ii dpi = dpi−1 + step I1 ri = ri − 1 − step L End loop i pk = pk − 1 + b dpk fk = || d − F(pk)||2/N If fk b fk − 1 ; break loop j Else : b = b/1.5; α = α/2 End loop j If fk b treshold ;break loop k End loop k
For k = 1 to ninv J calculation and Update W Weight J r = F(pk − 1) − d Io = JTr + α W (pk − 1 − pref) r0 = Io; b = 1 For j = 1 to line search rk = JTr + α W (pj − 1 − pref)
ð13Þ
where subscript k controls the weighting level. Since the inverse of the diag(JTJ) affects the partial derivatives, k values controls the weighting level. If k N 1, the parameters distant to the well are overweighted. On the other hand, if k → 0, those parameters closer to the well are over-weighted. In this study k was set at 0.5. 4. Conjugate gradient solution of the inverse problem As stated before, two different approaches were considered to carry out the solution to Eq. (6). The respective algorithms given in the later part are modified from Mackie and Madden (1993) and Zhdanov (2002). The ninv and the nrex define maximum number of
beta = rTkrk/rkT − 1rk − 1 Ik = rk + beta Ik step = rTkrk/ITk[JTJ + α I] Ik
pk = pk − 1 − b step Ik fk = || d − F(pk)||2/N If fk b fk − 1 ; break loop j Else b = b/1.5; α = α/2 End loop j If fk b treshold ; break loop k End loop k
Innermost loop can end when they reach their predefined number of iteration (nrex). Outer loop can end either when they reach their predefined number of iterations (ninv) or the misfit reaches to predefined threshold value or insignificant improvement in the model. The data fit, fk is measured by the quantity T
misfit = fk = ΔG ΔG = N
ð14Þ
where N is number of the observations. ΔG is a difference vector of logarithm of the observed and calculated data. Being a stopping criterion, target misfit value needs to be defined by the user. Unfortunately there is no way to select appropriate value since our problem is severely underdetermined and a variety models can provide the same fit to the data being inverted. Any model with a small misfit value (b1.E−06) may be accepted as satisfactorily explaining the inverted data. Insignificant improvement in the model means that the variation of successive misfit values, dmisfit =
if Δp N 2.0 then Δp = 2.0 if Δp b − 2.0 then Δp = −2.0.
133
misfitk−1 −misfitk misfitk
ð15Þ
is less than 1.E−5. Any priori information about any part of the resistivity model structure can reduce the number of iterations and move the result towards the expected model. Because the finer cell widths are used near the well, to recover their conductivities from inversion will be very difficult, thus, the borehole itself is defined with its mud conductivity. Knowledge of the mud conductivity can be incorporated into the initial model and then the left-most cells as priory information. 5. Regularization parameter The solution of Eqs. (12) and (13) requires the determination of regularization and line-search parameters. Each of them needs to be sought during the iterative process, the initial value for the regularization parameter is given, αo = ‖ JΔp−Δd‖ = ‖Wm po −pref ‖:
ð16Þ
If it does not reduce the misfit it is reduced by two, αn = αo/2.
134
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
A line-search is employed in relaxation solution only and performed in a simple manner as given in Pidlisecky et al. (2007). Initially, the value of b is set to unity and subsequently decreased using the relation bi = bi−1 = 1:5:
ð17Þ
6. Sensitivity matrix Calculation of the sensitivity matrix. J requires the derivative of Eq. (3), which, in the general case, is (Rodi, 1976), G
∂ν ∂G =− ν: ∂σ ∂σ
ð18Þ
The right-hand side (RHS) in the Eq. (18) is a new source term for the forward solution. Without employing the principle of reciprocity the most time-consuming part of the inversion process would be the calculation of partial derivatives. The method given by deLugao and Wannamaker (1996) reduces the calculations to two forward solutions per reading; one for current sources and another one for unit source. deLugão et al. (1997) define the forward solution with unit source as, GR = U1
ð19Þ
where U1 is a vector which has a value of unity at the receiver location. The solution vector, R contains the model responses in all nodes. Multiplying the appropriate responses by the cell parameters gives the derivative values. In the adopted FEM, the cell parameters are the element matrix obtained from the shape functions of each node. The derivatives are multiplied by the ratio of apparent resistivity to conductivity of relevant cell to obtain a logarithmic quantity. Fig. 3a and b shows the sensitivity zone for both the short- and long-normal tools. The receiver electrode (M) is at 862 in. (≈9.7 m) while a current electrode (A) is at 878 in. (≈22.3 m) and 926 in. (≈ 23.5 m) for short- and long-normal tools, respectively. The derivatives are normalized by their maximum value and presented for a selected receiver–transmitter position. The absolute values were used in the graphs. Selecting 0.03 arbitrarily, the maps indicate that 16 in. toll can comb almost twice the tool length while the 64-inch tool can only comb approximately 1.5 times the tool length. Note that this limit does not apply to the investigation distance (δ) but it should also adhere to a similar range ratio. One may stipulate this condition as, δ16 δ ≥ 64 : 16 64
ð20Þ
If one defines the investigation distance as a radius that data begins to sense the outer unit then, two-unit model can be used to examine the radius of investigation (ROI). Fig. 4 presents short (+) and long (triangular) normal curves for two-unit model. Horizontal axis is the radius of the interface between first and second unit. The number indicates the resistivity ratio of inner unit to outer unit (ρ1/ ρ2). Low contrast case (1/0.1 to 1/10) 64 in. toll sense ~ 24 in. (~60 cm) radius while 16 in. tool reaches to ~6 in. (~ 16 cm). Behavior of the curves becomes cumbersome in case of high contrast but ROI remains same. The results are in accord with the Eq. (20). 7. Application of the inversion algorithm to synthetic and field data The recovered geo-electrical sections from the joint inversion of short- and long-normal logging data are presented here for both the synthetic example and for real data. Although the Eq. (2) depends
Fig. 3. Normalized sensitivity map for short-normal tool (left) and long-normal tool (right). The resistivity of the medium is 1 Ω m.
upon the variation of conductivities, the sections are presented as resistivities. The first example considers a test model which is given in Fig. 5. The model considers the thin layer problem. The anomalous section has two blocky units of 0.0033 S/m (300 Ω m) and a layer of 0.033 S/m (30 Ω m) immersed in a 0.01 S/m (100 Ω m) host medium. The thickness of the each layer is 100 cm (~39.5 in.) while the width of each resistive block is 100 cm. The conductive layer (mid unit) extends to infinity. The ratio of layer thickness to tool length is ~2.4 for short-normal tool but ~0.6 for long-normal tool. Such combinations result in deviations in apparent resistivity curves. The long-normal curve has very little signatures for the layers while short normal curve represents all units. The well has 20 cm diameter and is filled with a mud of 10 Ω m (0.1 S/m) resistivity. The depth position of the first reading is 178.5 m in the test models. The number of data is 142, i.e. 71 readings for each tool. The noise adding is tricky in single reading data set. Short wave length variations (i.e. random noise adding to each reading) cause a recovery of artificial features in the vicinity of the borehole. Thus a trend is added to both data set. Minima and maximum values of the linear trend is 0.1 and 7.1, respectively. Once synthetic data were created, they were fed to the program which creates the initial model automatically depending on the data and survey parameters. The
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
135
are plotted versus location of their current electrodes (A in Fig. 1). The right panel is the recovered model. The vertical scale is the same as the left panel and white lines define the real boundaries of the sought model. The recovered model presents an expected result at the border of the radial exploration range (~ 1 m); thereafter, the image of model deteriorates in terms of the cell conductivities. Both resistive and conductive units appear over focused since the algorithm uses fixed value for focusing parameter. The over focusing caused separation in the horizontal interfaces of layers. Additionally, the inversion could not provide the exact conductivity values for each layer. If one defines the ratio for real and calculated and conductivities as r=
Fig. 4. Short (+) and long (triangular) normal curves for two-unit model. The number indicates the resistivity ratio of inner unit to outer unit (ρ1/ρ2). Horizontal axis is the radius of the interface between first and second unit.
initial model has 54 and 296 cells in x and z directions, respectively. The conductivity of leftmost 20 cm of the model is assumed to be known as priory information and set 0.1 S/m (10 Ω m) but not excluded from inversion. Only the anomalous zone of 3150 cells in total, between the first and last reading points, was considered in the inversion and rest of the model was kept fixed and excluded from inversion. After each iteration, the conductivity values at the border of anomalous zone are extended to mesh borders not to create any discontinuity between shoulder and calculation zones. For the CG routines the focusing parameter was fixed to 0.003 (Zhdanov and Tolstaya, 2004). The misfit threshold set 1.e−5. First test was run for the relaxation approach. The inversion of the test model data started with misfit value of 0.0063357 and stopped after 12 iterations. The misfit for the result was 6.3.e−6 . All the processes, including plotting for inter-steps took 192.35 s CPU time on the Intel Core2 @ 3 GHz processor. The left panel in Fig. 5 presents the observed (marked symbol) and calculated (solid line) resistivity curves for short normal (+) and long normal (o) tools. Both curves
ðσreal Þ ðσcalc Þ
it varies from 0.2 to 3.4 for the anomalous zone. Second test was run with the conventional solution approach. The inversion stopped after 10 iterations. The misfit for the result was 6.3. e−6. All the processes, including plotting for inter-steps took 164 s. Note that both approaches end the iterations with the misfit reaching to threshold value. This case, the ratio for real and calculated conductivities varies from 0.2 to 2.8 for the anomalous zone. Same as before the horizontal interfaces of layers are determined. Recovery is slightly better than previous case (Fig. 6). Variation of the objective functional (Eq. (6)) vs. CPU time is given in Fig. 7. The relaxation approach decreases objective functional very rapidly in first few iterations but then it appears that it slows down. On the other hand, the conventional solution decreases it steadily and the misfit functional (Eq. (9)) reaches to the threshold value earlier than the relaxation approach. Third tests were run on single log data. The long normal log presents smoothed curve, but the short normal data present detailed structural features, therefore, the inversion with the short normal data set results in better recovery than the long normal data (not presented). Thus, the long normal data were selected for testing the programs. Conventional approach end in 14 iterations and reduce the misfit from 0.0059 to 9.5574e−6 in 116.5 s of CPU time. The relaxation approach, on the other hand, did 8 iterations and reached the misfit value of 9.52e−6 in 69 s of CPU time. Both approaches could not recover the structural pattern, but only delineated the anomalous zone (Figs. 8 and 9). The visual inspection of both results indicates that relaxation approach recovers slightly better structural information since the trace of the conductive layer is clearer. Overall result of these tests points outs that the relaxation approach produces slightly better result than the conventional solution approach does and that the focusing parameter needs to be obtained through the iterations. One suggestion is to use percent of the minimum conductivity value in the model. ε = minðσ Þ = 100
Fig. 5. Observed (marked) and calculated (solid line) curves for short normal (+) and long normal (o) tools for two-blocks and a layered model (a, left panel), and the recovered resistivity section from combined inversion conventional solution approach (b, right panel). White lines are the real boundaries.
ð21Þ
ð22Þ
The smaller values could cause over focusing while larger values could lead to smooth conductivity distribution. The last example employed real logging data. The logs were recorded over a coal mine in Turkey. The structural setting is as follows. The geological sequence starts with marl unit (175.00– 182.00 m). A thin clay layer (182.00–182.30 m) separates it from a first lignite layer (182.30–183.30 m). A second clay (183.30– 184.50 m) layer overlies the second lignite (184.50–186.00 m) layer and logging ends in clay which also is the bottom of the well. To resolve these units using solely direct current resistivity logs is next to impossible. In order to evaluate the modeling result required depth information is obtained from Gamma logging and mechanical soundings. There is no clear indication of the clay–lignite boundaries
136
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
Fig. 6. Observed (marked) and calculated (solid line) curves for short normal (+) and long normal (o) tools for two-blocks and a layered model (a, left panel), and the recovered resistivity section from combined inversion with relaxation solution approach (b, right panel). White lines are the real boundaries.
Fig. 8. Observed (marked) and calculated (solid line) curves for long normal tool for two-blocks and a layered model (a, left panel), and the recovered resistivity section from relaxation approach (b, right panel). White lines are the real boundaries.
in any of the resistivity log curves. On the other hand, the lower boundary of the marl unit can be discerned in both curves. Thus it may be expected that the data inversion will delineate only the marl–clay interface. The result of an attempt to recover a model that fits the data is presented in Fig. 10. In this figure, the presentation format is the same as for Fig. 5. The number of observed data is 112, i.e. 56 readings for each tool. The automatically produced initial model has 54 and 281 cells in x and z directions, respectively. The first 6 cm into the borehole wall was assumed as mud-filled and assigned a conductivity of 0.1 S/m (10 Ω m) but not excluded from the inversion. The rest of the model consists of 0.01 S/m (100 Ω m). The inversion process with relaxation approach started with a misfit value of 0.02186, decreased to 0.0005199 after 11 iterations in 136 s of CPU time (Fig. 11). The curve ‘%rate’ plots the value of Eq. (15) and presents the variation in misfit functional iteration to iteration. It appears that after the sixth iteration improvement is very low in model below 0.1%. The recovered model delineates the marl–lignite boundary in radius of
1 m (Fig. 10). The low resistivity contrast limited model resolution to a large extent. The test with the short normal tool delineates the marl boundary with the sign of thin layer (not presented) but model appeared fairly layered. In such a case, a priori information would be needed to obtain a more detailed model. One could, for example, imposed known boundaries into the model with approximate conductivities and then let the inversion process adjust the model to provide a fit to the logging data. Such information may be obtained from either other wells exist in the local or the curves themselves by means of some transformation. Next test case considers obtaining the initial model from a conductivity transform. The calculation zone and then the rest of the mesh are filled with the apparent conductivity values obtained from short- and long-normal observed data (details are given in the later part). Once again the first 6 cm into the borehole wall was assumed as mud-filled and assigned a conductivity of 0.1 S/m (10 Ω m) but not excluded from the inversion. The inversion process started with a misfit value of 0.00901, decreased to 0.00035 after 30 iterations. The inversion took 405 s of CPU time. Addition to previous test case thin clay layer
Fig. 7. Objective function versus CPU time resulting from both relaxation and conventional approach.
Fig. 9. The result for Long normal tool with conventional solution (see Fig. 6 for the explanations).
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
Fig. 10. Observed (marked) and calculated (solid line) curves for short normal (+) and long normal (o) tools (a, left panel), and the recovered resistivity section from joint inversion without any priory information (b, right panel). White lines are the marl, clay and lignite layers (see text).
recovered clearly (Fig. 12). The variation of objective function and % rate are given in Fig. 13. %rate indicates that the model updates reduce the objective function very slowly while the number of iterations increases. In fact, 1e−4 can be used as trash hold and program can be terminated. Since after twenty-second iteration no significant improvement is provided by the rest of the iterations. Approximate conductivities may be calculated through a conductivity transform even if there is no information, about the structural setting. A simple conductivity transformation, σt can be obtained via −
σt = e
lnðρa Þ +
∂ρ = ∂z maxðj∂ρ = ∂zjÞ
ð23Þ
where derivatives of apparent resistivities are calculated via derivative of second order polynomial approximation to three successive data points over the logging curves. Then, transformed conductivities obtained from short normal curve fill the first 30 cm of the mesh while transformed conductivities of long normal curves fill the cells after the
Fig. 11. Objective function and variation of objective function versus CPU time resulting from relaxation approach in inversion without any priory model.
137
Fig. 12. Result from joint inversion with priory information (see Fig. 10 for explanation).
first 70 cm of the mesh. A linear variation is assumed between the two zones.
8. Discussion The success of 2D inversion for surface surveys, as described in many papers (e.g. Candansayar, 2008b; Dahlin and Zhou, 2004), depends largely on the use of different electrode spacings to provide the necessary vertical resolution. Thus, it is expected that the use of just two spacings for the borehole survey would have much poorer resolution in the radial direction. On the other hand, this approach still provides better resolution than usage of single spacing. Test runs and knowledge from surface surveys indicate that utilization of different tool lengths can improve the success of the recovery. Different types of data, both normal and focusing, can also help to obtain better subsurface resolution. Although this paper has considered a fairly simple test model, the overall results indicate that inversion can recover the structural information of thin layers, providing there is enough resistivity contrasts between the adjacent layers. On the other hand the ratios of
Fig. 13. Objective function and variation of objective function versus CPU time resulting from relaxation approach in inversion with a priory model.
138
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
the recovered cell conductivities vary over a wide range (0.10–3.0). Recovery of true conductivity of the layers still needs additional work. The solution of Eq. (6) may be obtained using different techniques. During the development stage three solvers were used; the SVD based solver, MATLAB built-in function and the conjugate gradients. The SVD based solver demands high memory, depends on the selection of α, and could not provide a better resolution than the others in every test case considered. The MATLAB built-in-function “\”, based on the Gaussian elimination, also demands high memory since both process the full matrices. The conjugate gradients solver demands less memory than previously mentioned solvers but may fail to find the model when the misfit is small. Memory restriction was a major problem resulting in the conjugate gradients being the solver of choice in this paper. The current version of the program can handle a limited number of readings since the longer the data set the larger the model size that needs to be considered. It is certain that parallel computation can provide the necessary CPU and memory capacity for handling large size problems. Two approaches were considered in CG solver; conventional and relaxation solution. The results of the various tests indicate that neither is superior to other. The relaxation approach reduce the misfit very rapidly in very first iterations therefore if the threshold (i.e. noisy data) is high it may save time. In the mean time, the conventional solution approach reduces the misfit functional steadily even if the threshold is extremely low. The focusing parameter was defined as model depended. The percent of the minimum conductivity serves for this purpose. Besides the inversion, the adopted axisymmetric approach also has some pitfalls; non symmetric structures like dipping layers require a three-dimensional modeling approach. This will be a topic of future study. 9. Conclusion A method for undertaking combined inversion of both the shortand long-normal direct current resistivity data in two-dimensional borehole environments has been developed and tested using synthetic and real data. Joint inversion of direct current resistivity data for both short and long normal tools recovered an image of the subsurface. Tests on synthetic data indicate that each tool spans a different zone and jointly inverting multi-tools data improves the resolution of subsurface targets. Application of the inversion scheme to field data from a coal exploration survey was less successful because of the low contrast between the coal and confining clay units. It was found that the problem may only be solved if it incorporates available reliable a priori information about geological structure. Acknowledgments The author thanks M. A. Meju (PETRONAS) for reading and making several helpful suggestions that improved the clarity of the draft manuscript. The author is grateful to U. Zaman (MTA) for providing the coal data and geological information. E. Peksen (KOU) and M. E. Candansayar (AU) are also thanked for their help in programming. Thanks extend to Klaus Holliger (Editor-in-chief) and the revivers for their patients and constructive comments.
References Candansayar, M.E., 2008a. Two-dimensional inversion of magnetotelluric data with consecutive use of conjugate gradient and least-squares solution with singular value decomposition algorithms. Geophysical Prospecting 56 (1), 141–157. Candansayar, M.E., 2008b. Two-dimensional individual and joint inversion of threeand four-electrode array dc resistivity data. Journal of Geophysics and Engineering 5, 290–300. Cozzolino, K., da Conceição da Silva, J., 2007. Synthetic focusing and simulation of dual laterolog tool in axisymmetric subsurface models. Journal of Applied Geophysics 61, 102–110. Dahlin, T., Zhou, B., 2004. A numerical comparison of 2D resistivity imaging with 10 electrode arrays. Geophysical Prospecting 52, 379–398. Daniels, J.J., 1978. Interpretation of buried electrode resistivity data using the layered earth model. Geophysics 43, 988–1001. Dey, A., Morrison, H.F., 1979. Resistivity modelling for arbitrarily shaped twodimensional structures. Geophys. Prospecting 27 (1), 106–136. Dutta, D.J., 1994. TRANS4: a FORTRAN program for computing apparent resistivity departure curves for an infinitely thick bed with transitional invaded zone in borehole geophysics. Computers & Geosciences 20 (3), 293–311. deLugao, P.P., Wannamaker, P., 1996. Calculating the two-dimensional magnetotelluric Jacobian in finite elements using reciprocity. Geophysical Journal International 127, 806–810. deLugão, P.B., Portniaguine, O., Zhdanov, M.S., 1997. Fast and stable two-dimensional inversion of magnetotellurics data. J. Geomag. Geoelectr. 49, 1469–1497. Haber, E., Oldenburg, D., 1998. A GCV based method for non-linear ill-posed problems. Computational Geosciences 4, 41–63. Kwon, Y.W., Bang, H., 2000. The Finite Element Method using MATLAB, 2nd Edition. CRC Press. LaBrecque, D.J., Oldenborger, G.A., Sharpe, R., Knoll, M.D., 2006. Axi-symetric inversion of electrical resistivity tomography data to monitor the movement of fluids injected in wells. Symposium on the Application of geophysics to Environmental and engineering problems, p. 1505-1503. Liu, Z., Lin, H., 2002. Joint inversion of induction/lateral/normal logs, case studies at Shenli field site, China. Journal of Petroleum Science and Engineering 34, 55–64. Mackie, R.L., Madden, T.R., 1993. Three dimensional magnetotelluric inversion using conjugate gradients. Geophysical Journal International 115, 215–229. Mufti, I.R., 1976. Finite-difference resistivity modeling for arbitrarily shaped twodimensional structures. Geophysics 41, 62–78. Pardo, D., Paszynski, M., Torres-verdin, C., Demkowicz, L., 2007. Simulation of 3D DC borehole resistivity measurements with a goal-oriented hp finite-element method, Part I: laterolog and LWD. Journals of the Serbian Society for Computational Mechanics 1, 62–73. Pidlisecky, A., Haber, E., Knight, R., 2007. RESINVM3D: a 3D resistivity inversion package. Geophysics 72, H1–H10. Portniaguine, O., Zhdanov, M.S., 1999. Focusing geophysical inversion images. Geophysics 64, 874–887. Rodi, W.L., 1976. A technique for improving the accuracy of finite element solutions for magnetotelluric data. Geophysical Journal of the Royal Astronomical. Society 44, 483–506. Rodi, W.L., Mackie, R.L., 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics 66, 174–187. Roy, K.K., Dutta, D.J., 1997. Relative performance of focussed and unfocussed devices in a two-dimensional borehole environment. Journal of Applied Geophysics 36, 181–194. Tikhonov, A.N., Arsenin, V.Y., 1977. Solutions of Ill-Posed Problems. V.H. Winston and Sons, Washington, D.C. Towle, G.H., Whitman, W.W., Kim, J.H., 1988. Electric log modeling with a finite difference method. The Log Analyst 184–195 (May–June). Yang, F.W., Ward, S.H., 1984. Inversion of borehole normal resistivity logs. Geophysics 49, 1541–1548. Zhang, G.J., Shen, L.C., 1984. Response of a normal resistivity tool in a borehole crossing a bed boundary. Geophysics 49 (2), 142–149. Zhang, J., Mackie, R.L., Madden, T.R., 1995. 3-D resistivity forward modeling and inversion using conjugate gradients. Geophysics 60, 1313–1325. Zhdanov, M.S., 2002. Geophysical Inverse Theory and Regularization Problems. Elsevier Science Pub\.Co., Inc. Zhdanov, M.S., Tolstaya, E., 2004. Minimum support nonlinear parameterization in the solution of 3D magnetotelluric inverse problem. Inverse Problems 20, 937–952.
Contents lists available at ScienceDirect
Journal of Applied Geophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j a p p g e o
Two dimensional combined inversion of short- and long-normal dc resistivity well log data Emin U. Ulugergerli ⁎ Çanakkale Onsekiz Mart University, Faculty of Engineering and Architecture, Department of Geophysical Engineering, Çanakkale, Turkey
a r t i c l e
i n f o
Article history: Received 19 October 2009 Accepted 8 December 2010 Available online 23 December 2010 Keywords: Well-logging Direct current resistivity 2D inversion Finite element method
a b s t r a c t A method is presented for the two-dimensional combined inversion of short- and long-normal tool direct current resistivity data with symmetry. The forward problem is solved using the finite element method in the cylindrical coordinates system. The inverse problem is solved using a conjugate gradient technique with the partial derivatives obtained using reciprocity. The parameters were obtained by means of both conjugate gradient relaxation and conventional conjugate gradient method. The solution of this highly underdetermined inverse problem is stabilized using Tikhonov regularization and the scheme yields a blurred image of the subsurface. The scheme is tested using synthetic data and field data. Tests using synthetic data suggest that traces of the horizontal boundaries are delineated in the range of the exploration distance while the resolution of vertical boundaries depends upon the solution regularization. Application to field data shows that additional information is necessary for resolving the resistivity structure when there are low resistivity contrasts between formation units. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Well logging applications produce useful information about geological formations where the thicknesses of formations are larger than the tool's size. Yang and Ward (1984) generalized the investigation results of Daniels (1978) to show the relation between electrodes spacing and layer thickness; if the layer is thicker than ten times the electrode spacing, the apparent resistivity will be closer to the true layer resistivity. If the thickness is less than six times the electrode spacing, estimating true resistivity will be cumbersome and departure curves (e.g. Dutta, 1994) may be employed to estimate the true resistivity values. Departure curves are similar to normalized apparent resistivity curves in surface application of the direct current resistivity (DCR) method. In reality, features such as borehole fillings, invaded zone, layering, and shoulder beds create a complex environment. The apparent resistivity curves represent filtered and simplified images of geo-electrical structures. Multi-dimensional modeling is a convenient way to image complex subsurface structures. As well known for surface measurements, two-dimensional (2D) modeling of the DCR method requires the solution of a boundary value problem in Fourier transformed space (e.g. Dey and Morrison, 1979). However, this requirement is not necessary when radial symmetry is assumed. The axisymmetric subsurface model indicates that the well
⁎ Tel./fax: +90 286 2180018x2207, +90 286 2180541. E-mail address: [email protected]. 0926-9851/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2010.12.004
itself together with the invasion zones and the resistive structures have radial symmetry in the range of exploration and the resistivity variations may occur in both vertical and radial directions. This condition assumes that the structures have concentric ring-like shapes, centered on the well axis. Axisymmetric layout of a well allows one to employ a cylindrical coordinate system. In such a case, either analytic solutions or numerical techniques can be used in the evaluation of electric tool responses. If the concentric ring-like resistivity structures have infinite thickness, the problem is reduced to one-dimension (1D) and an analytical solution will provide the required solution (e.g. Daniels, 1978; Yang and Ward, 1984). If the thicknesses are finite then, either the finite element method (FEM, e.g. Cozzolino and da Conceição da Silva, 2007), the finite differences method (FDM, e.g. Roy and Dutta, 1997; Towle et al., 1988), the integral equation method (IEM, e.g. Zhang and Shen, 1984) or a combination of all of these methods can be employed. In this paper, the FEM with bilinear rectangular element has been used for problem discretization and to solve the axisymmetric earth model. The modeling scheme takes into account only one half of an infinite number of planes in the domain. The use of non-uniform grid spacing allows for modeling large scale structures. The inversion of DCR logging data for subsurface resistivity distribution is a highly under-determined problem. While only the surveyed zone is considered in the inversion process, equations for hundreds of model parameters are solved per DCR reading. Hence, the problem requires regularization. The resulting regularized problem is then solved using the conjugate gradient (CG) method. Two conventional approaches are considered in this paper. First one is
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
CG relaxation (e.g. Mackie and Madden, 1993; Zhang et al., 1995) and will be called as “relaxation” approach in the text. Second solution is yielded via conventional CG (e.g. Zhdanov, 2002; Zhdanov and Tolstaya, 2004) and will be called as “conventional” approach in the text. The definitions are used solely for the purpose of abbreviations. For computational efficiency, the principle of reciprocity is used for the calculation of partial derivatives. Calculating the partial derivatives of potentials with respect to the model parameters for each reading point requires thousands of forward solutions for each inversion iteration. The reciprocity method given by deLugao and Wannamaker (1996) reduces this necessity to performing two forward solutions per reading; one for current sources and another for a unit source at the receiver site. As mentioned previously, simulation of DC logs in axiallysymmetric media was studied extensively. Many of the publications usually focus on departure curves. Some publications related to inversion may be found either in proceedings (e.g. LaBrecque et al., 2006) or special journals (e.g. Pardo et al., 2007). Unlike the induction logs (e.g. Liu and Lin, 2002), interpretation stage of DC logging usually consists of plotting and evaluating the apparent resistivity curves with using departure curves or performing single-log inversion only. The innovative part of the paper is the combined inversion of short and long normal DC data. This approach, then, tested on low-contrast case which is the bane of inversion. The resulting inversion schemes will be evaluated using synthetic and field data. Observable data for short normal (16 in.) tools and long normal (64 in.) tools are calculated for simple but representative borehole environments. The field data considered are from a low-contrast geoelectrical situation. Both data sets are used in the inversion routine. Tests on synthetic data indicated that the single log inversion of short normal data recovers structural pattern in the range of the exploration radius. On the other hand, the inversion of long normal data delineates only the anomalous zone. Besides the single log inversion, it is also demonstrated that joint inversion of short- and long-normal data produces a blurred picture of the subsurface, and that while horizontal boundaries can be delineated in the vicinity of the well, the recovery of structures with low conductivity contrast is a hard task requiring some a priori information. 2. Solution of the forward problem The governing equation for the flow of a steady state electric current in a non-uniform medium (Dey and Morrison, 1979; Mufti, 1976) is given by ∇ ðσ ðx; y; zÞ∇νðx; y; zÞÞ =
∂Q ðx; y; zÞ ; ∂t
ð1Þ
where σ is the conductivity, υ is the electrical potential and Q is the charge density. For the axisymmetric case, Eq. (1) can be rewritten as ∂ ∂νðr; zÞ ∂ ∂νðr; zÞ 1 ∂νðr; zÞ σ ðr; zÞ + σ ðr; zÞ + σ ðr; zÞ r ∂r ∂r ∂z ∂z ∂r =
131
measurement zone of the mesh has an equal spacing in the vertical (z) direction up to the end of the surveyed zone followed by a gradually increase of cell thicknesses. Similarly, the width of the cells gradually increases in the radial (r) direction. The presented test models used 2 in. (≈5 cm) spacing in the vertical direction which mimic borehole logging surveys with 2 in. intervals. The radial intervals vary between 0.125 and 100 in. Application of the FEM to Eq. (2) produces a set of equations, ð3Þ
Gν = Si
where G is the global or stiffness matrix, v is the vector of potentials and Si is the current source vector for the ith measurement. Eq. (3) is solved for each source location. These computations were implemented in MATLAB. Since the program was developed in MATLAB, a builtin-function of back slash (\), based on Gaussian elimination with partial pivoting is used to solve Eq. (3). The code uses four electrodes in the computation to obtain a complete simulation. Therefore, the possible effect of the shoulder bed or upper layers may also be considered. The distance between N and M electrodes is set at 15.24 m (50 ft) and B electrode placed at the surface. The geometric factor is set alike to the conventional four electrodes system to handle the deviation in the near surface measurements, k=
4π 1 1 1 AM − AN − BN
+
1 BM
ð4Þ
Due to the fact that all the distances are very large compared to the AM distance (Fig. 1), the operative distance is AM and is set to either 0.4 m (≈16 in.) or 1.6 m (≈64 in.), representing the short and long normal tools, respectively. Once the potential and geometric factors are calculated, the apparent resistivities are obtained as ρa = k
∇VMN : I
ð5Þ
Fig. 2a compares the forward solution from the present code with the analytical result of Dutta (1994). The 1D model considers mud and two zones; the mud, flushed and uncontaminated zones which have resistivities of 1, 10 and 0.5 Ω m, respectively. The first two zones are 1 and 1.5 m thick, respectively. A 6.5 m thick invaded zone has linear variation with radial distance. Details of the analytical calculation can be seen in Dutta (1994). Additionally, Fig. 2b presents the apparent resistivity curves of short-normal logging for the case where there is no borehole. The target layer is 80 in. (≈2 m) thick and has a resistivity of 10 or 100 Ω m in models 1 and 2 respectively. Models 1 and 2 have background resistivities of 100 and 10 Ω m, respectively. In both cases the borehole radius is 8 in. (≈0.2 m). Model 2 also has a
ð2Þ
Io δðr Þδðz−zc Þ 2π r
where r is the radial distance from well and z defines the depth. r = 0 is the centre of the well while z = 0 is the surface. Current, Io is injected at depth zc in the well (r = 0). v is the potential created by the current. FEM with bilinear rectangular elements with four nodes (e.g. Kwon and Bang, 2000) is used to solve Eq. (2). Various types of boundary conditions can be applied to edges of the mesh. Since the finite element method uses all nodes in calculations, I did not impose any external values to boundaries but set all borders away from the region of interest. Hence, the
Fig. 1. Basic model for borehole (not scaled). Right side is the calculation domain. Rings on toll indicate the electrodes.
132
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
optimization problem. The problem is underdetermined and needs to be regularized. The objective function to minimize therefore incorporates stabilizing functional (Tikhonov and Arsenin, 1977), viz: F = ϕðpÞ + α sðpÞ = min
Fig. 2. (a) Comparison of analytical (Dutta, 1994) and numerical (FEM) solution of Eq. (2) for 1D model. Resistivity of invasion zone obeys linear variation with radial distance. Axes are resistivity (Ω m) and transmitter–receiver spacing (AM). (b) Sample curves for the case with a 80″ (~ 2 m) thick layer. Both cases the borehole radius is 8″ (~0.2 m). Model 2 also has a transition zone of 50 Ω m, radius of which is 40″ (~1 m).
transition zone of 50 Ω m, the radius of which is 40 in. (≈1 m). Since the ratio of layer thickness to tool length (80/16) is smaller than 10 in Fig. 2b, the apparent resistivity curve deviates from the true value. An inversion technique will be employed to obtain the formation resistivity value in such a case instead of using the traditional correction charts.
3. Methodology for 2D joint inversion of short- and long-normal tool responses The inverse problem is to recover the conductivity distribution around the well from the typically few measured dc resistivity welllogging data. Since the mesh defines the electrical resistivity structure permeated by the well, the logarithms of conductivities of each cell in the 2D mesh are taken as the model parameters (p) while the logarithms of apparent resistivities serve as data (d) in our nonlinear
ð6Þ
where ϕ(p) is a measure of fit between observed and calculated data. α is a regularization parameter which provides a balance between ϕ(p) and stabilizing functional s(p). The solution to Eq. (6) is subject to various publications. The outlines of the researches focuses on that the selection of regularization functional, s(p) and that the investigation of algorithms to compute regularized solution. Starting from the first one, a wide range of stabilizing functional is available in the literature (e.g. Zhdanov, 2002). Depending on the nature of the inverse problem taken into consideration, different stabilizer can be employed in Eq. (6). Main expectations from a stabilizer functional can be summarized as follows (e.g. Zhdanov, 2002; Zhdanov and Tolstaya, 2004); it should (a) stabilize the illposed inverse problem, (b) reduce the number of possible models responds of which fit the data, (c) generate a compact image of the geo-electrical model and (d) not smear out the small features. Equivalent expressions in practice are all fallows, in all multidimensional modeling studies, the lack of enough information content in data, which leads to solution to an underdetermined problem (statement a), urges a priory model to be used. Then, the Tikhonov regularization theory seeks a model close to that priory (or reference) model (statement b). The most cases a reference model cannot be enough to get an appropriate solution. Thus some additional constraints can be provided by, for instance, the minimum norm of the second order finite difference operator, also called maximum smoothness stabilizing, through which neighbor parameters are forced to be correlated with each other. This constraint results in a model presents smooth resistivity distribution. As an alternative to smooth variation, focusing or minimum support (MS) or minimum gradient support stabilizing functionals could be considered where detailed structural information (statements c and d) is required (e.g. Portniaguine and Zhdanov, 1999). In our case, the sharp boundaries are needed to describe the complex resistivity distribution. Among a large family stabilizers, MS were proved to be effective to recover a model consists of discontinuous images (Zhdanov and Tolstaya, 2004). The stabilizing functional, s(p) is given as 2 sðpÞ = ‖Wm Ws p−pref ‖
ð7Þ
where ||.|| represent the two-norm, Wm is a model weighting matrix, diagonal matrix Ws is the stabilizing operator. p and pref are the logarithms of the current model and the reference model, respectively. The weighting for MS stabilizer, 0
1
B C 1 B C Ws = diagB 1 = 2 C; 2 @ A 2 +ε pi−1 −pref
ð8Þ
where ε is the focusing parameter. It is calculated, at each iteration, by dividing the minimum conductivity value in the mesh by a hundred. Note that Ws depends on the parameters from previous iteration. Once decided the stabilizer, the second step requires the selection of inversion algorithm. Rodi and Mackie (2001) compared nonlinear conjugate gradient (NLCG), Gauss-Newton (GN) and linear conjugate gradient (LCG) algorithms to compute regularized solution of the 2D magnetotellurics inverse problem. They concluded that the conjugate gradients (CG) based algorithms are superior to the GN algorithm. Later, Candansayar (2008a) included the comparison of different stabilizing functionals to the same inverse problem. He stated
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
weakness and superiority of the algorithms in terms of stabilizing functionals. Both studies declare that the CG based algorithms is favorable for the multi-dimensional inverse problems. Before going further, first term in Eq. (6) is also misfit functional is given as 2
ϕðpÞ = ‖Wd ðd−F ðpÞÞ‖
ð9Þ
where F(p) is the forward solution given by Eqs. (3) and (5). Diagonal matrix, Wd consists of positive weightings which are a function of the measure of the uncertainties of the observed data. Since we do not have any information on data errors, Wd set identity and is dropped from Eq. (9) for simplicity. Linearization of F(p) yields (e.g. Haber and Oldenburg, 1998), 2
ϕðpÞ = ‖ J∇p−Δd‖ ;
ð10Þ
where J is the sensitivity matrix, Δp are the model updates, the vector Δd is the discrepancy between the logarithm of observed and logarithm of calculated apparent resistivities. For the joint inversion problem, the partitioned matrix, J = ½ JSN JLN
T
ð11aÞ
and the partitioned vector, Δd = ½ΔdSN ΔdLN
T
ð11bÞ
represent the short-normal and long-normal contributions to the system of equations. There are many procedures to obtain the solution to Eq. (6) but the CG approach (e.g. Pidlisecky et al., 2007; Zhang et al., 1995; Zhdanov, 2002) is adopted in this study. Basic idea is that one needs to find Δp which, in turn, minimize the Eq. (6). The model updates are sought either in one step (e.g. Zhdanov, 2002) or by means of relaxation at each iteration of the inversion (e.g. Mackie and Madden, 1993). These procedures, in this study, will be named as the conventional and the relaxation solution, respectively and details of each will be given later in the text. In any case the solution to Eq. (6) yields perturbation (or updates) for parameter pi. pi +1 = pi + b Δp
ð12Þ
where b is a line-search parameter (e.g. Pidlisecky et al., 2007) which controls the magnitude of the update vector. Besides the scaling by b, the Δp values are truncated as:
These bounds are selected arbitrarily and prevent large steps in the model update process. Finally, the model weighting matrix Wm is defined as (e.g. Portniaguine and Zhdanov, 1999), k T Wm = diag J J
inversion iteration and maximum number of relaxation iteration, respectively. The empty lines at right panel indicate that the corresponding lines at left panel are excluded. Relaxation solution
Conventional solution
For k = 1 to ninv J calculation and Update W Weight J r = d − F(pk − 1) Io = JTr − αWTW(pk − 1 − pref) r0 = Io; dp0 = 0; b = 1 For j = 1 to line search rk = JTr − αWTW(pk − 1 − pref) For i = 1 to nrex beta = rTi ri/rTi − 1ri − 1 Ii = r i + beta Ii − 1 step = rTi ri/ITi [JTJ + αW] Ii dpi = dpi−1 + step I1 ri = ri − 1 − step L End loop i pk = pk − 1 + b dpk fk = || d − F(pk)||2/N If fk b fk − 1 ; break loop j Else : b = b/1.5; α = α/2 End loop j If fk b treshold ;break loop k End loop k
For k = 1 to ninv J calculation and Update W Weight J r = F(pk − 1) − d Io = JTr + α W (pk − 1 − pref) r0 = Io; b = 1 For j = 1 to line search rk = JTr + α W (pj − 1 − pref)
ð13Þ
where subscript k controls the weighting level. Since the inverse of the diag(JTJ) affects the partial derivatives, k values controls the weighting level. If k N 1, the parameters distant to the well are overweighted. On the other hand, if k → 0, those parameters closer to the well are over-weighted. In this study k was set at 0.5. 4. Conjugate gradient solution of the inverse problem As stated before, two different approaches were considered to carry out the solution to Eq. (6). The respective algorithms given in the later part are modified from Mackie and Madden (1993) and Zhdanov (2002). The ninv and the nrex define maximum number of
beta = rTkrk/rkT − 1rk − 1 Ik = rk + beta Ik step = rTkrk/ITk[JTJ + α I] Ik
pk = pk − 1 − b step Ik fk = || d − F(pk)||2/N If fk b fk − 1 ; break loop j Else b = b/1.5; α = α/2 End loop j If fk b treshold ; break loop k End loop k
Innermost loop can end when they reach their predefined number of iteration (nrex). Outer loop can end either when they reach their predefined number of iterations (ninv) or the misfit reaches to predefined threshold value or insignificant improvement in the model. The data fit, fk is measured by the quantity T
misfit = fk = ΔG ΔG = N
ð14Þ
where N is number of the observations. ΔG is a difference vector of logarithm of the observed and calculated data. Being a stopping criterion, target misfit value needs to be defined by the user. Unfortunately there is no way to select appropriate value since our problem is severely underdetermined and a variety models can provide the same fit to the data being inverted. Any model with a small misfit value (b1.E−06) may be accepted as satisfactorily explaining the inverted data. Insignificant improvement in the model means that the variation of successive misfit values, dmisfit =
if Δp N 2.0 then Δp = 2.0 if Δp b − 2.0 then Δp = −2.0.
133
misfitk−1 −misfitk misfitk
ð15Þ
is less than 1.E−5. Any priori information about any part of the resistivity model structure can reduce the number of iterations and move the result towards the expected model. Because the finer cell widths are used near the well, to recover their conductivities from inversion will be very difficult, thus, the borehole itself is defined with its mud conductivity. Knowledge of the mud conductivity can be incorporated into the initial model and then the left-most cells as priory information. 5. Regularization parameter The solution of Eqs. (12) and (13) requires the determination of regularization and line-search parameters. Each of them needs to be sought during the iterative process, the initial value for the regularization parameter is given, αo = ‖ JΔp−Δd‖ = ‖Wm po −pref ‖:
ð16Þ
If it does not reduce the misfit it is reduced by two, αn = αo/2.
134
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
A line-search is employed in relaxation solution only and performed in a simple manner as given in Pidlisecky et al. (2007). Initially, the value of b is set to unity and subsequently decreased using the relation bi = bi−1 = 1:5:
ð17Þ
6. Sensitivity matrix Calculation of the sensitivity matrix. J requires the derivative of Eq. (3), which, in the general case, is (Rodi, 1976), G
∂ν ∂G =− ν: ∂σ ∂σ
ð18Þ
The right-hand side (RHS) in the Eq. (18) is a new source term for the forward solution. Without employing the principle of reciprocity the most time-consuming part of the inversion process would be the calculation of partial derivatives. The method given by deLugao and Wannamaker (1996) reduces the calculations to two forward solutions per reading; one for current sources and another one for unit source. deLugão et al. (1997) define the forward solution with unit source as, GR = U1
ð19Þ
where U1 is a vector which has a value of unity at the receiver location. The solution vector, R contains the model responses in all nodes. Multiplying the appropriate responses by the cell parameters gives the derivative values. In the adopted FEM, the cell parameters are the element matrix obtained from the shape functions of each node. The derivatives are multiplied by the ratio of apparent resistivity to conductivity of relevant cell to obtain a logarithmic quantity. Fig. 3a and b shows the sensitivity zone for both the short- and long-normal tools. The receiver electrode (M) is at 862 in. (≈9.7 m) while a current electrode (A) is at 878 in. (≈22.3 m) and 926 in. (≈ 23.5 m) for short- and long-normal tools, respectively. The derivatives are normalized by their maximum value and presented for a selected receiver–transmitter position. The absolute values were used in the graphs. Selecting 0.03 arbitrarily, the maps indicate that 16 in. toll can comb almost twice the tool length while the 64-inch tool can only comb approximately 1.5 times the tool length. Note that this limit does not apply to the investigation distance (δ) but it should also adhere to a similar range ratio. One may stipulate this condition as, δ16 δ ≥ 64 : 16 64
ð20Þ
If one defines the investigation distance as a radius that data begins to sense the outer unit then, two-unit model can be used to examine the radius of investigation (ROI). Fig. 4 presents short (+) and long (triangular) normal curves for two-unit model. Horizontal axis is the radius of the interface between first and second unit. The number indicates the resistivity ratio of inner unit to outer unit (ρ1/ ρ2). Low contrast case (1/0.1 to 1/10) 64 in. toll sense ~ 24 in. (~60 cm) radius while 16 in. tool reaches to ~6 in. (~ 16 cm). Behavior of the curves becomes cumbersome in case of high contrast but ROI remains same. The results are in accord with the Eq. (20). 7. Application of the inversion algorithm to synthetic and field data The recovered geo-electrical sections from the joint inversion of short- and long-normal logging data are presented here for both the synthetic example and for real data. Although the Eq. (2) depends
Fig. 3. Normalized sensitivity map for short-normal tool (left) and long-normal tool (right). The resistivity of the medium is 1 Ω m.
upon the variation of conductivities, the sections are presented as resistivities. The first example considers a test model which is given in Fig. 5. The model considers the thin layer problem. The anomalous section has two blocky units of 0.0033 S/m (300 Ω m) and a layer of 0.033 S/m (30 Ω m) immersed in a 0.01 S/m (100 Ω m) host medium. The thickness of the each layer is 100 cm (~39.5 in.) while the width of each resistive block is 100 cm. The conductive layer (mid unit) extends to infinity. The ratio of layer thickness to tool length is ~2.4 for short-normal tool but ~0.6 for long-normal tool. Such combinations result in deviations in apparent resistivity curves. The long-normal curve has very little signatures for the layers while short normal curve represents all units. The well has 20 cm diameter and is filled with a mud of 10 Ω m (0.1 S/m) resistivity. The depth position of the first reading is 178.5 m in the test models. The number of data is 142, i.e. 71 readings for each tool. The noise adding is tricky in single reading data set. Short wave length variations (i.e. random noise adding to each reading) cause a recovery of artificial features in the vicinity of the borehole. Thus a trend is added to both data set. Minima and maximum values of the linear trend is 0.1 and 7.1, respectively. Once synthetic data were created, they were fed to the program which creates the initial model automatically depending on the data and survey parameters. The
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
135
are plotted versus location of their current electrodes (A in Fig. 1). The right panel is the recovered model. The vertical scale is the same as the left panel and white lines define the real boundaries of the sought model. The recovered model presents an expected result at the border of the radial exploration range (~ 1 m); thereafter, the image of model deteriorates in terms of the cell conductivities. Both resistive and conductive units appear over focused since the algorithm uses fixed value for focusing parameter. The over focusing caused separation in the horizontal interfaces of layers. Additionally, the inversion could not provide the exact conductivity values for each layer. If one defines the ratio for real and calculated and conductivities as r=
Fig. 4. Short (+) and long (triangular) normal curves for two-unit model. The number indicates the resistivity ratio of inner unit to outer unit (ρ1/ρ2). Horizontal axis is the radius of the interface between first and second unit.
initial model has 54 and 296 cells in x and z directions, respectively. The conductivity of leftmost 20 cm of the model is assumed to be known as priory information and set 0.1 S/m (10 Ω m) but not excluded from inversion. Only the anomalous zone of 3150 cells in total, between the first and last reading points, was considered in the inversion and rest of the model was kept fixed and excluded from inversion. After each iteration, the conductivity values at the border of anomalous zone are extended to mesh borders not to create any discontinuity between shoulder and calculation zones. For the CG routines the focusing parameter was fixed to 0.003 (Zhdanov and Tolstaya, 2004). The misfit threshold set 1.e−5. First test was run for the relaxation approach. The inversion of the test model data started with misfit value of 0.0063357 and stopped after 12 iterations. The misfit for the result was 6.3.e−6 . All the processes, including plotting for inter-steps took 192.35 s CPU time on the Intel Core2 @ 3 GHz processor. The left panel in Fig. 5 presents the observed (marked symbol) and calculated (solid line) resistivity curves for short normal (+) and long normal (o) tools. Both curves
ðσreal Þ ðσcalc Þ
it varies from 0.2 to 3.4 for the anomalous zone. Second test was run with the conventional solution approach. The inversion stopped after 10 iterations. The misfit for the result was 6.3. e−6. All the processes, including plotting for inter-steps took 164 s. Note that both approaches end the iterations with the misfit reaching to threshold value. This case, the ratio for real and calculated conductivities varies from 0.2 to 2.8 for the anomalous zone. Same as before the horizontal interfaces of layers are determined. Recovery is slightly better than previous case (Fig. 6). Variation of the objective functional (Eq. (6)) vs. CPU time is given in Fig. 7. The relaxation approach decreases objective functional very rapidly in first few iterations but then it appears that it slows down. On the other hand, the conventional solution decreases it steadily and the misfit functional (Eq. (9)) reaches to the threshold value earlier than the relaxation approach. Third tests were run on single log data. The long normal log presents smoothed curve, but the short normal data present detailed structural features, therefore, the inversion with the short normal data set results in better recovery than the long normal data (not presented). Thus, the long normal data were selected for testing the programs. Conventional approach end in 14 iterations and reduce the misfit from 0.0059 to 9.5574e−6 in 116.5 s of CPU time. The relaxation approach, on the other hand, did 8 iterations and reached the misfit value of 9.52e−6 in 69 s of CPU time. Both approaches could not recover the structural pattern, but only delineated the anomalous zone (Figs. 8 and 9). The visual inspection of both results indicates that relaxation approach recovers slightly better structural information since the trace of the conductive layer is clearer. Overall result of these tests points outs that the relaxation approach produces slightly better result than the conventional solution approach does and that the focusing parameter needs to be obtained through the iterations. One suggestion is to use percent of the minimum conductivity value in the model. ε = minðσ Þ = 100
Fig. 5. Observed (marked) and calculated (solid line) curves for short normal (+) and long normal (o) tools for two-blocks and a layered model (a, left panel), and the recovered resistivity section from combined inversion conventional solution approach (b, right panel). White lines are the real boundaries.
ð21Þ
ð22Þ
The smaller values could cause over focusing while larger values could lead to smooth conductivity distribution. The last example employed real logging data. The logs were recorded over a coal mine in Turkey. The structural setting is as follows. The geological sequence starts with marl unit (175.00– 182.00 m). A thin clay layer (182.00–182.30 m) separates it from a first lignite layer (182.30–183.30 m). A second clay (183.30– 184.50 m) layer overlies the second lignite (184.50–186.00 m) layer and logging ends in clay which also is the bottom of the well. To resolve these units using solely direct current resistivity logs is next to impossible. In order to evaluate the modeling result required depth information is obtained from Gamma logging and mechanical soundings. There is no clear indication of the clay–lignite boundaries
136
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
Fig. 6. Observed (marked) and calculated (solid line) curves for short normal (+) and long normal (o) tools for two-blocks and a layered model (a, left panel), and the recovered resistivity section from combined inversion with relaxation solution approach (b, right panel). White lines are the real boundaries.
Fig. 8. Observed (marked) and calculated (solid line) curves for long normal tool for two-blocks and a layered model (a, left panel), and the recovered resistivity section from relaxation approach (b, right panel). White lines are the real boundaries.
in any of the resistivity log curves. On the other hand, the lower boundary of the marl unit can be discerned in both curves. Thus it may be expected that the data inversion will delineate only the marl–clay interface. The result of an attempt to recover a model that fits the data is presented in Fig. 10. In this figure, the presentation format is the same as for Fig. 5. The number of observed data is 112, i.e. 56 readings for each tool. The automatically produced initial model has 54 and 281 cells in x and z directions, respectively. The first 6 cm into the borehole wall was assumed as mud-filled and assigned a conductivity of 0.1 S/m (10 Ω m) but not excluded from the inversion. The rest of the model consists of 0.01 S/m (100 Ω m). The inversion process with relaxation approach started with a misfit value of 0.02186, decreased to 0.0005199 after 11 iterations in 136 s of CPU time (Fig. 11). The curve ‘%rate’ plots the value of Eq. (15) and presents the variation in misfit functional iteration to iteration. It appears that after the sixth iteration improvement is very low in model below 0.1%. The recovered model delineates the marl–lignite boundary in radius of
1 m (Fig. 10). The low resistivity contrast limited model resolution to a large extent. The test with the short normal tool delineates the marl boundary with the sign of thin layer (not presented) but model appeared fairly layered. In such a case, a priori information would be needed to obtain a more detailed model. One could, for example, imposed known boundaries into the model with approximate conductivities and then let the inversion process adjust the model to provide a fit to the logging data. Such information may be obtained from either other wells exist in the local or the curves themselves by means of some transformation. Next test case considers obtaining the initial model from a conductivity transform. The calculation zone and then the rest of the mesh are filled with the apparent conductivity values obtained from short- and long-normal observed data (details are given in the later part). Once again the first 6 cm into the borehole wall was assumed as mud-filled and assigned a conductivity of 0.1 S/m (10 Ω m) but not excluded from the inversion. The inversion process started with a misfit value of 0.00901, decreased to 0.00035 after 30 iterations. The inversion took 405 s of CPU time. Addition to previous test case thin clay layer
Fig. 7. Objective function versus CPU time resulting from both relaxation and conventional approach.
Fig. 9. The result for Long normal tool with conventional solution (see Fig. 6 for the explanations).
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
Fig. 10. Observed (marked) and calculated (solid line) curves for short normal (+) and long normal (o) tools (a, left panel), and the recovered resistivity section from joint inversion without any priory information (b, right panel). White lines are the marl, clay and lignite layers (see text).
recovered clearly (Fig. 12). The variation of objective function and % rate are given in Fig. 13. %rate indicates that the model updates reduce the objective function very slowly while the number of iterations increases. In fact, 1e−4 can be used as trash hold and program can be terminated. Since after twenty-second iteration no significant improvement is provided by the rest of the iterations. Approximate conductivities may be calculated through a conductivity transform even if there is no information, about the structural setting. A simple conductivity transformation, σt can be obtained via −
σt = e
lnðρa Þ +
∂ρ = ∂z maxðj∂ρ = ∂zjÞ
ð23Þ
where derivatives of apparent resistivities are calculated via derivative of second order polynomial approximation to three successive data points over the logging curves. Then, transformed conductivities obtained from short normal curve fill the first 30 cm of the mesh while transformed conductivities of long normal curves fill the cells after the
Fig. 11. Objective function and variation of objective function versus CPU time resulting from relaxation approach in inversion without any priory model.
137
Fig. 12. Result from joint inversion with priory information (see Fig. 10 for explanation).
first 70 cm of the mesh. A linear variation is assumed between the two zones.
8. Discussion The success of 2D inversion for surface surveys, as described in many papers (e.g. Candansayar, 2008b; Dahlin and Zhou, 2004), depends largely on the use of different electrode spacings to provide the necessary vertical resolution. Thus, it is expected that the use of just two spacings for the borehole survey would have much poorer resolution in the radial direction. On the other hand, this approach still provides better resolution than usage of single spacing. Test runs and knowledge from surface surveys indicate that utilization of different tool lengths can improve the success of the recovery. Different types of data, both normal and focusing, can also help to obtain better subsurface resolution. Although this paper has considered a fairly simple test model, the overall results indicate that inversion can recover the structural information of thin layers, providing there is enough resistivity contrasts between the adjacent layers. On the other hand the ratios of
Fig. 13. Objective function and variation of objective function versus CPU time resulting from relaxation approach in inversion with a priory model.
138
E.U. Ulugergerli / Journal of Applied Geophysics 73 (2011) 130–138
the recovered cell conductivities vary over a wide range (0.10–3.0). Recovery of true conductivity of the layers still needs additional work. The solution of Eq. (6) may be obtained using different techniques. During the development stage three solvers were used; the SVD based solver, MATLAB built-in function and the conjugate gradients. The SVD based solver demands high memory, depends on the selection of α, and could not provide a better resolution than the others in every test case considered. The MATLAB built-in-function “\”, based on the Gaussian elimination, also demands high memory since both process the full matrices. The conjugate gradients solver demands less memory than previously mentioned solvers but may fail to find the model when the misfit is small. Memory restriction was a major problem resulting in the conjugate gradients being the solver of choice in this paper. The current version of the program can handle a limited number of readings since the longer the data set the larger the model size that needs to be considered. It is certain that parallel computation can provide the necessary CPU and memory capacity for handling large size problems. Two approaches were considered in CG solver; conventional and relaxation solution. The results of the various tests indicate that neither is superior to other. The relaxation approach reduce the misfit very rapidly in very first iterations therefore if the threshold (i.e. noisy data) is high it may save time. In the mean time, the conventional solution approach reduces the misfit functional steadily even if the threshold is extremely low. The focusing parameter was defined as model depended. The percent of the minimum conductivity serves for this purpose. Besides the inversion, the adopted axisymmetric approach also has some pitfalls; non symmetric structures like dipping layers require a three-dimensional modeling approach. This will be a topic of future study. 9. Conclusion A method for undertaking combined inversion of both the shortand long-normal direct current resistivity data in two-dimensional borehole environments has been developed and tested using synthetic and real data. Joint inversion of direct current resistivity data for both short and long normal tools recovered an image of the subsurface. Tests on synthetic data indicate that each tool spans a different zone and jointly inverting multi-tools data improves the resolution of subsurface targets. Application of the inversion scheme to field data from a coal exploration survey was less successful because of the low contrast between the coal and confining clay units. It was found that the problem may only be solved if it incorporates available reliable a priori information about geological structure. Acknowledgments The author thanks M. A. Meju (PETRONAS) for reading and making several helpful suggestions that improved the clarity of the draft manuscript. The author is grateful to U. Zaman (MTA) for providing the coal data and geological information. E. Peksen (KOU) and M. E. Candansayar (AU) are also thanked for their help in programming. Thanks extend to Klaus Holliger (Editor-in-chief) and the revivers for their patients and constructive comments.
References Candansayar, M.E., 2008a. Two-dimensional inversion of magnetotelluric data with consecutive use of conjugate gradient and least-squares solution with singular value decomposition algorithms. Geophysical Prospecting 56 (1), 141–157. Candansayar, M.E., 2008b. Two-dimensional individual and joint inversion of threeand four-electrode array dc resistivity data. Journal of Geophysics and Engineering 5, 290–300. Cozzolino, K., da Conceição da Silva, J., 2007. Synthetic focusing and simulation of dual laterolog tool in axisymmetric subsurface models. Journal of Applied Geophysics 61, 102–110. Dahlin, T., Zhou, B., 2004. A numerical comparison of 2D resistivity imaging with 10 electrode arrays. Geophysical Prospecting 52, 379–398. Daniels, J.J., 1978. Interpretation of buried electrode resistivity data using the layered earth model. Geophysics 43, 988–1001. Dey, A., Morrison, H.F., 1979. Resistivity modelling for arbitrarily shaped twodimensional structures. Geophys. Prospecting 27 (1), 106–136. Dutta, D.J., 1994. TRANS4: a FORTRAN program for computing apparent resistivity departure curves for an infinitely thick bed with transitional invaded zone in borehole geophysics. Computers & Geosciences 20 (3), 293–311. deLugao, P.P., Wannamaker, P., 1996. Calculating the two-dimensional magnetotelluric Jacobian in finite elements using reciprocity. Geophysical Journal International 127, 806–810. deLugão, P.B., Portniaguine, O., Zhdanov, M.S., 1997. Fast and stable two-dimensional inversion of magnetotellurics data. J. Geomag. Geoelectr. 49, 1469–1497. Haber, E., Oldenburg, D., 1998. A GCV based method for non-linear ill-posed problems. Computational Geosciences 4, 41–63. Kwon, Y.W., Bang, H., 2000. The Finite Element Method using MATLAB, 2nd Edition. CRC Press. LaBrecque, D.J., Oldenborger, G.A., Sharpe, R., Knoll, M.D., 2006. Axi-symetric inversion of electrical resistivity tomography data to monitor the movement of fluids injected in wells. Symposium on the Application of geophysics to Environmental and engineering problems, p. 1505-1503. Liu, Z., Lin, H., 2002. Joint inversion of induction/lateral/normal logs, case studies at Shenli field site, China. Journal of Petroleum Science and Engineering 34, 55–64. Mackie, R.L., Madden, T.R., 1993. Three dimensional magnetotelluric inversion using conjugate gradients. Geophysical Journal International 115, 215–229. Mufti, I.R., 1976. Finite-difference resistivity modeling for arbitrarily shaped twodimensional structures. Geophysics 41, 62–78. Pardo, D., Paszynski, M., Torres-verdin, C., Demkowicz, L., 2007. Simulation of 3D DC borehole resistivity measurements with a goal-oriented hp finite-element method, Part I: laterolog and LWD. Journals of the Serbian Society for Computational Mechanics 1, 62–73. Pidlisecky, A., Haber, E., Knight, R., 2007. RESINVM3D: a 3D resistivity inversion package. Geophysics 72, H1–H10. Portniaguine, O., Zhdanov, M.S., 1999. Focusing geophysical inversion images. Geophysics 64, 874–887. Rodi, W.L., 1976. A technique for improving the accuracy of finite element solutions for magnetotelluric data. Geophysical Journal of the Royal Astronomical. Society 44, 483–506. Rodi, W.L., Mackie, R.L., 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics 66, 174–187. Roy, K.K., Dutta, D.J., 1997. Relative performance of focussed and unfocussed devices in a two-dimensional borehole environment. Journal of Applied Geophysics 36, 181–194. Tikhonov, A.N., Arsenin, V.Y., 1977. Solutions of Ill-Posed Problems. V.H. Winston and Sons, Washington, D.C. Towle, G.H., Whitman, W.W., Kim, J.H., 1988. Electric log modeling with a finite difference method. The Log Analyst 184–195 (May–June). Yang, F.W., Ward, S.H., 1984. Inversion of borehole normal resistivity logs. Geophysics 49, 1541–1548. Zhang, G.J., Shen, L.C., 1984. Response of a normal resistivity tool in a borehole crossing a bed boundary. Geophysics 49 (2), 142–149. Zhang, J., Mackie, R.L., Madden, T.R., 1995. 3-D resistivity forward modeling and inversion using conjugate gradients. Geophysics 60, 1313–1325. Zhdanov, M.S., 2002. Geophysical Inverse Theory and Regularization Problems. Elsevier Science Pub\.Co., Inc. Zhdanov, M.S., Tolstaya, E., 2004. Minimum support nonlinear parameterization in the solution of 3D magnetotelluric inverse problem. Inverse Problems 20, 937–952.