1. Field of the Invention
This invention is related to seismic data processing. More specifically, the invention is related to a system for processing seismic data to detect azimuthal velocity variations.
2. Description of Related Art
Seismic surveys are routinely used in the search for oil and gas reservoirs in the earth's subsurface. Seismic surveys are performed by imparting acoustic energy into the earth, either at a land surface or in a marine environment, and then detecting the reflected and refracted acoustic energy. The delay time between the imparting of the acoustic energy wave at the source location and detection of the same wave at a receiver location indicates the depth of reflecting geological interfaces.
Until recently, only two-dimensional (“2D”) seismic surveys were conducted, with the seismic source locations being collinear with a line of receivers. Recent advances in technology have enabled three-dimensional (“3D”) seismic survey data to be gathered and analyzed. Typically, in 3D surveys arrays of seismic receivers are deployed which receive reflected acoustic energy imparted at varying locations that may be specifically selected to provide a rich assortment of azimuths for common midpoints.
A technique frequently used in seismic survey analysis is AVO analysis, which is the amplitude variation with offset, and is also referred to herein as amplitude variation with incidence angle. According to the AVO approach, attributes of a subsurface interface are determined both from the normal-incidence amplitude of reflected seismic energy, and also from the dependence of the detected seismic reflections on the angle of incidence of the seismic energy at a subsurface reflecting interface relative to the vertical. In conventional AVO analysis, multiple seismic traces having a common reflection point, commonly referred to as a common mid point or common depth point(CMP or CDP) gather, are collected. From the CMP (or CDP) gather, one may derive the amplitude R of a reflected seismic wave from an interface (i.e., the “target horizon”) as a function of the angle of incidence θ from the normal according to the following relationship:
R(θ)=A+B sin2 θ.
In this case, the coefficient A is the zero-offset response (also referred to as the AVO intercept), while the coefficient B is referred to as the AVO slope, or gradient, as it is representative of the rate of change of amplitude with the square of the angle of incidence. Analysis of the AVO slope and intercept can provide indicators of interesting formations, from an oil and gas exploration standpoint. For example, variations in the A and B values from a theoretical A-versus-B trend line for the expected stratigraphic sequences can indicate the location of hydrocarbon reserves.
While simple models of subsurface geology assume azimuthal isotropy in the propagation of acoustic energy it has been observed that azimuthal anisotropy is in fact present in many survey regions, such that the velocity of acoustic energy depends upon the azimuth of the source-receiver path. If azimuthal anisotropy is present, the conventional normal moveout correction may not adequately align the seismic traces in the gather, which can result in degraded AVO analysis.
Normal moveout correction of the seismic data, both for offset-dependent delays and also for azimuthal anisotropy caused by the overburden, is therefore typically performed in producing stacked traces of improved signal-to-noise ratio for use in a 3D seismic survey. For example, U. S. Pat. No. 5,532,978 describes a method of deriving and applying azimuthal anisotropy corrections to seismic survey signals.
The detection of a preferred azimuthal direction at a reflecting interface can also provide important information regarding geological features. For example, a preferred azimuthal reflection direction can indicate the presence of aligned vertical fractures. For moderately far offsets (25°-35° incidence angles), the P-wave traveling in the plane wave parallel to aligned vertical fractures has a higher velocity than the P-wave traveling in the plane perpendicular to the fractures.
Traditionally, azimuthal velocity analysis has been performed using azimuth-sectored supergathers and picking semblance maxima at various azimuths. This reduces the problem to a series of 2-D solutions, rather than solving the complete 3-D solution. In some cases as few as two sectors may be chosen, perpendicular and parallel to the (average) principal axes of the azimuthal anisotropy. If more than two sectors are used, an ellipse is fitted to the picked velocities to give fast and slow velocity magnitudes and the azimuth of the fast velocity. These procedures suffer from several drawbacks:
Picking semblance, by hand, on azimuth sectored data is processor/interpreter dependent and extremely time consuming.
Semblance works well for data which do not show amplitude variation with offset (“AVO”), however, if the data contain significant AVO, particularly if there is a polarity reversal, semblance can fail. In this case automatic picking of semblance maxima will be erroneous.
If the subsurface has azimuthal velocity variation (“AVV”) then this will appear as an offset-dependent static viewed on offset-sorted CMP gathers. This will reduce the effectiveness of any surface consistent statics solution, thus the azimuth-sectored supergathers will most likely be contaminated with statics. This will significantly degrade the semblance analysis and may result in several semblance maxima.
The semblance is based on giving the greatest stack power. However, for AVV analysis it is the actual subsurface velocity that is of interest, not simply the velocity that gives the best stack. For instance, if a higher amplitude occurs at a particular azimuth within the sector, then the velocity at that azimuth will be picked. In addition, if those high amplitudes are at the mid to near offsets and are contaminated with residual statics then a completely erroneous velocity could give the highest semblance.
Sectoring and partial stacking of the data means that it is extremely difficult to obtain error estimates. Not only is it difficult to attribute a picking error from picking semblance, but errors due to the acquisition geometry are not represented. In any analysis of this type it is important to compute the errors associated with the obtained results. For instance a weighted least squares approach has been used to compute the errors in a technique for inverting azimuthal variation of amplitude for shear wave data. It has also been observed that the reliability of the amplitude variation with azimuth analysis has been assessed by looking for an absence of the acquisition geometry being mirrored in the anisotropy maps.
It should be noted that the description of the invention which follows should not be construed as limiting the invention to the examples and preferred embodiments shown and described. Those skilled in the art to which this invention pertains will be able to devise variations of this invention within the scope of the appended claims.
The invention comprises a system for processing seismic data to estimate time shift resulting from velocity anisotropy in the earth's subsurface. A gather of seismic data traces is formed and selected seismic data traces included in said gather within selected time windows are cross-correlated to estimate the time shift in the seismic data traces included in said gather resulting from velocity anisotropy in the earth's subsurface.
Also shown in
As is well known in the art, normal moveout corrections are typically made to traces in a common midpoint gather to correct for the additional delay time for longer offset traces, so that the travel time of each seismic signal is effectively normalized to a zero-offset travel time. In situations where the earth exhibits no azimuthal anisotropy, the azimuthal variation will not introduce a variation in the data set. However, when azimuthal anisotropy is present in seismic survey signals, the standard normal moveout correction may not adequately correct for variations in delay times for traces from source-receiver pairs of different directions.
Variations in the seismic data traces will appear to be variations in amplitude as a function of the azimuthal variation, when in reality the variations in the seismic data traces are a result of azimuthally dependent time shifts. A velocity variation with azimuth of only a few percent can cause ten milliseconds or more of time difference in the location of an event in a seismic data trace. Accordingly, velocity variations below the resolution of conventional semblance based velocity analysis can distort the amplitude variation with incidence angle (AVO) and the amplitude variation with azimuth (“AVOA”), resulting in an incorrect computation.
In addition to the standard normal moveout correction, typically, for 3-D land acquisition, deconvolution, refraction statics, and two passes of surface consistent residual statics processes are applied to the data. As a preliminary to the performance of the invention as disclosed herein, noise reduction, trace scaling and any other single-trace process may be applied, however, the application of multi-trace processing should be avoided.
Flatten Events Within a Moving Time Window
If the velocity of the seismic signals in the subsurface vary with the azimuthal direction in which the seismic signals travel, the conventional normal moveout correction will not align the traces properly, and in accordance with a preferred embodiment of the invention, the following process may be utilized to achieve proper trace alignment. FIG. 2A and
The spatial window is then moved across the CMP gather within the selected time window by one trace and a new pilot trace is generated for this new spatial window. The new pilot trace is then correlated with the new input trace within the time window. This process is continued for the selected time window, with each successive data trace within the gather being designated as the “input trace” and cross-correlated with a pilot trace which comprises a plurality of nearby traces, to complete a time window trace correlation sequence. In this way the pilot trace represents the local phase and amplitude characteristics of the data for each “input” trace. Accordingly, a decision is made in step 60 as to whether a time window trace correlation sequence has been completed. If the answer is no, steps 52, 54, 56, 58 and 60 are repeated for a successive input trace. If the answer is yes, then in step 62 a decision is made as to whether the cross-correlation sequence just completed is the first performed cross-correlation process for the time window. In one implementation of the invention, if the answer is yes, the calculated time shift that provided maximum cross-correlation between the input trace and the pilot trace for each trace is applied to each trace in step 64 and steps 52, 54, 56, 58, 60 and 62 are repeated. If the answer in step 62 is no, then in step 66 a decision is made as to whether any reflection event in the time window is substantially aligned across the traces in the gather. If the answer is no steps 64, 52, 54, 56, 58, 60, 62 and 66 are repeated. This decision in step 66 is normally based on whether or not the additional time shifts calculated for the input traces in the just completed time window trace correlation sequence are significant. If the answer in step 66 is yes, then in step 68, the total computed time shifts for each trace are stored and the traces are returned to their form they were in prior to beginning the process outlined in FIG. 3. Typically, only two iterations of the process described in steps 52, 54, 56, 58, 60, 62 and 64 are performed, but further iterations may be performed if, in the judgment of the processor, data quality may be improved by further iterations.
A decision is made in step 70 as to whether all time windows of interest in the seismic data gather have been selected, and if the answer in no a new time window is selected in step 50 and the cross-correlation procedure described above with respect to steps 52, 54, 56, 58, 60, 6264, 66, 68 and 70 is repeated for all time windows of interest. The successive time windows selected may occupy successive time positions on the seismic data traces or the time windows may overlap, depending on the quality of the data.
Following completion of the cross-correlation procedure for all time windows and all traces within each time window, in step 72, the amount of time shift which achieved the maximum correlation for each trace within each time window is applied to each trace at the center point within each time window, and, in step 74, time shifts for the remainder of the data traces are interpolated between these center points.
Each pilot trace may comprise, for example, eleven traces. The trace in the center of the selected spatial window, i.e. the sixth trace, may be designated as the “input” trace and cross-correlated with the pilot trace to obtain a time shift. The spatial trace window is then moved one trace across the gather and a new pilot trace formed. Again, the trace in the center of this window, i.e. the next trace in the gather, is designated as the “input” trace and cross-correlated with the pilot trace to obtain a time shift for that trace. At the edges of the gather, in this example the first through the sixth traces, the spatial window comprising the pilot trace may be shortened so that, if the first trace is the “input” trace, the first through sixth traces are stacked to form the pilot trace, and for the second trace, the first through the seventh traces are stacked and so on, until the full number of traces desired in the spatial window is reached (in this example, eleven). The number of traces selected to form the pilot trace may be selected on the basis of the magnitude of amplitude variation with offset in the data and the magnitude of noise such as multiple contamination.
In performing the cross-correlation procedure, amplitude and phase variations with offset, including the case where an event reverses polarity at some offset are taken into account. The pilot trace represents the “local” data characteristics. If an event reverses polarity at far offsets then the pilot trace at near offsets should not include the far offset traces. Similarly, the pilot trace at far offsets will not include traces at near offsets.
Other trace attributes, including absolute values, RMS values, or trace envelope may be used for performing the cross-correlation in addition to the raw trace reflection amplitude.
At such time as the time shifts have been applied to the seismic data traces in the gather, AVO analysis as well as AVOA analysis, such as discussed herein with reference to
Computing the Amplitude Variation with Incidence Angle (AVO), and the Amplitude Variation with Azimuth (AVOA)
It is known to those of ordinary skill in the art that amplitude variation with incidence angle (also referred to as amplitude variation with offset), as well as the amplitude variation with azimuth for a reflection from a horizontal transverse isotropic layer of the earth which is overlain with an isotropic overburden can be approximated by the following equation:
R(θ,φ)=I+G1 Sin2(θ)+G2 Sin2(θ)Cos2(φ−β) (Eq. 1)
I is the P-wave impedance contrast between the subsurface layers from which the signal is reflected;
G1 is the isotropic AVO gradient;
G2 is the azimuthal or anisotropic term; and
β (with reference to
G1 and G2 are given by:
{overscore (ρ)}, {overscore (Vp )} and {overscore (Vs )} are the average density, the average P-wave velocity and the average S-wave velocity respectively,
the average P-wave velocity divided by the average S-wave velocity,
Δδ(ν) is the change in δ(ν) across the reflecting boundary, and
Δγ is the change in the shear wave splitting parameter γ across the reflecting boundary, where:
It is known to those of ordinary skill in the art that for a linearly elastic material, each component of stress σij is linearly dependent on every component of strain εkl, where i, j, k and l are directional indices that may assume values of 1, 2 or 3. The stress-strain dependency is given by Hooke's Law:
σij=Cijklεkl
where Cijkl is the elastic modulus tensor and completely characterizes the elasticity of the medium. The relationship between δ(ν) and the elastic modulus tensor is given by:
However, without knowing β, Eq. 1, cannot be solved using a least squares approach. However, Eq. 1 can be rewritten as:
R(θ,φ)=I+{G1+G2 Cos2(θ−β)} Sin2(θ) (Eq. 4)
which can be rewritten as:
R(θ,φ)=I+{G1*+(G2*−G1*)Cos2(φ−β)} Sin2(θ) (Eq. 5)
so that G1=G1 * and G2=G2*−G1*.
Utilizing the equality:
G1*+(G2*−G1*)Cos2(φ−β)=G2* Cos2(φ−β)+G1* Sin2(φ−β) (Eq. 6)
then
R(θ,φ)=I+[G2* Cos2(φ−β)+G1* Sin2(φ−β] Sin2(θ). (Eq. 7)
It is known to those of ordinary skill in the art that:
G1* Cos2(φ−β)+G1* Sin2(φ−β)=W11 Cos2(φ)+2W12 Cos(φ)Sin(φ)+W22 Sin2(φ) (Eq. 8)
which is linear in the unknowns W11, W12 and W13, which can be related back to the unknowns G1*, G2* and β, as follows:
G2*=0.5(W11+W22+√{square root over ((W11−W22)2+4W122)}) (Eq. 9)
G
1*=0.5(W11+W22−√{square root over ((W11−W22)2+4W122)}) (Eq. 10)
Thus, combining Equations 1 with Equations 4-11, Equation 1 can be written as:
R(θ,φ)=I+[W11 Cos2(φ)+2W12 Cos(φ)Sin(φ)+W22 Sin2(φ)] Sin2(θ) (Eq. 12)
with
G1=0.5(W11+W22−√{square root over ((W11−W22)2+4W122)}) (Eq. 13)
G2=√{square root over ((W11−W22)2+4W122)} (Eq. 14)
and
Values of the reflection coefficient R(θ,φ) for specific values of the incidence angle θ and the source-receiver azimuth φ can be obtained from the recorded seismic data for each reflection event by extracting the amplitudes of the seismic data traces as a function of offset and azimuth. With reference to
where:
X is the source to receiver offset;
T0 is the zero offset two way travel time;
Vrms is the RMS velocity; and
Vint is the interval velocity at the time of interest.
In step 88, a least squares method is used to compute reflection coefficient as a function of azimuthal angle and incidence angle for the seismic traces comprising the CMP gather. Eq. 12 is solved in a straightforward least squares manner, known to those of ordinary skill in the art, for the unknowns, I, W11, W22, and W12. Values for G1 (the isotropic AVO gradient), G2 (the azimuthal or anisotropic term) and β may then be computed from values computed for W11, W22, and W12. Accordingly, it is demonstrated that Eq. 1 is linear in I, G1, G2 and the direction β.
Note that as indicated above in Eqs. 2 and 3, the derived gradients G1 and G2 are related to physical rock properties
Δδ(ν) and Δγ.
Computing The Azimuthal Velocity Variation (AVV)
Because the process outlined in
where:
T=total travel time
T0=two way zero-offset traveltime
X=offset
Vnmo(φ)=the azimuthally varying velocity as a function of the azimuth φ,
and
Accordingly, the total traveltime may be written as:
T2=T02+[W11 Cos2(φ)+2W12 Cos(φ)Sin(φ)+W22 Sin2(φ)]X2. (Eq. 19)
In step 92, Eq. 19 may be solved by using a linear least squares method known to those of ordinary skill in the art, using the time shifts picked from the cross-correlation process described with reference to
the slowest velocity is given by
and the azimuth of the slowest velocity is given by β.
Because the travel times are being fitted by the least squares solving of Eq. 19, the azimuth β that is computed using Eq. 19 is the azimuth of the greater travel time. Accordingly, if the travel time is greater, the velocity is slower. The azimuthal velocity variation may then be computed in step 94 from the following relationship:
The amount of time shift resulting from azimuthal anisotropy may be determined for each reflection event as a function of azimuthal angle and the appropriate time shift may be applied to each trace to adjust for the azimuthal time shift. At such time as the time shifts have been applied to the seismic data traces in the gather, AVO analysis as well as AVOA analysis, such as discussed herein with reference to
Computing of Errors Associated with Calculation of Time Shift Variation, Velocity Variation and Amplitude Variation with Azimuth
Errors associated with the calculation of the time shift variation with azimuth, the velocity variation with azimuth, and the amplitude variation with azimuth may be estimated utilizing a least squares approach. In step 100, the least squares approach to estimating the errors associated with the calculation of the time shift variation with azimuth is formulated in matrix notation, and may be written:
A{right arrow over (x)}={right arrow over (b)} (Eq. 21)
Where {right arrow over (b)} is a 1×N matrix (i.e. a vector) containing the data (e.g. travel times or amplitudes), A is an M×N matrix of coefficients (e.g. Sin2 (θ)) and {right arrow over (x)} is a 1×M matrix (i.e. a vector) of the parameters to be solved for. For instance for Eq. 19:
Where Tn are the observed travel times for data with source-receiver offsets Xn and azimuth φn. The matrix equation is equivalent to N simultaneous equations for M unknowns. In the example shown, M is equal to four. For a least squares formulation, N, the number of data points, must be greater than M, the number of unknowns. There are various standard numerical methods, known to those of ordinary skill in the art, for computing the M unknowns.
In step 102, the standard error in the unknowns is computed by taking the square root of the diagonals of the matrix E:
E=(ATA)−1σ2, (Eq. 22)
where the superscript T refers to the matrix transpose, σ2 is the variance (i.e. sum of squares of the differences between the data and the computed fit divided by n−4). The diagonals of the matrix E represent the errors in each unknown so that the square root of the mth diagonal element (i.e. √{square root over (Emm)})is the standard error in the mth unknown.
The least squares method allows for the computation of errors contained in the matrix A, which include both the variance (error) caused by poor data quality or random noise as well as the expected error resulting from the data distribution. Note that the elements of the matrix A are functions of the offsets and azimuths in the data, which are then combined and used to compute the errors. In addition, since the variance caused by poor data quality or random noise can be computed independently, the variance caused by poor data quality or random noise,
E=(ATA)−1, (Eq. 23)
and the error resulting from the data distribution,
E=σ, (Eq. 24)
can be separated in step 104 and either one or both compared in step 106 to the AVV or AVOA result to confirm whether or not the result has an acquisition ‘footprint’—a pattern that is caused by the acquisition geometry. Typically, this comparison is done visually, although those of ordinary skill in the art would know to make the comparison mathematically. An indicia of the accuracy of the obtained results is the absence of the acquisition geometry (such as fold, maximum offset and minimum offset) being present in the obtained results.
Errors associated with the calculation of the amplitude variation with azimuth may be calculated in an analogous manner to the calculation of the errors associated with the calculation of the time shift variation with azimuth, with the process being applied to Eq. 12 rather than Eq. 19. Because the velocity variation is calculated from the travel time variation, the errors calculated for the time shift variations with azimuth are applicable to the velocity variations with azimuth.
Computing a New Surface Consistent Statics Solution.
An azimuthal velocity variation (AVV) will cause offset dependent and time dependant statics in CMP sorted gathers. Surface consistent statics solutions for a CMP gather typically use a pilot trace comprising a stack which includes all offsets and azimuths in the gather. If AVV exists, the far offsets will not align with the near offsets and will not stack coherently, thus the pilot trace will be representative of the near offset traces, which are relatively unaffected by the AVV. To obtain a single static per trace, a cross-correlation between each trace and the pilot trace is performed over a large time window (for example, two seconds), typically around a target horizon. If the time delay due to AVV changes with time, then the pilot trace will not be representative of the far offset traces, since the far offset traces will be stretched and squeezed relative to the near offsets and thus also relative to the pilot trace. For instance, a particular azimuth may be in the slow direction (resulting in a shift to later times) at one traveltime, but change to the fast direction (resulting in a shift to earlier times) at a later traveltime. Thus when this trace is correlated with the pilot trace, a poor cross-correlation may be obtained, possibly resulting in an incorrectly picked static for that trace. Even if a static is picked that represents the ‘average’ time delay due to AVV, it will appear as noise in the surface consistent statics computation (the “SCSC”) because this computation assumes the statics are surface consistent, when in fact the time delays due to AVV are not.
In accordance with an embodiment of the present invention, as outlined in
At such time as the time shifts resulting from azimuthal anisotropy have been applied to the seismic data traces in the gather, AVO analysis as well as AVOA analysis, such as discussed herein with reference to
The process of the invention disclosed herein is most conveniently carried out by writing a computer program to carry out the steps described herein on a work station or other conventional digital computer system of a type normally used in the industry. The generation of such a program may be performed by those of ordinary skill in the art based on the processes descried herein.
The results of the calculations according this invention may be displayed with commercially available visualization software. Such software is well known to those of ordinary skill in he art and will not be further described herein. It should be appreciated that the results of the methods of the invention can be displayed, plotted or both
While the invention has been described and illustrated herein by reference to certain preferred embodiments in relation to the drawings attached hereto, various changes and further modifications, apart from those shown or suggested herein, may be made herein by those skilled in the art, without departing from the spirit of the invention, the scope of which is defined by the following claims.
This application is a divisional of U.S. patent application Ser. No. 09/855,925 filed on May 15, 2001, now U.S. Pat. No. 6,681,184.
Number | Name | Date | Kind |
---|---|---|---|
3414370 | Baumgartner | Dec 1968 | A |
3611278 | Guinzy | Oct 1971 | A |
4206509 | Ruehle | Jun 1980 | A |
4570246 | Herkenhoff et al. | Feb 1986 | A |
4613960 | Hinkley et al. | Sep 1986 | A |
4755972 | Hanson et al. | Jul 1988 | A |
4779237 | Bodine | Oct 1988 | A |
4970696 | Crews et al. | Nov 1990 | A |
5258960 | Swan | Nov 1993 | A |
5508973 | Mallick et al. | Apr 1996 | A |
5532978 | Corrigan | Jul 1996 | A |
5610875 | Gaiser | Mar 1997 | A |
5737220 | Miller | Apr 1998 | A |
5764516 | Thompson et al. | Jun 1998 | A |
5933789 | Byun et al. | Aug 1999 | A |
6041018 | Roche | Mar 2000 | A |
6061301 | Corrigan | May 2000 | A |
6263284 | Crider et al. | Jul 2001 | B1 |
Number | Date | Country | |
---|---|---|---|
20040109387 A1 | Jun 2004 | US |
Number | Date | Country | |
---|---|---|---|
Parent | 09855925 | May 2001 | US |
Child | 10689423 | US |