Enhancing time resolution by stabilized inverse filter and Q estimated on instantaneous spectra

RESUMEN

1 Notwithstanding, this model does not satisfy the minimum delay condition of dispersive medium, rigorously stated by the Kramers-Krönig relation.To correct phase distortions, and considering the narrow seismic broadband, equation 1 is rewritten as (Wang and Guo, 2004).(2) Here, ω h is the maximum frequency of the seismic broadband and ᵧ= Q r Equation 2 leads to Hence, wave equation ( ∂u (ω, r) ⁄ ∂r −ik(ω)u(ω, r) = 0 admits as a solution u(ω, r + ∆r) = u(ω, r)exp exp i (4) The first exponential compensates for the amplitude spectra, and the second one the k(ω) phase spectra, and downward after continuing from the surface t = 0 to t u(ω, t) = u(ω, 0)exp (5) In consequence where β (ω, t) = exp − Equation 5 can be solved like an inversion problem with a stabilization factor o introduced through u(ω, t) = i (7) The final solution is obtained by integrating equation 13 along the broadband For a horizontal-layer model, the interval quality factor of the nth-layer and the average quality factor of n layers have a cumulative relation given by: ( With the above equation, Q n is estimated where T n and T n-1 are the times where the average quality factors Q n and Q n-1 are measured.

Matching Pursuit decomposition
The MP algorithm (Mallat and Zhang, 1993) decomposes a signal into a linear superposition of wavelets selected from a redundant functions dictionary.The dictionary is built with previously calculated wavelets that are estimated and chosen because of their best fitting of the seismic trace.A time-frequency wavelet family is generated by scaling, translating, and modulating a continuous window function m(t) L² (R) such that ‖m‖ = 1,m(0) ≠ 0, ∫ m(t)dt ≠ 0 defined by With scale factor s, modulating frequency ƒ m and translation t d .Due to m(t) being even mᵧ is centered at the t d abscissa concentrating the energy in the vicinity of t d .By Fourier transforming equation 10 becomes

Introduction
A seismic wave that propagates in inelastic media is transformed attenuating its amplitude due to frequency-dependent absorption and delaying its phase due to different phase velocities.Many approaches have been applied to increase time resolution using spectral inversion (Castaño et al., 2011) or by spectral blueing with relative success (Kazemeini at al., 2010).On the other hand, the quality factor Q quantifies the attenuation as the percentage of energy lost in a cycle at a certain frequency (Kjartansson, 1979), and methods to apply inverse Q filtering have been developed and applied to seismic data (Margrave, 1998;Wang, 2002), and their performances compared (Montana and Margrave, 2005).The spectral ratio method (SRM) is the most common approach used to calculate the interval attenuation factor.Initially, it was widely used to estimate attenuation in acoustic VSP data, and recently in shear VSP data (Chuandong & Stewart, 2006).In addition, new methodologies have been implemented to calculate P-wave and S-wave attenuation (Qp and Qs) from standard well log data, relating these to the presence of oil or natural gas (Walls et al., 2006).It is known that attenuation in seismic data is related to gas-bearing reservoirs, which can cause significant changes in seismic amplitude with offset providing erroneous AVO results in deep targets.In this case, attenuation is a valuable alternative tool as a direct indicator of hydrocarbons (Kumar et al, 2003).Seismic attenuation has been used to distinguish changes in spectra in order to detect the presence of fluids and fractures in rocks (Ramirez-Cruz et al., 2005).
Therefore, a stable inverse Q filter might compensate attenuation in seismic data if a reliable pseudo Q log was available, otherwise trying to recover amplitude spectra and correct phase spectra by Q filtering could lead to possible artifacts on seismic data.The stable inverse Q filter must achieve the previous goals without boosting ambient noise and increasing S/N ratio.
The accurate determination of Q may provide a more accurate Amplitudeoffset analysis and reliable seismic inversion for reservoir characterization.Therefore, the accurate estimation of Q is a core challenge to overcome before applying inverse Q filtering, considering that Q does not depend linearly on frequency (Kolsky, 1956).The robustness of Q estimation through transforms that use fixed and variable time-frequency windows have been analyzed, pointing out that the precision and accuracy achieved depends on the transform used to supply the spectrum (Reine et al, 2009).The WFT cannot describe signal structures with varying size, because all wavelets have a constant scale that is proportional to the window size, whereas in MP, the signal structures are represented by wavelets that match their time frequency signature (Wang, 2007).In this case, the MPD uses a Morlet wavelet dictionary to decompose a seismic trace into a time-frequency spectrum with a higher resolution than those ones provided by GT and WFT.The matching algorithm is based in Morlet wavelets because they have shown high performance along with computational advantages (Liu & Castagna, 2005).
MP, windowed Fourier (WF) and Gabor transforms (GT) were coded in Matlab and tested on synthetic seismograms, which pointed out the best performance of MP to provide accurate Q logs.The t-ω spectrum supplied by MP is brought by a change of variable, in which Q is calculated using a linear regression.Finally, the estimated Q log is input to a Matlab program that implements the stable inverse Q filtering (Wang, 2006).Applied to a stacked offshore Colombian seismic line, the above two applications corrected the wavelet distortion, compensated the energy losses, and increased the S/N ratio.The trace-by-trace approach guarantees the reliability of the improvement in continuity of seismic events, being applicable to pre-stack data as well.
When ω is much greater than ω r (an arbitrarily low reference frequency), the following expression becomes true Where, m(ω) is the Fourier transform of m(t) with its energy concentrated around ƒ m.Henceforth a trace x(t) might be expressed by a linear combination of wavelets mᵧ which belong to the redundant mᵧ -family dictionary.The Morlet wavelet was selected due to its best fit of the seismic wavelets and may include quantitative attenuation and dispersion (Wang & Pann, 1996).In each step, the MP algorithm looks for the wavelet that provides most energy to the signal and removes this energy from the signal.The energy of the signal is iteratively reduced until reaching the threshold energy previously established the remainder energy is considered noise.

Methodology and results
A synthetic seismogram was created by superposing three traces.The first trace by convolving a 0°-20 Hz Morlet wavelet with 1 and -1 reflection coefficients at 0.3 and 0.9 S, the second by convolving a 45°-30 Hz Morlet with 1, 1 and 1 reflection coefficients at 0.45, 0.93 and 1.2 S, and the third by convolving a 90°-50 Hz Morlet with 1 and -1 reflection coefficients at 0.6 and 1.23 S. In figure 1A, three Morlet wavelets at 0.3, 0.45 and 0.6 S are clearly separated, and also two wavelets at approximately 0.9 and 1.2 S where it is hard to discern the two close events, at 0.9 and 0.93 S, and 1.2 and 1.23 S.
To compare the performances of Matching Pursuit (MP), Continuous wavelet (CW), Gabor (G) and windowed Fourier (WF) transforms, the synthetic seismogram was spectrally decomposed by the four methods.The supplied MP, CW, WF and G spectra are shown in figures 1B, 1C, 1D and 1E.As a result, the 20, 30 and 50 Hz bandwidths are best defined in the MP spectrum, besides a better time resolution of the reflectors at 0.3, 0.45, 0.6, 0.9, 0.93, 1.2 and 1.23 S. On the other hand, the WF spectrum strongly depicts Gibbs phenomena as a consequence of the windowed effects, even at 0.3, 0.45 and 0.6 S. χₐ is associated with the highest energy value in the new χ −spectrum, on which amplitudes decrease along the χ axis with a slope equal to Qֿ ֿ ֿ ¹ .A map of isovalue χ-curves over the t-ω spectrum is shown in Figure 2B.The amplitudes along each χ curve are summed and associated to A (χ) in the χ spectrum, as shown in Figure 2C.As a clarification, the abscissa of Figure 2B represents frequency (Hz) and not angular frequency (ω) and χ represents ft instead of ωt, so the χ =10 curve cuts the abscissa at 2.0 S and 5 Hz and achieves 100 Hz at 0.1 S.So henceforth, the χ abscissa, as in the case of Figure 2C, must be scaled by a factor of 2π.The maximum amplitude in the spectrum of Figure 2C is set to χ 0 = 4. Ahead of the highest peak of the χ spectrum, a very low amplitude zone around χ =10 is observed, which agrees with the low energy content in the vicinity of the χ =10 curve as is shown by the t−ω spectrum of Figure 2B.The χ spectrum is normalized in a new one, see Figure 3, where all values are negative and with the highest value around χ =5.The slope is steepest after χ b = 80, a value that defines the superior limit of the interval along which Q is estimated.Tiny values after χ b = 80, are rejected because they are associated to numerical noise observed in Figure 2C.A linear regression provided a value of 16.1 for the inverse of the slope that scaled by 2π gives a Q value of 101.16, which is close to 100.
The same seismogram was decomposed by WFT and GT supplying the respective χ spectra, observed in Figures 4 and 5. WFT and GT χ spectra depict a similar pattern that differs from the MP χ spectrum.In both new spectra, the low amplitude zone around χ =10 is not observed, and instead is full of spurious amplitudes.In addition, the absence of a linear decreasing trend in both spectra increases the uncertainty to estimate Q.In an attempt to estimate Q, some χ values were excluded to achieve a linear trend.In the WFT case, a Q value of 193.9 was obtained, and in the GT case a value of 134, both considered very inaccurate.The preceding procedure was applied to the migrated offshore section RC2103-129 of the Caribbean Colombia basin, formed by a siliciclastic Neogene sequence that overlays a carbonate and siliciclastic Paleogene sequence.A trace of this section was decomposed by MP, and its spectrum shown in Figure 6A.The 0-2S, 0-3.5 S and 0-5 S intervals were considered to estimate their associated average quality factors.Their respective χ -spectra are shown in figures 6B, 6C and 6D, on which linear regressions provided ‹Q› values of 39.9, 55.2 and 73.1, respectively.Using equation 9, interval Q values of 39.9, 113.5 and 296.3 related to the 0-2S, 2-3.5 and 3.5-5S intervals were obtained, with these values being in agreement with Q measurements reported for this lithology.
Next, Q values were estimated at 0.5, 1.5, 2.5, 3.0, 3.5, and 4.5 seconds in the migrated section, and in order to consider the Q lateral changes it was estimated in positions separated forty CDPs.The final Q model was achieved by smoothing laterally the obtained attenuation factors.An inverse Q-filter was developed and implemented in Matlab, which includes a stabilization factor σ compensating for the amplitudes according to equation 13.The filter was applied to the above section, and as a result, a higher resolution image was achieved, as observed in Figure 7B.The seismic events in the Q filtered image look sharper, with a better lateral continuity than those in the former section and are reliable due to the trace by trace approach of the applied inverse Q filter.An area, from 1000 to 1300 ms and between the CDPs 1870 and 1890, was selected to observe the Q-filtering performance.
Figures 8A and 8B show the zoomed areas before and after applying the inverse filter.The seismic events between 1050-1100 ms and 1150-1200 ms look clearly separated after Q inversion, whereas the weak event between 1100-1150 ms becomes stronger with an observable lateral continuity not visible before.After Q filtering, the image depicts higher amplitudes and better lateral continuity than before Q filtering.A spectral analysis was performed in the limited area.Figures 9A and 9B show the amplitude's spectra before and after applying the Q filter.After Q filtering, the spectrum increased in one order of magnitude, and the dominant amplitude moved from 10 Hz in the stacked section, to 28 Hz.Moreover, frequency content above 30 Hz was recovered increasing the frequency bandwidth.

Conclusions
Q filters, which are widely used to recover seismic frequency bandwidth, may create artifacts in seismic images due to phase distortions, so to increase S/N ratio, an inverse Q filter that compensates for the amplitude spectrum without introducing phase distortion was implemented.Therefore, the strict use of the Q filter implies an accurate estimation of the quality factor.To accomplish this, Q must be estimated on the instantaneous spectra provided by different methods.The attenuation factors are estimated in a normalized spectrum making the achieved value reliable.Spectral decompositions supplied by Matching Pursuit, windowed Fourier, and Gabor transforms were compared, and the analysis of the estimated attenuation factors noted the high performance achieved by the Matching Pursuit method.Applied to a stacked offshore seismic line off the Colombia basin, the inverse Q filtering provided an enhanced image, where the reflectors depict sharp reflectors with better continuity than the stacked image.A zoomed zone of the inverted image showed separated events, not clearly seen in the stacked section.An analysis of the spectra before and after inversion indicated a one-order increase in the amplitude scale, and the increase of the frequency bandwidth.The trace-by-trace approach makes the method suitable to pre-stack inversion.All the decomposition methods and Q filter algorithms were specially implemented for this research.Figure 9 Spectral analysis of the zoomed area A) before and B) after Q filtering where the spectrum is improved at least in one order of magnitude and the dominant frequency moves from 10 to 25 Hz, increasing the frequency bandwidth.

Figure 1 A
Figure 1 A) Synthetic trace and spectra obtained by B) Matching pursuit C) Continuous wavelet transform D) Gabor and E) Windowed Fourier transforms.Next, a blocky model with Q = 100 constant in depth was made, along which a plane wave traveling generates the attenuated synthetic seismogram of Figure 2A.This seismogram was MP decomposed, providing the t-ω spectrum shown in Figure 2B.Considering exclusively the attenuation phenomena, the wave amplitude decays as A (ω, t ) = A 0 exp [ − ωt / 2Q ]

Figure 2 A
Figure 2 A) A trace propagated along a simple Q=100 model, B) the χ curve map superimposing the t − f spectrum provided by MP and C) the χ spectrum obtained by replacing ft with χ and summing along constant χ curves in the t − f spectrum.

Figure 3
Figure 3 Normalized χ -spectrum with χ b limiting the zone where Q is estimated through a linear regression.

Figure 4 Figure 5
Figure 4The χ -spectrum supplied by WFT and B) zooming of the interval along which the linear regression provides Q estimation.

Figure 6
Figure 6 A) The MP spectrum of a stacked trace of the stacked section RC2103-129, B) the 0-2 S spectrum, C) the 0-3.5 S spectrum and D) the 0-5 S spectrum.

Figure 7 A
Figure 7 A) three χ -spectra used to calculate Q at 2, 3.5 and 5 seconds, B) Stacked section of the Colombia basin, C) Stacked section after Q estimation and Q filtering.The seismic events look sharper than before Q filtering.

Figure 8
Figure 8 A) The zoomed area before Q filtering, and B) the same area after applying Q filtering with lateral continuity increased and new discerned reflectors.