US20020040274A1 - Method for 2D inversion of dual laterolog measurements - Google Patents

Method for 2D inversion of dual laterolog measurements Download PDF

Info

Publication number
US20020040274A1
US20020040274A1 US09/961,590 US96159001A US2002040274A1 US 20020040274 A1 US20020040274 A1 US 20020040274A1 US 96159001 A US96159001 A US 96159001A US 2002040274 A1 US2002040274 A1 US 2002040274A1
Authority
US
United States
Prior art keywords
earth model
tool
inversion
tool response
right arrow
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
US09/961,590
Other versions
US6430509B1 (en
Inventor
Hezhu Yin
Hanming Wang
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
ExxonMobil Upstream Research Co
Original Assignee
ExxonMobil Upstream Research Co
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by ExxonMobil Upstream Research Co filed Critical ExxonMobil Upstream Research Co
Priority to US09/961,590 priority Critical patent/US6430509B1/en
Assigned to EXXONMOBIL UPSTREAM RESEARCH COMPANY reassignment EXXONMOBIL UPSTREAM RESEARCH COMPANY ASSIGNMENT OF ASSIGNORS INTEREST (SEE DOCUMENT FOR DETAILS). Assignors: WANG, HANMING, YIN, HEZHU
Publication of US20020040274A1 publication Critical patent/US20020040274A1/en
Application granted granted Critical
Publication of US6430509B1 publication Critical patent/US6430509B1/en
Anticipated expiration legal-status Critical
Expired - Lifetime legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction

Definitions

  • This invention relates generally to the field of geophysical prospecting. More particularly, the invention relates to an automated two-dimensional inversion method for true resistivity of reservoir formations from a dual laterolog (DLL) tool's borehole measurements.
  • DLL dual laterolog
  • chartbooks provided by logging service companies.
  • chartbooks only contain a limited number of charts with strict assumptions (e.g., borehole diameter, mud resistivity, and Rt/Rxo ratio) that a chart can seldom match real world examples. Therefore, chartbook corrections may only serve to make a qualitative estimation.
  • the nonlinear resistivity tool response due to borehole diameter, mud resistivity, invasion, and bed thickness or shoulder bed effects all together) can not be corrected from the chartbooks' corrections without assumptions of linear superposition.
  • such an inversion process attempts to fit the computed tool response under a set of earth parameters (e.g., bed thickness, Rt, Rm, Rxo, borehole diameter and invasion depth) to an actual field resistivity log, or a set of actual field logs.
  • the parameters in the earth model can be refined by solving least-squares problems through the iterative process to minimize the sum of the squares of the errors between the computed tool response and the measured field log. The iteration may continue until the fit between the computed and field logs reaches predetermined criteria.
  • the present invention is a method for 2D inversion of true resistivity from dual laterolog tool measurements, using a pre-calculated look-up table.
  • a 2D tool response is calculated in each interval using the earth model. Matching is checked between the calculated 2D tool response and the tool measurements. The following steps are iterated until the match is satisfactory.
  • a 1D radial tool response is derived in each interval using the pre-calculated look-up table. Shoulder bed effects are approximated in each interval by subtracting the 1D radial tool response from the 2D tool response.
  • a non-linear least square optimization is applied at boundaries of the intervals and local maximum and minimum values in the intervals to update the earth model.
  • FIG. 1 is a flowchart illustrating an embodiment of the method of the present invention for 2D inversion of true formation resistivity from dual laterolog measurements;
  • FIGS. 2 a and 2 b are diagrams schematically illustrating the step in the present invention for the shoulder bed effect approximation to reduce matrix size and computation time;
  • FIG. 3 is a diagram schematically illustrating the 1D radial formation and borehole model used in a test of the present invention
  • FIG. 4 is a diagram schematically illustrating the 3-layer formation model used in a test of the present invention.
  • FIG. 5 a is a diagram schematically illustrating the 7-layer formation model used in a test of the present invention.
  • FIGS. 5 b and 5 c are graphs comparing the formation resistivity and invasion depth from the test of FIG. 5 a as measured and as obtained by the method of the present invention.
  • the present invention is a computationally efficient method for 2D inversion of true formation resistivity from the dual laterolog logging tool measurements.
  • the 2D dual laterolog tool response is decomposed into 1D radial and 2D responses in intervals.
  • the difference between the 1D radial response and the initial 2D response is attributed to the shoulder bed effect.
  • the dimensions of the Jacobian matrix in solving the least-square problem are significantly reduced, thereby drastically reducing the computational effort and time required to perform the inversion.
  • the results from theoretical benchmark tests and field log inversions have proven the method of the present invention to be a computationally efficient 2D inversion technique with sufficient accuracy for practical use.
  • the dual laterolog tool injects low-frequency alternating current into the formation. This process is described in Doll, H, 1951, “The laterolog: A new resistivity logging method with electrodes using an automatic focusing system”, AIME Trans. Petroleum, Vol. 192, PP. 305-315 and Suau, J., Grimadli, P., Poupon, A., and Souhaite, P., 1972, “The dual laterolog: Rxo tool,” SPE 4018.
  • the physical earth model can be simplified as a two-dimensional mathematical model where the model parameters vary symmetrically about borehole axis. The parameters are a function of radial direction and borehole vertical direction.
  • a strict 2D inversion method can address the coupled effects of shoulder beds (vertical), and borehole and invasion (radial) on the nonlinear tool response, but the computation is too intensive for a practical formation evaluation purpose, such as inversion.
  • One of the objectives of the 2D inversion method of the present invention is to de-couple the vertical and radial responses to reduce computation during preliminary iterations.
  • the nonlinear coupled-effects are calculated in strict 2D only at the very beginning of the iteration, and near the end of the inversion.
  • FIG. 1 is a flowchart illustrating an embodiment of the method of the present invention for 2D inversion of true formation resistivity from dual laterolog measurements.
  • an initial earth model is derived from available logging data, including the dual laterolog measurements. If there are N horizontal layers in the earth formation model, a vector ⁇ right arrow over (X) ⁇ j ⁇ db j , di j , Rm j , Rxo j , Rt j . . . ⁇ 1 ⁇ j ⁇ N represents the parameters of the j th layer of the earth model.
  • db j , di j , Rm j , Rxo j , Rt j are the borehole diameter, depth of invasion, mud resistivity, invasion zone resistivity, and uninvaded formation true resistivity of the j th layer of the earth model, respectively.
  • the initial earth model is preferably pre-processed. Pre-inversion processing of the parameter data develops an initial earth model that honors geological and geophysical constraints on defined bed boundaries, invasion depth, resistivity of invaded zones, and initial resistivity of uninvaded zone for each formation. Constraints from other logging data rather than the dual laterolog log itself are not only useful to reduce matrix size and computation time, but also to derive a more robust inversion result and to reduce opportunities for non-uniqueness.
  • the present invention preferably includes, but is not restricted to, the following preprocessing steps to set constraints.
  • the bed boundaries and number of intervals for the inversion process are derived with higher resolution logs, borehole images, or core descriptions or photographs.
  • invasion zone resistivity, Rxo is derived by averaging the log measurements from spherical focused laterologs or micro-spherical focused laterologs over each derived bed by simple arithmetic mean or weighted mean.
  • shaliness of each interval, Vsh is defined by using a normalization of gamma ray log measurements,
  • V sh ( GR ⁇ GR min )/( GR max ⁇ GR min ),
  • step 102 the earth model from step 101 is subdivided into intervals.
  • the intervals will be used below to make the inversion process more computationally efficient.
  • the number of sub-divided intervals is the same as the number of layers in the earth model. Then, as will be discussed below, the method of the present invention will be most efficient, since the problem of calculating a tool response in every interval becomes a 1D radial problem instead of a 2D problem.
  • step 103 a 2D dual laterolog tool response is calculated using the earth model derived in step 101 .
  • the calculation is done for each of the intervals selected in step 102 .
  • This calculation is discussed below in the mathematical treatment concerning Eq. (1).
  • a finite element method forward modeling code is preferably used for the 2D dual laterolog tool response calculation.
  • the calculation assumes a complex 2D borehole environment and is based on rigorous solution of electrostatic potential equations, efficient discretization, and automatic mesh generation for up to 500 layers.
  • the present invention preferably uses fast front-solver software for reducing the size of matrix operations.
  • step 104 a comparison is made to determine whether the calculated 2D tool response from step 103 and the original field resistivity measurements of the dual laterolog tool substantially match. If the determination is that the match is satisfactory, then the process continues to step 108 below for post-processing and output of the earth model. Otherwise, if the determination is that the match is not satisfactory, then the process continues to step 105 through 107 to update the earth model before checking the match again.
  • This matching process is described in Druskin, V., and Knizhnerman, L., 1987, “About one iterative algorithm for solving two-dimensional inverse problem of logging by lateral sounding” Geology and Geophysics, Vol., 9 , P111-116.
  • a 1D radial tool response is derived from a pre-calculated look-up table. The derivation is done for each of the intervals from step 102 . A pseudo-inversion of the dual laterolog tool response may be thus quickly approximated. Due to the complexity of the electrode configurations and focusing conditions for the deep measurements from a dual laterolog tool (commercial tools are currently available with 11 or 13 electrodes), a numerical solution of the 1D radial tool response for the dual laterolog tool may not be readily achieved.
  • a look-up table is pre-calculated once from the code to be used during the preliminary iterations to quickly approximate the dual laterolog response in such a 1D radial case.
  • the variables of the look-up table include the diameter of the borehole (db), the ratio of mud resistivity and invaded zone resistivity (Rm/Rxo), the ratio of formation resistivity and invaded zone resistivity (Rt/Rxo) and depth of invasion (di).
  • db the diameter of the borehole
  • Rm/Rxo the ratio of mud resistivity and invaded zone resistivity
  • Rt/Rxo the ratio of formation resistivity and invaded zone resistivity
  • di depth of invasion
  • db is varied from 6′′ to 32′′ with 2′′ increments
  • Rm/Rxo takes values of 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, and 10.
  • Rt/Rxo takes values of 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 2000, 5000, and 10000.
  • the total number of the 1D radial responses in a look-up table is 383,240 (13 ⁇ 11 ⁇ 67 ⁇ 40) for each type of dual laterolog tool, i.e., DLS-B, DLS-D, or DLS-E tools. Numerical interpolation is performed when the parameters fall between any two numbers in the look-up table.
  • step 106 shoulder bed effects are approximated by subtracting the 1D radial tool response interpolated from the look-up table in step 105 from the 2D tool response calculated in step 103 . The approximation is done for each of the intervals from step 102 . This shoulder bed approximation is discussed below in the mathematical treatment concerning Eq. (6).
  • step 107 the earth model is updated by solving the normal equation linking the measurements from the dual laterolog tool, the 2D dual laterolog tool response calculated in step 103 , and the shoulder bed effects calculated in step 106 .
  • This normal equation is discussed below in the mathematical treatment concerning Eq. (10).
  • the calculation is done for each of the intervals from step 102 . Preferably, this calculation is done by non-linear least square optimization. This calculation is schematically illustrated in FIGS. 2 a and 2 b.
  • FIG. 2 a shows a borehole 201 , the shoulder beds 202 above and below the interval, the invaded portion 203 of the interval, and the uninvaded portion 204 of the interval.
  • FIG. 2 a shows a borehole 201 , the shoulder beds 202 above and below the interval, the invaded portion 203 of the interval, and the uninvaded portion 204 of the interval.
  • 2 b shows the shoulder bed effect 210 and the combined borehole effect and invasion effect 211 relative to the 1D radial tool response 209 and the 2D dual laterolog tool response 208 .
  • this calculation is only done at the defined bed-boundaries 205 , local maximum 206 and minimum points 207 in the intervals, instead of calculating at every sampling point in the dual laterolog measurements 208 . This procedure decreases the size of the computations considerably.
  • step 103 After the earth model has been updated by the procedure in steps ( 105 ) through ( 107 ), the process then returns to step 103 to recalculate the dual laterolog 2D tool response with the new earth model parameters.
  • step 108 the results of the dual laterolog tool inversion calculated in step 103 are post-processed.
  • Post-inversion processing is preferred to avoid non-physical and non-logical inversion results which may be yielded from the poor sensitivity of the dual laterolog tool to changes of the invasion depth, di, and invasion zone resistivity, Rxo, in some cases.
  • the present invention preferably includes, but is not restricted to, the following post-processing steps. First, if
  • step 109 the parameters comprising the earth model are output.
  • the N horizontal layers of the earth model are represented by the vector ⁇ right arrow over (X) ⁇ j ⁇ db j , di j , Rm j , Rxo j , Rt j . . . ⁇ 1 ⁇ j ⁇ N , where db j , di j , Rm j , Rxo j , Rt j are the borehole diameter, depth of invasion, mud resistivity, invasion zone resistivity, and uninvaded formation true resistivity of the j th layer of the earth model, respectively.
  • the well logging measurement may be expressed as the product of the mapping operator G and earth model vector ⁇ right arrow over (X) ⁇ to form a tool response ⁇ right arrow over (R) ⁇ in vector space:
  • the operator G may be a nonlinear as function of ⁇ right arrow over (X) ⁇ , and EQ. (1) can be solved numerically for ⁇ right arrow over (R) ⁇ using finite element method forward modeling code.
  • ⁇ 2 ⁇ right arrow over (R) ⁇ right arrow over (D) ⁇ ⁇ (2)
  • ⁇ right arrow over ( ⁇ X) ⁇ is the increment step in modifying ⁇ right arrow over (X) ⁇ in the iteration process.
  • H is the Jacobian matrix and an element of H is a set of partial derivatives:
  • ⁇ ⁇ right arrow over (R) ⁇ / ⁇ right arrow over (X) ⁇ j ⁇ right arrow over (R) ⁇ / ⁇ db j + ⁇ right arrow over (R) ⁇ / ⁇ di j + ⁇ right arrow over (R) ⁇ / ⁇ Rm j + ⁇ right arrow over (R) ⁇ / ⁇ Rxo j + ⁇ right arrow over (R) ⁇ / ⁇ Rt j + . . . (4)
  • ⁇ right arrow over (R) ⁇ i ′ represents the predicted logging tool responses in the ⁇ right arrow over (X) ⁇ i ′ interval when ⁇ right arrow over (X) ⁇ i ′ exists independently.
  • the difference between ⁇ right arrow over (R) ⁇ i and ⁇ right arrow over (R) ⁇ i ′ is defined as the shoulder bed effect ⁇ right arrow over (S) ⁇ i ′ over the ⁇ right arrow over (X) ⁇ 1 ′ , ⁇ right arrow over (X) ⁇ 2 ′ , . . . , ⁇ right arrow over (X) ⁇ x ⁇ 1 ′ , ⁇ right arrow over (X) ⁇ i ′ , ⁇ right arrow over (X) ⁇ i+1 ′ , . . . , ⁇ right arrow over (X) ⁇ m ′ division.
  • x i ′2
  • H i ′ (or H) is nearly singular. This is because the tool response corresponding to the changes of the formation resistivity or invasion radius has poor sensitivity.
  • a damping parameter ⁇ can be introduced. This is described in Levenberg, K., 1944, “A method for the solution of certain nonlinear problems in least square, Quarterly of Applied Mathematics, 2, p164-168, and in Marquardt, D., W., 1963, “An algorithm for least square estimate of non-linear parameters, Journal of the Society of Industrial and Applied Mathematics, 11, p 431-441. Then equations (3) and (8), respectively, can be modified as:
  • equation (9) or equation (10) reverts to equation (3) or (8), respectively, and can be solved by a Newton-Gauss elimination method.
  • ⁇ is large, the matrix (H T H+ ⁇ I) ⁇ right arrow over ( ⁇ X) ⁇ (or (H ′T H ′ + ⁇ I) ⁇ right arrow over ( ⁇ X) ⁇ ′ ) is nearly diagonal and can be solved by the gradient method.
  • the singular value decomposition (SVD) method is chosen to solve the above normal equations for the computational benefits, as described in Lines, L. R., and Treitel, S., 1984, “A review of least-square inversion and its application to geophysical problems”, Geophysical Prospecting, Vol., 32, p.
  • the initial value of ⁇ is set to be 0.001.
  • the cost function is computed and compared with the one in previous step. If it becomes smaller than the one in the previous step, ⁇ is decreased by a factor of 10. Otherwise, ⁇ is increased by a factor of 10.
  • ⁇ 2 merit function is smaller than a pre-set value, a new set of parameters ⁇ right arrow over (X) ⁇ i ′ will be derived.
  • the above described method will be most efficient because the dimension of H i ′ reaches its minimum.
  • the tool response solution becomes a 1D radial solution in every sub-interval, i.e., the parameters ⁇ right arrow over (X) ⁇ i ′ are 1D in the radial direction.
  • the Jacobian matrix H i ′ described by equation (10) in the present invention may have four elements in maximum at a given depth and bed:
  • H i ′ ( ⁇ ( LLD )/ ⁇ Rt , ⁇ ( LLS )/ ⁇ Rt , ⁇ ( LLD )/ ⁇ di , ⁇ ( LLS )/ ⁇ di ). (11)
  • H i ′ ( ⁇ ( LLD )/ ⁇ Rt , ⁇ ( LLS )/ ⁇ Rt ). (12)
  • Benchmark model inversion is useful to evaluate how accurate the inversion method of the present invention can be. Using benchmark models, the formation parameters (such as Rt) are given in the model and hence, are known and the inversion results can be compared with the “true” value.
  • FIG. 3 shows a borehole 301 with diameter db 302 , surrounded by an invaded zone 303 with a diameter of invasion di 304 , sitting in an uninvaded bed 305 .
  • the LLD and LLS “log data” are computed under various Rt models using the 2D finite element method forward modeling code and used as “measured data” in the inversion. Under given Rxo, Rm, and db conditions, Rt and di are inverted from LLD and LLS laterolog responses.
  • the look-up table and interpolation technique is used to approximate the forward modeling process.
  • Table 1 The results show excellent agreement between the given Rt and di in the model and the inverted Rt and di from LLD and LLS responses with the look-up table.
  • the relative error is less than ⁇ 0.1% for Rt and ⁇ 2% for di, which is much more accurate than the laterolog tool measurement accuracy ( ⁇ 5% in resistivity) defined in the dual laterolog tool specification.
  • These 1D radial cases demonstrate that the look-up table and interpolation method is an excellent approximation of dual laterolog responses in a thick formation. Inversion with the look-up table and interpolation is thousands of times faster than a strict 2D tool response forward modeling-based inversion.
  • FIG. 4 shows a borehole 401 with diameter db 402 , surrounded by an invaded formation 403 with diameter of invasion di 404 and height h 406 , inside an uninvaded portion 405 and sandwiched by two infinitely thick shoulder beds 407 .
  • Rxo 5 ohm-m
  • Rs 1 ohm-m
  • Rm 0.1 ohm-m
  • db 8 inch
  • center bed resistivity Rt 20 ohm-m
  • the center bed thickness is varied from one to three meters.
  • a full 2D inversion, described by equation (9), is performed to compare with the present invention, described by equation (10). The results are shown in Table 2.
  • the resulting Rt and di from the two inversion methods are fairly similar, but the computation time required is quite different.
  • the full 2D inversion takes more than 5 minutes on the SUN SPARC-20 workstation, while the 2D inversion of the present invention needs less than 30 seconds for the same cases.
  • the present invention is about ten times faster than the full 2D inversion method.
  • the present invention can be hundreds of times faster than the full 2D inversion when inversion is performed over hundreds of beds.
  • FIG. 5 a shows a 7-layer formation benchmark with the 3 resistive and invaded beds is tested with the 2D inversion code of the present invention.
  • FIG. 5 a shows a borehole 501 with diameter db 502 , surrounded by 3 invaded formations 503 of diameter of invasion di 504 and height h 506 , inside an uninvaded portion 505 , 2 uninvaded formations 507 of height h 508 and sandwiched by two infinitely thick shoulder beds 507 .
  • 5 b and 5 c show that the inverted formation resistivity Rt 510 and invasion depth di 512 , respectively, agree well with the original formation resistivity Rt 509 and invasion depth di 511 set in the model.
  • the iteration is stopped when the relative difference between the computed LLD and LLS responses in the inversion process and the field dual laterolog log is less than 5% in a given bed, which is dual laterolog tool measurement accuracy. It can be seen that the relative difference between the inverted formation resistivity Rt 510 and the original formation resistivity Rt 509 set in the model is also less than 5% in this example.
  • the above benchmark inversion examples demonstrate that the 2D dual laterolog inversion method of the present invention is a fast and solid inversion technique with sufficient accuracy.
  • the advantages of the present invention include, but are not restricted to, computationally efficient dual laterolog resistivity inversion, more accurate hydrocarbon pore volume estimation analysis, better quantitative formation evaluation, more accurate petrophysical input (such as hydrocarbon pore volume) into reserve assessments, and improved identification of potentially by-passed hydrocarbon accumulations.

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Measurement Of Resistance Or Impedance (AREA)
  • Investigation Of Foundation Soil And Reinforcement Of Foundation Soil By Compacting Or Drainage (AREA)

Abstract

A 2D inversion of true formation resistivity from dual laterolog tool measurements is accomplished using a pre-calculated look-up table. An initial earth model is derived and divided into intervals. A 2D tool response is calculated in each interval using the earth model. Matching is checked between the calculated 2D tool response and the tool measurements. The following steps are iterated until the match is satisfactory. A 1D radial tool response is derived in each interval using the pre-calculated look-up table. Shoulder bed effects are approximated in each interval by subtracting the 1D radial tool response from the 2D tool response. A non-linear least square optimization is applied at boundaries of the intervals and local maximum and minimum values in the intervals to update the earth model.

Description

  • This application claims the benefit of U.S. Provisional Application No. 60/237,749 filed Oct. 3, 2000.[0001]
    • FIELD OF THE INVENTION
  • This invention relates generally to the field of geophysical prospecting. More particularly, the invention relates to an automated two-dimensional inversion method for true resistivity of reservoir formations from a dual laterolog (DLL) tool's borehole measurements. [0002]
    • BACKGROUND OF THE INVENTION
  • The accuracy of hydrocarbon saturation calculated from well logs, an essential component of petrophysical resource assessment, is fundamentally dependent upon accurate determination of the true formation resistivity (Rt). Apparent formation resistivity (Ra), as measured by a logging tool, however, is not equal to true formation resistivity (Rt) in most logging environments because of the limitations of tool physics and non-ideal borehole conditions. It is known in the art that deep-reading resistivity tools cannot resolve formations less than a few feet thick, and cannot make accurate true resistivity (Rt) measurements when the borehole diameter is variable (rugosity), and when the borehole fluid with a different resistivity than formation fluids has seeped into the formation (invasion), thereby altering the resistivity of the invaded zone (Rxo). [0003]
  • The traditional method of correcting these environmental effects on resistivity logs has been to use chartbooks provided by logging service companies. However, chartbooks only contain a limited number of charts with strict assumptions (e.g., borehole diameter, mud resistivity, and Rt/Rxo ratio) that a chart can seldom match real world examples. Therefore, chartbook corrections may only serve to make a qualitative estimation. Furthermore, the nonlinear resistivity tool response (due to borehole diameter, mud resistivity, invasion, and bed thickness or shoulder bed effects all together) can not be corrected from the chartbooks' corrections without assumptions of linear superposition. [0004]
  • Computer inverse modeling of resistivity tool response can be conducted to convert apparent resistivity from logs into a response profile that may closely approximate reality. In fact, modern environmental correction charts provided by service companies are the result of computer forward modeling. In general, the inverse modeling involves replicating the observed field log by numerically solving the mathematical boundary value problems of the electrical or electromagnetic fields generated by a specific resistivity tool under a predefined layered-earth model. To the degree that the field log and the computed tool response are in acceptable agreement through iterative forward modeling, the underlying earth model may be considered as one possible representation of the formation's true resistivity profile. Mathematically, such an inversion process attempts to fit the computed tool response under a set of earth parameters (e.g., bed thickness, Rt, Rm, Rxo, borehole diameter and invasion depth) to an actual field resistivity log, or a set of actual field logs. The parameters in the earth model can be refined by solving least-squares problems through the iterative process to minimize the sum of the squares of the errors between the computed tool response and the measured field log. The iteration may continue until the fit between the computed and field logs reaches predetermined criteria. [0005]
  • With the advent of modeling codes and the significant increase in computing power, resistivity tool response modeling has become a feasible option for formation evaluation. Strict 2D inversion of resistivity logging tool measurements based on iterative 2D forward modeling with finite element or hybrid methods are described in Gianzero, S., Lin, Y., and Su, S., 1985, “A new high speed hybrid technique for simulation and inversion of resistivity logs”, SPE 14189; Liu, Q., H., 1994, “Nonlinear inversion of electrode-type resistivity measurements”, IEEE on Geoscience and Remote Sensing, Vol. 32, No. 3, pp 499-507; and Mezzatesta, A. G., Eckard, M. H., Strack, K. M., 1995, “Integrated 2D Interpretation of Resistivity Logging Measurements by Inversion Methods”, Paper-E, SPWLA 36th Annual Logging Symposium Transactions. However, these methods are very computer time consuming. The massive computation required by some of the above methods may only be performed on super computers, making its application to well logging interpretation impractical. [0006]
  • Thus, the need exists for a computationally efficient method for 2D inversion of resistivity from dual laterolog measurements. [0007]
    • SUMMARY OF THE INVENTION
  • The present invention is a method for 2D inversion of true resistivity from dual laterolog tool measurements, using a pre-calculated look-up table. First, an initial earth model is derived and divided into intervals. A 2D tool response is calculated in each interval using the earth model. Matching is checked between the calculated 2D tool response and the tool measurements. The following steps are iterated until the match is satisfactory. A 1D radial tool response is derived in each interval using the pre-calculated look-up table. Shoulder bed effects are approximated in each interval by subtracting the 1D radial tool response from the 2D tool response. A non-linear least square optimization is applied at boundaries of the intervals and local maximum and minimum values in the intervals to update the earth model. [0008]
    • BRIEF DESCRIPTION OF THE DRAWINGS
  • The present invention and its advantages may be more easily understood by reference to the following detailed description and the attached drawings in which: [0009]
  • FIG. 1 is a flowchart illustrating an embodiment of the method of the present invention for 2D inversion of true formation resistivity from dual laterolog measurements; [0010]
  • FIGS. 2[0011] a and 2 b are diagrams schematically illustrating the step in the present invention for the shoulder bed effect approximation to reduce matrix size and computation time;
  • FIG. 3 is a diagram schematically illustrating the 1D radial formation and borehole model used in a test of the present invention; [0012]
  • FIG. 4 is a diagram schematically illustrating the 3-layer formation model used in a test of the present invention; [0013]
  • FIG. 5[0014] a is a diagram schematically illustrating the 7-layer formation model used in a test of the present invention; and
  • FIGS. 5[0015] b and 5 c are graphs comparing the formation resistivity and invasion depth from the test of FIG. 5a as measured and as obtained by the method of the present invention.
  • While the invention will be described in connection with its preferred embodiments, it will be understood that the invention is not limited thereto. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the spirit and scope of the invention, as defined by the appended claims. [0016]
    • DETAILED DESCRIPTION OF THE INVENTION
  • The present invention is a computationally efficient method for 2D inversion of true formation resistivity from the dual laterolog logging tool measurements. In the present invention, the 2D dual laterolog tool response is decomposed into 1D radial and 2D responses in intervals. The difference between the 1D radial response and the initial 2D response is attributed to the shoulder bed effect. Under this inversion scheme, the dimensions of the Jacobian matrix in solving the least-square problem are significantly reduced, thereby drastically reducing the computational effort and time required to perform the inversion. The results from theoretical benchmark tests and field log inversions have proven the method of the present invention to be a computationally efficient 2D inversion technique with sufficient accuracy for practical use. [0017]
  • The dual laterolog tool injects low-frequency alternating current into the formation. This process is described in Doll, H, 1951, “The laterolog: A new resistivity logging method with electrodes using an automatic focusing system”, AIME Trans. Petroleum, Vol. 192, PP. 305-315 and Suau, J., Grimadli, P., Poupon, A., and Souhaite, P., 1972, “The dual laterolog: Rxo tool,” SPE 4018. In a vertical wellbore with horizontal bed boundaries, the physical earth model can be simplified as a two-dimensional mathematical model where the model parameters vary symmetrically about borehole axis. The parameters are a function of radial direction and borehole vertical direction. [0018]
  • Fast and accurate forward modeling computer software plays an important role in the forward modeling-based iterative inversion process. The efficiency of an inversion algorithm is partially dependent on how much faster the forward modeling part is. In the present work, a robust finite element method for electrical logging simulation is adapted to produce a strict 2D dual laterolog tool response. In the finite element code, an automatic mesh generation method is introduced and an efficient linear equation solver from Irons, B., M., 1970, “A frontal solution program for finite element analysis”, Int. J. for Num. Methods in Eng., 2, 5-32, is used. The code can handle up to 500 layers in the vertical direction and 10 invasion zones in the radial direction. The computation time is almost independent of the total number of layers. Calculations of both the deep (LLD) and shallow (LLS) dual log responses take 3 seconds per logging station on a SUN SPARC 20 platform. Details of the finite element method for dual laterolog forward modeling are discussed by Zhang, G., J., 1986, “Electrical logging (II),” Publication House of petroleum Industry, China. [0019]
  • A strict 2D inversion method can address the coupled effects of shoulder beds (vertical), and borehole and invasion (radial) on the nonlinear tool response, but the computation is too intensive for a practical formation evaluation purpose, such as inversion. One of the objectives of the 2D inversion method of the present invention is to de-couple the vertical and radial responses to reduce computation during preliminary iterations. The nonlinear coupled-effects are calculated in strict 2D only at the very beginning of the iteration, and near the end of the inversion. [0020]
  • FIG. 1 is a flowchart illustrating an embodiment of the method of the present invention for 2D inversion of true formation resistivity from dual laterolog measurements. First, in [0021] step 101, an initial earth model is derived from available logging data, including the dual laterolog measurements. If there are N horizontal layers in the earth formation model, a vector {right arrow over (X)}j∈{dbj, dij, Rmj, Rxoj, Rtj . . . }1≦j≦N represents the parameters of the jth layer of the earth model. Here, dbj, dij, Rmj, Rxoj, Rtj are the borehole diameter, depth of invasion, mud resistivity, invasion zone resistivity, and uninvaded formation true resistivity of the jth layer of the earth model, respectively.
  • The initial earth model is preferably pre-processed. Pre-inversion processing of the parameter data develops an initial earth model that honors geological and geophysical constraints on defined bed boundaries, invasion depth, resistivity of invaded zones, and initial resistivity of uninvaded zone for each formation. Constraints from other logging data rather than the dual laterolog log itself are not only useful to reduce matrix size and computation time, but also to derive a more robust inversion result and to reduce opportunities for non-uniqueness. [0022]
  • Based on available geological and geophysical knowledge, the present invention preferably includes, but is not restricted to, the following preprocessing steps to set constraints. First, the bed boundaries and number of intervals for the inversion process are derived with higher resolution logs, borehole images, or core descriptions or photographs. Second, invasion zone resistivity, Rxo, is derived by averaging the log measurements from spherical focused laterologs or micro-spherical focused laterologs over each derived bed by simple arithmetic mean or weighted mean. Third, shaliness of each interval, Vsh, is defined by using a normalization of gamma ray log measurements, [0023]
  • V sh=(GR−GR min)/(GR max −GR min),
  • and then setting an initial invasion depth, di, with a user defined cut-off value of the shaliness, V[0024] cutoff. Thus, if Vsh(j)≧Vcutoff, then there is no invasion for shale formation. If Vsh(J)≦Vcutoff, then the initial invasion depth, di(j) for the jth interval is set to be 0.3 m, (about 11.8″). This third step is very effective in reducing the computation time because there will be no derivatives for the depth of invasion (∂{right arrow over (R)}/∂dij) to be calculated for the Jacobian matrix when there is no invasion. Inversion becomes a 1D issue with borehole effect considered when there is no invasion. Fourth, initial true resistivity Rt is derived for each interval with the deep laterolog from the center-point of each bed.
  • In [0025] step 102, the earth model from step 101 is subdivided into intervals. The intervals will be used below to make the inversion process more computationally efficient. Preferably, the number of sub-divided intervals is the same as the number of layers in the earth model. Then, as will be discussed below, the method of the present invention will be most efficient, since the problem of calculating a tool response in every interval becomes a 1D radial problem instead of a 2D problem.
  • In [0026] step 103, a 2D dual laterolog tool response is calculated using the earth model derived in step 101. The calculation is done for each of the intervals selected in step 102. This calculation is discussed below in the mathematical treatment concerning Eq. (1). A finite element method forward modeling code is preferably used for the 2D dual laterolog tool response calculation. The calculation assumes a complex 2D borehole environment and is based on rigorous solution of electrostatic potential equations, efficient discretization, and automatic mesh generation for up to 500 layers. In particular, the present invention preferably uses fast front-solver software for reducing the size of matrix operations.
  • In [0027] step 104, a comparison is made to determine whether the calculated 2D tool response from step 103 and the original field resistivity measurements of the dual laterolog tool substantially match. If the determination is that the match is satisfactory, then the process continues to step 108 below for post-processing and output of the earth model. Otherwise, if the determination is that the match is not satisfactory, then the process continues to step 105 through 107 to update the earth model before checking the match again. This matching process is described in Druskin, V., and Knizhnerman, L., 1987, “About one iterative algorithm for solving two-dimensional inverse problem of logging by lateral sounding” Geology and Geophysics, Vol., 9, P111-116.
  • In [0028] step 105, a 1D radial tool response is derived from a pre-calculated look-up table. The derivation is done for each of the intervals from step 102. A pseudo-inversion of the dual laterolog tool response may be thus quickly approximated. Due to the complexity of the electrode configurations and focusing conditions for the deep measurements from a dual laterolog tool (commercial tools are currently available with 11 or 13 electrodes), a numerical solution of the 1D radial tool response for the dual laterolog tool may not be readily achieved. Instead of running 2D finite element method forward modeling code in each iteration for the 1D radial response, a look-up table is pre-calculated once from the code to be used during the preliminary iterations to quickly approximate the dual laterolog response in such a 1D radial case.
  • The variables of the look-up table include the diameter of the borehole (db), the ratio of mud resistivity and invaded zone resistivity (Rm/Rxo), the ratio of formation resistivity and invaded zone resistivity (Rt/Rxo) and depth of invasion (di). In generation of the table, the following ranges of variable have been used in the look-up table generation: [0029]
  • (1) db is varied from 6″ to 32″ with 2″ increments; [0030]
  • (2) Rm/Rxo takes values of 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, and 10. [0031]
  • (3) Rt/Rxo takes values of 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 2000, 5000, and 10000. [0032]
  • (4) di is varied from borehole radius to 78.74 inches (2 m) with 1.97 inch (5 cm) increments. [0033]
  • The total number of the 1D radial responses in a look-up table is 383,240 (13×11×67×40) for each type of dual laterolog tool, i.e., DLS-B, DLS-D, or DLS-E tools. Numerical interpolation is performed when the parameters fall between any two numbers in the look-up table. [0034]
  • In [0035] step 106, shoulder bed effects are approximated by subtracting the 1D radial tool response interpolated from the look-up table in step 105 from the 2D tool response calculated in step 103. The approximation is done for each of the intervals from step 102. This shoulder bed approximation is discussed below in the mathematical treatment concerning Eq. (6).
  • In [0036] step 107, the earth model is updated by solving the normal equation linking the measurements from the dual laterolog tool, the 2D dual laterolog tool response calculated in step 103, and the shoulder bed effects calculated in step 106. This normal equation is discussed below in the mathematical treatment concerning Eq. (10). The calculation is done for each of the intervals from step 102. Preferably, this calculation is done by non-linear least square optimization. This calculation is schematically illustrated in FIGS. 2a and 2 b. FIG. 2a shows a borehole 201, the shoulder beds 202 above and below the interval, the invaded portion 203 of the interval, and the uninvaded portion 204 of the interval. FIG. 2b shows the shoulder bed effect 210 and the combined borehole effect and invasion effect 211 relative to the 1D radial tool response 209 and the 2D dual laterolog tool response 208. Preferably, this calculation is only done at the defined bed-boundaries 205, local maximum 206 and minimum points 207 in the intervals, instead of calculating at every sampling point in the dual laterolog measurements 208. This procedure decreases the size of the computations considerably.
  • After the earth model has been updated by the procedure in steps ([0037] 105) through (107), the process then returns to step 103 to recalculate the dual laterolog 2D tool response with the new earth model parameters.
  • If the results of the check for substantial matching in [0038] step 104 are satisfactory, then the process continues with step 108. In step 108, the results of the dual laterolog tool inversion calculated in step 103 are post-processed. Post-inversion processing is preferred to avoid non-physical and non-logical inversion results which may be yielded from the poor sensitivity of the dual laterolog tool to changes of the invasion depth, di, and invasion zone resistivity, Rxo, in some cases. The present invention preferably includes, but is not restricted to, the following post-processing steps. First, if
  • di(j)final ≦di cutoff,
  • then [0039]
  • di(j)final≡0,
  • where the user defines the cutoff invasion depth, di[0040] cutoff, based on the tool's maximum depth of investigation under the given mud resistivity and worst borehole condition of the well. Second, if the ratio of the Rxo vs. the final inverted Rt for the jth interval is less than 1.1 or greater than 0.9, i.e.,
  • 1.1≧R xo(j)/R t(j)final≧0.9,
  • then set [0041]
  • di(j)final≡0.
  • The depth of invasion is no longer physically meaningful when invasion zone resistivity Rxo is approximately equal to final true resistivity Rt. [0042]
  • Finally, in [0043] step 109, the parameters comprising the earth model are output.
  • Mathematical Foundations [0044]
  • The mathematical basis for the method of the present invention will be discussed in more detail. The N horizontal layers of the earth model are represented by the vector {right arrow over (X)}[0045] j∈{dbj, dij, Rmj, Rxoj, Rtj . . . }1≦j≦N, where dbj, dij, Rmj, Rxoj, Rtj are the borehole diameter, depth of invasion, mud resistivity, invasion zone resistivity, and uninvaded formation true resistivity of the jth layer of the earth model, respectively. The well logging measurement may be expressed as the product of the mapping operator G and earth model vector {right arrow over (X)} to form a tool response {right arrow over (R)} in vector space:
  • G{right arrow over (X)}={right arrow over (R)}  (1)
  • The operator G may be a nonlinear as function of {right arrow over (X)}, and EQ. (1) can be solved numerically for {right arrow over (R)} using finite element method forward modeling code. [0046]
  • If observed field data (e.g., resistivity measurements from a field dual laterolog logging tool) is expressed as {right arrow over (D)}, the difference between the computed and field log response may be represented by a χ[0047] 2 merit function:
  • χ2 =∥{right arrow over (R)}−{right arrow over (D)}∥  (2)
  • where ∥ ∥ represent the L[0048] 2norm of the misfit vector. Minimization of equation (2) can be accomplished by solving the following normal equation:
  • H T H{right arrow over (ΔX)}=H T({right arrow over (R)}−{right arrow over (D)})  (3)
  • where {right arrow over (ΔX)} is the increment step in modifying {right arrow over (X)} in the iteration process. H is the Jacobian matrix and an element of H is a set of partial derivatives: [0049]
  • ∂{right arrow over (R)}/∂{right arrow over (X)} j =∂{right arrow over (R)}/∂db j +∂{right arrow over (R)}/∂di j +∂{right arrow over (R)}/∂Rm j +∂{right arrow over (R)}/∂Rxo j +∂{right arrow over (R)}/∂Rt j+ . . .   (4)
  • Computing H and solving for {right arrow over (ΔX)} costs the most computation time in the inversion process (about 99%). Furthermore, it is difficult to solve equation (3) when the number of dimensions of H is large. [0050]
  • In order to reduce the dimensions of H and the computation time, the earth formation model is divided into m intervals whose parameters are represented by {right arrow over (X)}=({right arrow over (X)}[0051] 1 ′, {right arrow over (X)}2 ′, . . . {right arrow over (X)}i−1 ′, {right arrow over (X)}i ′,{right arrow over (X)}i+1, . . . {right arrow over (X)}m ′). Let {right arrow over (R)}i ′ represents the predicted logging tool responses in the {right arrow over (X)}i ′ interval when {right arrow over (X)}i ′ exists independently. Then
  • G{right arrow over (X)}i ′={right arrow over (R)}i ′.  (5)
  • The difference between {right arrow over (R)}[0052] i and {right arrow over (R)}i ′ is defined as the shoulder bed effect {right arrow over (S)}i ′ over the {right arrow over (X)}1 ′, {right arrow over (X)}2 ′, . . . , {right arrow over (X)}x−1 ′, {right arrow over (X)}i ′, {right arrow over (X)}i+1 ′, . . . , {right arrow over (X)}m ′ division. Thus,
  • {right arrow over (S)}i ′={right arrow over (R)}i−{right arrow over (R)}i ′ (1≦i≦m).  (6)
  • Then the χ[0053] 2 merit function of ith interval with the consideration of this shoulder bed effect can be expressed as:
  • x i ′2=|G{right arrow over (X)}i ′+{right arrow over (S)}i ′−Di| (1≦i≦m).  (7)
  • The normal equation for the i[0054] th interval is:
  • Hi ′THi ′{right arrow over (ΔX)}i ′=Hi ′T({right arrow over (R)}i−{right arrow over (D)}i) (1≦i≦m)tm (8)
  • The dimension of the Jacobian matrices H[0055] i ′ (1≦i≦m) is less than the H in equation (3), and the total computation will decrease drastically by dividing {right arrow over (X)} into {right arrow over (H)}i ′.
  • In the cases of deep invasion or thin beds, H[0056] i ′ (or H) is nearly singular. This is because the tool response corresponding to the changes of the formation resistivity or invasion radius has poor sensitivity. To overcome these deficiencies, a damping parameter λ can be introduced. This is described in Levenberg, K., 1944, “A method for the solution of certain nonlinear problems in least square, Quarterly of Applied Mathematics, 2, p164-168, and in Marquardt, D., W., 1963, “An algorithm for least square estimate of non-linear parameters, Journal of the Society of Industrial and Applied Mathematics, 11, p 431-441. Then equations (3) and (8), respectively, can be modified as:
  • (HTH+λI){right arrow over (ΔX)}=HT ({right arrow over (R)}−{right arrow over (D)})  (9)
  • (Hi ′THi ′+λI){right arrow over (ΔX)}i ′=Hi ′T({right arrow over (R)}i−{right arrow over (D)}i) (1≦i≦m).  (10)
  • where I is unit matrix. [0057]
  • When λ is small, equation (9) or equation (10) reverts to equation (3) or (8), respectively, and can be solved by a Newton-Gauss elimination method. When λ is large, the matrix (H[0058] TH+λI){right arrow over (ΔX)} (or (H′TH′+λI){right arrow over (ΔX)}′) is nearly diagonal and can be solved by the gradient method. In the present invention, the singular value decomposition (SVD) method is chosen to solve the above normal equations for the computational benefits, as described in Lines, L. R., and Treitel, S., 1984, “A review of least-square inversion and its application to geophysical problems”, Geophysical Prospecting, Vol., 32, p. 159-186. In the inversion loop, the initial value of λ is set to be 0.001. At the end of each iteration, the cost function is computed and compared with the one in previous step. If it becomes smaller than the one in the previous step, λ is decreased by a factor of 10. Otherwise, λ is increased by a factor of 10. When the χ2merit function is smaller than a pre-set value, a new set of parameters {right arrow over (X)}i ′ will be derived.
  • When m=N, i.e., the number of sub-divided intervals is the same as the number of layers in the earth model, the above described method will be most efficient because the dimension of H[0059] i′ reaches its minimum. In such a case, the tool response solution becomes a 1D radial solution in every sub-interval, i.e., the parameters {right arrow over (X)}i ′ are 1D in the radial direction.
  • To simultaneously invert both deep (LLD) and shallow (LLS) responses of the dual laterolog tool, the Jacobian matrix H[0060] i′ described by equation (10) in the present invention may have four elements in maximum at a given depth and bed:
  • H i′=(∂(LLD)/∂Rt,∂(LLS)/∂Rt,∂(LLD)/∂di,∂(LLS)/∂di).  (11)
  • When no invasion is defined in the pre-processing step, the above Jacobian matrix will be further reduced down to 2 elements: [0061]
  • H i′=(∂(LLD)/∂Rt,∂(LLS)/∂Rt).  (12)
  • The 1D radial tool response from the look-up table and the small Jacobian matrix by equation (11) or (12) drastically reduce computation time needed in the preliminary iteration during the inversion process. The difference between the 1D radial response and the initial 2D response is attributed to the shoulder bed effect. [0062]
  • When solving for equation (9) in the nonlinear least square optimization process, the 2D dual laterolog tool responses are calculated only at the defined bed boundaries, local minimum and maximum points of the interval in the preliminary iterations. This step for reducing both matrix size and computation time in the present invention is schematically illustrated in FIGS. 2[0063] a and 2 b.
    • EXAMPLES
  • Benchmark model inversion is useful to evaluate how accurate the inversion method of the present invention can be. Using benchmark models, the formation parameters (such as Rt) are given in the model and hence, are known and the inversion results can be compared with the “true” value. [0064]
  • The accuracy of the 1D radial inversion was tested. An invaded formation and borehole model with infinite bed thickness (1D radial model) is shown in FIG. 3. FIG. 3 shows a borehole [0065] 301 with diameter db 302, surrounded by an invaded zone 303 with a diameter of invasion di 304, sitting in an uninvaded bed 305. The LLD and LLS “log data” are computed under various Rt models using the 2D finite element method forward modeling code and used as “measured data” in the inversion. Under given Rxo, Rm, and db conditions, Rt and di are inverted from LLD and LLS laterolog responses. In the inversion process, the look-up table and interpolation technique is used to approximate the forward modeling process. Four cases were tested and the results are shown in Table 1. The results show excellent agreement between the given Rt and di in the model and the inverted Rt and di from LLD and LLS responses with the look-up table. The relative error is less than ±0.1% for Rt and ±2% for di, which is much more accurate than the laterolog tool measurement accuracy (±5% in resistivity) defined in the dual laterolog tool specification. These 1D radial cases demonstrate that the look-up table and interpolation method is an excellent approximation of dual laterolog responses in a thick formation. Inversion with the look-up table and interpolation is thousands of times faster than a strict 2D tool response forward modeling-based inversion.
  • The speed of the 2D inversion was tested. FIG. 4 shows a borehole [0066] 401 with diameter db 402, surrounded by an invaded formation 403 with diameter of invasion di 404 and height h 406, inside an uninvaded portion 405 and sandwiched by two infinitely thick shoulder beds 407. Given Rxo=5 ohm-m, Rs=1 ohm-m, Rm=0.1 ohm-m, db=8 inch, and center bed resistivity Rt=20 ohm-m, and the center bed thickness is varied from one to three meters. A full 2D inversion, described by equation (9), is performed to compare with the present invention, described by equation (10). The results are shown in Table 2. The resulting Rt and di from the two inversion methods are fairly similar, but the computation time required is quite different. For this 3 beds case, the full 2D inversion takes more than 5 minutes on the SUN SPARC-20 workstation, while the 2D inversion of the present invention needs less than 30 seconds for the same cases. Thus, the present invention is about ten times faster than the full 2D inversion method. Thus, the present invention can be hundreds of times faster than the full 2D inversion when inversion is performed over hundreds of beds.
  • A multi-beds inversion was tested. FIG. 5[0067] a shows a 7-layer formation benchmark with the 3 resistive and invaded beds is tested with the 2D inversion code of the present invention. In particular, FIG. 5a shows a borehole 501 with diameter db 502, surrounded by 3 invaded formations 503 of diameter of invasion di 504 and height h 506, inside an uninvaded portion 505, 2 uninvaded formations 507 of height h 508 and sandwiched by two infinitely thick shoulder beds 507. FIGS. 5b and 5 c show that the inverted formation resistivity Rt 510 and invasion depth di 512, respectively, agree well with the original formation resistivity Rt 509 and invasion depth di 511 set in the model. The iteration is stopped when the relative difference between the computed LLD and LLS responses in the inversion process and the field dual laterolog log is less than 5% in a given bed, which is dual laterolog tool measurement accuracy. It can be seen that the relative difference between the inverted formation resistivity Rt 510 and the original formation resistivity Rt 509 set in the model is also less than 5% in this example.
  • The above benchmark inversion examples demonstrate that the 2D dual laterolog inversion method of the present invention is a fast and solid inversion technique with sufficient accuracy. The advantages of the present invention include, but are not restricted to, computationally efficient dual laterolog resistivity inversion, more accurate hydrocarbon pore volume estimation analysis, better quantitative formation evaluation, more accurate petrophysical input (such as hydrocarbon pore volume) into reserve assessments, and improved identification of potentially by-passed hydrocarbon accumulations. [0068]
  • It should be understood that the invention is not to be unduly limited to the foregoing which has been set forth for illustrative purposes. Various modifications and alternatives will be apparent to those skilled in the art without departing from the true scope of the invention, as defined in the following claims. [0069]

Claims (6)

What is claimed is:
1] A method for 2D inversion of true resistivity from dual laterolog tool measurements, using a pre-calculated look-up table, comprising the steps of:
deriving an initial layered earth model;
dividing the earth model into intervals;
calculating a 2D tool response in each interval from the earth model;
checking whether the 2D tool response and the tool measurements substantially match; and
iterating the following steps until the match is satisfactory:
deriving a 1D radial tool response in each interval using the pre-calculated look-up table;
calculating shoulder bed effects in each interval by subtracting the 1D radial tool response from the 2D tool response; and
applying a non-linear least square optimization at boundaries of the intervals and local maximum and minimum values in the intervals to update the earth model.
2] The method of claim 1, wherein the step of deriving an initial earth model comprises the further step of:
pre-processing the initial earth model to honor geological and geophysical constraints.
3] The method of claim 1, wherein the number of selected intervals equals the number of layers in the derived earth model.
4] The method of claim 1, wherein the step of calculating the 2D tool response comprises the step of:
applying 2D finite element method forward modeling computer software.
5] The method of claim 1, wherein the step of deriving a 1D radial tool response comprises the step of:
interpolating values from the pre-calculated look-up table.
6] The method of claim 1, wherein the step of updating the earth model comprises the further step of:
post-processing the updated earth model.
US09/961,590 2000-10-03 2001-09-24 Method for 2D inversion of dual laterolog measurements Expired - Lifetime US6430509B1 (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US09/961,590 US6430509B1 (en) 2000-10-03 2001-09-24 Method for 2D inversion of dual laterolog measurements

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US23774900P 2000-10-03 2000-10-03
US09/961,590 US6430509B1 (en) 2000-10-03 2001-09-24 Method for 2D inversion of dual laterolog measurements

Publications (2)

Publication Number Publication Date
US20020040274A1 true US20020040274A1 (en) 2002-04-04
US6430509B1 US6430509B1 (en) 2002-08-06

Family

ID=22894999

Family Applications (1)

Application Number Title Priority Date Filing Date
US09/961,590 Expired - Lifetime US6430509B1 (en) 2000-10-03 2001-09-24 Method for 2D inversion of dual laterolog measurements

Country Status (8)

Country Link
US (1) US6430509B1 (en)
EP (1) EP1332381B1 (en)
AU (2) AU2001294682B2 (en)
CA (1) CA2423567C (en)
DE (1) DE60126311T2 (en)
MY (1) MY128395A (en)
NO (1) NO20031495L (en)
WO (1) WO2002029442A1 (en)

Cited By (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006067539A1 (en) * 2004-12-20 2006-06-29 Schlumberger Technology B.V. Computation of sensitivity for resistivity measurements
US7324898B2 (en) 2005-03-09 2008-01-29 Baker Hughes Incorporated System and method for determining a more accurate resistivity model of a geological formation using time-lapse well logging data
US20080133139A1 (en) * 2005-03-09 2008-06-05 Baker Hughes Incorporated System and method for determining a more accurate resistivity model of a geological formation using time-lapse well logging data
US20090082969A1 (en) * 2007-09-20 2009-03-26 Baker Hughes Incorporated Adaptive Borehole Corrections Accounting for Eccentricity for Array Laterlogs
US20110041051A1 (en) * 2000-05-26 2011-02-17 Libredigital, Inc. Method and system for replacing content in a digital version of a printed paper
US20110087459A1 (en) * 2009-10-09 2011-04-14 Zazovsky Alexander F Cleanup prediction and monitoring
CN102116869A (en) * 2011-02-12 2011-07-06 中国石油大学(华东) High-precision prestack domain least square migration seismic imaging technology
US20120194351A1 (en) * 2011-01-28 2012-08-02 Raymond E. Floyd Downhole Sensor MODBUS Data Emulator
EP2943817A4 (en) * 2014-04-11 2016-03-09 Halliburton Energy Services Inc Estimating subsurface formation and invasion properties
US20160327676A1 (en) * 2014-01-22 2016-11-10 Halliburton Energy Services, Inc. Cross-coupling compensation via complex-plane based extrapolation of frequency dependent measurements
NO339912B1 (en) * 2005-04-18 2017-02-13 Schlumberger Technology Bv Elimination of shoulder layer effects
EP3234657A4 (en) * 2014-12-19 2018-10-24 Baker Hughes, A Ge Company, Llc Hybrid image of earth formation based on transient electromagnetc measurements
US20180364390A1 (en) * 2017-06-16 2018-12-20 Pgs Geophysical As Electromagnetic Data Inversion
US10488547B2 (en) 2014-04-11 2019-11-26 Halliburton Energy Services, Inc. Estimating subsurface formation and invasion properties
CN111581791A (en) * 2020-04-22 2020-08-25 中国海洋石油集团有限公司 Resistivity forward simulation method and device and computer readable storage medium
CN112147706A (en) * 2019-06-26 2020-12-29 中国石油化工股份有限公司 Gravel cave double-laterolog response calculation method and system
US10914859B2 (en) 2015-10-28 2021-02-09 Baker Hughes Holdings Llc Real-time true resistivity estimation for logging-while-drilling tools
US10928540B2 (en) * 2015-06-26 2021-02-23 Halliburton Energy Services, Inc. Systems and methods for characterizing materials external of a casing
US11367248B2 (en) * 2018-03-06 2022-06-21 Halliburton Energy Services, Inc. Formation resistivity evaluation system
US20230228901A1 (en) * 2022-01-19 2023-07-20 Halliburton Energy Services, Inc. Correlating true vertical depths for a measured depth

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7363159B2 (en) * 2002-02-28 2008-04-22 Pathfinder Energy Services, Inc. Method of determining resistivity and/or dielectric values of an earth formation as a function of position within the earth formation
GB2449497A (en) * 2007-05-25 2008-11-26 Statoil Asa Method and apparatus for processing electromagnetic response data
DE112015005897T5 (en) 2015-01-07 2017-10-12 Halliburton Energy Services, Inc. Functional Earth Model Parameterization for Resistance Inversion
CN107305600A (en) * 2016-04-21 2017-10-31 新疆维吾尔自治区煤炭科学研究所 Least square method resistivity three-dimensional approximate inversion technology
WO2019203791A1 (en) 2018-04-16 2019-10-24 Halliburton Energy Services, Inc. Deconvolution-based enhancement of apparent resistivity and bed boundary identification in borehole resistivity imaging
WO2020040742A1 (en) 2018-08-21 2020-02-27 Halliburton Energy Services, Inc. Multi-step inversion using electromagnetic measurements

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4484278A (en) * 1981-12-22 1984-11-20 Schlumberger Technology Corporation Well logging: utilizing superposition of step-profile responses of logging tools to improve logs
US5355088A (en) * 1991-04-16 1994-10-11 Schlumberger Technology Corporation Method and apparatus for determining parameters of a transition zone of a formation traversed by a wellbore and generating a more accurate output record medium
US5675147A (en) * 1996-01-22 1997-10-07 Schlumberger Technology Corporation System and method of petrophysical formation evaluation in heterogeneous formations
US5781436A (en) * 1996-07-26 1998-07-14 Western Atlas International, Inc. Method and apparatus for transverse electromagnetic induction well logging
US6047240A (en) * 1998-01-16 2000-04-04 Schlumberger Technology Corporation Method and apparatus for evaluating the resistivity of invaded formations at high apparent dip angle

Cited By (28)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110041051A1 (en) * 2000-05-26 2011-02-17 Libredigital, Inc. Method and system for replacing content in a digital version of a printed paper
US20080097732A1 (en) * 2004-12-20 2008-04-24 Jean-Marc Donadille Computation of Sensitivity for Resistivity Measurements
US7756641B2 (en) 2004-12-20 2010-07-13 Schlumberger Technology Corporation Computation of sensitivity for resistivity measurements
WO2006067539A1 (en) * 2004-12-20 2006-06-29 Schlumberger Technology B.V. Computation of sensitivity for resistivity measurements
US7324898B2 (en) 2005-03-09 2008-01-29 Baker Hughes Incorporated System and method for determining a more accurate resistivity model of a geological formation using time-lapse well logging data
US20080133139A1 (en) * 2005-03-09 2008-06-05 Baker Hughes Incorporated System and method for determining a more accurate resistivity model of a geological formation using time-lapse well logging data
US8116979B2 (en) 2005-03-09 2012-02-14 Baker Hughes Incorporated System and method for determining a more accurate resistivity model of a geological formation using time-lapse well logging data
NO339912B1 (en) * 2005-04-18 2017-02-13 Schlumberger Technology Bv Elimination of shoulder layer effects
US20090082969A1 (en) * 2007-09-20 2009-03-26 Baker Hughes Incorporated Adaptive Borehole Corrections Accounting for Eccentricity for Array Laterlogs
US8775084B2 (en) * 2007-09-20 2014-07-08 Baker Hughes Incorporated Adaptive borehole corrections accounting for eccentricity for array laterologs
US9121263B2 (en) * 2009-10-09 2015-09-01 Schlumberger Technology Corporation Cleanup prediction and monitoring
US20110087459A1 (en) * 2009-10-09 2011-04-14 Zazovsky Alexander F Cleanup prediction and monitoring
US8781807B2 (en) * 2011-01-28 2014-07-15 Raymond E. Floyd Downhole sensor MODBUS data emulator
US20120194351A1 (en) * 2011-01-28 2012-08-02 Raymond E. Floyd Downhole Sensor MODBUS Data Emulator
CN102116869A (en) * 2011-02-12 2011-07-06 中国石油大学(华东) High-precision prestack domain least square migration seismic imaging technology
US10317562B2 (en) * 2014-01-22 2019-06-11 Halliburton Energy Services, Inc. Cross-coupling compensation via complex-plane based extrapolation of frequency dependent measurements
US20160327676A1 (en) * 2014-01-22 2016-11-10 Halliburton Energy Services, Inc. Cross-coupling compensation via complex-plane based extrapolation of frequency dependent measurements
EP2943817A4 (en) * 2014-04-11 2016-03-09 Halliburton Energy Services Inc Estimating subsurface formation and invasion properties
US10488547B2 (en) 2014-04-11 2019-11-26 Halliburton Energy Services, Inc. Estimating subsurface formation and invasion properties
EP3234657A4 (en) * 2014-12-19 2018-10-24 Baker Hughes, A Ge Company, Llc Hybrid image of earth formation based on transient electromagnetc measurements
US10928540B2 (en) * 2015-06-26 2021-02-23 Halliburton Energy Services, Inc. Systems and methods for characterizing materials external of a casing
US10914859B2 (en) 2015-10-28 2021-02-09 Baker Hughes Holdings Llc Real-time true resistivity estimation for logging-while-drilling tools
US20180364390A1 (en) * 2017-06-16 2018-12-20 Pgs Geophysical As Electromagnetic Data Inversion
US10871590B2 (en) * 2017-06-16 2020-12-22 Pgs Geophysical As Electromagnetic data inversion
US11367248B2 (en) * 2018-03-06 2022-06-21 Halliburton Energy Services, Inc. Formation resistivity evaluation system
CN112147706A (en) * 2019-06-26 2020-12-29 中国石油化工股份有限公司 Gravel cave double-laterolog response calculation method and system
CN111581791A (en) * 2020-04-22 2020-08-25 中国海洋石油集团有限公司 Resistivity forward simulation method and device and computer readable storage medium
US20230228901A1 (en) * 2022-01-19 2023-07-20 Halliburton Energy Services, Inc. Correlating true vertical depths for a measured depth

Also Published As

Publication number Publication date
DE60126311D1 (en) 2007-03-15
CA2423567A1 (en) 2002-04-11
EP1332381B1 (en) 2007-01-24
NO20031495D0 (en) 2003-04-02
US6430509B1 (en) 2002-08-06
WO2002029442A1 (en) 2002-04-11
NO20031495L (en) 2003-06-02
MY128395A (en) 2007-01-31
AU9468201A (en) 2002-04-15
EP1332381A1 (en) 2003-08-06
CA2423567C (en) 2012-05-01
AU2001294682B2 (en) 2006-04-13
DE60126311T2 (en) 2007-06-28
EP1332381A4 (en) 2005-06-15

Similar Documents

Publication Publication Date Title
US6430509B1 (en) Method for 2D inversion of dual laterolog measurements
AU2001294682A1 (en) Method for 2D inversion of dual laterolog measurements
US8095345B2 (en) Stochastic inversion of geophysical data for estimating earth model parameters
EP2810101B1 (en) Improving efficiency of pixel-based inversion algorithms
US7565244B2 (en) Method and system for removing effects of conductive casings and wellbore and surface heterogeneity in electromagnetic imaging surveys
US7991553B2 (en) Method and system for removing effects of conductive casings and wellbore and surface heterogeneity in electromagnetic imaging surveys
Eaton 3D ELECTROMAGNETIC INVERSION USING INTEGRAL EQUATIONS1
Li et al. Fast modeling and practical inversion of laterolog-type downhole resistivity measurements
Rovetta et al. Petrophysical joint inversion of multi-geophysical attributes and measurements for reservoir characterization
Dyes Inversion of induction log data by the method of maximum entropy
Domenzain et al. Efficient inversion of 2.5 D electrical resistivity data using the discrete adjoint method
US7756641B2 (en) Computation of sensitivity for resistivity measurements
Butler Research Note: The mean sensitivity depth of the electrical resistivity method
Li et al. 2D cross-hole electromagnetic inversion algorithms based on regularization algorithms
Yin et al. A practical 2D dual laterolog (DLL) inversion method and its impacts on HPV estimation
Cosmo et al. Fast deconvolution of laterologs by direct inverse filtering
Yan et al. 2D Pixel Based Inversion of Ultra-deep Electromagnetic Logging Data for Look-ahead Applications
Costa et al. 1D inversion of electromagnetic logs in the characterization of carbonate formations from Santos basin–Brazil
Yang et al. Forward and inversion of azimuthal lateral resistivity logs
Israil Delineation of layer boundaries from smooth models obtained from the geoelectrical sounding data inversion
Passey et al. AAPG Archie Series, No. 1, Chapter 10: Modeling Log Responses in Thinly Bedded Reservoirs
Terpolilli et al. A Method for Computing Fast and Accurate Analytical Gradient: the Key for Routinely Performing Efficient 2D and 3D Resistivity Inversions
Bordakov et al. Robust Log Sharpening with Unknown Tool Response Function
Yang Modeling and inversion of azimuthal lateral resistivity logging
Reddig et al. Imaging of magnetotelluric data

Legal Events

Date Code Title Description
AS Assignment

Owner name: EXXONMOBIL UPSTREAM RESEARCH COMPANY, TEXAS

Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNORS:YIN, HEZHU;WANG, HANMING;REEL/FRAME:012210/0587

Effective date: 20010920

STCF Information on status: patent grant

Free format text: PATENTED CASE

FPAY Fee payment

Year of fee payment: 4

FPAY Fee payment

Year of fee payment: 8

FPAY Fee payment

Year of fee payment: 12