Wave-equation trace interpolation
GEOPHYSICS.
VOL. 52. NO. 7 (JULY 1987); P. 973-984, 15 FIGS
Wave-equation
trace
interpolation
Joshua Ronen* These questions, in a general context, have been discussed by several authors (Shannon, 1949; Linden, 1959; Papoulis, 1977; Brown, 1981), and in particular they found that if a function is filtered by some independent filters before sampling, as ~&own in Figures I, a certairr anmurrt of aiiasing is allowed and one can still recover the original signal. The same idea can be applied to seismic data (Bolondi et al., 1982; Claerbout, 1984, p. 219). The earth’s image is the signal that is filtered by some different filters before subsampling. The Hj filters, based on the wave equation, in Figure 2 are the operators that produce seismic data for a given earth model and survey geometry. The data they produce may be aliased in each channel. but the combination of all channels can give a high-resolution image.
ABSTRACT
Spatial aliasing in multichannel seismic data can be overcome by solving an inversion in which the model is the section that would be recorded in a well sampled zero-offset experiment, and the data are seismic data after normal moveout (NMO). The formulation of the (linear) relation between the data and the model is based on the wave equation and on Fourier analysis of aliasing. A processing sequence in which one treats missing data as zero data and performs partial migration before stacking is equivalent to application of the transpose of the operator that actually needs to be inverted. The inverse of that operator cannot be uniquely determined, but it can be estimated using spatial spectral balancing in a conjugate-gradient iterative scheme. The first iteration is conventional processing (including prestack partial migration). As shown in a field data example in which severe spatial aliasing was simulated, a few more iterations are necessary to achieve significantly better results.
FORMULATION
Multichannel abased one-dimensionalsignal When a model function #z(x) is sampled with a certain sampling interval AX, its Fourier transform is replicated: if @I(X) has the Fourier transform m(k), then the Fourier transform of the sampled signal is FT{sampled[ril(ll]~
= njcm(k
- nk,J,
(1)
INTRODUCTION
where the sampling wavenumber k, is The Nyquist (1928) condition for sampling a wave field is that the sampling interval may not exceed half the smallest wavelength. While this is achievable in two-dimensional (2-D) surveys, practical difficulties arise in three-dimensional (3-D) surveys. The result is a problem of missing data; one does not have enough data to sample the seismic wave field adequately. On the other hand, it seems that one has too much data in multichannel seismic surveys; there are many traces in every common-midpoint (CMP) gather, and summing (stacking) is used to merge what seems to be redundant information from different channels. The questions are whether the information is redundant and whether summing is the optimal way to merge it.
k, = 2nlAx. Nyquist (1928) showed that if the function m is band-limited, m(k) = 0
for Ikl 2 W,
then the summation in equation (1) has a finite number of terms N, where N is the smallest integer that is larger than 2 W/k,. In particular, if k, is more than 2 W, then N = 1 and FT (sampled [K(x)]}
= m(k),
for (k ( < W. In this case, m(k) can be easily restored from its samples simply by extracting the spectrum for 1k 1< kJ2.
Manuscriptreceivedby the Editor March 7, 1986;revisedmanuscriptreceivedSeptember29, 1986. *Formerly Departmentof Geophyscs,StanfordUniversity;presentlyCenterfor Wave Phenomena,ColoradoSchoolof Mines,Golden,CO 80401. al 1987Societyof Exploration Geophysicists. All rights reserved. 973
Ronen
974
Recovering from aliasing
Analogfilters
Sub-
sampling
Aliasedsequences
FIG. 1. Overcomimg aliasing of a one-dimensional signal. The signal at the left is filtered by five filters H,, .., If,, before sampling by the analog-to-digital (A/D) converters. The data are five different aliased sequences.The original signal cannot be recovered from any one sequencealone, but the combination of them may be suffkient.
Multioffset data
Wave propagjation
Constantoffset sections [aliased 1
FIG. 2. Overcoming aliasing in reflection seismology. The earth reflectivity model is the signal sampled after it has
propagated through the filters Hj.
975
Wave-equaiion Trace interpolation The point is that aliasing can be overcome if the function 6i(.x) is filtered by some convolution filters hj before sampling. The data in channel j are samples from the function h,(x) * G(x). In the wavenumber domain, one has d,(k) = C h,(k - nk,)m(k -
nk,s),
forj=1,2
k) a matrix-vector multiplication, d,(k) = Q,? (k)m(k).
(3)
The vectors are the common-offset section
,._., J. dj(k) =
In a matrix form, for five channels (J = 5) and three-fold aliasing (N = 3), (nt is the number of samples in a trace) and the zero-offset section
m(k) =
The solution of this system of equations interpolates on the x-axis because it extrapolates on the k-axis. The data are known up to the wavenumber kJ2, the model is found (in this case) for ) k ) c 3kJ2, because m(k), m(k - k,), and m(k + k,J are found for )k ) < kJ2.
(no is the number of frequencies in the Fourier transforin of a trace). The matrix Q,’ is the inverse DMO operator (an nt by nw matrix). The tth row and the wth column element are A-
Multichannel abased seismic data To apply the idea of equation (2) to reflection seismology, one must first describe the relations between the earth (the model) and the seismicdata. For the purpose of trace interpolation, I chose the model to be the section that would be recorded if a well-sampled, zerooffset survey were performed. The relationship between the zero-offset section and the earth’s reflectivity is approximately linear; the linear operator is poststack migration. Non-zero offset data are also approximately linearly related to the reflectivity; the linear operator is prestack migration. Both the zero-offset section and the non-zero offset section are linear with reflectivity; therefore, they are linearly related to each other; the linear operator is a prestack-partial migration also known as dip moveout (DMO). The prestack migration process can be performed in three steps: normal moveout (NMO), DMO, and poststack migration (Figure 3). Choosing the model as the unmigrated zerooffset section and taking the data after NM0 leave only the DMO and the spatial sampling between the data and the model (Figure 4). To set up the relation between the model and the data, the inverse DMO, and not the DMO itself, is needed. It is shown in the Appendix that the relation between the time-space Fourier transform of the zero-offset section m(k, w) and the space Fourier transform of a common-offset section dj(t, k) is
One can now merge sampling theory [equation (i)] and migration [equation (3)]. When spatial aliasing is introduced, equation (3) becomes d,(k) = c Qf (k - nk,)m(k - nk,). n
Ptestack
migration
ie
dw A
where rl = ,/I + [hj k/w]’ form of this integral is
~ ‘ -‘oartn(~,
k),
and llj is the half-offset. A discrete
dj(t, k) = 1 A-‘e
~i2noAr~nrtn(w, k),
w (nt is the number of time samples in a trace). This is (for each
(4)
Prestack partial migration
mid oint
Constant offset 22 section $J
.
NM0 $: Bi DMO r-i
Model d, (t, k) =
1 e - ilnodrlnt
m
1
f
,$ 5
FIG. 3. Decomposition of prestack full migration into NMO, DMO, and poststack migration. An impulse on a commonoffset section can be migrated to the expected ellipse with the shot and the receiver in the foci, either in one process (ptestack full migration), or in three steps (NMO, DMO, and zero-offset migration. which can be done as a poststack process).
976
Ronen
For J = 5 offsets and an aliasing of N = 3, the system is
~~~=~~~i._~~~]~~~~~.
(5) This is the same as equation (2), except d, and m are vectors and &),Tis a matrix.
Equation (5) is a system of J x nt equations (one equation for every time sample of every common-offset section), with N x no unknowns for each wavenumber k. The matrix has J x N blocks with each block Q,? (k - nk,) being an nt x nw matrix. Solving this equation extrapolates the data in the offset direction (to zero offset), and also in the wavenumber direction from the low wavenumbers (k 1< k,/2, on which the aliased data dj are given, to the full range (k 1 < Nk,/2 required to describe the unaliased zero-offset section. This formulation is similar to the multichannel deconvolution given by Davies and Mercado (1968). Both formulations involve inversion of large block matrices, except that here the blocks are inverse-DMO operators instead of convolutions.
Wave-equation trace interpolation in two steps Conventional data processing does not involve inversion of huge matrices, such as the one in equation (5), and missing data are often assumedto be zero. Nevertheless,it works quite well when the data are not severely aliased. If there is no spatial aliasing, then N = 1 in equation (5), and one has
(6)
m-
4
D;
H A/D
k
dj If DMO is applied to each common-offset section, stack = i Qjdj j=l
-
.
. = CQjLlfm ,:
1
.I x
=
1
DMO- I*
sub-
DMO
sampling
--
Commonoffset section
-__I b r2; Ii
Prestack migration
Jm.
(7)
The processing in equation (7) is the application of the Hermitian transpose of the matrix in equation (6),
D,(k)
Q,(k)
0 E
Common offset m section I
Cm j=1
Dh(k)
Q,(k)
1
(the asterisk denotes taking the Hermitian transpose). This application usesthe result p+ *,(,4) = Q, {k) I
FIG. 4. The choice of the model and the data. The model m is the ideal zero-offset section. The data d. are the commonoffset sections after NMO. The Q,’ are &verse-DMO operators.
(8)
given in the Appendix. When the data are not spatially aliased, equation (6) can be solved by applying the Hermitian transpose matrix. Since this is adequate processing for the nonaliased case, it is interesting to see what the results of applying the Hermitian transpose are in general.
Wave-equation
Trace
Interpolation
977
400
600
800
N
FIG. 5. Synthetic example. (a) The model: a well-sampled, zero-offset section. (b) The data: four midpoint gathers after NMO. The midpoint interval is 200 m. Each gather has five traces (at the same midpoint but with different offsets). Note that the dipping reflector is overcorrected by the NMO. (c) Transpose processing: the result of DMO and stacking after addition of zero traces. (d) The result of wave-equation trace interpolation.
I
I
ORIGINRL I
OQTFI I
I
I I
r
SRtiPLED 1
DRTG I
I
I
1
I
1
I
I
I I
1
1
4
shot (4
I
9
5oom
shot
I
E
8
I
@I
6. (a) Stacking chart of the original data. The original data do not have an aliasing problem: 25 m shot interval, 12.5 m group interval (one offset is missing because of a dead receiver). (b) Stacking chart of the data used for the trace interpolation. Less than 1 percent of the data were used, simulating severe spatial aliasing. FIG.
978
Ronen
If equation (5) is written as
Synthetic data example
d=Gm,
(9)
(where G is the J x nt by N x nw matrix), then the vector ti can be defined as m = G*d.
(10)
To find 61,consider the example (N = 3 and J = 4). Then
_+E;/;
;;
;;;;I.
The transpose of G is D,(k - k,) D, (k - k,) !I, (k - k,) D,(k - k,) G* =
Dz (k) !I, (4 D,(k) [@(k + A,, ” &Z(k+!e) 5 &(k+k,)
124(k) ?,(k+kj
The two-step method was applied to a synthetic data set of four CMP gathers, 200 m apart. Each gather had five traces with 0 to 800 m offset (Figure 5b). The earth model includes a 90 degree dipping reflector and a flat reflector. The ideal solution, the well sampled (25 m midpoint interval) zero-offset section, is shown in Figure 5a. This section is not part of the simulated data but is included for comparison with the results of the trace interpolation. The results of applying the transpose operator (DMO and stacking) are shown in Figure SC. Poststack correction, using QR factorization to invert G*G + C, (with a diagonal c), generated the section shown in Figure 5d. The improvement by the poststack correction is substantial, in particular for the Rat reflector. The conjugate gradient method with spectral balancing
1 I “‘J
The two-step method is not easily applicable to field data because the formation and the inversion of G*G are too com-
where equation (8) was used again. KI is therefore
d,(k)
D,(k
- ks)
P,(k)
Dz (k + k,)
Dx (k + k,)
D,(k
for - 3k,/2 < k =c - kJ2
i Dj(k)d,(k) j= 1
for - kJ2 < k < kJ2
i Qj(k)dj(k-ks) j= I
for k,/2 < k < 3k,/2.
(11)
! Equation (1 I) describes DMO (application of Dj) and stacking
using the (unjustified) assumption d(k -
k,) = d(k) = d(k + k,).
P,(k
Da(k)
j=l k(k)=
- k,J
9, (k)
for ) k ) < kJ2. This can be written as
i D,(k)dj(k+kJ
b(k
(12)
This assumption of replicating the Fourier transform is equivalent to adding zero traces between data traces in the space domain. An important conclusion is that G* is equivalent to DMO and stacking, treating missing data as zero data. Applying C* [equation (lo)] can be the first step in solving equation (9). The second step is
[Actually G*G is singular, so the pseudoinverse (Strang, 1980) should be used.] To summarize, the two steps are (1) DMO and stacking with zero traces in place of missing traces [equation (lo)], and (2) a poststack process [equation (13)].
-
IiI
k,)
d, (k) d,(k)
’
+ kJ
4 (k)
putationally intensive. Fortunately, there are practical iterative techniques, such as the conjugate gradient method, that produce good results at a reasonable cost. Conjugate-gradient inversion is an iterative method in which each iteration involves the application of a forward operator and its transpose (Golub and Van-Loan, 1983; Luenberger, 1984; Paige and Saunders, 1982). Conjugate-gradient inversion of an M x M matrix converges within M iterations. Often, just a few iterations (much less than M) are actually needed. To solve equation (9) each conjugate-gradient iteration involves DMO stacking (G*) and transpose DMO stacking (G). Therefore, the total cost in four iterations is the cost of eight DMOs. Actually, the cost is less than that because Fourier transforms are not repeated, data are not resorted, and operators are re-used. When applied to field data, straightforward conjugategradient inversion converged reasonably quickly, but to a disappointing result. The reason for the poor results is that solving equation (9) involves inversion of singular matrices that have no exact inverse. Moreover, even the approximate solution is not unique. There are many models that reasonably fit the data. Those models differ by components of what is called the null space. Null-space components are the eigenvectors that correspond to zero singular values. [See Strang (1980) for a diszussiorr of singular-vaiue~decomposition and null space.] To achieve a unique solution, the null-space components in the model must be estimated. A standard estimation procedure for the null space in leastsquares problems is damping: no null-space components in the solution. Damping was used in the synthetic example to gcneratc Figue 5d, which is of a reasonable quality. For the
Wave-equation
150
300
Interpolation
979
field data, one has to be more careful, because in practice the null space is composed of the eigenvectors that correspond to small (not necessarilyzero) singular values. How small is small is determined by the noise level. The synthetic example is free of noise, so the null space is small and damping produces a reasonable model. Field data, however, have enough noise to require a more sophisticated estimation of the null space. The estimation of the null-space components is an optimization problem, the objective of which is to fit a priori assumptions. A safe a priori assumption is that the spatial spectrum of the earth model should be balanced. A study of spatial abasing in multichannel seismic data (Rocca and Ronen, 1984) predicted that in the presence of spatial aliasing, when performing DMO and stacking there will be resonance wavenumbers for which the spectrum will have high amplitude. An a priori assumption is that these are artifacts. This a priori assumption can be used in an iterative optimization scheme (such as the conjugate gradient method). Each iteration can involve, in addition to the usual conjugategradient routine, the following spectral balancing procedure:
(m) 0
Trace
450
(1) Transforming the raw model into the wavenumber-time (k, r) domain. (2) Calculating the input envelope (a leaky integration of the absolute values will do). The expected envelope is calculated based on the input envelope monotonously decreasing with wavenumber. (3) Dividing the input by its envelope, and multiplying by the expected envelope. Field-data
Fto. i. Tne data used for the trace interpolation. Every trace is corrected for gain and NMO, and is plotted at its midpoint position. As seen on the stacking chart of Figure 6b, there are four groups of five near traces and four groups of far traces. This figure illustrates what would be obtained by NM0 and stacking.
RAW 200
1
RAW 400
200
2
RAW 400
200
3
A small part of a marine line from the Gulf of Mexico was used as a test case. The original data do not have any aliasing because a shot interval of 25 m and 240 channels, 12.5 m apart, provides a 6.25 m midpoint interval on the 6000 percent stack. Only a small number of these data were used: 150 m
RAW 400
example
200
4
RAW 400
200
5
RAW 400
200
0 400
FIG. 8. The models of the first six iterations. Convergence is achieved within four iterations, but the resulting model is disappointing; it is very different from the true zero-offset section of Figure 12.
(m)
980
Ronen
RAW
1
0
RAW
0.05 0
2
RAW
0.05
3
RAW
0.05
4
RAW
0.05
5
RAW
0.05
6
0.05
(cycles/m)
N
FIG. 9. The raw (unbalanced) spectra of the first six iterations. As predicted (Rocca and Ronen, 1984), there are some resonance wavenumbers where the amplitude is too high compared to the true spectrum of Figure 13. BALANCED 1 BALANCED 2 200 400 200 400
BALANCED 3 200 400
BALANCED 4 200 400
BALANCED 200
5
BALANCED
400
200
6 400
0
0-n)
N
FIG. 10. The models of the first six iterations with spectral balance. Compared to Figure 8, these sections are much more similar to the near-offset section in Figure 12. BALANCED
0
0.05
1 BALANCED
2 BALANCED
0.05
3 BALANCED
0.05
4 BALANCED
0.05
5 BALANCED
0.05
6
0.05
(cycles/m)
0
FIG. 11. Spectra of the first six iterations with balancing. The spectral balancing routine used the knowledge that the raw spectra of Figure 9 may have resonance wavenumbers in which the amplitude is too high.
Wave-equation
0-W 200
Trace
981
Interpolation
shot interval (killing five shots out of every six), and ten channels per shot (killing 230 channels out of 240) in two groups, using five near traces and five far traces. This particular selection of data simulates the aliasing that would occur in the cross-line direction of data collected in a two-boat survey with a 150 m line interval where the boats alternate shooting and only one boat is recording (Ronen? 1985). A stacking chart of the data is shown in Figure 6. The data that were used in the inversion are shown in Figure 7. As in many 3-D surveys, the coverage is not uniform. Nonuniform coverage is natural for the data; it is the model, not the data, which should be uniformly and adequately sampled. The models obtained by the first six iterations (Figure 8) are the disappointing results mentioned earlier in this paper. The spatial spectra are clearly unbalanced (Figure 9). The models for the first six iterations, with spectral balancing, are shown in Figure 10, and their spatial spectra are shown in Figure 11. The results with spectral balancing are a much better estimate for the near-offset section (Figure 12) and its spectrum (Figure 13) compared to the results without spectral balancing (Figures S-9). The final result (Figure 15), although not precisely the near-
0-N 200
400
0
FIG 12.The near-offset section. This section at offsets of about 220 m is considered the “true solution.” Only four traces from this section were included in the data of Figure 7.
TRUE
SPECTRUM
0
FIG. 13. The spectrum of the near-offset section of Figure 12. This spectrum was not used in the trace interpolation, but is shown here for comparison with the estimated spectra of Figures 9 and 11.
FIG. 14. The model obtained by the first iteration, without spectral balancing. This would be the result of adding zero traces in place of missing data, then performing DMO and stacking.
R0n0ll
982
offset section (Figure 12), is significantly better than what is achieved by DMO stacking (Figure 14), and tremendously better than the result of NM0 and stacking (Figure 7). CONCLUSIONS
The trace-interpolation method is based on the wave equation and on an a priori assumption of a smooth spatial spectrum. The method involves an iterative inversion of the linear relations d = Gm, where the model vector m is the ideal zerooffset section, and the data vector d is the seismic data after NMO. Treating missing data as zero data and performing DMO and stacking are equivalent to application of the transpose operator I?I = G*d. This processing is adequate only in the absenceof spatial aliasing. 0-n) 200
400
Two alternative methods of trace interpolation presented are a poststack process (on synthetic data) and an iterative, conjugate-gradient scheme (on field data). The poststack process involves large matrix multiplication and inversion and is too costly to be applied to field data. The iterative method, a combination of conjugate-gradient and spectral balancing, was successfully applied to field data. In both methods, the first step, or first iteration, is DMO and stacking. ACKNOWLEDGMENTS
This work would not have been possible without the guidance and support of the staff and sponsors of the Stanford Exploration Project. Special thanks are due to Fabio Rocca who started me on this project; to Jon Claerbout who taught me how to make it work; to Einar Kjartansson, Stew Levin, Rick Ottolini, Chuck Sword, and John Toldi with whom I often discussedraw ideas, and whose software was very useful; to Paige and Saunders for writing LSQR; and to Chevron who donated the field data. Thanks are also due to Jim Black and Dave Hale who carefully reviewed the manuscript.
REFERENCES
Bolondi,G.. Loinger,E., and Rocca,F., 1982,Offset continuationof seismicsections:Geophys.Prosp.,30, 813-828.
FIG. 15. The model obtained at the sixth conjugate-gradient iteration with spectral balancing. This is the estimation for the zero-offset section (Figure 12) that the wave-equation trace interpolation produces. Flat as well as dipping reflectors are well-interpolated.
Brown, J. L., 1981, Multichannel sampling of lowpass signals: Inst. Electr. Electron. Eng., Trans. Circuit Theory, CAS 28, 101-106. Claerbout, J. F., 1984. Imaging the Earth’s interior: Blackwell Scientific Publ. Ltd. Davies E. B., and Mercado, E. J., 1968, Multichannel deconvolution filtering of field recorded seismicdata: Geophysics, 33,711-722. Dereeowski. S. M.. and Rocca, F.. 1981, Geometrical optics and wave theory of constant offset sections in layered media: Geophys. Prosp., 29, 384-406. Golub, G. H., and Van-Loan, C. F., 1983, Matrix computations: The Johns Hopkins Univ. Press. Hale, D., 1983, Dip moveout by Fourier transform: Ph.D. thesis, Stanford Univ. Judson, D. R., Schultz, P. S.. and Sherwood, J. W. C., 1978, Equalizing the stacking velocities of dipping events via Devilish: Presented at the 48th Ann. Internat. Mtg., Sot. Explor. Geophys., San Francisco. Linden, D. A., 1959, A discussion of sampling theorems: Proc., I.R.E., 47, 1219. Luenberger, D. G., 1984, Linear and non linear programming: Addison-Wes!w. Nyquist, H., 192& Certain topics in telegraph transmission theory: AIEE Trans., 47,617-644. Ottolini, R. A., 1982, Migration of reflection seismic data in anglemidpoint coordinates: Ph.D. thesis,StanfordUniv. Paige, C. C., and Saunders,M. A., 1982,LSQR: An algorithm for sparse linear equations and sparse least squares: Assn. Comp. Math., Trans. Math. Soft., 8, 195-209. Papoulis, A., 1977, Generalized sampling expansion: Inst. Electr. Electron. Eng., Trans. Circuit Syst., CAS 24, 652-654. Rocca, F.. and Ronen, J., 1984, Improving resolution by dip moveout: Presented at the 54th Ann. Internat. Mtg., Sot. Explor. Geophys., Atlanta. Ronen, J., 1985, Multichannel inversion: Ph.D. thesis, Stanford Univ. Shannon, C. E., 1949, Communication in the presence of noise: Proc. I.R.E., 37, l&21. Strang, G., 1980, Linear algebra and its applications: Academic Press, Inc. Yilmaz, O., 1979, Prestack partial migration: Ph.D. thesis, Stanford Univ.
Wave-equation Trace Interpolation
983
APPENDIX INVERSE DMO
Different common-offset sections are not independent because they are all reflected from the same reflectors. Many authors (Judson et al., 1978; Yilmaz, 1979; Deregowski and Rocca, 1981; Bolondi et al., 1982; Ottolini, 1982; Hale, 1983) have shown how to get the zero-offset section from commonoffset sections. For the trace interpolation, the inverse operator is needed and it is derived in this Appendix.
Zero-offset migration is described by the dispersion relation
(A-5) (k, , li,, kZ, and WI,are the wavenumbers that correspond to X, I’, Z, and t,,) Inserting the Fourier transform of the change of variables equation (A-4)
Prestack full and partial migration The presence of a single spike in the data set, at midpoint (x, ,3 = (0, 0), and traveltime t,, with a shot at (x, y, z) = (h, 0. 0) and a receiver at (-h, 0, 0), means there is an ellipsoidal reflector such that the distance from any point (x. I’. z) on the ellipsoid to the source plus the distance from that point to the receiver equals the trayeltirne multiplied by the velocity, ~~(x-~,2+~~+Z*+~~~++)~+~~+z~=uth.
(A-l)
Points (x, y, Z) in equation (A-l) are on an ellipsoid which can be written as
J
into equation (A-5) gives
which is the dispersion relation for prestack, post-NMO, full migration. Both t, and w, appear in equation (A-6) because the operator is time-variable. Prestack full migration can be done in three steps (Figure 3): (1) NMO: mapping,
Perform the well-known trace-by-trace
a, can be found by setting 2 = 0 and r‘ = 0 in the last two equations, i.e.,
t; = th’- (2h/a)Z. (2) DMO:
(A-2)
Define (A-7)
Similarly, a, = a, = V/(ut,/2)2 - h2 = ut,/2,
co0is the frequency that corresponds to the time on the zero-offset section. (3) Migration: Using
where t, is the normal moveout (NMO) time r; = t,: - (2h/U)Z.
(A-2)
[gyk:+k;+k:]=lu:.
Therefore, prestack migration of the impulse produces the ellipsoid
(A-8)
obtain the migrated section. Inverse DMO or
(A-3)
x2 1 + (2/l,/nr,)Z
t 4.2f 22 = L
vt,/2 a, 1
d(t,, k,, h) =
The change of variables
1 done-‘“+@,,
k,, h)
J
and using equation (A-7) to change variables o,, to o,, ,
x2
x2 =
1+
(2h,ILQ
(A-4)
compressesthe ellipsoid (A-3) to the sphere X2 + $ + z2 = (ot,/2)2, which is the zero-offset migration of a spike at the NM0 t”
Equation (A-7) can be used to derive the relation between the zero-offset section m and the common-offset section d. Starting from the inverse Fourier transform
time
984
Ronen
Defining
with A = ,/l
+ (hkJw, t,)‘,
A= Jl
one obtains (A-10)
W” = Ao,
To return to the original (x, y) coordinates, recall that a rotation in the space domain corresponds to a rotation in the same direction in the wavenumber domain:
and
4l a0 = A-‘, do, - a,
cos 9 kX’ z-z kY’ -sin0 (1[
(A-11)
sin 0 cos8
k,
I I k,
Therefore,
Substituting equations (A-10) and (A-11) into equation (A-9),
d(t,, k,, h) =
(A-14)
+ (hk,./o, t,)‘.
dw, A-le-‘“OA’nm(oo,
k,),
hk,, = h(cos Ok, +
(A-12)
sin
Ok,)
= (h cos O)k, + (h sin O)k,
s = h,k, + h,ky
which implies a method of generating common-offset sections d(t,, x, h) from the zero-offset m(t, , x).
=h.k. Equations (A-13) and (A-14) become
Inverse DMO in three dimensions In equation (A-l), I made the assumption that the shot and receiver are on the x-axis. In general they are not, and the half-offset vector has two components:
d(t,, k,, k,, I?)=
do, A-le-‘“UA”m(oo,
k,, k,,)
s and A = \/l
+ (h - k/w,t,)‘,
h=
where h. k is the linear product of the offset and the wavenumber vectors.
Let 0 be the angle of rotation such that
Summary
One can rotate the (x, y) axis to (x’, y’) by
[I [ XI
The inverse DMO operator is given by
=
-sin0
Y’
d(r, k, h) =
cost3
so that the shot and receiver are on the x’-axis, and one can make the same development of equations (A-l) through (A- 12) in the (x’, y’) coordinates, obtaining
dw A- e‘ -iWA’m(w, k).
This is the Hermitian transpose to the DMO operator, m(w, k) =
dt Am’e-‘OA’d(t, s
d(t,, k,. , k,, ,
11)=
dw, As
‘e-iooALm(wO
k,,
, k,.),
(A-13)
that was derived by Hale (1983).
k, h),
VOL. 52. NO. 7 (JULY 1987); P. 973-984, 15 FIGS
Wave-equation
trace
interpolation
Joshua Ronen* These questions, in a general context, have been discussed by several authors (Shannon, 1949; Linden, 1959; Papoulis, 1977; Brown, 1981), and in particular they found that if a function is filtered by some independent filters before sampling, as ~&own in Figures I, a certairr anmurrt of aiiasing is allowed and one can still recover the original signal. The same idea can be applied to seismic data (Bolondi et al., 1982; Claerbout, 1984, p. 219). The earth’s image is the signal that is filtered by some different filters before subsampling. The Hj filters, based on the wave equation, in Figure 2 are the operators that produce seismic data for a given earth model and survey geometry. The data they produce may be aliased in each channel. but the combination of all channels can give a high-resolution image.
ABSTRACT
Spatial aliasing in multichannel seismic data can be overcome by solving an inversion in which the model is the section that would be recorded in a well sampled zero-offset experiment, and the data are seismic data after normal moveout (NMO). The formulation of the (linear) relation between the data and the model is based on the wave equation and on Fourier analysis of aliasing. A processing sequence in which one treats missing data as zero data and performs partial migration before stacking is equivalent to application of the transpose of the operator that actually needs to be inverted. The inverse of that operator cannot be uniquely determined, but it can be estimated using spatial spectral balancing in a conjugate-gradient iterative scheme. The first iteration is conventional processing (including prestack partial migration). As shown in a field data example in which severe spatial aliasing was simulated, a few more iterations are necessary to achieve significantly better results.
FORMULATION
Multichannel abased one-dimensionalsignal When a model function #z(x) is sampled with a certain sampling interval AX, its Fourier transform is replicated: if @I(X) has the Fourier transform m(k), then the Fourier transform of the sampled signal is FT{sampled[ril(ll]~
= njcm(k
- nk,J,
(1)
INTRODUCTION
where the sampling wavenumber k, is The Nyquist (1928) condition for sampling a wave field is that the sampling interval may not exceed half the smallest wavelength. While this is achievable in two-dimensional (2-D) surveys, practical difficulties arise in three-dimensional (3-D) surveys. The result is a problem of missing data; one does not have enough data to sample the seismic wave field adequately. On the other hand, it seems that one has too much data in multichannel seismic surveys; there are many traces in every common-midpoint (CMP) gather, and summing (stacking) is used to merge what seems to be redundant information from different channels. The questions are whether the information is redundant and whether summing is the optimal way to merge it.
k, = 2nlAx. Nyquist (1928) showed that if the function m is band-limited, m(k) = 0
for Ikl 2 W,
then the summation in equation (1) has a finite number of terms N, where N is the smallest integer that is larger than 2 W/k,. In particular, if k, is more than 2 W, then N = 1 and FT (sampled [K(x)]}
= m(k),
for (k ( < W. In this case, m(k) can be easily restored from its samples simply by extracting the spectrum for 1k 1< kJ2.
Manuscriptreceivedby the Editor March 7, 1986;revisedmanuscriptreceivedSeptember29, 1986. *Formerly Departmentof Geophyscs,StanfordUniversity;presentlyCenterfor Wave Phenomena,ColoradoSchoolof Mines,Golden,CO 80401. al 1987Societyof Exploration Geophysicists. All rights reserved. 973
Ronen
974
Recovering from aliasing
Analogfilters
Sub-
sampling
Aliasedsequences
FIG. 1. Overcomimg aliasing of a one-dimensional signal. The signal at the left is filtered by five filters H,, .., If,, before sampling by the analog-to-digital (A/D) converters. The data are five different aliased sequences.The original signal cannot be recovered from any one sequencealone, but the combination of them may be suffkient.
Multioffset data
Wave propagjation
Constantoffset sections [aliased 1
FIG. 2. Overcoming aliasing in reflection seismology. The earth reflectivity model is the signal sampled after it has
propagated through the filters Hj.
975
Wave-equaiion Trace interpolation The point is that aliasing can be overcome if the function 6i(.x) is filtered by some convolution filters hj before sampling. The data in channel j are samples from the function h,(x) * G(x). In the wavenumber domain, one has d,(k) = C h,(k - nk,)m(k -
nk,s),
forj=1,2
k) a matrix-vector multiplication, d,(k) = Q,? (k)m(k).
(3)
The vectors are the common-offset section
,._., J. dj(k) =
In a matrix form, for five channels (J = 5) and three-fold aliasing (N = 3), (nt is the number of samples in a trace) and the zero-offset section
m(k) =
The solution of this system of equations interpolates on the x-axis because it extrapolates on the k-axis. The data are known up to the wavenumber kJ2, the model is found (in this case) for ) k ) c 3kJ2, because m(k), m(k - k,), and m(k + k,J are found for )k ) < kJ2.
(no is the number of frequencies in the Fourier transforin of a trace). The matrix Q,’ is the inverse DMO operator (an nt by nw matrix). The tth row and the wth column element are A-
Multichannel abased seismic data To apply the idea of equation (2) to reflection seismology, one must first describe the relations between the earth (the model) and the seismicdata. For the purpose of trace interpolation, I chose the model to be the section that would be recorded if a well-sampled, zerooffset survey were performed. The relationship between the zero-offset section and the earth’s reflectivity is approximately linear; the linear operator is poststack migration. Non-zero offset data are also approximately linearly related to the reflectivity; the linear operator is prestack migration. Both the zero-offset section and the non-zero offset section are linear with reflectivity; therefore, they are linearly related to each other; the linear operator is a prestack-partial migration also known as dip moveout (DMO). The prestack migration process can be performed in three steps: normal moveout (NMO), DMO, and poststack migration (Figure 3). Choosing the model as the unmigrated zerooffset section and taking the data after NM0 leave only the DMO and the spatial sampling between the data and the model (Figure 4). To set up the relation between the model and the data, the inverse DMO, and not the DMO itself, is needed. It is shown in the Appendix that the relation between the time-space Fourier transform of the zero-offset section m(k, w) and the space Fourier transform of a common-offset section dj(t, k) is
One can now merge sampling theory [equation (i)] and migration [equation (3)]. When spatial aliasing is introduced, equation (3) becomes d,(k) = c Qf (k - nk,)m(k - nk,). n
Ptestack
migration
ie
dw A
where rl = ,/I + [hj k/w]’ form of this integral is
~ ‘ -‘oartn(~,
k),
and llj is the half-offset. A discrete
dj(t, k) = 1 A-‘e
~i2noAr~nrtn(w, k),
w (nt is the number of time samples in a trace). This is (for each
(4)
Prestack partial migration
mid oint
Constant offset 22 section $J
.
NM0 $: Bi DMO r-i
Model d, (t, k) =
1 e - ilnodrlnt
m
1
f
,$ 5
FIG. 3. Decomposition of prestack full migration into NMO, DMO, and poststack migration. An impulse on a commonoffset section can be migrated to the expected ellipse with the shot and the receiver in the foci, either in one process (ptestack full migration), or in three steps (NMO, DMO, and zero-offset migration. which can be done as a poststack process).
976
Ronen
For J = 5 offsets and an aliasing of N = 3, the system is
~~~=~~~i._~~~]~~~~~.
(5) This is the same as equation (2), except d, and m are vectors and &),Tis a matrix.
Equation (5) is a system of J x nt equations (one equation for every time sample of every common-offset section), with N x no unknowns for each wavenumber k. The matrix has J x N blocks with each block Q,? (k - nk,) being an nt x nw matrix. Solving this equation extrapolates the data in the offset direction (to zero offset), and also in the wavenumber direction from the low wavenumbers (k 1< k,/2, on which the aliased data dj are given, to the full range (k 1 < Nk,/2 required to describe the unaliased zero-offset section. This formulation is similar to the multichannel deconvolution given by Davies and Mercado (1968). Both formulations involve inversion of large block matrices, except that here the blocks are inverse-DMO operators instead of convolutions.
Wave-equation trace interpolation in two steps Conventional data processing does not involve inversion of huge matrices, such as the one in equation (5), and missing data are often assumedto be zero. Nevertheless,it works quite well when the data are not severely aliased. If there is no spatial aliasing, then N = 1 in equation (5), and one has
(6)
m-
4
D;
H A/D
k
dj If DMO is applied to each common-offset section, stack = i Qjdj j=l
-
.
. = CQjLlfm ,:
1
.I x
=
1
DMO- I*
sub-
DMO
sampling
--
Commonoffset section
-__I b r2; Ii
Prestack migration
Jm.
(7)
The processing in equation (7) is the application of the Hermitian transpose of the matrix in equation (6),
D,(k)
Q,(k)
0 E
Common offset m section I
Cm j=1
Dh(k)
Q,(k)
1
(the asterisk denotes taking the Hermitian transpose). This application usesthe result p+ *,(,4) = Q, {k) I
FIG. 4. The choice of the model and the data. The model m is the ideal zero-offset section. The data d. are the commonoffset sections after NMO. The Q,’ are &verse-DMO operators.
(8)
given in the Appendix. When the data are not spatially aliased, equation (6) can be solved by applying the Hermitian transpose matrix. Since this is adequate processing for the nonaliased case, it is interesting to see what the results of applying the Hermitian transpose are in general.
Wave-equation
Trace
Interpolation
977
400
600
800
N
FIG. 5. Synthetic example. (a) The model: a well-sampled, zero-offset section. (b) The data: four midpoint gathers after NMO. The midpoint interval is 200 m. Each gather has five traces (at the same midpoint but with different offsets). Note that the dipping reflector is overcorrected by the NMO. (c) Transpose processing: the result of DMO and stacking after addition of zero traces. (d) The result of wave-equation trace interpolation.
I
I
ORIGINRL I
OQTFI I
I
I I
r
SRtiPLED 1
DRTG I
I
I
1
I
1
I
I
I I
1
1
4
shot (4
I
9
5oom
shot
I
E
8
I
@I
6. (a) Stacking chart of the original data. The original data do not have an aliasing problem: 25 m shot interval, 12.5 m group interval (one offset is missing because of a dead receiver). (b) Stacking chart of the data used for the trace interpolation. Less than 1 percent of the data were used, simulating severe spatial aliasing. FIG.
978
Ronen
If equation (5) is written as
Synthetic data example
d=Gm,
(9)
(where G is the J x nt by N x nw matrix), then the vector ti can be defined as m = G*d.
(10)
To find 61,consider the example (N = 3 and J = 4). Then
_+E;/;
;;
;;;;I.
The transpose of G is D,(k - k,) D, (k - k,) !I, (k - k,) D,(k - k,) G* =
Dz (k) !I, (4 D,(k) [@(k + A,, ” &Z(k+!e) 5 &(k+k,)
124(k) ?,(k+kj
The two-step method was applied to a synthetic data set of four CMP gathers, 200 m apart. Each gather had five traces with 0 to 800 m offset (Figure 5b). The earth model includes a 90 degree dipping reflector and a flat reflector. The ideal solution, the well sampled (25 m midpoint interval) zero-offset section, is shown in Figure 5a. This section is not part of the simulated data but is included for comparison with the results of the trace interpolation. The results of applying the transpose operator (DMO and stacking) are shown in Figure SC. Poststack correction, using QR factorization to invert G*G + C, (with a diagonal c), generated the section shown in Figure 5d. The improvement by the poststack correction is substantial, in particular for the Rat reflector. The conjugate gradient method with spectral balancing
1 I “‘J
The two-step method is not easily applicable to field data because the formation and the inversion of G*G are too com-
where equation (8) was used again. KI is therefore
d,(k)
D,(k
- ks)
P,(k)
Dz (k + k,)
Dx (k + k,)
D,(k
for - 3k,/2 < k =c - kJ2
i Dj(k)d,(k) j= 1
for - kJ2 < k < kJ2
i Qj(k)dj(k-ks) j= I
for k,/2 < k < 3k,/2.
(11)
! Equation (1 I) describes DMO (application of Dj) and stacking
using the (unjustified) assumption d(k -
k,) = d(k) = d(k + k,).
P,(k
Da(k)
j=l k(k)=
- k,J
9, (k)
for ) k ) < kJ2. This can be written as
i D,(k)dj(k+kJ
b(k
(12)
This assumption of replicating the Fourier transform is equivalent to adding zero traces between data traces in the space domain. An important conclusion is that G* is equivalent to DMO and stacking, treating missing data as zero data. Applying C* [equation (lo)] can be the first step in solving equation (9). The second step is
[Actually G*G is singular, so the pseudoinverse (Strang, 1980) should be used.] To summarize, the two steps are (1) DMO and stacking with zero traces in place of missing traces [equation (lo)], and (2) a poststack process [equation (13)].
-
IiI
k,)
d, (k) d,(k)
’
+ kJ
4 (k)
putationally intensive. Fortunately, there are practical iterative techniques, such as the conjugate gradient method, that produce good results at a reasonable cost. Conjugate-gradient inversion is an iterative method in which each iteration involves the application of a forward operator and its transpose (Golub and Van-Loan, 1983; Luenberger, 1984; Paige and Saunders, 1982). Conjugate-gradient inversion of an M x M matrix converges within M iterations. Often, just a few iterations (much less than M) are actually needed. To solve equation (9) each conjugate-gradient iteration involves DMO stacking (G*) and transpose DMO stacking (G). Therefore, the total cost in four iterations is the cost of eight DMOs. Actually, the cost is less than that because Fourier transforms are not repeated, data are not resorted, and operators are re-used. When applied to field data, straightforward conjugategradient inversion converged reasonably quickly, but to a disappointing result. The reason for the poor results is that solving equation (9) involves inversion of singular matrices that have no exact inverse. Moreover, even the approximate solution is not unique. There are many models that reasonably fit the data. Those models differ by components of what is called the null space. Null-space components are the eigenvectors that correspond to zero singular values. [See Strang (1980) for a diszussiorr of singular-vaiue~decomposition and null space.] To achieve a unique solution, the null-space components in the model must be estimated. A standard estimation procedure for the null space in leastsquares problems is damping: no null-space components in the solution. Damping was used in the synthetic example to gcneratc Figue 5d, which is of a reasonable quality. For the
Wave-equation
150
300
Interpolation
979
field data, one has to be more careful, because in practice the null space is composed of the eigenvectors that correspond to small (not necessarilyzero) singular values. How small is small is determined by the noise level. The synthetic example is free of noise, so the null space is small and damping produces a reasonable model. Field data, however, have enough noise to require a more sophisticated estimation of the null space. The estimation of the null-space components is an optimization problem, the objective of which is to fit a priori assumptions. A safe a priori assumption is that the spatial spectrum of the earth model should be balanced. A study of spatial abasing in multichannel seismic data (Rocca and Ronen, 1984) predicted that in the presence of spatial aliasing, when performing DMO and stacking there will be resonance wavenumbers for which the spectrum will have high amplitude. An a priori assumption is that these are artifacts. This a priori assumption can be used in an iterative optimization scheme (such as the conjugate gradient method). Each iteration can involve, in addition to the usual conjugategradient routine, the following spectral balancing procedure:
(m) 0
Trace
450
(1) Transforming the raw model into the wavenumber-time (k, r) domain. (2) Calculating the input envelope (a leaky integration of the absolute values will do). The expected envelope is calculated based on the input envelope monotonously decreasing with wavenumber. (3) Dividing the input by its envelope, and multiplying by the expected envelope. Field-data
Fto. i. Tne data used for the trace interpolation. Every trace is corrected for gain and NMO, and is plotted at its midpoint position. As seen on the stacking chart of Figure 6b, there are four groups of five near traces and four groups of far traces. This figure illustrates what would be obtained by NM0 and stacking.
RAW 200
1
RAW 400
200
2
RAW 400
200
3
A small part of a marine line from the Gulf of Mexico was used as a test case. The original data do not have any aliasing because a shot interval of 25 m and 240 channels, 12.5 m apart, provides a 6.25 m midpoint interval on the 6000 percent stack. Only a small number of these data were used: 150 m
RAW 400
example
200
4
RAW 400
200
5
RAW 400
200
0 400
FIG. 8. The models of the first six iterations. Convergence is achieved within four iterations, but the resulting model is disappointing; it is very different from the true zero-offset section of Figure 12.
(m)
980
Ronen
RAW
1
0
RAW
0.05 0
2
RAW
0.05
3
RAW
0.05
4
RAW
0.05
5
RAW
0.05
6
0.05
(cycles/m)
N
FIG. 9. The raw (unbalanced) spectra of the first six iterations. As predicted (Rocca and Ronen, 1984), there are some resonance wavenumbers where the amplitude is too high compared to the true spectrum of Figure 13. BALANCED 1 BALANCED 2 200 400 200 400
BALANCED 3 200 400
BALANCED 4 200 400
BALANCED 200
5
BALANCED
400
200
6 400
0
0-n)
N
FIG. 10. The models of the first six iterations with spectral balance. Compared to Figure 8, these sections are much more similar to the near-offset section in Figure 12. BALANCED
0
0.05
1 BALANCED
2 BALANCED
0.05
3 BALANCED
0.05
4 BALANCED
0.05
5 BALANCED
0.05
6
0.05
(cycles/m)
0
FIG. 11. Spectra of the first six iterations with balancing. The spectral balancing routine used the knowledge that the raw spectra of Figure 9 may have resonance wavenumbers in which the amplitude is too high.
Wave-equation
0-W 200
Trace
981
Interpolation
shot interval (killing five shots out of every six), and ten channels per shot (killing 230 channels out of 240) in two groups, using five near traces and five far traces. This particular selection of data simulates the aliasing that would occur in the cross-line direction of data collected in a two-boat survey with a 150 m line interval where the boats alternate shooting and only one boat is recording (Ronen? 1985). A stacking chart of the data is shown in Figure 6. The data that were used in the inversion are shown in Figure 7. As in many 3-D surveys, the coverage is not uniform. Nonuniform coverage is natural for the data; it is the model, not the data, which should be uniformly and adequately sampled. The models obtained by the first six iterations (Figure 8) are the disappointing results mentioned earlier in this paper. The spatial spectra are clearly unbalanced (Figure 9). The models for the first six iterations, with spectral balancing, are shown in Figure 10, and their spatial spectra are shown in Figure 11. The results with spectral balancing are a much better estimate for the near-offset section (Figure 12) and its spectrum (Figure 13) compared to the results without spectral balancing (Figures S-9). The final result (Figure 15), although not precisely the near-
0-N 200
400
0
FIG 12.The near-offset section. This section at offsets of about 220 m is considered the “true solution.” Only four traces from this section were included in the data of Figure 7.
TRUE
SPECTRUM
0
FIG. 13. The spectrum of the near-offset section of Figure 12. This spectrum was not used in the trace interpolation, but is shown here for comparison with the estimated spectra of Figures 9 and 11.
FIG. 14. The model obtained by the first iteration, without spectral balancing. This would be the result of adding zero traces in place of missing data, then performing DMO and stacking.
R0n0ll
982
offset section (Figure 12), is significantly better than what is achieved by DMO stacking (Figure 14), and tremendously better than the result of NM0 and stacking (Figure 7). CONCLUSIONS
The trace-interpolation method is based on the wave equation and on an a priori assumption of a smooth spatial spectrum. The method involves an iterative inversion of the linear relations d = Gm, where the model vector m is the ideal zerooffset section, and the data vector d is the seismic data after NMO. Treating missing data as zero data and performing DMO and stacking are equivalent to application of the transpose operator I?I = G*d. This processing is adequate only in the absenceof spatial aliasing. 0-n) 200
400
Two alternative methods of trace interpolation presented are a poststack process (on synthetic data) and an iterative, conjugate-gradient scheme (on field data). The poststack process involves large matrix multiplication and inversion and is too costly to be applied to field data. The iterative method, a combination of conjugate-gradient and spectral balancing, was successfully applied to field data. In both methods, the first step, or first iteration, is DMO and stacking. ACKNOWLEDGMENTS
This work would not have been possible without the guidance and support of the staff and sponsors of the Stanford Exploration Project. Special thanks are due to Fabio Rocca who started me on this project; to Jon Claerbout who taught me how to make it work; to Einar Kjartansson, Stew Levin, Rick Ottolini, Chuck Sword, and John Toldi with whom I often discussedraw ideas, and whose software was very useful; to Paige and Saunders for writing LSQR; and to Chevron who donated the field data. Thanks are also due to Jim Black and Dave Hale who carefully reviewed the manuscript.
REFERENCES
Bolondi,G.. Loinger,E., and Rocca,F., 1982,Offset continuationof seismicsections:Geophys.Prosp.,30, 813-828.
FIG. 15. The model obtained at the sixth conjugate-gradient iteration with spectral balancing. This is the estimation for the zero-offset section (Figure 12) that the wave-equation trace interpolation produces. Flat as well as dipping reflectors are well-interpolated.
Brown, J. L., 1981, Multichannel sampling of lowpass signals: Inst. Electr. Electron. Eng., Trans. Circuit Theory, CAS 28, 101-106. Claerbout, J. F., 1984. Imaging the Earth’s interior: Blackwell Scientific Publ. Ltd. Davies E. B., and Mercado, E. J., 1968, Multichannel deconvolution filtering of field recorded seismicdata: Geophysics, 33,711-722. Dereeowski. S. M.. and Rocca, F.. 1981, Geometrical optics and wave theory of constant offset sections in layered media: Geophys. Prosp., 29, 384-406. Golub, G. H., and Van-Loan, C. F., 1983, Matrix computations: The Johns Hopkins Univ. Press. Hale, D., 1983, Dip moveout by Fourier transform: Ph.D. thesis, Stanford Univ. Judson, D. R., Schultz, P. S.. and Sherwood, J. W. C., 1978, Equalizing the stacking velocities of dipping events via Devilish: Presented at the 48th Ann. Internat. Mtg., Sot. Explor. Geophys., San Francisco. Linden, D. A., 1959, A discussion of sampling theorems: Proc., I.R.E., 47, 1219. Luenberger, D. G., 1984, Linear and non linear programming: Addison-Wes!w. Nyquist, H., 192& Certain topics in telegraph transmission theory: AIEE Trans., 47,617-644. Ottolini, R. A., 1982, Migration of reflection seismic data in anglemidpoint coordinates: Ph.D. thesis,StanfordUniv. Paige, C. C., and Saunders,M. A., 1982,LSQR: An algorithm for sparse linear equations and sparse least squares: Assn. Comp. Math., Trans. Math. Soft., 8, 195-209. Papoulis, A., 1977, Generalized sampling expansion: Inst. Electr. Electron. Eng., Trans. Circuit Syst., CAS 24, 652-654. Rocca, F.. and Ronen, J., 1984, Improving resolution by dip moveout: Presented at the 54th Ann. Internat. Mtg., Sot. Explor. Geophys., Atlanta. Ronen, J., 1985, Multichannel inversion: Ph.D. thesis, Stanford Univ. Shannon, C. E., 1949, Communication in the presence of noise: Proc. I.R.E., 37, l&21. Strang, G., 1980, Linear algebra and its applications: Academic Press, Inc. Yilmaz, O., 1979, Prestack partial migration: Ph.D. thesis, Stanford Univ.
Wave-equation Trace Interpolation
983
APPENDIX INVERSE DMO
Different common-offset sections are not independent because they are all reflected from the same reflectors. Many authors (Judson et al., 1978; Yilmaz, 1979; Deregowski and Rocca, 1981; Bolondi et al., 1982; Ottolini, 1982; Hale, 1983) have shown how to get the zero-offset section from commonoffset sections. For the trace interpolation, the inverse operator is needed and it is derived in this Appendix.
Zero-offset migration is described by the dispersion relation
(A-5) (k, , li,, kZ, and WI,are the wavenumbers that correspond to X, I’, Z, and t,,) Inserting the Fourier transform of the change of variables equation (A-4)
Prestack full and partial migration The presence of a single spike in the data set, at midpoint (x, ,3 = (0, 0), and traveltime t,, with a shot at (x, y, z) = (h, 0. 0) and a receiver at (-h, 0, 0), means there is an ellipsoidal reflector such that the distance from any point (x. I’. z) on the ellipsoid to the source plus the distance from that point to the receiver equals the trayeltirne multiplied by the velocity, ~~(x-~,2+~~+Z*+~~~++)~+~~+z~=uth.
(A-l)
Points (x, y, Z) in equation (A-l) are on an ellipsoid which can be written as
J
into equation (A-5) gives
which is the dispersion relation for prestack, post-NMO, full migration. Both t, and w, appear in equation (A-6) because the operator is time-variable. Prestack full migration can be done in three steps (Figure 3): (1) NMO: mapping,
Perform the well-known trace-by-trace
a, can be found by setting 2 = 0 and r‘ = 0 in the last two equations, i.e.,
t; = th’- (2h/a)Z. (2) DMO:
(A-2)
Define (A-7)
Similarly, a, = a, = V/(ut,/2)2 - h2 = ut,/2,
co0is the frequency that corresponds to the time on the zero-offset section. (3) Migration: Using
where t, is the normal moveout (NMO) time r; = t,: - (2h/U)Z.
(A-2)
[gyk:+k;+k:]=lu:.
Therefore, prestack migration of the impulse produces the ellipsoid
(A-8)
obtain the migrated section. Inverse DMO or
(A-3)
x2 1 + (2/l,/nr,)Z
t 4.2f 22 = L
vt,/2 a, 1
d(t,, k,, h) =
The change of variables
1 done-‘“+@,,
k,, h)
J
and using equation (A-7) to change variables o,, to o,, ,
x2
x2 =
1+
(2h,ILQ
(A-4)
compressesthe ellipsoid (A-3) to the sphere X2 + $ + z2 = (ot,/2)2, which is the zero-offset migration of a spike at the NM0 t”
Equation (A-7) can be used to derive the relation between the zero-offset section m and the common-offset section d. Starting from the inverse Fourier transform
time
984
Ronen
Defining
with A = ,/l
+ (hkJw, t,)‘,
A= Jl
one obtains (A-10)
W” = Ao,
To return to the original (x, y) coordinates, recall that a rotation in the space domain corresponds to a rotation in the same direction in the wavenumber domain:
and
4l a0 = A-‘, do, - a,
cos 9 kX’ z-z kY’ -sin0 (1[
(A-11)
sin 0 cos8
k,
I I k,
Therefore,
Substituting equations (A-10) and (A-11) into equation (A-9),
d(t,, k,, h) =
(A-14)
+ (hk,./o, t,)‘.
dw, A-le-‘“OA’nm(oo,
k,),
hk,, = h(cos Ok, +
(A-12)
sin
Ok,)
= (h cos O)k, + (h sin O)k,
s = h,k, + h,ky
which implies a method of generating common-offset sections d(t,, x, h) from the zero-offset m(t, , x).
=h.k. Equations (A-13) and (A-14) become
Inverse DMO in three dimensions In equation (A-l), I made the assumption that the shot and receiver are on the x-axis. In general they are not, and the half-offset vector has two components:
d(t,, k,, k,, I?)=
do, A-le-‘“UA”m(oo,
k,, k,,)
s and A = \/l
+ (h - k/w,t,)‘,
h=
where h. k is the linear product of the offset and the wavenumber vectors.
Let 0 be the angle of rotation such that
Summary
One can rotate the (x, y) axis to (x’, y’) by
[I [ XI
The inverse DMO operator is given by
=
-sin0
Y’
d(r, k, h) =
cost3
so that the shot and receiver are on the x’-axis, and one can make the same development of equations (A-l) through (A- 12) in the (x’, y’) coordinates, obtaining
dw A- e‘ -iWA’m(w, k).
This is the Hermitian transpose to the DMO operator, m(w, k) =
dt Am’e-‘OA’d(t, s
d(t,, k,. , k,, ,
11)=
dw, As
‘e-iooALm(wO
k,,
, k,.),
(A-13)
that was derived by Hale (1983).
k, h),
