1 INTRODUCTION

Seismic data is nonstationary series whose amplitude and phase both change with travel time. The time-varying character is significant to predict the subsurface geological features. Moreover, as the target depth and complexity increases, the difficulty of reservoir prediction increases also in turn asks for higher resolution and inversion accuracy. The accuracy of seismic wavelet estimation affects greatly both the seismic section resolution and inversion result. However, the most of existing seismic wavelet estimation focus on the stationary wavelet. We studied the time-varying seismic wavelet estimation method.

The effect of inaccurate wavelet estimation on seismic interpretation has been studied. For instance, the minor error of wavelet phase estimation can cause severe inversion errors (Yuan et al., 2011; Luo et al., 2014). White (2003) also pointed out that a small wavelet time-lapse causes a large error in amplitude attribute slice, which brings difficulties in interpretation. An accurate wavelet estimation result can not only improve the accuracy of well-tie and formation calibration but more importantly, can provide the base of deconvolution and inversion. Gao (2009) introduced the mathematic expression of the varying wavelet model of the seismogram and improved the seismic resolution by adaptive time-frequency decomposition. Han et al. (2015) showed that the wavelet phase changes with the existence of fluid because of attenuation by examining the relation between wavelet phase and fluid prediction with both simulation and areal data example. Moreover, studying the time-varying wavelet can make a contribution to the study of attenuation and compensation of seismic data to improve the resolution of the seismic profile. The lateral continuity of seismic profile is influenced by wavelet's side-lobe (Yuan et al., 2013). Deconvolution can eliminate the wavelet effects, which in turn corrects the phase distortion and increases the seismic profile's bandwidth (Wang et al., 2015). Other than deconvolution, the inverse Q-filtering can also increase the resolution by compensating the attenuation. Van der Baan (2012) compared the results of deconvolution and inverse Q-filtering, found that phase only inverse Q-filtering is good for dispersion correction while time-varying deconvolution is better for compensation of attenuation. He suggested that both methods should be used successively in real data processing.

The existing seismic wavelet estimation methods can be classified into two categories (Tygel et al., 2000): (1) Deterministic methods estimate wavelet by matching seismic data to the well-log data. The resulting wavelet can be used to correct seismic profile formation calibration. However, in most real synthetic profile calibration processing, the so-called "stretching" or "squeezing" methods are not sufficiently scientific and lack quality control (Li et al., 2008). Edgar et al. (2011) introduced the following steps. Firstly, estimate the seismic wavelet amplitude spectrum and then add zero-phase to it. Secondly make synthetic seismic data based on the reflectivity from well-log data and estimate zero-phase wavelet. Finally match the synthetic data to the real one and adjust both the amplitude and phase till the best matching result. Based on the above steps, Wang et al. (2015) introduced the ant colony optimization algorithm to determine wavelet phase effectively.

Walden and White (1998) proposed a coherency matching technique, which used auto-and cross-power spectra of well-log and seismic data to estimate wavelet without any hypothesis on the seismic wavelet. However, this method estimates the stationary wavelet. White et al. (2014) also used this method on broadband seismic data and studied the low-frequency spectral decay of the seismic well-tie wavelet especially. Based on the coherency matching technique, Ma et al. (2010) put forward the wavelet estimation by well-log, VSP and seismic data.

(2) Statistical methods analyze the reflectivity's statistic characters to estimate wavelet with a certain hypothesis, for example assuming the white reflectivity and minimum-or zero-phase. Wiggins (1978) introduced a classical statistic method, the kurtosis maximization. Kurtosis is used to evaluate the series non-Gaussianity in statistics. Other than kurtosis criterion, there are lots of criterions for deconvolution result evaluation (Yu et al., 2012), for instance, the parsimony criterion (Claerbout, 1977) and exponential transform criterion (Ooe et al., 1979). Wang et al. (2014) pointed out that the kurtosis maximization makes inaccurate phase estimation of complex seismic data such as low dominant frequency data. They then constructed a new Rxenyi divergence-based criterion which showed better phase detection results than kurtosis for low dominant frequency data. Liu et al. (2014) combined the kurtosis maximization and local similarity maximization with the envelope to estimate the time-varying phase of well-log and seismic data. Ma et al. (2015) compared the results of wavelet phase estimation between kurtosis maximization and skewness maximization and found that the latter has higher dynamic range under the same phase variation and is more robust for noisy data even both methods are influenced by the distribution character of reflectivity.

Based on the studies by Levy et al. (1987), White (1988) and Longbottom et al. (1988), Van der Baan(2008, 2009, 2010) improved Wiggins' method and proposed the non-minimum phase estimation method. Hypothesizing white and non-Gaussian reflectivity, rotating the phase and the wavelet phase can be determined when the seismic data reaches the maximum non-Gaussianity. The wavelet estimation can be completed by combining the estimated phase with the amplitude spectrum estimated by seismic data autocorrelation. Fomel and Van der Baan (2014) proposed a more stable time-varying locally observed phase estimation method called the local skewness maximization. This method uses local skewness, a smooth nonstationary measure (Fomel, 2007a, 2007b), to test Gaussianity. Local skewness is defined by using regularized least-squares inversion. In this article, we tried to introduce skewness maximization into the time-varying propagating wavelet phase estimation. The locally observed wavelet is the response of the local geological features, while the propagating wavelet is the comprehensive response of the subsurface media, geometric spreading, attenuation and scattering.

We proposed a time-varying wavelet hybrid estimation method by combining the advantages of the co-herency matching technique and skewness maximization, using coherency matching technique to estimate the time-varying wavelet's amplitude from well-log and seismic data, while using skewness maximization to de-termine the time-varying phase from seismic data directly. Comparing to the conventional Wiener-Levinson method (Robinson et al., 2000) and constant phase rotation method (Levy et al., 1987; Longbottom et al., 1988; White, 1988; Gao et al., 2010; Wang et al., 2014), we don't have any hypothesis on wavelet. The synthetic examples showed the correct estimation result of both wavelet's amplitude and phase. Then we applied this method to the real seismic data and used quality control, including the qualitative analysis by well-tie results and quantitative analysis by PEP and NMSE calculation.

2 METHOD 2.1 Time-Varying Wavelet Amplitude Estimation

From the well-log data and seismic data, we have the following relations.

(1)

(2)

(3)

The true seismic reflection signal s_t is the convolution of true reflectivity r_t and seismic wavelet w_t. The synthetic reflectivity from well-log data x_t contains certain random noise u_t, and so is the real seismic signal y_t. We assumed that the random noise u_t and v_t are independent and smooth. We used the spectral theory of stationary and stochastic processes to solve this noisy system. Based on the Eqs.(1)–(3), we can get the auto-and cross-power spectra of x_t and y_t.

(4)

(5)

(6)

where Φ_xx(f), Φ_yy (f), Φ_uu(f), Φ_vv (f), Φ_rr (f) and Φ_ss(f) are the auto-power spectra of x_t, y_t, u_t, v_t, r_t and s_t, respectively, while Φ_xy (f) is the cross-power spectra of x_t and y_t.

From Eqs.(4)–(6), three variables are known from well-log and seismic data, while four are still unknown which are Φ_uu(f), Φ_vv (f), Φ_rr (f) and W (f). Therefore, we introduced one more variable v(f) to solve the wavelet amplitude W (f).

(7)

v(f) is defined as the noise power ratio between the well-log and seismic data and can be calculated by cross-correlation between multi-traces of seismic data (Walden and White, 1998).

During the W (f) estimation, the smoothing window is needed to minimize the error caused by spectra estimation of finite time series. However, the smoothing window can also cause some distortion. If the smoothing window is not long enough, it results in severe distortion and deviation, otherwise, it results in not enough error minimizing. Hence, tests on the balance point are necessary.

2.2 Time-Varying Wavelet Phase Estimation

The conventional phase correction method is kurtosis maximization method, which was first applied to blind deconvolution by Wiggins (1978). Then Levy et al. (1987) and White et al. (1988) made the constant-phase assumption that reduced the number of free parameters to only one; hence ensure the stableness of the inversion. The statistical rationale behind the Wiggins algorithm is that convolution of any filter with a time series that is white with respect to all statistical orders makes the outcome more Gaussian (Donoho, 1981). Based on this theory, wavelet phase can be determined by phase rotation while evaluating the Gaussianity. The seismic data rotation in time-domain can be expressed as the Eq.(8).

(8)

where ϕ is the rotation angle and H[·] is the Hilbert transform.

The third-order statistics, skewness, evaluates the non-Gaussianity by evaluating asymmetry. Comparing with the kurtosis, skewness has great advantages on phase detect application, including much higher dynamic range and fewer sample points (Fomel et al., 2014).

(9)

where skew(y) is the skewness of seismic data y(t) and N is the sample point. The bigger the skewness absolute value is, the higher the non-Gaussianity of series is. When applying a series of constant phase rotation from –180° to 180° to the seismic data, the angle that produces the most skewed distribution is the wavelet phase.

Based on the correlation coefficient of series a_n and b_n,

(10)

We can get skew²(y),

(11)

Here, we followed the local attribute estimation method (Fomel, 2007a), which uses regularized least-squares inversion and has better robustness than conventional sliding window method. The long or short radius of the regularization smoothing determines whether the estimated wavelet phase is global or local, respectively. Long regularization smoothing averages the local geological characters and results in the propagating wavelet phase.

2.3 Time-Varying Wavelet Estimation

Finally, combine the estimated amplitude and phase to get the estimated time-varying wavelet.

(12)

where |W_i(f)| is the final estimated wavelet in window i. |W_i(f)| and ϕ_{skew, i} are the estimated wavelet amplitude and average phase in the window i. sgn(f) is sign function.

2.4 Quality Control

We evaluated the estimated result from both quality and quantity aspect. On one hand, convolving the estimated wavelet with the synthetic reflectivity from well-log data makes the synthetic seismic data ŝ(t). We can then evaluate the match between the synthetic ŝ(t) and real seismic trace y(t) to evaluate the wavelet estimation result.

On the other hand, there are two factors, the PEP and NMSE, that can quantitatively evaluate the wavelet estimation result. The proportion of energy in the surface seismic data segment y(t) predicted by the well-log synthetic ŝ(t) (PEP) is a measure of goodness-of-fit (Ma et al., 2010). Higher PEP means better match between ŝ(t) and y(t), hence better wavelet estimation result.

(13)

Normalized mean square error (NMSE) means the error between synthetic ŝ(t) and real reflection signal s(t). Lower NMSE means more accurate wavelet estimation.

(14)

However, the real reflection signal s(t) is unable to know. Based on the prediction theory, NMSE has an approximate expression (Walden et al., 1984).

(15)

where T is the length of matched seismic data and b is the analysis bandwidth (White, 1984). From Eq.(15), multiple tests of smoothing window width can be made to get the most appropriate value based on the higher PEP and lower NMSE.

2.5 Time-Varying Deconvolution

The deconvolution is straightforward after wavelet estimation. Based on the least-square inverse filtering, filter coefficients F in frequency domain is

(16)

where c² is a stable factor that means the noise level of seismic data and superscript * means complex conjugate.

Time-varying deconvolution contains filter coefficients that change with time. According to Eqs.(12) and (16), time-varying filter coefficients can be achieved. Then applying these coefficients to the entire seismic trace and put all these deconvolved results together using the following function:

(17)

d_i(t) is the deconvolved result by estimated wavelet i and d(t) is the final time-varying deconvolution result (van der Baan, 2008).

3 MODEL TEST

We tested the time-varying characters of wavelet's amplitude and phase separately with different models.

3.1 Time-Varying Dominant Frequency and Fixed Phase Wavelet Model

The first synthetic seismic trace is made by convolution of reflectivity and time-varying dominant frequency and fixed phase wavelet. The reflectivity is a random time series. The wavelets have a dominant frequency of 60 Hz initially and transform linearly to 15 Hz at the end of the travel time. The 0° and 90° phase wavelets are considered separately. The total travel time is 2.1 s with a sample rate of 2 ms.

Figure 1a is the random time series, which convolves with the 0° and 90° phase wavelets and produces the synthetic seismic trace of Fig. 1b and 1c, respectively. Both traces show the reducing resolution caused by reducing dominant frequency wavelets, while 0° phase trace shows better symmetric character than 90° phase trace, for instance at 0.4 s. After time-varying wavelet phase estimation, both estimated wavelet phases match well with the true phases with less than 20° error (Figs. 1d and 1e). In order to evaluate the accuracy of wavelet estimation, we convolved the estimated wavelet with the true reflectivity and matched it with real synthetic seismic trace. Moreover, we calculated the PEP and NMSE. Fig. 2 is the result of 90° wavelet estimation and the match evaluation. Over the interval of 0.33~0.45 s, 0.93~1.15 s and 1.65~1.77 s, the three estimated wavelets clearly show the reduction of dominant frequency. The PEP for three match evaluations is 86.6%, 99% and 92%, respectively. Similarly, the result of 0° wavelet estimation and the evaluation also comes well.

3.2 Time-Varying Dominant Frequency and Phase Wavelet Model

Next, the synthetic seismic trace is made by convolution of reflectivity and time-varying dominant fre-quency and phase wavelet. The wavelet has dominant frequency and phase of (60 Hz, 0°) and transforms linearly to (15 Hz, 90°) at the end of the travel time. This trace simulates the real seismic data suffered from all the attenuation factors underground. The total time is 2.1 s with sample rate of 2 ms. Fig. 3 shows this synthetic trace result, which also shows the reducing resolution caused by the wavelet's reducing dominant frequency.

Figure 4 shows the time-varying wavelet estimation results at 9 instant times from 0.3~2 s. For the time-varying amplitude estimation, the smoothing window is Papoulis window of 120 ms. The result wavelets have 60 ms length. For the time-varying phase estimation, the radius of the regularization smoothing was 200 samples or 0.4 s. The estimated phase matches well with the true phase. The estimated wavelets' waveforms accord well with the true waveforms. In order to compare with the existing wavelet estimation method, we used coherency matching technique on the total trace of 2 s and got the result of stationary wavelet in Fig. 4d. Because the true wavelet varies with time, the stationary wavelet estimation result has a certain deviation, especially for the phase.

Figure 5 shows the time-varying deconvolution result. After deconvolution, the trace improves the reso-lution clearly, which helps thin layer identification, for example at 0.3 s and 0.6 s. Moreover, the time-varying deconvolution improves seismic trace's symmetry by rotating to zero phase, for instance at 1.1 s.

4 REAL DATA

The field data is from well TZ72 at Tazhong low uplift area. Seismic data is a time migrated post-stack section with sample rate of 4 ms. For the well-log data, certain pre-processing is necessary. For instance, CheckShot from VSP can correct acoustic log data which greatly improves the wavelet estimation result. Fur-thermore, the density log data is missing at shallow layers. The missing part can be made by acoustic log data. Then we can use the corrected acoustic and density log data to construct the synthetic reflection series. Possibly resampling is necessary to match the sample rate of real seismic data. Fig. 6 is the result of resampled reflectivity from TZ72.

Using well-log data and seismic data to estimate the time-varying amplitude, while estimating the time-varying phase directly from seismic data. Fig. 7 shows the result of time-varying phase searching. The radius of the regularization smoothing was 0.5 s. The result shows the phase of seismic data varies greatly in this area. Therefore, the time-varying deconvolution is necessary in order to meet the demand of fine reservoir interpretation. Fig. 8 shows the result of time-varying wavelet estimation over the interval of 0.94~1.44 s, 1.83 ~2.33 s and 2.56~3.06 s. The evaluation of estimated wavelet over 2.56~3.06 s shows the satisfactory result, with PEP of 72.9% and NMSE of 0.11, considering the real seismic data having low SNR and coming from complex geological features.

5 DISCUSSION

The method proposed in this article has two hypotheses: (1) reflectivity is white and (2) reflectivity is non-Gaussian. Based on the analysis of well-log data, reflectivity underground can basically meet the non-Gaussianity demand (Walden et al., 1986). Skewness attribute needs asymmetric reflectivity. The real reflec-tivity underground generally increases with depth; hence the possibility of positive numbers is much higher, which in turn results in the asymmetric reflectivity distribution and then the effectiveness of skewness maxi-mum method in practice. However, some part of reflectivity series may appear to be symmetric or Gaussian in practice, which can be verified by well-log data.

Basically, the non-Gaussianity character of reflectivity can be met in practice. However, the reflectivity is not white but blue instead, which means the power spectrum increases with frequency in limited bandwidth (Walden et al., 1985). This issue can also be solved by well-log data (Saggaf et al., 2000). Van der Baan (2008) made more detailed analysis on this issue.

Fomel and Van der Baan (2014) introduced the local skewness maximization method to detect the phase variation of the locally observed wavelet. They pointed out the difference between locally observed wavelet and propagating wavelet. The propagating wavelet is the physical wavelet that propagates through the subsurface media and reflects comprehensively the geology, geometric spreading, intrinsic and apparent attenuation and scattering. The locally observed wavelet is the instant wavelet as observed at a certain point in space and time. The locally observed wavelet is subject to the current shape of propagating wavelet and current reflectivity. The propagating wavelet is the input of deconvolution and follow-up inversion, while the locally observed wavelet is the supplementary of interpretation for its reflection of local geological structure and layer. In this article, we tried to determine the propagating wavelets' phase.

This time-varying wavelet estimation method needs pre-processing of well-log and seismic data. The real well-log data includes noise and errors which influence the wavelet estimation result (White, 1997; Yang and Yin, 2008). Therefore, it is necessary to make corrections on well-log data and de-noise and migration on seismic data.

6 CONCLUSION

We have presented a hybrid method to estimate time-varying seismic wavelet, without any hypothesis on wavelet, using coherency matching technique to estimate the time-varying amplitude and skewness maximization to estimate time-varying phase. In order to evaluate the estimation accuracy, we then make synthetic matching from convolution of the estimated wavelet with reflectivity made by well-log. The two quality control factors, PEP and NMSE, can also be calculated to quantify the evaluation. Furthermore, the estimated wavelet can be considered as the input of the time-varying deconvolution to improve the seismic section's resolution. In synthetic and field data examples, our approach exhibits effective time-varying wavelets estimation on the nonstationary seismic data.

ACKNOWLEDGMENTS

This work was supported by the National Basic Research Program of China (2013CB228602) and National Science and Technology Major Project (2016ZX05004003-002).