Method for determination of earth stratum elastic parameters using seismic energy
This is a continuation of copending application Ser. No. 707,838, filed on Mar. 4, 1985, and now abandoned. 
    
    
     BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The invention relates generally to computer-aided processing of seismic energy and, more particularly, but not by way of limitation, it relates to a specific method of processing wherein the elastic parameters of an earth substratum, such as pressure modulus, density and shear modulus, are determinable in relative value using conventionally acquired seismic survey data. 
     2. Description of the Prior Art 
     Applicants are not aware of any prior art that relates to the procedures of the present invention. Prior inquiry indicates that U.S. Pat. Nos. 4,316,267 and 4,316,268 are of interest due only to the fact that they relate generally to identification of gas-bearing hydrocarbon zones from seismic data; however, these teachings are directed to a different form of compressional wave evaluation wherein a progressive comparison of Common Depth Point event indications is examined for high intensity amplitude anomalies in adherence to selected criteria which will tend to verify presence of gas in the formation. The primary teaching of the patent is of little interest relative to the present invention as it is directed to an entirely different form of evaluation procedure. 
     SUMMARY OF THE INVENTION 
     The present invention relates to a method for processing common depth point seismic data to derive relative changes in the elastic parameters of substratum. Resolution of pressure modulus, density and shear modulus relative change is effected by inversion of acoustic reflection data for a selected terrain. The method processes plural groups of common depth point data using the Born approximation for the obliquity of downward propagating energy to determine the relative changes in the elastic parameters at downwardly displaced formations. The CDP seismic data, a function of surface offset and energy travel time, is processed along with an estimate of selected interval velocities along the depth of survey in order to determine a function of the angle of incidence of energy at each reflector depth. A linear system relating to the obliquity factor and including relationships for pressure modulus, density and shear modulus is then solved by matrix inversion for each depth of interest. 
     Therefore, it is an object of the present invention to provide a method for deriving direct gas production likelihood indication utilizing conventional common depth point seismic data. 
     It is also an object of the present invention to provide a method for deriving data indicative of gas bearing formations with greater reliability. 
     Finally, it is an object of the invention to provide a method that is applicable to conventional seismic file data to provide yet another interpretation tool relative to hydrocarbon presence. 
    
    
     Other objects and advantages of the invention will be evident from the following detailed description when read in conjunction with the accompanying drawings which illustrate the invention. 
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is an idealized drawing of an earth section showing selected strata and seismic energy propagation; 
     FIG. 2 is a flow diagram of the computational program of the present invention for deriving elastic parameter data output; 
     FIG. 3 is a photographic reproduction of a seismic survey line showing relative amplitude stack; 
     FIG. 4 is a photographic reproduction of a seismic survey line showing relative change in pressure modulus; 
     FIG. 5 is a photographic reproduction of a seismic survey line showing relative change in shear modulus; and 
     FIG. 6 is a photographic reproduction of a seismic survey line illustrating relative change in density. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The present invention accomplishes inversion of seismic reflection data, using a programmed digital computer, to determine selected elastic parameters for earth formations at various depths. The method functions with input of multiple offset seismic data of the type commonly collected in field surveys and processed as common depth point or common mid-point data. This is a well-known seismic energy evaluation technique as used extensively in the industry, such common depth point data being preprocessed in conventional manner to correct for static and dynamic signal variations. The present method treats the data still further to determine and output the relative changes of the elastic parameters pressure modulus, density and shear modulus for subsurface interfaces at selected depths (having known average seismic interval velocity). For best results, the process is limited to locally layered media, i.e. substrata that exhibit locally negligible lateral change in elastic parameters along the length of the seismic survey line. 
     In the usual and conventional form of seismic reflection data analysis, the reflectors obliquity is ignored. That is, the process is not concerned with the dependence of the seismic reflection strength to the angle of incidence of the seismic wave front to a respective reflector. The result of such conventional analysis is to define a reflection strength at each depth or interface that is averaged over the different angles of incidence that correspond to the different offsets. Thus, composited reflection data for that stratum tends to produce an average event for each two-way vertical travel time, (i.e. depth). 
     The present method uses the Born approximation of the obliquity factor in order to determine the relative changes in the elastic parameters at each depth. FIG. 1 illustrates a typical common depth point &#34;CDP&#34; under test, or at least a portion of one, wherein a seismic source 10, such as a vibrator or other energy source, is engaged at earth surface 12 to impart energy to the subsurface. To illustrate by means of a single raypath 14 that behaves in accordance with Snell&#39;s law which states that ##EQU1## FIG. 1 illustrates the change in the angle of incidence θ as raypath 14 proceeds deeper through the substratum Z 1  -Z n  relative to common depth point 15. The interfaces 16, 18, 20 and 22 each reflects seismic energy in accordance with the angle of incidence at that depth. 
     At stratum 16, raypath energy 24 (dash-line) would reflect at angle θ&#39; 1  to point 26 on surface 12 and, similarly, reflected energy from each of stratum 18-22 would reflect at the respective angle θ for that particular stratum. Energy passing through the successive stratum are refracted successively at different angles θ 2  -θ n . 
     If D is designated as the reflection data, i.e. the seismic impulse response minus the direct wave, then a result of a far field WKBJ approximation of the acoustic wave equation is that D is a product of the impulse response for reflectors having no obliquity A and an obliquity factor 0. Thus, 
     
         D=A.O                                                      (2) 
    
     Divergence A depends only upon the depth z and the shotreceiver offset x. The obliquity factor 0 depends upon the angle of incidence of the raypath 14 at z (θ) and relative changes in the elastic parameters at the depth z. These relative changes are: ##EQU2## In the Born approximation, the obliquity factor 0 is ##EQU3## where the velocities shear and compressional, V s  and V p , may be taken as the locally averaged shear and compressional velocities, respectively. Thus, it is only necessary to obtain a gross average of the velocities and these may be determined by standard velocity estimation methods using reflection travel times from data over a thousand feet or more. 
     FIG. 2 illustrates a computer processing flow diagram for carrying out the present process. Seismic data is input at flow stage 30 in pre-processed form and flow stage 32 orders the data into common depth point sort groups, i.e. source-receiver data having the same mid-point or common depth point. Flow stage 34 selects a group D, that is a CDP group that is a function of selected source-receiver offset x and the source-reflector-receiver travel time T. Flow stage 36 then estimates. the interval velocities as a function of depth z. This estimate can be input by the operator in accordance with previous knowledge of the field, or if prior data is not available, then specific velocity processing of the data can be performed to obtain the necessary input data. 
     Proceeding then to flow stage 38, the program functions to determine each of the offset x, divergence A, and travel time T as a function of angle of incidence θ and depth z for the depth point of interest. This involves solution of the far field WKBJ relationship for each of offset x, divergence A and two-way travel time T. For a layered, piecewise constant velocity media, the solutions are as follows: ##EQU4## given by Snell&#39;s law, Equation (1). 
     Having determined x as a monotone function of θ and z, flow stage 40 functions to express θ as a function of x and z and, subsequently, using θ as a function of x and z to express A and T as functions of x and z. Then, in flow stage 42 with T as a function of x and z, the data for each group D can be expressed as a function of x and z. With all of T, A, 0 and D expressed as functions of x and z, flow proceeds to stage 44 wherein for each z, the linear system found by combining equations (1) and (2) is solved. Thus, ##EQU5## Equation (10) is solved by matrix inversion of a linear system, a well-known computational routine, to determine the relative moduli and density as set forth in Equations (3), (4) and (5). This solution is effected by the method of least squares, and since seismic data is often contaminated, a two-pass procedure is used. If, after the first pass, some of the data is badly at variance with the first least square&#39;s estimate of the changes in the elastic parameters, the bad data is rejected and a new estimate is made. Suitable comparison criteria may be input by the operator. 
     After solution for elastic parameters at a selected z or depth on the seismic section in flow stage 44, the program recycles as indicated by flow line 46 to process data for a next selected depth z. As the elastic parameters for each depth z are determined they are available at the data output stage 48 for output in selected representation. 
     FIGS. 3-6 illustrate output data for a seismic survey line wherein surface offset distance is represented on the abscissa while depth in travel time (seconds) is illustrated on the ordinate. The output is a variable amplitude-type of presentation with white indicating negative values and black indicating positive values. The output presentations of FIGS. 3-6 clearly illustrate that a great number of intermediate depth strata have been processed in the operation; however, a stratum of particular interest lies at 1.5 seconds time (depth) where the location of existing wells No. 1 and No. 2 are indicated. Well No. 2 was a gas producing well and well No. 1 could not be produced in commercial quantity. 
     FIG. 3 is an output display of the relative amplitude stack for the seismic data of the section. A strong event response can be noted along the 1.5 second line with the white amplitude upper limits underlain by the counterpart black limits. FIG. 4 illustrates for the same section a relative change in pressure modulus display 52. Here again, one can discern the strong response at the 1.5 second time level, since in-phase with the relative amplitude stack of FIG. 3 and exhibiting nearly as strong response. 
     A section 54 of FIG. 5 for the same section illustrates the data output for relative change in shear modulus. At the 1.5 second time line, the black response indicates a positive relative change in shear modulus, i.e. the shear modulus of the sandstone gas reservoir is greater than the shear modulus of the shale overlaying the gas reservoir. The display output 56 of FIG. 6 illustrates the same section line of the relative change in density, and the response is again in-phase with the relative amplitude stack display 50 (FIG. 3), the extremely bright response variation for the well No. 2 ordinate further bears out the favorable density and capability for gas production as is the fact. 
     The foregoing discloses a novel computer implemented data handling technique wherein common depth point seismic data is effectively linearized to define an obliquity factor for selected depths that enables determination of relative indications for certain elastic parameters. The method provides output data that is more nearly determinative of hydrocarbon production capability than conventional seismic sounding. Seismic data undergoing the present method may be variously pre-processed prior to entry for inversion processing. 
     Changes may be made in combination and arrangement of steps or routines as heretofore set forth in the specification and shown in the drawings; it being understood that changes may be made without departing from the spirit and scope of the invention.