Introduction

Because of the anelasticity and heterogeneity of earth medium, seismic waves experience absorption when traveling through the earth (Wang et al., 2018b; Kneib and Shapiro, 1995; Zhang and Ulrych, 2002; Wang, 2002; Ricker, 1953; Futterman, 1962; Ruan and Zhou, 2012; Wright and Hoy, 1981; Lawrence et al., 2006; Ruan and Zhou, 2010; Gan et al., 2016). The attenuation of seismic waves is a property of the earth, and is quantified by Q, the quality factor. Seismic attenuation affects both amplitude and waveforms of the recorded seismic data. The accuracy and stability when estimating the quality factor have effects on processing quality, seismic imaging, amplitudes interpretation and reservoir description (Wang et al., 2018a; Guo et al., 2014). Q can be used to enhance seismic resolution by inverse Q filtering or Q-compensated migration (Delle Piane et al., 2014; Dutta and Schuster, 2014; Kjartansson, 1979; Zhang and Ulrych, 2007; Zhu, 2014; Li et al., 2015; Wang, 2011; Carcione, 2007; Matsushima et al., 2015; Zhu and Sun, 2017; Toverud and Ursin, 2005; Wang et al., 2015; Prudencio et al., 2013; Carcione et al., 2002; Zhu, 2015; Li et al., 2018), and can also be used to aid AVO analysis because the attenuation effects are superimposed on AVO signature and can be corrected if accurate Q values are available for compensation. Besides, Q can be potentially used as a hydrocarbon indicator (Hedlin et al., 2001; Wang, 2014,2004). In addition, Q can be used to provide information about lithology, porosity, and fluid or gas saturation of subsurface to facilitate seismic interpretation (Cheng and Margrave, 2012).

The algorithm for estimating Q from recorded seismic data has been investigated exclusively in the literature. Q can be estimated in time domain, frequency domain and time-frequency domain, each domain has its own advantage and limitation. In time domain, the pulse amplitude decay (PAD) (Brzostowski and McMechan, 1992), pulse rising time (PRT) (Kjartansson, 1979), and pulse broadening (PB) (Wright and Hoy, 1981) are the most common approaches. In the frequency domain, those approaches include spectral ratio method (SRM) (Hauge, 1981; Raike and White, 1984), centroid frequency shifting (CFS) (Quan and Harris, 1997) and peak frequency shifting (PFS) (Zhang and Ulrych, 2002; Gao and Yang, 2007). All these frequency domain approaches require Fourier transforms to calculate the frequency spectra of seismic records sampled within a time window (Yang et al., 2009). The time-frequency domain approach includes peak scale variations (PSV) (Li et al., 2006) in the wavelet domain to estimate Q by assuming an idealized pulse as the seismic source wavelet. However, the difference between the real source wavelet and an idealized pulse may be substantial, which is the limitation of the PSV method. There also exists a type of Q estimation approaches based on the instantaneous frequency (IF) of seismic record. In this type of method, Q can be estimated by inversion from the relationship between the measured instantaneous spectra and the seismic attenuation (Tonn, 1991; Engelhard, 1996; Barnes, 1991). Matheney and Nowack (1995) proposed a IF matching (IFM) method, which can iteratively modifies the causal attenuation operator (Aki and Richards, 2002) so that the IF at the envelope peak of the reference pulse is the closest to that of the target pulse and thus Q can be estimated. Although IFM requires a time-window and an iterative procedure, it does not have the requirement of selecting variable frequency bands as required by SRM.

Several researchers compared different Q estimation methods in the literature. Jannsen et al. (1985) compared PRT, SRM, spectrum modeling method and wavelet modeling method, and concluded that the wavelet modeling method is relatively better. Tonn (1991) compared the Q estimation approaches based on zero-offset VSP data, and concluded that no single method is generally superior. The performance of different Q estimation approaches depends on the specific situation where they are used. For example, analytical signal method is better for true amplitude recording, and SRM is the best for noise free data. Yang et al. (2009) compared four Q estimation methods, including SRM, CFS, PFS and wavelet envelope peak instantaneous frequency (WEPIF) methods according to their performance regarding wavelet independence, noise resistance and resolution of thin beds. Their comparison results showed that SRM is the best in both wavelet independence and noise resistance. WEPIF is the best in obtaining a highest resolution for a wedge model. Reine et al. (2009) compared four time-frequency transforms used for SRM. Their results show that time-frequency transforms with a systematically varying time window, such as the S transform and continuous wavelet transform (CWT), can obtain more robust estimation of Q (Du et al., 2010).

By reformulating the SRM in the time domain, Cheng and Margrave (2012) proposed a matching filtering method (MFM). Instead of calculating a spectral ratio, the Q estimation turns to finding a convolution operator in the time domain. The MFM is theoretically same to the SRM because the inverse Fourier transform of a spectral ratio is a matching filter. However, the time domain MFM is more robust in the presence of noise than direct spectral division in the frequency domain. Based on the concept of frequency down-shift, Raji and Rietbrock (2013) proposed a new mathematical model to achieve a reliable Q estimation, which shows less sensitivity to noise interference and change of frequency bandwidth. Their inversion process mainly consist of two stages: (1) computation of the centroid frequency for the individual signal using variable window length and fast Fourier transform; (2) estimation of the difference in the centroid frequency and traveltime for a paired incident and transmitted signals.

From the literature review, we can conclude that the noise level existing in the seismic records is the main effect on the Q estimation performance. In this paper, we propose a multi-channel Q estimation approach using a regularized spectral ratio method, which can utilize the spatial coherence of useful reflection signals to alleviate the influence of strong random noise. The proposed method can help obtain an accurate Q estimation result even from the raw seismic data without the need of denoising. We formulate the spectral division as a regularized inversion problem and use shaping regularization with a local smoothing constraint to solve the optimization problem. The smoothness constraint is applied both along the frequency and space directions. According to the conclusion from Reine et al. (2009) and Du et al. (2010), we propose to use the S transform (Stockwell et al., 1996) as the optimal time-frequency transform for the Q estimation. We organize the paper as follows: we first introduce the state-of-the-art single-channel Q estimation method based on the spectral ratio method. Then, we introduce the multi-channel Q estimation method based on the regularized spectral division framework, and we also introduce the shaping regularization method to solve the inversion-based regularized spectral division problem. Next, we test the proposed multi-channel Q estimation method on several synthetic and real data examples in the section of EXAMPLES. We also discuss the issues regarding the influence of noise and denoising, and computational complexity in the section of INTRODUCTION. Finally, we draw some key conclusions at the end of the paper.