Spectral-ratio tomography for seismic attenuation estimation

When measuring surface seismic data, an accurate attenuation estimation method is necessary to compensate for the energy loss and phase distortion of seismic waves, and is also beneficial for further quantitative amplitude analyses and reservoir parameter predictions. For conventional Q-estimation methods (such as the log spectral-ratio (LSR) method and attenuated traveltime tomography), accuracy may be affected by the differences between the overburden ray paths of two selected reflections (we call it the overburden effect). In this study, we design a more accurate Q-tomography method to estimate Q-values (both in the overburden and target layer simultaneously) without overburden assumptions. We address the overburden effect by using an inversion method, which allows us to separate attenuation effects from the overburden through the traveltime differences in the tomography grid cells. We test the method on synthetic data and prove its feasibility and effectiveness by applying it to field data.


Introduction
It is necessary to estimate accurate Q-values to compensate for the phase distortion and the energy loss of seismic waves propagating through viscoelastic subsurface media. Ignoring compensation for attenuation may decrease the seismic data resolution and challenge quantitative amplitude studies. Many methods for Q-value estimation have been developed (Quan & Harris 1997;Brzostowski & McMechan 1992;Liao & McMechan 1997;Wang 2004Wang , 2014Rakshit et al. 2017;Jin & Sun 2017;Tary et al. 2017;Ilyina et al. 2018). Based on accurate Q-estimation results, Q-compensation migration (Mittet et al. 1995;Bano 1996;Wang & Guo 2004;Liu et al. 2017) and inverse Q-filtering (Wang 2002(Wang , 2006Wang et al. 2013;Zhang et al. 2015) can be applied. Q-estimation can also be beneficial for quantitative amplitude analysis. Q-tomography could be interpreted as a combination of tomography (Liao et al. 2017) and frequency-domain methods for interval Q-value estimation. Instead of estimating only one interval Q-value of a single layer, it can estimate multiple subsurface Qs through inversion (Li et al. 2016;Jin & Sun 2017). According to how the frequency-domain information is used, two types of existing Q-tomography method can be classified: indirect Q-tomography and direct Q-tomography. Quan & Harris (1997) and Liao & McMechan (1997) indirectly used the frequency-domain information (center frequency downshift) of seismic waves to reconstruct the attenuation distribution tomographically. In contrast, Brzostowski & McMechan (1992) introduced the amplitude-tomography technique for Q-value inversion. This method directly estimated attenuation based on the relation between the change in frequency-domain information (seismic amplitude spectral ratio) and the Q-value model. However, their algorithm does require absolute amplitude input, which is easily contaminated by geometric spreading loss, transmission or reflection loss, or other factors.
For the log spectral-ratio (LSR) method and attenuated traveltime tomography (t* tomography), the accuracy of the conventional methods may be affected by the overburden effect when there is a difference between the ray paths of the reference and target reflections (Reine et al. 2012). To address this problem, Reine et al. (2012) selected reflections that shared identical overburden ray paths. However, this kind of reflection is not easy to identify, especially for complex subsurface structures. Zhang & Ulrych (2002) introduced the peak frequency shift method to estimate Q-values in multiple layers using a layer-stripping approach. The Qvalue results of deeper layers are affected by the Q-value error estimated in the shallower layers. The robustness of the method is constrained by using only one trace for calculation.
We design a new attenuation tomography method for Qvalue estimation. We can invert for the Q-values both in the overburden and target layer simultaneously and separate the attenuation effects from the target and the overburden layers during the inversion. This method makes no assumption of the overburden structure and attenuation properties. This method can be regarded as a combination of LSR and tomography. We test the feasibility and effectiveness of the proposed method on synthetic and field data and find that it is more robust than conventional frequency-domain methods as it uses multiple reflections; it is also more reliable and accurate as it does not require assumptions about the overburden.

Spectral-ratio tomography (SRT) for Q-estimation
The amplitude of a seismic wave traveling through a viscoelastic medium will be attenuated as : where is angular frequency, A is the amplitude spectrum of seismic wave after traveling distance s, t is traveltime and S is the amplitude spectra of seismic wave at distance s 0 . Here we use G to stand for the frequency-independent attenuation instead of the PG value of Wang et al. (2015), which contains energy loss from reflection and transmission and also geometric spreading. For amplitude tomography, the G is not considered, thus the Q-inversion is contaminated by the non-intrinsic attenuation. For t* tomography based on the spectral ratio or centroid frequency shift, both the t* from real data and model are calculated and the initial Q model is updated by minimizing the misfit between the two t*. Both the methods neglect the overburden effects on the Q-inversion.
The proposed SRT method chooses two well-defined reflections at different traveltimes as the input to the spectral ratio from the same source. For a two-layer model as shown in figure 1, we assume the grid height is the same as the layer thickness and the number of grid cells is two (vertical dimension). Each layer has m grid cells (horizontal dimension), and the total number of grid cells is 2m. We pick reference reflections A 1 , A 3 , … A 2n-1 in the overburden, and target reflections A 2 , A 4 , … A 2n in the target layer along the horizontal direction.
Amplitude spectra of seismic waves reflected from points A 1 and A 2 in figure 1 can be expressed as where Q j is the Q-value in the jth grid cell, t j 1 is the traveltime in the jth grid related to reflection A 1 and t j 2 is the traveltime in the jth grid cell related to reflection A 2 . The logarithm ratio of these two amplitude spectra is: 393 where B = ln(G 2 /G 1 ) contains frequency-independent attenuation and Δt j 12 is the traveltime difference between the two reflections in the jth grid cell. Conventional t* tomography method neglects attenuation in the overburden (related to 1/Q 1 ,1/Q 2 , … , and 1/Q m in equation (4)) and invert for Q-values in the second layer, which leads to errors in the Q results. We will further discuss this later in section 4 (Overburden effects analysis).
Next, we select the other reflections from different offsets as shown in figure 1. The logarithm ratio of the amplitude spectra of these reflections can be written as where Δt j 2n,2n−1 (the superscript is the number of grid cell) is the traveltime difference between reflections A 2n and A2n−1 in the jth grid cell. When we select reflections, if the target layer reflections are not available, reflections from a deeper layer with rays traveling through the target layer can be used instead. We need the target layer information but are not restricted to target layer reflections.
The kernel of our method is that the differences in compensation factors t * 2n,2n−1 in Q-values for the overburden, which are neglected in conventional methods, are considered: If the ray paths of two reflections from a trace at fixed overburden location do not match, the overburden effect is considered in our method by adding this compensation term into the inversion. If there is no overburden traveltime difference, the compensation term becomes zero and the above equation simplifies to the conventional Q-tomography method.
We will use components within frequency band which have good signal-to-noise ratio and assume that Q-values are frequency independent in the seismic frequency range. Equations (5) and (6) can be represented in matrix where vector d contains spectra-ratio information (b 1 b 2 … b n ), and vectors v = (1 1 … 1) T, = (f 1 f 2 … f l ) T and z = (0 0 … 0) T are of length l (the number of frequencies). The vector m to be estimated is B and Q −1 . Thus, the geometrical spreading loss and transmission or reflection loss that contaminate the absolute amplitudes can also be derived. The method can be expanded for the N-layers model using the log ratio of the amplitude spectra of reflections in each layer. We minimize the least-squares objective function to solve equation (7) (Sacchi 1997).

Feasibility test and non-intrinsic attenuation immunity analysis
A four-layer model was used for the model test (with parameters in Table 1). The structure of Q-values and velocities can be found in figure 2a. A Ricker wavelet was chosen as the source wave, with 60-Hz peak frequency. Through ray tracing, travel paths and traveltimes in the grid cells were calculated. The shot gathers were simulated by convolution in the frequency domain using equation (1). Velocity dispersion was also considered using the frequency-independent Q-value model by Futterman (1962). The non-intrinsic attenuation was not considered here, and the frequencyindependent term B in equation (4) is set to 0. Figure  b and c show the synthetic shot gathers. The offset ranges were from 20 to 800 m with an interval of 20 m. We chose two shot gathers to perform the four unknown Q-values inversion. The two shot sources (with 78 traces that travel through the area) were located on the left and right side of the area. This geometry with opposite source location was retained for the following 2D model test to ensure that all the grid cells could be covered by the rays. The SRT method (actually all the existing tomography methods) needs the rays to travel through the grid cells as much as possible, which carry the information of the traveltime in each cell and the overall amplitude loss, for Q-inversion. The more the coverage of the grid cells by the rays, the more accurate the Q-inversion. Each model layer forms a grid cell and the total number of grid cells is four. The initial Q −1 values in the model are set to 0.01. The Q-values in the grid cells here were estimated simultaneously by the SRT method using four reflections in the shot gathers. The estimated Q −1 -results (figure 2d) are close to the real values, which shows the feasibility of our method. Next, the model was added with Gaussian random noise of signal-to-noise ratios (SNR) of 6, 4 and 2. The mean value of Q −1 was calculated by averaging the results of 200 independent runs. Figure 3 shows the mean values and standard deviations of the inverted Q −1 with different Gaussian random noise values.
Since the inversion results are dependent on the chosen frequency ranges, we invert for Q −1 within a range of bandwidths. The lower bound frequency is set from 1 to 40 Hz, and the upper bound frequency is set from 41 to 80 Hz. The inversion results are shown in figure 4 with mean values and standard deviations. Although the lower frequency bound 1 Hz gave the maximum errors of Q −1 , the mean Q −1 values match well with the true values.
By keeping other parameters unchanged, we allowed the Q-value along the offset in the four-layer model to change as shown in figure 5. Each layer has four grid cells (horizontal dimension), and the total number of grid cells is 16. To test whether non-intrinsic attenuation could introduce error to the proposed method, angle-dependent reflection and transmission coefficients were included. The frequencyindependent factor G n in equation (1) can be derived as where R, T and T' are the reflection coefficient, upgoing transmission coefficients and downgoing coefficients, respectively. The Zoeppritz equation (Aki & Richards 2002) was used to calculate the reflectivity and transmission coefficients. Figure 5 shows the Q −1 inversion results from amplitude tomography and the SRT method. There was no result of 395 amplitude tomography in the first layer as there was no source wave knowledge, while the SRT method could invert the Q −1 values in all layers simultaneously (figure 5a). Because of the disadvantage of amplitude tomography (requiring absolute amplitudes), changing the reflection and transmission coefficients caused obvious errors in the results. The estimated Q −1 of SRT tomography showed good agreement with the true values (figure 5b-d), which confirms the validity of our method. The effect of non-intrinsic losses (changes in reflection and transmission coefficients, figure 5e) was addressed by introducing vector B in equation (2) and using multifrequency components.

Overburden effects analysis
In this section, we quantitatively derive and analyze the overburden effects on the Q-tomography methods. The t* calculation in the tomography method may introduce errors by using reflections from the top and bottom of the target layer (figure 6). The t* tomography method assumes that there is no attenuation from the overburden on the reflections from the top, which is not reliable.
The t* of the two reflections in figure 4 can be expressed as where Q 1 and Q 2 represent the overburden and the target layer, respectively; and t AF , t FE , t AB , t DE , t BC and t CD are the traveltimes with different ray paths. Taking the logarithm ratio of the two amplitude spectra, we obtain: Conventional t* tomography neglects the second term related to 1/Q 1 in equation (11), which will introduce errors as there is always traveltime difference, t AF + t FE ≠t AB + t DE.  Therefore, we rewrite equation (11) The error factor u for t* is expressed as where t 2 = t BC + t CD , Δt 1 = t AB + t DE − t AF − t FE, and t * 2 = t 2 Q 2 . If Δt 1 = 0, the error factor is 1 and there will be no error in the t* calculation. Except in this case, there will be errors introduced by attenuation in the overburden. The error can be affected by the Q −1 values in both overburden and target layer, the distance of the selected reflections from the source and the overburden thickness. For small offsets and deep layers, the traveltime difference is small, while for larger offsets and shallower layers, there will be a large error in the Q −1 inversion. The general form of error factor in a multilayered model can be expressed as where Δt j is the traveltime difference in the jth layer, and t * n is the t* in the nth layer.
The t* error factors of different layers in the 1D model are shown in figure 7a. Without source wave knowledge, t* and the Q −1 value in the first layer cannot be calculated by t* tomography, while the SRT method can invert the Q −1 values in all the layers simultaneously. For all the t* error factors, with the increase in the distance between the receiver and source, the overburden traveltime differences increased, resulting in larger t* errors. The fourth layer was deepest, which gave the least overburden traveltime difference, thus the t* error here was the smallest. The error factor of the third layer was larger than that of the second because of the relative larger attenuation in the overburden and the correspondingly larger error factor of the third layer. As the offset increases, the error factor will decrease to zero (where the ray path of the two reflections are the same) and then become negative. The inverted Q −1 results of both methods are shown in figure 7b. The SRT method produced more accurate Q −1 397 Figure 7. (a) The t* error factor variation with distance for different layers and (b) comparison of Q −1 -inversion results using the SRT method and t* tomography.
inversion results and Q −1 values in the first layer, while t* tomography contained errors. The fourth layer generated the smallest error in accordance with the error factor. In the third layer, because some traces with large source to receiver distances generated negative error factors, the overburden effect balanced out from the positive error of small distance traces; thus, the third layer contained a relatively smaller error.
For the 2D model with Q −1 change in the horizontal direction, assuming there are m cells in each layer, so the t* error factor is where Δt j represents the traveltime difference in the jth cell. The general form of error factor in a multilayered model can be expressed as The Q −1 structure in the four-layer model with horizontal changes is shown in figure 8a. The changes in Q −1 in the horizontal direction modify the error factor curves (figure 8b). Compared to the t* error factor in figure 7a, the curve vibrates, and the positions of slope change are related to the change in Q −1 values. The t* results of the fourth layer of t* tomography contained the smallest error because of the smallest overburden traveltime difference.
Q −1 inversion results of the two methods for the four layers are shown in figure 9a-d. SRT method produced good results in the first layer, while t* tomography gave no result in this layer. The second and third layers showed large t* tomography errors due to large error factors. The decrease in the error of cells at a large distance was because of the source gather input from the right-side source location. The SRT method showed good Q −1 inversion results for all layers, and the overburden effects were well handled despite the horizontal variation in Q −1 values.
The model overburden was not just the first layer, but comprised all the layers that overlapped above the target layer. To prove this, we assumed the first layer in the fourlayer model had no attenuation (very close to zero) as shown in figure 10a. We expected no t* error factor (figure 10b) in the second layer. The third and fourth layers still had errors as there was attenuation in the overlying layers. Thus, in the Q −1 inversion results shown in figure 11a-c, the Q −1 results of t* tomography in the second layer are as accurate as in the SRT method. There are still errors for t* tomography in the third and fourth layers, while the SRT method produces good results.
Then, we analyzed how the overburden Q-values affected the error factor. Here, we assumed that the attenuation was high in the second layer as shown in figure 12a, which may relate to a gas or oil reservoir. Due to the relative low attenuation in the overburden, the t* error in the second layer becomes small. This is similar to the no attenuation case in figure 12b, but the other two layers have more errors. Thus, the error in the Q −1 inversion results, as shown in figure 13a-c, is relatively small in the second layer, but large in the third and fourth layers. There was no overburden effect in the SRT method even with an embedded high-attenuation layer.
Next, we set the third layer as the high-attenuation layer as shown in figure 14a. The Q −1 values were the same as those in the second layer in figure 10a. The t* error factor in the third layer was small due to its high attenuation (figure 14b). The error factor of the second layer was not affected by this high attenuation, and was the same as that of the original model in figure 8b. Although the traveltime difference was small in the fourth layer, the error increased due to the high overburden attenuation. Figure 15 parts a-c show the Q −1 inversion results of the two methods. The t* tomography gave small errors in the third layer, but large errors in the other two layers, which corresponded to the t* error factors. The SRT method gave good results despite the high-attenuation layer set in a different location. The complex Q model (figure 16a) was used to test the ability of the proposed method to distinguish the embedded high-attenuation thin layer and four high-attenuation heterogeneous bodies below. The velocity model is shown in figure 16b. A Ricker wave with a 60-Hz dominate frequency is used for modeling. The model consisted of 40 × 40 grid   The SRT method and t* tomography were both applied for attenuation inversion. The inversion results are shown in figure 17. For t* tomography, the traces that offset are larger than 280 m, so they are not used for inversion to avoid unacceptable errors in the results. Thus, the resolution is reduced due to less trace input. Due to the overburden effect, the interface of the high-attenuation thin layer was hard to  distinguish using t* tomography. The four high-attenuation heterogeneous bodies could be distinguished but contained approximately 15-35% errors, especially for the deeper right two. An interesting abnormality of the t* tomography was that there were four symmetric errors right below those with small or even negative Q-values. Some linear errors were due to the geometry. The SRT method, which accounted for the overburden effect, gave high-quality results with a clear thin-layer boundary and satisfactorily located the four highattenuation heterogeneous bodies.

Field data test
The proposed method was then applied to field data acquired in northern China. Before Q-estimation, some essential denoising processes (ground roll suppression, random noise suppression) and geometrical spreading compensation were applied.
The inversion used 25 shots with 1990 traces. As shown in the stack gather (figure 18), the locations at 500, 950, 1150 and 1850 ms are used as the main interfaces to calculate the subsurface Q −1 values. The Q −1 results calculated by the t* tomography and SRT methods are shown in figure 19 parts a and b. For t* tomography, Q-values above 750 m are not calculated and set to be zero. The resulting Q-structures from the t* tomography contains more abnormal negative values due to the overburden effect.
Next, the Q −1 results (at location shown in figure 17) were calculated using logging data from the Koesoemadinata-McMechan rock physics model (Koesoemadinata and McMechan 2001). This model is used in the studied area as an empirical Q prediction model for Q −1 values in the frequency range from 1 to 100 Hz. Figure 20a shows the Q −1 results by this model at depths between 1500 and 2500 m. In addition, the Q −1 results by the t* tomography and the SRT method are also shown for comparison. The SRT method match better with the empirical model, especially   for depths between 1800 and 2300 m. The t* tomography contains more errors due to the overburden effect. Figure 20b shows the well-to-seismic calibration before and after using Q −1 results from the SRT method as guidance for stretching the events. The two events indicated by the red arrows at depths between 2200 and 2300 m are stretched 6 and 8 ms, respectively. The velocity dispersion calculated is roughly 1.2%. By setting the frequency of sonic logging frequency at roughly 12 500 Hz and the dominated frequency of seismic data at 40 Hz, the Q −1 value was calculated based 402 Figure 18. The stack gather shows the interfaces used for Q −1 inversion. on the velocity dispersion (based on equation 7 by Li et al. 2015). The calculated Q −1 result was roughly 0.0065, which was consistent with that of the empirical model (averaged 0.0066) and the SRT method (averaged 0.0059), while the t* tomography gave average Q −1 value of 0.0045, which is less than the SRT method. More accurate Q-estimation results are useful for well-to-seismic calibration. In turn, the effectiveness and robustness of the SRT method were proved.

Discussion
The Q-value resolution (the grid size) is determined by the number of effective traces used for inversion. The proposed method uses all the ray paths from all the shots simultane-ously and defines a smaller grid size for higher resolution. Meanwhile, we use multi-frequency components of the spectral ratio, which improves the inversion robustness. We do not need to choose an optimistic frequency band for the t* regression as in t* tomography. The accuracy of the results is also dependent on seismic wave raytracing accuracy.
This method addresses the overburden effects without requiring overburden assumptions. Although not mentioned, the conventional LSR is also affected by the overburden effects.
Equation (4) can be simplified as where error factor u is expressed as where t 2 = t BC + t CD, Δt 12 = t AB + t DE + t BC + t CD − t AF − t FE . When there is no difference in the ray paths in the overburden Δt 1 = 0 and the term Δt 12 = t 2 , the error factor is 1 and there will be no error in the Q calculation. The error of LSR can be affected by the Q-value ratios, the position of the selected reflections, and the overburden thickness. The general form of the error factor in a multi-layer model can be expressed as where Δt j is the traveltime difference in the jth layer and t n is the traveltime in the nth layer. The Q error factors are shown in figure 21.

Conclusions
A more accurate Q-estimation method, called SRT, is developed in this study, which is less affected by the overburden effect. The effect is addressed through inversion by separating the attenuation effects of the target layer from the overburden based on traveltime differences in the tomography grid cells. The proposed method can handle the overburden effect despite the vertical or horizontal change in Q −1 values, even if there are high-attenuation layers and heterogeneous bodies. Subsurface Q-values can be inverted simultaneously including those in the overburden. The proposed method can avoid the effect of non-intrinsic losses by introducing a frequency-independent term and using multi-frequency components. The components within frequency band that have good signal-to-noise ratios can be used as input. For thin interbed formation with seismic response of composite 404 waveforms, time-frequency analysis methods are suggested when windowing the reflections. The model and real data test prove that the proposed method is more robust compared to t* tomography by using multiple reflections. In addition, traces with large offset can be used for inversion, which improves the resolution of Q −1 results.
The resolution of SRT method is controlled by the size of the grid cell, which is dependent on whether there are enough traces that travel through the cells can be used for inversion. Further work will be carried out on the subject of frequencydependent Q-inversion by the SRT method. A Q-value can describe the difference of velocity between log and seismic data, so the knowledge about how the Q behaves between the high and low frequency band is essential. Since the SRT method considers multiple frequency components in the inversion process, it has advantages in handling Q variation versus frequency.