Pure Appl. Geophys. Ó 2016 Springer International Publishing DOI 10.1007/s00024-016-1238-7

Pure and Applied Geophysics

Source Wavelet Phase Extraction DIAKO HARIRI NAGHADEH1 and CHRISTOPHER KEITH MORLEY1 Abstract—Extraction of propagation wavelet phase from seismic data can be conducted using first, second, third and fourthorder statistics. Three new methods are introduced, which are: (1) Combination of different moments, (2) Windowed continuous wavelet transform and (3) Maximum correlation with cosine function. To compare different methods synthetic data with and without noise were chosen. Results show that first, second and third order statistics are not able to preserve wavelet phase. Kurtosis can preserve propagation wavelet phase but signal-to-noise ratio can affect the extracted phase using this method. So for data set with low signal-to-noise ratio, it will be unstable. Using a combination of different moments to extract the phase is more robust than applying kurtosis. The improvement occurs because zero phase wavelets with reverse polarities have equal maximum kurtosis values hence the correct wavelet polarity cannot be identified. Zero-phase wavelets with reverse polarities have minimum and maximum values for a combination of different-moments method. These properties enable the technique to handle a finite data segment and to choose the correct wavelet polarity. Also, the existence of different moments can decrease sensitivity to outliers. A windowed continuous wavelet transform is more sensitive to signal-tonoise ratio than the combination of different-moments method, also if the scale for the wavelet is incorrect it will encounter with more problems to extract phase. When the effects of frequency bandwidth, signal-to-noise ratio and analyzing window length are considered, the results of extracting phase information from data without and with noise demonstrate that combination of differentmoments is superior to the other methods introduced here. Key words: Propagation wavelet phase, combination of different moments, Kurtosis, Zero phasing.

1. Introduction Seismic data is one of the main source information for estimating reservoir properties in hydrocarbon reservoirs. Since information obtained from well logs is limited to a narrow zone around the

1 Petroleum Geophysics, Department of Geological Sciences, Chiang Mai University, Chiang Mai, Thailand. E-mail: [email protected]

well path, seismic data is used to project properties between wells. While well log data alone is used to predict rock properties between wells using geo-statistical methods, the resulting estimates can have a high degree of uncertainty depending upon geological conditions and data quality. In seismic sections amplitude variations occur due to impedance contrasts at layer boundaries. Acoustic impedance, which is related to subsurface rock properties are obtained by seismic inversion. Since acoustic impedance is the product of density and velocity, wells with logs of these properties can be used to generate impedance logs. To obtain an acoustic impedance section from a seismic cube a wavelet should be estimated from either well and seismic data or just seismic data. Inaccurate estimation of the wavelet can cause inaccurate reflection coefficients and acoustic impedance estimations. It is possible to extract the amplitude spectrum of the wavelet from the seismic data but estimation of the best phase property is not easy. Every mistake in phase estimation can cause a major error in reflection coefficient and acoustic impedance estimations. Different techniques for wavelet extraction and phase correction of seismic data have been proposed. These include: 1. Automatic phase correction of common-midpoint stacked data (LEVY and OLDENBURG 1987) 2. Maximum kurtosis phase correction (WHITE 1988) 3. Time-varying wavelet estimation and deconvolution by kurtosis maximization (VAN DER BAAN 2008) 4. Nonstationary phase estimation using regularized local kurtosis maximization (VAN DER BAAN and FOMEL 2009) 5. Short-time wavelet estimation in the homomorphic domain (HERRERA and VAN DER BAAN 2012)

D. H. Naghadeh, C. K. Morley

6. Local skewness attribute as a seismic phase detector (FOMEL and VAN DER BAAN 2014) Real seismic signals are not zero phase and they have to be rendered zero-phased to be suitable for interpretation and automatic tracking on workstations (MANSAR and RODRIGUEZ 1996). Perturbation of wavelet phase can be created by acquisition and processing, consequently it is necessary to process the data for zero phasing and to correct residual phase distortion (LEVY and OLDENBURG 1987). To do zero-phasing statistical approaches are commonly applied to find the best phase shift that will generate zero phased data for interpretation. These approaches need many considerations to obtain the best results. To determine local and propagation wavelet phase, a finite data segment and long data segments should be chosen, respectively. The main target of this paper is to compare different statistical approaches to find the best one to extract propagation wavelet phase. We introduce three methods which are combination of different moments (CDM), windowed continuous wavelet transform (WCWT) and maximum correlation with cosine function (MCWC). CDM is a technique that measures the sharpness of the probability distribution of a discrete time series. It can overcome the limitations of kurtosis especially for wavelet polarity detection and local seismic phase correction using short time windows. WCWT is a windowed continuous wavelet transform. It works based on obtaining the maximum wavelet coefficient for special phase shifted data set in each window. MCWC is based on calculating the correlation between the Cosine function and constant phase rotated signal. The maximum correlation coefficient shows the best phase rotation for attaining zero phasing. To compare these methods a synthetic data set is used. 1.1. First-Order Statistics (FOS) FOS contains mean, median and mode. Mean of a time series xt returns the mean value of time sample amplitudes of xt (DAVID 1970, 1981). The quantity of mean is defined as

Pure Appl. Geophys.

PN Meanðxt Þ ¼

t¼1 xt

N

ð1Þ

where N is the total number of the samples. To calculate the median of a series, each number is ordered by size so the number in the middle is the median (DAVID 1970, 1981). The mode is the number that occurs most often within a set of numbers (DAVID 1970, 1981). In Fig. 1b values of FOS for rotated Ricker wavelets from -90° to 90° (Fig. 1a) are shown. It is clear that the mean value for the zero phase wavelet is the largest one but it is very close to zero and it shows small random number. For noisy data which includes many wavelets with different polarities, the mean value goes to zero under this condition it is unsuitable to render zero phasing. The median value for shifted wavelets is unable to preserve phase because the value for zero phase wavelet is neither the largest nor smallest one. The mode value for the zero phase wavelet is the largest one. However, using mode as a FOS for phase extraction from seismic data is problematic because its value for a long window, that contains many wavelets, is near zero so it cannot preserve the phase of the propagation wavelet. To confirm the inability of FOS to preserve phase, a synthetic reflection coefficient and a synthetic seismogram (Fig. 2b) created by convolution of reflection coefficient (Fig. 2a) with a -60degree phase shifted Ricker wavelet, were chosen. In Fig. 3a, b and c mean, median and mode values for a seismogram rotated from -90° to 90° are shown, respectively. To rotate the phase of signal Eq. 2 is used. xrot ðtÞ ¼ cosðhÞ:xðtÞ þ sinðhÞ:H ½xðtÞ

ð2Þ

where, xrot is the rotated data with the phase rotation angle h and His the Hilbert transform. The results show that the FOS is not able to preserve the phase of the wavelet despite the maximum values of mode and mean corresponding to a zero phase wavelet. Seismic recording systems cannot record DC (zero frequency) components, hence their-FOS for a long data segment, whether based on the mean, median or mode, are all inherently zero and cannot yield information on wavelet phase.

Source Wavelet Phase Extraction

Figure 1 a Rotated Ricker wavelets from -90° to 90°; b values of mean, median and mode for rotated Ricker wavelets

D. H. Naghadeh, C. K. Morley

Pure Appl. Geophys.

Figure 2 a Synthetic reflection coefficient; b synthetic seismogram created by convolution of reflection coefficient (a) with a -60° phase shifted Ricker wavelet

1.2. Second-Order Statistics (SOS)

1.3. Third-Order Statistics (TOS)

Variance is a SOS and is a measure of dispersion. It is equal to the second moment of a time series around the mean (ROBERTS and RICCARDO 1999). The zero value of variance indicates that all the values of time series samples are identical and are equal to the mean. A small value of variance shows that the data points are very close to each other and to the mean value. A high value of variance indicates that the data points are distributed far from the mean and from each other (ROBERTS and RICCARDO 1999). PN ðxt  MÞ2 ð3Þ Varinaceðxt Þ ¼ t¼1 N

The target of this section is to discuss the use of skewness as a TOS to extract phase of wavelet. Skewness is discussed in terms of the second and third moments around the mean (DOANE and SEWARD 2011). P 1=N Nt¼1 ðxt  M Þ3 ð4Þ Skewnessðxt Þ ¼  3=2 1= PN ðxt  M Þ2 N t¼1

where xt is an input series, M is mean of xt and N is the total number of elements of xt. In Fig. 4 the values of variance for rotated wavelets (Fig. 1a) are shown. Since changing the phase of the data does not affect the variance value it is concluded that SOS is not a useful tool for phase correction.

where, xt, M and N are amplitude of seismic samples, mean of xt and the number of samples, respectively. The skewness value represents the asymmetry of the distribution of a variable (DOANE and SEWARD 2011). The skew value of a Gaussian distribution is zero. Positive and negative skew values show that the tail on the right side of the distribution is longer than the left side, and the bulk of the values lie to the right of the mean, respectively (KIM 2013). In Fig. 5 skewness values for rotated wavelets (Fig. 1a) are shown.

Source Wavelet Phase Extraction

Figure 3 a Mean values for rotated synthetic seismogram (Fig. 2b); b median values for rotated synthetic seismogram (Fig. 2b); c mode values for a rotated synthetic seismogram (Fig. 2b)

The skew value for a zero phase wavelet is maximum. The absolute values of skewness for zero phase wavelets with reverse polarities are equal. So the skew value for a signal with an equal number of wavelets with reverse polarity could be zero. Commonly for long data segments the use of skewness cannot preserve the phase of wavelet. The skew values for a phase rotated synthetic seismogram are shown in Fig. 6. The maximum value of skewness is for a 90° phase shifted synthetic trace. It is clear that TOS is not suitable for the wavelet phase extraction not only for noise free but also for noisy dataset.

1.4. Fourth (High)-Order Statistics (HOS) In this section kurtosis is introduced as a HOS. The kurtosis value of a variable equals to its fourth moment around the mean divided by the square of its variance (WHITE 1988) and it is calculated by Eq. 5.

Kurtosisðxt Þ ¼ 

PN

ðxt  M Þ4 2 1= PN ðxt  M Þ2 N t¼1 1=N

t¼1

ð5Þ

where, xt, M and N are amplitude of seismic samples, mean of xt and the number of samples, respectively. The kurtosis value should be three for a Gaussian distribution. HOS values above and below a value of three have a leptokurtic and platykurtic distributions, respectively (KIM 2013). Wavelet phase extraction using HOS is a two steps procedure. In the first step, a series of constant phase rotations are applied to the seismic data using Eq. 2. In the final step, the angle corresponding to the amount of rotated phase that yields the highest kurtosis value of the data set with the minus sign is chosen as the wavelet phase. In Fig. 7 kurtosis values for rotated Ricker wavelets from -180° to 180° are shown. It is clear that zerophase wavelets with reverse polarities have equal maximum kurtosis value. Also, it can be concluded

D. H. Naghadeh, C. K. Morley

Pure Appl. Geophys.

Figure 4 Variance value for phase rotated wavelets (Fig. 1a)

that external control is required to determine wavelet polarity using kurtosis. Kurtosis due to the higher powers is sensitive to outliers. Also fourth-order statistics are very difficult to estimate for short data sequences, but that preserve phase information for long data segments. The kurtosis values for a phase shifted seismogram (Fig. 2b) are shown in Fig. 8. The HOS is able to extract the correct phase of wavelet from long data segments with a high signalto-noise ratio (S/N). In Fig. 8, the amount of phase correction to yield maximum kurtosis is 60°, consequently the wavelet phase should be -60°. 1.5. Combination of Different Moments (CDM) We looked for a statistical approach to overcome kurtosis limitations (inability to detect local phase using short windows, inability to detect the correct wavelet polarity and it is sensitive to outliers). We found L-kurtosis (HOSKING 1990) which is combination of L-moments or probability weighted moments and calculated by Eq. 6:

L Kurtosis ¼

20b3  30b2 þ 12b1  b0 2b1  b0

ð6Þ

where, bi are probability weighted moments and defined by Z bi ¼ x½FðxÞi dFðxÞ; i ¼ 0; 1; 2; . . . ð7Þ where, F(x) is a cumulative distribution function of a probability distribution x. In Fig. 9, L-kurtosis values for rotated Ricker wavelets (from -180° to 180°) are shown. When L-kurtosis values for zero phase wavelets with diverse polarities are same, then it is not able to detect the correct polarity of wavelets (as same as kurtosis). It is essential to mention that L-kurtosis is less sensitive to outliers than kurtosis because of combination of L-moments, also it doesn’t have any problem with short window analysis to extract local phase. Finally we decided to use moments instead of L-moments to create the CDM method with special properties (ability of correct polarity detection, ability

Source Wavelet Phase Extraction

Figure 5 Skew value for rotated wavelets (Fig. 1a)

Figure 6 Skew value for phase rotated synthetic seismogram (Fig. 2b)

D. H. Naghadeh, C. K. Morley

Pure Appl. Geophys.

Figure 7 Kurtosis value for rotated Ricker wavelets from -180° to 180°

of zero phase rendering not only for long windows but also for short windows, less sensitivity to noise because of combination of different moments), to overcome all the limitations of kurtosis and L-kurtosis. CDM is a technique that measures the sharpness of the probability distribution of a discrete time series. The CDM of a sequence x is defined as

CDMðxÞ ¼

Second-Moment = 0.9842, Third-Moment(xt) = -5.99e - 17, Fourth-Moment(xt) = 1.467 and CDM(xt) = 2.491. The value of different moments for a window with 10 samples, are First-Moment(xt) = Mean(xt) = 0, Second Moment(xt) = 0.2, Third-Moment(xt) = -0.054, Fourth-Moment(xt) = 0.088 and CDM(xt) = 1.702.

momentðx; 4Þ  momentðx; 3Þ þ momentðx; 2Þ  momentðx; 1Þ momentðx; 2Þ  momentðx; 1Þ

To show variability between different moments, a harmonic wave is chosen with a 5 (Hz) predominant frequency [xt = 1.41 9 sin (10pt), t = 0:0.01:1(s)]. The value of different moments for a window with 101 samples, are First-Moment(xt) = Mean(xt) = 0,

ð8Þ

The values of different moments show that for a harmonic wave (with a long data segment) the first and third moments are zero. The most important parameter used to create a non-zero CDM value for a long data segment is the fourth moment. It is able to

Source Wavelet Phase Extraction

Figure 8 Kurtosis values for a rotated synthetic seismogram (Fig. 2b) from -90° to 90°

extract the propagating wavelet phase using long seismic data windows, which have near zero mean values. To conduct local phase correction, short window lengths are required. Different-moment values for the portion of harmonic wave with ten samples show that the third-moment and fourthmoment play important roles determining CDM value. Hence a combination of different moments can help not only to conduct phase correction but also for polarity detection. In Fig. 10 CDM values for rotated Ricker wavelets are shown. CDM values for zero phase wavelets with diverse polarities are different hence zero phase wavelets with the correct polarity can be detected. To extract a wavelet phase using CDM, the procedure contains two steps. First step is: to apply a series of constant phase rotations to the signal using Eq. 2. The second one is: to determine the angles corresponding to the amount of rotated phase to yield the highest and smallest CDM values. The absolute values of those angles are compared and the smallest one is chosen as the essential rotation to create zero phase signals. In

Fig. 11 CDM values for a phase shifted seismogram (Fig. 2b) are shown. The phase information is preserved correctly and the wavelet phase is -60°. 1.6. Windowed-Continuous-Wavelet-Transform (WCWT) The Continuous wavelet transform (CWT) provides a different approach to time–frequency analysis. Instead of producing a time–frequency spectrum, it produces a time-scale map called a scalogram (RIOUL and VETTERLI 1991). The CWT of a signal x(t) with the wavelet w is defined as Xw ðs; TÞ ¼  xðtÞ; ws;T ðtÞ    Zþ1 1 T t xðtÞ: pffiffi w ¼ dt s s

ð9Þ

1

where, s is scale value, T is translation value, and Xw(s, T)is the time scale map (SINHA et al. 2005). CWT also represents the correlation of the input signal with a time-reversed version of w rescaled by a

D. H. Naghadeh, C. K. Morley

Pure Appl. Geophys.

Figure 9 L-Kurtosis values for rotated Ricker wavelets from -180° to 180°

factor of s (MUNOZ et al. 2002). For a 1-D input signal, the result is a 2-D description of the signal with respect to time T and scales. The scale s is inversely proportional to thecentral frequency of the  rescaled wavelet ws ðtÞ ¼ w st which is typically a band pass function, where t represents the time location at which we analyze the signal. The larger the scales, the wider the analyzing function ws, and hence smaller the corresponding analyzed frequency. The output value is maximized when the frequency of the signal matches that of the corresponding dilated wavelet (MUNOZ et al. 2002). At this stage the target is using WCWT to obtain maximum correlation for each window. The procedure contains four steps. The first step is to choose the best scale for the Mexican hat wavelet (DAUBECHIES 1992; MALLAT 1989) which is used as a mother wavelet based on the predominant frequency of signal. The second step is windowing

the signal, it is essential to choose windows which include at least 200 samples. Third is constant phase shifting of the signal and calculation of the wavelet coefficient for each window. The last step is choosing the maximum Xw(s, T) for each window based on the special phase shifted data in that window. The angle corresponding to the amount of rotated phase to yield the maximum Xw(s, T) is the amount of angle to obtain zero phase wavelets. It is possible to use this approach for local phase correction as well. To clarify the reason of choosing Mexican hat wavelet to detect the phase of signal, many wavelets (Morlet, complex Gaussian4 and Daubechies2) are compared. In Fig. 12a, b, c and d the extracted phase of synthetic seismogram (Fig. 2b) for a window equals to the whole length of trace using Mexican hat, Morlet, complex Gaussian4 and Daubechies2 are shown, respectively. The results show that the extracted

Source Wavelet Phase Extraction

Figure 10 CDM values for rotated Ricker wavelets from -180° to 180°

phase of wavelet using Mexican hat is more precise than the others because of short tail which has it. In Fig. 13 phase corrected angles for different windows of synthetic seismogram (Fig. 2b) using Mexican hat as a mother wavelet, are shown. The estimated phase in Middle part of Fig. 13 is not correct because that window includes 100 samples and it is lower than essential number of samples to carry it out correctly, then using WCWT for small windows is not suitable. 1.7. Maximum Correlation with the Cosine Function (MCWC) In this part MCWC is introduced as an approach to extract wavelet phase. To conduct MCWC correlation between f ðtÞ ¼ A2 ðexpðiwtÞ þ expðiwtÞÞ and signal is calculated. Where A is max amplitude of signal, w is angular frequency (w = 2pf and f is predominant frequency of signal) and t is time samples. Correlation coefficients between two signals are calculated using Eq. 10.

Correlationðf ; xÞ N P P P N Ni¼1 fi xi  xi Ni¼1 fi i¼1 ffi ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN 2 PN PN 2 PN 2 x f  f N i¼1 xi  N i i i¼1 i¼1 i i¼1

ð10Þ where N is number of samples, f and x can be same or different signals. To extract correct propagation phase angle, a series of constant phase rotations are applied to the signal and for each rotation the correlation between f(t) and signal is calculated. Maximum correlation coefficient shows the best phase rotation to attain zero phasing. To extract the phase of the wavelet correctly the dominant frequency of signal is chosen as a dominant frequency of f(t) and the total length of trace should be chosen as a length of correlation window. In Fig. 14 correlation values for the phase shifted synthetic seismogram (Fig. 2b) from -90° to 90° are shown. Maximum correlation is

D. H. Naghadeh, C. K. Morley

Pure Appl. Geophys.

Figure 11 CDM values for phase shifted synthetic seismogram (Fig. 2b)

obtained by a 60° phase shift so the wavelet phase is -60°. 1.8. Effect of Noise on Extracted Phase of Wavelet Extracting propagation wavelet phase from noisy data is a big challenge using statistical methods. HOS are very sensitive to S/N ratio, so they will be unstable in low S/N situations. On the other hand, due to the higher powers, HOS are significantly less robust than lower-order statistics. CDM is a combination of different moments hence it should be more stable than kurtosis as a HOS. To compare phase extraction using HOS, CDM, WCWT and MCWC a synthetic seismogram (Fig. 2b) with different level of added random noise, is chosen. In Fig. 15 the amount of phase corrections for S/N ratios from 1 to 10 using different methods are shown. The results of phase correction show that MCWC and CDM (Fig. 15b, c) are just affected at S/N = 1, and for other values they can extract the correct phase angles. Phase corrected

angles using WCWT are incorrect for S/N ratios ranging from 4 to 1 (Fig. 15d). Consequently for short windows that are chosen to extract phase in regions with low S/N, the wrong results can be generated. Using kurtosis for phase extraction from noisy data encounters problems where the S/N ratios range from 5 to 1 (Fig. 15e). Based on the examples above, it is clear that MCWC and CDM are methods that can extract correct wavelet phase not only for data sets free of noise, but also for noisy data sets. To choose the best approach the effect of frequency bandwidth on two last methods needs to be clarified. 1.9. Effect of Frequency Bandwidth on Extracted Phase of Wavelet Using MCWC and CDM In this part, true reflectivity from well E_01 which is located at inline 892-823 of Songkhla basin in Thailand, is chosen to test phase extraction using MCWC and CDM in different frequency bandwidths. True reflectivity (Fig. 16a) is convolved with a -30°

Source Wavelet Phase Extraction

Figure 12 The extracted phase of synthetic seismogram (Fig. 2b) for a window equals to the whole length of trace using: a Mexican hat; b Morlet; c complex Gaussian4; d Daubechies2

phase shifted Ricker wavelet to create a trace which is shown in Fig. 16b. Then the trace is filtered using a trapezoidal filter with corner frequencies (f1, f2, f3, f4). The results of phase extraction from noise free trace using these methods are shown in Tables 1 and 2. The low and high frequency cutoffs do not play an important role in extracting the correct phase from noise free data using MCWC. But the amount of low frequency cutoff can affect CDM results especially when a high cutoff exceeding 50 hz is chosen. In this situation where the low frequency cutoff is high, the signal spectrum needs to be extended toward lower frequencies. To compare their results on noisy data, random noise is added to trace (Fig. 16b) which S/ N is equal to 2. The results of phase extraction from noisy trace are given in Tables 3 and 4 for MCWC and CDM, respectively. The MCWC phase extracted results from filtered noisy data are range between -18° and -25°, and have maximum and minimum errors of 40 and 16.6 %, respectively. CDM is more stable, with maximum and minimum errors of 10 and 3 %.

1.10. Effect of Window Length on Extracted Phase of Wavelet Using MCWC and CDM The length of window that is used to extract phase of wavelet plays an important role for all of statistical methods. As we know, long data segments are used to extract the phase of propagation wavelet. In regions with low S/N ratio the length of window should be increased to decrease the effect of noise. To compare the effect of this parameter on extracted results using both methods the real seismic data set from New Zealand entitled ‘Bo_Parihaka, inline 81SY-11’ was chosen which is located in the Taranaki basin. This basin is the only producing basin in New Zealand with commercial reserves of oil and gas, where it has many units of sedimentation and quite complex structures. To apply MCWC on a data set the procedure contains four steps as follows: (1) application of a series of constant phase rotations, (2) obtain the correlation between the cosine function and each trace for each rotation angle, (3) calculate the mean value for the obtained absolute correlation

D. H. Naghadeh, C. K. Morley

Pure Appl. Geophys.

Figure 13 Phase corrected angles for different windows of synthetic seismogram using WCWT

coefficients for each rotation angle, and (4) choose the maximum correlation coefficient that shows the best phase rotation to obtain zero phasing. In Table 5 the extracted phase angles for different windows are shown. Based on the results MCWC can extract the propagation wavelet phase correctly when the window length is equal to total length of record, but it is not suitable for use as a tool to extract the time variant wavelet phase angle. Conversely CDM can be used both for time invariant phase extraction and for time variant wavelet phase extraction. Since the amplitude spectrum of the wavelet can be obtained from the amplitude spectrum of the seismic signal, the extracted phase information can be used to construct a wavelet to conduct seismic inversion and determine the correct reflection coefficient. In Fig. 17a and c real seismic data and extracted wavelets using CDM, from different windows of data set are shown, respectively. The extracted

wavelets from different windows were used in sparse-spike deconvolution (ROBINSON and TREITEL 1980) process to obtain reflection coefficient from seismic data. Sparse-spike deconvolution assumes that the reflectivity is a sparse sequence of spikes. The main target of the Sparse-spike deconvolution method is to provide a significant increase in bandwidth content from band-limited seismic data (VELIS 2008). To obtain reflection coefficient from a data set, the source wavelet effects must be removed from the data using deconvolution. In Fig. 17b the reflection coefficient that obtained from a real data set was shown. 1.11. Lateral Variability of Estimated Phase Angle To obtain zero phasing of a data set both temporal and lateral phase correction should be considered. Zero phasing was conducted on Bo_Parihaka, inline

Source Wavelet Phase Extraction

Figure 14 Correlation coefficients for a rotated seismogram (Fig. 2b) from -90° to 90°

81SY-11. To correct phase distortion the analyzing window included 40 traces and 500 samples. For each window the optimum phase angle to do phase correction using CDM was estimated. The overlap of two adjacent windows is 39 traces. In Fig. 18a and b the phase corrected seismic section and phase angles used for phase correction for different windows are shown, respectively. To estimate the correct phase angle it is essential to choose at least 40 traces and the number of samples for each trace is equal to 500. To obtain zero phase data set, dealing with phase angle extraction in vertical and lateral directions is essential. The zero phase rendered data set is great to pick horizons because of high continuity of reflectors that is prepared by removing lateral phase distortion. In Fig. 17b which is the result of deconvolution, data set was phase corrected just in vertical direction and it is assumed that phase of wavelet in horizontal direction is constant without any changing. In reality

in both directions the phase of wavelet because of frequency absorption, is changed.

2. Conclusion To extract the propagation wavelet phase from seismic data statistical approaches are used. Choosing the FOS as a tool for zero phasing is not practical, because seismic recording systems cannot record zero frequency components. Hence the FOS for long data segments, whether based on the mean, median or mode, all are inherently zero and cannot yield information on wavelet phase. On the other hand changing the phase of the data does not affect the variance value, so it is concluded that SOS is not a tool for phase correction. The absolute values of skewness for zero phase wavelets with reverse polarities are equal. So the skew value for a signal with an equal number of

D. H. Naghadeh, C. K. Morley

Pure Appl. Geophys.

Figure 15 a Synthetic seismogram (Fig. 2b) where the S/N ratio is equal to 1; b phase correction angles for different S/N ratios using MCWC; c phase correction angles for different S/N ratios using CDM; d phase correction angles for different S/N ratios using WCWT; e phase correction angles for different S/N ratios using Kurtosis

wavelets with reverse polarities could be zero. Then commonly long data segments that use skewness will not preserve the phase of the wavelet. Kurtosis, due to use of higher powers, is sensitive to outliers. Additionally as a fourth-order statistic, it is very difficult to estimate phase for short data sequences. It cannot determine the correct wavelet polarities, but it preserves phase information for long data segments. Kurtosis is unable to extract the correct phase where the S/N ratio ranges between 5 and 1. WCWT coefficients are a correlation between wavelet and windowed data. The angle corresponding to the amount of rotated phase to yield the maximum CWT coefficient is the amount of angle required to create zero phase wavelets. By choosing one constant scale the method can extract the correct phase for some windows and the incorrect phase for others.

Also low S/N ratio affects many extracted phase angles using this method. Correlation between signal and cosine function was introduced as a method to extract correct phase not only for data without noise but also for noisy data. This method can be used just for extracting propagation wavelet phase using long data segments. To carry out the method correctly, short windows are avoided, and longer windows are used. These long windows are equal to the total length of the record, and are used to obtain correct maximum correlation. The extracted phase from noisy data without frequency filtering is the correct phase for different S/ N ratios, except where S/N = 1. CDM can overcome the limitations of kurtosis especially for wavelet polarity detection and local seismic phase correction using short time windows. The improvement occurs because zero-phase

Source Wavelet Phase Extraction

Figure 16 a True reflectivity; b trace which is created by convolution of (a) with -30° phase shifted Ricker wavelet

D. H. Naghadeh, C. K. Morley

Pure Appl. Geophys.

Table 1 Effect of frequency bandwidth on extracted phase of wavelet from free noise data using MCWC ðf3 ;f 4 Þ ðf1 ;f2 Þ ðHzÞ

30–35

40–45

50–55

60–65

1–5 5–10 10–15

-30 -30 -30

-30 -30 -30

-30 -30 -30

-30 -30 -30

Table 2 Effect of frequency bandwidth on extracted phase of wavelet from free noise data using CDM ðf3 ;f4 Þ ðf1 ;f2 Þ ðHzÞ

30–35

40–45

50–55

60–65

1–5 5–10 10–15

-30 -30 -30

-30 -30 -30

-30 -30 -28

-30 -30 -27

Table 3 Effect of frequency bandwidth on extracted phase of wavelet from noisy data using MCWC ðf3 ;f4 Þ ðf1 ;f2 Þ ðHzÞ

30–35

40–45

50–55

60–65

1–5 5–10 10–15

-23 -21 -24

-25 -21 -21

-18 -22 -23

-25 -22 -22

Table 4 Effect of frequency bandwidth on extracted phase of wavelet from noisy data using CDM ðf3 ;f4 Þ ðf1 ;f2 Þ ðHzÞ

30–35

40–45

50–55

60–65

1–5 5–10 10–15

-30 -30 -30

-30 -30 -30

-30 -30 -29

-30 -29 -27

Table 5 Extracted phase of wavelet for different windows using CDM and MCWC Window length (number of samples)

Extracted phase using CDM

Extracted phase using MCWC

1–1500 (total samples) 1–500 501–1000 1001–1500

60 60 70 75

62 58 50 54

Source Wavelet Phase Extraction

D. H. Naghadeh, C. K. Morley b

Figure 17 a Real seismic data set from New Zealand entitled ‘Bo_Parihaka, inline 81SY-11’; b obtained reflection coefficient from data set (a) using time variant sparse-spike deconvolution method; c extracted wavelets from different windows using CDM

wavelets with reverse polarities have equal maximum kurtosis value hence the correct wavelet polarity cannot be identified. Zero-phase wavelets with reverse polarities have minimum and maximum val-

Pure Appl. Geophys.

ues for combination of different-moments method. Those enable the technique to handle a finite data segment and it can also choose the correct wavelet polarity. Also, the existence of different moments can decrease sensitivity to outliers. When the effects of frequency bandwidth and analyzing window length are considered, the results of extracting phase information from data without and with noise demonstrate

Figure 17 continued

Source Wavelet Phase Extraction

Figure 18 a Phase corrected real seismic data set from New Zealand entitled ‘Bo_Parihaka, inline 81SY-11’; b phase corrected angles for different analyzing windows, each window includes 40 traces which include 500 samples

D. H. Naghadeh, C. K. Morley

that CDM is superior to the other methods discussed here.

Acknowledgments The authors thank GNS science (Institute of Geological and Nuclear Sciences of New Zealand, www.gns. cri.nz) for permission to use and show the real data examples. REFERENCES David, H.A. (1970), Order Statistics, John Wiley and Sons, New York. David, H.A. (1981), Order Statistics, 2nd edition, John Wiley and Sons, New York. DAUBECHIES, I. (1992), Ten Lectures on Wavelets, Philadelphia, PA, Society for Industrial and applied Mathematics. DOANE, D.P., SEWARD, L.E. (2011), Measuring skewness, A Forgotten Statistic? Journal of statistics education Volume 19. FOMEL, S., and VAN DER BAAN, M. (2014), Local skewness attribute as a seismic phase detector, Interpretation 2, no. 1, SA49–SA56. HOSKING, J. R. M. (1990), L-Moments: analysis and estimation of distributions using linear combinations of order statistics, J. Roy. Statist, Soc. Ser. B 52(2), 105–124. HERRERA, H.R., and VAN DER BAAN, M. (2012), Short-time wavelet estimation in the homomorphic domain, AAPG, GeoConvention, Calgary. AB. Canada.

Pure Appl. Geophys.

KIM, H.Y. (2013), Statistical notes for clinical researchers, assessing normal distribution (2) using skewness and kurtosis, Open lecture on statistics. LEVY, S., and Oldenburg, D.W. (1987), Automatic phase correction of common-midpoint stacked data, Geophysics. 52, No. 1, 51–59. ROBERTS, M.J., and RICCARDO, R.A. (1999), Student’s Guide to Analysis of Variance, London, Routledge. MALLAT, S. (1989), A theory of multi-resolution signal decomposition, the wavelet representation, IEEE Trans, Pattern Anal, Machine Intell, PAMI-11,674–693. MANSAR, S., and RODRIGUEZ, J.M. (1996), Analyzing and zero phasing seismic data using wavelet transform, Proc. EAEG. 96, Amsterdam, The Netherlands, 3–7 June. MUNOZ, R., RAPHAEL, E., and UNSER, M. (2002), Continuous wavelet transform with arbitrary scales and O(N) complexity, Signal Processing. 82, 749–757. ROBINSON, E.A., and TREITEL, S. (1980), Geophysical signal analysis, Prentice-Hall, Englewood Cliffs. N. J. RIOUL, O., and VETTERLI, M. (1991), Wavelets and signal processing, IEEE Signal processing, 8, no. 4, 14–38. SINHA, S., PARTHA, S.R., PHIL, A.D., and CASTAGNA, J.P. (2005), Spectral decomposition of seismic data with continuous-wavelet transform, Geophysics. 70, No 6. VAN DER BAAN, M. (2008), Time-varying wavelet estimation and deconvolution by kurtosis maximization, Geophysics. 73, NO. 2 MARCH–APRIL 2008; P. V11–V18, 5 FIGS. 10.1190/1. 2831936. VAN DER BAAN, M., and FOMEL, S. (2009), Non-stationary phase estimation using regularized local kurtosis maximization, Geophysics. 74, NO. 6, P. A75–A80, 5 FIGS. 10.1190/1.3213533. VELIS, R.D. (2008), Stochastic sparse-spike deconvolution, Geophysics. 73, NO. 1, P. R1–R9, 14 FIGS. 10.1190/1.2790584. WHITE, R.E. (1988), Maximum kurtosis phase correction, Geophysical Journal International. 95, 371–389.

(Received June 1, 2015, revised October 29, 2015, accepted January 8, 2016)