AU735894B2 - Method of analyzing capabilities of migration and DMO computer seismic data processing - Google Patents
Method of analyzing capabilities of migration and DMO computer seismic data processing Download PDFInfo
- Publication number
- AU735894B2 AU735894B2 AU65506/98A AU6550698A AU735894B2 AU 735894 B2 AU735894 B2 AU 735894B2 AU 65506/98 A AU65506/98 A AU 65506/98A AU 6550698 A AU6550698 A AU 6550698A AU 735894 B2 AU735894 B2 AU 735894B2
- Authority
- AU
- Australia
- Prior art keywords
- impulse
- interest
- attribute
- data
- response
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Ceased
Links
- 230000005012 migration Effects 0.000 title claims description 60
- 238000000034 method Methods 0.000 title claims description 50
- 238000013508 migration Methods 0.000 title claims description 50
- 238000012545 processing Methods 0.000 title claims description 20
- 230000004044 response Effects 0.000 claims description 88
- 238000003384 imaging method Methods 0.000 claims description 17
- 238000001914 filtration Methods 0.000 claims description 9
- 238000011158 quantitative evaluation Methods 0.000 claims description 4
- 230000000052 comparative effect Effects 0.000 claims description 3
- 238000004519 manufacturing process Methods 0.000 claims description 3
- 238000011160 research Methods 0.000 claims 1
- 230000006870 function Effects 0.000 description 19
- 230000010363 phase shift Effects 0.000 description 6
- 238000004364 calculation method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 4
- 239000000284 extract Substances 0.000 description 3
- 238000004445 quantitative analysis Methods 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 239000007787 solid Substances 0.000 description 3
- 230000000007 visual effect Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000009472 formulation Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 230000010355 oscillation Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 230000007480 spreading Effects 0.000 description 2
- 241000364021 Tulsa Species 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000003796 beauty Effects 0.000 description 1
- JLQUFIHWVLZVTJ-UHFFFAOYSA-N carbosulfan Chemical compound CCCCN(CCCC)SN(C)C(=O)OC1=CC=CC2=C1OC(C)(C)C2 JLQUFIHWVLZVTJ-UHFFFAOYSA-N 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 238000003780 insertion Methods 0.000 description 1
- 230000037431 insertion Effects 0.000 description 1
- 238000003908 quality control method Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 239000013598 vector Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Description
WO 98/43189 PCT/US98/04798 TITLE: METHOD OF ANALYZING CAPABILITIES OF MIGRATION AND DMO COMPUTER SEISMIC DATA
PROCESSING
BACKGROUND OF INVENTION: 1. FIELD OF INVENTION: The present invention relates to quality control of seismic data processing. More particularly, the present invention involves computerized quantitative evaluation of 3-D seismic imaging software, specifically migration and DMO or dip moveout software.
2. DESCRIPTION OF PRIOR ART: It has been widely known that no processing step alters the seismic data more than imaging with 3-D migration and DMO software. Imaging of seismic data in computers generally changes the position, amplitude, and phase of reflections as it converts them into images of reflectors. Many different types of migration have existed (time, depth, 2-D, 3-D, poststack and prestack) and there have been algorithms for each type (Kirchhoff, phase-shift, finite-difference, for example). There has been really only one type of DMO (time, in 2-D or but there have been different algorithms for it k-w, for example). Also, the fidelity or processing accuracy of a particular imaging algorithm is implementation-dependent, because of ever-present approximations. All of the foregoing factors have made it difficult to quantify the accuracy of particular imaging software.
Imaging software has commonly been tested by creating accurate synthetic data from some known earth model, passing these data through a migration or DMO program, and then comparing the output with the known answer from the model.. This approach could be very accurate, but it was model-dependent, and creating accurate and proper 3-D synthetic data has been difficult.
Huygens' principle allowed for decomposition and imaging df data, after migration or DMO, into what in effect was a superposition of a number of impulse responses. Alternatively then, imaging software could be understood by analyzing WO 98/43189 PCT/US98/04798 2 impulse responses. An impulse response was formed by processing a special impulse data set with a migration or DMO program which was to be evaluated. The special impulse data set contained any number of "dead" traces (zero amplitudes), and just one "live" trace, having several wavelets separated in time.
Impulse responses were simple to make, and they were actual migration or DMO operators. Impulse responses of this type were thus used to qualitatively test imaging software. An impulse response was formed by processing the impulse data set with the program being evaluated. The results were then plotted and visually examined.
Visual comparisons of impulse responses, from a variety of software being tested, allowed qualitative ranking of the software. But, this type of side-by-side visual comparison, much like beauty contests, was subjective and imprecise. Visual comparison and evaluation did not reveal how good or how bad the software really was.
SUMMARY OF INVENTION: Briefly, a computer-implemented method is performed according to the present invention to enable quantitative evaluation of the functionality or accuracy of computer software programs for imaging seismic data. The present invention thus permits evaluation of the imaging programs to be quantitative, so that it reveals precisely how good or bad the imaging software really is. An impulse data set is created in a programmed digital computer. The impulse data set contains a single live impulse data trace of one wavelet and a plurality of dead data traces. The impulse data set is then subjected to inverse filtering in the computer to compensate for natural filtering resulting from processing it with the imaging software program under evaluation.
The impulse data set is then processed in the computer using the computer imaging software program and a seismic velocity function which is representative of the seismic velocity of a portion of the earth's subsurface. This processing produces a set of processed impulse-response data traces. A set of corresponding theoretical impulseresponse events is also obtained in the computer using the seismic velocity functions.
The impulse-response formulas are derived from integral equations for migration and DMO. Output records or plots are formed of attributes such as vertical traveltime (or WO 98/43189 PCT/US98/04798 3 depth), amplitude, and phase, as predicted by these theoretical impulse responses against corresponding values measured from actual impulse-response data sets. The method of the present invention thus permits comparison of actual imaging software against theory. The plots so formed reveal just how closely imaging software agrees with theory.
For example, if an imaging operator exhibits an amplitude result, at 600 true-dip, which is ten times weaker than it should be, then the migration or DMO software will image a 60° reflector with an amplitude error of 1/10: Also if an operator extends only out to 70° dip, then the software will not image any dips steeper than 70°. The method of the present invention produces such quantitative predictions. Impulse responses can thus be used to precisely understand an imaging program's operator, and thus generally predict its imaging performance.
The method of the present invention outlined generally above for migration/DMO impulse-response software evaluation is according to the sequence set forth below: 1. An impulse data set is created that contains one live trace and any number of dead traces. The live trace is located near the center of this data set, and the data set includes enough dead traces so that the resulting impulse response does not interact with the edges of the data set. The live trace contains one or more bandlimited, zero-phase wavelets that are separated in time so that they do not interfere. The live trace is then filtered with for 3-D migration, or (i w)" 2 for 2-D migration and DMO, where o is frequency.
2. An appropriate velocity function is then created. The velocity function need not be of any exact form and to this extent may be somewhat arbitrary so long as it is realistic. Generally-variable velocity functions, V(xy,z), may be used to evaluate depth migration software, and depth-only variable velocity functions, are most appropriate for evaluating time migration and DMO software.
3. The impulse data set is then processed with migration or DMO software for selected evaluation.
WO 98/43189 PCT/US98/04798 4 4. At least two of the actual impulse-response events in the processed result are then preferably selected. An estimate is made of the vertical traveltime (or depth), amplitude, and phase attributes along these events as functions of lateral distance, measured relative to the location of the live trace in the original impulse data set.
Preferably, instantaneous phase is used for the phase attribute, but other measures may be acceptable.
Vertical traveltime (or depth), amplitude, and phase are then computed for the corresponding theoretical impulse-response events.
6. The theoretical amplitude curves are allowed to be out of scale with the imaging software by an arbitrary constant. A scale factor is computed and all amplitudes are scaled, estimated from the actual response events. One actual amplitude curve is selected, along with its theoretical counterpart. The scale factor is set to equal the ratio of the two curves at the location of the original impulse trace. This value, SCALE, is set to equal Atheory x-O Aatal, x=o, where x 0 refers to the location of the impulse trace. All actual amplitude curves are then multiplied by the SCALE value.
7. Vertical traveltime (or depth) is plotted for the actual and theoretical impulse-response events together, then the plotting process is repeated separately for the amplitude and phase attributes.
8. From the plots so formed, quantitative evaluation of the imaging software is now possible. This is done by observing how well the "actual" vertical traveltime (or depth), amplitude, and phase curves agree with their "theoretical" counterparts. In this manner, the performance of an imaging program's operator can be quantitatively assessed and evaluated.
BRIEF DESCRIPTION OF DRAWINGS: Fig. 1 is a diagram illustrating the geometry for exploding-reflectors migration.
Fig. 2 is a diagram illustrating the geometry for zero offset migration.
Fig. 3 is a diagram illustrating the geometry for common offset migration.
Fig. 4 is a diagram illustrating the geometry for common offset DMO.
Figs. 5A, 5B, and 5C are plots of an example impulse data set used according to the present invention.
WO 98/43189 PCT/US98/04798 Fig. 6 is a plot of output data traces representing a computed impulse response formed according to the present invention.
Figs. 7, 8, 9 and 10 are comparative plots of data processing test results against values obtained from an integral equation accurately representative of migration.
Fig. 11 is another plot of output data traces representing a computed impulse response formed according to the present invention.
Figs. 12, 13 and 14 are comparative plots of data processing test results against values obtained from an integral equation accurately representative of migration.
Fig. 15 is another plot of output data traces representing a computed impulse response formed according to the present invention.
Fig. 16 is a plot of vertical traveltime from the data of Fig. Fig. 17 is a plot of amplitude from the data of Fig. Fig. 18 is a schematic diagram or flow chart of the process steps of the present invention.
DESCRIPTION OF PREFERRED
EMBODIMENT
1. PHYSICAL PRINCIPLES AND UNDERLYING
THEORY
The physical principles of seismic data processing performed during imaging are expressed by a set of theoretical, integral equations for 3-D poststack or zero offset migration, 3-D prestack migration, and DMO. The equations for each of these are discussed separately below, since the theory and equations differ for each. The equations are first reviewed from the natural perspective of imaging. Thereafter, an implementation made according to the present invention is discussed, using the equations to generate theoretical impulse responses to quantitatively evaluate imaging software.
Poststack seismic data represent approximate zero-offset data, but most poststack migration software relies upon the artificial, exploding-reflectors model, as disclosed for example in Sheriff, 1989, Encyclopedic Dictionary of Exploration Geophysics, 2nd ed.: Soc. ofExpl. Geophys., Tulsa, OK. It is helpful to recall that the exploding-reflectors model assumes that all reflectors explode at time t 0, waves propagate upward with the half-velocity, V/2, and receivers, located at the surface, WO 98/43189 PCT/US98/04798 6 record the wavefield. The exploding-reflectors model does not give amplitudes consistent with the true zero-offset model, and thus the exploding-reflectors is somewhat deficient. The kinematics for both models are almost always identical, however.
The present section explains how separate theory equations are developed representing 3-D poststack migration for the exploding-reflectors and zero-offset, amplitude models. The actual poststack migration software being tested almost always conforms to exploding-reflectors, with the possible exception of poststack Kirchhoffmigration programs, these being relatively expensive, however. Thus, the present invention permits users to test processing software against that model in a quantitative manner to learn how well the migration software meets its assumptions. Alternatively, users can test against the zero-offset model to see how much error the software really makes to their approximate zero-offset data. Equation is discussed below while Equations were obtained from available public literature for the purposes of "trueamplitude" imaging. True-amplitude migration programs remove geometric spreading from the data as they image, and true-amplitude DMO programs do not (by definition) alter the amplitude of planar reflections.
A. EXPLODING REFLECTORS MIGRATION EQUATION MODEL True-amplitude, variable-velocity, exploding-reflectors migration in 3-D (Fig. 1) may be represented by: P(x,y,z, 0) f dxdyo. (1) 2t zo detQ at where Q 12 where means complex conjugate.
Q
2 1
Q
22 The integration surface is the plane, zo constant, containing the receivers. The image point is and the surface location of some receiver (representing a seismic WO 98/43189 PCTIUS98/04798 7 trace) is (xo,yo,z 0 Variable T represents traveltime for a ray that propagates between the surface location and the image point, and Q represents geometric spreading incurred on this ray. Factor a /az o is "obliquity," equaling 2cos6/V(xo,y 0 ,z 0 where 6 is the angle that the ray makes with the vertical at the surface location. The values T and aTr/zo are computed with kinematic ray tracing and Q with dynamic ray tracing. Q is a 2x2 matrix, and det Q represents the determinant of matrix Q (see below).
Kinematic ray tracing involves solving the well-known differential equations describing 3-D rays for the coordinates, direction vectors, and traveltimes along rays.
See, for example (Equations 3.1, of Cerveny, 1985, The application of ray tracing to the propagation of shear waves in complex media; in Seismic shear waves, ed. Dohr, G., Geophysical Press Ltd.) Dynamic ray tracing involves solving other well-known differential equations (Equations 5.22 and 6.27' of Cerveny, 1985) for the amplitudes everywhere along these computed rays. These differential equations are then solved numerically by integrating them with the well-known Runge-Kutta method (as discussed in Press, et al., 1992, Numerical Recipes in Fortran, 2nd ed.: Cambridge Univ. Press. The general Kirchhoff integral of Equation is validated by Berkhout and Wapenaar (1989), One-way versions of the Kirchhoff Integral; Geophysics, 54, 460-467 for variable velocity. Classic derivations of Kirchhoff integrals assume constant velocity, as described by Schneider, W.A. (1978), Integral formulation for migration in two and three dimensions: Geophysics, 43:49-76.65 B. ZERO OFFSET MIGRATION EQUATION
MODEL
True-amplitude, variable-velocity, zero-offset migration in 3-D may be represented by, P y,,t 0) f dxodyo P(xoY (2) 2 dzo at Equation lacks the eQ factor of Equation but is otherwise identical.
The simpler Equation is more accurate for migrating conventional seismic data WO 98/43189 PCT/US98/04798 8 because it does zero-offset amplitude corrections. Kirchhoff 3-D poststack migration is not common in practice, however, because other methods are more efficient. Equations and are both used in processing with the present invention to evaluate impulse responses of production 3-D poststack migration software. Equation is disclosed by Schleicher, et al. (1993) 3-D true-amplitude finite-offset migration: Geophysics, 58:1112-1126 (his equations 5 and 21).
C. COMMON OFFSET MIGRATION EQUATION MODEL True-amplitude, variable-velocity, common-offset migration in 3-D may be represented by: Ph(x, y,z, t 0) f dy dy 2 w( y, 2 zo, h, x, y,z) O 'P(y 2 h,zo t (3) 2 7 1 at
T
wcos scosOG Idet(Ns NGM) where w s jldet(NsM) det(NGM)l The integration surface is the plane, z 0 constant, containing the sources and receivers. This equation is disclosed in Schleicher, et al. (1993) 3-D true-amplitude finite-offset migration: Geophysics, 58:1112-1126 (his equations 5 and 22) for the trueamplitude migration of common-offset 3-D data. Quantity h represents the constant, half-offset value, and y, and y 2 represent midpoint coordinates, in two orthogonal directions on the recording surface (inline and crossline directions). The image point is M source location is S +h,y 2 ,z 0 and receiver location is G (y,-h,y 2 ,Zo).
It is to be noted that S and G may be defined in other ways, but the distance between S and G must always equal 2h, as required for common-offset data. Variables Ts and
T
G
represent traveltimes for rays that connect points S with M, and G with M, respectively, and Os and 0 G are angles that these rays make with the vertical at points S and G. V s is the velocity at point S. Finally, NsM and NGM are 2x2 Hessian matrices of the traveltimes T s and
T
G as defined in Schleicher, et al. (1993) 3-D true-amplitude finiteoffset migration: Geophysics, 58:1112-1126 (his equations 17 and 18). The components of these matrices are second derivatives of travel times, rs and r, with SUBSTITUTE SHEET (RULE 26) WO 98/43189 PCT/US98/04798 9 respect to lateral spatial coordinates. Equation is used to evaluate production 3-D common-offset migration software.
Traveltimes,
T
s and rG, and angles Os and 0 may be computed as indicated by Schleicher et al (1993) by kinematic ray tracing. In this case, matrices, NsM and NGM may be computed by dynamic ray tracing. These quantities may alternatively be computed by specialized 3-D traveltime-calculation methods (such as Schneider Jr., W.A, 1995 Robust and efficient upwind finite-difference traveltime calculations in three dimensions: Geophysics, 60:1108-1117). In this case, the various traveltime derivatives may be computed by finite differencing (see Press, et al., 1992, Numerical Recipes in Fortran, 2nd ed.: Cambridge Univ. Press.) the computed 3-D traveltime fields. When the velocity varies with depth only, as assumed when using Equation to evaluate 3-D prestack time migration software, then it can be seen that the weight function, w, and the traveltimes, -r and may all be computed with simple 2-D traveltime-calculation methods. An example of this is described, for example, in van Trier, J. and Symes, 1991, Upwind finite-difference calculation of traveltimes: Geophysics, 56:812-921, with the 2-D stability condition published in Schneider Jr., W.A.,1995, Robust and efficient upwind finite-diference traveltime calculations in three dimensions: Geophysics, 60:1108-1117. In this case, it can also be seen that the computation of the matrices NsM and NGM is much simplified-they are easily computed from 2-D finite-difference traveltime derivatives.
D. COMMON OFFSET DMO EQUATION
MODEL
True-amplitude, variable-velocity, common-offset DMO in 3-D or 2-D may be represented by: Py 0 th) d w(y,ty,tf,h) Pj(y 0 j h (4) where do (i)H e is the inverse Fourier transform of WO 98/43189 PCT/US98/04798 This expression has been developed from equation (31) of Black et al (1993).
Here represents temporal convolution, w represents frequency, and the y direction is the direction of the line connecting the source and receiver. Then y, is midpoint coordinate measured along this line before DMO (where half-offset equals and yo is midpoint coordinate measured along this line after DMO (where half-offset equals zero). Variables t, and t o represent time on a seismic trace at half-offsets, h and 0, respectively, and t, t. (to, T) to (1-l 1 2 where n y, Thus, y, y,(yo, Yo rl. Plus, P, and Po represent the common-offset data before and after DMO, at halfoffsets h and zero, respectively. It is to be noted that the familiar (i w) 1 filter is represented through the time-domain convolution operator 8 but the filter may still be applied in the frequency domain, as is usually preferable.
This general integral formulation for DMO is well-known, but the weight function, w, for true amplitude is attributed to Black, et al. (1993). True-amplitude imaging and dip moveout: Geophysics, 58:47-66. Their weight function, modified to accommodate Equation above is w(yoyo,,,h) (2C)to 2a( 1 27T h where a(r) 2/h2)-' Equations 2, and 3) above may be converted from 3-D to 2-D by making assumptions about how the source wavefield varies in the y-direction and that the earth model does not vary in the y-direction, then performing the y, and y, integrals analytically. The integrals then represent 2-D migration, and can be used to evaluate 2- D imaging software. This embodiment of the present invention focuses on 3-D imaging, but the methods herein described easily extend to 2-D imaging software.
WO 98/43189 PCT/US98/04798 11 2. THEORETICAL IMPULSE RESPONSES FROM THE INTEGRAL
EQUATIONS
According to the present invention, the theoretical impulse responses described above are used to evaluate actual, computed impulse responses. A theoretical impulse response is extracted from each of equations (1 through 4) above, and these analytical representations are used to generate theoretical impulse response curves.
This is done for Equation by inserting the wavefield, P(xoYz~t) 8(x xi) 68(y o yl) R(t t 1 (6) where 6 represents the dirac delta function, and R(t) is a zero-phase wavelet, with R(0) representing the central wavelet peak. Such an insertion yields the result, U I-I P(xyx,t 0) az 0 eQ, (ta (7) where T r(x, y, z; zo), ar/az o 2 cosN/V(x,y,zo), and Q Q(x,y,z; xl,yl,zo) are computed by tracing rays and performing dynamic ray tracing between the fixed surface point yl, zo) and all possible image points Equation is nonzero only for such that r= where the wavelet assumes its peak, so Equation extracts one 3-D surface from the volume of possible image points and this represents the impulse response of Equation For a fixed impulse time, the data values for a surface so defined are computed and this surface is saved in memory. Also, the amplitude and phase associated with it are similarly obtained and saved. The associated phase is developed in a manner explained below. Curves extracted from this impulse response are then plotted against curves estimated from actual impulse-response events, as described previously.
WO 98/43189 PCT/US98/04798 12 The theoretical impulse response for Equation is Equation with the det Q factor removed. The theoretical impulse response is extracted from Equation by inserting the wavefield, y2, h, zo, 8(y, (y 2 2) R(t- (8) and obtaining -1 Ph(x, y, z, t= 2 zo, h, x, y, z) R(t (9) where in Equation T, h, 2 ,zo) and TG To(x,y,z; h, 2 ,zo)are traveltimes for rays that travel between the fixed source, S h, 2 ,z 0 and receiver, G i h, 2 locations respectively, and all possible image points M M The weight function, w( 2 ,zo,h, x, y, z) is computed for each pair of these traveltimes, -r and rG. Equation is nonzero only for such that T, Tz t,, where the wavelet assumes its peak, so Equation extracts one 3-D surface, from the volume of possible image points and this represents the impulse response of Equation For fixed impulse time, the data values for a surface so defined are computed and this surface is stored in memory. The weight function and phase associated with it are similarly obtained and stored. Curves extracted from this impulse response are again plotted against curves estimated from actual impulse-response events, as described previously.
The theoretical impulse response from Equation for DMO, is extracted by inserting the wavefield, SUBSTITUTE SHEET (RULE 26) WO 98/43189 PCT/US98/04798 13 and using the previous-defined relations y, yo ir and t, to (1-12/h 2 1 2 to get Po(Yo,to,h)=w(yo, ,to,a(no)to,h)0 R(a(no)to-T) (11) a(T) to 2 2a(ro) .2 to) R (a( 0 o) to T) where a(lo) (1 12/h 2 )-1 2 Here y y Equation (10) is nonzero only for (yo, to) such that a(rio)to-T, where the wavelet assumes its peak, so Equation (10) extracts one 2-D curve, from the range of possible image points (yo, to), and this represents the impulse response of Equation For a fixed impulse time, T, the data values for a surface so defined are computed again and saved or stored along with the amplitude and phase associated with it. Curves in extracted form of this impulse response are then plotted against curves estimated from actual impulse-response events, as described previously.
3. IMPULSE-RESPONSE PHASE The -1 scaling and time derivatives in Equations and may be interpreted as frequency-domain filtering operations because a/at Fourier transforms to io, where co is frequency, and -1 scaling represents a 1800 phase shift. The convolution with the "half-derivative" operator, 81/2, in Equation also implies a frequency-domain filter of Generally, -io filters associate with 3-D imaging software and filters associate with 2-D imaging software, including all forms of DMO. Interpretation of these filters is easy: io scales the Fourier components of the data by co and introduces a constant 900 phase shift, and (io)1 2 scales the Fourier components of the data by (o) 1 /2 and introduces a constant 450 phase shift.
The analysis is simplified and these filters are accounted for by applying their inverses to the original impulse trace, when creating the impulse data set (see Section 1 above). That way all wavelets on computed impulse responses are expected to be the original zero-phase Ricker wavelets, and zero is used for the theoretical phase. Thus, SUBSTITUTE SHEET (RULE 26) WO 98/43189 PCT/US98/04798 14 the derivative operators and -1 scalings in Equations 7, 9 and 11 are ignored when computing theory curves for comparison with actual impulse responses generated by imaging software.
Imaging software will apply a time- and dip- varying, "stretch" to the wavelet.
The theory here discussed does not account for this, but the phenomenon is easily recognized and understood from the QC plots (see subsequent examples).
4. VELOCITY AND ITS IMPACT The method and theory described above, with the exception of DMO theory, is valid for variable velocity models, V V(xy,z). Wave propagation quickly becomes very complicated when V becomes complicated computed impulse responses may contain triplications or "bow ties" making any impulse-response event multivalued in the (xy) coordinates. Impulse responses may also contain diffractions for complicated V. Such things make computing the theoretical impulse-response curves difficult.
Also, ray theory may differ considerably from wave theory when V becomes complicated, which may invalidate the theoretical operators. There are acceptable solutions to these difficulties, such as computing the theoretical impulse-responses with the method of Gaussian beams (see Cerveny, 1985, The application of ray tracing to the propagation of shear waves in complex media: in Seismic shear waves, ed. Dohr, G., Geophysical Press Ltd.), or even with wave theory, although these alternatives could be expensive. With the present invention, it is preferable to restrict the velocity model. In this way ray theory still approximates wave theory well, and the theoretical curves are easily computed. Such a restriction is a smooth velocity field, or even one that varies only in depth, V For most purposes, such models are adequate for software evaluation, especially for 3-D time-migration and DMO software, the most commonlyused imaging software today.
The method may also be applied to 2-D migration software, although this is believed to be less important today than 3-D imaging, so the preferred mode is to operate in 3-D. Any wavelet can be used to construct the original impulse trace, but it is believed that a zero-phase wavelet is most appropriate, because it simplifies the WO 98/43189 PCT/US98/04798 subsequent phase analysis. The preferred zero-phase wavelet is of the type known as the Ricker wavelet (Figs. 5A, 5B, and
EXAMPLE
The general method of the present invention can be illustrated for 3-D poststack time migration using a smooth, V(z) velocity function that was derived from North Sea data. The impulse data set contained 360 traces in the x and y directions, respectively, trace spacing were Ax=12.5m and Ay=25m, representing realistic sampling. Center frequency for the impulse trace of the Ricker wavelet was 40 Hz. The impulse trace was the center trace in this data volume, and it contained three wavelets, at times 1, 2, and 3 s, respectively. The origin of coordinates was defined to be the surface location of the original impulse trace. The data set was processed with a 3-D poststack timemigration program, based upon the phase-shift algorithm (see Gazdag, 1978, Wave equation migration with the phase-shift method: Geophysics, 43:1342-1351). Impulseresponses for V(z) velocity should be symmetric about the impulse trace. Thus, single vertical slices were analyzed through the 3-D impulse response. Figure 6 shows traces extracted from the 3-D impulse-response volume, along the entire positive y-axis (the leftmost trace represents the original impulse trace, where x, Trace spacing in Figure 6 is 25 m.
A computer QC program was written that takes an actual impulse response, automatically digitizes selected impulse-response events, and tabulates vertical traveltime, amplitude, and instantaneously phase for the events, all as functions of the horizontal coordinate in the case of Figure The functions of that program are illustrated in Figure 18. The sequence of steps of how those functions are performed has been set forth above. The program computes theoretical curves, corresponding to the impulse-response events, using the appropriate equation or Equation applies with respect to Fig. 1. Finally, the program plotted the theoretical curves against those derived from the actual impulse response.
Figures 7, 8, 9, and 10 represent these plots for the data of Figure 6. Figure 7 shows vertical traveltime for the events corresponding to original impulse times of one and two seconds. Theory curves are solid, and curves derived from the actual impulse WO 98/43189 PCT/US98/04798 16 response are dashed. Figure 7 shows that vertical traveltimes agree well, although the dashed curves do not extend out as far as the solid curves. Figure 8 is the comparable amplitude plot, which reveals that the actual amplitudes oscillate up to 15% about the correct amplitudes. The actual amplitudes can be seen to agree well for small y distances (0-200m on the top curve and 0-600 m on the bottom curve), while the actual amplitudes drop off quickly with y-distance. The two dashed curves were scaled by one scale factor, SCALE (defined previously), to match them to the solid curves. Good agreement is noted, especially at the left edge of the plot. The oscillation is an implementation artifact, caused by setting a user-parameter to save costs at the expense of accuracy. The drop off however, is the natural and expected consequence of spatial aliasing. The amplitude oscillation is not desirable because it represents coherent noise that will appear in migrated images. Figure 9 is the related phase plot, where phase of both events is near zero, except for slight deviations at large y-distance because of wavelet stretch, as mentioned previously. Figure 10 shows true-dip angle versus ydistance, and this plot may be used to convert the horizontal axes of the other plots from distance to true-dip angle. For example, the amplitude curve corresponding to the original impulse time of 2 s, in Figure 8, agrees well with theory for dip angles under 250. The corresponding vertical traveltime curve in Figure 7 extends out to dip angles of almost 70°, but amplitude at these high dips has been severely attenuated.
The migration program was modified by resetting the user-parameter to maintain accuracy and incur higher costs, and then a new impulse response created.
Figure 11 represents the same vertical slice as in Figure 6, from this new run. Figure 12 is the amplitude display computed from Figure 11, showing that now the amplitudes vary smoothly, as desired. The differences between the amplitudes in Figures 8 and 12 are significant and important, but they cannot be discerned through the traditional approach of visually comhparing the impulse responses (Figures 6 and 11). This demonstrates the resolution and power of the quantitative method of the present invention.
Decreasing the trace spacing would reduce spatial aliasing and the amplitude drop off with distance (or dip angle) indicated in Figures 8 and 12. Vertical slices were WO 98/43189 PCT/US98/04798 17 again extracted from the impulse response volumes, but this time the slice direction was parallel to the x-axis, where trace sampling was 12.5 m. The slices were passed through the program, which produced the amplitude plots of Figures 13 and 14. The only difference between the plots Figures 13 and 14 and those in Figures 8 and 12, is the reduced trace spacing (from 25 m to 12.5 The plots of Figs. 13 and 14 show the improved amplitude match obtained by reducing the trace spacing, and by using appropriate user-parameter settings that maintain accuracy.
This again demonstrates the usefulness of the quantitative method of the present invention. Here it was used to understand how changes in the trace sampling change the amplitude fidelity of migrated images. Clearly it can be understood, therefore, that to image steep dips, with this velocity model, fine trace spacing (12.5 m or less) is needed. It is also to be noted that the 3-D impulse response was computed with coarser trace spacing in the y-direction than in the x-direction (25 m vs. 12.5 This is common for 3-D seismic marine processing, and through the quantitative method of the present invention described here, one can learn just how much better imaging will be in the x-direction (12.5 m) than in the y-direction (25 m).
Finally, another way of evaluating 3-D impulse responses through horizontal slices is shown. A horizontal slice, at time 1.2s, was extracted from the original impulse-response volume and plotted in Figure 15. This time slice now shows the 3-D intersection of the 2.0 s impulse-response event with the horizontal plane at time T 1.2 s (compare Figure The vertical direction is the x-direction and the horizontal direction is the y-direction in Figure 15. This slice was evaluated with the program of the present invention, and vertical traveltime, amplitude, and phase were computed as before. Figures 16 and 17 show the vertical traveltime and amplitude plots, respectively. Poststack 3-D time migration operators are azimuthally symmetric, so in time slices, impulse responses should have circular shape, with constant amplitude and phase. Figure 16 shows the theory event and the actual impulse-response event where it intersects the 1.2 s horizontal plane. Good agreement is exhibited, since theory and actual curves overlie one another. Figure 17 shows amplitude, along the bottom and top halves of the impulse-response event of Figure 15, plotted as functions ofy. It can be WO 98/43189 PCT/US98/04798 18 seen that at y=0, which is on the x-axis, strong amplitude is present, but for large y, which is near x=0, weak amplitude is present. This is again the consequence of spatial aliasing and trace spacings of dx 12.5 m and dy 25 m. These conclusions support what was seen when comparing Figures 8 and 12 for vertical slices.
From the foregoing, it can be seen how 3-D poststack time migration can be used to demonstrate various ways in which one may quantitatively analyze computed 3- D impulse responses. This information may be used to accurately predict how imaging software will perform on real data. It should be understood that 3-D prestack time migration software may be analyzed similarly, as well as 2-D or 3-D DMO. Also, depth migration software may be revaluated analogously, although the procedure becomes more difficult as the velocity function is allowed more degrees of freedom.
The foregoing disclosure and description of the invention are illustrative and explanatory thereof, and various changes in the data acquisition and processing sequences, the parameters selected and analyzed, as well as the processing sequences described, as well as in the details of the process of operation described above, may be made without departing from the scope of the present invention.
Claims (19)
1. A computer-implemented method for forming comparative plots to enable quantitative evaluation of the functionality of a computer software program for imaging seismic data, said method including the steps of: creating an impulse data set in the computer containing a single live impulse data trace and a plurality of dead data traces, said live impulse data trace containing at least one wavelet; inverse filtering said impulse data set in the computer to compensate for natural filtering resulting from processing said impulse data set in the computer with said imaging software program; processing said impulse data set in the computer using said imaging software program and a seismic velocity function representative of the seismic o: velocity of a portion of the earth's subsurface to produce a set of processed impulse-response data traces; o:oo° forming an output plot representative of an attribute of interest in said set of processed impulse-response data traces; computing a set of corresponding theoretical impulse-response events using said seismic velocity function; and forming an output plot representative of an attribute of interest in said theoretical impulse-response events for comparison with the output plot of the attribute of interest in said processed impulse-response data traces to enable quantitative assessment of the accuracy of said imaging software program.
2. The method of claim 1, further including the step of: determining the value of the attribute of interest for a selected location in said processed impulse-response data traces.
3. The method of claim 2, wherein said attribute of interest is vertical traveltime. WO 98/43189 PCT/US98/04798 1
4. The method of claim 2, wherein said attribute of interest is depth. 1 5. The method of claim 2, wherein said attribute of interest is amplitude. 1
6. The method of claim 2, wherein said attribute of interest is phase. 1 7. The method of claim 1, further including the step of: .2 determining the value of an attribute of interest for a selected location in 3 said processed computed theoretical impulse-response events. 1
8. The method of claim 7, wherein said attribute of interest is vertical 2 traveltime. 1
9. The method of claim 7, wherein said attribute of interest is depth. 1
10. The method of claim 7, wherein said attribute of interest is amplitude. 1
11. The method of claim 7, wherein said attribute of interest is phase. 1
12. The method of claim 1, further including the step of: 2 determining the value of an attribute of interest for a selected location in 3 said processed impulse-response data traces; and 4 determining the value of an attribute of interest for a selected location in said processed computed theoretical impulse-response events. 1
13. The method of claim 12, wherein said attribute of interest is vertical 2 traveltime. 1
14. The method of claim 12, wherein said attribute of interest is depth. WO 98/43189 PCT/US98/04798 21 1 15. The method of claim 12, wherein said attribute of interest is amplitude. 1 16. The method of claim 15, further including the step of: 2 computing a scale factor between the amplitude of impulse-response data 3 traces and the theoretical impulse-response events. 1 17. The method of claim 16, further including the step of: 2 applying the scale factor to the amplitude of the impulse-response data 3 traces. 1 18. The method of claim 12, wherein said parameter of interest is phase. 1 19. The method of claim 1, wherein said step of computing a set of 2 corresponding theoretical impulse-response events is performed based on an exploding 3 reflectors migration model. 1 20. The method of claim 19, wherein the exploding-reflectors migration 2 model is represented by the following mathematical expression: P(xy,t 0) ro U dy. [xaP Y 27t JJ z detQ) at 1 21. The method of claim 1, wherein said step of computing a set of 2 corresponding theoretical impulse-response events is performed based on a zero offset 3 migration model. 1 22. The method of claim 21, wherein the zero offset migration model is 2 represented by the following mathematical expression: WO 98/43189 PCT/US98/04798 22 1 1 P(x,y,z,t=O) f f dxody. *I P(xoo,zo 0 t 2x a zo at
23. The method of claim 1, wherein said step of computing a set of corresponding theoretical impulse-response events is performed based on a common offset migration model.
24. The method of claim 23, wherein the common offset migration model is represented by the following mathematical expression: Ph(x, y,z,t 0) I dy dy 2 w (l _2,Zo,h,x yz) P( 2,hzo,t The method of claim 1, wherein said step of computing a set of corresponding theoretical impulse-response events is performed based on a common offset DMO model.
26. The method of claim 25, wherein the common offset DMO model is represented by the following mathematical expression: h Po(yo,to,h)= dl w(yo, p, 7l=-h SUBSTITUTE SHEET (RULE 26)
27. The method of claim 1, wherein said step of computing a set of corresponding theoretical impulse-response events is performed based on a prestack migration model.
28. A method for evaluating the functionality of a computer software program for imaging seismic data, said method including the steps of: creating an impulse data set containing a single live impulse data trace and a plurality of dead data traces, said live impulse data trace containing at least one wavelet; inverse filtering said impulse data set to compensate for natural filtering resulting from processing said impulse data set with said imaging software program; selecting a seismic velocity function representative of the seismic velocity of a portion of the earth's subsurface; processing said impulse data set using said imaging software S program and said seismic velocity function to produce a set of processed impulse-response data traces; computing the corresponding theoretical impulse-response events using said seismic velocity function; and comparing said theoretical impulse-response events computed in step to said processed impulse-response data traces from step to assess the accuracy of said imaging software program.
29. The method of claim 28, wherein said impulse data set is a two- dimensional (2D) data set and said step of inverse filtering is done with (iw)- 1 2 The method of claim 28, wherein said impulse data set is a three- dimensional (3D) data set and said step of inverse filtering is done with
31. The method of claim 28, wherein said step of comparing said theoretical impulse-response events to said processed impulse-response data traces includes quantitatively comparing attributes of said theoretical impulse-response events to those of said processed impulse-response data traces. DATED this 24 th day of April 2001 EXXON PRODUCTION RESEARCH COMPANY WATERMARK PATENT TRADEMARK ATTORNEYS LEVEL 21 "ALLENDALE SQUARE" 77 ST GEORGES TERRACE PERTH WA 6000 0 0o
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US08/832329 | 1997-03-26 | ||
US08/832,329 US5812963A (en) | 1997-03-26 | 1997-03-26 | Method of analyzing capabilities of migration and DMO computer seismic data processing |
PCT/US1998/004798 WO1998043189A1 (en) | 1997-03-26 | 1998-03-03 | Method of analyzing capabilities of migration and dmo computer seismic data processing |
Publications (2)
Publication Number | Publication Date |
---|---|
AU6550698A AU6550698A (en) | 1998-10-20 |
AU735894B2 true AU735894B2 (en) | 2001-07-19 |
Family
ID=25261343
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
AU65506/98A Ceased AU735894B2 (en) | 1997-03-26 | 1998-03-03 | Method of analyzing capabilities of migration and DMO computer seismic data processing |
Country Status (7)
Country | Link |
---|---|
US (1) | US5812963A (en) |
EP (1) | EP0970432A4 (en) |
AU (1) | AU735894B2 (en) |
CA (1) | CA2276272A1 (en) |
MY (1) | MY132838A (en) |
NO (1) | NO994168L (en) |
WO (1) | WO1998043189A1 (en) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP3929152B2 (en) * | 1997-01-08 | 2007-06-13 | 株式会社ジオトップ | Recording medium describing method and program for calculating array Y (I) by numerical calculation by computer |
US6016462A (en) * | 1997-08-29 | 2000-01-18 | Exxon Production Research Company | Analysis of statistical attributes for parameter estimation |
US6625543B1 (en) | 2002-09-05 | 2003-09-23 | 3Dgeo Development, Inc. | Output based azimuth moveout re-gridding of seismic data |
US7133777B2 (en) * | 2004-04-26 | 2006-11-07 | Schlumberger Technology Corporation | Method for transmitting wellbore data acquired in the wellbore to the surface |
US20100067328A1 (en) * | 2008-09-17 | 2010-03-18 | Andrew Curtis | Interferometric directional balancing |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US4380059A (en) * | 1980-08-20 | 1983-04-12 | Mobil Oil Corporation | F-K Filtering of multiple reflections from a seismic section |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US4878204A (en) * | 1988-10-28 | 1989-10-31 | Geophysical Service, Inc. | Method for true-amplitude dip moveout correction |
-
1997
- 1997-03-26 US US08/832,329 patent/US5812963A/en not_active Expired - Fee Related
-
1998
- 1998-03-03 CA CA002276272A patent/CA2276272A1/en not_active Abandoned
- 1998-03-03 AU AU65506/98A patent/AU735894B2/en not_active Ceased
- 1998-03-03 EP EP98911580A patent/EP0970432A4/en not_active Withdrawn
- 1998-03-03 WO PCT/US1998/004798 patent/WO1998043189A1/en not_active Application Discontinuation
- 1998-03-24 MY MYPI98001258A patent/MY132838A/en unknown
-
1999
- 1999-08-27 NO NO994168A patent/NO994168L/en not_active Application Discontinuation
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US4380059A (en) * | 1980-08-20 | 1983-04-12 | Mobil Oil Corporation | F-K Filtering of multiple reflections from a seismic section |
Also Published As
Publication number | Publication date |
---|---|
EP0970432A4 (en) | 2002-01-02 |
WO1998043189A1 (en) | 1998-10-01 |
US5812963A (en) | 1998-09-22 |
MY132838A (en) | 2007-10-31 |
EP0970432A1 (en) | 2000-01-12 |
CA2276272A1 (en) | 1998-10-01 |
NO994168D0 (en) | 1999-08-27 |
NO994168L (en) | 1999-11-26 |
AU6550698A (en) | 1998-10-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US7376539B2 (en) | Method for simulating local prestack depth migrated seismic images | |
CA2775561C (en) | Migration-based illumination determination for ava risk assessment | |
EP1047957B1 (en) | Method of prestack 3-d migration | |
US6594585B1 (en) | Method of frequency domain seismic attribute generation | |
EP2601542B1 (en) | Method for separating independent simultaneous sources | |
EP2356492B1 (en) | Method for separating independent simultaneous sources | |
US8560242B2 (en) | Pseudo-analytical method for the solution of wave equations | |
US7525873B1 (en) | Seismic inversion of conditioned amplitude spectra | |
US6263284B1 (en) | Selection of seismic modes through amplitude characteristics | |
US7701804B2 (en) | Method and constructing a model of a heterogeneous medium described by several parameters from data expressed in different time scales | |
US20060056272A1 (en) | Methods for earth modeling and seismic imaging using interactive and selective updating | |
US20130289879A1 (en) | Process for characterising the evolution of a reservoir | |
Gjøystdal et al. | Improved applicability of ray tracing in seismic acquisition, imaging, and interpretation | |
AU735894B2 (en) | Method of analyzing capabilities of migration and DMO computer seismic data processing | |
Bourgeois et al. | 3D seismic modeling: an accurate and fast tool for reservoir characterization | |
Kaschwich et al. | Advanced ray-based synthetic seismograms | |
Milkereit et al. | A new method for the migration velocity analysis of noisy seismic data | |
Huang et al. | Surface-consistent decomposition of seismic data | |
Gibson et al. | Wavefront construction modeling and inversion of VSP data for anisotropy parameters | |
Tetalovic | Numerical and physical modeling studies of the resolution of thin beds for thickness calculation |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
TC | Change of applicant's name (sec. 104) |
Owner name: EXXONMOBIL UPSTREAM RESEARCH COMPANY Free format text: FORMER NAME: EXXON PRODUCTION RESEARCH COMPANY |
|
FGA | Letters patent sealed or granted (standard patent) | ||
MK14 | Patent ceased section 143(a) (annual fees not paid) or expired |