CN101116008A - Method and aparatus for true relative amplitude correction of seismic data for normal moveout stretch effects - Google Patents

Method and aparatus for true relative amplitude correction of seismic data for normal moveout stretch effects Download PDF

Info

Publication number
CN101116008A
CN101116008A CNA2006800045862A CN200680004586A CN101116008A CN 101116008 A CN101116008 A CN 101116008A CN A2006800045862 A CNA2006800045862 A CN A2006800045862A CN 200680004586 A CN200680004586 A CN 200680004586A CN 101116008 A CN101116008 A CN 101116008A
Authority
CN
China
Prior art keywords
stretched
wavelet
spectrum
destretched
spectra
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.)
Granted
Application number
CNA2006800045862A
Other languages
Chinese (zh)
Other versions
CN101116008B (en
Inventor
E·F·赫肯霍夫
理查德·B·阿尔弗德
哈里·L·马丁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chevron USA Inc
Original Assignee
Chevron USA Inc
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Chevron USA Inc filed Critical Chevron USA Inc
Publication of CN101116008A publication Critical patent/CN101116008A/en
Application granted granted Critical
Publication of CN101116008B publication Critical patent/CN101116008B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362Effecting static or dynamic corrections; Stacking
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/50Corrections or adjustments related to wave propagation
    • G01V2210/52Move-out correction
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/50Corrections or adjustments related to wave propagation
    • G01V2210/52Move-out correction
    • G01V2210/522Dip move-out [DMO]

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geology (AREA)
  • Environmental & Geological Engineering (AREA)
  • Acoustics & Sound (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)
  • Radar Systems Or Details Thereof (AREA)
  • Indication And Recording Devices For Special Purposes And Tariff Metering Devices (AREA)

Abstract

The present invention provides a method and apparatus for arriving at true relative amplitude destretched seismic traces from stretched seismic traces. The method compensates for offset varying reflection interference effects due to normal moveout. Stretch factors ss and also input spectra are determined for NMOR stretched seismic traces. Estimates are then made of stretched wavelet spectra from the input spectra. A destretched wavelet spectra is then obtained. Shaping correction factors are determined by taking the ratio of the destretched wavelet spectra to the stretched wavelet spectra and are applied to the input spectra of the stretched traces to arrive at a destretched trace spectra. True relative amplitude scaling factors are computed by taking the ratio of a true relative amplitude property of the destretched wavelet spectra to a corresponding true relative amplitude property of the stretched wavelet spectra. Finally, the true relative amplitude scaling factors are applied to the destretched trace spectra to arrive at true relative amplitude destretched seismic traces.

Description

Method and apparatus for true relative amplitude correction of seismic data for normal moveout stretch effects
Technical Field
The present invention relates generally to methods for analysis of seismic reflection data for the purpose of obtaining subsurface properties, and more particularly to methods for compensating for offset-varying reflection interference effects caused by normal moveout correction (NMOR) present in Common Midpoint (CMP) or Common Reflection Point (CRP) seismic trace gathers.
Background
Seismic data acquired in field surveys are typically recorded using a Common Midpoint (CMP) field technique as shown in fig. 1. Acoustic energy propagating in the form of a waveform is introduced into the earth from a series of "detonation" sources S spaced from a Common Midpoint (CMP) location. The energy from each source S strikes a Common Reflection Point (CRP) in the subsurface, a portion of which is returned to a series of spaced-apart receivers R. Using this acquisition technique, the recorded gather is characterized by having a known common surface point (CMP) or common subsurface reflection point (CRP) with increased offset between the shot and the receiver. The gathers contain recordings of desired signals at various reflection angles θ r Reflected from a Common Reflection Point (CRP) in the subsurface and/or refracted from subsurface formations. Furthermore, these traces that are recorded contain other undesirable components, such as noise, in addition to the desired signal.
The reflection coefficient is a measure of the ratio of the reflected wave amplitude to the incident wave amplitude, indicating how much energy is reflected back from the subsurface interface. The reflection coefficient is a function of the elastic properties of the subsurface formation, including changes in compressional wave velocity, shear wave velocity, and density at each interface. In reflection seismic techniques, the earth's reflection coefficients below a known coplanar position are recovered from the recorded seismic amplitude responses or seismic traces. The actual seismic disturbance from a single reflecting interface is characterized by a time-varying response or wavelet that is related to the overburden properties within the earth and also related to the reflection seismic data acquisition equipment.
Wavelets are one-dimensional pulses characterized by their amplitude, frequency and phase. The wavelet originates as a packet of energy from a source S with a particular origin time and returns to the receiver R as a series of events or reflected wavelets distributed in time and energy. This distribution depends on the velocity and density of the subsurface and the relative positions of the sources S and receivers R.
The field recorded CMP gather traces are typically processed at several steps in a processing sequence to separate the desired signal from noise, reduce the effects of time and offset varying wave trains, and align and compare the amplitude responses from the common interface. An important step in trace alignment is to apply normal moveout correction NMOR directly to the data in NMOR applications or indirectly through pre-stack imaging steps. The travel time to the common subsurface interface is calculated using the CMP gather acquisition geometry and estimates of the subsurface propagation velocity of seismic energy from the shot location to the common subsurface reflection point (CRP) and back to the receiver location to distinguish the offset from the shot to the recorder group. The travel time difference between the zero and non-zero offsets of the shot to the recorder is used to map the trace amplitude from the field recording time coordinate to the zero offset time coordinate. Whether NMOR is applied directly to CMP gather traces or NMOR is applied indirectly to pre-stack migration to produce CRP gather traces, the amplitudes of the signals in the gather of traces after NMOR application may (1) be added together to form the stacked traces; (2) Comparing them to each other in an Amplitude Versus Offset (AVO) analysis; or (3) a transformation to obtain amplitude properties from which detailed interface characteristics can be derived from changes in the amplitude response.
Figures 2A-C show the effect of wavelets and normal moveout correction (NMOR) on a single time-shot offset CMP gather consisting of identical constant amplitude reflections from a layered earth model consisting of multiple subsurface interfaces with random spacing. FIG. 2A shows an interface reflectance (RC series) CMP gather illustrating the moveout effect (time convergence) of reflections from different interfaces. FIG. 2B shows this same CMP gather, but each reflection coefficient is replaced by a wavelet whose amplitude is proportional to the reflection coefficient. For a common event, offset-dependent interference effects are manifested in the form of offset-dependent amplitudes. Fig. 2C presents the results of the data shown in fig. 2B after application of NMOR, indicating that the NMOR moveout has produced offset-varying wavelets that result in offset-varying amplitudes for equal reflection coefficients. Notably, there are variations in the reflected amplitude and bandwidth in fig. 2C due to wavelet interference prior to NMOR and NMOR correction. As a result, the amplitudes in traces from different offsets are different from each other even when the local reflectivity is the same. Therefore, these NMOR corrected amplitudes cannot be considered "true relative amplitudes".
Particularly because of ongoing efforts to explore and develop deep water, AVO analysis and inversion is now being applied to CRP trace gathers containing processed seismic amplitudes that are reflected from subsurface interfaces at reflection angles from 0 ° to 60 ° or more. Fig. 3 depicts an amplitude spectrum of a single event reflected from one interface at an angle theta. As shown in FIG. 3, applying NMOR will map the amplitude and phase spectra of seismic wavelets to be equal to COS θ r Multiplying by a series of frequencies prior to the original NMOR processing while also amplifying the amplitude spectrum of the data by a factor (COS θ) relative to zero angle reflection events r ) -1 . Thus, for a 60 ° reflection event, NMOR will shift a 40Hz amplitude response to 20Hz while simultaneously vibratingThe intensity of the amplitude spectrum doubles. Wavelets have an amplitude spectrum and also a phase spectrum. For the purposes of this specification, the term "spectrum" hereinafter refers to both the amplitude spectrum and the phase spectrum of a wavelet.
When multiple reflection events are present, NMOR stretches the interference event response differently for each offset, resulting in more complex offset-dependent interference, as shown in FIG. 2C. This NMOR stretching effect makes it difficult to directly compare the amplitudes of common events from different offset traces to each other. Another complication is that even after multifaceted processing, traces in a CMP gather will typically also contain embedded wavelets that vary with both time and offset. These wavelet variations are due to residual acquisition and propagation effects as well as NMOR stretch effects. The velocity analysis required to align near offset events with far offset events becomes problematic when the amplitude response of a common event varies significantly from near offset to far offset. Furthermore, at high frequencies, NMOR stretching will reduce the signal-to-noise ratio improvement normally expected as a result of stacking seismic traces.
U.S. patent No. 5,684,754 to Byun et al teaches a method for removing NMOR stretch from CMP gather traces. This method relies on prior knowledge of the embedded wavelet and a measure of the NMOR stretch factor from a similar analysis of the seismic data. This technique does not provide true relative amplitude compensation for NMOR induced amplitude effects and is therefore not an ideal technique for AVO analysis.
An article entitled "removing offset dependent tuning in AVO analysis" (int. Seg mtg. Page 175-178 of the 67 th annual meeting detailed abstract) by h.w.swan (1997) proposes a method for reducing NMOR stretch effects from AVO attributes (e.g., AVO intercept and gradient) calculated from NMOR processed traces that do not compensate for NMOR stretch effects. As a result, this method has the disadvantage of not being applicable to the correction of CMP or CRP gather traces.
U.S. patent No. 6,516,275 to Lazaratos describes removing wavelet stretch effects from seismic traces prior to performing operations such as stacking or computing AVO attributes. This patent teaches a method of removing stretch effects from a single trace in which a filter that varies with time and offset is used to match the response of a stretched non-zero offset trace to a zero offset (and unstretched) trace. Because this method involves matching each non-zero offset trace to a zero offset trace by designing and applying an equalization filter, the method can change the relative amplitude relationship between traces as the reflection intensity changes. This method does not restore the trace amplitude to a relative value consistent with the reflection coefficient of each trace that has been over-convolved with the wavelet prior to NMOR processing. To obtain true amplitudes, this method must assume that all traces before NMOR processing have the same wavelets as the zero offset NMOR pre-processing trace. Furthermore, this method also implicitly assumes that the time-averaged reflection coefficient for each trace has the same value for all offsets as for a zero offset trace — this assumption is generally not satisfied for large offset or reflection angle variations.
Therefore, there is a need for a method and apparatus to overcome the disadvantages of previous methods and apparatus to recover true relative amplitudes of seismic reflections between traces of different offsets. Previous methods and devices have not been able to remove stretch from seismic traces. More specifically, these methods cannot compensate for offset-varying reflection interference due to normal moveout. The present invention provides a solution to these disadvantages.
Disclosure of Invention
Seismic traces may be stretched due to direct normal moveout correction (NMOR) processing, or indirectly through pre-stack imaging processing steps. The present invention provides a method for obtaining destretched seismic traces having true relative amplitudes from such stretched seismic traces. In particular, the method compensates for offset-varying reflection interference effects due to normal moveout.
In a preferred embodiment of the method, stretch factors β are determined for NMOR stretched seismic traces, and input spectra are also determined. The stretched wavelet spectrum is then estimated from the input spectrum. A destretched wavelet spectrum is then obtained, which may be the same wavelet embedded in the seismic data traces prior to NMOR processing, or may utilize an externally specified target wavelet. The shape correction factor is then determined by taking the ratio of the destretched wavelet spectrum to the stretched wavelet spectrum. The shape correction factor is applied to the input spectrum of the stretched trace to obtain a destretched trace spectrum.
Then, the true relative amplitude scale factor is calculated by taking the ratio of the true relative amplitude characteristic of the destretched wavelet spectra to the corresponding true relative amplitude characteristic of the stretched wavelet spectra. Examples of true relative amplitude characteristics may include, but are not limited to, the zero time value of the wavelet, the area under the wavelet amplitude spectrum, or the time-averaged absolute value of the stretched trace. Finally, true relative amplitude scale factors are applied to the destretched trace spectra to obtain true relative amplitude destretched seismic traces, thereby substantially preserving the relative amplitude characteristics of the stretched wavelet spectra.
The stretched wavelet spectrum is mapped to a destretched wavelet spectrum using the similarity theorem and the stretching factor, whereby the destretched wavelet spectrum is obtained. Furthermore, the target wavelet spectrum can be modified to correct for non-white reflectance. Preferably, these stretch factors are deterministically calculated as a function of such variables as offset, time, rms velocity, interval velocity, overburden anisotropy, and geologic dip, examples of trace sets for which a stretch factor may be calculated include CDP, DMO, or CRP trace sets.
It is an object of the present invention to correct for variations in reflection amplitude and bandwidth due to trace NMOR processing so that amplitudes from traces of different offsets are substantially proportional to the subsurface reflection coefficient and are identical to each other (true relative amplitudes) when the subsurface reflection coefficients are equal.
Drawings
These and other objects, features and advantages of the present invention will become better understood with regard to the following description, appended claims and accompanying drawings where:
FIG. 1 is a schematic diagram showing the geometric layout of a Common Midpoint (CMP) gather of traces acquisition system in which energy generated by an explosive source S is reflected back from a Common Reflection Point (CRP) and recorded by a receiver R;
FIGS. 2A-C show a CMP gather of interface reflection coefficients (RC series), the same CMP gather resulting from replacing each reflection coefficient with a wavelet whose amplitude is proportional to the respective reflection coefficient, and the CMP gather after removing normal time difference (NMO);
FIG. 3 shows the stretching effect on wavelet amplitude spectra; the maximum amplitude is amplified and shifted towards lower frequencies due to NMOR processing and the wavelet is stretched;
FIGS. 4A-D show the time domain response when a single event is decomposed into multiple frequency bands, where FIG. 4A is the input pulse, FIG. 4B shows the data after decomposition into several narrow frequency bands, FIG. 4C depicts the smoothed envelopes of these narrow frequency bands, and FIG. 4D shows a reconstructed spike;
FIGS. 5A-B illustrate an example of destretching applied to a common angle CMP gather containing the effects of normal moveout stretching, in which the stretched and destretched gathers are compared for a multiple event model;
FIG. 6 is a flowchart depicting steps taken to transition from stretched seismic traces to destretched seismic traces having true relative amplitudes in accordance with the invention;
FIG. 7 is a flowchart showing more specifically the steps in FIG. 6 when time-offset trace data (upper half of the flowchart) and time-angle trace data (lower half of the flowchart) are used;
FIG. 8 is a schematic diagram showing the geometry and terminology used in calculating the stretch factor β; and
figure 9 shows the shape of a single filter used in constructing a complementary band track.
Detailed Description
1. Convolution model applied to normal moveout stretching
An accepted model for the processed seismic data amplitudes indicates that these amplitudes represent the convolution of the wavelet excited by the seismic source with the subsurface reflection coefficients derived from the changes in the elastic properties of the subsurface interfaces. The time domain form of this model shows that the processed seismic data traces can be represented as a convolution of the wavelet with the earth reflection coefficient function prior to NMOR processing, or as:
d(t,t 0j ,Δt j )=∫w(τ)r(t-t 0j -Δt j -τ)dτ (1)
wherein d (t, t) 0j ,Δt j ) Is a seismic data trace, w (T) is a wavelet, r (t-t) 0j -Δt j ) Is the sum of the discrete subsurface reflection coefficient, δ (pulse) functions, given by:
r(t-t 0j -Δt j )=∑r j δ(t-t 0j -Δt j ). (2)
in this expression, r j δ(t-t 0j -Δt j ) Is of size r j T is the time domain representation of the jth reflection coefficient of (a) 0j Is zero offset time, Δ t j Is the time offset before NMOR processing.
In the frequency domain, this convolution is expressed as the product of the wavelet and the Fourier transform of the function of the earth's reflection coefficients, or as
Figure A20068000458600151
In this expression, D (f), W (f) and Σ r j e -2πif(t0j-Δtj) Data d (t, t) before NMOR processing 0j ,Δt j ) Wavelet function w (t) (hereinafter referred to simply as wavelet) before NMOR processing and reflection coefficient function r (t-t) 0j -Δt j ) The fourier transform of (d).
The key insight with respect to removing the effect of NMOR on stretched traces comes from this understanding: the effect of NMOR on the reflection coefficient function is different from its effect on the wavelet. First, the model supporting the application of NMOR to the reflection coefficient function is the time of reflection coefficient and the zero offset time t 0j Aligned without changing their size r j . Conceptually, NMOR changes the travel time between the reflection coefficients of the top and bottom surfaces of a layer to the vertical travel time through the layer. Using the migration theorem (r.m. bracewell, fourier transform and its applications, mcGraw-Hill,1965, P104-107), NMOR transforms the reflection coefficient function as follows:
Figure A20068000458600152
wherein r is j δ(t-t 0j -Δt j ) Is to apply the NMOR time offset Δ t j Before zero offset time t 0j Size r j ⊃ represents Σ r j δ(t-t 0j ) To frequencyFourier transform of the domain. Second, the effect of applying NMOR to the wavelet in the time domain is to stretch the wavelet by a time-varying factor β, which is related to the travel-time difference to an interface, which is dependent on the offset. The similarity theorem (R.M. Bracewell, fourier transform and its applications, mcGraw-Hill,1995, P101-104) can be used to describe when pulling a constantThe effect of the extension factor β in the frequency domain when applied to the time domain wavelet w (t) preceding NMOR is shown in fig. 3. The similarity theorem gives:
w(t)stretchw(t/β)⊃W(fβ)|β| (5)
where W (t/β) is a wavelet stretched in time, ⊃ shows the fourier transform of W (t/β) to the frequency domain, W (f β) | β | is the fourier transform of the stretched wavelet, β is the stretch factor, and β typically varies with time and offset in a CMP or CRP trace. Conceptually, the similarity theorem controls the effect of NMOR on wavelets, while the offset theorem controls the effect of NMOR on reflection coefficients. NMOR changes the time difference between the reflection coefficients but does not change their size.
Finally, because of the different effects of NMOR on reflection coefficients and wavelets, the Fourier transform Dnmo (f) of a typical NMOR corrected seismic trace yields a spectrum of the form:
Figure A20068000458600161
it is shown that the spectral components of the wavelet before NMOR are scaled by the sum of the earth reflection coefficients that are frequency dependent, in addition to being shifted and scaled by the stretch factor β.
In summary, the convolution model predicts that applying NMOR to CMP or CRP gather traces will cause frequency-dependent changes in both the amplitude and phase spectra of the convolved wavelet, but will not change the amplitude of the subsurface reflection coefficients. To remove the stretch caused by NMOR, the wavelet spectrum needs to be estimated to obtain and compensate for the stretch factor β, to produce an output trace amplitude that scales proportionally with the local subsurface reflection coefficient.
Principle of true relative amplitude correction of stretched traces
The invention provides a method and a device for correcting true relative amplitude of seismic traces aiming at NMOR stretching influence. The de-NMOR stretching (de-stretching) is performed in the frequency domain, i.e., the wavelet stretch factor β is estimated and the factor β is removed from the estimated value W (f β) | β | of the stretched wavelet, so that the resulting data (convolution of the embedded wavelet and the reflection coefficient function) becomes:
Figure A20068000458600171
wherein the wavelet stretch factor beta is equal to the interface reflection angle theta r W (f) is the destretched or externally specified target wavelet.
The stretched wavelet spectrum W (f β) | β | can be estimated using trace spectral averaging techniques. These techniques reduce the effect of the earth's reflection coefficient on the trace amplitude spectrum to the effect of a constant scalar multiplier. Spectral averaging techniques rely on earth reflection coefficients averaged over a frequency band taken from a large time gate to be statistically constant or of known spectral shape, and include fourier transformation of the in-band average of the trace spectrum or of the finite delayed autocorrelation function of the trace. Application of the reflection whitening filter described herein can remove the effect of the non-whitened reflection spectrum on the estimated wavelet spectrum. In a preferred embodiment of the invention, the effect of varying earth reflection coefficients in the trace spectrum is minimized by averaging the amplitude and phase sample values within a frequency band. Ideally, the number of spectral samples participating in the averaging is the inverse of the expected maximum duration of the wavelet, and is also larger than the ratio of the transform time window to the wavelet duration (preferably, at least 10 spectral samples should be included in each band). The trace amplitude and phase spectra averaged from these bands will be the estimate Wn (f) of the stretched wavelet spectrum scaled by the reflection coefficient, in the form shown below:
W n (f)=(r c )W(fβ)|β| (8)
where W (f) is the wavelet spectrum before NMOR, r c Is a frequency independent, time and offset dependent scaling factor that depends on the local subsurface reflection coefficient. As in fig. 7Alternatively, wavelet spectral estimates can be generated from frequency domain averaging of small time gated fourier transform samples, as shown in steps 130B, 132B, and 140B.
"true relative amplitude" destretching is defined as an operation that preserves the zero-valued values of the stretched (pre-NMO) wavelet convolved with an isolated reflection coefficient. If it is stretched wavelet zero time value w (t/beta) converter t=0 And removing zero time value w (t/beta) of stretched wavelet t=0 Constrained such that:
w(t)| t=0 =w(t/β)| t=0 = constant (9)
This occurs.
In the frequency domain, this constraint would require:
∑W(f i β)|β|=∑W(f i ) (10)
wherein W (f) i ) Are discrete values of the stretched wavelet fourier transform.
If the phase spectrum of the wavelet is assumed to be zero, or only the amplitude spectrum of the wavelet is known, then there may be another definition of true relative amplitude destretching of the data: it is an operation of not changing the area under the small wave amplitude spectrum so that
∑(W(f i β)|β| W*(f i β)|β|) 1/2 =∑(W(f i )W*(f i )) 1/2 (11)
Wherein W (f) i ) And W (f) i β) are each W (f) i ) And W (f) i β), (W [ (] f) ] i )W(f i )) 1/2 Are amplitude spectrum samples. For the case where the stretch factor β or wavelet varies over time, this true relative amplitude condition should be satisfied on an instantaneous basis.
In the most general form of this method, tracks stretched by a normal moveout are subjected to a discrete Fourier transform and separated into tracks of overlapping frequency bands having a center frequency f i . These frequency bands should be complementary, fromSuch that the sum of each frequency band in the frequency domain is equal to the frequency domain representation of the input track data. Alternatively, the wavelet spectral estimate may be generated by averaging samples of the fourier transform of the small time gate in the frequency domain.
By utilizing the principle, the real relative amplitude destretching of the seismic traces is completed in 3 steps. First, the wavelet spectrum estimate is preferably non-white reflection coefficient spectrum corrected to produce a stretched wavelet spectrum Ws (f) corrected for time varying reflection coefficients i ,t-t 0j ) Given by:
W s (f i ,t-t 0j )=W n (f i ,t-t 0j )R(f i )/R(f i β) (12)
wherein Wn (f) i ,t-t 0j ) Is at time t-t 0j And frequency f i Stretched wavelet spectrum of R (f) i ) Is at a frequency f i Beta is a stretch factor. Then, the stretched (pre-NMO) target wavelet spectrum O (f) will be removed using the similarity theorem i ) Is defined as:
O(f i )=W s (f i /β,t-t 0j )/|β| (13)
wherein Ws (f) i ,t-t 0j ) Is a stretched wavelet spectrum (wavelet spectrum after NMOR) corrected by the reflection coefficient, and β is a stretching factor.
Alternatively, for use of a user-specified target wavelet Wd (f) i ) Desired output in the case of true relative amplitude destretching (instead of NMOR's wavelets)The wavelet is defined as:
O(f i )=W d (f i ,t-t 0j )R(f i ) (14)
wherein R (f) i ) Is user-specified at frequency f i Non-white vertical time reflection coefficient spectrum. Second, for the frequency at the center f i And time t-t 0j Output amplitude Aout (f) for each input amplitude value of (1) i ,t-t 0j ) Given by:
A out (f i ,t-t 0j )=A in (f i ,t-t 0j )O(f i )/W s (f i ,t-t 0j ). (15)
in general, the desired and estimated wavelet samples may be complex, having both amplitude and phase components.
Finally, at each time t-t 0j The sum of the output band data ∑ A out (f i ,t-t 0j ) Is constrained so that the spectrally lower area products of the input and desired output wavelets are equal. For minimum and maximum output frequency f min And f max The desired output wavelet is stretched to have a spectral value of Ws (f) i ,t-t 0j ) From beta to -1 f min To beta -1 f max Summing at O (f) i ) From f min To f max The sum is applied as follows to produce true relative amplitude trace data A with the stretch removed out (t- t 0j ):
A out (t-t 0j )=∑A out (f i ,t-t 0j )∑(W s (f i ,t-t 0j )W s *(f i ,t-t 0j )) 1/2 /∑(O(f i )O*(f i )) 1/2 (16)
Detailed procedure for implementing true relative amplitude correction of stretched traces
FIG. 6 is a flow diagram of a preferred embodiment of a "destretch" method for obtaining true relative amplitude destretched seismic traces from stretched seismic traces. FIG. 7 expands the steps in FIG. 6 to identify each step by a corresponding profile descriptor at the top of the flowchart. Two examples of de-stretching methods are shown for processing (a) time-offset data in the upper flow path and (b) time-angle data in the lower flow path. The steps in each flow path are identified by the addition of the letter "a" to the reference number in the time-offset data flow path, and by the letter "b" in the time-angle data. In both example flow paths, the destretching of traces is accomplished by replacing the wavelet spectra in stretched seismic traces with the true relative amplitude NMOR pre-wavelet spectra to obtain true relative amplitude destretched seismic traces.
Fig. 6 summarizes the general steps taken in the present de-stretching process. Stretched seismic traces are acquired in step 110. A stretch factor β is then determined for each stretched seismic trace in step 120. The input spectrum for each stretched trace is determined in step 130. An estimate of the stretched wavelet spectrum is determined from the input spectrum in step 140. In step 150, the destretched wavelet spectrum is obtained by calculation using the stretching factor β determined in step 120 or the desired wavelet spectrum after destretching input by the user.
In step 160, a shape correction factor is determined by comparing the destretched wavelet spectra to the stretched wavelet spectra. The shape correction factor is applied to the input spectrum of the stretched trace in step 170 to obtain a destretched trace spectrum. The true relative amplitude scale factor is calculated in step 180 by taking the ratio of the true relative amplitude characteristic of the destretched wavelet spectrum to the corresponding true relative amplitude characteristic of the stretched wavelet spectrum. Finally, true relative amplitude scale factors are applied to the destretched trace spectra in step 190 to obtain destretched seismic traces of true relative amplitude, thereby substantially preserving true relative amplitude characteristics of the stretched wavelet spectra.
Referring now to fig. 7, the above-described de-stretching method may be performed in the time-offset domain or the time-reflection angle domain. The method of de-stretching applied to the time-offset domain will be described first below.
A. Method for destretching in time-offset domain
Stretched seismic trace data is acquired at step 110. The input data is processed in such a way that the true relative amplitudes are maintained. Again, the data is time corrected so that the time now seen in each trace represents the time as if the source and receiver positions were coincident, i.e., the time at zero offset. This time correction can be done either by normal moveout correction or by pre-stack imaging processing. The source to receiver distance needs to be known and is constant over the length of the trace. This time correction process is an operation for recording tracks one by one, and any track order can be used. However, the preferred track data organization is to combine data gathers into Common Depth Point (CDP), common Midpoint (CMP), or Common Reflection Point (CRP) gathers. In addition, the phase spectrum of an embedded waveform in each track is constant over the available bandwidth of the data (with sufficiently strong amplitude). The auxiliary data required to calculate the stretch factor β are the velocity model as a function of position and time, information about the traces relative to the velocity model, and the distance from the source to the receiver (or offset).
Then, in step 120A, a stretch factor β is determined for each stretched seismic trace. The stretch factor beta is defined as the reflection angle theta of the track r The inverse of the cosine.
Fig. 8 shows an explosive source S and a receiver R. The energy wavelet emitted by the explosive source S will generally be refracted as it traverses the subsurface formations. Interfaces containing Common Reflection Points (CRP) may be oriented at geological dips θ d from the normal. The angle of reflection alpha of a straight ray is shown as the angle between the normal of the interface at the CRP and the ray extending from the CRP to the receiver R. Angle of reflection theta r Is the angle from which the wavelet was transmitted from the CRP. As shown in FIG. 8, the wavelet will refract further after reflection from the CRP and before being received at receiver R due to changes in the subsurface velocity.
A deterministic estimate of the stretch factor beta as a function of time and offset can be generated from the rms velocity function and the interval velocity at the interface.
β=(1-sin 2 θ r ) -1/2 (17)
Wherein
sinθ r =xv i cosθ d /(v a 2 (t 0 2 +(xcosθ d ) 2 /v a 2 ) 1/2 ) (18)
v a =v rms (1+2εsin 2 α(2-sin 2 α)) 1/2 (19)
sin 2 α=x 2 /(x 2 +v rms t 0 ) 2 ) (20)
For a CDP gather, the CDP gather,
sin 2 θ r =(xv i cosθ d ) 2 /(v a 4 (t 0 2 +(xcosθ d ) 2 /v a 2 )) (18a)
for Dip Moveout (DMO) gathers:
sin 2 θ r =(xcosθ d ) 2 /((t dmo 2 v a 4 /v i 2 )+(cos 2 θ d v a 2 /v i 2 ))
for PSTM data:
sin 2 θ r =(xcos 2 θ d ) 2 /((t mig 2 v a 4 /v i 2 )+(cos 4 θ d v a 2 /v i 2 ))
wherein, sin θ r Is the sine of the angle of reflection at the interface and x is the offset of the detonation point to a set of receivers. V i Is at t 0 Velocity of subterranean zone, V rms Is the rms velocity, t, of the overlying layer 0 Zero offset travel time, t dmo Zero offset travel time, t, for DMO gather mig Is the zero offset time of the migrated gather, ε is the overburden anisotropy parameter, θ d Is the geological dip at the interface.
An input trace spectrum for the stretched trace is determined in step 130A, where the stretched trace is decomposed into a plurality of narrowband traces, as shown in fig. 4B. The property of these narrowband tracks is that summing them will produce the original input track.
Fig. 9 shows the design of the band filters required to produce the band traces.
Number of frequency bands N b Preferably, it is selected as the minimum value satisfying the following conditions:
N w ≤N b ≤N t /10 (21)
wherein N is w Is the number of time samples in the expected wavelet, N t Is the number of samples in the fourier transform time window. N is a radical of b And N t A typical selection value of/10 implies a time averaging window of several hundred milliseconds.
The preferred practice is to select the number of frequency bands to be approximately equal to the number of sample points in the embedded waveform.
The user-defined band filter has the characteristic of passing data completely at the center frequency and tapering off from the center frequency to completely rejecting data. The adjacent bands are fully complementary such that the all-pass center frequency is a zero-pass frequency in the adjacent bands. The corresponding frequency samples of each frequency filter are multiplied by the fourier transform of the input trace. The result of this operation will form filtered fourier transform frequency data. The product is then inverse fourier transformed to produce the desired narrowband trace, as shown in fig. 4B.
The stretched wavelet spectrum estimate is caused by the input spectrum of step 130A in step 140A.
FIGS. 4A-D show the use of N satisfying the foregoing constraints b The single event shown in FIG. 4AUsing the envelope or rms time average of each of these time traces, according to equation (8), for each center frequency fi and time t-t 0j Generating a time and frequency independent post-NMOR wavelet spectral estimate Wn (f) i ,t-t 0j ):
W n (f i ,t-t 0j )=(r c )W(f i β,t-t 0j )|β| (22)
Wherein beta is the stretch factor, r c Is a local reflection coefficient scalar.
The amplitude estimate for each band trace is calculated by taking the filtered fourier transform frequency data (which is an intermediate product from step 130A) and applying a 90 phase rotation, and then calculating an inverse fourier transform to form a hilbert transform for each corresponding band trace from step 130A. To form the trace envelope corresponding to the time samples of each narrowband trace, the Hilbert transform pair is squared together and summed, with the square root being the final sum. This constitutes the trace envelope (fig. 4C) which is filtered to reduce the effect of noise. The resulting trace amplitude set represents an estimate of the pre-correction wavelet amplitude spectrum determined time-by-time.
The destretched output spectrum is obtained in step 150A. In a preferred manner, the stretched wavelet spectra determined in step 140A are summed using the similarity theorem equation (5)The estimated value is corrected. This correction uses the stretch factor (β) (equation (17-20)) calculated in step 120A, and the stretched wavelet spectrum estimate (equation (8)) determined in step 140A, to correct the stretch in the track data. For each time and average frequency in step 140A, the amplitude of the stretched wavelet is interpolated at a frequency that is the product of β and the current center frequency. The interpolated amplitude is multiplied again by the inverse of the stretch factor (β) -1 ) Becomes a corrected sample of the destretched output spectrum. The product of this operation is completedThe product is a corrected wavelet amplitude spectrum that constitutes a function of time.
As an alternative to calculating a destretched output spectrum, a user-defined target waveform may be used to generate a destretched output spectrum. Preferably, the destretched output spectra are characterized by high and low frequency characteristics that extend over the widest offset trace over the available data frequencies.
In step 162A, a shape (or corrected shape) correction factor is determined by taking the ratio of the destretched wavelet spectrum to the stretched wavelet spectrum, using equation (12-14). A stability factor is preferably added to the denominator (stretched wavelet spectrum) to prevent the possible divide-by-zero condition from occurring. If necessary, calculations can be performed in step 160A using equations (12) and (14) to make optional non-white reflection coefficient corrections to the destretched wavelet spectra. This correction is achieved by modeling the subsurface earth reflection coefficient amplitude spectrum. This modeled earth reflectance spectrum is decomposed into a stretched wavelet spectrum in step 150A.
In step 170A, a shape (or corrected shape) correction factor is applied to the input spectrum of the stretched track using equation (15), resulting in a destretched track spectrum.
The true relative amplitude scale factor is calculated in step 180A. The ratio of the corresponding true relative amplitude characteristics of the stretched wavelet spectrum and the destretched wavelet spectrum is taken to determine the true relative amplitude scale factor. Examples of such true relative amplitude characteristics include the zero time value of the wavelet, the area under the wavelet amplitude spectrum, or the time-averaged absolute value of the stretched trace.
In step 192A, the true relative amplitude scale factor and shape (or corrected shape) factor are applied to the destretched trace spectra using equation (16) to obtain true relative amplitude destretched seismic traces, thereby substantially preserving the true relative amplitude characteristics of the stretched wavelet spectra. In step 192A, the input spectrum of step 130A is corrected to form a corrected output spectrum containing the traces of each band. The corresponding samples given in time and center frequency (i.e., correction factors) output by step 180A and the corresponding samples given in time and center frequency (i.e., shape (or corrected shape) factors) from step 162A are multiplied by those samples of step 130A (the input data represented as frequency) to form corrected frequency bands. In step 195A, the corrected output spectrum is converted to destretched seismograms. The corrected individual band traces of step 190A are added together to form an unstretched output trace.
Finally, at each time t-t 0j The sum of the output band data ∑ A out (t i ,t-t 0j ) Constrained so that the areas under the input and desired output wavelet spectra are equal. For desired minimum and maximum output frequency f min And f max Is output wavelet, the value of the stretched spectrum is W s (f i ,t-t 0j ) From β f min To β f max Sum at O (f) i ) From f min To f max Summed and applied as follows:
A out (t-t 0j )=∑A out (f i ,t-t 0j) ∑(W s (f i ,t-t 0j )W s *(f i ,t-t 0j )) 1/2 /∑(O(f i )O*(f i )) 1/2 (16)
thereby producing true relative amplitude trace data A of destretched out (t-t 0j )。
B. Time-angle domain de-stretching method
Stretched seismic traces are acquired in step 110. As previously described, the stretched traces may have been processed by normal moveout correction (NMOR) or by a prestack imaging process.
Then, a stretch factor beta is determined for the stretched seismic traces in step 120B,the stretch factor β is defined as the reflection angle θ r The inverse of the cosine of (b), due to this reflection angle theta r A constant value on a track, so that the reflection angle theta can be directly determined by searching the track head section of each track r
The input spectrum is determined for the stretched trace window in step 130B. The data is partitioned into overlapping windows. Each window is again converted to the frequency domain by taking its fourier transform and calculating its amplitude spectrum in step 132B.
An estimate of the stretched wavelet (amplitude) spectrum is calculated in step 140B. These estimates are determined by smoothing the input spectrum of step 132B by low pass filtering to form estimates of the stretched wavelet spectrum.
The destretched wavelet spectrum is obtained in step 150B. Preferably, the destretched wavelet spectra are calculated as follows: the stretched wavelet spectrum estimate determined in step 140B is now corrected using the similarity theorem. This correction corrects the stretch in the data using the stretch factor (β) determined in step 120B and the estimate of the stretched wavelet spectrum. The amplitude of each frequency sample is interpolated at step 140B, at a frequency that is the product of β and the current frequency. The amplitude of the samples at that frequency is again scaled by a stretch factor β to a destretched wavelet (amplitude) spectrum.
As an alternative to calculating a destretched output spectrum, a user-defined target waveform may be used to generate a destretched output spectrum. Preferably, this destretched output spectrum is characterized by high and low frequency characteristics that overlay the available data frequencies on the maximum offset trace.
If necessary, additional adjustments are made to the track correction factor in step 160B to accommodate the non-white reflection coefficient. This correction is calculated by modeling the amplitude spectrum of the subsurface earth's reflection coefficients. This modeled earth reflection coefficient spectrum is assigned to the destretched wavelet spectrum determined in step 150B.
In step 162B, a shape correction factor (or corrected shape factor) is determined by comparing the estimated values of the stretched wavelet spectra to the desired destretched wavelet spectra. The quotient (or ratio) of the destretched wavelet spectra and each corresponding frequency sample of the stretched wavelet spectra is calculated to obtain a shape correction factor. Preferably, a stabilizing factor is added to the denominator to prevent division by zero.
The shape correction factor (or adjusted shape factor) is applied to the stretched trace input spectrum in step 170B to obtain a destretched trace spectrum.
Then, the true relative amplitude scale factor is calculated in step 180B. The ratio of the corresponding true relative amplitude characteristics of the stretched and unstretched wavelet spectra is calculated to determine a true relative amplitude scale factor. As before, these true relative amplitude characteristics may include the zero time value of the wavelet, the area under the wavelet amplitude spectrum, or the time averaged absolute value of the stretched trace.
True relative amplitude scale factors and shape correction factors (or corrected shape correction factors) are applied to the destretched trace spectra in step 190B to obtain true relative amplitude destretched seismic traces, thereby substantially preserving true relative amplitude characteristics of the stretched wavelet spectra. More specifically, corresponding frequency samples of both the shape correction factor (or corrected shape correction factor) and the true relative amplitude scale factor are applied to the input spectrum of the stretched trace window to obtain an output spectrum with the stretched trace window removed. In step 192B, a destretched time window of data is formed by taking the inverse fourier transform of the product of step 190B. Traces are reconstructed by summing the windows in step 194B to form destretched seismic traces.
FIG. 5 is an example of destretching a common angle CMP gather. The set of stretched input traces on the left side of the graph has residual-of-view moveout induced by interference and reflection angle-dependent amplitude variations, regardless of the fact that the reflection coefficient for each event is equal at different reflection angles. In the destretched data on the right side of the figure, all events have almost exactly the same amplitude at different reflection angles and reduced (reduced) induced apparent residual moveout.
The present invention also includes a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to obtain true relative amplitude destretched seismic traces from stretched seismic traces, the method comprising the steps of:
a. acquiring a stretched seismic trace;
b. determining a stretching factor beta of the stretched seismic traces;
c. determining an input frequency spectrum of the stretched seismic trace;
d. determining an estimate of the stretched wavelet spectrum from the input spectrum;
e. obtaining a destretched wavelet spectrum;
f. determining a shape correction factor by taking off the ratio of the stretched wavelet spectrum to the stretched wavelet spectrum;
g. applying a shape correction factor to the input spectrum of the stretched trace to obtain a destretched trace spectrum;
h. calculating a true relative amplitude scale factor by taking the ratio of the true relative amplitude characteristic of the stretched wavelet spectrum to the corresponding true relative amplitude characteristic of the stretched wavelet spectrum; and
i. true relative amplitude scaling factors are applied to the destretched trace spectra to obtain true relative amplitude destretched seismic traces, thereby substantially preserving true relative amplitude characteristics of the stretched wavelet spectra.
While in the foregoing specification this invention has been described in relation to certain preferred embodiments thereof, and many details have been set forth for purpose of illustration, it will be apparent to those skilled in the art that the invention is susceptible to alteration and that certain other details described herein can vary considerably without departing from the basic principles of the invention.

Claims (19)

1. A method for obtaining true relative amplitude destretched seismic traces from stretched seismic traces, the method comprising the steps of:
a. acquiring a stretched seismic trace;
b. determining a stretching factor beta of the stretched seismic traces;
c. determining an input frequency spectrum of the stretched seismic trace;
d. determining an estimate of the stretched wavelet spectrum from the input spectrum;
e. obtaining a destretched wavelet spectrum;
f. determining a shape correction factor by taking off the ratio of the stretched wavelet spectrum to the stretched wavelet spectrum;
g. applying a shape correction factor to the input spectrum of the stretched trace to obtain a destretched trace spectrum;
h. calculating a true relative amplitude scale factor by taking the ratio of the true relative amplitude characteristic of the stretched wavelet spectrum to the corresponding true relative amplitude characteristic of the stretched wavelet spectrum; and
i. true relative amplitude scaling factors are applied to the destretched trace spectra to obtain true relative amplitude destretched seismic traces, thereby substantially preserving true relative amplitude characteristics of the stretched wavelet spectra.
2. The method of claim 1, wherein:
the following numerical formula is used:
W(f)=W n (f/β)/|β|
in the formula
W (f) = destretched wavelet spectra;
W n (f) = stretched wavelet spectra;
β = the wavelet stretch factor of the wavelet,
the stretched wavelet spectra are mapped with a stretch factor to destretched wavelet spectra using the similarity theorem, resulting in destretched wavelet spectra.
3. The method of claim 2, wherein the stretched wavelet spectra are modified for non-white reflection coefficient correction using the following mathematical formula:
W s (f)=W n (f)R(f)/R(f/β)
in the formula
W s (f) = stretched wavelet spectra after non-white reflection coefficient spectra correction;
W n (f) = stretched wavelet spectra;
r (f) = reflection coefficient spectrum at normal incidence specified by user; and
β = wavelet stretch factor.
4. The method of claim 1, wherein:
deriving a destretched wavelet spectrum from the user-specified target wavelet spectrum, wherein:
W(f)=W d (f);
in the formula
W (f) = destretched wavelet spectra;
W d (f) = user-specified target wavelet spectrum.
5. The method of claim 4, wherein the target wavelet spectra are modified for non-white reflection coefficient correction using the following numerical formula:
W(f)=W d (f)R(f)
in the formula
W (f) = destretched wavelet spectra;
W d (f) = user specified target wavelet spectrum;
r (f) = user-specified vertical time reflection coefficient spectrum.
6. The method of claim 1, wherein:
the stretch factor β is calculated deterministically as a function of offset, time, rms velocity and interval velocity.
7. The method of claim 1, wherein:
the stretch factor β is deterministically calculated as a function of offset, time, rms velocity, interval velocity, and at least one of overburden anisotropy and geological dip.
8. The method of claim 7, wherein:
the stretched seismic traces are CDP gathers, and beta -1 Is calculated according to the following mathematical expression:
β=(1-sin 2 θ r ) -1/2 =1/cosθ r
in the formula
sin 2 θ r =(x v i cosθ d ) 2 /(v a 4 (t 0 2 +(xcosθ d ) 2 /v a 2 ));
v a =v rms (1+2εsin 2 α(2-sin 2 α)) 1/2
sin 2 α=x 2 /(x 2 +v rms t 0 ) 2 );
In the formula
θ r Is the angle of reflection at the interface;
x is the offset from the shot to the receiver group;
v i is at t 0 The subterranean zone velocity of;
v rms is the rms velocity of the overburden;
t 0 zero shot-detection travel time;
ε is the upper cladding anisotropy parameter;
θ d is the geological dip at the interface.
9. The method of claim 7, wherein:
the stretched seismic traces are Dip Moveout (DMO) gathers, and beta -1 Is calculated according to the following mathematical expression:
β=(1-sin 2 θ r ) -1/2
in the formula
sin 2 θ r =(x cosθ d ) 2 /((t dmo 2 v a 4 /v i 2 )+(cos 2 θ d v a 2 /v i 2 ))
v a =v rms (1+2εsin 2 α(2-sin 2 α)) 1/2
sin 2 α=x 2 /(x 2 +v rms t 0 ) 2 );
In the formula
θ r Is the angle of reflection at the interface;
x is the offset from the shot to the receiver group;
v i is at t 0 The subterranean zone velocity of;
v rms is the rms velocity of the overburden;
t dmo zero offset travel time of the DMO gather is obtained;
ε is the upper cladding anisotropy parameter;
θ d is the geological dip at the interface.
10. The method of claim 7, wherein:
the stretched seismic traces are CRP gathers, and beta -1 Is calculated according to the following mathematical expression:
β=(1-sin 2 θ r ) -1/2
in the formula
sin 2 θ r =(xcos 2 θ d ) 2 /((t mig 2 v a 4 /v i 2 )+(cos 4 θ d v a 2 /v i 2 ))
v a =v rms (1+2εsin 2 α(2-sin 2 α)) 1/2
sin 2 α=x 2 /(x 2 +v rms t 0 ) 2 );
In the formula
θ r Is the angle of reflection at the interface;
x is the offset from the shot to the receiver group;
v i is at t 0 The subterranean zone velocity of;
v rms is the rms velocity of the overburden;
t mig is the zero offset time of the migrated gather;
ε is the upper cladding anisotropy parameter;
θ d is the geological dip at the interface.
11. The method of claim 1, wherein:
the true relative amplitude characteristic maintained in the true relative amplitude destretched trace is selected from one of the following characteristics: the zero time value of the wavelet, the area under the wavelet amplitude spectrum, and the time-averaged absolute value of the stretched trace.
12. A method of calculating true relative amplitude scale factors, the method comprising the steps of:
determining a stretched wavelet spectrum for the stretched record gather after NMOR processing;
determining a destretched wavelet spectrum; and
the ratio of the true relative amplitude characteristic of the destretched wavelet spectra to the corresponding true relative amplitude characteristic of the stretched wavelet spectra is taken to calculate a true relative amplitude scale factor.
13. A method of deriving a destretched wavelet spectrum from a stretched gather of records, the method comprising:
using the following mathematical formula
W(f)=W n (f/β)/|β|
In the formula
W (f) = destretched wavelet spectra;
W n (f) = stretched wavelet spectra;
β = the wavelet stretch factor of the wavelet,
the stretched wavelet spectra are mapped with a stretch factor to destretched wavelet spectra using the similarity theorem.
14. A method of calculating a stretch factor β for a stretched seismic trace, wherein:
the stretch factor β is calculated deterministically as a function of offset, time, rms velocity and interval velocity.
15. The method of claim 14, wherein:
the stretch factor β is deterministically calculated as a function of offset, time, rms velocity, interval velocity, and at least one of overburden anisotropy and geological dip.
16. The method of claim 14, wherein:
the stretched seismic traces are CDP gathers, and β -1 Is calculated according to the following mathematical expression:
β=(1-sin 2 θ r ) -1/2
in the formula
sin 2 θ r =(x v i cosθ d ) 2 /(v a 4 (t 0 2 +(xcosθ d ) 2 /v a 2 ));
v a =v rms (1+2ε sin 2 α(2-sin 2 α)) 1/2
sin 2 α=x 2 /(x 2 +v rms t 0 ) 2 );
In the formula
θ r Is the angle of reflection at the interface;
x is the offset from the shot to the receiver group;
v i is at t 0 The subterranean zone velocity of;
v rms is the rms velocity of the overburden;
t 0 zero shot-detection travel time;
ε is the overburden anisotropy parameter;
θ d is the geological dip at the interface.
17. The method of claim 14, wherein:
the stretched seismic traces are Dip Moveout (DMO) gathers, and beta -1 Is calculated according to the following mathematical expression:
β=(1-sin 2 θ r ) -1/2
in the formula
sin 2 θ r =(x cosθ d ) 2 /(t dmo 2 v a 4 /v i 2 )+(cos 2 θ d v a 2 /v i 2 ))
v a =v rms (1+2εsin 2 α(2-sin 2 α)) 1/2
sin 2 α=x 2 /(x 2 +v rms t 0 ) 2 );
In the formula
θ r Is the angle of reflection at the interface;
x is the offset from the shot to the receiver group;
v i is at t 0 The subterranean zone velocity of;
v rms is the rms velocity of the overburden;
t dmo zero offset travel time of the DMO gather;
ε is the upper cladding anisotropy parameter;
θ d is the geological dip at the interface.
18. The method of claim 14, wherein:
the stretched seismic traces are CRP gathers, and beta -1 Is calculated according to the following mathematical expression:
β=(1-sin 2 θ r ) -1/2
in the formula
sin 2 θ r =(x cos 2 θ d ) 2 /((t mig 2 v a 4 /v i 2 )+(cos 4 θ d v a 2 /v i 2 ))
v a =v rms (1+2εsin 2 α(2-sin 2 α)) 1/2
sin 2 α=x 2 /(x 2 +v rms t 0 ) 2 );
In the formula
θ r Is the angle of reflection at the interface;
x is the offset from the shot to the receiver group;
v i is at t 0 The subterranean zone velocity of;
v rms is the rms velocity of the overburden;
t mig is the zero offset time of the migrated gather;
ε is the upper cladding anisotropy parameter;
θ d is the geological dip at the interface.
19. A program storage device readable by a machine, tangibly embodying a program of instructions executable by the machine to perform obtaining true relative amplitude destretched seismic traces from stretched seismic traces using the steps of:
a. acquiring a stretched seismic trace;
b. determining a stretching factor beta of the stretched seismic traces;
c. determining an input frequency spectrum of the stretched seismic trace;
d. determining an estimate of the stretched wavelet spectrum from the input spectrum;
e. obtaining a destretched wavelet spectrum;
f. determining a shape correction factor by taking off the ratio of the stretched wavelet spectrum to the stretched wavelet spectrum;
g. applying a shape correction factor to the input spectrum of the stretched trace to obtain a destretched trace spectrum;
h. calculating true relative amplitude scale factors by taking the ratio of true relative amplitude characteristics of the stretched wavelet spectra to corresponding true relative amplitude characteristics of the stretched wavelet spectra; and
i. true relative amplitude scaling factors are applied to the destretched trace spectra to obtain true relative amplitude destretched seismic traces, thereby substantially preserving true relative amplitude characteristics of the stretched wavelet spectra.
CN2006800045862A 2005-02-12 2006-02-09 Method and aparatus for true relative amplitude correction of seismic data for normal moveout stretch effects Expired - Fee Related CN101116008B (en)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US11/056,640 2005-02-12
US11/056,640 US7230879B2 (en) 2005-02-12 2005-02-12 Method and apparatus for true relative amplitude correction of seismic data for normal moveout stretch effects
PCT/US2006/004688 WO2006088729A2 (en) 2005-02-12 2006-02-09 Method and aparatus for true relative amplitude correction of seismic data for normal moveout stretch effects

Publications (2)

Publication Number Publication Date
CN101116008A true CN101116008A (en) 2008-01-30
CN101116008B CN101116008B (en) 2011-11-23

Family

ID=36916933

Family Applications (1)

Application Number Title Priority Date Filing Date
CN2006800045862A Expired - Fee Related CN101116008B (en) 2005-02-12 2006-02-09 Method and aparatus for true relative amplitude correction of seismic data for normal moveout stretch effects

Country Status (8)

Country Link
US (1) US7230879B2 (en)
EP (4) EP2680041A1 (en)
CN (1) CN101116008B (en)
AU (1) AU2006214552B2 (en)
CA (2) CA2849608A1 (en)
EA (1) EA013009B1 (en)
NO (1) NO20074608L (en)
WO (1) WO2006088729A2 (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102016494B (en) * 2008-02-28 2012-12-19 斯塔特伊公司 Improved interferometric methods and apparatus for seismic exploration
CN103217707A (en) * 2012-01-18 2013-07-24 中国石油天然气集团公司 Method of directly extracting longitudinal wave time domain converted wave angle gather
CN103842853A (en) * 2011-08-05 2014-06-04 沙特阿拉伯石油公司 Correcting time lapse seismic data for overburden and recording effects
CN107132577A (en) * 2017-07-05 2017-09-05 西安交通大学 A kind of seismic attenuation method of estimation changed based on area under spectrum

Families Citing this family (62)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060227662A1 (en) * 2005-03-04 2006-10-12 Richard Foy Stretch free trace processing using block move sum and phase-based move out corrected data
US7848220B2 (en) * 2005-03-29 2010-12-07 Lockheed Martin Corporation System for modeling digital pulses having specific FMOP properties
US20060256657A1 (en) * 2005-05-11 2006-11-16 Prism Seismic, Inc. Method for improving the time-depth tie of well log data and seismic data
US7405997B2 (en) * 2005-08-11 2008-07-29 Conocophillips Company Method of accounting for wavelet stretch in seismic data
ES2652413T3 (en) * 2006-09-28 2018-02-02 Exxonmobil Upstream Research Company Iterative inversion of data from simultaneous geophysical sources
CN102112894B (en) * 2008-08-11 2015-03-25 埃克森美孚上游研究公司 Estimation of soil properties using waveforms of seismic surface waves
US20100149912A1 (en) * 2008-12-17 2010-06-17 Luren Yang System and method for reducing signature variation of seismic sources
US8537638B2 (en) * 2010-02-10 2013-09-17 Exxonmobil Upstream Research Company Methods for subsurface parameter estimation in full wavefield inversion and reverse-time migration
US8223587B2 (en) * 2010-03-29 2012-07-17 Exxonmobil Upstream Research Company Full wavefield inversion using time varying filters
RU2558013C2 (en) 2010-05-05 2015-07-27 Эксонмобил Апстрим Рисерч Компани Q tomography method
US8694299B2 (en) 2010-05-07 2014-04-08 Exxonmobil Upstream Research Company Artifact reduction in iterative inversion of geophysical data
US8756042B2 (en) 2010-05-19 2014-06-17 Exxonmobile Upstream Research Company Method and system for checkpointing during simulations
US8767508B2 (en) 2010-08-18 2014-07-01 Exxonmobil Upstream Research Company Using seismic P and S arrivals to determine shallow velocity structure
US8437998B2 (en) 2010-09-27 2013-05-07 Exxonmobil Upstream Research Company Hybrid method for full waveform inversion using simultaneous and sequential source method
WO2012047378A1 (en) 2010-09-27 2012-04-12 Exxonmobil Upstream Research Company Simultaneous source encoding and source separation as a practical solution for full wavefield inversion
US20120099396A1 (en) * 2010-10-22 2012-04-26 Chevron U.S.A. Inc. System and method for characterization with non-unique solutions of anisotropic velocities
RU2568921C2 (en) * 2010-11-11 2015-11-20 Шлюмбергер Текнолоджи Б.В. Pulse form inversion and frequency-domain whitening-inversion of seismic data
WO2012074592A1 (en) 2010-12-01 2012-06-07 Exxonmobil Upstream Research Company Simultaneous source inversion for marine streamer data with cross-correlation objective function
US8892413B2 (en) 2011-03-30 2014-11-18 Exxonmobil Upstream Research Company Convergence rate of full wavefield inversion using spectral shaping
WO2012134609A1 (en) 2011-03-31 2012-10-04 Exxonmobil Upstream Research Company Method of wavelet estimation and multiple prediction in full wavefield inversion
US9140812B2 (en) 2011-09-02 2015-09-22 Exxonmobil Upstream Research Company Using projection onto convex sets to constrain full-wavefield inversion
US9176930B2 (en) 2011-11-29 2015-11-03 Exxonmobil Upstream Research Company Methods for approximating hessian times vector operation in full wavefield inversion
EP2795366A4 (en) * 2011-12-20 2016-06-22 Conocophillips Co Critical reflection illumination analysis
EP2823335A4 (en) 2012-03-08 2016-01-13 Exxonmobil Upstream Res Co Orthogonal source and receiver encoding
CN102636813B (en) * 2012-05-04 2013-03-20 郭平 Normal moveout stretch cutting method for processing geophysical exploitation seismic data
US10577895B2 (en) 2012-11-20 2020-03-03 Drilling Info, Inc. Energy deposit discovery system and method
US10317548B2 (en) 2012-11-28 2019-06-11 Exxonmobil Upstream Research Company Reflection seismic data Q tomography
US20140254321A1 (en) * 2013-03-08 2014-09-11 Chevron U.S.A. Inc. Methods and systems for determining clathrate presence and saturation using simulated well logs
US10853893B2 (en) * 2013-04-17 2020-12-01 Drilling Info, Inc. System and method for automatically correlating geologic tops
US10459098B2 (en) * 2013-04-17 2019-10-29 Drilling Info, Inc. System and method for automatically correlating geologic tops
SG11201508195PA (en) 2013-05-24 2015-12-30 Exxonmobil Upstream Res Co Multi-parameter inversion through offset dependent elastic fwi
US10459117B2 (en) 2013-06-03 2019-10-29 Exxonmobil Upstream Research Company Extended subspace method for cross-talk mitigation in multi-parameter inversion
US9702998B2 (en) 2013-07-08 2017-07-11 Exxonmobil Upstream Research Company Full-wavefield inversion of primaries and multiples in marine environment
WO2015026451A2 (en) 2013-08-23 2015-02-26 Exxonmobil Upstream Research Company Simultaneous sourcing during both seismic acquisition and seismic inversion
US10036818B2 (en) 2013-09-06 2018-07-31 Exxonmobil Upstream Research Company Accelerating full wavefield inversion with nonstationary point-spread functions
US9910189B2 (en) 2014-04-09 2018-03-06 Exxonmobil Upstream Research Company Method for fast line search in frequency domain FWI
SG11201608175SA (en) 2014-05-09 2016-11-29 Exxonmobil Upstream Res Co Efficient line search methods for multi-parameter full wavefield inversion
US10185046B2 (en) 2014-06-09 2019-01-22 Exxonmobil Upstream Research Company Method for temporal dispersion correction for seismic simulation, RTM and FWI
AU2015280633B2 (en) 2014-06-17 2018-07-19 Exxonmobil Upstream Research Company Fast viscoacoustic and viscoelastic full-wavefield inversion
US20160018541A1 (en) * 2014-07-18 2016-01-21 Chevron U.S.A. Inc. System and method for rock property estimation of subsurface geologic volumes
US10838092B2 (en) 2014-07-24 2020-11-17 Exxonmobil Upstream Research Company Estimating multiple subsurface parameters by cascaded inversion of wavefield components
US10422899B2 (en) 2014-07-30 2019-09-24 Exxonmobil Upstream Research Company Harmonic encoding for FWI
US10386511B2 (en) 2014-10-03 2019-08-20 Exxonmobil Upstream Research Company Seismic survey design using full wavefield inversion
EP3210050A1 (en) 2014-10-20 2017-08-30 Exxonmobil Upstream Research Company Velocity tomography using property scans
US11163092B2 (en) 2014-12-18 2021-11-02 Exxonmobil Upstream Research Company Scalable scheduling of parallel iterative seismic jobs
US9989660B2 (en) * 2014-12-18 2018-06-05 Pgs Geophysical As Correcting time shifts
US10520618B2 (en) 2015-02-04 2019-12-31 ExxohnMobil Upstream Research Company Poynting vector minimal reflection boundary conditions
US10317546B2 (en) 2015-02-13 2019-06-11 Exxonmobil Upstream Research Company Efficient and stable absorbing boundary condition in finite-difference calculations
MX2017007988A (en) 2015-02-17 2017-09-29 Exxonmobil Upstream Res Co Multistage full wavefield inversion process that generates a multiple free data set.
CA2985738A1 (en) 2015-06-04 2016-12-08 Exxonmobil Upstream Research Company Method for generating multiple free seismic images
US10838093B2 (en) 2015-07-02 2020-11-17 Exxonmobil Upstream Research Company Krylov-space-based quasi-newton preconditioner for full-wavefield inversion
CN108139499B (en) 2015-10-02 2020-02-14 埃克森美孚上游研究公司 Q-compensated full wavefield inversion
MX2018003495A (en) 2015-10-15 2018-06-06 Exxonmobil Upstream Res Co Fwi model domain angle stacks with amplitude preservation.
US10908316B2 (en) 2015-10-15 2021-02-02 Drilling Info, Inc. Raster log digitization system and method
US10768324B2 (en) 2016-05-19 2020-09-08 Exxonmobil Upstream Research Company Method to predict pore pressure and seal integrity using full wavefield inversion
US10067252B2 (en) 2016-07-25 2018-09-04 Chevron U.S.A. Inc. Methods and systems for identifying a clathrate deposit
CN107807390B (en) * 2016-09-09 2019-08-23 中国石油化工股份有限公司 The processing method and system of seismic data
US10852450B2 (en) * 2017-05-03 2020-12-01 Saudi Arabian Oil Company Refraction-based surface-consistent amplitude compensation and deconvolution
WO2020227327A1 (en) * 2019-05-06 2020-11-12 RS Energy Group Topco, Inc. CA System and method for well interference detection and prediction
CN112305612B (en) * 2019-07-23 2022-05-20 中国海洋石油集团有限公司 High-resolution complex spectrum decomposition time-frequency space domain amplitude variation correction method along with offset distance
CN112394393B (en) * 2019-08-13 2023-03-21 中国石油化工股份有限公司 CRP gather data volume reconstruction method
CN113267808B (en) * 2020-02-14 2023-04-25 中国石油天然气集团有限公司 Amplitude compensation method and device

Family Cites Families (26)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB1306588A (en) 1969-12-23 1973-02-14 Messerschmitt Boelkow Blohm Nozzle assembly for jet propulsion engines
US3613071A (en) * 1969-12-24 1971-10-12 Petty Geophysical Eng Co Simultaneous dual seismic spread configuration for determining data processing of extensive seismic data
US3705382A (en) * 1970-02-26 1972-12-05 Petty Geophysical Eng Co Methods for improved deconvolution of seismic or similar data
US3735337A (en) * 1971-07-22 1973-05-22 Amco Prod Co Minimizing multiple reflecting
US4503527A (en) * 1981-03-30 1985-03-05 Mobil Oil Corporation Method for enhancement of the signal-to-noise ratio in seismic reflection signals
US4747054A (en) * 1984-11-28 1988-05-24 Conoco Inc. Method for non-linear signal matching
US4894809A (en) * 1985-05-23 1990-01-16 Mobil Oil Corporation Method for bin, moveout correction and stack of offset vertical seismic profile data in media with dip
US4759636A (en) * 1985-12-16 1988-07-26 Amoco Corporation Method and system for real-time processing of seismic data
US4759836A (en) 1987-08-12 1988-07-26 Siliconix Incorporated Ion implantation of thin film CrSi2 and SiC resistors
US5197039A (en) 1988-03-29 1993-03-23 Shell Oil Company Methods for processing seismic data
US4964103A (en) * 1989-07-13 1990-10-16 Conoco Inc. Three dimensional before stack depth migration of two dimensional or three dimensional seismic data
US4995007A (en) * 1989-12-21 1991-02-19 Shell Oil Company Method for processing seismic data
US5083297A (en) 1990-06-26 1992-01-21 Chevron Research And Technology Company Method of improving the seismic resolution of geologic structures
US5150331A (en) * 1991-03-25 1992-09-22 Amoco Corporation Method for enhancing seismic data
US5500832A (en) * 1993-10-13 1996-03-19 Exxon Production Research Company Method of processing seismic data for migration
US5596546A (en) * 1994-12-21 1997-01-21 Western Atlas International, Inc. Spatially distributed signal sampling method
US5684754A (en) * 1995-12-13 1997-11-04 Atlantic Richfield Company Method and system for correcting seismic traces for normal move-out stretch effects
US5661697A (en) * 1995-12-18 1997-08-26 Atlantic Richfield Company Method and apparatus for detection of sand formations in amplitude-versus-offset seismic surveys
US5784334A (en) * 1996-03-13 1998-07-21 Atlantic Richfield Company Method and system for detecting hydrocarbon reservoirs using amplitude versus offset analysis of seismic signals
US5862100A (en) * 1996-05-28 1999-01-19 Atlantic Richfield Company Method and system for detecting hydrocarbon reservoirs using statistical normalization of amplitude-versus-offset indicators based upon seismic signals
CA2186297A1 (en) 1996-09-24 1998-03-25 R. Daniel Wisecup Spatially distributed signal sampling method
US6021379A (en) * 1997-07-29 2000-02-01 Exxon Production Research Company Method for reconstructing seismic wavefields
GB9726928D0 (en) * 1997-12-19 1998-02-18 Geco Prakla Uk Ltd Method of stacking seismic signals
GB9813851D0 (en) * 1998-06-27 1998-08-26 Geco Prakla Uk Ltd Seismic data acquisition and processing method
MY129095A (en) * 2001-02-13 2007-03-30 Exxonmobil Upstream Res Co Method for spectral balancing of near-and far-offset seismic data.
US6798714B1 (en) * 2003-04-30 2004-09-28 Kelman Technologies Inc. Method of performing stretch-free normal moveout (NMO) and stacking of seismic traces

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102016494B (en) * 2008-02-28 2012-12-19 斯塔特伊公司 Improved interferometric methods and apparatus for seismic exploration
CN103842853A (en) * 2011-08-05 2014-06-04 沙特阿拉伯石油公司 Correcting time lapse seismic data for overburden and recording effects
CN103842853B (en) * 2011-08-05 2017-05-31 沙特阿拉伯石油公司 Time shift geological data is corrected for coating and record effect
CN103217707A (en) * 2012-01-18 2013-07-24 中国石油天然气集团公司 Method of directly extracting longitudinal wave time domain converted wave angle gather
CN107132577A (en) * 2017-07-05 2017-09-05 西安交通大学 A kind of seismic attenuation method of estimation changed based on area under spectrum

Also Published As

Publication number Publication date
EP2511735B1 (en) 2014-05-14
CA2597598C (en) 2014-11-04
WO2006088729A2 (en) 2006-08-24
EP2680040A1 (en) 2014-01-01
EP1849025A4 (en) 2012-04-04
CA2597598A1 (en) 2006-08-24
US20060193205A1 (en) 2006-08-31
EA013009B1 (en) 2010-02-26
CA2849608A1 (en) 2006-08-24
AU2006214552B2 (en) 2011-10-06
EP2680041A1 (en) 2014-01-01
EP1849025B1 (en) 2013-01-16
EP1849025A2 (en) 2007-10-31
NO20074608L (en) 2007-09-11
EP2511735A2 (en) 2012-10-17
WO2006088729A3 (en) 2007-07-05
US7230879B2 (en) 2007-06-12
AU2006214552A1 (en) 2006-08-24
CN101116008B (en) 2011-11-23
EA200701709A1 (en) 2008-08-29
EP2511735A3 (en) 2013-03-20

Similar Documents

Publication Publication Date Title
CN101116008A (en) Method and aparatus for true relative amplitude correction of seismic data for normal moveout stretch effects
EP1360635B1 (en) Method for spectral balancing offset seismic data
US7477992B2 (en) Method for combining seismic data sets
AU2002243981A1 (en) Method for spectral balancing seismic data
US7433265B2 (en) Converted wave energy removal from seismic data
CA2357604A1 (en) Method for identifying and removing multiples from seismic reflection data
US7616524B1 (en) Wavelet based intercept attribute for seismic exploration
WO1998037437A1 (en) A method of processing seismic data signals
US5684754A (en) Method and system for correcting seismic traces for normal move-out stretch effects
US6393366B1 (en) Deconvolution of seismic data based on fractionally integrated noise
AU710992B2 (en) Method for removing non-geologic amplitude variations from seismic data
AU2011224012A1 (en) Method and apparatus for true relative amplitude correction of seismic data for normal moveout stretch effects
Lines et al. A model-based comparison of modern velocity analysis methods
Carrion et al. Source wavelet and its angular spectrum from plan-wave seismograms
Simmons Jr et al. A matched-filter approach to impedance estimation
Zhou et al. Surface diffraction noise attenuation for marine seismic data processing with mathematical morphological filtering
Przebindowska et al. Acoustic imaging away from the borehole using coherence-based migration
Guerrero et al. VTI anisotropy parameter estimation in the tau-p domain: An example from the North Sea
Huang Uncertainty analysis of P-wave receiver functions
Amundsen et al. Estimation of phase velocities and Q-factors from VSP data

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20111123

Termination date: 20180209