Method for Migration of Seismic Reflection Waves

This invention relates to a method of determining the presence and location of seismic reflections from the earth's subsurface formations.

In seismic exploration, seismic energy is generated at a shotpoint at or near the surface of the earth, is reflected from subsurface interfaces between layers of the earth, and is received by a spread of geophones on the surface of the earth. The geophone signals are recorded in the form of a seismic record section. This seismic time section contains information which can be used to represent the charateristics of the subsurface formations.

Primary reflection signals on the seismic record section indicate the presence of a subsurface reflecting interface, and time occurances of these primary reflection signals represent the depth of the reflecting interface from the earth's surface. Furthermore, the time-shift or attitude of a primary reflection from trace-to-trace indicates the dip or slope of the subsurface interface.

Skilled seismic interpreters can sometimes determine the location and dip of subsurface interfaces directly from the seismic record section. However, this requires a large amount of judgment due to various distortions introduced into the record section. Such distortions may take the form of multiples, ghosts, reverberations, ground roll, shot noise, etc. Another type of distortion is from the physical phenomena, commonly known as migration. Such migration distortion is introduced by assuming that the reflection events come from a subsurface interface directly below the receiver and have been displayed on the record section as though this were the case. The impedance discontinuity producing the reflection is commonly termed the "reflector surface" and its apparent orientation on a record section is commonly termed the "record surface". The two surfaces as illustrated in FIGS. 1 and 2 will coincide only when the reflector surface is flat and horizontal. Migration is, therefore, a process by which the reflector surface is constructed on the record section as illustrated in FIG. 3. This is carried out by picturing the reflection points A', B' and C' as being moved to an arc to A", B" and C" directly under the shotpoints A, B and C.

In a typical seismic data processing sequence, migration is preceded by conventional common-depth-point stacking. Even though a large fraction of data is processed in this way, there are data sets where the depth-point stack does not enhance all the reflecting events of interest. In areas of complex geometry, criss-cross events may not stack well. Also fault plane reflections are lost during the common-depth-point stack. Diffractions, which help delineate the termination of reflectors, are generaly not preserved during the common-depth-point stack. Further, conventional migration techniques are desinged to operate in a homogeneous medium and are not satisfactory for continuing a seismic wave field in a medium with lateral velocity variations.

It is therefore, an object of the present invention to apply a prestack migration operation in such a manner as to produce a migrated seismic record section which overcomes the difficulties inherent in conventional migration processing.

In accordance with the present invention, there is produced a migrated seismic depth section from seismic data, representing the occurance times of seismic energy from a plurality of spaced sources and reflected from subsurface interfaces. A source of seismic energy transmits seismic energy waves downward into the earth. Seismic energy detectors on the surface of the earth record the upcoming seismic energy waves reflected from subsurface interfaces. Both the known downgoing and the recorded upcoming seismic energy waves are Fourier transformed into the frequency domain. Each frequency component in both the downgoing and upcoming seismic energy waves are downward continued. All the downward continued frequency components of the product of upcoming wave and the complex conjugate of downgoing wave are summed to produce a migrated seismic record section representing the actual depth of the subsurface interfaces.

In a further aspect, the downward continuation of the frequency components of the downgoing and upcoming seismic energy waves employs a finite difference operation applied to the X-Z variables for the frequency components of said seismic energy waves.

In a still further aspect, the summation of the downward continued frequency components is carried out in accordance with arbitrary source-to-reciver imagery.

FIGS. 1-3 illustrate the known principle of migration of seismic reflection signals from subsurface interfaces.

FIG. 4 illustrates a conventional seismic exploration technique.

FIG. 5 is a flow sheet depicting the migration method of the present invention.

FIG. 6 is a water bottom model used to generate a common shot gather of seismic traces.

FIG. 7 is a seismic record section obtained from the water bottom model of FIG. 6 without migration processing enhancement.

FIG. 8 is the seismic record section of FIG. 7 corrected in accordance with the migration enhancement method of FIG. 5.

Referring to FIG. 4, seismic energy, produced of each of a plurality of spaced shotpoints SP.sub.1 -SP.sub.n, is reflected from a subsurface interface at a common-depth-point (CDP) and is received at a plurality of spaced locations X.sub.1 -X.sub.n for the production of a seismic time section. This seismic energy can be described as having two components, a source generated downgoing wave and a field recorded upcoming wave. Even though the source generated downgoing wave is never directly measured, the location and geometry of the source are known. Hence, the surface distribution of the amplitude of the downgoing wave can be estimated. When the velocity of the earth medium is known, both the measured upcoming wave and the estimated source generated downgoing wave can be downward continued to arbitrary depth. Downward continuation is a computation operation by which wave field at any depth level can be computed from a known field specified at a different level, in a manner consistent with the wave equation. The upcoming wave is generated at the moment when the downgoing wave is impinging on a subsurface reflector. Therefore, at the reflector, the upcoming wave and the downgoing wave exist simultaneously and their amplitude ratio is the reflector strength. Conversely, reflectors exist only where there are simultaneous arrivals of downgoing and upcoming waves, the ratio of their amplitudes being the reflector strength.

Referring now to FIG. 5, there is illustrated a flow diagram of the prestack migration operation of the present invention utilizing both the downgoing and upcoming components of the recording seismic signals. The field recorded upcoming wave is initially tranformed at step 10 into the frequency domain. This is carried out by applying the Fourier transform to the differential equation with respect to the time variable t and only the signal frequency band is kept. A finite difference operation is then applied to the X-Z variables for each of the kept frequency components. Such operation includes a split-step time-sharing finite difference algorithm for efficient computations.

More particularly, the downward continuation of upcoming wave is carried out in accordance with the following described expressions. From the upcoming wave measured on the surface U.sub.s (x,z=o, t), where s represents shot location, the upcoming wave over the whole depth is computed by use of the following total velocity gradient upcoming wave equation ##EQU1## where V, the velocity of the earth medium may vary in both x, the lateral distance variable; and z, the depth variable. To accommodate for lateral velocity variations, a split-step algorithm is implemented in the frequency domain wherein equation (1) is split into two equations: ##EQU2## Numerically, equation (2) is implemented with a Crank-Nicolson type finite difference scheme, and equation (3) is directly solved for each depth increment. For a more detailed discussion of splitting-up techniques, reference is made to Marchuk, G. I. (1975), Methods of Numerical Mathematics, Springer-Verlag and to Kjartansson, E. (1978), "Modeling and Migration with the Monochromatic Wave Equation--Variable Velocity and Attenuation", Stanford Exploration Project (Stanford, Calif.: Leland Stanford Junior University) Report No. 15.

In like manner, the source generated downgoing wave is initially transformed at step 12 into the frequency domain. The equation which governs the source generated downgoing wave in step 12 can be obtained from equations (1), (2), and (3) by changing the sign of variable z. ##EQU3## where D(x,z,w) is the fourier transform of the downgoing wave field.

At step 11 each frequency component in the signal band for the upcoming wave is downward continued, while at step 13 each frequency component in the signal band for the downgoing wave is downward continued.

The relevant quantity in the frequency domain is not the downgoing wave itself, but rather the complex conjugate of it. Taking the complex conjugate of equations (5) and (6) from steps 10 and 12 the complex conjugate of the downgoing wave D satisfies the same equation which the upcoming wave U satisfies. ##EQU4## The operation of complex conjugation in the frequency domain corresponds to the time reversal in the time domain, and a time reversed downgoing wave will behave in every respect like an upcoming wave.

After U(x,z,w) is thus computed at all depth levels, the migration process is finished by summing over the real part of the D*(x,z,w)U(x,z,w) for all the frequencies under consideration. This summing is carried out at step 14 in accordance with arbitrary source-to-receiver imagery:

where R(x,z) is the reflector map to be estimated. Both U(x,z,t) and D(x,z,t) are time dependent wave fields whereas R(x,z) is independent time. Thus, equation (9) for any instant of time provides a solution method for R(x,z). However, to take advantage of the available U(x,z,t) and D(x,z,t) for all time, this overly redundant set of equations as described in equation (9) should be solved in a least squares manner:

From which it follows that ##EQU5## In equation (11), the denominator behaves as a scaling factor. In the case of complex geology, it is possible that a certain region in depth is poorly illuminated by the downgoing wave, then equation (11) gives undesirable results. By ignoring the denominator, the following expression is adequate for summing the frequency components to form the arbitrary source-to-receiver imagery of step 14:

In the frequency domain, ##EQU6## Equation (13) shows that only D*, not D is used in imaging. Thus, only D*, not D needs to be computed in the downward continuation process. Equation (13) is implemented with a numerical summing over the discrete frequency components.

To illustrate application of the foregoing described migration operation, reference is now made to FIGS. 6-8. A water bottom model which consists of an undulating water bottom overlying three flat reflectors, shown in FIG. 6, is used to generate a common shot gather of synthetic traces. The shot is located at x=5,000 ft. Ninety-six receiver locations, ranging from x=1,000 ft. to x=10,500 ft. are spaced at 100 ft. Based on ray tracing technique, this common shot gather of traces is generated and is shown in FIG. 7. The prestack migration of the present invention is then applied. The resulting depth section is shown in FIG. 8. All the flat reflectors are accurately mapped.

From the foregoing it will be apparent that the migration method of the present invention can be practiced either with the use of analog computing apparatus such as delay lines, shift registers, drum recorders, contour plotters, etc., or with the use of digital computing apparatus. In one particular embodiment of a digital computing system a Control Data Corporation Model 6600 General Purpose Digital Computer is utilized along with the following input/output components:

Control Computer, 65K Memory

6602 Console Display

6681 Data Channel Converter

405 Card Reader

3447 Card Reader Controller

501 Line Printer

3256 Line Printer Controller

The foregoing has described the migration method of the present invention. It is to be understood that various modifications to the disclosed embodiment, as well as alternative embodiments, may become apparent to one skilled in the art without departing from the scope and spirit of the invention as hereinafter defined by the appended claims.