2D phase unwrapping by DCT method Kui Zhang, Marcilio Castro de Matos, and Kurt J. Marfurt

1. A(1,1) B(1,1) 2D phase unwrapping by DCT method Kui Zhang, Marcilio Castro de Matos, and Kurt J. Marfurt ConocoPhillips School of Geology & Geophysics, University of Oklahoma Phase dispersion correction Seismic waves propagating through the subsurface undergo strong energy dissipation and velocity dispersion due to both anelasticity and heterogeneity within the earth. High frequency data components suffer more loss than low frequency components that traveled along the same ray path, resulting in a relatively narrow-band, low-frequency spectrum. In general, the high frequencies travel at different velocities than the low frequencies resulting in a significant change in waveform shape. Due to the above two effects, the seismic wavelet become noticeably stretched, and displays “ringing” characteristic as the travel-time increases. The traditional approach to compensate seismic wave attenuation and dispersion is inverse Q filtering by assuming the earth Q model to be a multilayedstructure. Inverse Q filtering is implemented in a layer stripping manner(Wang, 2002) by downward continuation. However, very few papers have mentioned about how to estimate Q values of the shallow earth effectively. Noted exceptions are work by Singleton and Taner (2006) who use logs and Chopra and Alexeev (2003) who use VSPs. In our study, we provide a different approach to eliminate velocity dispersion. Without considering attenuation, a plane wave at time t that has undergone dispersion can be formulated as: 2πft is the linear phase part caused by the time delay, φQ(f,t) is a nonlinear phase delay that depends on both intrinsic and geometric Q dispersion, and φG(f,t) is directly related to the underlying geology and impedance of the layers, including π phase changes associated with negative reflection coefficients,  π/2 phase changes associated with thin bed tuning and upward fining/coarsening, as well as more complicated phase changes due to stratigraphic layering. Our objective is to retain the linear part and φG(f,t), and remove the effects of φQ(f,t) to compensate for the dispersion to obtain φ(f,t): The phase dispersion corrected signal in time domain can be obtained by summing all plane waves. (c)The unwrapped result from (b) (using B(1,1) as the reference) (d)The unwrapped result from (b) (using A(1,1) as the reference) (e) 1D unwrapped result from (b) (notice its big difference with (c)) Figure 1. (a) 2D Gaussian model (b)The wrapped phases of (a) Real data example 2D phase unwrapping by Discrete Cosine Transform Ghiglia (1994; 1998) gave several discrete approaches to unwrap the phase by solving Poisson’s equation: The least squares approach requires Neumann boundary condition on Poisson’s equation, that is: Applying the two dimensional Fourier transform to the two sides of (8), we get Where Ψand P are the two-dimensional Fourier transform of ϕ and ρ respectively. The solution ϕcan then be obtained by inverse Fourier transform. To satisfy the Neumann boundary condition (11), we need either mirror the 2D Laplacianρi,j, and then apply the fast Fourier transform, or apply the two dimensional forward DCT (Discrete Cosine Transform) to ρi,j thereby eliminate the need for mirroring. In our study, we implemented the algorithm using the DCT method. To test the algorithm, we designed a 2D Gaussian model, wrapped it, and then unwrapped the phase and compared with the original model. The least squares solution described above is relative such that one reference point needs to be set by the user to calibrate it. Figure 1d shows that if we know only one true absolute phase of the model and use it as the reference, the unwrapped result is exactly the same as the original model (Figure 1a). Why phase unwrapping? In order to estimate φQ(f,t), we average phase over a finite time window of thickness T over the entire seismic survey so that equation 2 becomes: Where, the sign < > denotes the average operator. Assuming that the reflectivity character of the earth is white, Combining (4), (5), and assuming the invariablity of Q model horizontally, we get In order to calculate the average of equation (6), all the phases ψ(f,t) must be unwrapped . (b) Unwrapped result from 2(a) Figure 2. (a) Wrapped phase References Chopra, S., A., Alexeev, V. Sudhakar, 2003, High frequency restoration of surface seismic data: The leading Edge, 22, 730-738. Ghiglia, D. C., and M. D. Pritt, 1998, Two-dimensional phase unwrapping: Theory, algorithm, and software: John Wiley & Sons, Inc. Itoh, K., 1982, Analysis of the phase unwrapping problem: Applied Optics, 21, p2740. Lomask, J., 2007, Seismic volumetric flattening and segmentation, Ph.D. thesis, Stanford University. Singleton,. S., M. T., Taner, 2006, Q estimation using Gabor-Morlet joint time-frequency analysis techniques: 76th Annual meeting SEG, Expanded abstracts, 1442-1446. Wang, Y., 2002, A stable and efficient approach to inverse Q-filtering: Geophysics, 67, 657-663. Itoh’s 1D phase unwrapping In 1982 Itoh proposed that the unwrapped phase can be obtained by integrating wrapped phase differences (Itoh, 1982). The unwrapped phase will equal the true phase provided the true phase differences are less than πradians in magnitude everywhere. If W is the wrapped operator, Φis the true phase, and φis the wrapped phase (φ (n)=W{Φ(n)}, n=0,… N-1). The 1D phase unwrapping can be realized by: Since the spectral phase of a seismic trace after spectral decomposition will be a 2D panel, if we apply 1D phase unwrapping to every component, the output will have some vertical stripes because spectral phase also change laterally. Considering this fact, 2D or multi-dimensional phase unwrapping is essential. Conclusions • The statistical phase compensation after 2D phase unwrapping will be investigated using real seismic data. • The DCT method of 2D phase unwrapping provides a robust result, but it still needs to be tested by more models and real data examples. • Some attributes based on the unwrapped phase will be analyzed to map structural and stratigraphic information. Acknowledgement We thank all the sponsors of Attribute Assisted Seismic Processing and Interpretation Consortium for their support. We also thank Tim Kwiatkowski, Jesse Lomask for their helpful discussions. 11/18/2008