Time-frequency analysis maps an 1D time signal into 2D time and frequency domains, which can capture the non-stationary character of seismic data. Time-frequency analysis is a fundamental tool for seismic data analysis and geological interpretation (Liu et al., 2016; Reine et al., 2009; Chen et al., 2014; Castagna et al., 2003). Conventional time-frequency methods, such as short time Fourier transform (Cohen, 1989), wavelet transform (Mallat, 1989) and S-transform (Stockwell et al., 1996) are under the control of Heisenberg/Gabor uncertainty principle, which states that we cannot have the energy arbitrarily located in both time and frequency domains (Mallat, 2009). Moreover, short time Fourier transform, wavelet transform and S-transform are using a windowing process, which often brings smearing and leakage (Tary et al., 2014). Therefore spurious frequencies are often generated, which will "color" the real time-frequency map and affects the interpretation. In recent years, many new methods were proposed such as matching pursuit (Mallat and Zhang, 1993), basis pursuit (Chen et al., 1998), empirical mode decomposition (Huang et al., 1998; Chen and Fomel, 2015), the synchrosqueezing wavelet transform (Daubechies et al., 2011) and its variants such as the synchrosqueezing short time Fourier transform (Oberlin et al., 2014), the synchrosqueezing S-transform (Huang et al., 2015). The matching pursuit and basis pursuit methods represent the energies of the input signal by the energies of atoms found using different methods, which prevents smearing and leakage, creating highly localized time-frequency decompositions. The efficiency of these two methods depend on the predefined dictionary (Tary et al., 2014). The empirical mode decomposition method decomposes a signal into symmetric, narrow-band waveforms named intrinsic mode functions to compress artificial spectra caused by sudden changes and to therefore improve the time-frequency resolution (Han and van der Baan, 2013). However, the empirical mode decomposition also suffers from mode mixing and splitting problems. In order to solve the above problems, alternative methods were developed based on empirical mode decomposition: ensemble empirical mode decomposition(Wu and Huang, 2009), complete ensemble empirical mode decomposition (Torres et al., 2011). However, these two methods, like the empirical mode decomposition, are still "empirical" because their sketchy mathematical justifications. The synchrosqueezing wavelet transform (Daubechies et al., 2011) and its variants capture the philosophy of empirical mode decomposition, but use a different method to compute the intrinsic mode functions providing a rigorous mathematical framework.

Similar to the Fourier transform, the Prony method (Prony, 1795) decomposes a signal into a series of damped exponentials or sinusoids in a data-driven manner, which allows for the estimation of frequencies, amplitudes, phases and damping components of a signal. Fomel (2013) proposed the non-stationary Prony method (NPM) based on regularized non-stationary auto-regression. The NPM decomposes a signal into intrinsic mode functions with controlled smoothness of amplitudes and frequencies like the empirical mode decomposition does, but uses NPM instead. Unlike Fourier transform, the coefficients of the extracted intrinsic mode functions for the Prony method do not clearly define the energy distribution for the input signal in the time-frequency domain. Therefore, the NPM used by Fomel does not clearly define a "real" time-frequency map but a "time-component" map. In this paper, we couple the NPM (Fomel, 2013) and the Hilbert transform to give a time-frequency decomposition. The proposed method has a rigorous mathematical framework. Furthermore, synthetic and real data tests confirm that the intrinsic mode functions derived by the proposed method are more smooth with respect to the amplitudes and frequencies than the intrinsic mode functions of ensemble empirical mode decomposition (Wu and Huang, 2009). Synthetic and real data tests also confirm that the proposed method has a higher time-frequency resolution than the ensemble empirical mode decomposition. The proposed method can be used to facilitate seismic interpretation.