The traditional convolution model of a seismic trace is often stated as (Robinson, 1967) y s t a t t = w t * r t + e t = ∫ − ∞ ∞ w t − τ r τ d τ + e t , (1) where t and τ are the time indexes, y s t a t t is the stationary seismic trace, w t is the seismic wavelet, r t is the seismic reflectivity, and e t is the random noise term. We can write a matrix equivalent expression for Eq. 1, as y = W r + e , (2) where W is a Toeplitz matrix formed from seismic wavelet; y , r and e are the vector forms of the seismic trace, the seismic reflectivity series and the random noise term respectively.

We first convolute the reflectivity r t and delta function δ t to derive the mathematical formula of the non-stationary convolution model (Gholami, 2015) x t = ∫ − ∞ ∞ δ t − τ r τ d τ . (3)

Using x t to replace r t in Eq. 1, we have y s t a t t = w t * x t + e t . (4)

We can write a matrix equivalent expression for Eq. 4, as y = W x + e . (5)

Considering the attenuation effect of seismic wave, the attenuation coefficient (Aghamiry and Gholami, 2018) is introduced into the delta function to obtain δ ∧ f , τ = α f , τ δ ∧ f , τ − 1 = ∏ k = 1 τ α f , k δ ∧ f , 0 , (6) where f is the frequency index, δ ∧ f , τ is the frequency spectrum of the delta function, δ ∧ f , τ = F t → f δ t , τ , and F t → f denotes the Fourier transform with respect to t. δ ∧ f , 0 is the frequency spectrum of the delta function at the initial time. Using the constant Q theory (Aki and Richards, 1980), the attenuation coefficient can be defined as follows (Aki and Richards, 1980; Gholami, 2015): α f , τ = exp − π β f Q τ − i 2 π f N , (7) where α f , τ is the attenuation coefficient on the propagation time τ , β f is a function of frequency and its role is to impose some physical constraints, N is the number of data samples, and Q τ is the seismic quality factor which is a function changing with the propagation time. The average quality factor can be defined as follows (Margrave et al., 2011): t Q a v e t = ∑ k = 1 K Δ t k Q k , (8) where Q a v e t is the average of the seismic quality factor as a function of the propagation time, k = 1, … , K, K is the total number of assumed Q layers, Δ t k and Q k are the traveltime and Q value through the kth layer, respectively.

In Eq. 7, β f = f − i H f , where H is the Hilbert transform operator. We can rewrite Eq. 7 into the following matrix form (Aghamiry and Gholami, 2018) Δ ^ = α f , τ 0 α f , τ 1 . . . α f , τ N − 1 , (9) where Δ is the attenuation operator in the time domain and Δ ^ is the frequency domain representation of it. Figure 1A shows the attenuation matrix Δ for a constant Q = 30, and Figure 1B shows some columns of this matrix. This figure suggests that the shape of the propagated delta functions is severely affected by the attenuation effect. According to Eqs 3, 9, we have x = F − 1 Δ ^ r = Δ r , (10) where F − 1 is the inverse Fourier transform matrix. Substitute Eq. 10 into Eq. 5 and obtain the matrix form of the non-stationary convolution model y = W Δ r + e . (11)

We construct a dictionary matrix corresponding to the Q model. Let G = W Δ . (12)

Substitute Eq. 12 into Eq. 11 and obtain y = G r + e . (13)

The dictionary matrix G is related to Q value. The sparsity of the seismic reflectivity series can be used to determine the corresponding dictionary matrix G in the optimal case (Aghamiry and Gholami, 2018), so as to obtain the quality factor. We first assume that G is known, then the solution to the seismic reflectivity series r can be converted to an optimization problem, which is expressed in the following form (Gholami, 2015; Aghamiry and Gholami, 2018) a r g m i n r 1   s u b j e c t   t o y − G r 2 2 ≤ ξ , (14) or a r g m i n y − G r 2 2 + λ r 1 , (15) where ξ is an error bound of the data, λ is the regularization parameter and p denotes the p-norm. Due to the sparseness of the reflectivity series, it is necessary to use the sparse norm as much as possible to constrain it. The value of the p-norm in [0, 1] can ensure its sparsity. The value of the reflectivity series always lies between −1 and 1, so we use p=1/2 to measure the sparsity of the reflectivity in this paper.

In the inverse Q filtering method, the first step is usually to estimate the Q value, and then use the estimated Q value to perform inverse Q filtering, so the inverse Q filtering result depends on the estimated Q value. In this paper, we adopt a similar matching pursuit algorithm in dictionary learning to solve Eq. 14. The traditional matching pursuit algorithm needs to construct the dictionary matrices in Hilbert space, and each column of the dictionary matrices is dealed with normalization. In this paper, however, we construct a complete dictionary based on the theory of the non-stationary convolution model, and compute the reflectivity series under different dictionary matrices with the corresponding referencing Q values. Finally, we select the optimal dictionary matrix by comprehensively analyzing the residual and reflectivity sparsity, so as to obtain the seismic reflectivity and Q concurrently. The complete dictionary can describe changes in seismic amplitude over time, which is more in line with the laws of physics.

In the specific implementation process, the first step is to construct a three-dimensional dictionary matrix G d i c as shown in Figure 2. The first dimension of this matrix is related to the value of the quality factor. Slicing along the first dimension is a two-dimensional N×N matrix G , with the N th column representing the seismic wave at the propagation time N ms under the quality factor corresponding to the slice. According to the theory of the average Q value, we construct the dictionary matrix G by extracting the corresponding column from the matrix G d i c for the multi-layered signal (multi-segment Q value), which avoids the reduplicate calculation for constructing the matrices in every loop and significantly improves the computational efficiency. The second step is to initialize the reference Q values. The exact approach is to select a certain number of Q values in a reasonable range in logarithm scale as the reference Q values. The third step is to perform the subsection process for the seismic signals. The seismic signal is divided to several segments according to the rough seismic horizon to estimate the reflectivity and Q. For each seismic signal segment, the similar matching pursuit algorithm is used to compute the reflectivity series under the different referencing Q values, so as to obtain the corresponding residual ε = y − G r 2 2 and the sparsity of the seismic reflectivity series σ = r 1 / 2 . ε * and σ * are defined to represent the normalization of ε and σ respectively, and we can construct the following formula to evaluate the estimation results ξ = 1 − ϕ ε * + ϕ σ * , (16) where ξ is a synthetic criterion, and ϕ is an adjustable constant, measuring the weight of sparsity versus residual when evaluates the estimation results. A pseudo-code for this procedure is presented in Algorithm 1.

Algorithm 1

A pseudo-code for seismic reflectivity and Q estimation.

1 Input: y , the number of signal segment n s , the number of the reference Q n q

We first test the performance of the new method for seismic reflectivity and Q estimation of the constant-Q model. Some synthetic traces have been generated by convolving the reflectivity with a Ricker wavelet of dominant frequency 35 Hz and then contaminated by the random noises. The recording time is 1 s with the sample interval at 1 ms. These wavelet and recording parameters are used for all the synthetic examples in this paper.

Figures 3A–C shows the simple synthetic traces attenuated by using the constant-Q model with different Q values (Q = 30, Q = 60 and Q = 110) and contaminated by the 0.1% random noises. Here, the reflectivity series of 0.5, −0.6 and 0.7 are set at 250, 500, and 750 ms respectively, and no interference occurred between the reflections. The new method is applied to these simple traces for the reflectivity and constant Q estimation. Thirty logarithmically spaced Q values in the range [10, 120] are selected to run Algorithm 1. The estimated reflectivity series are shown in Figures 3D–F. As seen from these figures, the estimated reflectivity series are very close to the true values. The corresponding estimations of the Q values are showed in Figures 3G–I, showing a relatively good match with the true Q values, as shown in these figures by the red lines for comparison.

Figures 4A–C shows the complex synthetic traces attenuated by using the constant-Q model (Q = 60) and contaminated by random noises of different levels (0.1%, 5% and 10%). Here, the high reflections are set at 333, 500, and 666 ms, and the corresponding reflection coefficients are 0.75, −0.75 and 0.75, respectively. At other times, the reflection coefficients are set randomly. There are some reflection interferences from the adjacent reflections in the complex synthetic traces. These complex traces are used to test the new method for the reflectivity and constant-Q estimation. Thirty logarithmically spaced Q values in the range [10, 120] are also selected to run Algorithm 1. The estimated reflectivity series are shown in Figures 4D–F. As seen from these results, the new method improves the resolution and gives satisfactory results in the case of complex waves. The corresponding estimations of the Q values are showed in Figures 4G–I, which show good consistency with the true Q values.

We next test the performance of the new method for seismic reflectivity and Q estimation of the interval-Q model. Similar to the previous examples, some synthetic traces have been generated by convolving the reflectivity with a Ricker wavelet of dominant frequency 35 Hz, and then different random noises have been added to the resulting traces. The new method is applied to these synthetic traces for the reflectivity and interval-Q estimation. Algorithm 1 is performed for thirty logarithmically spaced Q values in the range [10, 120].

Figures 5A–C shows the simple synthetic traces attenuated by using the different interval-Q models and contaminated by the 0.1% random noises. Here, the reflectivity series of 0.3, 0.4, −0.6, 0.7, −0.6, 0.4 and 0.3 are set at 111, 222, 333, 500, 666, 777, and 888 ms respectively, and no interference occurred between the reflections. The estimated reflectivity series are shown in Figures 5D–F. The corresponding estimations of the Q values are showed in Figures 5G–I.

Figures 6A–C shows the complex synthetic traces attenuated by using the interval-Q model (Q = 60, 30, 110) and each contaminated by random noises of different levels (0.1%, 5% and 10%). Here, the high reflections are set at 333, 500, and 666 ms, and the corresponding reflection coefficients are 0.7, −0.7 and 0.7, respectively. At other times, the reflection coefficients are set randomly. The complex synthetic traces contain some reflection interferences from the adjacent reflections. The final obtained reflectivity series are shown in Figures 6D–F. The corresponding estimations of the Q values are showed in Figures 6G–I.

As seen from these results, the estimations of the reflectivity series are very accurate and the estimated Q values are relatively close to the true values in all the cases. Furthermore, the accuracy of the estimates is related to the noise level, which decreases as the noise level increases under normal circumstances. The new method improves the resolution by removing the wavelet and attenuation effects even in the case of high-level noises, especially at places with great attenuation effects.

For comparison, we also perform the Gabor deconvolution (Margrave et al., 2011) for the reflectivity estimation corresponding to the synthetic traces shown in Figures 6A–C. The Gabor deconvolution estimates and corrects for the effects of source wavelet and anelastic attenuation (Margrave et al., 2011). The corresponding estimated reflectivity series are shown in Figure 7. It is clearly seen from Figures 6, 7 that the accuracy of the reflectivity estimation by the Gabor deconvolution is not better than that by our proposed method.

In addition, we compare the proposed method with the conventional spectral ratio method (Hauge, 1981) for the Q estimation. The results of the Q estimation are reported in Table 1. It can be seen from Table 1 that the proposed method provides more accurate Q values than the spectral ratio method in all cases. The new method uses both the amplitude and phase information of the attenuation operator, and regularizes the estimated Q by the sparsity information of the seismic reflectivity series. However, the conventional spectral ratio method only uses the amplitude information to determine the Q value, and the spectral interference from the adjacent reflections limits it (Xue et al., 2020). Hence, the Q values estimated by the new method are more accurate than the results provided by the conventional spectral ratio method.