AU2003278717B2

AU2003278717B2 - A vibrator sweep shaping method

Info

Publication number: AU2003278717B2
Application number: AU2003278717A
Authority: AU
Inventors: Rainer Moerig
Original assignee: Westerngeco Seismic Holdings Ltd
Current assignee: Westerngeco Seismic Holdings Ltd
Priority date: 2002-09-04
Filing date: 2003-08-21
Publication date: 2007-04-26
Anticipated expiration: 2023-08-21
Also published as: NO336868B1; CN1285921C; NO20051672L; AU2003278717A1; CN1678925A

Description

WO 2004/023163 WO 204/03163PCT/US2003!0261 INTERNATIONAL PATENT APPLICATION VIBRATOR SWEEP SHAPING METHOD TITLE: INVENTORS: INVENORS: RAINER MOERIG A VIBRATOR SWEEP SHAPING METHOD BACKGROUND OF THE INVENTION In the acquisition of seismic data, seismic waves are used to determine the nature, orientation and location of subsurface geological formations. In the late 1950s, Conoco, Inc., pioneered the development of a 00 new type of geophysical prospecting technique, now generally known as t "vibroseis" prospecting. Vibroseis prospecting employs a seismic vibrator used to generate a controlled wavetrain that propagates through the earth to be t

C

detected by seismic detectors. The vibrator operator selects an energization (sweep) interval, along with a subsequent period during which the vibrator is not energized. Reflected signals are received both during the time when the vibrator is energized and when it is not. Typically, the energization takes the form of a sinusoidal vibration of continuously varying frequency applied to the surface of the earth (or in a body of water) during a sweep period lasting from about two to about 20 seconds or even more.

Various types of sweeps are possible, each employing some sort of amplitude taper applied at the beginning and at the end of the sweep to insure the amplitude of the sweep goes to zero smoothly at its endpoints. The frequency of the sweep may be varied linearly or nonlinearly with time. The standard signal is a linear sweep. A linear sweep is a sinusoidal-type signal having an essentially constant amplitude envelope where the frequency of the signal varies linearly with time, either increases or decreases monotonically within a given frequency range, and yields a constant sweep rate. A non-linear sweep is a sinusoidal-type signal where the frequency varies nonlinearly with time.

Recently, a new type of signal known as a "shaped-sweep" has been developed for use in vibroseis prospecting. Shaped-sweep technology is disclosed In, US Pat. No. 5,347,494 to Andersen. One benefit of using a shaped-sweep is that the sweep is designed to have an autocorrelation pulse length and an impulse response spectrum that facilitates subsequent data processing activities.

;In vibroseis, data that is generated from the vibratory source is correlated with a reference sweep to produce a. correlated record. A reference sweep signal is, generally, an ideal signal that the vibrator is programmed to apply. The correlated record resembles a conventional seismic record, as one would receive with an explosive or impulsive seismic source.

_It is well known in seismic art that an undesirable byproduct in 17 vibration-generated seismic signals is "side lobe energy." Side lobes are 17 byproducts of the correlation process and lengthen and complicate the desired Scorrelated seismic wavelet. Visually, side lobe energy appears as small oscillations to either side of the central three lobes of the seismic wavelet. Side lobe energy degrades data quality and affects adversely the ability to estimate and control the seismic wavelet in processing. There is, therefore, a need to generate vibrator correlation data that have a simple wavelet shape with minimal side lobe energy. Such data would reduce seismic signal distortion and enhance seismic resolution.

In attempting to solve the side lobe problem, Rietsch, "Vibroseis Signals with Prescribed Power Spectrum," Geophysical Prospecting, Vol. pp. 613-620 (1977), developed a relationship between a sweep's instantaneous phase function and its power spectral density for sweeps having a constant amplitude envelope. This relationship is based on the fact that the power spectrum of a sweep is inversely related to the rate of frequency change of the sweep. Rietsch proposed a method for determining an appropriate phase function for a sweep that has a certain predetermined power spectrum, noting that the method could be used to design sweeps with autocorrelation functions that had low side lobes. Sweeps having predefined power spectra could, theoretically, have been designed using this relationship; however, vibrator electronic control systems of that time could not accurately reproduce or follow a user-defined sweep.

Post Rietsch, the advent of new-generation vibrator control instruments based on advanced microprocessor technology allowed for tight control of the vibrator output force (both amplitude and phase). This advancement enabled user-defined sweeps to be more accurately reproduced and followed by the vibrator. For example, Andersen in US Pat. No 5,347,494, incorporated herein by reference in its entirety, adds to Rietsch's method for producing

(N

;improved wavelet shapes with minimal side energy by using a feedback loop that compensates G for the effects of a non-constant amplitude envelope (taper). Andersen further proposes certain power spectra that produce a substantially three-lobe wavelet. Both Rietsch and Andersen employ methods that relate a sweep's phase function to its power spectrum.

Recently developed cascaded sweep techniques suffer in part from the side lobe problem, described above. In these techniques, multiple sweeps (segments) are linked end-to-end, output by the vibrator, and recorded with a single listen period at the end of the 00 cascaded sweep. This significantly reduces the field time required to record the data as compared to conventional methods. The recorded cascaded sweep data are then correlated with a cascaded reference sweep sequence. The correlated data show, due to the cascading of (N the sweeps, a repetitive (segment) structure. As a result of this structure, side lobes of the first breaks in one segment may extend into a previous segment masking weak reflections and, thus, degrading the data quality.

There is a need in the art for effective and efficient methods to enhance the quality of vibroseis data. Such methods produce wavelets with minimal side lobe energy, solve the side lobe problem in cascaded sweep data, and provide shaped data when cascaded sweeps are used in simultaneous shooting with two or more vibratory sources.

SUMMARY OF THE INVENTION According to the present invention there is provided a method of reducing side lobe energy of a wavelet representing an autocorrelation of a pilot sweep, comprising: shaping an instantaneous frequency of the pilot sweep, wherein the instantaneous frequency is a nonlinear function of time and-autocorrelating the pilot sweep.

BRIEF DESCRIPTION OF THE DRAWINGS So that the manner in which the above recited features of the present invention can be understood in detail, a more particular description of the invention, briefly summarized above, may be had by reference to the embodiments of the invention which are described in the present specification and illustrated in the appended drawings. It is to be noted, however, that the specification and appended drawings illustrate only certain embodiments of this invention and are therefore not to be considered limiting of its scope, for the invention may admit to other equally N\Melbourne\Cases\Patent\55000-55999\P55761.AU\Specis\P55761 AU Specification 2007-3-15.doc 15/03/07

O

ct 0 1^

(N

mq effective embodiments.

FIG. 1 depicts a flow diagram showing one embodiment of the basic method of the present invention.

FIG. 2 depicts a flow diagram for one embodiment of the method of the present invention where the difference between a computed and desired autocorrection wavelet is minimized.

FIG. 3A is a graph of the instantaneous frequency at the beginning of a shaped four-second sweep.

N:\Melbourne\Cases\Patent\55000-55999\P55761.AU\Specis\PS5761.AU Specification 2007-3-15.doc 15/03/07 FIG. 3B Is a graph of the instantaneous frequency at the end of a Sshaped four-second sweep.

Q, FIG. 4 depicts error development during the shaping of a sweep

C

with a 50 Hz Ricker autocorrelation wavelet.

FIG. 5A is a graph of the instantaneous frequency of a linear sweep.

FIG. 5B is an autocorrelation wavelet of a linear sweep.

00 FIG. 5C is a graph of the Instantaneous frequency of a shaped t' sweep.

FIG. 5D is an autocorrelation wavelet of a shaped sweep.

FIG. 6A is the first portion of raw data of a cascaded linear sweep.

FIG. 6B is the first portion of raw data of a cascaded shaped sweep for the same vibrator point as in FIG. 6A.

FIG. 7A is the seismic data of FIG. 6A after correlation using a linear sweep.

FIG. 7B is the seismic data of FIG. 6B after correlation using a shaped sweep.

FIG. 8A is the last portion of correlated sweep data of a linear sweep.

FIG. 8B is the last portion of correlated sweep data of a shaped sweep.

FIG. 9A is raw data of cascaded shaped sweeps acquired with two simultaneous sources.

FIG. 9B is separated and compressed data for the left source of the cascaded sweep data of FIG. 9A.

FIG. 9C Is separated and compressed data for the right source of the cascaded sweep data of FIG. 9A.

DETAILED DESCRIPTION OF THE INVENTION The goal of shaped-sweep technology is to provide methods to shape a vibrator sweep to generate a simple wavelet shape with minimal side lobe energy after correlation. The pilot sweeps used to drive vibrators are, aside from the taper regions, constant amplitude sweeps. The most commonly Sused sweep in data acquisition is a linear sweep with a constant sweep rate S(linear instantaneous frequency). However, non-linear sweeps also are Q commonly employed and are created by changing the sweep rate during sweeping. In so doing, the spectral content and the autocorrelation wavelet can be modified while maintaining a constant amplitude envelope. Non-linear sweeps are usually employed to compensate for the attenuation of the signal Q that occurs as the signal travels through geological formations. Yet regardless of whether the sweeps are linear or non-linear, the autocorrelation wavelets of Cthese sweeps are multiple-lobe wavelets with, generally, only slowly decaying 0side lobes.

Unlike prior art methods, the methods of the current invention do not use the relationship between the pilot sweep's phase function and its power spectrum to minimize side lobe energy. Instead, in the methods of the current invention, the instantaneous frequency of a pilot sweep is shaped by an optimization process in order to obtain an autocorrelation wavelet that is virtually a three-lobe wavelet with minimal side lobe energy.

Figure 1 depicts a flow diagram of the method of the present invention. Generally, in step 10, an instantaneous frequency of a pilot sweep is shaped by application of an algorithm. Next, this pilot sweep is autocorrelated.

In so doing, the autocorrelation wavelet produced has reduced side lobe energy.

Figure 2 depicts one embodiment of the method shown in Figure 1. Referring to Figure 2, a polynomial with a degree n on the order of, for example, about 10 is chosen (step 100). In step 102, the instantaneous frequency f(t) of the pilot sweep is computed. Instantaneous frequency f(t) is a non-linear function of time t and Is represented by an nth degree polynomial n f(t) (1) i=0 where the coefficients ai, the time shift T, and the degree n are parameters that are used to shape the sweep. Initial values for the parameters may be T= a f/=f (3) l^where T is the sweep length and fo and fi are the start and end frequency of a 0linear sweep that otherwise would have been used in the survey area. The value of ao will be approximately the center frequency fc of the correlation wavelet and should be chosen such that it is within the band of the recoverable Sfrequencies in the survey area. The initial value for a, may be a, o (4) which Is, as is well known, the sweep rate of a linear sweep. The initial values for the other coefficients (I are a, 10 for 2< i<n which is small enough to avoid numerical instabilities in the computation of the polynomial.

In addition, it is preferable to obtain an instantaneous frequency function: 1) that is positive for 0 <t and 2) where f increases monotonically with time. Both requirements, however, generally are not fulfilled at all points in time during the sweep. Particularly for polynomials of very high degree, f(t) in equation 1 tends to take on extreme large positive or negative values at both ends of the polynomial's support, within the taper regions.

Thus, in order to achieve the above properties, the instantaneous frequency function must be damped within the taper regions (step 104) of the pilot sweep. This is done by computing a weighted sum of the polynomial in equation 1 and another linear frequency function d(t): f (1 du(t) wu(t) for 0 t< f(t) for t1. t, m with d, t* r tW)=1- O0 alpr ter qp where du and dd are the linear frequency functions in the up- and down-taper regions. wu and Wd are the corresponding weighting functions, fmin and fm~ are user-specified frequency limits, and t,,pe is the taper length assumed in this embodiment to be of equal length in the up- and down-taper regions. Equation provides the damped instantaneous frequency using a linear damping function Instead of using linear functions, other embodiments might employ quadratic functions.

Equation 5 guarantees that the instantaneous frequency function fs(t) starts at fmin and ends at fm,. However, it does not guarantee that is indeed a monotonic function. Therefore, equation 5 is applied repeatedly until monotonic behavior of the instantaneous frequency is obtained (step 105).

Instantaneous frequency was determined for a shaped four-second sweep employing equations 1 and 5 of the present invention. The results obtained within the taper regions of the shaped sweep can be seen in Figure 3. The taper region of the sweep is 0.2 seconds long. The curve labeled "original polynomial" is the instantaneous frequency according to equation 1. The curve labeled "damped polynomial" Is the Instantaneous frequency after repeated applications of equation 5 with fmn 3 Hz and fax 150 Hz.

Once monotonic behavior is attained, the instantaneous phase is Scomputed (step 106). The instantaneous phase 4(t) is the integral of the instantaneous frequency fs(t): (N(6) 0 2a fs(t)dt (6) 00 C- The integration in equation 6 is performed numerically. Next, the shaped sweep is computed (step 108), given by: S(t) sinj(t) (7) where A(t) is the amplitude envelope and P 0 is the initial phase.

The autocorrelation of S(t) is computed (step 110) and compared to the desired autocorrelation wavelet. The error, the sum of the squares of the difference between the computed and desired autocorrelation wavelet, is computed and stored (step 114). Once the error is computed, a determination Is. made as to whether the error falls below a threshold value (or is not decreasing in subsequent iterations) (step 114). If the answer to the determination at step 114 is yes, the error is below the threshold value (or is not decreasing in subsequent iterations), the process is stopped (step 116). If the error is not below the threshold value, an additional iteration of the process is performed.

The process employs an optimization process such as the Levenberg-Marquart method (Press, Teukolsky, Vetterling, Flannery, Numerical Recipes, Cambridge University Press) to continue to minimize the difference between the computed and desired autocorrelation wavelet. Instead of using the Levenberg-Marquart method, one with skill in the art may select any other optimization algorithm suitable to minimize the error. The optimization process involves reiterating steps 100-110. In this Iterative process, the degree of the polynomial of step 100 is increased by, for example, 1 (n 1) (step 118).

Additionally, derivatives of the autocorrelation function at step 110 are determined (step 119) and are used to estimate new parameters (step 120).

Finally, the new polynomial (n 1) (step 118) and the new parameters (step S120) are used to calculate another, further optimized, instantaneous frequency (steps 102-114).

Ci Figure 4 shows the development of error during the sweep shaping process of the present invention with a desired autocorrelation wavelet that is a 50 Hz Ricker wavelet. In this example, the error essentially stops 17- decreasing when the degree of the polynomial approaches 00 t' Example 1 O Figures 5 A-D compare the instantaneous frequency functions and the autocorrelation wavelets of a conventional linear 10-100 Hz sweep (Figures 5A and 5B) and a sweep shaped using the methods of the present invention and having a center frequency of 50 Hz (Figures 5C and 5D). Both sweeps were four seconds in duration and the taper length was 0.2 seconds.

The start and end frequency fmin and fmax used in equation 5 for the shaped sweep were 3 Hz and 150 Hz, respectively. At the end of the up-taper 0.2 seconds) the frequency of the shaped sweep is about 15 Hz and at the beginning of the down-taper (t 3.8 seconds) the frequency of the shaped sweep is about 110 Hz. The side lobe level in the autocorrelation of the shaped sweep is about 60 dB below the main lobe of the three-lobe wavelet.

Example 2 Figures 6A and 6B show the first portion of raw cascaded sweep data using a linear sweep (Figure 6A), and.a sweep shaped by the methods of the present invention (Figure 6B) for the same geological model. The data contain fifth order harmonic distortion, five segments plus a listen period in length, where each segment and listen period length is four seconds long. The initial phases of the segments are 0, 144 288, 72, and 216 degrees.

Figures 7A and 7B display correlated data (of Figures 6A and 6B) after correlation with a linear (Figure 7A) and a shaped (Figure 7B) sweep. The improvement in the quality of data employing a shaped sweep is apparent. The linear sweep data look "noisy" compared to the shaped sweep data. The onset ;of the first breaks in the shaped sweep data is much clearer and better defined Sthan in the linear sweep data, and may serve to improve statics. Furthermore, the high level side lobes in the linear sweep data mask weak reflections that are clearly visible in the shaped sweep data (marked by arrows).

Figures 8A and 8B show the last portion of correlated data of a _linear (Figure 8A) and a shaped (Figure 8B) sweep. The high amplitude noise 17- with hyperbolic move-out seen in the linear sweep data covering the late in- 00 17 coming reflections (seen in Figure 8A) is the side lobe energy of phase-shifted C' first breaks in the subsequent segment. This effect is usually pronounced and often affects one or two seconds of the correlated data. When a shaped sweep is used, however, this effect is much less severe due to the strong suppression of the side lobes (Figure 8B).

Example 3 Figures 9 A-C show cascaded shaped sweep data acquired with two simultaneous sources. When using a single source or vibrator set, a result similar to those shown in Figures 7B and 8B for the shaped sweep can be obtained using a linear sweep for data acquisition and applying a shaping filter in data processing. The shaping filter converts the measured ground force into the desired wavelet. How well this works, however, depends on how well the measured ground force represents the actual ground force put into the ground.

Moreover, such post-acquisition data shaping does not work well on correlated cascaded sweep data, due to the first breaks' side lobes creeping in from the next segment.

Furthermore, post-acquisition data shaping fails completely on raw cascaded sweep data that are acquired with two or more sources simultaneously. The simultaneous sources are to be separated first without compressing the sweeps (maintaining raw data) before a shaping filter could be applied to accomplish the data shaping. The source separation is possible if certain phase rotation schemes are used, but would slow down production significantly. Therefore, the post-acquisition shaping of multiple-source data is considered to be impractical.

On the other hand, pre-acquisition shaping--the deployment of shaped sweeps- -is not affected by the number of simultaneous sources, as can be seen in Figures 9 A-C.

SFigure 9A shows the raw data from cascaded sweeps acquired with two simultaneous sources using shaped sweeps. Figure 9B shows shaped sweep data from the left source that has been S 5 separated and compressed by correlation. Figure 9C shows the corresponding data from the right source.

The equations set forth above are illustrative examples of equations that can be used to shape the instantaneous frequency of a sweep in order to obtain a three-lobe 00 0autocorrelation wavelet with minimal side lobe energy, according to the methods of the present invention. However, these equations should not be considered as limiting the scope of the Sinvention. Those skilled in the art will recognize that the particular equations used to implement Sshaping instantaneous frequency are a matter of convenience and personal preference, and derivatives thereof can be used effectively as well.

In the claims which follow and in the preceding description, except where the context requires otherwise due to express language or necessary implication, the word "comprise" or variations such as "comprises" or "comprising" is used in an inclusive sense, i.e.

to specify the presence of the stated features but not to preclude the presence or addition of further features in various embodiments of the invention.

It is to be understood that, if any prior art publication is referred to herein, such reference does not constitute an admission that the publication forms a part of the common general knowledge in the art, in Australia or any other country.

N:\Melbourne\Cases\Patent\55000-55999\P55761.AU\Specis\P55761.AU Specification 2007-3-15.doc 15/03/07 13

Claims

1. A method of reducing side lobe energy of a wavelet representing an autocorrelation of a pilot sweep, comprising: C 5 shaping an instantaneous frequency of the pilot sweep, wherein the instantaneous frequency is a non-linear function of time and autocorrelating the pilot sweep. 00

2. The method of claim 1, wherein the shaping step comprises computing the instantaneous (N frequency of the pilot sweep. 0 c

3. The method of claim 2, wherein the instantaneous frequency is represented by an nth degree polynomial.

4. The method of claim 3, wherein the nth degree polynomial is greater than about The method of claim 2, wherein the instantaneous frequency of the pilot sweep is computed by n t) aT, (1) .i=0

6. The method of claim 5, wherein initial parameters are (2) f2 (3)

7. The method of claim 2, wherein the shaping step further comprises damping the instantaneous frequency of the pilot sweep within taper regions. N:\Melbourne\Cases\Patent\S5000-55999\P55761.AU\Specis\P55761.AU Specification 2007-3-15.doc 15/03/07 14

8. The method of claim 7, wherein the instantaneous frequency of the pilot sweep within the taper regions are damped by t f (1 w, du(t) wu(t) for 0 5 t t e f(t) for t T- t wd forT t, T Oo with f f f(Tt) and dd(t)= f(T

9. The method of claim 8, comprising the further step of applying equation 5 repeatedly until the instantaneous frequency of the pilot sweep within the taper regions behaves monotonically.

10. The method of claim 2, wherein the shaping step further comprises computing an instantaneous phase of the pilot sweep.

11. The method of claim 10, wherein the instantaneous phase of the pilot sweep is computed by 2 Jf,(t)dt (6) 0

12. The method of claim 2, wherein the shaping step further comprises computing a shaped pilot sweep.

13. The method of claim 12, wherein the shaped pilot sweep is computed by S(t) A(t) sin[(t) (7)

14. The method of claim 7, wherein the shaping step comprises: computing an instantaneous frequency of the pilot sweep; damping an instantaneous frequency of the pilot sweep within the taper regions; computing an instantaneous phase of the pilot sweep; and computing a shaped pilot sweep. N:\Melbourne\Cases\Patent\55000-55999\P55761.AU\Specis\P55761.AU Specification 2007-3-15.doc 15/03/07 The method of claim 14, further comprising after the damping step, continuing damping until the instantaneous frequency within the taper regions behaves monotonically.

16. The method of claim 1 wherein the autocorrelating step further comprises computing an error between the computed and desired autocorrelation wavelet.

17. The method of claim 16, wherein the autocorrelating step further comprises determining whether the computed error falls below a threshold value.

18. The method of claim 17, further comprising increasing the degree of the polynomial if the computed error does not fall below a threshold value.

19. The method of claim 1, further comprising computing derivatives of the autocorrelation function of the pilot sweep. A method as claimed in any one of the preceding claims, and substantially as herein described with reference to the accompanying drawings. N\Melbourne\Case\Patent\55000-55999\P55761.AU\Specis\P55761 .AU Specification 2007-3-l5.doc 15/03/07