Matching pursuit based on the rotated time-frequency atomic dictionary and its applications to seismic non-stationary data

The time-frequency spectrum of an attenuated seismic wave rotates due to the influence of intrinsic quality factor. Traditional time-frequency analysis methods are either limited by time-frequency resolution or incapable of rotated time-frequency representation. To solve it, a matching pursuit (MP) method based on the rotated time-frequency atomic dictionary (MPFR) was proposed. By introducing a frequency rotation factor to the four-parameter dictionary in the conventional MP method, a five-parameter dictionary with frequency rotation characteristics was constructed, and the calculation process of each parameter was given. We tested the method on the synthetic traces and field data sets. The results showed that the proposed method can effectively realize rotated time-frequency characterization of attenuated seismic waves. Compared with the conventional atoms, the rotated time-frequency atoms are more compliant to the local features of non-stationary seismic data. Moreover, the effective description of rotated time-frequency characteristics can be used as an auxiliary technique to predict the reservoirs related to the attenuation. A care has to be taken when thin layers are encountered, since wavelet interference may cause rotated time-frequency characteristics independent of dispersion and attenuation in MPFR method.


Introduction
Numerous petrophysical evidence has shown that, due to the friction between rock particles or fracture surfaces and relative movement of fluid, such as squirt flow, Biot fluid flow and fluid diffusion in porous media ( Johnston et al. 1979;Müller et al. 2010), attenuation occurs during the propagation of seismic waves. The attenuation can be partitioned to a phase and amplitude component, which leads to the variation of amplitude, frequency and phase characteristics, results in the non-stationarity of seismic waves and introduces challenges to the processing and interpretation of seismic data (Yuan et al. 2016).
Non-stationary signal analysis has been widely used to accurately describe the non-stationary features of seismic data and achieve the precise characterization of small-scale objects (Li et al. 2016b;Yuan et al. 2019a). By transforming the one-dimensional signal in the time domain into the two-dimensional time-frequency domain, local anomalies related to the target of interest were obtained (Oliveira et al. 2010). The representative methods include short time Fourier transform (STFT) (Zhong & Huang 2010), continuous wavelet transform (CWT) (Kazemeini et al. 2009;Chakraborty & Okaya 1995), S transform (ST) (Stockwell et al. 1996;Li et al. 2016a), Wigner-Ville distribution (WVD) (Claasen & Mecklenbräuker 1980) and matching pursuit (MP) (Mallat & Zhang 1993;Wang et al. 2016).
Assuming that the non-stationary signal is stationary within short time intervals, STFT describes the local features of the signal by introducing a sliding analysis window. However, due to the fixed window function, the resolution of STFT is invariant, which means, once the window function is determined, the time-frequency resolution is the same irrespective of frequency (Kwok & Jones 2000). In real applications, it is necessary to have a good transcendental understanding of the signal, so as to select appropriate time-window parameters and obtain relatively satisfied time-frequency resolution.
CWT decomposes the signal into the sum of a series of wavelet functions. By implementing the scale factor, the multi-resolution analysis is realized. However, previous studies (Puryear et al. 2012) have shown that CWT has highfrequency resolution and low temporal resolution at low frequencies, while it has low-frequency resolution and high temporal resolution at high frequencies. The mutual restriction between time-domain and frequency-domain resolution conforms to the uncertainty principle (Yang et al. 2014).
ST is an alternative multi-resolution analysis method, which is developed on the basis of CWT. Because of the timeshift and scale-change factors in the basic wavelet of CWT, its direct calculation is a time-scale result. To get the timefrequency result conveniently, ST defines the basic wavelet as the combination of Gaussian functions and complex exponential harmonics, which can establish direct connection with the Fourier spectrum (Stockwell et al. 1996). However, the time-frequency resolution is also limited by the uncertainty principle.
WVD is a quadratic transformation, defined as the Fourier transform of the instantaneous correlation function. By using the conjugation of the signal as the analysis window, WVD can effectively eliminate the influence of external time-window and avoid the mutual restrictions between the time and frequency resolution in the uncertainty principle. Thus WVD can achieve the optimal performance of timefrequency energy concentration for the single-component signal. However, because WVD makes use of the autocorrelation of the signal itself, cross-term energy interference was introduced in analysis of multi-component signals (Li & Zheng 2007;Yuan et al. 2019a). To suppress the cross terms, a time-window function, such as pseudo-WVD (Lerga & Sucic 2009;Wu & Liu 2009), is usually adopted. The cross-term interference can be solved to some extent at the sacrifice of resolution.
MP belongs to the adaptive method with parameterized time-frequency analysis, which decomposes the signal into the linear combinations of time-frequency atoms associated with the local structure of the signal (Mallat & Zhang 1993). Since the atoms are independent on each other and the WVD of a single atom does not arouse a cross-term, the timefrequency distribution of the signal can be restored by superimposing the WVD of each atom. As a result, it greatly avoids the influence of the cross terms and solves the timefrequency resolution problem of multi-component signals (Li & Zheng 2007).
Based on the time-frequency spectrum obtained by MP, Castagna et al. (2003) and Wang (2007), respectively, showed the low-frequency shadows associated with sandstone and carbonate gas-bearing reservoirs and applied it to gas-bearing reservoir detection. Based on the different variation features of spectrum at different time-space locations, Wang (2012) exploited the spectrum variation pattern for the first time to predict coalbed methane reservoirs. Chen & Song (2008) extracted the information about instantaneous frequency and bandwidth to predict the submarine reservoirs. Zhang et al. (2012) identified the full pattern of the deltaic sandstone based on instantaneous spectral features. Liu & Marfurt (2007) extracted the peak frequency and amplitude, and combined them with the coherence volume, which not only characterized the weak lateral variation in thickness of channel, but also predicted the width successfully.
Based on the decomposition results of seismic data through MP, Wang (2010) extracted the strong reflection waveforms of coal seam and stripped them from the seismic profile, which effectively enhanced the weak reflection information related to the reservoirs beneath the coal seam. Cao et al. (2017) identified micro-seismic signals during fracturing reservoirs and improved signal-to-noise ratio (SNR) by means of local correlation spectral constraints, which reduced the probability of misjudgment in micro-seismic events. Deng et al. (2017) improved the accuracy of fault recognition by separating the continuous reflection from the discontinuous residuals in seismic section. Zhang et al. (2017) realized attenuation compensation and dispersion correction of the prestack gathers with the amplitude variation of matching wavelets. Zhang et al. (2013) solved the dynamic stretching problem in the NMO (Normal Moveout) correction of prestack gathers by replacing the original oneby-one sampling point correction with one-by-one matching wavelet correction. Özbek et al. (2009) characterized seismic data by defining the sinusoidal basis functions related to wave number and realized the regular processing in the irregular space sampling. Yuan et al. (2017) found that wavelet interference gives rise to the equivalent Q-filtering effect and proposed to build the wavelet library with Q-dependent timevariant wavelets. Then they adopted sparse Bayesian theory to sequentially add or delete time-variant wavelets and applied to high-quality reflectivity imaging successfully.
The validity of the MP depends on the similarity between the time-frequency atoms and the local features of the analyzed signals. If the similarity is high enough, the sparse 823 Figure 1. Time-frequency relationship of the attenuated full-band signal with the frequency-independent Q and the frequency-dependent Q by amplitude and isolines. (a) Amplitude representation with the frequency-independent Q = 100. (b) Contour representation with the frequency-independent Q = 100. (c) The relationship of the Q in terms of frequency in the modified Kolsky model when Q r = 100. (d) Comparison of attenuation amplitude contours when the Q is frequency-independent (blue) and frequency-dependent (red), respectively. decomposition can be achieved by using a small number of time-frequency atoms. Otherwise, a large number of atoms are usually needed to reconstruct the signal. Moreover, the decomposed atoms cannot characterize the local characteristics of the signal effectively. Consequently, the definition and selection of the time-frequency atomic dictionary are of importance. The conventional time-frequency atoms usually consist of four parameters: scale factor, center time, center frequency and phase factor. For the seismic data with relatively gentle dispersion and attenuation, the four-parameter atoms can usually meet the requirements of seismic trace decomposition. However, in the presence of strong attenuation, such as gas-bearing reservoirs, the time-frequency variation becomes more intense. Thus, four-parameter atoms may fail to effectively describe the time-frequency characteristics in seismic data. To solve this issue, a five-parameter time-frequency atomic dictionary is constructed in this paper by introducing the frequency rotation factor. Meanwhile, MP based on the rotated time-frequency atomic dictionary is proposed to realize the sparse decomposition and timefrequency characterization of non-stationary seismic signals. Wang (2002Wang ( , 2006 gave the expression of seismic wave propagation and attenuation compensation in the frequency domain,

The time-frequency relation of the attenuated seismic waves
where U(t, ) denotes the seismic wave field at time t and radial frequency , Δt denotes the travel time interval between t and (t + Δt), h denotes the highest frequency in seismic frequency band, Q r denotes the intrinsic quality factor of the strata related to reference frequency and i is the imaginary unit. The first and second exponential terms in equation (1 exponential terms characterize the non-stationarity of seismic wave propagation. According to equation (1), every frequency component of seismic wave will be affected by amplitude attenuation and phase change, particularly for the high-frequency component. Assuming the Q is frequency-independent, figure 1a and b show the time-frequency relationship of the full-band signals at varying times in the case of the Q as 100. For the same travel time, the attenuation of the high-frequency component is stronger than that of the low-frequency component. For the same frequency component, the attenuation goes up with the travel time increases. Combining these two factors, the time-frequency relationship presents the contour map in figure 1b. On the whole, the time-frequency attenuation relationship is curved. But within a certain local timefrequency window, the relationship can be approximately linear and the rotation angle versus the coordinate axis varies with the window location: (i) Relative to time axis, the rotation angle in the high-frequency zone is larger than that in the low-frequency zone for a fixed time window; (ii) Relative to frequency axis, the rotation angle of time-frequency in the high time-window area is larger than that in low time-window area for a fixed frequency window. Overall, the longer the travel time, the more obvious the frequency rotation feature. The higher the frequency, the more obvious the time rotation feature.
For the frequency-dependent Q, several models have been proposed, such as the Kolsky model (Kolsky 1953), the Strick-Azimi model (Strick 1967;Azimi et al. 1968), Azimi's second and third models (Azimi et al. 1968) and Müller's power-law model (Müller 1983). Wang & Guo (2004) analyzed these different models and proposed the modified Kolsky model, which can more accurately describe the velocity dispersion in seismic frequency band. The expression is Corresponding to the frequency-independent Q = 100, the Q r in equation (2) is set to 100 and the relationship of the Q in terms of frequency is shown in figure 1c. As the frequency goes down, the Q goes up. As the frequency goes up, the Q gradually goes down with a more gentle rate of decrease. In the frequency range of 0-500 Hz, the Q changes from 102 to 100. Figure 1d shows the comparison of attenuation amplitude contours between the frequency-independent Q and the frequency-dependent Q. The two are very close. The reason is that although the Q varies with frequency, the magnitude of change is very small. Compared with the frequency-independent Q, the largest difference in the frequency-dependent Q occurs at low frequencies with a difference of only 2, which is much less than the Q values, and the difference decreases with the increasing frequency ( figure 1c). With such a small change of the Q value, the variation of attenuation and dispersion is also very small. So is the disparity in the difference of amplitude attenuation isolines between the frequency-dependent Q and the frequency-independent Q. Therefore, for the frequencydependent Q, the time-frequency relationship can also be approximately linear in a small time-frequency window within a seismic frequency band.
To further analyze the time-frequency relationship under the influence of attenuation, two kinds of wavelets commonly used in seismic exploration, Morlet and Ricker (Liu & Marfurt 2007;Wang 2007), are taken as examples to simulate the attenuation of seismic wave propagation when the Q is ∞ (infinity), 200, 100 and 50. Taking into account the advantages of WVD for single-component signal analysis, figure 2 parts b and d, respectively, calculate and superimpose WVD of the independent waveforms in corresponding traces and time positions in figure 2 parts a and c, forming the timefrequency distributions of seismic waves under the influence of different Q values. Meanwhile, to ensure the effective display of all waveforms and spectra within the same amplitude range, only the time range of 0.6-2 s is shown. In addition, the amplitudes of spectra are treated using decibels (Huang et al. 2018), which is expressed as 20 × log 10 (Amp), where the Amp represents the amplitude.
When the Q = ∞, the amplitude and phase terms related to the Q in equation (1) are both 1; that is, all frequency components of seismic waves are not affected by attenuation, so the waveform will not change during propagation. As the Q value decreases, the influence of attenuation gradually appears. From the perspective of waveform, (i) for the same Q value, the amplitude of seismic wave decreases as the travel time increases and the waveform gradually varies from symmetry to asymmetry; (ii) for the same travel time, the smaller the Q value, the stronger the attenuation of the waves. Corresponding to the waveform, the time-frequency spectrum is also affected. When the Q = ∞, the spectrum remains unchanged. As the Q value decreases, the attenuation becomes much clearer: (i) for the same Q value, the spectral amplitude decreases with the increasing propagation time and the spectral morphology rotates gradually; (ii) for the same travel time, the smaller the Q value is, the more serious the attenuation and the greater the rotation will be. In summary, the longer the propagation time and the smaller the Q value, the stronger the amplitude attenuation and waveform variation of the seismic wave, and the greater the amplitude attenuation and rotation of the time-frequency spectrum (figure 2).

Construction of the rotated time-frequency atomic dictionary
To characterize non-stationary signals, Mallat & Zhang (1993) proposed the concept of time-frequency atoms, which are expressed as where u is the time delay, is the spread in the time axis, is the phase shift, w is the time window and g (t) denotes time-frequency atom, which is characterized by four parameters, = (u, , , ). By changing the parameters, a series of time-frequency atoms can be generated and then combined to form the time-frequency atomic dictionary. Compared with the basis functions of STFT, CWT and ST, the parameters in equation (3) are independent from each other and the way to generate atoms is relatively flexible. Therefore, parameters can be designed and optimized according to the features of signals in real applications. So, a time-frequency atomic dictionary more suitable for the signals can be established to realize the best description of the local structure.
Due to the advantages of the Gaussian function in timefrequency resolution, Gaussian window is usually used to replace the window function in equation (3), forming the expression of Gabor atoms (Wang 2007) (4) To better characterize seismic signals, Morlet et al. (1982) proposed a form of wavelet similar to equation (4), forming the Morlet wavelet, whose expression is A series of Morlet wavelets can be generated by frequency variation, time translation, scale variation and phase rotation from equation (5). Thus, the Morlet time-frequency atomic dictionary can be formed to represent seismic data (Wang, 2007).
For the Ricker wavelet, its basic expression is (Liu & Marfurt 2007) Since the Ricker wavelet belongs to the zero-phase wavelet, non-zero phase is usually introduced by means of Hilbert transform ( to enrich the wavelet series, where H[⋅] denotes the Hilbert transform. Similarly, the Ricker wavelet dictionary can be generated by changing the parameters in equation (7). Unlike the Morlet wavelet, the Ricker wavelet is not constructed by the combination of a window function and a complex exponential function, but is equivalent to the second derivative of a Gaussian function (Wang 2015). Therefore, Ricker wavelet dictionary has no scale factor. Instead, it adjusts waveform expansion in time direction by the frequency, .
In the MP method, the time-frequency atomic dictionary generated by a Morlet or Ricker wavelet is usually used to realize the basic characterization of seismic data. spectrum of Morlet and Ricker wavelets after phase rotation, respectively. One can note that different phases correspond to different wavelet shapes, but their time-frequency spectra are the same, as expected. In addition, the phase change does not cause the rotation of time-frequency spectrum. Because of the symmetry of the Gaussian function, other parameters such as frequency change, time translation and scale change do not contribute to the time-frequency rotation either. Therefore, it is difficult to represent the time-frequency spectrum of attenuated signals in figure 2 by the atoms using conventional wavelet dictionary. A new wavelet dictionary needs to be constructed to realize the effective characterization of the rotated time-frequency spectrum.
On the basis of the linear frequency-modulation wavelet transform, Yang et al. (2014) proposed the concept of frequency rotation (FR) operator, which is expressed as where c denotes the FR factor and Φ R (t) denotes the FR function. In fact, by calculating the first derivative of its phase, the instantaneous frequency IF(t) of Φ R (t) can be obtained as follows Equation (9) indicates that the instantaneous frequency of Φ R (t) is a linear function of time and c represents the time-changing rate of frequency. When c = 0, Φ R (t) becomes a constant and the instantaneous frequency does not change with time. When c ≠ 0, the instantaneous frequency changes linearly with time and its graph will have an angle with the time axis, which corresponds to the FR. Figure 4 shows the real part of Φ R (t) and its time-frequency spectrum when c > 0. The oscillation of the waveform increases with time, and the spectrum rotates relative to the axis. Because c > 0, both the rotated angle of the time-frequency spectrum in equation (8) and the frequency axis are negative (Yang et al. 2014). When c < 0, both the rotated angle and the frequency axis will be positive.
To establish an atomic dictionary with rotated timefrequency characteristics, equation (8)  conventional wavelet function, so that equation (3) becomes (10) Accordingly, equations (5) and (7) become and respectively, where = (u, , , c, ) represents the basic parameters of the Morlet wavelet dictionary. For the Ricker wavelet dictionary, the parameter is neglected.
Therefore, time-frequency atoms with different characteristics can be generated by changing the parameters. Figure 5a, from the left to the right, shows the time-domain waveforms and time-frequency spectra of the initial Morlet wavelet after scale change, time shift, frequency change, phase rotation and FR. Analogously, figure 5b demonstrates the waveforms and spectra of the Ricker wavelet with different parameters. Note that the waveform and time-frequency spectrum are changed in shape and position due to the variation of conventional parameters, but they do not induce time-FR. After involving FR, the rotated time-frequency characteristics are shown. Based on this, the rotated time-frequency atomic dictionary is generated and the effective representation of attenuated seismic data is realized.

Matching pursuit based on the rotated time-frequency atomic dictionary (MPFR)
Once the atomic dictionary is determined, seismic data decomposition can be accomplished by using MP. The where s(t) denotes the seismic trace, g n (t) denotes the wavelet generated by equations (11) or (12), n = (u n , n , n , c n , n ), a n denotes the amplitude and R (N) s(t) denotes the residual trace after the s(t) is reconstructed by N wavelets. To obtain the optimal decomposition, the wavelet g n (t) should represent the trace optimally and the residual energy reach minimum in each iteration. Therefore, the wavelet parameters and amplitude are determined by the following equations (Wang 2007): where P denotes the wavelet parameters of the atomic dictionary and <,> denotes the inner product. Equation (13) suggests that MP belongs to an iterative algorithm. The main challenge is the accurate determination of wavelet parameters. For the rotated time-frequency atoms, there are more parameters (a n , u n , n , n , c n , n ), which need to be calculated step by step in each iteration: (1) Calculate the wavelet parameters (u n , n , n ). By using Hilbert transform to calculate the complex seismic trace, then the time u t , instantaneous frequency t and instantaneous phase t at the maximal envelope are extracted as the time delay u n , frequency n and the phase factor n , respectively. For a Ricker wavelet, special frequency processing (Liu & Marfurt 2007) is required, which is n = ( √ ∕2) t . (2) Calculate the wavelet parameters (c n , n ). Set the range of FR factor c s and scale factor s , respectively. The optimal parameters are obtained by scanning. The basic principle is (c n , n ) = arg max For the Ricker wavelet, the scaling factor n in the above equation can be ignored. As for the range of s , its minimum value should be slightly less than 1 and maximum value should be slightly greater than 1 according to Wang (2007). About the range of c s , the minimum value can be set to zero. In this situation, the time-frequency atom becomes the atom of conventional MP method, which does not have the rotated time-frequency characteristic. For the maximum value of c s , it can be set according to the specific research problems. However, the larger the value, the greater the numbers of calculations.
(3) Wavelet parameters optimization. Since steps (1) and (2) are carried out independently, the calculated wavelet parameters may not reach the local optimum. So, parameter optimization is needed to make the results satisfy equation (14). Generally, the results of steps (1) and (2) are taken as the initial values and the local optimization of parameters is achieved by the Newton method (Mallat & Zhang 1993). (4) Calculate the amplitude of wavelet a n . After determining the fundamental parameters of the wavelet, the amplitude can be calculated based on equation (15).
Going through these four steps, one iteration of MPFR will be completed. By repeating these steps until the maximum number of iterations or residual energy reaches a predefined level, signal decomposition can be accomplished. Based on this, WVD of each decomposed wavelet is computed and superimposed, and the time-frequency distributions of the original seismic data are obtained.  Compared with the Ricker wavelet dictionary, the Morlet wavelet dictionary has more basic parameters. Therefore, taking a Morlet wavelet for instance, the MPFR is demonstrated and compared with the common time-frequency analysis methods. Figure 6a is a synthetic trace, the constituent wavelets and parameters are shown in figure 6b and Table 1, respectively, where the wavelets at 0.7 and 0.9 s are superimposed. Figure 6c shows the WVD superposition of wavelets in figure 6b, where the resolution of time-frequency spectrum is quite high. Besides the first, fifth and seventh wavelets, the spectra of other wavelets rotate to different degrees. Figure 6  of figure 6a calculated by STFT, CWT and ST, respectively. It is demonstrated that the time-frequency resolution of STFT is the lowest. Yet, in comparison with figure 6c, energy aggregation still exists in the results of CWT and ST where the wavelets at 0.7 and 0.9 s are overlapped and difficult to distinguish. However, it is gratifying to notice that the rotated time-frequency characteristics of wavelets can be expressed by these three transforms from the spectrum at 0.2, 0.35 and 0.5 s. The reason is that these methods actually belong to the average fitting process of local signals under the continuous translation of the time window. To better describe the local signals, the kernel function in the transformation method keeps changing as the time-window position changes. Taking the wavelet at 0.2 s as an example, when the center of the window locates in the upper half of the waveform, the frequency in the time window is higher and the corresponding kernel function also belongs to the high-frequency characteristic. When the window moves downward, the frequency of the waveform in the time window decreases gradually and the frequency characteristic of the kernel function also decreases correspondingly. As a result, the frequency decreases with time in the time-frequency spectrum. Figure 6g is the time-frequency spectrum of figure 6a obtained using MPFR, where the superimposed wavelets at 0.7 and 0.9 s are well separated and the rotated time-frequency feature is also accurately expressed. Figure 6h is the decomposed wavelets cor-responding to figure 6g. Only the results of the first 10 iterations are shown here. The wavelets in the first eight channels are consistent with those in figure 6b.

Description of rotated time-frequency characteristics of attenuated seismic waves
To describe the rotated time-frequency characteristics of attenuated seismic waves, the Morlet and Ricker wavelets at 1.9 s corresponding to the Q = 50 in figure 2 are adopted to compute the wavelet decomposition and time-frequency spectrum based on MPFR, which is compared with the conventional MP method. Targeting the attenuated Morlet wavelet, figures 7 and 8 show the results of MP based on the rotated time-frequency atomic dictionary (MPFR) and the conventional atomic dictionary (MP), respectively. In both cases, 10 iterations are used. The blue lines in figures 7a and 8a are the original attenuated signal, and the red lines are the superposition of the decomposed wavelets. Both of them are in good agreement with the original signal. Figures 7d and 8d are curves with error display, which show that the error levels are 2 × 10 −5 and 5 × 10 −4 , respectively. Figures 7b and 8b  the attenuation rate of wavelet amplitude by MPFR is faster than that by the conventional MP, which indicates that the rotated time-frequency atomic dictionary is more competent to represent attenuated signals than the conventional atomic dictionary. Figures 7c and 8c are the time-frequency spectra associated to the two methods. The long or short axis of the spectrum by the conventional MP method is parallel to the coordinate axis, which cannot express the rotated timefrequency features evoked by the attenuation. However, the long axis direction of the time-frequency spectrum obtained using MPFR has a certain angle to the time or frequency axis, which effectively shows the rotated time-frequency characteristics caused by the attenuation. Similar to figures 7 and 8, figures 9 and 10 show the results of 10 iterations of MPFR and conventional MP methods targeting an attenuated Ricker wavelet, respectively. From the curves with error display (figures 9d and 10d) and the wavelet amplitudes (figures 9e and f, figures 10e and f), it is demonstrated that after the same number of iterations, the error between the reconstructed and the attenuated signals of MPFR is less than that of the conventional MP. Meanwhile, the rotated time-frequency characterization of the attenuated Ricker wavelet can be observed from the spectrum of MPFR (figure 9c).
To further evaluate the MPFR method, a velocity model is constructed in figure 11a. The Q model (figure 11b) is ob-tained using the empirical formula about the Q value and velocity (Li 1994). Figures 12a and 13a simulate the noise-free and noise-added traces with a noise level of 5%, respectively. In the noise-added trace, the noise is filtered into the seismic frequency band to make the simulation more realistic. The time-frequency spectra of noise-free and noise-added traces are shown in figure 12b, c, d and e and figure 13b, c, d, e, respectively. The results indicate that STFT, CWT, ST and MPFR methods can describe the rotated time-frequency characteristics of the attenuated seismic trace and have certain anti-noise ability. However, MPFR method has a higher time-frequency resolution than the other three methods.

Seismic data reconstruction
Reconstruction of seismic data is an important application of MP (Wang, 2010). Seismic traces are decomposed into the superposition of a series of constituent wavelets. The smaller the number of constituent wavelets, the more the wavelets can represent the basic characteristics of seismic traces. Figure 14 shows a portion of seismic profile, which is reconstructed by the conventional MP and MPFR methods. Figure 15 shows the total residual energy illustrated in the normal and the decibel values when the reconstructed 835 Figure 14. A field seismic profile.
wavelets of every trace gradually increases from 1 to 500. The blue dotted line denotes the conventional MP method and the red circled line denotes the MPFR method. The reconstruction errors of both methods decrease with the increasing wavelet number. The residual energy of MPFR method is always lower than that of the conventional MP method with the same number of reconstructed wavelets. Figures 16a, d and g show the results of seismic data reconstruction using the conventional MP method (blue) and the MPFR method (red) when the number of decomposed wavelets in each trace is 10, 20 and 40, respectively. To compare the difference in reconstruction results between the two methods more clearly, the seismic traces were extracted with an interval of 50, and the two results were displayed together with the original seismic traces (black) in a way of trace comparison. When the number of reconstructed wavelets is small, only part of the waveforms in the original seismic trace can be restored by these two methods. In particular, there is a large error in the interval 2.8-3.1 s. With the increase of the wavelet number, the reconstructed results are gradually close to the original seismic trace and the mismatch between 2.8 s and 3.1 s is also improved. Figure 16b shows the residual section of the two methods corresponding to figure 16a, which adopts the same scale to analyze the results distinctly. Likewise, figure 16 parts e and h are residuals corresponding to figure 16d and g. Figure 16 parts c, f and i, respectively, show the residual energy values of the reconstructed traces relative to the original traces by the two methods. With the increase of the wavelet number, the relative residual energy shows a decreasing trend, and under the same number of wavelet reconstruction, the relative residual energy of each trace based on MPFR method is smaller than that based on conventional MP method, which indicates the advantages of the rotated time-frequency atomic dictionary in seismic data reconstruction.

Time-frequency analysis of a gas-bearing reservoir
The spectral decomposition technique has been widely used in seismic exploration (Partyka et al. 1999;Wang 2007). By extracting several common frequency profiles, formation structure or reservoir information can be displayed to carry out reservoir prediction (Castagna & Sun 2006). Figure 17a shows a section of a seismic profile, the region marked by the ellipse has been confirmed as a gas-bearing reservoir. Based on the spectral decomposition, the common frequency profiles of 10, 20, and 30 Hz are extracted in figure 17 parts b, c and d, respectively. One can note that the reservoir response is more obvious in the low-frequency section of 10 Hz attributed to the absorption and attenuation effect of the reservoir. As the frequency moves toward higher values, the response of the strata gradually enhances and becomes confused with the gas-bearing region. Therefore, the location of gas-bearing reservoir can be effectively predicted in the low-frequency profile. Based on this, the time-frequency spectra of gas-bearing and non-gas-bearing regions are extracted by means of MPFR (figure 18) and conventional MP ( figure 19). For the 200th and 700th traces in the non-gas region, the rotated time-frequency characteristics are both not obvious by the two methods. However, for the 400th and 500th traces in the gas region, the spectra in figure 18 show a certain characteristic of time-FR, which is invisible in figure 19. Therefore, MPFR technology can be used as an auxiliary technology, combined with spectral decomposition, to achieve the prediction of reservoirs with attenuated characteristics (such as gas-bearing reservoirs).

Discussion
In real applications, how to understand the results of different time-frequency analysis methods is a problem worthy of attention (Castagna et al. 2003). Figure 20 shows synthetic data from the wedge model with a Morlet wavelet, whose dominant frequency is 30 Hz. The WVD superposition of constituent wavelets (WVDS), STFT, conventional MP and MPFR of each seismic trace are calculated, respectively. Based on this, three representative time-frequency spectra are shown in figure 20b.
(1) When the thickness of the stratum is thicker than the tuning thickness (trace number = 70), the two wavelets are highly discernible in the time axis and the spectrum is well separated on the four sub-graphs, just the STFT has a lower resolution.   the seismic trace corresponds to two waveforms, the time-frequency spectra of STFT, MP and MPFR are not well separated since they are whole-trace calculations. Especially for MP and MPFR methods, it is impossible to distinguish the two wavelets effectively. Instead, the overlapped waveforms are treated as a single wavelet and a single spectrum is obtained. In addition, the time-frequency spectrum of MPFR is rotated, which is caused by wavelet interference and is independent of dispersion and attenuation. So, in this situation, MPFR could produce misleading results. (3) When the thickness of the stratum is less than the tuning thickness (trace number = 5), the two wavelets are intensively interfered with. The time-frequency spectrum obtained by the four methods all shows one energy group, so it is impossible to distinguish the two wavelets. However, the frequency center in the spec-trum of WVDS corresponds to 30 Hz, which is the same as the dominant frequency of the wavelet. While the frequency center of STFT, MP and MPFR deviates from 30 Hz, this is because the latter three methods take the overlapped wavelets as a single wavelet. The tuning effect changes the dominant frequency, leading to the variation in the frequency center of the spectrum.
Therefore, different time-frequency analysis methods contribute to varying results for the same seismic data. For the synthetics of forward modeling, the real time-frequency spectrum, namely WVDS, can be obtained as the wavelet is known. However, it is hard to get the true spectrum for the real data. In this case, Castagna et al. (2003) pointed out that the key issue is not to distinguish right or wrong of the time-frequency analysis, but whether the results can capture the features important in target interpretation. 839 Therefore, the practical problems should be combined with specific analysis. For example, for the trace number = 5 in figure 20b, although the frequency center calculated by the last three methods deviates from the dominant frequency of the wavelet, it is also a true hint of the frequency when the two wavelets overlap to form a single seismic waveform. In practice, it is necessary to distinguish whether the frequency variation is a response of a thin layer. For the MPFR calculation of trace number = 25 in figure 20b, we should analyze the reasons of time-FR. If there is other evidence to indicate the existence of attenuation (such as gas-bearing reservoirs), to some extent, it can be considered that the attenuation causes the time-FR. Otherwise, it is more likely to be the result of wavelet interference.

Conclusion
Absorption and attenuation affect seismic wave propagation, resulting in the rotated time-frequency spectrum. Traditional time-frequency analysis methods, such as STFT, CWT and ST, can describe the rotated characteristics, but they are limited by resolution and cannot meet the requirements of subtle interpretation. The conventional MP based on a Morlet or Ricker wavelet dictionary can improve the resolution of the time-frequency spectrum, but cannot effectively describe the rotated time-frequency characteristics of attenuated seismic waves. The MPFR method not only guarantees the timefrequency resolution, but also effectively solves the characterization problem of the rotated time-frequency by introducing the FR factor. Real data applications show that the rotated time-frequency atomic dictionary are more effective in characterizing non-stationary seismic data. Meanwhile, the description of rotated time-frequency characteristics can be used as an auxiliary technique to predict the attenuated reservoirs, such as gas-bearing reservoirs, by combining with spectral decomposition.