US 4,398,273 AGrant
Method for Interpreting Seismic Records to Yield Indications of Gas/Oil in an Earth Formation
Issue Date:1983-08-09
•25 Claims
•7 Drawing Sheets
Abstract
The present invention indicates that acoustic characteristics associated with gas/oil-containing strata of an earth formation, including reflectivity coefficients can be normalized (and favorably compared) with similar characteristics calculated and displayed by means of a machine-implemented data processing method in which well logging and geologic data are fed thereto to calculate such characteristics without the need for shear-wave velocities. In more detail, in accordance with the invention brine-saturated bulk and shear moduli, (i.e., Kw* and Gw*) of the strata of interest can be predicted as a function of, say, brine-saturated P-wave modulus (Pw*) alone (independent of shear-wave velocity). In that way, resulting acoustic values including seismic velocities and amplitudes (also, reflectivities) as a function of saturation can ultimately be provided. Such values, when compared to actual field-generated characteristics, are surprisingly accurate predictors of the amount of gas/oil saturation in the zone of interest. The method has particular applicability to designating gas-sand zones within formations of interest.
Metadata
Assignee
- Chevron Research Company
Inventors
- Don D. Thompson
- Robert J. S. Brown
- Richard J. Runge
Application Information
Application Number:US 3225097
Filing Date:1981-11-18
Priority Date:1979-10-02
Art Unit:222
Classifications
IPC:
G01V 132G01V 134
Field of Search:
36736473;75421
Patent Drawings (7 sheets)
Description
DISCLAIMER
While the Abstract, supra, has been carefully written, the purpose of such statements is to provide a non-legal description of the contents of this application as a searching, scanning and classification aid for technical persons. Accordingly, all hereinbefore-presented statements are not intended to be used in understanding or otherwise comprehending the principles of the invention hereinafter described in detail and are not, more particularly, to be used in interpreting or in any way limiting the scope or fair interpretation of the claims appended hereto.
FIELD OF THE INVENTION
The present invention pertains to the art of seismic prospecting for petroleum reservoirs, and more particularly to the art of interpreting high-intensity reflection amplitude anomalies observed on acoustic record traces into diagnostic indicators of the presence and amounts of gas/oil in the underlying subsurface strata.
BACKGROUND OF THE INVENTION
For several decades, seismic prospecting for petroleum has involved the creation of acoustic disturbances above, upon, or just below the surface of the earth, using explosives, air guns, or large mechanical vibrators. Resulting acoustic waves propagate downwardly in the earth, and are partially reflected back toward the surface when acoustic impedance changes within the earth are encountered. A change from one rock type to another, for example, may be accompanied by an acoustic impedance change, so that the reflectivity of a particular layer depends on the velocity and density content between that layer and the layer which overlies it, say according to the formula ##EQU1## where AR is the amplitude from the reflected signal and Ai is the amplitude of the incident signal; V.sub.1 is the compressional velocity of the wave in the overlying medium 1; V.sub.2 is the compressional velocity in the medium layer below the contact line; d.sub.1 is the density of the overlying medium 1; and d.sub.2 is the density of the underlying medium.
Today's seismic interpretors have made good use of ultra-high amplitude anomalies in seismic traces to infer the presence of natural gas in situ. So-called "bright-spot" analysis has been used to good advantage to indicate several large gas reservoirs in the world, especially in the Gulf Coast of the United States. Such analysis is now rather common in the oil industry, but it is not without its critics, especially in the area of predicting gas saturation based on the characteristics of the amplitude anomalies of the traces alone.
The present invention improves the ability of the seismologist to correctly interpret gas as well as gas/oil content of a formation normalized to (and comparable with) a series of patterned acoustic characteristics associated with zones of similar mineralogy and determinable gas or gas/oil saturations.
OBJECT OF THE INVENTION
An object of the invention is the provision of a novel method of correctly predicting gas and/or gas/oil saturations of hydrocarbon-bearing structures using at least in part high-intensity anomalies provided by traces, including seismic traces.
SUMMARY OF THE INVENTION
The present invention indicates that acoustic characteristics associated with gas/oil-containing strata of an earth formation can be normalized (and favorably compared) with similar synthetic values calculated and displayed by means of a machine-implemented data processing method in which well logging and geologic data are fed thereto to calculate such values.
In more detail, in accordance with the invention brine-saturated bulk and shear moduli, (i.e., Kw* and Gw*) of the strata of interest can be predicted, say, as a function of brine-saturated P-wave modulus (Pw*) (independent of shear-wave velocity), thereafter a surprisingly accurate estimation of bulk moduli of the fluid-filled pore spaces in the zone of interest can occur. In that way, various synthetic values (say velocities, amplitudes, etc.) as a function of depth and of gas/oil saturation can be provided. Ultimately, such values when compared to actual field-generated values have accurate prognostic capabilities as to the amount of gas/oil saturation in the zone of interest.
DESCRIPTION OF THE DRAWINGS
Further features of the invention will become more apparent upon consideration of the following detailed description of the invention when taken in connection with the accompanying drawings, wherein:
FIG. 1 is a plan view of a grid of centerpoints produced in the field by the systematic positioning and energization of an array of seismic sources and detectors whereby a series of locational traces associated with individual centerpoints between respective source-detector pairs are ultimately generated;
FIG. 2 is a model of typical reflecting horizons within an earth formation that can be associated with the characteristics of the locational traces of FIG. 1;
FIG. 3 is a plot of percentage of gas saturation vs. compressional velocity as a function of change in pressure (depth) of a particular subsurface material;
FIG. 4 is a diagram of process steps for carrying out the method of the present invention using, in at least part, programmed digital computing system;
FIGS. 5 and 6 are schematic diagrams of elements within the digital computing system of FIG. 4;
FIGS. 7 and 8 illustrate the diagnostic capability of the method of the present invention in predicting the presence of and amounts pertaining to gas strata in an actual field example.
PREFERRED EMBODIMENTS OF THE INVENTION
Before discussion of an embodiment of the invention with an actual field environment, a brief description of the mathematical and theoretical concepts behind the discovery may prove beneficial and are presented below.
Firstly, it may be of interest to indicate lithology limitations associated with the present invention. Anomalies associated with gas or gas/oil sands over shale rock cap are one example where the method of the present invention offers surprising predictive capabilities. Another example relates to gas-saturated limestone or gas/oil-saturated limestone over shale. Secondly, it should be pointed out that certain relationships of various petrophysical parameters are also of great importance, viz., bulk and shear moduli of a rock formation as a function of the P-wave modulus in the zone of interest.
Now, in more detail, attention should be directed to the Figures, particularly FIG. 1, in which the field-collection phase of the present invention is shown. Note that, inter alia, FIG. 1 illustrates in some detail how seismic multi-centerpoint data are collected in the field, and derivations of terms of interest in this application occur: e.g., the term "centerpoint" is a geographical location located midway between a series of sources S.sub.1,S.sub.2 . . . Sn of a geophysical field system 9 and a set of detectors D.sub.1,D.sub.2 . . . Dm at a datum horizon near the earth's surface. The centerpoints are designated C.sub.1,C.sub.2 . . . Cp in the Figure, and are associated with a trace derived by placement of a source at that centerpoint location followed immediately by relocating a detector thereat.
I.e., if the sources S.sub.1 . . . Sn are excited in sequence at the source locations indicated, traces received at the different detector locations shown can be related to common centerpoints therebetween. If such traces are summed, a gather or group of traces is formed. I.e., if the reflecting interface is a flat sand horizon containing gas, or oil/gas fractions, the depth point where reflection occurs will define a vertical line which passes through the centerpoint of interest. Applying static and dynamic corrections to the field traces is equivalent (under the above facts) to placing the individual sources S.sub.1,S.sub.2 . . . Sn at the centerpoint in sequence followed by replacement with the detectors D.sub.1 . . . Dm of interest at the same locations. If the traces associated with a common centerpoint are summed, a series of enhanced traces, sometimes called CDPS traces, is provided.
FIG. 2 illustrates reflection phenomena of a three-layer model typical of a young, shallow geologic section 10, such as found in the Gulf Coast, illustrating how reflection phenomena associated with traces of the field system 9 of FIG. 1 can be related to the presence of gas alone.
Section 10 includes a gas sand 11 embedded in a shale stratum 12. A 20% velocity reduction can occur at interface 13, say from 10,000'/sec to 8000'/sec, and a 10% density reduction from 2.40 g/cc to 2.16 g/cc, can likewise happen.
Actual P-wave reflection coefficient Apr of gas sand 11 can be calculated by Equation (1) below: also, P-wave transmission displacement amplitude coefficient Apt can similarly be related in accordance with Equation (2) below. ##EQU2##
Equations (1) and (2) are, of course, the two basic equations of wave travel in an earth formation and are for isotropic media with the layer index being i=1 for the incident medium and i=2 for the underlying medium. Equations (3) through (17) simple define intermediate variables.
As an example of calculations associated therewith, if .theta.-=0.degree. (normal incidence), the P-wave reflection coefficient Apr is equal to about -0.16 and +0.16, respectively, assuming a Poisson's ratio of 0.1 for sand 11 and 0.4 for the shale 12 of FIG. 2. Such values are larger, of course, than reflections from barren strata.
And how various gas fractions affect elastic moduli and seismic velocity, which can be of monumental interest to seismologists, is depicted in FIG. 3.
FIG. 3 also illustrates to some degree the complexity that is involved in relating various petrological parameters one to the other to the elastic moduli in the presence of gas. E.g., FIG. 3 shows that very small gas saturations cause nearly as much reduction in moduli of compression as do larger gas saturations. In highly porous sandstone a very small gas saturation may lead to a substantial reduction in compressional wave velocities, especially at low differential pressures. However, an increased gas saturation leads to a density decrease which may cause a velocity increase, partly offsetting the initial velocity decrease. Also, effects of gas saturation are greatest for rocks of highly compressable elastic framework: Thus, in most low-porosity carbonates of high differential pressure, the effect of gas saturation is minimal.
Briefly, the present invention shows that appropriate expressions can, in fact, be provided for calculating the elastic moduli of competent and incompetent strata based on reasonable assumptions, followed by a determination of synthetic characteristics as a function of oil/gas saturations. Comparison techniques with actual field-collected data then can be utilized. Result: accurate estimations of gas/oil saturations for the zones of interest occur. The expressions for providing such data are derived, in part, from surprising laboratory measurements to be described below. Also, well logging and other geologic data of the zone of interest are involved to some degree.
GENERAL PRINCIPLES
Expressions shown below relate densities, velocities and elastic moduli. First, in accordance with theory:
and
Here, .rho.b is the bulk density (including contributions of pore fluids; Vp is the P-wave velocity, or plane-compressional-wave velocity; Vp is the shear-wave velocity; p* is the P-wave elastic modulus; K* is the compressional modulus (reciprocal of the compressibility) and G* is the shear modulus, or modulus of rigidity.
General Definitions
The rocks in these equations are assumed to be macroscopically homogeneous and isotropic. In fact, macroscopic isotropy and homogeneity are assumed to exist when petrological parameters are represented by only density and two other parameters, such as Vp and Vs or else G* and K*. The asterisk superscript (*) applied to an elastic modulus indicates reference to porous rock with fluid (any mixture of brine, oil and gas) in the pores. By measurement of .rho.b, Vp and Vs, other parameters such as P*, K*, G*, Young's modulus or Poisson's ratio can be imparted.
Two additional macroscopic elastic moduli, the solid-material elastic modulus Km or Gm and the solid-frame modulus Ka, can be used.
If the solid material of a porous rock is reassembled in such a manner that there is no pore space, the resulting elastic moduli depend on the manner in which the material is assembled. However, upper and lower limits to the moduli can be computed. For such hypothetical rocks without pore space, upper and lower limits are not widely separated, and the means between these limits will be indicated by subscripts "m". Thus, Km and Gm are solid-material moduli. These quantities are set forth below in Table I.
But if the pore space in a porous rock is all connected, and if the differential pressure in the fluid in the pores is held constant during a small change in the externally applied pressure, the ratio of the pressure change to the relative change in the bulk volume is the solid-frame modulus of compressibility, Ka. This modulus is sometimes referred to as the modulus of compressibility of the rock with the pore system empty. However, many oilfield rocks contain clays which may be altered by actually drying the rocks to empty the pore system. It is also possible that certain surface forces may be different when all water is removed from the pore system of a rock. It is also possible that there are several different possibilities for finding values of Ka.
If all pore space is connected, and if the rock is macroscopically homogeneous and isotropic, a change in the pore pressure cannot lead to macroscopic shear strains, since excluded from the specification of the system is anything which could define direction. Therefore, the compressibility of the fluid cannot influence the macroscopic shear modulus of the rock. Thus, G*=Ga for a rock satisfying the above conditions.
Note, however that the dynamic elastic moduli of rocks usually depend on the static stresses and therefore on the fluid pressure in the pore system. This, in turn, affects the differential pressure, or the pressure supported by the elastic framework of the rock. Furthermore, note the requirement that there be time for the fluid pressure to equalize throughout the pore system. If the pore space is not connected, this cannot happen at all. If there are localized regions of low permeability, there may not be enough time for pressure equalization in the pore system if pressure changes are applied at high frequency. The time required for pressure equalization depends not only on the local permeability of the material but on the microscopic distribution of the pore fluid components. This consideration applies, for example, if the pore fluid is a mixture of brine and gas. The higher the seismic or sonic frequency, the finer must be the mixture of brine and gas.
However, in estimating gas/oil-saturation effects in accordance with the present invention, input parameters can be limited to compressional velocity, density and porosity, as set forth below.
The Effect of Gas Saturation
Gassmann derived an equation (Gassmann, F., "Elasticity of Porous Media" in Vierteljahrschrift der Naturforshenden Gesellschaft in Zurich, Heft I, 1951) and given also in a book by White (J. E. White, "Seismic Waves: Radiation, Transmission and Attenuation", McGraw-Hill, N.Y., 1965), relating K* and Ka. ##EQU3## or the equivalent, ##EQU4## or another equivalent, ##EQU5## or still another equivalent, ##EQU6## where a=Ka/Km; b=K*/Km; g=G*/Gm=Ga/Gm; f=(Km/Kf)-1.
Here .phi. is the fractional porosity, Kf is the compressibility modulus of the fluid in the pores, and the primes (') indicate values for a different pore fluid system. Thus, if Kf'=0, the K*' becomes Ka. If the pore fluid is a multiphase system, such as any mixture of brine, oil and gas, Kf can be calculated by noting that the compressibility of a fluid mixture is the volume-weighted average of the compressibilities of the components. Thus, for a gas-brine mixture, ##EQU7## where Sgas is the fractional gas saturation.
It can be seen from Equations 19 and 20a-20e that, if G*, P*, Km, Kf and .phi. are known, calculations are straightforward. But since G* and Km are not known, simultaneous solutions of equations is a distinct possibility and is set forth in detail below. If K* for two different values of Kf is measured, two equations with two unknowns are provided and both Ka and Km can be determined. Km may be computed with adequate accuracy from estimates of the composition of the solid framework of the rock, as set forth below, in which case measurements of K* and .phi., together with knowledge or estimates of Kf and Km give values of Ka.
Equations 18, 19 and 20a-20e are rigorously true if the rock is macrohomogeneous, macroisotropic, all pore space connected, time is allowed for pressure equalization in the pore fluid, and the material of the solid framework is microhomogeneous and microisotropic. The assumptions of microisotropy and microhomogeneity are not satisfied rigorously for real rocks. The lack of microisotropy is probably well accounted for by the spatial averaging employed in the determination of the "isotropized" values of the compressional moduli of the rock-forming minerals given in Table I, supra. It is highly probable that anisotropic mineral grains are deposited more or less randomly. Even if they are not, the directions of the stress vary greatly in a granular rock, justifying the spatial averaging.
The effects of microhomogeneity of the solid material of the rock may affect the strict validity of Equations 20a-20e supra if one or the components of the solid framework does not carry its share of the load and if that component has a compression modulus substantially different from those of other components. A candidate would be a sand consisting of a quartz framework with a large amount of clay or other materials with higher compressional moduli and attached in such a manner that it does not significantly share the stress.
It should be noted that if shear-wave logs would be available for the basin of interest, the bulk modulus can be calculated from the P- and S-wave transit times in the density of the zone of interest in accordance with Equations 18 and 19, supra. Unfortunately, the shear-wave logs needed for direct measurement of K* are not available for most basins.
Hence, a method is needed for accurately determining bulk and shear moduli of rock formations in such basins, these moduli being needed for use in the Biot-Gassmann relationship, i.e., the Equations 20a-20d supra, to model the influences of pore fluids on the bulk modulus, and therefore on the compression wave velocity. While these fluids may be either brine, liquid, hydrocarbon or gas or any combination thereof, the present invention relates to brine-gas mixtures principally but not exclusively.
In accordance with the present invention, only compressional velocity of the brine-saturated material, its brine-saturated density and porosity measurements are needed to estimate bulk and shear moduli for the zone of interest as functions of pressure, temperature and fluid content (gas saturation).
These equations which are set forth below and above permit an in-situ estimate of K* and G* as defined above, which may be based on measurement of the P-wave modulus of the aggregate (i.e., Pw*) only, and do not require the need for shear-wave velocity. Also, a gas correction term applies only to well-consolidated sandstones, not to sand packs or loosely or poorly consolidated sands, as set forth in detail and discussed below.
A key to the above discovery that Pw* can be used to estimate the elastic moduli of the formation of interest on the basis of empirical expressions of the form
where G* is the shear modulus, and the bulk modulus by the form
where Po1* is the P-wave modulus. The subscript "ol" means "observed in the laboratory". Of course, the subscript "ow" is also operative in this context: in circumstances where P-wave modulus is "observed from well log data", as explained below.
The empirical coefficients for the above Equations 21 and 22 have been found for sands without performing undue experimentation and, moreover, produce correlation coefficients in a range of 0.927 to 0.962.
In general, Equations 21 and 22 are proven excellent predictors. E.g., the above equations have been found to give correlation coefficients of 0.96 and 0.93 respectively at the highest pressure conditions. A standard acoustic measurement system (AMS) can be used to determine densities, porosities and three moduli (P*, K* and G*) as a function of pressure of known materials to predict the coefficients of the above equations.
Note also that values of P-wave modulus estimated from well logs is in accordance with
where .rho.b is the density and Vp is the compressional velocity from such logs, as explained below.
However, mineralogy can also be a determining factor.
The limestone equations have practical validity in areas well represented by that particular sampling used above. Experience has shown that such data are particularly good predictors of velocity and porosity values in the Ardmore and Arcoma basins in the Midcontinent of the United States of America. However, the Biot-Gassmann relationship as set forth in Equation 20b should be used: ##EQU8## where Ka=frame modulus (the bulk modulus with empty pores)
Km=bulk modulus of the solid material
Kf=bulk modulus of the fluid (including the gas) in the pores.
The second term of the right-hand member of Equation 20d is referred to as the fluid correction term, Kfc, so that K*=Ka+Kfc, as previously indicated.
Mathematical Modeling of Gas in Sandstone Formations
Equations 21 and 22 are assumed to be applicable in the area under study, as determined from well logs via equation 35. Without this assumption, calculation of fluid saturation effects on the P-wave modulus is not possible.
Parenthetically, there must be available from well logs +.rho.b, .phi. and Vp, from which the observed P-wave modulus, Pow*, is calculated ("ow" means "observed from well logs") by equation 35, supra. These data must come from a section of the well known to be 100% brine saturated. From these steps, the elastic moduli can be calculated.
FIG. 4 is a flow diagram of the data transformation process of the present invention.
Briefly, in the initial step, predicted shear modulus, Gw*, for the depth (pressure) is calculated from the predictor Equation 21 supra, i.e.,
In step (2), there is the calculation of the predicted bulk modulus, Kw*, for the depth (pressure) from the predictor Equation 22 supra, i.e.,
In step (3), to model gas-saturated formations, there is next calculated Ka via Equations 20a-20d. In these statistical terms, Equations should be modified to read ##EQU9## where Kbr is the bulk modulus of brine. Assume a pressure-independent value of 25 Kbar. The bulk modulus of the solid material, Km, is also assumed to be pressure-independent.
From the several sandstone samples, a means Km of 393 Kbar is noted with a standard deviation of less than 10%. This value can be used as a starting point if one has no other data.
It is physically impossible for Kw* to exceed Km. However, the predicted modulus, Kw*, can exceed the mean value of Km. In that event, either:
1. abort the method for that set of data, or
2. independently determine the true value of Km.
Errors in Ka are apt to be largest for porosities less than about 5%. Above that value, we had no problems in this step of the calculation.
In step (4), calculation of Kgas*, the predicted bulk modulus with gas present, in accordance with: ##EQU10## where Kfg equals the bulk modulus of a liquid-gas mixture. It is equal to
where Sq is the fractional gas saturation and Kg is the bulk modulus of methane. The latter is both temperature- and pressure-dependent.
Estimation of the dynamic bulk modulus of natural gas can be more rigorous. Approximate methods to estimate the fluid modulus, Kf, for use in the Gassmann relations can be used, and are set forth supra.
Attempts to estimate the dynamic bulk modulus of a gas component, Kdgas, from either the isothermal or adibatic perfect gas laws are inadequate for some applications, but results set forth below in Table III can be useful to obtain the speed of sound in methane over a range of pressure and temperature conditions of interest in accordance with:
where .rho.gas is the density and Vp is the compressional velocity.
Methane compressibilities can be obtained from solutions of the Benedict, Webb and Rubin equations of state (see Journal of Petroleum Technology, Vol. 22, 1970, pp 889-895), or from conventional empirical gas equations equivalent thereto.
Table V shows a constituent analysis of gas from a typical South Louisiana well. The gas gravity is 0.6995 with a 44,000 G.O.R. Table V indicates that methane and ethane are the predominant components and, hence, the use of the methane equations presented here seems justified for use in Gulf Coast estimations of reflectivities.
In step (5), equations 21 and 23 are used to obtain Pgas*. Since Gw*=Ggas*, according to the Biot-Gassmann theory,
Pgas* is the predicted P-wave modulus with gas present.
In step 6, to convert the moduli to velocities, use:
and the density equations, assuming a value for Sg,
where .rho..sub.fg =S.sub.g .rho..sub.g +(1-S.sub.g).rho..sub.br, to obtain the predicted value of V.sub.p, i.e., V.sub.p with gas present.
In step 7, seismic characteristics of interest are calculated as a function of different gas saturations and compared with actual values until a pattern match occurs. The amount of gas saturation is then evident.
DESCRIPTION OF EMBODIMENT
Having now established a firm mathematical and theoretical basis for the process of the present invention, perhaps a description for a system for carrying out the invention is in order. In this regard, reference again should be made to FIG. 4, which is a flow diagram for such a system. From the above descriptions, it is evident that displays ultimately will be provided which indicate both presence and amounts of gas/oil-bearing strata within a zone of interest.
It should be pointed out that while theoretical acoustic characteristics can be provided using the expressions and descriptions set forth above, actual acoustic data must be available from the field; such data are initially collected, say, using one of today's conventional techniques, e.g., involving seismic multiple-area coverage. One of the better field-gathering techniques is a conventional form of recording seismic data using what is called common-depth-point gathers (CDP). From such data, characteristics of events associated with gas/oil-bearing strata can be normalized and then compared with the theoretical data predicted by the method involving the present invention, such values being variable as a function of gas/oil saturation, inter alia.
Note further in returning to FIG. 4 that it can be assumed that a section of data has been analyzed for "bright spots" related to, say, gas-bearing strata alone, such events being known by geographic location and/or time-depth basis and the traces have been dynamically and statically corrected.
The steps of FIG. 4 include calculation of synthetic acoustic values of interest, as a function of gas saturation, and comparing these values with actual field values obtained using conventional field-gathering techniques.
Such comparison will include calculations of the theoretical data values as a function of gas saturation; then the pattern of such values as a function of gas saturation can be analyzed so that an interpreter can compare actual field data therewith and therefore estimate the gas saturation in the zone of interest.
Comparison Techniques
The comparison of such values can occur within (or without) the particular computing system to be adopted. If internal comparison techniques are used, the nature of the results (to be equated in a comparison sense) must be taken into account:
(i) for velocity comparisons, an array comparison technique can be used within the system to be adopted in which dimensions of the array include: velocity amplitudes (both synthetic and actual) as a function of geometrical location, depth (or equivalent time) values, and gas/oil saturation.
While either the synthetic or actual values can serve as the primary selector values, it is perhaps better to use the latter to compare with the former. In that way, the array addresses associated with individual velocity values can be incremented in sequence (in one or more dimensions) and the results compared. The best match becomes apparent from a reading or listing of all comparisons made, or from a designation of the best match on a GO-NO GO basis. In this regard, some preliminary processing analysis of the seismic record is in order to establish the active seismic velocity deducible from the record as a function of time or depth. A technique well known in the art to reconcile theory and data, uses an iterative technique in which a stacking velocity for a set of traces (corrected) is associated. Hereafter, the set of corrected traces (all associated with a common centerpoint) is stacked and coherency of the data noted. Then the process is repeated with a different stacking velocity; the most coherent of the stacked data allows a "best" stacking velocity determination to be deduced.
(ii) for amplitude comparisons, an array technique similar to that previously described can be used; but conversely an auxiliary comparitor could also be implemented to function with disk data provided from the computing system to be described hereinafter. In the former, dimensions of the array would include: actual and synthetic amplitude values as a function of geometrical centerpoint location, time (or migrated equivalent depth) values and gas/oil saturation values. Incrementing of the array (in one or more dimensions) allows for a determination of "best" match of the data. In the latter, "best" match statistics can be developed after a synthetic amplitude-versus-time trace has been generated. Of course, such a synthetic trace requires amplitude values to be established, say based on the relationship ##EQU11## where Ar is the amplitude from the reflected signal and Ai is the amplitude of the incident signal; V.sub.1 is the compressional velocity of the wave in an overlying medium 1; V.sub.2 is the compressional velocity d.sub.2 is the density of the underlying medium below the contact reflection line; d.sub.1 is the density of the overlying medium 1. Also, the above data must be migrated to the time dimension as such information is initially associated with established depth values within the earth. Thereafter the original seismic record (usually in a composite trace basis) is cross-correlated with the various synthetic traces (each of the latter being associated with different gas/oil saturation values so that the best match therebetween becomes easily established.
FIG. 5 illustrates particular elements of a computing system for carrying out the steps of FIG. 4, including certain aspects related to the comparison of data.
The System
While many computing systems are available to carry out the process of the present invention, to best illustrate operations a microcomputing system 50 of FIG. 5 is didactically best and is presented in detail below.
System 50 of FIG. 5 can be implemented on hardware provided by many different manufacturers, and for this purpose elements provided by MOS Technology, Norristown, Pa. 19401, may be preferred.
Such a system 50, which can be an MCS-6500 system from the above corporation, can include a CPU 51 controlled by a control unit 52 in conjunction with signals from external devices (one or more), not shown. Two memory units 53 and 54 connect to CPU 51 through bus 55. Program memory unit 53 (ROM) stores instructions for directing the activities of the CPU 51, while data memory 54 (RAM) contains data as data words related to logging and other geologic data provided by field acquisition systems. Since the logging and geologic data banks comprise large amounts of big-bit data, an auxiliary memory unit 56 may also be useful in carrying out the method of the present invention. Inputting and accessing such information is usually via peripheral interfaces 57. To rapidly input or access data, the CPU 51 addresses a particular port (say at input port 58, output ports 59). Such ports are addressed to receive or access information as required, such information being fed say via keyboard 60 or outputted or storable on usual external equipment known in the art, e.g., printer 61, floppy disks 62, paper-tape readers, etc.
The number and type of peripheral devices needed to carry out the method of the invention can generally be estimated very accurately. However, it is important to keep in mind that such estimates must be subject to review after a full analysis of system performance is completed. The seismic-interpreter-designer may find it necessary to use a special-purpose interface part or to redesign the I/O structure if the evaluation of total system performance reveals that the system performance reveals that the system cannot operate at the required speed. Use of special-purpose peripheral interface parts will reduce the number of tasks which must be handled by the CPU 51 and consequently can increase the over-all system speed, but this generally involves additional component cost.
Likewise, the use of a fully vectored interrupt can lead to increased performance at increased cost. The goal of any program must be to meet all the system performance at the minimum possible cost.
After the various peripheral devices in the system have been evaluated to determine the number of inputs and outputs required, the total required by all peripherals can be divided by 16 to determine the number of devices required. This is a good first approximation which will be reevaluated as the system development progresses.
Evaluation of the amount of RAM required by the system is a somewhat more difficult problem than estimation of peripheral devices. This is due primarily to the fact that much of the RAM is required by the system software as working storage, such as storage of immediate results in arithmetic operations. Since the machine-implemented program will probably not be written when these estimates are first attempted, the probability of error in this portion of the estimate may be fairly high.
In addition to working storage, the RAM must provide storage for:
1. The Stack (to be described below);
2. Peripheral input data storage;
3. Peripheral output data storage.
Items 2 and 3, above, can be evaluated quite accurately, since a detailed analysis of the peripheral devices has usually been completed when these estimates are first attempted. In general, a block of RAM must be made available for each peripheral device. The amount of RAM required for each is a function of the type of peripheral device being interfaced and just how the device is to be controlled.
The amount of RAM required by the stack is a function of both the interrupt structure and the system software. As a result, an estimate of this requirement must be based on the system programmer's best estimates of his requirements. This should be combined with an estimate of the required working storage and the peripheral data storage requirements to obtain an estimate of the total system RAM.
Amount of ROM required in a system cannot be determined accurately until the system program is completed. However, by partitioning the system program into definable pieces, an estimate can be made of each task and the total can be obtained of the ROM required by each section.
Most programs consist of easily defined sections such as the software for each peripheral device, arithmetic routines, etc. These are the pieces which should be examined separately to estimate the ROM required by each.
For outputting information, system 50 can be included with printer 61, again accessed through a particular interface 57 at an output port 59, whereby the operations in accordance with the present invention are printable. Of more use as an outputting unit, however, would be a disk unit which would temporarily store the accumulated data after the latter had again been accessed through a different interface 57 at another output port 59. In that way, the information on the disk unit can be used in conjunction with an off-line digital plotter capable of generating proper side-by-side displays of the data for use in conjunction with typical seismic data sections for comparisons and interpretation, as explained below. Such plotters are available in the art. One proprietary model that we are familiar with uses a computer-controlled CRT for optically merging onto photographic paper as a display mechanism the data provided by the present invention. Briefly, in such a plotter, the data are converted to CRT deflection signals; the resulting beam is drawn on the face of the CRT and the optically merged record of the events recorded, say via photographic film. After a predetermined number of side-by-side lines have been drawn, the film is processed in a photography laboratory and hard copies returned to the interpreters for their review.
Of course, the CPU 51 of FIG. 5 is the key to the operation in accordance with the present invention. As a general matter, it controls the functions performed by the other components, including addressing all input and output ports. It also fetches instructions from program memory unit 53, decodes their binary content and executes them. It also references data memory unit 54 (and/or other auxiliary memory units) and other interface ports as required in the execution of the instructions. It also recognizes and responds to certain control signals as set forth below.
FIG. 6 illustrates CPU 51 in more detail.
As shown, the CPU 51 includes an array of registers surrounding but operatively tied to an arithmetic logic unit (ALU) 63 through an internal data bus 64, all under control of a control unit generally indicated at 65. As indicated, the control unit includes clock generator 66, timing control unit 67 and instruction decoder 68. Note that the instruction decoder 68 is also operatively timed to the following units for control purposes: an instruction register 69, processor status register 70 and interrupt logic circuit 71, operating as set forth below. While instruction register 69 and program counters 72 have dedicated uses, an array of other registers, such as accumulator 73 and index registers 77 have more general uses.
The accumulator 73 usually stores one of the seismic operands to be manipulated by the ALU 63. E.g., in the summation of two or more numbers, the instruction may direct the ALU 63 to not only add in sequence the contents of the temporary registers containing predetermined values of logging or geologic data in the accumulator, but also store the result in the accumulator itself. Hence the accumulator 73 operates as both a source (operand) and a destination (result) register. Additional registers such as input data latch 74, data bus buffer 75 and address latches and buffers 76 are useful in manipulation of other seismic results and data, since they eliminate the need to shuffle results back and forth between the external memory units of FIG. 5 and accumulator 73.
The internal organization of the CPU 51 can be split into two sections. In general, the instructions obtained from the program memory 54 (FIG. 5) are executed by implementing a series of data transfers in the register section indicated at 80. The control lines which actually cause the data transfers to take place are generated in the control section generally indicated at 81. Instructions enter the CPU 51 on the bus 55 (FIG. 5), are latched into instruction register 69 and are then decoded along with timing signals at instruction decoder 68 to generate the register control signals to the various other elements of the system.
The control unit 65, including timing controller 67, keeps track of the specific cycle being executed. These units are set to "T0" for each instruction fetch cycle and are advanced at the beginning of each "Phase One" clock pulse. Each instruction starts in T0 and goes to T1, T2, T3, etc., for as many cycles as are required to complete execution of the instruction. Each data transfer, etc., which takes place in the register section 80 is caused by decoding the contents of both the instruction register 69 and the timing control unit 65.
Additional control lines which affect the execution of the instructions are derived from the interrupt logic unit 71 and from the process status register 70. The interrupt logic unit 71 controls the CPU interface to the interrupt inputs to assure proper timing, enabling, sequencing, etc., which the CPU 51 recognizes and services.
The processor status register 70 also can contain a set of latches which serve to control certain aspects of the CPU operation, say to indicate the results of CPU arithmetic and logic operations, and to indicate the status of data either generated by the CPU or transferred into the CPU from outside.
Since the real work of the CPU is carried on in the register section 80 of the CPU, perhaps a detailed study is in order. The components of section 80 include:
Data Bus Buffers 75
Input Data Latch (DL) 74
Program Counter (PCL, PCH) 72
Accumulator (A) 73
Arithmetic Logic Unit (ALU) 63
Stack Pointer (S) 82
Index Registers (X, Y) 77
Address Bus Latches/Buffers 76
Process Status Register (P) 70
Again referring to FIG. 5, at 1 mHz the data which come into the CPU from the program memory 54, the data memory 53 or from peripheral devices appears on the bus 55 during the last 100 nanoseconds of Phase Two. No attempt is made to actually operate on the data during this short period. Instead, it is simply transferred into the input data latch 74 (FIG. 6) for use during the next cycle. The data latch 74 serves to trap the data on the bus during each Phase Two pulse. It can then be transferred onto one of the internal buses, say bus 64, and from there into one of the internal registers. For example, data being transferred from memory into the accumulator 73 will be placed on the internal data bus 64 and will then be transferred from the internal data bus 64 into the accumulator 73. If an arithmetic or logic operation is to be performed, say using the data from a memory unit and the contents of the accumulator 73, data in the input data latch 74 will be transferred onto the internal data bus 64 as before. From there, it will be transferred into the ALU 63. At the same time, the contents of the accumulator 73 will be transferred onto a bus in the register section 80 and from there into the second input to the ALU 63. The results of the arithmetic or logic operation will be transferred back to the accumulator 73 on the next cycle by transferring first onto the bus and then into the accumulator 73. All of them data transfers take place during the Phase One clock pulse.
The program counter (PCL, PCH) 72 provides the addresses which step the CPU through sequential instructions to the program. Each time the CPU fetches an instruction from program memory, the contents of PCL 72 are placed on the low-order 8 bits of the address bus and the contents of PCH 72 are placed on the high-order 8 bits. The counter 72 is incremented each time an instruction or data is fetched from program memory.
Accumulator 73 is a general-purpose 8-bit register which stores the results of most arithmetic and logic operations. In addition, the accumulator 73 usually contains one of the two data words used in these operations.
All logic and arithmetic operations take place in the ALU 63. This includes incrementing and decrementing of internal registers (except PCL and PCH 72). However, the ALU 63 cannot store data for more than one cycle. If data are placed on the inputs to the ALU 63 at the beginning of one cycle, the result is always gated into one of the storage registers or to external memory during the next cycle. Each bit of the ALU 63 has two inputs. These inputs can be tied to various internal buses or to a logic zero; the ALU 63 then generates the SUM, AND, OR, etc., function using the data on the two inputs.
Stack pointer (S) 82 and the two index registers (X and Y) 77 each consist of 8 simple latches. These registers store data which are to be used in calculating addresses in data memory.
Address bus buffers 76 consist, in part, of a set of latches and TTL compatible drivers. These latches store the addresses which are used in accessing the peripheral devices (ROM, RAM, and I/O's of FIG. 5.
Table VI illustrates the total CPU instructions which must be used to carry out the method of the present invention.
In brief, the CPU instruction set is divided into three basic groups. The first group has the greatest addressing flexibility and consists of the most general-purpose instructions such as Load, Add, Store, etc. The second group includes the Read, Modify, Write instructions such as Shift, Increment, Decrement and the Register X movement instructions. The third group contains all the remaining instructions, including all stack operations, the register Y, compares for X and Y and instructions which do not fit naturally into Group One or Group Two.
There are eight Group One instructions, eight Group Two instructions, and all of the 39 remaining instructions are Group Three instructions.
The three groups are obtained by organizing the OP CODE pattern to give maximum addressing flexibility (16 addressing combinations) to Group One, to give 8 combinations to Group Two instructions and the Group Three instructions are basically individually decoded.
Group One Instructions
These instructions are: Add With Carry (ADC), (AND), Compare (CMP), Exclusive Or (EOR), Load A (LDA), or (ORA), Subtract With Carry (SBC), and Store A (STA). Each of these instructions has a potential for 16 addressing modes, although usually in practice only 8 of the available modes are used.
Addressing modes for Group One are: Immediate, Zero Page, Zero Page Indexed by X, Absolute, Absolute Indexed by X, Absolute Indexed by Y, Indexed Indirect, Indirect Indexed. The unused 8 addressing modes are used to allow, possibly, addressing of additional on-chip registers, of on-chip I/O ports, and to allow two-byte word processing.
Group Two Instructions
Group Two instructions are primarily Read, Modify, Write instructions. There are really two subcategories within the Group Two instructions. The components of the first group are shift and rotate instructions and are: Shift Right (LSR), Shift Left (ASL), Rotate Left (ROL), and Rotate Right (ROR).
The second subgroup includes the Increment (INC) and Decrement (DEC) instructions and the two index register X instructions, Load X (LDX) and Store X (STX). These instructions would normally have 8 addressing modes available to them because of the bit pattern. However, to allow for upward expansion, only the following addressing modes have been defined: Zero Page, Zero Page Indexed by X, Absolute, Absolute Indexed by X, and a special Accumulator (or Register) mode. The four shift instructions all have register A operations; the incremented or decremented Load X and Store X instructions also have accumulator modes, although the Increment and Decrement Accumulator has been reserved for other purposes. Load X from A has been assigned its own mnemonic, TAX. Also included in this group are the special functions of Decrement X which is one of the special cases of Store X. Included also in this group in the X decodes are the TXS and TSX instructions.
All group One instructions have all addressing modes available to each instruction. In the case of Group Two instructions, another addressing mode has been added; that of the accumulator and the other special decodes have also been implemented in this basic group. However, the primary function of Group Two instructions is to perform some memory operation using the appropriate index.
It should be noted for documentation purposes that the X instructions have a special mode of addressing in which register Y is used for all indexing operations; thus, instead of Zero Page Indexed by X, Y instructions have Zero Page Indexed by Y, and instead of having Absolute Indexed by X, X instructions have Absolute Indexed by Y.
Group Three Instructions
There are really two major classifications of Group Three instructions; the modify Y registers instructions, Load Y (LDY), Store Y (STY), Compare Y (CPY), and Compare X (CPX), instructions actually occupy about half of the OP CODE space for the Group Three instructions. Increment X (INX) and Increment Y (INY) are special subsets of the Compare X and Compare Y instructions and all of the branch instructions are in the Group Three instructions.
Instructions in this group consist of all of the branches: BCC, BCS, BEQ, BMI, BNE, BPL, BPC and BPS. All of the flag operations are also devoted to one addressing mode--they are: CLC, SEC, CLD, SED, CLI, SEI and CLV. All of the push-and-pull instructions and stack operation instructions are Group Three instructions. These include: BRK, JSR, PHA, PHP, PLA and PLP. The JMP and BIT instructions are also included in this group. There is no common addressing mode available to members of this group. Load Y, Store Y, BIT, Compare X and Compare Y have Zero Page and Absolute, and all of the Y and X instructions allow Zero Page Indexed operations and Immediate.
Returning to FIG. 5, the total CPU instruction set of Table VII which controls CPU 51 is stored in the program memory unit 53. The CPU 51 references the contents of the memory unit 53 in order to determine appropriate processing steps.
In operations, the addresses of the instruction set, in order, are numerically adjacent, the lower addresses containing the first instructions to be executed and the higher addresses containing later instructions.
The only time the program violates the sequential rule previously mentioned is when the instruction in one section of memory is a "jump" instruction to another section of memory.
A jump instruction can also contain the address of the instruction which is to follow it. The next instruction may be stored in any location in memory unit 53, FIG. 5, as long as the programmed jump specifies the correct address.
A special kind of program jump occurs when the stored program "calls" a subroutine. In this kind of jump, the CPU 51 is required to "remember" the contents of the program counter at the time that the jump occurs. This enables the CPU 51 to resume execution of the main program which it is finished with the last instruction of the subroutine.
A "subroutine" is a program within a program. Usually it is a general-purpose set of instructions that must be executed repeatedly in the course of a main program.
The CPU 51 has a special way of handling subroutines in order to insure an orderly return to the main program. When the CPU receives a call instruction, it increments the program counter and notes the counter's contents in a reserved memory area of the memory unit known as the stack area, which includes Such an instruction need specify no address. When the CPU 51 fetches a return instruction, it simply replaces the current contents of the program counter with the address on the top of the stack. This causes the CPU 51 to resume execution of the calling program at the point immediately following the original call instruction.
Subroutines are often "nested", that is, one subroutine will sometimes call a second subroutine. The second may call a third, and so on. This is perfectly acceptable, as long as CPU 51 has enough capacity to store the necessary return addresses, and the logical provision for doing so. In other words, the maximum depth of nesting is determined by the depth of the stack. If the stack has space for storing three return addresses, then three levels of subroutines may be accomplished.
CPU's have different ways of maintaining stack contents. Some have facilities for the storage of return addresses built into the CPU itself. Other CPU's use a reserved area of external memory as the stack area and simply maintain a stack "pointer" register, such as pointer register 82, FIG. 6, which contains the address of the most recent stack entry. The external stack hence allows virtually unlimited subroutine nesting. In addition, if the processor provides instructions that cause the contents of the accumulator 73 and other general-purpose registers to be "pushed" onto the stack or "pulled" off the stack via the address stored in the stack pointer 82, multilevel interrupt processing (as described later) is possible. The status of CPU 51 (i.e., the contents of all the registers) can be saved in the stack when an interrupt is accepted and then restored after the interrupt has been serviced. This ability to save the CPU's status at any given time is possible, even if an interrupt service routine itself is interrupted.
Each operation that the CPU can perform is identified by a unique byte of data known as an "operation code". An 8-bit word used as an operation code can distinguish between 256 alternative actions, more than adequate for most processors.
In general, all instructions begin with T0 and the fetch of the OP CODE and continue through the required number of cycles until the next T0 and the fetch of the next OP CODE. Basic terminology used in conjunction with the Figures and is as defined below, for ease of reference.
OP CODE--The first byte of the instruction containing the operator and mode of address.
OPERAND--The data on which the operation specified in the OP CODE is performed.
BASE ADDRESS--The address in Indexed addressing modes which specifies the location in memory to which indexing is referenced. The high order of byte of the base address (AB08 to AB15) is BAH (Base Address High) and the low order of the base address (AB00 to AB07) is BAL (Base Address Low).
EFFECTIVE ADDRESS--The destination in memory in which the data are to be found. The effective address may be loaded directly as in the case of page Zero and Absolute Addressing or may be calcined as in Indexing operations. The high order byte of the effective address (AB08 to AB15) is ADH and the low order byte of the effective address (AB00 to AB07) is ADL.
INDIRECT ADDRESS--The address found in the operand of instructions utilizing (Indirect), Y which contains the low order byte of the base address. IAH and IAL represent the high and low order bytes.
JUMP ADDRESS--The value to be loaded into Program Counter as a result of a Jump Instruction.
During operation, the CPU 51 fetches an instruction in two distinct operations. First, CPU 51 transmits the address in its program counter 72 of FIG. 6 to the memory. Then the memory returns the addressed byte to CPU 51, viz. to the instruction decoder 68, and uses it to direct activities during the remainder of the execution.
The mechanism by which the CPU translates an operation code into specific processing actions requires some elaboration. The 8 bits stored in the instruction decoder 68 can be decoded and used to selectively activate one of a number of output lines, in this case up to 256 lines. Each line represents a set of activities associated with execution of a particular operation code. The enabled line can be combined with selected timing pulses, to develop electrical signals that can then be used to initiate specific actions. This translation of code into action is also aided by the control circuitry associated therewith, as previously explained.
An 8-bit operation code is often sufficient to specify a particular processing action. There are times, however, when execution of the instruction requires more information than 8 bits can convey, e.g., say when the instruction references a memory location. The basic instruction code identifies the operation to be performed, but cannot specify the object address as well. In a case such as this, a 2- or 3-byte instruction must be used. Successive instruction bytes are stored in sequentially adjacent memory locations, and the CPU performs two or three fetches in succession to obtain the full instruction. The first byte retrieved from memory is placed in the CPU's instruction 66 and subsequent bytes are placed in temporary storage; the CPU then proceeds with the execution phase.
CPU 51 usually uses register-pairs to hold the address of a memory location that is to be accessed for data. If the address register is programmable (i.e., if there are instructions that allow the programmer to alter the contents of the register), the program can "build" an address in the address register prior to executing a memory reference instruction (i.e., and instruction that reads data from memory, writes data to memory or operates on data stored in memory).
ALU 63, as previously mentioned, is that portion of the CPU performing the arithmetical and logical operations on the binary data.
The ALU 63 must contain an adder which is capable of combining the contents of registers of array, the register of the CPU and accumulator 73 in accordance with the logic of binary arithmetic. This provision permits the CPU to perform arithmetical manipulations on the data it obtains from memory and from its other inputs.
Using only the basic adder, routines can subtract, multiply and divide, giving the machine complete arithmetical capabilities. In practice, however, most ALU's provide other built-in functions, including hardware subtraction, boolean logic operations and shift capabilities in conjunction with index registers 77.
As previously indicated, interrupt provisions are included on many central processors as a means of improving the CPU's efficiency. Consider the case of a computer that is processing a large volume of data, portions of which are to be output to a printer. The CPU 52 can output a byte of data within a single machine cycle, but it may take the printer the equivalent of many machine cycles to actually print the character specified by the data byte. The CPU 51 could then remain idle, waiting until the printer can accept the next data byte. If an interrupt capability is implemented to the computer, the CPU can output a data byte, then return to data processing. When the printer is ready to accept the next data byte, it can request an interrupt. When the CPU 51 acknowledges the interrupt, it suspends main program execution and automatically branches to a routine that will output the next data byte. After the byte is output, the CPU 51 continues with main program execution. Note that this is, in principle, quite similar to a subroutine call, except that the jump is initiated externally rather than by the program.
More complex interrupt structures are possible, in which several interrupting devices share the same processor but have different priority levels. Interruptive processing is an important feature that enables maximum utilization of a processor's capacity for high system throughput.
Another important feature than improves the throughput of a CPU is the hold. The hold provision enables direct memory access (DMA) operations.
In ordinary input and output operations, the CPU 51 supervises the entire data transfer. Information to be placed in memory is transferred from the output device to the CPU, and then from the CPU to the designated memory location. In similar fashion, information that goes from memory to output devices goes by way of the CPU.
Some peripheral devices, however, are capable of transferring information to and from memory much faster than the CPU itself can accomplish the transfer. If any appreciable quantity of data must be transferred to or from such a device, then system throughput will be increased by having the device accomplish the transfer directly. The CPU 51 must temporarily suspend its operation during such a transfer, to prevent conflicts that would arise if the CPU and peripheral device attempted to access memory simultaneously. It is for this reason that a hold provision is provided.
MODIFICATION
To repeat, in estimating gas-liquid saturation effects in sands having any degree of consolidation, inputting only compressional-wave velocity, density and porosity is all that is required in order to provide estimates of the compression (or bulk) moduli, K*, of fluid-saturated sands; the compression moduli, Ka, of sands with empty pores (or with pore-fluid pressures held constant Kfc=K*-Ka; and the shear moduli, G* of fluid-saturated sands (or Ga, of empty sands).
And for an isotropic and homogeneous medium, the velocity, density and moduli are as related above:
where .rho..sub.b is the bulk density of the formation.
If both compressional- and shear-wave velocities are known, of course, Equation 18 gives K* directly. But no shearwave logs exist for most basins under study; furthermore, if they exist, the results are unreliable in loose (unconsolidated) sands.
For convenience and brevity, a program listing has been formulated to ease the solution of the above equations and is set forth below in the Appendix. In such listing, several lower-case symbols will represents dimensionless ratios of moduli where a=Ka/Km, b=K*/Km, g=G*/Gm=Ga/Gm, and F=(Km/Kf)-1. Gassmann's Equation 20d results: ##EQU12##
If G*, P*, Km, Kf and .phi. are known, calculation of Ka through 18 and 20d is of course possible. To calculate P* for a pore fluid other than brine only, i.e., for gas (or oil) and brine, a value of K with the new value of Kf in Equations 18 and 20d is used. But since G* and Km are not usually known, the algorithm provided in the listing uses density, porosity, and compressional-wave velocity data from well logs in the manner previously discussed. It has been proven to be surprisingly accurate in computing the change in compressional-wave velocity due to a change in pore fluids.
EXAMPLES
Let us assume P* and .phi. for a brine-saturated sand without large amounts of carbonate or high-elastic-modulus iron minerals, i.e., P* is available from compressional-wave velocity, density, and Equation 18 above. To estimate the other elastic parameters and the compressional-wave velocity, variations in brine-oil-gas saturations must be accounted for.
For the assumed conditions, the modulus Kf=25 Kbar (to be within the accuracy needed) and is pressure-independent. However, an estimate of Km, Gm, and a relationship between a and g are also needed.
To estimate Km and Gm, all information on rock composition of a given area (basin, formation, well, etc.) must be used, i.e., using samples or logs (nuclear, etc.). But mean values for Km and Gm are fairly accurate for many sands, consolidated or not, viz:
The standard deviations of the above quantities for several samples measured were 7% and 12% respectively. Extensive calculations with and without physical models, with sandstone data, produced an acceptable empirical relationship between a and g:
where c is an adjustable empirical parameter, a and g are as described above, and f.sub.3 is a functional operator which relates to the variables of interest in an empirical fashion. It was surprising found by minimizing the rms error between predicted and measured moduli values and is pressure-independent in the above approximation. The correlation coefficients for estimating K* and G* using Equation 34 are both over 0.95.
Final data base includes suites of "hard" sandstones and "soft" sands, with densities and porosities measured and with all compressional-wave and shear-wave velocities measured as functions of pressure. In addition, compositions were determined by an overlapping suite of measurements to give estimates of Km and Gm. Equation 18 supra gives G*, P* and K*; with values of Km and Gm, a and g can then be provided. Equation 34 then gives a c-value for each sample. Such c-values were found not to be significantly correlated with either G, .phi., Km/Gm, Gm, or some other pertinent parameters.
Samples containing substantial amounts of carbonates or high-elastic-modulus iron minerals, e.g., had somewhat negative correlations of c with Km. Deleting such samples from the data base reduced the suite of hard sandstones slightly. A slight negative correlation also existed between c and pressure.
By adopting Equations 33 and 34 and the value Kf=25 Kbar, all the information needed to use Equations 18 and 20d become available. The computation, however, is easy to do only "backwards"; that is, by starting with some value of (a+g) and computing (a-g), a and g. Then Equation 20d gives b, and finally, Equation 18 gives P*. Of course, use of the definitions of the various quantities as indicated above is mandatory.
To use P* and .phi. to compute K*, Ka, G*, etc., one approach is to program an interative solution as set forth below. But also, another approach would be to start with various values of (a+g) and construct plots or of the desired quantities as functions of P* for various values of .phi., as also indicated below.
The complete algorithm of the program listing set forth below was used to compute Ka and G*, using P* brine values only.
The observed or calculated quantities in Table X can then be used to compute seismic velocities in gas and water sands with good accuracy, using steps and equations previously mentioned. Such operations indicate that Equations 18, 20d, 33 and 34 are valid to estimate gas-sand velocities, if Vp, .phi. and .rho.b for a corresponding 100% brine-saturated section are known, in accordance with above assumptions.
Data banks for the process of the invention, of course, must come from a 100%-brine-saturated section of a well. If Km and Gm are known from mineralogical analysis, such values are used. If not, the assumed values of 393 and 361 Kbar, respectively, are available. The solution to Gassmann's equation is very sensitive to errors in Km only when the porosity is low. As the major application of these equations is intended to be in high-porosity "soft" rocks, errors in Km are not apt to be a major drawback.
The process by which the above calculation of the aspects of the present invention were used is as set forth below in steps I-III. Thereafter, velocities, amplitudes, reflectivities, etc. of the strata of interest were calculated in the manner previously mentioned.
Step I
Calculated the observed P-wave modulus from the well, Pow*, in kilobars. 1 Kbar=10.degree. dynes/cm.sup.2, or 1.45.times.10.sup.4 psi
STEP II
Compute Ka and G*, using the above algorithm. (The caret sign refers to a predicted quantity, as previously indicated).
STEP III
Use Equations 18, 20d, 33 and 34 to obtain a general expression for Vp for any fluid saturant: ##EQU13## where Kfg equals the bulk modulus for a brine-gas mixture. It is given by Equation 24, supra,
where Sg is the fractional gas saturation, Kg is the bulk modulus of the gas, and Kbr is the bulk modulus of the brine. In this approximation, Kbr=25 Kbar and may be assumed pressure- and temperature-independent; Kg is both temperature- and pressure-dependent. Estimate for both Kg and .rho.g for methane, as set forth above.
Correction must be made to .rho.b for gas saturation and density:
where
Use
for most cases. Incorrect estimation of .rho.g and Kg may cause serious errors.
There is some evidence of changes in seismic response for some "soft" rocks when the fluid changes from brine to highly compressible liquid hydrocarbons. The range of Kf is from about 25 Kbar for brines to about 15 for most crude oils to about 5 or 3 for liquid hydrocarbons containing much dissolved gas (but no free gas); and from about 2 to 0.001 Kbar for gases, depending on pressure and temperature. Thereafter, synthetic seismic values can be established and compared with field records in the manner previously described.
Iterative Method in Detail
The objective is to start with assumed values of Km, Gm, and Kf, and the measured values of P* and .phi.; and then to use the empirical relationship in Equation 34 and Gassmann's Equation 20d supra to estimate Ka, K*, Kfc and G*. A simple iterative calculation with interval-halving to update the variable (a+g) of Equation 34 is available, and such sequence is as follows:
(1) From the current value of (a+g), calculate (a-g); then solve for a and g, and use Gassmann's equation 20d to get b. Use Equations 33 and 18 and the definitions of the quantities to get trial values of K*, G* and finally .rho.*. Let the trial value of P* be P'; the error in P* is then P'-P*.
(2) If the sign of the error is opposite to that of the preceding error, divide the increment for updating (a+g) by -2. That is, if (P'-P*) has crossed zero, reverse directions and use steps only half as large.
(3) Update (a+g) by adding the current size and sign of the increment, as mentioned above.
(4) If the increment [or (P'-P*)] is satisfactorily small, terminate the iteration. Otherwise, repeat from Step (1).
Suitable starting values are (a+g)=0, increment=+0.4, and negative (P'-P*). The results for (a+g), (a-g) and K* are easily used to calculate the various parameters.
APPENDIX
Having now briefly described the operation of the CPU 51 and of an instruction set for operation of a particular type microprocessor system, e.g. MGS 6500, the program listing for carrying out, in the main, the method of the present invention is of consequence and is set forth below in detail.
The listing is divided into two sections. Section A contains comments as to the purpose and scope of the method. E.g., Steps 00015-00020 detail the purpose and limitations of the method; Steps 00022-00083 indicate the form and nature of the input data; Steps 00085-00098 state the form the nature of the output information; Steps 00100-00105 indicate what error messages are to be expected; Step 00107 indicates the nature and scope of subroutines, etc.
Section B describes the method in detail for calculating various petrophysical parameters of interest. Note that each instruction number that is followed by a ";" describes the particular operation that is carried by the subsequent instructions. E.g., the Step 00116 is "SG(1)=FLOAT(0)" and the Steps 00117-00130 show the detailed steps for carrying out that command; similar circumstances occur at Steps 00131; 00146; 00161, etc. These commands provide in sequence form the operations previously described with reference to FIGS. 5 and 6. Comparison of actual versus synthetic values then occurs on a separate internal (or external) basis as previously described. ##SPC1## ##SPC2## ##SPC3## ##SPC4## ##SPC5## ##SPC6## ##SPC7## ##SPC8##
CONCLUSION
Diagnostic capability of the method of the present invention can be perhaps better illustrated by an example of actual use.
Seismic data were obtained in a field environment and are depicted in FIG. 7. The geologic section was limestone embedded in shale. The gas-related anomaly is located over the indicated area.
From well logging and geologic data, the method of the present invention was used to provide the synthetic depthogram of FIG. 8. Note the close similarity of relative amplitudes for the events of interest. From FIG. 8, it was determined that the area of interest was 70% gas-saturated. Field data affirmed the estimate.
From the above, it is apparent that the method of the present invention as described hereinbefore provides a geophysicist with a strong tool for determining amounts of gas/oil-filled reservoirs in a variety of structural combinations, e.g., in sand, limestone, etc. However, the invention is not limited to the above structural combinations but is applicable to other anomalous circumstances as known to those skilled in the art. E.g., while the above presentations were directed in the most part to creating synthetic gas/oil-indicating characteristics for sandstones at different saturations, other strata can be used, as previously indicated, e.g., limestone. Other additional strata could be mapped in similar fashion, viz. dolostone, without departing from the intended scope of the claims appended hereto. Additionally, while gas-brine characteristics were mentioned in some detail in the specification, multi-phase systems would also be implemented, i.e., a mixture of brine, oil and gas. In such case, the factor Kf can be adjusted to take into account such combinations by noting that compressibility of any fluid mixture is the volume-weight average of the compressibilities of the components. Similarly, such characteristics could relate to estimations of oil-brine mixtures or to oil alone without departing from the intended scope of the present invention. Still further, actual data from logged wells (as provided, e.g., by acoustic well-logging tools) may serve as the actual data base to which the synthetic information as set forth herein can be compared. It should thus be understood that the invention is not limited to any specific embodiments set forth herein, as variations are readily apparent, and thus the invention is to be given the broadest possible interpretation within the terms of the following claims.
Claims
We claim:
1. System for resolving event characteristics of an acoustic log of a first well so as to determine gas/oil fractions in a selected zone of an earth formation by creating synthetic gas-indicating characteristics from petrophysical parameters corresponding to each selected zone of interest, comprising: (A) a digital programmable computer for: (i) predicting brine-saturated bulk and shear moduli (Kw* and Gw*) using empirical equations containing coefficients and a P-wave modulus value (Pw*), said coefficients varying in empirical fashion as a function of at least porosity, pressure and temperature, said P-wave modulus value (Pw*) being a P-wave modulus value which corresponds to said zone of interest based on well logging geologic data; (ii) estimating a frame modulus value (Ka) of said zone of interest using a relationship involving Kw*, Kb, Km and .phi., where Kw* is brine-saturated bulk density of the aggregate, Kb is the bulk modulus of brine, and Km is the bulk modulus of the solid material making up the zone of interest and .phi. is porosity; (iii) predicting bulk modulus values Kg* with the gas/oil present in differing amounts in said zone of interest using a relationship involving Ka of step (ii), Km, .phi., Kfg, where Kfg is the bulk modulus of a gas/oil mixture; (iv) assuming Gw* is equal to Ggas/oil*, predicting P-wave modulus of the zone with gas/oil present (Pgas/oil*) in differing amounts based on the relationship (v) predicting at least acoustic velocities within said zone from said P-wave modulus values (Pgas/oil*) of step (iv) as a function of gas/oil saturation; and (B) comparator means comparing said predicted velocities or derivations thereof with actual velocities or derivations thereof, from said acoustic record to predict presence and amounts of gas/oil fractions in said zone of interest.
2. System of claim 1 in which (B) repeats in sequence velocity comparisons for different zones of said earth formation and at different depths within said formation so as to provide a multiplicity of compared and predicted gas/oil-saturated zones of said formation as a function of depth.
3. System of claim 1 in which said empirical equations of (i) are o the form where Kw* is the brine-saturated bulk modulus of said zone of interest; Gw* is the brine-saturated shear modulus of said zone of interest; Pol is P-wave modulus observed in a laboratory environment; and f.sub.1 and f.sub.2 describe operators which relate the variables of interest.
4. System of claim 1 in which said frame modulus Ka of (i) is estimated based on the relationship ##EQU14## where Ka is the frame modulus; Km is the bulk modulus of solid material; Kf is the bulk modulus of the fluid and .phi. is the porosity.
5. System for resolving event characteristics of a seismic record so as to determine gas/oil fractions in a selected zone of an earth formation by creating synthetic gas-indicating characteristics from brine-saturated petrophysical parameters corresponding to each selected zone of interest, comprising a programmable digital computer including a CPU-control means, program and data memories connected to said CPU-control means and display-storage means also connected to said CPU-control means, said CPU-control means operating on petrophysical data within said data memory control of program statements initially held in said program memory but sequentially fetched to said CPU-control means so as: (i) from brine-saturated bulk and shear moduli (Kw* and Gw*), calculated by empirical equations including a P-wave modulus value (Pw*) which corresponds to said zone of interest based on well logging geologic data and estimated frame modulus values (Ka) of said zone of interest, to predict bulk modulus values Kg* with the gas/oil present in differing amounts in said zone of interest; (ii) to determine P-wave modulus values of the zone with gas/oil present (Pgas/oil*) in differing amounts based on the relationship (iii) to predict from seismic velocities within said zone (Pgas/oil*) calculated from step (B) as a function of gas/oil saturation; (iv) to compare said predicted seismic velocities or derivations thereof with actual velocities or similar derivatives thereof from said seismic record to predict amounts of gas/oil fractions in said zone of interest.
6. System of claim 5 whereby (iv) is an indexable multi-dimensional array whereby said actual velocities are sequentially matched with said predicted velocity values as a function of gas/oil saturation until a best match is obtained.
7. System of claim 5 in which said CPU-control means additionally operates to provide: amplitude reflection values calculated based on a relationship involving velocity and density of said zone of interest and incident amplitude information (Ai) and in which (iv) is an indexable multi-dimensional array whereby said actual reflection values from said seismic record are sequentially matched with said predicted reflection values as a function of gas/oil saturation until a best match is obtained.
8. System of claim 7 in which (iv) includes first means for providing a plurality of synthetic amplitude-v.-time seismic traces from said synthetic reflection data as a function of gas/oil saturation and second means operatively connected to said first means for cross-correlating said synthetic traces with actual traces of said record until a best match is obtained.
9. System of claim 7 in which said estimated amplitude reflection values of (iv) are estimated in accordance with ##EQU15## where Ar is the amplitude from the reflected signal at said zone of interest and Ai is the amplitude of the incident signal; V.sub.1 is the compressional velocity of the wave in an overlying medium 1; V.sub.2 is the compressional velocity in the zone of interest below the contact line; d.sub.1 is the density of the overlying medium 1; and d.sub.2 is the density of the underlying medium 2.
10. System for resolving event characteristics of a seismic record so as to determine gas/oil fractions in a selected zone of an earth formation by creating synthetic gas/oil-indicating characteristics from brine-saturated petrophysical parameters corresponding to each selected zone of interest, comprising: (A) a programmable digital computer for: (i) from calculated normalized brine-saturated petrophysical parameters involving a laboratory and/or well-log generated brine-saturated P-wave modulus value for a brine-saturated rock corresponding to said zone of interest, predicting compressional velocities of said zone with gas/oil in differing amounts present in said zone, in accordance with the relationship ##EQU16## where .rho.b is the bulk density, .phi. is the porosity, G* is the shear modulus at the saturations of interest, Kfg equals the bulk modulus for the liquid-gas mixture, Ka is the frame modulus, Km is the bulk modulus of a hypothetical rock, i.e., same composition but perfectly cemented with zero porosity; (ii) from seismic velocities of said zone calculated from (i), predicting synthetic seismic amplitude information with gas/oil present in said zone; and (iii) comparing said predicted amplitude information with gas/oil present with actual field-acquired amplitude data associated with said seismic record so as to determine both presence and amounts of gas/oil within said zone of interest.
11. System of claim 10 in which (i) is further characterized by: simultaneously solving through iteration the following equations: ##EQU17## to provide estimations of Ka and G* as a function of fluid saturation, where: .rho.b is the bulk density of the zone of interest, Vp is the compressional velocity, P* is the P-wave modulus at said saturations of interest, K* is the bulk modulus at said saturations of interest, G* is the shear modulus at said saturations of interest, Ka and Ga are frame moduli, Kfc is the fluid correction equal to K*-Ka, Km and Gm are the bulk and shear moduli of a hypothetical rock (i.e., same composition but perfectly cemented with zero porosity), .phi.=porosity a=Ka/Km b=K*/Km g=G*/Gm=Ga/Gm f=(Km/Kf)-1 C is an empirical parameter and f.sub.3 describes a functional operator which relates to the variables of interest of Equation (3).
12. System of claim 11 in which said simultaneous iteration involves: (i)' predicting brine-saturated bulk and shear moduli (Kw* and Gw*) using empirical equations containing coefficients and a P-wave modulus value (Pw*), said coefficients varying in empirical fashion as a function of at least porosity, pressure and temperature, said P-wave modulus value (Pw*) being a P-wave modulus value for a brine-saturated rock which corresponds to said zone of interest based on well logging geologic data, said shear modulus Gw* being equal to Ggas/oil* where Ggas/oil* is the shear modulus of said zone with gas/oil present; (ii)' estimating frame modulus values (Ka) of said zone of interest using a relationship involving Kw*, Kb, Km and .phi., where Kw* is brine-saturated bulk density of the aggregate, Kb is the bulk modulus of brine, and Km is the bulk modulus of the solid material making up the zone of interest and .phi. is porosity.
13. System of claim 10 in which (iii) is provided by a separate comparator means connected to said computer.
14. Method for resolving event characteristics of an acoustic record so as to determine gas/oil fractions in a selected zone of an earth formation by creating synthetic gas-indicating characteristics from brine-saturated petrophysical parameters corresponding to each selected zone of interest, comprising: (a) predicting brine-saturated bulk and shear moduli (Kw* and Gw*) using empirical equations containing coefficients and a P-wave modulus value (Pw*), said coefficients varying in empirical fashion as a function of at least porosity, pressure and temperature, said P-wave modulus value (Pw*) being a P-wave modulus value which corresponds to said zone of interest based on well logging geologic data; (b) estimating a frame modulus value (Ka) of said zone of interest using a relationship involving Kw*, Kb, Km and .phi., where Kw* is brine-saturated bulk density of the aggregate, Kb is the bulk modulus of brine, and Km is the bulk modulus of the solid material making up the zone of interest and .phi. is porosity; (c) predicting a bulk modulus value Kg* with the gas/oil present in said zone of interest in a predetermined amount using a relationship involving Ka of step (b), km, .phi., Kfg, where Kfg is the bulk modulus of a gas/oil mixture; (d) assuming Gw* is equal to Ggas/oil*, predicting P-wave modulus value of the zone with gas/oil present (Pgas/oil*) in said predetermined amount based on the relationship (e) predicting at least an acoustic velocity within said zone from said P-wave modulus value (Pgas/oil*) of step (d); (f) repeating steps (c)-(e) for gas/oil-brine mixtures of different gas/oil saturations to provide a suite of synthetic acoustic velocities for said zone of interest as the function of a plurality of gas/oil saturations; and (g) comparing said suite of predicted velocities with actual velocity information produced from said actual acoustic record to determine presence and amounts of gas/oil fractions in said zone of interest.
15. Method of claim 14 in which steps (a)-(g) are repeated in sequence for different zones of said earth formation and at different depths within said formation so as to provide a multiplicity of compared and predicted gas/oil-saturated zones of said formation as a function of depth.
16. Method of claim 14 in which step (g) is further characterized by (i) indexing predicted and actual velocity values as a function of gas/oil saturations and (ii) determining best match of all indexed values so as to accurately determine amounts of gas/oil within said zone of interest.
17. Method for resolving event characteristics of an acoustic record so as to determine gas/oil fractions in a selected zone of an earth formation by creating synthetic gas/oil-indicating characteristics from brine-saturated petrophysical parameters corresponding to each selected zone of interest, comprising the steps of: (a) from calculated normalized brine-saturated petrophysical parameters involving a laboratory or otherwise annotated brine-saturated P-wave modulus value corresponding to said zone of interest, predicting P-wave modulus values with gas/oil in differing amounts present in said zone, in accordance with the relationship where Pgas/oil* is the predicted P-wave modulus with gas/oil present in predetermined amounts; Kgas/oil* is the predicted bulk modulus of said zone with gas/oil present in said predetermined amounts; Gw* is the brine-saturated shear modulus of said zone; (b) from acoustic velocities of said zone calculated from (a), predicting synthetic acoustic values associated with different gas/oil fractions present in said zone; and (c) comparing said synthetic acoustic values with actual field-acquired values so as to determine both presence and amounts of gas/oil within said zone of interest.
18. Method of claim 17 in which said resolution of said acoustic record is limited to estimations of gas fractions only.
19. Method of claim 17 in which said resolution of said acoustic record is limited to estimations of oil fractions only.
20. Method for resolving event characteristics of an acoustic record so as to determine gas/oil fractions in a selected zone of an earth formation by creating synthetic gas/oil indicating characteristics from brine-saturated petrophysical parameters corresponding to each selected zone of interest, comprising the steps of: (a) from calculated normalized brine-saturated petrophysical parameters involving laboratory and/or well-log generated brine-saturated P-wave modulus values corresponding to said zone of interest, predicting velocities of said zone with gas/oil in differing amounts present in said zone, in accordance with the relationship ##EQU18## where .rho.b is the bulk density, .phi. is the porosity, G* is the shear modulus at the saturations of interest, Ka is the frame modulus, Kfg equals the bulk modulus for the liquid-gas mixture, Km is the bulk modulus of a hypothetical rock, i.e., same composition but perfectly cemented with zero porosity; (b) from velocities calculated from (a), predicting synthetic acoustic information as a function of gas/oil present in said zone; and (c) comparing said predicted acoustic information with actual field-acquired data of said acoustic record so as to determine both presence and amounts of gas/oil within said zone of interest.
21. Method of claim 20 in which said frame modulus Ka of step (a) is estimated based on the relationship ##EQU19## where Ka is the frame modulus; Km is the bulk modulus of solid material; Kf is the bulk modulus of the fluid and .phi. is the porosity.
22. Method of claim 21 in which bulk modulus Kfc is equal to the second term of said relationship whereby said relationship becomes
23. Method of claim 20 in which step (a) is further characterized by: simultaneously solving through iteration the following equations: ##EQU20## to provide estimations of Ka and G* as a function of fluid saturation, where: .rho.b is the bulk density of the zone of interest, Vp is the compressional velocity, P* is the P-wave modulus at said saturations of interest, K* is the bulk modulus at said saturations of interest, G* is the shear modulus at said saturations of interest, Ka and Ga are frame moduli, Kfc is the fluid correction equal to K*-Ka, Km and Gm are the bulk and shear moduli of a hypothetical rock (i.e., same composition but perfectly cemented with zero porosity), .phi.=porosity a=Ka/Km b=K*/Km g=G*/Gm=Ga/Gm f=(Km/Kf)-1 C is an empirical parameter and f.sub.3 describes a functional operator which relates to the variables of interest of Equation (3).
24. Method for resolving events of a seismic record associated with a selected zone of an earth formation so as to determine gas/oil fractions therein by creating synthetic gas/oil indicating characteristics from brine-saturated petrophysical parameters corresponding to said selected zone of interest, comprising the steps of: (a) from calculated normalized brine-saturated petrophysical parameters involving a laboratory or otherwise annotated brine-saturated P-wave modulus value corresponding to said zone of interest, predicting P-wave modulus values of said zone with gas/oil present in differing amounts in said zone in accordance with the relationship where Pgas/oil* is a predicted P-wave modulus with gas/oil present in said predetermined amounts; Kgas/oil* is a predicted bulk modulus of said zone with gas/oil present in said predetermined amounts; Gw* is the brine-saturated shear modulus of said zone; (b) from seismic velocities calculated from (a), predicting synthetic seismic amplitude values with gas/oil present in said zone in said predetermined amounts, and (c) comparing predicted amplitude values with actual values of said seismic record whereby amounts of gas/oil within said zone of interest can be accurately estimated.
25. Method of claim 24 in which said predicted amplitude values of step (b) are determined in accordance with ##EQU21## where Ar is the amplitude from the reflected signal at said zone of interest and Ai is the amplitude of the incident signal; V.sub.1 is the compressional velocity of the wave in an overlying medium 1; V.sub.2 is the compressional velocity in the zone of interest below the contact line; d.sub.1 is the density of the overlying medium 1; and d.sub.2 is the density of the underlying medium 2.
Patent Citations (2)
| Patent | Date | Inventor | Cited By |
|---|---|---|---|
| US3241102 | 1966-03-01 | Peterson | |
| US4293933 | 1981-10-01 | Park et al. |
Non-Patent Literature (1)
- Levitsova et al., "Razvedochnaya Geofizika", 1966, pp. 3-17, #15, G-222.