Time-Space Varying Spectra for Seismic Processing

• US 20090292475A1
• Filed: 05/21/2008
• Published: 11/26/2009
• Est. Priority Date: 05/25/2007
• Status: Active Grant

• Current Assignee: Prime Geoscience Corporation
• Original Assignee: Prime Geoscience Corporation
• Inventors: Taylor, James D., Alam, Aftab
• Granted Patent: US 8,185,316 B2
• Time in Patent Office: Days
• US Class Current: 702/14
• CPC Class Codes: G01V 1/34 Displaying seismic recordings or visualisation of seismic data or attributes
• CPC Class Codes: G01V 1/36 Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy

Abstract

A method and visualization apparatus for spectral analysis of time-and-space varying signals enables high resolution investigation of 3D seismic data for the exploration of oil and gas. The method extrapolates multi-resolution short windows into an average long window then computes its FFT. Extrapolation uses the continuity relationship between data inside and outside of short windows. Applications of the method are illustrated with graphical screen 3D volume displays of amplitude spectra, dip and azimuth, curvature and faults (figure below). Aside from high resolution these displays improve the productivity of a seismic interpreter.

Claims

1. A method for estimating amplitude and phase spectra of data consisting of time-varying sensor signals sampled continuously in time with uniform high frequency-resolution, comprising:

(a) obtaining a representation of said data sampled at an input (Δ

T) time-rate;

(b) specifying an output sample time-rate (Δ

T′

) that is an integer multiple of Δ

T;

(c) specifying a data analysis window Δ

t that is larger than the output sample time rate (Δ

T′

), wherein the analysis window is an integer multiple of Δ

T;

(d) specifying a frequency sample rate (Δ

f);

(e) computing a time window (2/Δ

f) that would produce the frequency sample rate and that is an integer multiple of the input time rate (Δ

T);

(f) generating a model that extrapolates data inside the data analysis window to fill the frequency sample rate (2/Δ

f) window, wherein

(g) the model uses a continuity relationship between data within and data outside of the data analysis (Δ

t) window,

(h) wherein model has the property that forward and backward extrapolated values decrease in amplitude with respect to time from edges of the data analysis (Δ

t) window,

(j) capturing data within the data analysis (Δ

t) window at a start of data time;

(k) extrapolating in forward and backward directions the captured data to fill the (2/Δ

f) time window;

(l) computing fast Fourier transform (FFT) coefficients for the data in the (2/Δ

f) time window;

(m) incrementing by the output sample time (Δ

T′

) from the start time of data capture for a next FFT computation;

(n) repeating (k) through (m) until an end of data time is reached; and

at least one of

(o) saving in a storage medium the FFT coefficients sampled at (Δ

T′

,Δ

f) intervals, and

(p) saving in the storage medium a transformation of the FFT coefficients into amplitude and phase spectra sampled at (Δ

T′

,Δ

f) intervals.

2. The method of claim 1 wherein:

(a) the mathematical model is a statistical Linear Prediction Filter (LPF) model;

(b) the continuity relationship is an autocorrelation function with a maximum lag, Lmax, of said data in a time window of length, (2/Δ

f), that is symmetrically centered relative to a moving data time window of length Δ

t; and

(c) the forward extrapolation is done with repeated application of a minimum phase filter derived via Levinson recursion from the autocorrelation function and backward extrapolation is done with repeated application of the time reversed version of the minimum phase filter.

3. The method of claim 2 further comprising:

(a) computing the auto-correlation function of the signal in a (2/Δ

f)-window to determine statistical correlation between the signal within and outside of a signal analysis window,

(b) using the auto-correlation function to compute LP coefficients via Levinson recursion and linearly combining said coefficients with analysis window in a recursive loop for extrapolations in the forward direction and then repeating said linear combination with the reversed sequence of said coefficients for extrapolations in the backward direction,

(c) specifying a plurality of multi-resolution data analysis windows, where each such multi-resolution data analysis window has a different length, the length being in a range from said Δ

t to a selected fraction of said (2/Δ

f);

(d) initializing the start of spectral computation time;

(e) capturing data within said plurality of multi-resolution signal analysis windows, each centered at the spectral computation time;

(f) computing a sequence of LP coefficients for each multi-resolution window;

(g) extrapolating captured data in each multi-resolution window in forward and backward directions to fill said (2/Δ

f)-window;

(h) averaging data over plurality of multi-resolution extrapolated (2/Δ

f)-windows to compute an average (2/Δ

f)-window;

(j) computing FFT coefficients for the average (2/Δ

f)-window;

(k) incrementing spectral computation time by Δ

T′

,

(l) repeating (f) through (k) until an end of spectral computation time;

(m) at least one of storing said FFT coefficients as 2-dimensional arrays sampled at

(Δ

T′

, Δ

f) intervals, saving in said storage medium a transformation of said FFT coefficients into amplitude and phase spectra as 2-dimensional arrays sampled at (Δ

T′

, Δ

f) intervals, and

(p) displaying the amplitude and phase spectra as 2-dimensional images on an interactive graphical screen,

9. The method of claim 1 further comprising displaying the amplitude and phase spectra as 2-dimensional images on an interactive graphical screen.

4. A method to estimate time delay and dispersion between seismic signals recorded after propagation between spatially separated receivers, where said time delay is an intercept at a reference frequency and dispersion is the slope of a linear relationship of timed delay with respect to frequency, the method comprising:

(a) specifying a reference frequency (ω

ref);

(b) computing amplitude and phase spectra of the recorded seismic signals within windows centered at specified times t1 and t2 wherein t2>

t1;

(c) computing normalized amplitude spectra A with respect to frequency (ω

) namely, (A(t1, ω

)) and (A(t2, ω

)), such that;

A(t1,ω

ref)=A(t2,ω

ref);

(d) computing frequency-by-frequency mean amplitudes

Mean A(ω

)=0.5*[A(t1,ω

)+A(t2,ω

)]

(e) identifying frequency pass bands wherein Mean A(ω

) amplitudes are above a specified threshold;

(f) computing phase difference with respect to frequency;

dφ

(ω

)=φ

(t2,ω

)−

φ

(t1,ω

),

(g) computing time delay with respect to frequency;

dt(ω

)=dφ

(ω

)/ω

,

(h) plotting dt(ω

)-vs-ω

scatter on a graph,

(j) fitting a straight line that passes through points on said graph in (h) that lie within said pass bands in (e) such that a Mean A(ω

) weighted average of absolute difference between said points and said line is a minimum,

(k) extracting a residual time delay, δ

t, as the value on said line at a user-specified frequency and dispersion as the slope of said line;

(l) computing a final time delay;

τ

=t2−

t1+δ

t

(m) at least one of storing and displaying the final time delay and dispersion parameter.

5. The method of claim 4 to estimate the Earth'"'"'s attenuation between reflectors by using only the amplitude spectra, the method further comprises:

(n) performing (a) through (e) wherein the specified times (t1, t2) correspond to two different reflectors on a same seismic trace;

(p) computing ln(A(t1, ω

) and ln(A(t2, ω

) for data within the pass bands;

(q) computing an energy absorption time rate with respect to frequency

E(ω

)=2[ln(A(t1,ω

)−

ln(A(t1,ω

)]/(t2−

t1);

(r) fitting a linear least squares line that is constrained to pass through zero at (ω

ref) and also through E(ω

) data points within the pass bands;

(s) storing the slope of the fitted line as an attenuation parameter called ‘

Inverse Q’

between the reflectors at the specified times t1 and t2;

(t) repeating (n) through (s) for a plurality of specified times that belong to a same reflector pair in a plurality of seismic traces, wherein Mean A(ω

) is computed for the plurality of traces.

7. The method of claim 6 further comprising:

(a) computing from the output in step (l) of claim 5 the quantities, kx(ω

)/ω and

ky(ω

)/ω

,

(b) computing the amplitude, A(ω

), weighted mean values of kx(ω

)/w and ky(ω

)/ω

,

(c) the weighted mean of kx(ω

)/ω

gives the apparent time dip, dt/dx, and weighted mean of ky(ω

)/ω

gives the apparent time dip, dt/dy,

(d) the true time dip=sqrt[(dt/dx)̂

2+(dt/dy)̂

2] and azimuth=arc tan[(dt/dy)/(dt/dx) relative to the y-direction, and

(e) storing at least one of the apparent time dip vector and dip and azimuth vector in a medium,

6. A method to estimate wavenumber spectrum continuously at each point in a three dimensional seismic volume of traces associated with an area of Earth'"'"'s subsurface, the method comprising:

(a) computing amplitude and phase spectra of said traces;

(b) specifying a lateral perimeter enclosing a plurality of traces;

(c) placing said perimeter at a center point in said volume;

(d) computing a mean amplitude with respect to frequency A(ω

), over said perimeter at said center time and then computing an average amplitude for all frequencies within the seismic traces;

(e) defining a set of frequency passbands wherein a mean amplitude therein is greater than a specified multiplier of the average amplitude;

(f) initializing a first frequency within each passband;

(g) computing for the first frequency phase differences between the center trace and each other trace in said perimeter at the time of said center point;

(h) fitting a plane such that the plane is constrained to pass through zero at the center trace and minimizes a least squares phase differences with respect to scatter of points in said perimeter;

(j) computing a root mean squared (rms) error (err(ω

)), of phase difference scatter relative to the fitted plane;

(k) comparing err(ω

) with a specified threshold;

if the error exceeds the threshold then reducing the size of the perimeter and repeating (h) and (j) for the reduced perimeter until either the error does not exceed the threshold or a minimum perimeter size is reached, in which case re-fitting the plane is performed as in (h) after removing any outliers

(l) storing a gradient of the final fitted plane that gives horizontal wavenumbers kx(ω

), ky(ω

), and associated err(ω

), the mean amplitude A(ω

) and variance, var[A(ω

)], over the final perimeter,

(m) incrementing to the next frequency in the passband and repeating (g) through (l) until the last frequency in the passband is processed,

(n) incrementing the center time, resetting the perimeter to initial size and repeating the steps (e) through (n) until the end of data is reached,

8. A method to estimate curvature of reflection time surfaces continuously in a three dimensional seismic volume, the method comprises:

(a) obtaining a three dimensional seismic data volume comprising apparent time dips (dt/dx);

(b) specifying a lateral perimeter enclosing a plurality of unaliased seismic data traces;

(c) initializing time at a start of data recording and placing said perimeter at a center point at a beginning of survey area in said volume;

(d) computing (dt/dx) apparent time dip differences between the center trace and each other trace in said perimeter at the time of said center point;

(e) fitting a plane through the scatter of said (dt/dx) differences over said perimeter such that the plane is constrained to pass through zero at the center trace and minimizes a least squares phase differences with respect to the scatter of points in said perimeter;

(f) computing root mean square (rms) error (err) of (dt/dx) difference scatter relative to the fitted plane;

(g) comparing err with a specified threshold;

if the error exceeds the threshold then reducing the size of the perimeter and repeating (c) through (g) for the reduced perimeter until either err does not exceed the threshold or a minimum perimeter size is reached, in which case re-fitting the plane after removing any outlier is performed as in (e);

(h) storing a gradient of the final fitted plane that gives horizontal second derivatives (d²t/dx²), (d²t/dydx), and err, over the final perimeter,

(j) repeating (d) through (h) for all time and space points in the seismic volume to provide two sub volumes, (d²t/dx²) and (d²t/dydx),

(k) replacing the volume of time dips (dt/dx) with the volume of (dt/dy) and repeating (b) through (j) to provide two additional subvolumes, (d²t/dy²) and (d²t/dxdy), and

(l) computing and storing K+ and K−

measures of curvature as defined by equations (2) and (3) in paragraph [0050] and where coefficients (a, b, c) are given by equations (15), (17) and (19) in paragragraphs [0139] and [0140].

10. A method for estimating amplitude and phase spectra of seismic data sampled continuously in time with uniform high frequency-resolution, comprising:

(a) obtaining a representation of the seismic data sampled at an input time rate (Δ

T);

(b) specifying an output time-rate (Δ

T′

) that is an integer multiple of the input time rate Δ

T;

(c) specifying a data analysis window (Δ

t) that is larger than the output time rate (Δ

T′

), wherein the analysis window is an integer multiple of said Δ

T;

(d) specifying a frequency sample rate (Δ

f);

(e) computing a time window (2/Δ

f) that would produce the frequency said sample and is an integer multiple of the input time rate (Δ

T);

(f) generating a model that extrapolates data inside the data analysis window to fill the (2/Δ

f) time window, wherein

(g) the model uses a continuity relationship between data within and outside of the data analysis(Δ

t) window,

(h) the model has the property that forward and backward extrapolated values decrease in amplitude with respect to time from edges of the data analysis (Δ

t) window;

(j) capturing data within the data analysis (Δ

t) window at a start of data time;

(k) extrapolating in forward and backward directions the captured data to fill the (2/Δ

f) time window;

(l) computing fast Fourier transform (FFT) coefficients of data in the (2/Δ

f) time window;

(m) incrementing by the output sample time (Δ

T)′

the start time of data capture for a next FFT computation;

(n) repeating (k) through (m) until and end of data time is reached; and

at least one of

(o) saving in a storage medium the FFT coefficients sampled at (Δ

T′

,Δ

f) intervals, and

(p) saving in the storage medium a transformation of the FFT coefficients into amplitude and phase spectra sampled at (Δ

T′

,Δ

f) intervals.

11. The method of claim 10 wherein:

(a) the mathematical model is a statistical Linear Prediction Filter (LPF) model,

(b) the continuity relationship is an autocorrelation function with a maximum lag, Lmax, of said seismic data in a time window of length, (2/Δ

f), that is symmetrically centered relative to a moving said data time window of length Δ

t, and

(c) the forward extrapolation is done with repeated application of a minimum phase filter derived via Levinson recursion from the autocorrelation function and backward extrapolation is done with repeated application of the time reversed version of the minimum phase filter.

12. The method of claim 11 further comprising:

(a) computing an auto-correlation function of the data in a (2/Δ

f) window to determine statistical correlation between the data within and outside of the signal analysis window,

(b) using the auto-correlation function to compute LP coefficients via Levinson recursion and linearly combining said coefficients with analysis window in a recursive loop for extrapolations in the forward direction and then repeating said linear combination with the reversed sequence of said coefficients for extrapolations in the backward direction,

(c) specifying a plurality of multi-resolution data analysis windows, where each of the multi resolution data analysis windows has a different length, wherein the length thereof is in a range from Δ

t to a fraction of (2/Δ

f),

(d) initializing a start of spectral computation time,

(e) capturing data within the plurality of multi-resolution signal analysis windows, each centered at the spectral computation time,

(f) computing a sequence of LP coefficients for each multi-resolution window,

(g) time extrapolating captured data in each of the multi-resolution windows in forward and backward directions to fill the (2/Δ

f)-window,

(h) averaging data over plurality of multi-resolution extrapolated (2/Δ

f)-windows,

(j) computing FFT coefficients of data in said average (2/Δ

f)-window,

(k) incrementing spectral computation time by said Δ

T′

,

(l) repeating steps (l) through (q) until an end of spectral computation time.

(m) at least one of storing said FFT coefficients as 2-dimensional arrays sampled at (Δ

T′

,Δ

f) intervals and saving in said storage medium a transformation of said FFT coefficients into amplitude and phase spectra as 2-dimensional arrays sampled at (Δ

T′

,Δ

f) intervals, and

(p) displaying the amplitude and phase spectra as 2-dimensional images on an interactive graphical screen.

13. A method to estimate time delay and dispersion between un-aliased time-windowed seismic data, and said time delay is an intercept at a reference frequency and dispersion is the slope of an optimal line through a time delay-with respect to frequency graph, the method comprising:

(a) computing amplitude and phase spectra of said seismic data by a high resolution time-frequency method,

(b) computing frequency-by-frequency phase differences between data windows of said time windowed seismic data,

(c) dividing said phase differences by associated frequencies to compute a time delay-vs-frequency scatter,

(d) computing frequency-by-frequency mean amplitude over said data windows,

(e) identifying frequency pass bands where amplitudes in said data windows are above a specified threshold,

(f) fitting a straight line through points on said graph in (c) that lie within said pass bands in (e) such that a normalized amplitude weighted mean difference between said points and said line is a minimum,

(g) extracting time delay estimate as the value on said line at a user-specified reference frequency and dispersion as the slope of said line.