GB2385435A - Amplitude and phase surface consistent deconvolution - Google Patents

Abstract

Surface consistent deconvolution is commonly applied to seismic data, based on the assumption that the seismic trace consists of the convolution of wavelets which are consistent for sources, receivers, offsets and common mid-points for all the data. These wavelets are usually separated in the Log-Fourier domain, where the convolution becomes a summation. This produces an amplitude spectra. But the phase spectra is usually derived by other means as the phase spectra is wrapped between Ò f . Assumptions about the phase spectra may be made to derive the phase spectra from the amplitude spectra, or some ad hoc process of unwrapping the phase so that the phase can be partitioned between the four components. A method is described here where the phase is not unwrapped but transformed so that it can be operated on by the same matrix operator generated for the amplitude spectra. Thus, solutions for the phase spectra for the four components are produced without relying on assumptions and with little extra computational expense. The method is demonstrated on some example ultrasonic data on a concrete block.

Description

Amplitude and Phase, Surface Consistent Deconvolution Introduction

Surface consistent deconvolution is based on the commonly used premise that the seismic trace is the convolution of four components, a source operator, a receiver operator, a sub-surface reflectivity operator and an offset related operator. These operators are consistent between different traces to which they contribute. This is a very similar concept to surface consistent statics e.g., Taner et al (1974). With sufficient data the equations can be inverted to obtain the amplitude spectra. The usual approach, then, is to assume minimum phase and derive operators on this assumption, e.g., Morley and Claerbout (1983), Cambois and Stoffa, (1992). The phase function is discontinuous and the inversion process needs continuity to converge on the correct solutions. Unwrapping the phase by various means has been attempted (Poggiagliolrni et al, 1982, Spagnolini, 1993, Cambois and Stoffa, 1993) but it is a computationally costly and ad hoc process. This is a similar problem to side-scan sonar interferometers where an ambiguous phase is used to calculate a beam angle and hence depth.

Four component decomposition The four-component surface-consistent hypothesis can be written as dij (a)) = s; (o).g; (lo) hk (do) YE (lo) Equation 1 where dij(a7)is the Fourier transformed date for source land receiver j, sj(a>) is the source spectra, gj(a>) is the receiver spectra, hk (a), k = li - Oil is the offset spectra, and Yl (a>), l = (i + j)/2 are the residual CMP (common mid-point) spectra. This can be derived from consideration of the source and receiver spectra convolved with the Green's function of the response. In seismic surveys there is some question of reciprocity of source and receivers due to near surface anomalies and non-linear responses to sources. The convolutions are linearised by taking logarithms of the spectra, (Cambois and Stoffa, 1993, Cary and Lorentz, 1993). Making substitutions like Si = ln(si(a')) Equation becomes Dj; (is) = Si (is) + Gj (is) + Hk (a') + Yl (m) Equation 2

The real part.of the logarithm of the spectra.is the amplitude, and the imaginary part is the principal part of the phase; (i.e. the phase wraps around between +).

Let NS, NG, NH, NC and NY be the number of receivers, sources, offsets, channels and common midpoints, respectively.:Equation then represents a system of NS x NC linear equations in NS + NG + NTI + IVY unknowns. I have also dropped the subtraction of the average spectra introduced by Morley and Claerbout (1983).

FOI a modern seismic survey this becomes a very large system of equations and Cary and Lorentz (op. cit.) offer a Gauss-Seidel method requiring only two passes of the data tapes (as the volume of prestack data is usually too large to fit onto the available disk space).

For my test data set I considered this approach as unnecessary. I reformulate Equation 2 as a mask or selector matrix operating on the component spectra to produce the observed data.

Si (1) À Si (n) :::: G1(1) Gj(77) try(l) Dy(n): ...... . JIG (1) Hk (n) =. Equation 3 ....:::: D,j(l) D'j(77) Ye (1). Fi(n) :::: Each row of the component matrix and data matrix is a spectra for n frequencies from O to COmaX Making the substitution of the matrices S. G. H and Y for Si, Gj, Hk and F' respectively, A for the selector matrix, and the data as D gives S Equation 4 G A _ =n Y The augmented matrix is assembled from the matrices S. G. lI and. Making this substitution gives AX = D Equation 5 Damped least squares solution Some elementary matrix algebra on Equation gives = (ATA) ATI) Equation6 There is much discussion in the literature on whether this is the best inverse, but I use it here, it is the least squares solution. The matrix is independent of the data and consists only of ones and zeros so numerical accuracy should not be a problem. The matrix A contains (NS x NC) x (NS + NG + NH + NY) elements. For a large seismic survey this is rather a large matrix, though the data could be easily processed in patches. It is derived only from the geometry so may be preferable to reading all the data. For even my

- very largest data set the matrix is only 1421 by 231 elements. To find a solution we need at least the condition (l\rS + NG + N2ir + NY) c (NS x NC), or the matrix to have sufficient rank, or the smallest eigenvalue of the matrix A can be checked. The square matrix (ArA) is rank deficient because the fold of data drops at the begirming and end of the data set. This is equivalent to saying there are insufficient independent equations for the number of unknowns here, which results in zeros or near zeroes (to the computational tolerance) on the diagonal of the square matrix. Additional constraints can be provided on the value that elements of X ( xi) can take. Damping of the solution is achieved by adding equations Hi = 0, where is a small number. Equation then becomes -I]X = -O] Equation 7 Forming the least squares solution as before gives the damped least squares solution, X = (ATA +,52iJ ATD Equation 8 This is equal to adding a small value to the diagonal of the previous square matrix. I found that a value of = 0.01 gave satisfactory results. This is similar to adding white noise for stability in deconvolution. As this data contains noise, adding ú2 to the diagonal limits the resolution of the results to this small value and so this should be estimated to be small enough to give a stable solution.

The generalized matrix inverse Another candidate method for obtaining a solution is the Generalised Matrix Inversion method. The matrix equation y = Ax can be transformed by the substitutions y = UY, = VX, and so BY = AYX.

If U is orthogonal, Y = UtAV = Equation 9 where A is a diagonal matrix with the eigenvalues of A as the elements of the main diagonal. The factoring of A into U. V and is referred to as singular value decomposition. Now a suitable inverse operator H would be H = VA TUT Equation lO Since U and V are orthogonal square matrices HA = Y.-iU7U.AYT = I Equation 11 However as in the simple least squares case, problems arise when tile main diagonal elements of 1\ are zero or small so that 1cannot be formed. If the matrix of eigenvalues is reduced in size to a matrix p x p by dropping the singular eigenvalues then corresponding columns of U and V are also dropped to give rectangular matrices of p columns each. The generalised inverse operator then becomes H[ -'U T Equabon 12

where the suffix p denotes reduced to p non-singular eigenvalues. I used the estimated rank of as the estimate for p which gave very satisfactory results. Very similar results to the damped least squares solution were obtained, despite the extra complication and computation time. (e.g., Hatton et al, 1986).

For my largest test data set, a damped least squares method gives an acceptable solution in about 2 minutes on a 400Iz PC, with the singular value decomposition method it takes about twice as long.

Treatment of the phase spectra The previous sections are for the amplitude spectra. The usual approach to take for the phase spectra is to calculate it from the amplitude spectra on the assumption of minimum phase or to unwrap the phase for the calculations. The example data is recorded with a chirp which when correlated with the data should give a zero phase wavelet. Examples from the previous sections show that the wavelet is usually mixed phase. If the partitioned phases can be calculated directly from the data's phase spectra then no assumptions need be made. The phase given by the log spectra is the principal phase ( dip) i.e. = lip _ n Fir where 7 is an integer.

Taking the logarithm of the spectra might appear to corrupt the phase spectra, but a random integer number of or radians added to the phase spectra has no effect on the result of the inverse Fourier transform. Phase unwrapping is usually achieved by adding an integer number of r radians to the phases to smooth out the phase spectra, judged by some criteria such as minimising the steps in the result, e.g., Spagnolini (1993).

I take the phases, (the imaginary part of the log spectra), and transform them by the function em. This is equivalent to constructing a complex number coS O + i. sin. This fimction is now smooth and cyclic, period 2. Their magnitudes are always one, and so operations such as multiplication, mean and weighted means (the weights must sum to one) are valid operations on the phases. As the matrix has been derived for four components the matrix need only be scaled by 1/4 to produce the same phases, i.e. the same matrix inverse used to calculate the amplitude spectra is applied to the exponentiated phases ( eiR) and the argument of the result gives the phase spectra of the components. This is mathematically similar to a process used by McGlymn and Ioup (19S5) to calculate a phase coherency filter to enhance seismic data.

The complete spectra can then be assembled from the amplitude spectra (the real part) and the phase spectra (the imaginary part). The inverse Fourier spectrum then gives the impulse response or transfer function of each component.

An example data set Equation is not specific about the nature of the reflection data and could quite easily apply to the surface wave response as well, for example. Signals such as a high amplitude surface wave are best removed before this analysis so that they do not dominate desired data. Most of the data sets recorded with our self built equipment have very significant electrical crossfeed which would dominate the Fourier transform. Cross-

correlation with the sweep should remove it as it has zero delay, but the side-lobes of the sweep correlation are still significant. Spherical divergence gain compensation is best applied to make the wavelet statistically stationary, and therefore compensation for attenuation in the processing can be included here. For seismic data NMO would also be applied. NIvIO is derived from the approximation of a horizontally layered earth model, which obviously does not apply to ultrasonic data in concrete. An equivalent operation to NMO would be a prestack migration process for point reflector sources. Not applying the operator to the data will mean that it is part of the solution, appearing as part of the offset operator. Potentially, the data could be left uncorrelated with the sweep and this surface consistent deconvolution will extract the sweep as part of the source component (the transducer transfer function and coupling are also part of the component). However, this would preclude time domain, data-conditioning operations like gain compensation as the late part of the sweep will overlap the early part of the sweep from later reflections.

A data set was acquired with the Transeise equipment on a slab of the TRL sample concrete jOO mm by 500 trim by 150-165 mm (i.e. slightly wedge shaped). The full number OI available transmit/receive transducers were used. The coverage produced is displayed below. The missing zero offset traces are evident as the gap in the coverage.

Figure 1. A plot of the coverage of the data acquired on a slab of the sample concrete. Each dot marks one trace, horizontal lines are common shot gathers, vertical lines of dots are common receiver gathers. The ordinals of the common CMP gathers and common offset gathers are at 45 .

The data was correlated with the sweep and severe electrical cross-talk is evident on the near traces. The cross z-correlation should move the cross-talk to time zero, but the side lobes of the correlated sweep are still of high amplitude. Careful windowing of the data removes much of this. Deconvolution applied to the data whitens the spectra and also removes some more of this cross-talk.

Figure 2 A plot of one correlated shot gather of the input data. Severe electrical cross-talk is evident on the near traces.

Figure 3 A plot of one shot gather of the input data with windowing, scaling and deconvolution applied.

Results and conclusions

Fourier transformation of the input data traces gives the phases and log amplitudes shown in Figure 4. The frequency content of the 0.1-0.3 MHz sweep is evident. The amplitude of the transmitted sweep was reduced for the later part as it was thought the system input was overdriven. Better phase coherency can be seen in this data. A damped least squares solution was applied to these spectra to give the surface consistent spectra shown in Figure 5. The surface consistent components have been resolved very well. The weaker sources have been resolved with their better phase coherency. The frequency spectra of the receivers is similar to the sources as expected with a dead receiver resolved too. The offset spectra, show the expected higher amplitudes at near offsets with the spectra become noisier towards the edges as the fold of data drops, similar to the common midpoint spectra. The time domain equivalents of the four components are shown in Figure to Figure 9.

Figure 4 A plot of the input data amplitudes (dB) and phases. The horizontal scale is an arbitrary trace number. Note that the transmit power was reduced for the later traces and gives better phase coherency.

Figure 5 A plot of the output surface consistent component amplitudes (dB) and phases. The first 31 traces are the sources, the next 47 traces are the receivers, the next 76 are the offset components (from) and the last 77 are the common midpoint components.

The time domain equivalents of the source components shown in Figure 6 are hard to interpret. The amplitude spectra are as expected but the phases are very noisy. Ideally the source components would be a zero phase wavelet. Coda generated near the sources would be seen as a common source component. The phase would in turn be rapidly changing spatially, so that the phase has been aliased by the sampling. It could just be that the cross-talk has swamped the phases or resonances of the transducer with the couplant to the concrete. The weaker sources have been resolved.

Figure 6 A plot of the source components. Some of the transducers are dead, and the power was reduced for the last part as it was thought clipping may be occurring.

The time domain versions of the receiver components shown in Figure 7 are similarly difficult to interpret.

The regular bands of impulses across the traces are unexpected. They may be from the cross-talk or reflections from the bottom of the slab. From a similar experimental set up on the stone slab, shear waves at approximately half the velocity were observed.

Figure 7 A plot of the receiver components, one or two may be weak.

The time domain offset components shown in Figure 8 also have regular bands across the traces and has resolved the linear moveout of the surface waves as an offset component. Note that the null trace for the zero offset, as there is no zero offset data in the data set. Again, the far offset traces become noisy as the fold of data drops.

Figure 8 A plot of the offset components from -336mm to 272mm. There are no zero offset traces in the data set.

The the domain common midpoint traces shower in Figure 9 also have regular bands across the traces.

These could be the cross-talk again. But the corornon midpoint traces are equivalent to the zero offset stack obtained in seismic processing, even though NMO has not been applied. The signal at about 50,u; and again at 100 Mu could be the reflection and multiple from the bottom of the slab. The slab is slightly wedge shaped and the signal does appear to curve across the plot. The traces are noisier at the edges due to the drop in fold, as with the other components.

Figure 9 A plot of the common midpoint components at 4mm spacing.

Conclusion

The quality of the example data set does not coznpare very well with a typical exploration seismic data set but is sufficient to demonstrate the technique here. The separation of the data into four components with amplitudes and phases has worked very well. The new method of solution for the phase has been shown to work well, without assumptions (e.g. minimum phase) and using a similar method to that used for the amplitudes. The phase solution has been found despite the mixed phase waveless of the source. The biggest problem appears to be the quality of the sample data. Better pre-conditioning of the data to remove noise and known effects like diffraction, should help. Also interpolation of the data to a finer trace interval to remove phase aliasing effects (equivalent to spatial aliasing), would be beneficial. A seismic data set will probably benefit from suitable scaling and NMO and other conditioning which might be applied pre-stack.

Even pre-stack migration could be applied before this process demonstrated here.

The large matrix to be solved might appear to be a problem, but it is only built from the geometry of the data. Also, once the matrix is solved, it is used for the phase solution as well as the amplitude solution.

Experiments with a normal seismic data set should be performed.

(alit References Cambois, G. and Stoffa P., 1992, Surface-consistent deconYolution in the log/Fourier domain, Geophysics, 57, now, 823-840.

Cambois, G. and Stoma P., 1993, Surface-consistent phase decomposition in the log/Fourier domain, Geophysics, 58, no.8, 1099-1111.

Cary P. W. and Lorentz G.A., 1993, Surface-Consistent deconvolution in the log/Fourier domain, Geophysics, 58, no.8, 383-392.

Hatton, L., Worthington, M. H., Makin, J., 1986, Seismic data processing, theory and practice, Blacksvell Science Ltd. McGlynn, J. D. and soup, G. E:, 1985, Phase coherency filtering of reflection seismic data (short note): Geophysics, 0, no. 09, 1505-1509.

Morley, L. and Claerbout, J., 1983, Predictive deconvolution in shotreceiver space, Geophysics, 48, 515-531.

Poggagliolmi, E., Berkhout, A. J., Boone, M. M., 1982, Phase unwrapping, possibilities and limitations, Geoph. Prosp., 30, 281-291.

Spagnolini, U., 1993, 2-D Phase unwrapping and phase aliasing, Geophysics, 58, now, 1324-1334.

Taner, M.T., Lee, L., Baysal, E., 1974, Estimation and correction of nearsurface time anomalies, Geophysics, 39, 441-463.

_.

Claims (3)

1. A method or process is described Hat will find a solution to a set of linear equations, where He data are phasors, complex numbers or numbers representing angles as well as magnitude.

2. A method or process is described that will find the solution to the phase spectra as well asithe amplitude spectra, for a composite set of input data, where a solution exists.

3. A method or process is described Hat will deconvolve a set of signals into their component parts, where He signals share component parts g