System and Method of Petrophysical Formation Evaluation in Heterogeneous Formations

This invention relates to the field of wireline logging and logging-while-drilling (LWD) and interpretation, and in particular, relates to a system and method for interpreting geological formations comparing a multiplicity of actual log measurements with log measurements predicted using computational tool response models. A specific embodiment relates to a system and method for determining bed boundaries and associated bedding formation properties in a horizontally stratified geological formation.

An oil & gas company, after drilling a well borehole, is particularly interested in the subsurface geological formations surrounding the well borehole. The oil & gas company would like to understand the geological formation to determine if the formation contains hydrocarbon resources (oil and/or gas), the amount and producibility of hydrocarbon the formation contains, and the best options for completing the well, as well as anticipate any problems that might be encountered during completion and production. A significant aid in this evaluation is the use of wireline logging and/or LWD measurements of the formation surrounding the borehole (sometimes referred to herein collectively as "logs" or "log measurements"). Typically, one or more logging tools are lowered into the borehole and the tool readings or "measurement logs" are recorded as the tools traverse the borehole. These measurement logs are used to infer the desired formation properties.

In evaluating the hydrocarbon production potential of a subsurface formation, the formation under study is described in terms of a set of "petrophysical properties." The petrophysical properties include (1) the lithology or the rock type, e.g. amount of sand, shale, limestone, or more detailed mineralogical description, (2) the porosity or fraction of the rock that is void or pore space, and (3) the fluid saturations or fractions of the pore space occupied by oil, water and gas. Wireline logging tools do not directly measure petrophysical properties, they measure "log properties", for example, bulk density, electrical resistivity and acoustic velocity. Log properties are related to petrophysical properties via a set of mathematical or statistical relations called "petrophysical mixing relations." In formation evaluation, a set of wireline logs are interpreted or processed, using knowledge of petrophysical mixing relations, to estimate the formation petrophysical properties such as porosity, saturation and lithology.

As no single physical measurement provides all the information needed to evaluate hydrocarbon production potential, many different types of tools have been developed for providing information regarding the petrophysical properties of a formation. See, e.g. U.S. Pat. Nos. 4,780,679, 4,722,220, 4,703,277, 4,646,240, 4,597,290, 4,531,188, 4,939,648, 5,355,088, and 5,350,925 (incorporated by reference). Typically, in a logging operation, several different logging tools are combined and used simultaneously to obtain an integrated set of measurements which are used to obtain the desired formation information. In this context, different tools may, in fact, be used to obtain information about the same set of formation properties using different techniques, or different tools may be used to obtain information about different formation properties. Due to differences in physical measurement modalities and sensor configurations, different logging tools have different volumes or zones of investigation, hence different measurement resolutions. For example, a sonic tool may look at the geological formation around the borehole extending for several meters away from the borehole, but its vertical resolution along the borehole is about a meter. In contrast, the Formation MicroScanner.TM. tool has a vertical resolution of several millimeters, but its depth of investigation is correspondingly small.

Subsurface formations are generally heterogeneous, so that porosity, saturation and lithology vary with position. Examples of heterogeneity include the presence of geological beds and the result of migration or invasion of a bank of borehole fluid into the formation. Because each logging tool has a nonzero volume of investigation, more than one heterogeneous region (for example, multiple beds) may lie within the volume of investigation of a tool. Consequently, as different logging tools have different volumes of investigation, in processing the logs from multiple tools during formation evaluation, the effect of formation heterogeneity must be corrected to avoid erroneous results. Otherwise, for example, the petrophysical evaluation of one region may be interfered and distorted by the presence of another region of heterogeneity lying within the larger volume of investigation of one tool, for example, while not in that of the others.

Prior to performing formation evaluation, the effects of formation heterogeneity are partially accounted for using a series of "corrections". These compensate for certain simplified types of heterogeneity--the presence of the borehole ("borehole correction"), the presence of a single layer ("shoulder correction"), the presence of borehole fluid which has migrated into the formation ("invasion correction") and the presence of a non-horizontal bedding orientation ("dip correction"). Despite the fact that these effects are mutually coupled in the instrumental response of the measurement, each of these corrections assumes that the only type of heterogeneity present is the one being corrected. These corrections comprise a sequence of independent preprocessing steps performed prior to formation evaluation, and during formation evaluation it is implicitly assumed that the formation is homogeneous.

Computational models exist for precisely simulating how a wireline logging or LWD tool responds to a specified heterogeneous formation, a process sometimes referred to as "forward-modeling." However, such forward models are typically of individual tool responses and are primarily used in research for designing tools. Use of forward model tool responses for formation evaluation has been limited or nonexistent. Generally, it would be a significant advance if formation evaluation used tool response forward models to better accommodate multiple logging tools with multiple volumes of investigation for improved petrophysical evaluation in heterogeneous formations.

Another objective and challenge in multiple tool formation evaluation is the ability to depth align or depth match different tools traversing a common depth in the borehole at different times. In the wireline case, due to the elasticity of the cable supporting downhole tools and a nonconstant drag friction along the borehole wall, tools record downhole information at true depths that are different from the "cable depth" reported by the uphole winch. In the LWD case, similar difficulties exist in tool depth alignment due, for example, to the temperature-dependent elasticity of the drilling pipe and length variations in drill pipe segments. For both situations, fine-scale depth alignment of logs from different tools using conventional algorithms fails because all tools have a nonzero resolution and blur the formation with a tool response signature which is characteristic to each tool. It follows that fine-scale registration of logs from multiple tools requires explicit use of tool response models.

A commonly encountered heterogeneous formation is a stack of thin beds. Formation evaluation in thin beds gives rise to the particular problem of identifying the bed or layer geometry (number of beds, their locations, sizes and orientation). Usually, bed boundaries are placed at the estimated inflection points of the acquired logs. Although it is computationally efficient, this approach of estimating bed boundary positions has several important drawbacks: the tool responses are not taken into account, and the inflection point estimates obtained from different types of measurements (e.g. acoustic and nuclear) of the same bedding structure will typically differ and cannot be easily combined. This problem is especially difficult when the logging tools used have a vertical resolution greater than the bed thicknesses, e.g. a vertical tool resolution of 2-3 meters with beds less than 1 meter thick.

In this disclosure, all patent references are incorporated by reference and all non-patent references are incorporated by reference for background.

The limitations of current formation evaluation techniques outlined above are addressed by the system and method of the present invention which uses a set of measured wireline and/or LWD logs and knowledge of the associated tool response models to estimate the heterogeneous distribution of earth formation properties. For the important particular case of a stack of horizontal beds and linear tool response models, the bed boundaries and the formation properties inside each bed are estimated in a computationally efficient manner.

The predicted logs are computed using a separate forward response model for each tool, but with a common input formation description. The formation description typically includes both a geometrical description of the heterogeneity (how the formation is divided into regions such as beds or invasion zones) as characterized by the set of parameters g and a description of the formation properties within each region (associated with the rock and fluid types present), as characterized by the set of parameters p. There is a hierarchy of geometrical descriptions, depending on the physical complexity of the heterogeneity. As shown in FIG. 1, this hierarchy is based on the formation's geometrical complexity. The "one-dimensional" model corresponds to a formation having a spatial variation only in the depth direction, being composed of a set of beds with otherwise homogeneous characteristics. In this case, the geometrical description g involves the depth locations of the bed boundaries. In the "two-dimensional" case, borehole fluid invasion is introduced into this model, providing spatial variation in the radial direction from the borehole. The geometrical parameters are the bedding boundaries augmented by the invasion boundary location in each bed. Correspondingly, the "three-dimensional" case, involves the full geometric generality, e.g. introducing tilt or dip into the bedding structure, with the attendant parametric representations of the bedding orientation added to the parameter set.

The formation properties in each region of the chosen geometrical model may be represented in a hierarchy of properties models, based on the elemental level of composition of the constitutive rock/fluid components as shown in FIG. 2. At the "macroscopic" level, the model involves a set of so-called log properties, p.sub.l which are those parameters associated with the direct log response, such as bulk density, electrical resistivity and acoustic velocity. A more detailed model involves using a set of the petrophysical properties p.sub.p such as porosity, fluid saturations and lithology, which are related to the log properties through petrophysical mixing laws. For example, using Archie's equation, the electrical resistivity is a function of the water saturation and porosity. At the "atomic" level, the elemental or mineralogical volumes such as the volume quartz, volume dolomite, volume oil etc. comprise the set of model parameters, p.sub.v. They are, in turn, related to the petrophysical properties through a forward model or mixing law. Different representations offer different advantages in solving the inverse problem, with the dimensionality of the models along with the statistical independence of the parameter set increasing from the log properties level through the hierarchy.

Generally speaking, the method estimates a formation description by comparing the actual log measurements with a corresponding set of predicted logs. The predicted logs generated from the chosen formation descriptions (g & p) and the tool response models (including mixing laws if appropriate) are compared with the actual measured logs, and a difference error between the predicted and actual logs is determined. An object function or mismatch criterion is defined, which is a prescribed function of this error. It is used to provide a measure of "goodness" of the predictive logs to represent the physical reality of the measured logs. This object function is combined with a constraint function (e.g. upper and lower bounds on the formation description parameters), and an optimization problem is formulated where the formation parameters are chosen to minimize the mismatch criterion, subject to the imposed constraints. Depending on the form of the object function, constraints imposed, and variational knowledge of the forward model, this minimization involves an iterative process whereby the formation parameter estimates are successively refined until the mismatch criterion is reduced to an acceptable level. Because, in general, the computed logs depend nonlinearly on both of the formation geometric and property parameters, the best estimates are obtained using an iterative nonlinear optimization procedure. See P. E. Gill, W. Murray, and M. H. Wright, Practical Optimization, San Diego, Academic Press, 1981.

The preferred method begins by specifying an initial set of formation geometry and property parameters, for example, layer boundaries from the position of log inflection points, or invasion front positions and properties using logs from tools having multiple depths of investigation. Beginning with this initial model, the procedure iteratively estimates the formation geometry and property parameters by comparing a number of actual log measurements with a number of predicted log measurements. The predicted logs are computed using a response model for each tool based on the initial estimate of the formation geometry and properties. The initial predicted logs are compared with the actual log measurements to generate a mismatch criterion or object function that depends on the difference error. A new estimate of the formation geometry and properties is determined by the optimization algorithm employed, so as to reduce the mismatch criterion. This process is iteratively repeated, the final formation representation being chosen as the "best" estimate of the formation geometry and properties when the value of the object function is sufficiently small. Further, the parameters may also be constrained to satisfy auxiliary conditions, for example, to lie within a predetermined range of values, or for a set of values to sum to unity. In this case, a constrained Gauss Newton optimization procedure is used.

In a preferred form, the system and method of the present invention is used in layered formations with one-dimensional linear tool response models, including convolutional models, and a weighted least-squares mismatch criterion. With linear tool response models, if the geometrical bed boundary parameters are specified, the optimum layer property values can be obtained directly without a search process. This leads to an efficient optimization procedure where only the set of geometrical parameters are estimated iteratively, and at each step of the iteration the optimum formation properties are identified noniteratively by matrix inversion. Compared to the more general optimization problem, this leads to searches in lower dimensional parameter spaces, with attendant improvements in computational efficiency and faster solutions.

With layered formations, fine-resolution "ground truth" information reflecting bed boundaries is often available, such as from core or log data from the Formation MicroScanner.TM. or Formation MicroImager.TM. (FMS, FMI) imaging services by Schlumberger Technology Corporation. See, M. P. Ekstrom, C. A. Dahan, M. Y. Chen, P. M. Lloyd, and D. J. Rossi, Formation Imaging with Micro Electrical Scanning Arrays, The Log Analyst, 28(3), 294-306, 1987. In this context, fine-resolution implies resolutions which are much better than those of the logging tools used to measure the formation parameters (e.g. vertical resolutions on the centimeter or millimeter scale, rather than meter scale). Such fine-scale information can be introduced into the algorithm in several ways. These data may be introduced into the algorithm as additional weighted, "measured logs." The weight attributed to this type of log must reflect (1) the influence of the disparity of resolution, as the fine-scale logs will resolve many more beds than the conventional logs and not all of these beds may be included in the formation model, and (2) the quality of the fine-scale data, which is sensitive to bad borehole conditions and/or mud conditions. Alternatively, the fine-scale data may be used to establish a bedded or zoned formation model, for example, by using an image edge detection or texture segmentation algorithm to process FMS or FMI imagery. Some or all of the detected beds or zones may be selected for inclusion in the formation model, and this is especially useful in establishing an initial model.

A corollary to the present method is the ability to evaluate the noise sensitivity and contribution of each log measurement to the formation description. For example, the effect of measurement noise can be incorporated into one or more tool response models and the sensitivity of the estimated formation description to the measurement noise determined. This, of course, can help define the quality and accuracy of the estimated formation description and the contribution of each tool to the quality of the estimated formation description. It may also be used to assess the value of additional measurements not yet acquired in the design of future logging programs. It may additionally be used to assess resolution-noise tradeoffs in modeling thinly layered formations. Specifically, the iterative computational method in this invention can be used several times beginning with different initial models, ranging from a thickly bedded model having few layers to a thinly bedded model having a large number of layers (e.g. established from core or an imaging logging tool). By performing noise sensitivity analysis on the resulting estimated formations, the relationship may be established between modeled bed thickness and accuracy of the property estimates. This information serves as a guide to select the minimum acceptable bed thickness in the formation model.

Another corollary to the present method is the ability to depth align or depth match logs from different tools traversing a common depth in the borehole at different times. One method of this invention makes explicit use of tool response models to account for nonzero tool resolution, and introduces auxiliary "depth shift" variables that are estimated during inversion and used for fine-scale registration of logs from multiple tools.

Tool response models are developed and used primarily in research and engineering in support of new tool development. A diversity of tool models exist, and even for a specific tool, several models may be available to compute the response to different types of heterogeneity. For example, several models might be available to compute the response of a tool in formations of varying geometric complexity, such as a vertical stack of layers (one-dimensional model), or a stack of layers with fluid invasion into the various layers (two-dimensional model), or a stack of layers that are dipping at some angle with respect to the borehole (three-dimensional model). Among the simplest tool response models are the linear models, but nonlinear models are required to handle mixing laws and more complicated heterogeneity. This invention makes use of multiple tool response models for each tool, and invokes the tool response model that is most appropriate for the modeling task being performed and the geometric complexity of the heterogeneity. In a preferred embodiment of the algorithm, a hierarchical formation evaluation approach is used, consistent with Tukey's principle of parsimony, which begins with the simplest tool response models, and invokes more complex models as they are needed to obtain data consistency and a suitably small mismatch criterion.

One embodiment invoking increasingly complex models begins with a formation model that is a stack of layers, and assumes linear tool response models to obtain a computationally efficient solution for the bed positions and properties. Because the formation model is assumed to be a stack of layers, the logs are preprocessed as necessary to correct for the presence of the borehole, fluid invasion and formation dip. Although the algorithm is designed under the assumption of a layered formation description, it functions well even with moderate amounts of fluid invasion and layer dip angle. If fluid invasion is present, in one embodiment it can be handled to a first-order approximation by assigning a separate layer parameter to each tool that is sensitive to invasion. For example, each tool with a different depth of investigation may provide its own estimate of layer resistivity, and the different resistivity estimates from different tools may be compared to infer fluid invasion.

If invasion is present beyond some limit, say with highly permeable sands between shale layers, a more complicated two-dimensional layered-invaded formation model is invoked. Not all of the layers will require a description with fluid invasion. For example, impermeable shales do not require additional invasion parameters, and a layer with deep invasion extending beyond the depth of investigation of the available tools does not require additional invasion parameters. Only for those layers which require additional invasion parameters to reconcile the predicted and measured logs, 3, 4 or 5 invasion parameters are introduced: bed positions, invasion radius, and properties in the invaded and uninvaded zones. For computational efficiency with certain tool models, the bed boundary estimation problem may be performed separately from the estimation of the radial invasion boundary and properties.

Alternatively, if bedding dip is present beyond some angular limit, say in highly deviated and horizontal well applications, a more complicated three dimensional formation model of a dipping stack of layers is invoked. In this case, the relative dip angle of the formation to the borehole is assumed to be known from external information such as from the Dipmeter, FMS or FMI imaging services from Schlumberger Technology Corporation, and the dipping tool response models are used for formation evaluation.

The system for interpreting a geological formation of the present invention ideally comprises a specially programmed computer having memory storing a set of actual measurement logs from an oil or gas well borehole in a geological formation of interest. A set of tool response models is also stored in memory where each response model predicts a measurement log based on an input formation description. The specially programmed computer includes a mechanism for inputting or changing the formation description. The processor of this specially programmed computer generates a predicted log based on an input formation description and said tool response model retrieved from memory. The processor operates to compare the predicted logs and the corresponding measurement logs and generates a difference error between the predicted logs and the measurement logs. A display communicates the formation description, and optionally the difference error between the predicted logs and measurement logs, to the interpreter. In the preferred embodiment, the processor operates to iteratively change the formation description, computing a mismatch criterion which includes the contribution of the difference errors and arriving at a formation description corresponding to a minimal mismatch criterion.

FIGS. 1A-1C is a vertical sectional view of a collection of geometrical models of the subsurface formation showing an increasing order of geometric complexity.

FIG. 2 is a chart displaying the tree structure of formation property models for an individual subsurface layer, showing log properties, petrophysical properties and volumetric properties.

FIG. 3 is a block diagram illustrating the method of the present invention with the formation properties given in terms of "log" parameters.

FIG. 4 is a block diagram illustrating the method of the present invention with the formation properties given in terms of "petrophysical" parameters.

FIG. 5 is a block diagram illustrating the method of the present invention with the formation properties given in terms of "volumetric" parameters.

FIG. 6 is an illustration of an interpretation workstation for evaluating geological formation in accordance with the present invention;

FIG. 7 is a side view of a wireline logging operation of a typical geological formation;

FIGS. 8(A)-8(C) are graphs showing actual wireline measurements of a test well;

FIGS. 9(A)-9(C) are graphs depicting predicted logs of the test well used in FIGS. 8(A)-8(C); and

FIG. 10 is a graph showing the total weighted fit error as a function of the number of beds of the test well used in FIGS. 8 and 9.

FIG. 11 is a flowchart illustrating a MATLAB routine for optimizing bed boundary positions.

I. General Approach

Turning to the drawings, FIG. 7 schematically illustrates a typical wireline logging operation where one or more tools 10 are assembled into a string and lowered into an oilwell 12. Tools 10 may be conveyed on a wireline cable 14, or on the end of drill pipe during logging while drilling (LWD), or in some special types of wells such as horizontal wells, the tools may be conveyed by coiled tubing. The oilwell 12 is bored through a geological formation 16 which in FIG. 7 is a layered formation containing a number of beds or layers 18 of different petrophysical properties. As previously mentioned, a number of different types of tools 10 are available for looking at different aspects of the formation. Measurements of the formation from the tools 10 are conveyed by the cable 14, or with LWD using a retrievable downhole memory or using mud column pulsing, to a computer acquisition system 20, such as residing in the truck shown in FIG. 7. These measurements are typically called "logs" and usually comprise a tool measurement recorded versus distance along the oilwell 12.

FIG. 6 shows a workstation 22 for interpreting the logs or measurements obtained from the acquisition system of FIG. 7. It should be appreciated that the interpretation function of FIG. 6 may actually occur in the acquisition truck 20 during well logging operation or shortly thereafter. Interpretation can occur long after the logging operation in the offices of the oil and gas company.

Log analysis is sometimes called a "geophysical inverse problem". Prediction of the response of a sensor, (e.g. an oilwell logging tool) from a known physical system (e.g. a geological formation) is called a forward modeling problem. Traditional log analysis or formation evaluation defines the physical properties of a geological formation from tool measurements, and is called an inverse problem. While a great deal of research and understanding is available relating to inverse problems, the commonly used processes of estimating a formation description from multiple tool measurements has not incorporated forward tool modeling nor used any process to insure a consistent formation model for multiple tool observations.

Many types of forward modeling codes have been developed and are available for commercial electromagnetic logging tools (e.g. linearized Born approximation for induction tools, numerical mode matching for induction tools in layered, invaded formations, 2D and 3D finite-element models for Laterolog resistivity tools), nuclear logging tools (e.g. 1D and 2D sensitivity functions based on Monte Carlo analysis), and sonic logging tools (e.g. ray tracing, real axis integration and finite-difference). Such forward modeling codes have as input a description of the geological formation and they compute a predicted logging tool response.

Turning to FIG. 3, the general operation of the method of the present invention is schematically illustrated. In FIG. 3, a formation description 30 is used as input into one or more tool response forward models 32. Forward models 32 may be constructed for practically any type of logging tool, with nuclear, electromagnetic, and acoustic forward models 32 being shown in FIG. 3. The models chosen should match the tools used to provide the set of measured logs. The tool response forward models 32 output the predicted response from the sensor, that is, a set of predicted logs 34.

The predicted logs 34 are compared with the actual measured logs 24 in a comparison step 36 as denoted in FIG. 3. The measured logs 24 are those obtained from the acquisition process in FIG. 7 and are illustrated in FIG. 6. An object function or mismatch criterion 38 is defined which is a function of the difference error obtained by comparing the predicted logs 34 and the measured logs 24. It is used to indicate the accuracy and consistency of the formation description 30. In the method of the present invention, the object function 38 is combined with constraint information 38c based on the formation description 30. A constrained optimization algorithm 38a is used to compute a correction or adjustment 38b to the formation description 30. The formation description at 30, consisting of both geometry and property parameters, is modified by the adjustment 38b so that the object function 38 is reduced, and the process of comparing the predicted logs 34 with the measured logs 24 iteratively repeated. When the object function 38 attains a minimum, or is below a certain predetermined threshold, the current formation description 30 is used as the estimate of the geological formation 16.

The object function 38 is general and may, for example, be the weighted sum of several factors. The most important contributing factor is the difference error between the predicted and measured logs. The difference error from each logging tool may be weighted independently to account for different logging units (e.g. density in grams/cc versus velocity in ft/sec), the relative importance of the log, and the quality and noise levels in the logs.

Constraint information 38c about the formation model 30 may also be incorporated as input to the constrained optimization 38a. It is included in order to favor formation estimates with desired properties. One constraint is boundedness, for example, the density and acoustic velocity are known to be positive. Another constraint is a statistical relationship between properties, for example, a specified statistical correlation between acoustic velocity and density. Other formation model constraints such as smoothness may be invoked using "regularization" terms, see, e.g. A. N. Tikhonov and V. Y. Arsenin, Solutions of ill-Posed Problems, New York, Wiley, 1977.

In more detail, the formation description 30 of FIG. 3 depends on a set of numerical parameters g which represent the geometry (e.g. bed boundary positions), and parameters p.sub.l which represent the formation log properties, such as the density, resistivity and acoustic velocity within each bed. The formation description 30 is denoted as F.sub.l (x, g, p.sub.l), where F.sub.l is the set of formation log properties at the spatial point x, which depends on the geometry parameters g and log property parameters p.sub.l. Note that the properties are not required to be constant-valued in each geometrical region; for example, there may be a ramp in property value within a bed.

F.sub.l (x, g, p.sub.l), the formation description 30 of FIG. 3, is based on formation log properties p.sub.l, but may be expressed in terms of alternative properties as shown in FIG. 2. One alternative is F.sub.p (x, g, p.sub.p), the formation description 42b of FIG. 4, which is based on formation petrophysical properties p.sub.p such as porosity, fluid saturations and lithology. Another form is F.sub.v (x, g, p.sub.v), the formation description 42a of FIG. 5, which is based on formation volumetric parameters p.sub.v, such as volume clay, volume quartz, etc.

As shown in FIG. 2, the formation property parameters p.sub.l, p.sub.p and p.sub.v in these representations are related to one another through statistical or deterministic petrophysical mixing relations 44, which may be considered as an additional nonlinear forward model. One example of such a nonlinear relation is the Archie saturation equation describing the electrical resistivity of a water-saturated clean sandstone: ##EQU1## where R.sub.l is the rock resistivity, .phi. is the porosity, m is the cementation exponent, R.sub.w, is the resistivity of the pore-filling brine water, and s.sub.w is the water saturation. The Archie relation is typical of many petrophysical mixing relations, which include additional auxiliary variables, such as coefficients, exponents, and connate water properties. These variables are sometimes assumed to be known a priori, but may also be considered as unknown variables which are estimated during the inversion procedure. Some of these petrophysical auxiliary variables may be constant within a geological zone containing many beds, and external constraint information such as FMS and FMI imagery may be used to divide the formation into zones of constant values for these variables.

A particular benefit to the method of FIG. 3 is the ability to characterize the stability, robustness, and quality of the resulting formation description. That is, known statistical techniques can be applied to evaluate the sensitivity of the method to measurement noise. For example, if the predicted nuclear log is changed 5% in a particular bed, the corresponding change in the estimated formation description 30 in FIG. 3 can be assessed. The availability of the tool response models 32 permits a mechanism for propagating such measurement noise to evaluate the impact on the estimated formation description. Using this method with the system of FIG. 3, the contribution of each measurement to the overall formation description may be assessed, which can lead to the development of a logic or methodology for deciding when and how to include diverse measurements in a given interpretation. For example, if the set of measured logs 24 does not include an acoustic log, but analysis using the tool response models 32 predicts that the presence of an acoustic log would contribute significantly to the quality of the formation description 30, then the acquisition of an acoustic log might be called for, particularly for future wells in the same field.

While logging tools have been discussed in general as a method of predicting a petrophysical formation description (see FIG. 7), it is understood that other types of tools and forward models may be similarly integrated using this approach. For example, fluid pressure and fluid flow measurements such as those available from well testing tools may be incorporated into the method of FIG. 3. With the use of well testing tools, the formation description can be extended from the lithology-porosity-saturation formation description generally used in formation evaluation to a permeability formation description. Measurement while drilling (MWD) and logging while drilling (LWD) tools which acquire data during the drilling process can be incorporated into the method of FIG. 3, as well as point measurements such as those obtained from core samples. Deep reading measurements such as seismic, gravity and electromagnetic (single-well, cross-well, well-surface and surface-surface) could be incorporated into the method of FIG. 3. As used herein, "tool" may incorporate logging tools, well test tools, core sample measurements, LWD/MWD and deep reading measurements.

II. Computational Efficiency in Layered Formations

An important application of petrophysical evaluation using tool response models is in a layered formation. Here, the problem is one of locating bed boundary positions and determining the petrophysical properties within each layer. In the case where the tool response models 32 in FIG. 3 are linear, this problem may be solved in a computationally efficient way.

The initial or starting formation description 30 is an initial estimate of bed boundary positions, obtained for example, from the inflection points of the measured logs 24.

Given M logs, the following three-step algorithm is used:

S0: Initialize the number of beds: N=N.sub.0, and their positions: z.sub.k,i,1.ltoreq.k.ltoreq.N-1,1.ltoreq.i<M, and fix z.sub.0.i to the depth of the beginning of the zone and z.sub.N,i to the depth at the end of the zone. Also initialize measured property values within each bed: .mu..sub.k,i which is the value of property p.sub.i (z) for z.sub.k-1,i .ltoreq.z<z.sub.k,i.

S1: Optimize bed boundary positions, z.sub.k,i, 1.ltoreq.k.ltoreq.N-1 and property values .mu..sub.k,i, 1.ltoreq.k.ltoreq.N-1 so that the object function (depending on the fit error between the measured and the predicted logs) is minimum.

S2: If the overall fit is acceptable STOP, otherwise introduce one more bed boundary at the depth where the error is the largest and GO TO S1.

The optimization step, S1, is preferably a constrained optimization to minimize an object function given by a weighted fit error as follows: ##EQU2## where z.sub.i and .mu..sub.i are sets of variables for the bed boundary positions and log property values for the i.sup.th log. The predicted i.sup.th log l.sub.i is computed by using the following 1-D linear measurement model:

where g.sub.i is the one-dimensional linearized tool response model for the i.sup.th tool, * denotes convolution operation and p.sub.i is the layered property model defined as:

The minimization of this object function with respect to the layer positions and the property values is performed in two steps: first choose the layer positions z.sub.k,i and then find optimum property values .mu..sub.k,i within these layers. Given the layer positions z.sub.k,i, the predicted logs are linear functions of the layer properties .mu..sub.k,i, so the optimum layer properties .mu..sub.k,i can be directly and noniteratively identified by solving an appropriate system of linear equations. Therefore, the only variables in this optimization problem which must be estimated iteratively are the layer positions. For linear tool response models, then, the above algorithm may be rewritten as follows:

S0: Initialize number of beds: N=N.sub.0 and their positions: z.sub.k,i, 1.ltoreq.k.ltoreq.N-1,1.ltoreq.i.ltoreq.M, and fix z.sub.0,i to the depth at the beginning of the zone and z.sub.N,i to the depth at the end of the zone.

S1: Optimize bed boundary positions, z.sub.k,i, 1.ltoreq.k.ltoreq.N-1 so that the following fit error between the measured and the predicted logs is minimized: ##EQU3## where the predicted i.sup.th log l.sub.i is computed by using the 1-D linear measurement model in equation (1) and p.sub.i is the layered property model in equation (2), where the .mu..sub.k,i are set equal to the optimum layer property values given the layer boundary positions.

S2: If the overall fit is acceptable STOP, otherwise introduce one more bed boundary at a depth where the error is the largest and GO TO S1.

Preferably the bed boundaries are constrained, otherwise the above process may produce different layered formation models for each measured log property. The following set of constraints are used, but not necessarily all used at the same time.

1. Minimum allowed bed thickness of T.sub.min :

2. No allowed depth shift across the measurements:

3. Static depth shift across the measurements relative to the first measurement:

with a maximum allowed static depth shift of .DELTA.z.sub.max which is greater than or equal to .DELTA.z.sub.i.

4. Dynamic depth shift across the measurements:

where .DELTA.z.sub.max is the maximum allowed dynamic depth shift at each layer boundary.

Appendix A contains MATLAB computer routines which implement the above algorithm and handle different sets of constraints. The main routine is called "beds". A flowchart for "beds" is included in Appendix A, as well as a list of the important variables in "beds". The "beds" program calls the built-in MATLAB constrained optimization routine called "constr". It also calls several functions which are included in Appendix A, such as:

"fcc", which computes the object function

"gcc", which computes the gradients of the object function

"fmu", which computes the optimal bed properties given the boundary positions

"fmcontr", which computes the predicted log by convolution

"fraugment", which combines two matrices into one matrix.

The "beds" routine and its associated functions include the following functionalities:

Multiple measurements of different properties with different sampling rates can be used for bed boundary identification. The various measurements and their sample positions are loaded in at the beginning of the "beds" routine.

Each measurement covers approximately the same zone, in order that the "beds" routine can effectively combine information from multiple logs traversing a common section of the formation. The common depth interval endpoints are stored in S(1) and S(2) in the "beds" routine.

Each measurement is treated as discrete samples of an underlying continuous signal which is found by third order polynomial interpolation between any two samples. This is achieved in the "fcc", "gcc" and "fmu" functions by computing a weighted combination of multiple predicted logs computed by repeat calls to the "fmcontr" function.

Linear 1-D vertical tool response functions are used in generating predicted logs. The convolution of the property variation along the well with the tool response function is performed in the "fmcontr" function.

User specifies the minimum bed thickness Tmin in the final description. This value is used at the end of the "beds" routine to insure that new bed boundaries are only introduced into beds thicker than 2Tmin. It is also used as an explicit constraint on the bed thickness at the end of the "fcc" function, where S(3)=Tmin.

Main routine calls a constrained iterative nonlinear optimization package for solution of optimal positions for a given number of beds. The nonlinear optimization program is "constr" and is called by the "beds" program.

Constrained optimization package has access to the function "fcc" which computes the weighted fit error and its partial derivatives with respect to bed boundary positions.

Constrained optimization package has access to the function "gcc" which computes the constraints and their partial derivatives with respect to the bed boundary positions.

A function called "fmcontr" exists for predicting the log in a layered formation by convolving the formation property with the tool response model.

A quantitative relative importance of each bed in a layered description is computed near the end of the "beds" routine to provide an indication of the significance of each bed boundary.

Additional Matlab routines to perform computations similar to those in the "beds" routine are possible. Such routines include the additional capability to estimate static depth shifts between different logs, and to estimate minor perturbations of the sample positions (so-called depth matching problem) to compensate for imperfect sampling alignment due to cable stretch and erratic tool movement in the borehole.

The approach implemented in the "beds" Matlab routine has been validated on real logs obtained by running the Schlumberger Maxis Express tool in a man-made formation consisting of a stack of rock layers, in which case the true rock layer boundary positions and layer properties are known. Three measurement channels were chosen for analysis: the first energy window count rates of the gamma-gamma back scattering detector, the second energy window count rates of the gamma-gamma short spacing detector, and the high-resolution resistivity measurement.

FIGS. 8(A), (B), and (C) show the three measurement channels as a function of depth. In this test case, the bed boundaries were identified using the "beds" routine of Appendix A. The corresponding predicted logs for the measurements of FIGS. 8(A)-(C) are shown in FIGS. 9(A)-(C). As can be seen, there is good agreement between the measurement channels (FIG. 8) and the predicted logs (FIG. 9) for a formation consisting of 11 beds. FIG. 10 depicts the weighted fit error versus the number of modeled beds. As can be seen, the weighted fit error decreases significantly in going from 8 beds to 11 beds, and it then levels off. The following table shows a comparison between actual bed positions and predicted bed positions for a final model with 11 beds.