testing techniques, and geophysics. The early concept of slant-stack processing for geophysical applications originated with Rieber’s construction of a new reflection system with controlled directional sensitivity (Rieber, 1936), which delayed and added recordings (analog sound tracks) at ten equally spaced ground locations to form steered traces with a range of directions. A similar method called the controlled directional reception (CDR) method was extensively used and developed in Russia (Gardner and Lu, 1991). transforms can be The general theory of Radon found in the book by Deans (1983). The fundamental properties of the Radon transform were examined by Durrani and Bisset (1984). For the linear transform, Chapman (1981) presented the exact forms of the generalized Radon transform pairs for a point source in Cartesian or spherical coordinates, and for a line source in cylindrical coordinates. Thorson and Claerbout (1985) imposed the least-squares error constraint on the reconstructed data and developed a method for velocity-stack and slant-stack stochastic inversion. Hampson (1986) implemented the Thorson-Claerbout method efficiently by a parabolic approximation to residual moveout for the NMO-corrected data in the f-x domain. Hampson (1987) identified his method as the discrete Radon transform (DRT) which was explored by Beylkin (1987). Kostov (1990) derived fast and accurate algorithms for computing the least-square inverse of the slant-stack. The classic and the discrete Radon transforms were examined by Gulunay (1990). Growing interest in the transform in the mid-1970s, especially the effort by Claerbout’s geophysics group at Stanford University, led to several applications of this technique in seismic exploration. Table 1 lists the main transform in seismic exploration. One applications of the of the primary purposes of transforming 2-D seismic records domain is to separate, in the transformed domain, to the coherent events that normally interfere with one another in the original domain of recording (i.e., the x-t domain) (Tatham, 1989). This can be enhanced by deconvolution in the transform domain as explained in a later section.
ABSTRACT New derivations for the conventional linear and transforms in the classic continuous parabolic function domain provide useful insight into the distransformations. For the filtering of uncrete wanted waves such as multiples, the derivation of the transform should define the inverse transform first, and then compute the forward transform. The forward transform usually requires a p-direction deconvolution to improve the resolution in that direction. It aids the wave filtering by improving the separation of events in domain. The p-direction deconvolution is the transrequired for both the linear and curvilinear formations for aperture-limited data. It essentially compensates for the finite length of the array. For the transform, the deconvolution is required parabolic even if the input data have an infinite aperture. For sampled data, the derived transform formulas are identical to the DRT equations obtained by other researchers. Numerical examples are presented to space after decondemonstrate event focusing in volution.
INTRODUCTION The Radon transform is a mathematical technique that ha seen popular usage in seismic data processing and analysis in recent years. It is more commonly known as the transform, the slant-stack technique, or the velocity-stack method. The original mathematical development relating to this technique, but not containing any immediate physical applications, was by Radon in 1917 (Deans, 1983). A similar mathematical method is the Hough transform (Hough, 1962). It has provided a unifying mathematical basis for a large class of image reconstruction (tomography) problems in physics, medicine, astronomy, molecular biology, material science, optics, nuclear magnetic resonance, nondestructive
In this paper, new derivations for the conventional linear transforms in the classic continuous doand parabolic main are presented. They provide powerful insight into the discrete transform. Numerical examples are used to deconvolution. explore event focusing by DERIVATIONS OF THE Conventional (linear)
TRANSFORMS
+ px. The seismic wavefield transforms to a compact space, and for broadband data, hyperbolic region in the reflection events map to ellipses and linear refraction events map to points (Schultz and Claerbout, 1978; Diebold and Stoffa, 1981; Treitel et al., 1982; Tatham, 1984). To obtain the inverse transform in equation the standard back-projection is used
transform dp
There are several ways to construct a transform pair: transform, then obtain the first define the forward transform in both the continuous domain (anainverse lytical) and the discrete domain; or first define the inverse transform, then derive the forward transform formula. It also can be worked in both the analytical (continuous) domain and the discrete domain for the latter case. In this transform pair for section, a new way to construct the continuous data is presented which provides useful insight into the discrete matrix version of the transform pair. Linear
dx’ dp
x’) J
J (3)
Introducing a new function
transform, version 1 (classic) =
Define the forward
e iopx
(4)
transform as follows: then equation (3) can be rewritten as +
=
dx.
(1)
I
dx’
=
In the temporal-frequency domain, this becomes =
dx.
(2)
I
Here and , are the Fourier transforms of and respectively. For seismic functions t) can be exploration applications, the quantity thought of as the common mid-point (CMP) or common is the slant stack of shot-point (CSP) gather; while the input gather with horizontal slowness p (or ray parameter, or Snell’s parameter, or the reciprocal of the horizontal phase velocity) and intercept time Slant stacking is a in equation (1). The description of the integral across integral stands for adding or stacking along the slant line t = Table 1. Applications of the Application Velocity analysis Migration and modeling Inversion Interpretation Plane-wave decomposition Wave separation and/or filtering Noise attenuation
Data interpolation Controlled directional reception & controlled directional source methods
(5) = Here, the “*” stands for convolution with respect to the spatial variable x. Equation (5) represents the standard normal equations for deconvolution (Robinson, 1978) which can be solved by the modified Levinson algorithm. The convolution operator is equivalent to the autocorrelation function in the standard deconvolution equation. It clearly shows that the conventional back-projection is not the inverse transform of the quantity p , [see defining equation If the exact inverse transform is desired, the spatial deconvolution must be . This is performed on the back-projected data similar to Beylkin’s inverse DRT (1987). Now let us find out p transform in seismic exploration. References
Schultz and Claerbout, 1978; Schultz, 1982; Gray and Golden, 1983 Chapman, 1978; Wenzel et al., 1982; Ottolini and Claerbout, 1984; Miller et al., 1987; Ruter, 1987; Hubral, 1991 Clayton and McMechan, 1981; Diebold and Stoffa, 1981; Thorson and Claerbout, 1985 McMechan and Ottolini, 1980; Kennett, 1981; Phinney et al., 1981 Stoffa et al., 1981; Treitel et al., 1982 Tatham, 1984; Harlan et al., 1984; Moon et al., 1986; Benoliel et al., 1987; Hu and McMechan, 1987; Greenhalgh et al., 1990 Alam and Austin, 1981; Tatham, 1984; Carrion, 1986; Noponen and Keeney, 1986; Hampson, 1986 & 1987; Yilmaz, 1989; Mitchell and Kelamis, 1990; Foster and Mosher, 1992; Zhou and Greenhalgh, 1993 Lu, 1985 Zavalishin, 1982; Taner et al., 1991; Sword, 1991
.
Transforms Revisited
the function defined by equation (4). There are two cases to be considered: (1) variable p occupies an infinite range of and (2) variable p is of a limited range of Case 2 is more practical because a computer can handle only limited digital number data. However, case 2 can be included in case 1 if it is assumed that there is no useful information for reconstructing the original data outside the range for parameter p. If this holds, then the two situations should produce equivalent results for the same data set. Letting = and substituting p = into equation for case 1, equation (4) becomes =
=
(6)
II The convolution operator respect to the spatial variable pair for the delta function
is a delta function with Here the Fourier transform
has been used. Substituting equation (6) into equation (5), we have =
I
=
(8)
II
(9) Therefore the conventional linear frequency domain is
transform pair in the
= I . =
II
(10)
dp
Equation (10) is widely referred to in the relevant literature. Now let us consider case 2; viz, the restricted domain for p. In this situation, equation (4) becomes =
max
e iopx
1135
1 = 0’ Pmax (11) In this situation, x, is no longer the delta function as shown by equation (6). When = then equation (11) becomes a function
In the limit as approaches infinity, equation (12) reverts to a delta function in x. Therefore, the normal equation (5) must be solved; i.e., the spatial deconvolution must be performed if an exact inverse is desired. From the above derivation of the convolution operators for the infinite and the following observations can be made: finite ranges of 1) The spatial deconvolution is required because the infinite range of p is truncated to a finite range. It is clear x, is of complex conjugate symmetry with that variable x. Therefore, if equation (5) is rewritten in will be of matrix form, the matrix formed by Hermitian-Toeplitz structure (see Appendix A). The result is the same as the inverse DRT of Beylkin (1987), Kostov (1990), Gulunay (1990) and Foster and Mosher (1992), but it is derived in a totally different way. The transform is derived in the continuous inverse while the function domain by limiting the range of other researchers computed the inverse in the discrete domain with the least-squares matrix inversion method. 2) The convolution operator for the finite range of p decreases in amplitude according to the factor l/x. This is clearly shown by equations (11) or (12). When approaches and becomes both equations (11) and (12) will be infinite. Therefore, the convolution operator for case 2 (limited p range) becomes a delta function like case 1 (infinite p range). This observation is similar to that of Gulunay but it is easily and clearly identified by equations (11) and (12). 3) Fast matrix solvers (e.g., the Levinson algorithm) can be used to solve the deconvolution equation (5) to obtain the exact inverse transform. It is worth noting that elements of the Toeplitz matrix for the Levinson algorithm are related to the minimum and maximum p values only and do not depend on how the p In p values are sampled in the range of [p transform, the inverse the conventional discrete transform can be derived directly by discretising equation (5) for the second case of limited p values. of the Toeplitz matrix will be The elements calculated by summation rather than the analytical solution (11). The computation of the summation depends on the sampling rate of the p value. Therefore, the analytical convolution operator (11) guarantees the efficiency in constructing the Toeplitz matrix for the
Zhou and Greenhalgh
11 36
Levinson algorithm to perform the x-direction deconvolution. 4) If the wavenumber in the original data lies out of the range 1 numerical instability arises in solving equation (5) (see Appendix B). This occurs especially for low-frequency components. To stabilize the Levinson algorithm for solving the deconvolution equation (5), prewhitening the convolution operator; i.e., adding a small value into equation (12) is necessary and is common practice in standard seismic deconvolution processing. When the prewhitening noise is large, the operator described by equations (11) and (12) will approach a delta function. Therefore, the inverse transform result for case 2 will approach that of case 1. It is possible to use other methods to stabilize the Levinson algorithm (Kostov, 1990). Linear
transform, version 2
In a similar but alternative fashion to the transform definition of version 1, we can start with the inverse transform and define it as
forward-projected data p, w). There are two situations to be considered for the function p, defined by equation (16): (1) the spatial variable is of infinite range (2) the spatial variable is of finite range [ As stated before, the second case can be included in the first case if there is no useful information for constructing the forward transform data outside of the range [ for variable x. In this situation, the two cases should be equivalent. Letting = and substituting = into equation (16) for case 1, we have 00
=
(18) Here, equation (7) has been used. Substituting equation (18) into equation (17), we obtain
00 =
(13)
I -CO
=
In the frequency domain, this is equal to dp.
(14)
Therefore, the transform pair in the frequency domain can be written as II
Here t) and p, are the input data and the transform of the input data, respectively, and x, and p , are the Fourier transforms of t) and (p, To obtain the forward transform p, in a forward-projection p, is defined as equation =
-al
II
dx I (21)
=
dp
which is quite similar to the transform pair (10). The only difference is that the frequencies attenuated in the inverse transform are compensated for in the forward transform. Now let us consider case viz, the finite aperture array. In this situation, equation ( 16) becomes:
dx
I (15) If a new function
J x min (16)
1 otherwise*
is introduced, equation (15) can be rewritten as =
(22)
(17)
stands for convolution with respect to variHere, the able slowness p. Again, equation (17) represents the standard normal equations for deconvolution which can be solved by the modified Levinson algorithm (Robinson, 1978). It clearly shows that the forward-projection equation (15) and the transform of the defining equation (14) are not a inverse transform pair. If the exact forward transform corresponding to the inverse transform (14) is required, then the p-direction deconvolution must be performed on the
In this situation, h(p, is no longer the delta function as function as in case 2 of shown by equation (18), but is a h( p, can version 1. This is because when be written as
= .2x max
1137
Transforms Revisited
Therefore, the normal equations (17) must be solved; i.e., the p-direction deconvolution must be performed. Equations (22) show that the filter is related to the aperture, or of but independent of the sample number in the range. Significant computer time savings are possible using equations (22) and (23) to construct the Toeplitz matrix for the Levinson algorithm to perform the deconvolution. Similar conclusions as were drawn for version 1 can be made from the above derivation. Comparing the two versions of the linear transforms, it is clear that the deconvolution for version 1 is performed in the inverse transform in the x-direction, while the deconvolution for version 2 is performed in the forward transform in the p-direction (Foster and Mosher, 1992). Deconvolution can improve resolution in the p- or x-direction. Therefore, transform will have better resolution or version 2 of the separation of waves than version 1 for the forward transform. It is preferable to use version 2 for noise attenuation in the domain. From the above derivation and discussion, it is clear that if the input array is aperture limited, deconvolution must be performed in the destination (forward or inverse transform) domain to improve the transform quality. In practice, both and p are of limited range. Therefore in theory, deconvolution should be carried out in both forward and inverse transforms. It transpires that we actually need to perform only the p-direction deconvolution in the forward transform because the range of p can be chosen in such a way that the major energy/information of the input data falls within this p range. Parabolic
transform
On common midpoint (CMP) or common shot-point (CSP) gathers, seismic events are more likely to be hyperbolic than linear (reflections and diffractions are hyperbolic, refractions and direct waves are linear). The standard linear transform converts hyperbolas in the x-t domain into ellipses in the domain. This kind of transform does not really help much for suppressing noise in the form of multiples or for separating reflected P- and S-waves which also have hyperbolic trajectories in the x-t space. Instead of the linear transform, a hyperbolic transform can be designed. It transforms the hyperbolic events on the seismic section into points in much the same way that the linear transform converts linear events into points. This focusing of energy improves the resolution of the section by separating events. In particular, the multiples can be more easily isolated from the signals of interest. However, the direct hyperbolic transform is too expensive to realize because the transform is time variable (Hampson, 1986). To make the transform (Yilmaz, 1989) or NMO cortime-invariant, a rection (Hampson, 1986) is performed on the seismic data to make the events become parabolic. Then a parabolic transform is applied to the -stretched or NMO-corrected data set, and a better focused section can be obtained for suppressing the coherent noise. Foster and Mosher (1992) achieved a generalized transform that provided good resolution in the transform domain by integrating/stacking data along a suboptimal hyperbolic surface. The following is transform. a derivation of the parabolic (curvilinear)
To obtain superior image resolution in the transform domain, a similar approach as was used with the linear transform, version 2, is taken in the following derivation of transform. First, we define the inverse the parabolic parabolic transform as (24)
=
In the frequency domain, this trans-
where forms =
(25) I
Equation (24) shows that the inverse parabolic T-p transform is computed by integrating along a straight line in the transform domain. Here p is no longer the slowness as transform. It can be defined by the conventional linear positive or negative in a mathematical sense. The physical meaning of the parameter p depends on the input data in the x-t domain. If the input data have been NMO-corrected, then p is taken to be the reciprocal of the square of the (Hampson, 1986). If the input is the residual velocity -stretched data, then p is taken as the square of the (Yilmaz, 1989). To obtain the forhorizontal slowness ward parabolic T-p transform, a forward projection function (26)
=
is introduced . Substituting equation (25) into the above equation, we obtain
I (27)
=
Here the convolution operator with respect to variable p is =
e -CO
(28)
which takes the form of a delta function = 0, In common with the second version of the linear transform described by equation (17), equation (27) represents the standard normal equations for deconvolution with respect to variable p which can be solved by the fast Levinson algorithm. To get an exact forward parabolic transform corresponding to the inverse transform (25), the p-direction deconvolution must be performed. Again, there are two cases that must be considered for the computation of the convolution operator (28): (1) the integration bounds are is in the range and (2) the infinite; i.e., integration interval is finite, or lies in the range [
Zhou and Greenhalgh
1138
For the first situation, the convolution operator defined by equation (28) becomes =2
=
Jo
[cos I0
+ =
sign (p) sin +
sign (p)]
(29) l This equation indicates that p-direction deconvolution is required even if the input data have an infinite spatial extent. transform, version 2. In This is different from the linear the Fourier transform domain, equation (27) becomes (30)
= or =
(31) Here the new function in equation (31) is defined as =
=
e
(35)
unless = 0, in which case 0) Note that other more general functions of (e.g., hyperbolic) can be used in place of the parabolic formula in the exponent of the integrand, but they will not be considered here (see Foster and Mosher, 1992). This integral must be approximated by numerical quadrature methods such as the rectangular, trapezoidal or Simpson rules. With the simple rectangular integration rule, equation (35) can be expressed as
+ i sin
0 2
max JX min
[cos
2
For case 2; i.e., the range of variable is finite, there is no analytical solution for the integral
[ 1 + i sign (p)]
dp
+ sign
(32) P The derivation is given in Appendix C. Here the notation for the Fourier transform is a problem. It is difficult to maintain the usual convention of using lower-case letters for the domain of physical space and upper-case letters for the Fourier domain because the convention cannot include the mixed objects with physical and transformed variables. Rather than invent some new notation, it seems best to let the arguments of the function aid in the naming of the function. It is easy to demonstrate that the Fourier transform variable corresponding to variable p is positive. Therefore, equation (32) can be rewritten as =
(33) P
and equation (31) becomes
e iopx .
Ax
(36)
+ (k If integral (35) is approximated by equation (36), the paratransform constructed here will be similar to the bolic algorithms described by Hampson (1986), Beylkin (1987), Gulunay (1990), Kostov (1990), and Foster and Mosher (1992). Indeed, Gulunay (1990) derived an expression similar to equation (29) using Fresnel integrals which was an upper-limit to equation (36) for the discrete case of the parabolic transform. Our approach is distinctly different from theirs. Figure 1 is a flow chart for the parabolic forward and transformations based on the above derivation. inverse transform based on equation (26) is called the The low-resolution transform. If the p-direction deconvolution of such a transform is performed based on equations (27) and (29), then a medium resolution result is obtained. It is referred to as analytical matrix deconvolution. If deconvotransform [equation (26)] is carried out by lution of the equation (34), it is called the Fourier transform (FT) deconvolution. The analytical matrix deconvolution and the FT deconvolution are equivalent in theory, but they are realized in different domains. The result obtained by the FT deconvolution method may differ from that of the analytical matrix deconvolution method because of differences arising from numerical discretization and truncation. The choice of the diagonal elements of the matrix also affects the result. In this = 0, in equation (29) was set to paper, the value of be much larger (ten times) than the neighboring values. As 0, the function a becomes like a delta function, and p deconvolution is no longer required. If the deconvolution is performed according to equations (27) and (34), then a high resolution result is obtained. It is referred to as the standard matrix deconvolution. Gulunay (1990) notes that the hightransforms resolution appearance of forward parabolic result from mainlobe sharpening when the problem at hand is not singular and sidelobe suppression when singularities occur. A similar flow chart to Figure 1 can be formed for the transformation. linear NUMERICAL EXAMPLES
(34) which shows that the p-direction deconvolution enhances the Fourier transform components. Therefore the resolution transform domain is increased. in the parabolic
The following numerical examples are used to demonstrate the necessity of p-direction deconvolution in the transforms when the input array has a relatively small aperture.
Transforms Revisited
Linear
transform
Figure 2 shows a horizontal, linear event at a time of 0.5 s for testing the linear transform. The offset of the first trace is -70 m. The trace interval is 10 m. The input array has a small aperture. The wavelet on each trace is a 25Hz, zero-phase Ricker wavelet with equal amplitude. The extransform of an idealized delta function pected linear version of this event, recorded on an infinitely long array, is a point at p = 0 and = 0.5 s. Figure 3a is the discrete linear transform of the seismic section shown in Figure 2. It was computed using equation (15). No p-direction deconvolution has been performed on this section. The number of p values is 50, over the p range [0, l/500]. The expected point on this section is smeared with linear artifacts because of the relatively small aperture of the input data. The smearing is dramatically decreased after applying the p-direction deconvolution to Figure 3a based on equations (17) and (22). The improvement of the p-direction deconvolution in event focusing is obvious. Figure 4 shows the inverse linear transforms of Figure 3. The inverse transform recovers the linear event in both cases. However, the recovered linear event in Figure 4a, extracted from the data of Figure 3a without p-direction deconvolution, is distorted in phase. By contrast, the event of Figure 4b, which was recovered from the data of Figure 3b after the p-direction deconvolution using a white noise level of 0.003 percent, closely matches the original seismic section (Figure 2). Therefore, the p-direction deconvolution not only improves the event focusing or resolution, but also helps in recovering the original waveshape. Certainly, the deconvolution operator compen-
FIG. 1. A flow chart for parabolic
1139
sates for the lost frequencies during the transform, as suggested by equation (22). Hyperbolic
transform
Figure 5b is a common shot-point (CSP) gather generated for the simple five-layer model shown in Figure 5a. The seismograms were calculated using a nine-point convolutional differentiator method (Zhou et al., 1993). The sampling rates for both space and time used in the modeling are AX = AZ = 10 m and At = 0.5 ms. The seismograms look complicated compared with the relatively simple model (Figure 5a) because of multiple reflections. From the traveltime, we can easily identify the four primary reflections from , and the sea-bottom B and subbottom reflectors The water-layer reverberations corresponding to the reflecand which are lst-, 2nd-, tion B are 3rd-, and 4th-order surface multiples, respectively. The symbols M , and on the seismic section stand for the lst-, 2nd-, and 3rd-order peg-leg multiples from the reflector . The peg-leg multiples from reflector are M and Only the lst-order multiple R M from can be observed in the recorded time the reflector window. It is worth mentioning that the model (Figure 5a) is specially designed to create reflections such that the following situations are combined on one section: (1) isolated primary reflections which do not interfere with multiples such as B and R (2) a primary reflection overlapping a and (3) a multiple on the near-offset traces (e.g., R primary reflection overlapping a multiple on the far-offset traces (e.g.,
processing. A high-resolution transform is achieved by p-direction deconvolution at each frequency.
1140
Zhou and Greenhalgh
Figure 6 shows the parabolic transform results for the seismic section given in Figure 5b. A t2-stretch operation (Yilmaz, 1989) was performed before the actual parabolic transforms were carried out to make the hyperbolic events parabolic. The number of p values is 150, over the range [-5E-8, 1/(1.5E6)]. This is a sufficient p sampling interval. Questions of spatial and slowness sampling and aliasing considerations are beyond the scope of this paper and are treated by Turner (1990). Figure 6a is the low-resolution transform result; i.e., no p-direction deconvolution was performed on the result. The energy for each event dies out slowly in the p-direction, which means the resolution or the focusing of energy in this direction is relatively poor. The resolution is improved by the p-direction FT deconvolution based on equation (34). The result is shown in Figure 6b. The deep reflector energy becomes relatively weak compared with Figure 6a because the resolution of the shallow reflectors is improved much more than that of the deep reflections. Therefore, the energy becomes much higher than prior to deconvolution. The energy of the shallow events is originally stronger than that of the deep events. The analytical and standard matrix deconvolutions in the p-direction for Figure 6a are shown in Figures 6c and 6d. The white noise level added for the matrix deconvolution is 0.1 percent. The result of the standard deconvolution is slightly noisy because the added white noise for the deconvolution is relatively low. In this domain, it is difficult to say which deconvolution method is best for the improvement of resolution. However, it is clear that the events are well separated from each other. Therefore, this domain provides a suitable space in which to remove unwanted waves such as multiples. An inverse parabolic transformation, followed by an inverse t2-stretch, was performed on the sections of Figure 6. The results are presented in Figure 7. The inverse transform result for the low-resolution section of Figure 6a is shown in Figure 7a. The high-frequency components are obviously lost during the transformation, although all the events in the original input data are preserved. The inverse transform result (Figure 7b) from the FT deconvolution section of Figure 6b recovers the high-frequency compo-
FIG. 2. A horizontal linear event at a time of 0.5 s for testing the linear transform. The desired linear transform is a point at p = 0 and = 0.5 s.
nents, but the energy for the deep reflections is attenuated by the transformation. The situation is improved by the inverse transform results for the analytical matrix (Figure 7c) and standard matrix (Figure 7d) deconvolution. Both yield almost identical results to the original input data of Figure 5b. The difference between these two sections is minor. However, the computer time needed in forming the HermitianToeplitz matrix for these two deconvolution methods differs slightly (the analytical matrix method uses about 21 percent less computer time than the standard matrix method). CONCLUSIONS In this paper, new derivations for the conventional linear and parabolic transforms were presented in the classic continuous function domain. The derivations provide useful insight into the discrete transformations. For the filtering of unwanted waves such as multiples, the derivation of the
FIG. 3. The linear transrorms of Figure 2. The range of p is [0, 1/500]. Diagram (a) is the result without p-direction deconvolution. The desired transform point is smeared with linear artifacts because of the relatively small aperture of the input data. Diagram (b) is the result with p-direction deconvolution. The point smearing in (a) has significantly decreased. Therefore, a better event focusing is achieved by the deconvolution.
Transforms Revisited
transform should define the inverse transform first and then compute the forward transform. The forward transform usually requires a p-direction deconvolution to improve the resolution in that direction and aids the wave filtering by achieving a better separation of events. The p-direction deconvolution is required for both the linear and curvilinear transformations for aperture-limited data to compensate for the finite length of the array. For the parabolic transform, the deconvolution is required even if the input data have an infinite aperture. In the discrete domain, the transform formulas derived in this paper are identical to the DRT explored by Beylkin (1987), Hampson (1986, 1987), Kostov (1990), Gulunay (1990) and Foster and Mosher (1992), as explained in Appendix A. However, our versions of the transform are derived in a totally different way. The convolution operators derived here guarantee the computational efficiency in constructing the Toeplitz matrix for the Levinson algorithm to perform the deconvolution.
FIG. 4. The inverse transforms of Figure 3 using equation (13). Diagram (a) is derived from the data in Figure 3a. The inverse transform recovers the linear event in Figure 2 but with serious phase-distortion at both sides of the section. Diagram (b) is derived from Figure 3b. The result is identical to the original input data in Figure 2.
1141
ACKNOWLEDGMENTS The research reported in this paper was supported by grants from the Australian Research Council and Western Mining Corporation. We wish to thank Dr G. Jackson, Dr C. Mosher, and the anonymous referees for their constructive comments. REFERENCES Alam, A., and Austin, J., 1981, Suppression of multiples using slant stacks: 5 1 st Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts. Benoliel, S. D., Schneider, W. A., and Shurtleff, R. N., 1987, Frequency wavenumber approach of the tau-p transform-Some applications in seismic data processing: Geophys. Prosp., 35, 5 17-538. Beylkin, G., 1987, Discrete Radon transform: IEEE transactions on acoustic, speech and signal processing, 35, 162-172. Carrion, P. M., 1986, A layer-stripping technique for the suppression of water-bottom multiple reflections: Geophys. Prosp., 34, 330-342. Chapman, C. H., 1978, A new method for computing synthetic seismograms: Geophys. J. Roy. Astr. Soc., 54, 481-518. - 1981, Generalized Radon transforms and slant stacks: Geophys. J. Roy. Astr. Soc., 66, 445-453. Clayton, R. W., and McMechan, G. A., 1981, Inversion of refraction data by wavefield continuation: Geophysics, 46, 860-868. Deans, S. R., 1983, The Radon transform and some of its applications: J. Wiley & Sons, Inc. Diebold, J. B., and Stoffa, P. L., 1981, The traveltime equation tau-p mapping and inversion of common midpoint data: Geophysics, 46, 238-254. Durrani, T. S., and Bisset, D., 1984, The Radon transform and its properties: Geophysics, 49, 1180-l 187. Foster, D. J., and Mosher, C. C., 1992, Suppression of multiple reflections using the Radon transform: Geophysics, 57, 386-395. Gardner, G. H. F., and Lu, L., 1991, Editors’ introduction to “Initial concepts’’ : Slant-stack processing, Gardner, G. H. F., and Lu, L., Eds., l-2: Soc. Expl. Geophys. Gray, W. C., and Golden, J. E., 1983, Velocity determination in a complex earth: 53rd Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts. Greenhalgh, S. A., Mason, I. M., Lucas, E., Pant, D., and Eames, R. T., 1990, Controlled direction reception filtering of P- and S-waves in space: Geophys. J. Int., 100, 221-234. Gulunay, N., 1990, F-X domain least-squares tau-p and tau-q: 60th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1607-1610. Hampson, D., 1986, Inverse velocity stacking for multiple elimination: J. Can. Soc. Expl. Geophys., 22, 44-55. - 1987, The discrete Radon transform: A new tool for image enhancement and noise suppression: 57th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts. Harlan, W. S., Claerbout, J. F., and Rocca, F., 1984, Signal/noise separation and velocity estimation: Geophysics, 49, 1869-1880. Hough, P. V. C., 1962, Method and means for recognizing complex patterns: U.S. Patent 3,069,654. Hu, L. Z., and McMechan, G. A., 1987, Wavefield transformations of vertical seismic profiles: Geophysics, 52, 307-321. Hubral, P., 1991, Slant-stack migration, in Gardner, G. H. F., and Lu, L., Ed., Slant-stack processing: Soc. Expl. Geophys., 111117. Kennett, B. L. N., 1981, Slowness techniques in seismic interpretation: J. Geophys. Res., 86, 11 575-11 584. Kostov, C., 1990, Toeplitz structure in slant-stack inversion: 60th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1618-1621. Lu, L., 1985, Application of local slant stack to trace interpolation: 55th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts. McMechan, G. A., and Ottolini, R., 1980, Direct observation of a curve in a slant-stacked wavefield: Bull. Seis. Soc. Am., 70, 775-789. Miller, D., Oristaglio, M., and Beylkin, G., 1987, A new slant on seismic imaging-Migration and integral geometry: Geophysics, 52, 943-964. Mitchell, A. R., and Kelamis, P. G., 1990, Efficient tau-p hyperbolic velocity filtering: Geophysics, 55, 619-625. Moon, W., Carswell, A., Tang, R., and Dilliston, C., 1986, Radon transform wavefield separation for vertical seismic profiling data: Geophysics, 51, 940-947. Noponen, I., and Keeney, J., 1986, Attenuation of waterborne
1142
Zhou and Greenhalgh
FIG. 5. (a) A simple five-layer shallow water structural model for numerical modeling computation of CSP seismograms. (b) A CSP seismic section computed using a nine-point convolutional differentiator. The direct wave from the source has been muted before display. Symbol definitions are: B = seabottom, R = reflection from sub-reflectors, M = multiple, number = reflector number or order of multiple.
Transforms Revisited
FIG. 6. The parabolic
transformation of the section in Figure 5b with various p-direction deconvolution 2 transform. (a) No deconvolution; (b) FT methods. The input data were t - stretched before the deconvolution; (c) analytical matrix deconvolution; (d) standard matrix deconvolution. The T-axis is in t2 sampling number.
1143
1144
Zhou and Greenhalgh
FIG. 6. (continued)
Transforms Revisited
1145
FIG 7 The inverse parabolic transformation of the sections in Figure 6. The inverse t2 -stretch was performed on the transformed data before being displayed. (a) from Figure 6a; (b) from Figure 6b; (c) from Figure 6c; (d) from Figure 6d.
1146
Zhou and Greenhalgh
FIG. 7. (continued)
1147
Transforms Revisited coherent noise by application of hyperbolic velocity filtering during the tau-p transform: Geophysics, 51, 20-33. Ottolini, R., and Claerbout, J. F., 1984, The migration of commonmidpoint slant stacks: Geophysics, 49, 237-249. Phinney, R. A., Chowdhury, K. R., and Frazer, L. N., 1981, Transformation and analysis of record sections: J. Geophys. Res., 86, 359-377. Rieber, F., 1936, A new reflection system with controlled directional sensitivity: Geophysics, 1, 97-106. Robinson, E. A., 1978, Multichannel time series analysis with digital computer programs: Holden-Day . Ruter , H., 1987, Migration and tomography-Radon migration: First Break, 5, 399402. Schultz, P. S., 1982, A method for direct estimation of interval velocities: Geophysics, 47, 1657-l671. Schultz, P. S., and Claerbout, J. F., 1978, Velocity estimation and downward-continuation by wavefront synthesis: Geophysics, 43, 691-714. Stoffa, P. L., Buhl, P., Diebold, J. B., and Wenzel, F., 1981, Direct mapping of seismic data to the domain of intercept time and ray parameter-A plane-wave decomposition: Geophysics, 46, 255267. Sword, C., 1991, The CDR method: in Gardner, G. H. F. and Lu, L., Eds., Slant-stack processing: Soc. Expl. Geophys., 365-389. Taner, M. T., Baysal, E., and Koehler, F., 1991, The controlled directional source method, in Gardner, G. H. F., and
Tatham, R. H., 1984, Multidimensional filtering of seismic data: Proc. IEEE, 72, 1357-1369. - 1989, Tau-p filtering, in Stoffa, P. L., Tau-p: A plane-wave approach to the analysis of seismic data: Kluwer Academic Publ., 35-70. Thorson, J. R., and Claerbout, J. F., 1985, Velocity stack and slant-stochastic inversion: Geophysics, 50, 2727-2741. Treitel, S., Gutowski, P. R., and Wagner, D. E., 1982, Plane-wave decomposition of seismograms: Geophysics, 47, 1375-1401. Turner, G., 1990, Aliasing in the tau-p transform and the removal of spatial aliased coherent noise: Geophysics, 55, 1496-1503. Wenzel, F., Stoffa, P. L., and Buhl, P., 1982, Seismic modeling in the domain of intercept time and ray parameter: IEEE, ASSP, 33, 407-423. Yilmaz, O., 1989, Velocity stack processing: Geophys. Prosp., 37, 357-382. Zavalishin, B. R., 1982, Improvements in constructing seismic images using the method of controlled directional reception: 52nd Annual Intemat. Mtg., Soc. Expl. Geophys., Expanded Abstracts. Zhou, B., Greenhalgh, S. A., and Zhe, J., 1993, Numerical seismogram calculations for inhomogeneous media using a short, variable-length convolutional differentiator: Geophys. Prosp., 41, 751-766.
APPENDIX A MATRIX REPRESENTATION OF THE LINEAR In this paper, the transformation pairs were derived analytically assuming a continuous function. If the results are presented in the discrete domain matrix form, they are identical to those of Beylkin (1987), Kostov (1990), Gulunay (1990), and Foster and Mosher (1992). This can be seen from the following derivation of the matrix representation of the linear transform, version 1. In the discrete domain, equation (2) can be expressed by
TRANSFORM, VERSION 1
(A-4) The combination of the discrete form for equations (3) and (5) is
(A-5) which can be rewritten in a concise matrix form as For j = 1, 2, . . . , Np,
(A-1)
where Nx and Np are the total sampling numbers along the x- and p- directions, respectively. The apertures for the x and p variables are limited in this case. Equation (A-l) can be written in matrix form
(A-6) or = where the elements for the new matrix H are =
=
(A-2)
where the elements of the Nx x Np dimensional phase-shift matrix are
(A-3)
The vectors V and U are
(A-7)
(A-8)
The back-projection operator in equations (A-6) and (A-7) is simply the Hermitian conjugate of There are two ways to form the matrix either using the discrete version of equation (4) or the integration result of equation (11). Based on the direct discretization of equation (4) for the limited range of p values, equation (A-8) becomes (A-9) which is the standard way to form the matrix. If the integration result, equations (5)-(11), is used, equation (A-8) can be rewritten as
Zhou and Greenhalgh
1148
(A-10)
=
Equation (A-10) is the accurate expression of equation (A-9), which can be understood from the derivation. It guarantees the computation efficiency in constructing the matrix while the efficiency of equation (A-9) depends on the way of sampling. The matrix has complex conjugate symmetry. has a HermitianTherefore the convolution matrix Toeplitz structure, and equation (A-6) can be solved by the fast Levinson algorithm. It is obvious that the relationship between and in equation (A-9) is
(A-l 1) which indicates that equation (A-6) or equation (A-7), based on equation (A-9), is identical to the result obtained by Beylkin (1987), Kostov (1990), Gulunay (1990), and Foster and Mosher (1992). However, the derivation is totally different. In a similar fashion, the matrix representations for the linear transform, version 2, and the parabolic transform can be derived.
APPENDIX B INSTABILITY OF THE SPATIAL DECONVOLUTION FOR EQUATION (5) An instability problem in solving the spatial deconvolution equation (5) arises if the Fourier transformation with respect to variable x of the convolution has any zero value. This occurs only for the operator situation of a limited range of p. In the Fourier transform becomes domain, the convolution operator
Therefore, the spatial deconvolution algorithm for solving equation (5) will be unstable if the wavenumber in the original data lies outside the range The evaluation of possible instabilities with the deconvolution algorithm arising in other transform constructions can be carried out in a similar fashion.
APPENDIX C DERIVATION OF EQUATION (32)
Transforms Revisited
1149
In the above derivation, the following formulas have been used
and
Comments
Report "Linear and parabolic τ- p transforms revisited "