Method for Estimating Shear Wave Reflection Data From Acquired Compressional Wave Reflection Data

1. Field of the Invention

The present invention relates to a method for estimating shear wave reflection data from compressional wave reflection data and, more particularly, to a method for producing pseudo-shear wave reflection seismic sections from compressional wave reflection data.

2. Discussion of the Prior Art

The use of compressional (P) wave reflection data in geophysics analysis is notoriously well known. Typically, seismic sections are produced from compressional wave subsurface reflections which provide extremely useful information concerning formation characteristics. Recently, interest has been growing in obtaining shear (S) wave seismic sections which would provide additional useful information concerning subsurface formation characteristics which can be used in conjunction with information obtained from compressional wave seismic sections. For example, P wave seismic sections can provide useful information on the compressibility of subsurface formations, while S wave seismic sections can provide useful information on subsurface formation rigidity. However, the detection of shear wave reflections is more difficult than is the case with compressional wave reflections. The shear wave reflections are typically of lower amplitude than compressional wave reflections, making detection difficult, and the direct propogation of a shear wave into a subsurface formation to induce a shear wave reflection requires special transducers and additional steps over and above those required for obtaining compressional wave reflection data. This makes obtaining shear wave reflection data difficult, more costly and time consuming.

One object of the present invention is the provision of a method for estimating shear wave reflection data from acquired compressional wave reflection data, thus removing the requirement for a "shooting" or other special data acquisition procedures to obtain shear wave reflection data.

Another object of the invention is the provision of a method for producing pseudo shear wave seismic sections from acquired compressional wave reflection data.

The method of the invention analyzes variations in the amplitude of recorded compressional wave traces with sourcereceiver offset. Many variables affect the amplitude of the compressional wave signal with offset, including the presence of a shear wave at a reflection boundary. Each of the factors affecting amplitude is quantified and assigned a value, either estimated or measured, allowing the effects of shear wave reflection on the amplitude changes with offset to be estimated and recorded, thus producing shear wave reflection data. Seismic sections based on the estimated shear wave reflection data can then also be produced.

The above and other objects, advantages and features of the invention will be more readily understood from the following detailed description of the invention, which is provided in connection with the accompanying drawings.

FIG. 1 illustrates a marine seismic exploration system employing an array of sources and receivers;

FIG. 2 illustrates a CDP gather of compressional reflection wave traces and graphs showing trace amplitude variations with offset; and

FIG. 3 illustrates in flow chart form a computer program implementing the method of the invention.

FIG. 1 illustrates source and receiver arrays which are typically used in marine seismic surveying. It should be understood that similar arrays are used for land seismic surveying and that the present invention is applicable to data gathered in both types of surveying. When land surveying is employed, however, an additional amplitude correction is used, as discussed below, which is not needed in a marine environment.

During data gathering, the sources S.sub.1 . . . S.sub.n are activated, typically in sequence, to generate descending acoustic waves which partially reflect off subsurface boundary layers, producing ascending reflection signals (traces) which are received by receivers R.sub.1 . . . R.sub.n and recorded. During subsequent processing to produce seismic sections, the recorded signals are often gathered for common depth point (CDP) stacking. A typical common depth point gather is shown in FIG. 2 for near and far coupled pairs of sources and receivers, as well as for paired sources and receivers in between. The signals in this gather have been preprocessed by application of normal moveout (NMO) corrections in a known manner to the gathered signals. If land shooting is involved, the recorded signals in the FIG. 2 gather have also been preprocessed for surface consistent statics. Both NMO corrections and surface consistent statics corrections are well known in the art and a detailed discussion of these techniques is not repeated herein. Normally, the signal gather in FIG. 2 is then CDP stacked to enhance the signal-to-noise ratio of the reflection signals and the stacked traces may then be used to produce seismic sections of compressional reflection wave data.

The method of the invention uses variations in the amplitude of the CDP gathered compressional waves with offset, denoted herein as A.sub.p (x), at selected points in time, to estimate shear wave reflection data, which can in turn be used to produce pseudo-shear wave seismic sections corresponding to CDP stacked compressional wave seismic sections. As known in the art, offset refers to the distance between source receiver pairs used to generate a signal trace. The amplitude variation with offset A.sub.p (x) for the gathered traces at two selected points in time t.sub.1, t.sub.1 +.DELTA. at is also shown in FIG. 2.

To understand the method of the invention, a brief discussion of the factors affecting amplitude variations with offset A.sub.p (x) is necessary. These are the P wave reflection coefficient R.sub.p (the coefficient of reflection of the compressional wave at a reflection interface), travel path factors (acoustic wave spreading and transmission characteristics), geometric characteristics of the source-receiver arrays and surface boundary conditions, and random variations such as transducer coupling (mainly for geophones in a land based shooting), and shallow attenuation. Of these factors, proper choice of shooting conditions can minimize the effects of some factors (except R.sub.p) on the value A.sub.p (x). However, all of these factors, save R.sub.p, cannot be completely removed and will be compensated for in processing steps of the invention, described in greater detail below. In the invention, the various factors affecting the value A.sub.p (x) are quantified so that the effects of changes in A.sub.p (x) due to shear wave reflections can be determined to provide information on the magnitude of the shear wave reflections.

The reflection coefficient R.sub.p is one of the important factors which must be determined in the method of the invention. The well know Zoeppritz equation relates the reflection coefficient R.sub.p to the variables -.sub..rho., V.sub.p, V.sub.s, -.sub..rho., V.sub.p ', V.sub.s ', where -.rho., V.sub.p and V.sub.s are the density, compressional wave velocity and shear wave velocity in the upper medium at a formation interface and -.rho.', and V.sub.p ' and V.sub.s ' are the density, compressional wave velocity and shear wave velocity in the lower medium at the interface.

Certain assumptions can be made concerning the relationships of these factors in the Zoeppritz equation. First, the fractional velocity/density changes across a reflection boundary are small and second, that the incident angles .theta. a between source and receiver pairs are small. Typically, angles up to 22.degree. from the normal for sources and receivers can be used in the gathered traces to ensure that the second assumption is accurate. With these simplifying assumptions the Zoeppritz power series equation can be expanded, with the second order of terms dropped, yielding the following simplification:

where

.theta. is the angle from the normal of the wave path between a source-receiver pair; and,

R.sub.po is the normal incident P wave reflection coefficient; and

R.sub.so is the normal incident S wave reflection coefficient (R.sub.po and R.sub.so are the coefficients which exists for P and S waves propagated normal to a reflection boundary).

For the first order factors discussed above, e.g., travel path factors, geometric characteristics random variations, etc., which have not been removed from the gathered signals by appropriate selection of shooting conditions and preprocessing and which still affect the amplitude variation with offset A.sub.p (x), each can be rCepresented by a function

where i=1 to the number of influencing factors N.

Thus, the total of all remaining influencing effects can then be represented as:

The value A.sub.p (.theta.) can also be represented by a general equation of the form:

with

and

where ##EQU1##

Equation (10) represents the totality of amplitude effects in the gathered traces caused by the various influencing factors discussed above, e.g., travel path factors, geometric characteristics, random variations, etc., which have not been removed from the gathered traces by shooting conditions or preprocessing.

For recorded traces at known offsets, equation (7) can be rewritten as:

where x represents known offset values determined by the spatial positioning of sources and receivers in the array (FIG. 1), A.sub.p (x) is a measured amplitude value for a determined offset x, and f represents a conversion factor for offset incident angle derived from MNO corrections

where V.sub.INT is an estimated interval velocity (P wave velocity through a medium of interest) and V.sub.st is an estimated stacking velocity (effective velocity used to line up traces in a stack), and t is a two-way travel time.

Equation (11) can be solved for P and Q at each selected time point using, for example, an iterative least mean squared procedure to find the best "fit" of P and Q values. From calculated P and Q values, equation (8) can be solved to yield R.sub.po values. Moreover, the determined value of Q can be used in equation (9) to determine R.sub.so (the normal shear wave reflection coefficient) if the values of T, G and .DELTA..rho./.rho. are known. The relationship of the value T (equation 4) to Q can be illustrated as follows:

As demonstrated, if T=Vp/Vs=2.0, then the term .DELTA..rho./.rho. is eliminated as a factor in the calculation of Q. In many instances, the value of T will indeed be 2.0. However, it is also now possible to obtain direct measurements of V.sub.p /V.sub.s using compressional and shear wave borehole logs, enabling the value R.sub.so to be determined from determined values of T, Q, and R.sub.po, if .DELTA..rho./.rho. and G can also be determined.

The value .DELTA..rho./.rho. can often be predicted with reasonable accuracy from known subsurface characteristics. Typically, .DELTA..rho./.rho..about.1/n R.sub.po, where n is an observed value for a particular area and is typically 5. It is also possible to obtain .DELTA..rho./.rho. directly as a measured value from actual well logs of a surveyed area (.DELTA..rho./.rho.=.DELTA.V.sub.p /V.sub.p).

The value of G can be estimated, as it typically has the value G.about.0 when land vibrator velocity geophones are used in the source-receiver array. It also typically has the value G.about.1 for marine explosive and hydrophone sources and hydrophone receivers. However, it is also possible to more accurately estimate the value of G. A description of a technique for doing this is as follows:

Recalling equations 8 and 9,

P and Q are estimated using, for example, a least mean squared technique, for each time point in the CDP trace gather. The density term .DELTA..rho./.rho. in equation 9 can be removed by assuming .DELTA..rho./.rho.1/n(R.sub.po). Equation 9 then reduces to: ##EQU2##

In simplifying the Zeoppritz equation initially, it was assumed that T and G are slowly varying functions of time (and, of course, depth). Accordingly, a smoothing or low-pass filtered version of their constituent parts is all that is needed to estimate their values. This can be obtained from P wave and S wave log data from boreholes in the near vicinity. Thus, if P and S wave velocities (V.sub.p ',V.sub.s ') are taken from corresponding logs, the function T can be estimated as ##EQU3## where <> designates a low pass filtering operation. Next, if synthetic P and S wave reflectives R.sub.po and R.sub.so are derived from P and S wave logs, then combining equations 8 and 9 and low pass filtering provides the following estimate of the function G: ##EQU4##

Once the value of G is known, equation (9) can be solved for R.sub.so.

Thus, by first solving equation (7) using a least mean squared procedure to determine P (where P=R.sub.po) and Q and then solving equation (9) to determine R.sub.so, both the compressional wave reflection coefficient R.sub.po and the shear wave reflection coefficient R.sub.so can be determined from measured amplitude variations with offset A.sub.p (x). These values can be used to construct both compressional wave and shear wave seismic sections for interpretation and analysis.

A representative computer program implementation of the method of the invention, which can be run on commercially available general purpose digital computers, is illustrated in FIG. 3.

In the first step 101, values for V.sub.RMS, T, S and acquisition (shooting) and physical parameters are input to the computers, where V.sub.RMS represents stacking velocities (V.sub.st),

which is a filtered (smoothed) version of V.sub.s /V.sub.p data obtained from actual borehole logs, and ##EQU5## which is filtered (smoothed) statistical data also obtained from actual borehole logs. The acquisition parameters include types of sources and receivers used and array parameters. The physical parameters include an elastic attenuation parameter q. In the next step 103, NMO corrected CDP trace gathers are read into the computer (these traces are also preprocessed to correct for near surface statics, if required, e.g., for land shootings). In step 105, a time T=0 is set in a time counter. Following this, in step 107, for the time then set in the time counter, a "depth"=z is calculated as ##EQU6## as is an estimated interval velocity V.sub.INT, and a deterministic component G.sub.det of the influencing factors on amplitude versus offset behavior. For example, certain deterministic factors are governed by known equations. These include spherical spreading (SS), which is defined as

where ##EQU7## where t=a one-way travel time, source receiver directivity in arrays (SRD), which is defined as

where s represents the frequency of the source signal and 1=array length; and attenuation (AT), which is defined as

Using these equations, G.sub.det can be determined for these influencing factors. Thereafter, in step 109, a counter is set, with the near offset value x for the gathered trace having the smallest incident angle .theta. and in step 111, this incident angle .theta. is estimated using the equation ##EQU8## In subsequent step 113, the offset counter x is incremented, with the next offset value x (the offset of the next trace in the gather moving toward the trace having the farthest offset), and in step 115, a determination is made of whether all offsets (all traces) have been processed. If not, the incident angle .theta. is estimated for the newest offset x incremented at step 113. When all offsets x have processed to estimate the corresponding incident angles .theta., the processing computer proceeds to step 117, where it calculates values for P and Q using Equation (7) and a least mean squared processing technique. Following this, the computer proceeds to step 119, where an estimate is made of the component of amplitude versus offset behavior which is not due to a shear wave velocity component and is not accounted for by the deterministic component G.sub.det calculated in step 107. This estimate G is made in accordance with the equation ##EQU9## Following this, in step 121, a shear wave reflectivity component R.sub.s is calculated from the equation ##EQU10## which is a rewritten version of Equation (16). The computer then proceeds to step 123, where the TIME counter is incremented. In the next step 125, a determination is made of whether the end of the processed traces has been reached. If not, the computer proceeds back to step 107, with a new trace time value and repeats steps 107 . . . 123. If the end of the traces is reached, it is determined in step 125, following which the computer proceeds to step 127, where it determines if all NMO corrected gathers have been processed. If they have, the program ends; if not, the computer proceeds to step 103 and begins processing a new trace gather.

In step 121, the computer calculates shear wave reflectivity coefficients which can then be applied to the gathered traces to produce traces representing shear wave reflection traces. These shear wave traces can then be used in conventional manner (CDP stached) to produce pseudo-shear wave seismic sections.

While a preferred embodiment of the method of the invention has been described and illustrated, it should be apparent that modifications can be made thereto without departing from the spirit and scope of the invention. Accordingly, the invention is not limited by the foregoing description, but is only limited by the scope of the appended claims.