Inverse spectral decomposition using an lp-norm constraint for the detection of close geological anomalies

An important application of spectral decomposition (SD) is to identify subsurface geological anomalies such as channels and karst caves, which may be buried in full-band seismic data. However, the classical SD methods including the wavelet transform (WT) are often limited by relatively low time–frequency resolution, which is responsible for false high horizon-associated space resolution probably indicating more geological structures, especially when close geological anomalies exist. To address this issue, we impose a constraint of minimizing an lp (0 < p < 1) norm of time–frequency spectral coefficients on the misfit derived by using the inverse WT and apply the generalized iterated shrinkage algorithm to invert for the optimal coefficients. Compared with the WT and inverse SD (ISD) using a typical l1-norm constraint, the modified ISD (MISD) using an lp-norm constraint can yield a more compact spectrum contributing to detect the distributions of close geological features. We design a 3D synthetic dataset involving frequency-close thin geological anomalies and the other 3D non-stationary dataset involving time-close anomalies to demonstrate the effectiveness of MISD. The application of 4D spectrum on a 3D real dataset with an area of approximately 230 km2 illustrates its potential for detecting deep channels and the karst slope fracture zone.


Introduction
Detecting anomalous bodies in the subsurface, such as channels, karst caves, faults and alluvial fans, is one of the key objectives in seismic exploration (e.g., Bacon et al. 2003;Qi et al. 2014;Kang et al. 2020). These geological anomalies will be conductive to the construction of traps and reservoirs in the Earth, since they can usually provide storage space and migration channel for oil and natural gas (e.g., Hart et al. 2002;Ben-Zion 2008;Mittempergher et al. 2009). Seismic data acquired by the dense sensors placed at the surface can carry the information associated with such geological anomalies. However, due to the attenuation related to near surface and fluids, the frequency band of seismic data is reduced from shallow to deep, and seismic waveforms are commonly interfered together. Therefore, it is sometimes difficult to directly use full-band and non-stationary seismic data to identify the subsurface geological anomalies. Decomposing seismic data into certain modes with different dominant frequencies and bandwidth can facilitate extracting intrinsic spectral features (e.g., Li et al. 2017;Liu et al. 2019).
Spectral decomposition (SD) is one of the widely used techniques for signal decomposition. SD that decomposes a seismic signal into a 2D function of time and frequency is an effective method to characterize full-band and nonstationary signals (Partyka et al. 1999). It is often useful for evaluating the frequency-dependent response of seismic data associated with layered earth structures. Moreover, spectral anomalies at different frequency components (e.g., at low frequencies and high frequencies) can illuminate geological anomalous bodies of different scales (e.g., Laughlin et al. 2002;Li et al. 2011) and indicate hydrocarbon (such as oil, gas or gas hydrate) in sandstone or carbonate reservoirs (e.g., Castagna et al. 2003;Li et al. 2012). Therefore, SD has been routinely used in thin-bed prediction (Marfurt and Kirlin 2001), reservoir characterization (Li et al. 2011), and hydrocarbon detection (Huang et al. 2016;Naseer and Asim 2018). As well, the amplitude-frequency spectral components, phase-frequency spectral components or both amplitude-frequency and phase-frequency spectral ones are widely adopted to calculate frequency division coherence (e.g., Li and Lu 2014) for detecting anomalous bodies clearer with more geological boundary details, in contrast to the broadband coherence.
In the last few decades, various SD methods have been proposed and widely developed for analyzing a great variety of signals, such as the short-time Fourier transform (STFT) (Gabor 1946), the continuous wavelet transform (CWT) (Sinha et al. 2005), and S-transform (ST) (Stockwell et al. 1996) with its generalized variations (Pinnegar and Mansinha 2003a, b;Liu et al. 2018Liu et al. , 2019a. The classical SD methods are simple and fast to be implemented. However, they suffer from a trade-off between the time resolution and the frequency resolution, due to the effect of the window function and the limitation of the Heisenberg uncertainty principle (Heisenberg 1927). In special, seismic signals interfere together when two or several geological bodies are close, which will also affect the frequency components of local signals, especially for the deep. In order to yield high time resolution, the isofrequency amplitude slice obtained by using the conventional SD methods readily suffers from strong spectral energy mixture issue across the frequency band or the spectral leakage phenomenon (Oyem and Castagna 2013). Consequently, the isofrequency amplitude slice is still the superposition of different frequency components, which is not easy to interpret close geological bodies. In addition, the multitude of these methods signifies the non-unique nature of spatiotemporal transformation (Castagna et al. 2003). The inverse spectral decomposition (ISD) method by solving an inverse problem with a sparse spike constraint (Portniaguine and Castagna 2004) has the potential to alleviate the above these issues. Without the constraint, ISD is equivalent to the conventional SD. After the ISD is proposed, its improvements are focused primarily on the appropriate selection of the wavelet library (such as the truncated sinusoid, the Ricker wavelet, the Morlet wavelet or the extracted wavelet), the type of the constraint or the prior information (such as the l 2 norm, the l 1 norm, the mixed l 1 -l 2 norm, the coherency-based constraint, or the hierarchical Gaussian prior), and the algorithm (such as the iterative soft thresholding algorithm, the iteratively reweighted least-squares algorithm, spectral projected gradient for l 1 minimization, Bregman iterative algorithm or sparse Bayesian learning) for solving ISD (Puryear et al. 2012;Han et al. 2012;Gholami 2013;Tary et al. 2014;Amosu et al. 2016;Ma et al. 2019;Yuan et al. 2019). In contrast to the development of methods, the fine application of the ISD results is rare (Gholami 2013;Oyem and Castagna 2013;Han et al. 2016;Li et al. 2016;Wang et al. 2018), especially for interpreting 3D seismic data more than tens of thousands of traces, which are very common in oil companies. Moreover, how to make full use of the ISD results is not unimportant.
In this paper, we investigate the potential application of a modified ISD method (named MISD for short) with an l p -norm (0 < p < 1) constraint to detect close geological anomalies along the time direction or those at the frequency along a horizon. Similar to the typical widely adopted l 1 -norm constraint, the l p -norm regularization with the proven convergence (Raskutti et al. 2011) can also promote the sparsity of the resulting time-frequency spectrum (Yuan et al. 2020). However, the l p norm has the better sparse representation than the l 1 norm, and overcomes the disadvantage that larger coefficients are penalized more heavily in the l 1 -norm than smaller coefficients to some extents (Candes et al. 2008). Furthermore, the sparsity of the l p -norm (1/2 ≤ p < 1) solution increases as p decreases, whereas the sparsity of the solution for l p norm (0 < p ≤1/2) does not overly change with respect to p (Fan and Peng 2004). Xu et al. (2010) have shown that the l 1/2 norm is an unbiased estimator which imposes strong sparsity upon the minimization problem at hand. Based on these advantages, the l p norm is applied to yield a more compact or focusing time-frequency spectrum for broadband and non-stationary seismic signal analysis, and can separate the frequency components of the wavelet interference and reduce spectrum leakage, which contribute to detect the distributions of close geological features from high-resolution spectral slices of horizons.
We begin this paper with the methodology of MISD. Then, a 3D synthetic seismic dataset involving frequencyclose geological anomalies along the horizon and the other 3D non-stationary synthetic dataset complicated by both the time-variant wavelets and thin-bed interference are adopted to demonstrate the effectiveness of MISD in detecting close channels and alluvial fans with different spatial distributions, and to illustrate its advantages over CWT and ISD using a typical l 1 norm. As well, we extracted the seismic traces from 3D synthetic data to further explain why MISD has the highest time-frequency resolution. At last, the application of 4D high-resolution time-frequency spectrum on a 3D field dataset from the Northwest China illustrate that MISD can be used to interpret deep channels and the karst slope fracture zone at two close frequencies, where each frequency presents the minimum number of geological structures than CWT and ISD using a typical l 1 norm.

3 2 Methodology
The inverse wavelet transform of a signal s(t) can usually be expressed as: where a is the scale factor (that can be converted to frequency), τ is the translation factor, C(τ, a) is the spectral coefficient of wavelet transform, ψ((t-τ)/a) is the wavelet library composed by the mother wavelet through stretching and translation, and t is time. In seismic exploration, the mother wavelet can usually be the Morlet wavelet, Ricker wavelet, Haar wavelet alone or in combination with characterize the recorded seismic data. As studies found, Ricker wavelets can be adopted to better fit seismic data that tend to be biased toward the lower frequencies (Ricker 1953;Liu and Marfurt 2007). Moreover, Ricker wavelets have been widely applied to seismic data modeling, processing, inversion and interpretation. Consequently, we choose the Ricker wavelet as the mother wavelet to construct the wavelet library in this paper.
As well, the integration of the translation factor τ in Eq. (1) can be derived to express as a convolution form: where * represents the convolution operator.
It is well known that the integral can be solved discretely, and a general expression form is obtained as follows: where f k (k = 1, 2…,K) represents the frequency corresponding to the scale factor a in Eq. (2), K represents the frequency sample number, w(t, f k ) corresponding to a −5/2 ψ(t/a) in Eq. (2) represents the frequency-dependent wavelet or basis function, r(t, f k ) corresponding to C(t, a) represents spectral coefficients (or the time-frequency pseudo reflectivity), and s k (t, f k ) represents the decomposed stationary signal. According to Eq. (3), a non-stationary signal can be described by using the superposition of the convolution results of the known wavelets with different dominant frequencies and the spectral coefficients corresponding to the dominant frequency. In this way, the inverse CWT of a nonstationary signal can be regarded as a process of quadratic linear superposition. Because the convolution operator in the time domain is equal to a product operation in the frequency domain, Eq. (3) can be implemented fast by using the Fourier transform. Besides the non-stationarity caused by the near surface and fluids, thin-bed interference can lead to the data non-stationarity (Yuan et al. 2017). In turn, we can decompose the non-stationary seismic data resulting from the interference to obtain the time-frequency pseudo reflectivity corresponding to a series of wavelets of varying frequencies. It provides a basis for interpreting subsurface close geological abnormalities from non-stationary data.
To facilitate the description of the subsequent inverse problem, Eq. (3) is written in a matrix-vector formulation: represents the wavelet dictionary, W(t, f k ) refers to the convolution matrix of the frequency-dependent wavelet [w(t 1 , represents the spectral coefficients corresponding to a certain dominant frequency f k , N represents the time sample number, [·] T represents the transposition, and the vector s represents a vertical trace of non-stationary signal. In this way, a non-stationary signal can be described by using the product of the wavelet dictionary matrix with a large column number and a long spectral coefficient column vector. Similarly, STFT and ST (Stockwell et al. 1996) can also be expressed in the form of Eq. (4), but only W(t, f k ) is different. In theory, the result of the conventional SD can be obtained by solving m in Eq. (4) in a least-squares method. Nevertheless, Eq. (4) is ill-posed since the row number of the matrix G is usually far smaller than its column number and the frequency-dependent wavelets are bandlimited, which lead to the solution of Eq. (4) not unique. It probably provides another insight why can the conventional SD obtain multiple time-frequency spectra. As expected, the general purpose of SD is to obtain high-resolution or more compact time-frequency spectrum. The physical meaning of the time-frequency focusing property is strongly relevant to sparsity. The mathematical expression of the physical meaning can be approximated in one non-convex l p -norm minimization (Chartrand 2007) as follows: , and the p value is limited between 0 and 1. Theoretically, minimizing the l 0 norm (i.e., the number of nonzero elements) can yield the sparsest solution (Candes et al. 2008). This is of little practical use, however, since minimizing the l 0 norm is a non-deterministic polynomial (NP) hard problem. Mathematically, all situations are required to exhaust to find the optimal solution. Especially for the large-scale problem, it is almost infeasible to be solved (Xu et al. 2010). It is intuitive that when p is arg min ‖ ‖ p , close to 0, the solution to the l p (0 < p<1) norm will be close to that of the l 0 norm. Moreover, the l p norm penalizes large coefficients smaller, while penalizes small coefficients more heavily, in contrast to the typical l 1 norm ‖ ‖ 1 = ∑ K×N j=1 � j � . Therefore, the l p norm can produce a sparser solution than the typical l 1 norm (Chartrand 2007). On the whole, under the premise of satisfying the conventional SD, we hope the time-frequency spectrum is sparse, that is, to say the resolution of both time-frequency amplitude spectrum and phase spectrum is high. Consequently, Eq. (5) is imposed on Eq. (4) to obtain an optimal sparse time-frequency spectrum by solving the following function: In general, the above Eq. (6) can be solved by transforming it into an l p (or hyper-Laplacians) regularization problem as follows: where ‖ − ‖ 2 2 represents the non-stationary data misfit, and λ represents a regularization parameter balancing the sparsity desired in the resulting time-frequency spectrum and the non-stationary data misfit. If p = 1, Eq. (7) is equivalent to the conventional ISD in the lasso (Tibshirani 2011) form. Therefore, the conventional ISD can be regarded as a special case of MISD in this way.
Adopting the basic gradient idea to the above equation leads to a following iterative scheme: where β > 0 is an appropriate step size that depends on the maximum eigenvalue of G H G or the Lipschitz constant (Beck and Teboulle 2009), m u represents the solution at the uth iteration, and [·] H represents the complex conjugate.
Equation (8) can be effectively solved by using a generalized iterated shrinkage algorithm (Zuo et al. 2013), where each iteration involves relatively cheap matrix-vector multiplication followed by a generalized shrinkage/soft-threshold step and a fast iterative strategy (Beck and Teboulle 2009). The main iteration expressions are as follows: where T GST p is the generalized shrinkage operator (Zuo et al. 2013), y u is a clever update of the model estimate considering a very specific linear combination of the two previous updates that can play an important role in accelerating convergence (Beck and Teboulle 2009), γ is usually initialized to 1, and τ p is the threshold related to both λ and p values (She 2009). When p = 1, τ p = λ, i.e., the typical threshold corresponds to the l 1 norm. The matrix-vector multiplication Gy u can be regarded as the forward modeling from the highdimension time-frequency domain to the low-dimension data domain, whereas G H (s − Gy u ) can be understood as the adjoint operator of data residual (i.e., the difference between the observed non-stationary data and the calculated data) from the data domain to the time-frequency domain. These both matrix-vector multiplications are similar to Eqs. (3) and (4), and thus can be implemented quickly by using the fast Fourier transform.
After iteratively solving the optimal m, we can obtain the corresponding time-frequency spectrum. Although both the sparsity of time-frequency spectrum and how to get the sparse time-frequency spectrum are meaningful, it is not unimportant to make full use of the sparse time-frequency spectrum. Moreover, it is helpful for seismic interpretation that the sparse time-frequency spectrum is applied reasonably. In this paper, we use the 4D sparse time-frequency spectrum to identify close geological anomalies at the frequency along a horizon or those along the time direction. We also emphasize that the horizon magnitude spectra at the predominant frequencies can present the less number of geological structures due to the sparsity constraint, and thus improving the spatial resolution of horizons and reducing false horizon-associated geological structures.

Examples
In this section, a 3D synthetic dataset with an area of 100 km 2 involving single-stage thin sand bodies of three close dominant frequencies as well as the other 3D nonstationary synthetic dataset with the same area involving two-stage close sand bodies along the time direction are designed to test the performances of MISD. These two synthetic examples are adopted to focus on illustrating that the 4D spectral results obtained by using MISD can favorably identify close geological bodies with different spatial distributions, as well as showing its advantages over CWT and conventional ISD. Furthermore, a 3D real dataset with an area of approximately 230 km 2 from the Northwest China is used to test its application potential for detecting channels and the karst slope fracture zone in the deep target reservoir. It will be emphasized for these three 3D data examples that MISD-based spectral results of horizons present the minimum number of geological structures than CWT-based and ISD-based those, and thus making seismic structural interpretation clear. For the three examples, MISD adopts the same p value of 0.5 and the total iteration number of 100. For each trial, ISD and MISD adopt the same algorithm flow and the wavelet library except for different p values.

3D synthetic data example
In order to clarify the ability of MISD to identify geological bodies with different spatial distributions and close dominant frequencies, we design a single-stage thin sand layer with a time range of 160 ms to 170 ms, which consists of one meandering channel complex and two alluvial fan structures. Then, a 3D seismic data cube with a size of 201 Inlines × 201 Crosslines × 300 time samples is synthesized by using three zero-phase Ricker wavelets with dominant frequencies of 25 Hz, 30 Hz and 35 Hz for three sand bodies from Northwest to Southeast, respectively. For the model, the space intervals along both Inline (North-South) and Crossline (East-West) directions are 50 m. Figure 1 shows the comparisons among CWT-, ISD-and MISD-based spectral amplitude slices at three close frequencies and 165 ms. It can be seen from Fig. 1a-c that no matter in 25-Hz, 30-Hz or 35-Hz CWT-based spectral amplitude slices, both channel complex and alluvial fan structures exist. In other words, CWT cannot distinguish these geological anomalies located at different spatial locations by using spectral amplitude attributes. It suggests that the frequency focus of CWT is poor at both relatively low frequency and relatively high frequency, as well as the frequency leakage is serious. It can be observed that the 25-Hz ISD-based spectral amplitude slice (Fig. 1d) shows only the Northwest alluvial fan, but both 30-Hz and 35-Hz ISD-based spectral amplitude slices (Fig. 1e, f) indicate the same channel complex and Southeast alluvial fan together. It means that although the frequency resolution of ISD has been improved at the relatively low frequency, it is low at middle and high Only MISD can indicate three 10-ms time thickness sand bodies with the true spatial distribution characteristics separately by using spectral amplitude slices at three close frequencies frequencies. As expected, MISD can image all sand bodies with the true spatial distribution characteristics separately by using spectral amplitude slices at three corresponding close dominant frequencies, as shown in Fig. 1g-i. The above comparisons demonstrate that the sparsity constraint represented by using an l p norm with a p value of 0.5 plays a role in obtaining accurate spectral segmentation, so the MISD-based spectral amplitude attribute gets rid of the frequency mixture issue. In order to test the ability of MISD to identify geological bodies with different spatial distributions and close time locations, we design a 3D synthetic data set with a size of 201 Inlines × 201 Crosslines × 300 time samples. The shallow upper layer with the top interface of 100 ms and the bottom interface of 120 ms, mainly consists of one meandering channel complex and two alluvial fan structures. The deep lower layer consists of one distributary channel complex and one alluvial fan structure, with a time range of 160 ms to 180 ms. Then, a 3D seismic data cube with frequency attenuation is synthesized by using two zero-phase Ricker wavelets with dominant frequencies of 40 Hz and 20 Hz for the shallow single-stage sand bodies and the deep single-stage sand bodies. Figure 2a, b shows the lateral distributions of the upper-layer and lower-layer sand bodies with the space intervals of 50 m along both Inline and Crossline directions. It can be observed that all five sand bodies share different lateral distributions. Figure 2c is the reflection coefficient corresponding to one common depth point (CDP) involving both upper-layer and lower-layer sand bodies.
It is difficult to distinguish the two-stage sand bodies by using the original full-band amplitude slices, which are complicated by both the time-variant wavelets and thin-bed interference related to the shallow single-stage sand bodies and the deep single-stage sand bodies. Time-frequency analysis is an effective way to characterize full-band nonstationary signals by decomposing a series of frequencies.
Figure 3a-f shows 40-Hz isofrequency magnitude slices at 125 ms and 131 ms, which are obtained by using CWT, ISD and MISD. It can be seen that CWT (Fig. 3a, b) has the strongest interference pattern with the imprints of all five sand bodies overlapped together. Although ISD can image three shallow sand bodies, there are some shadows of the deep sand bodies (indicated by the red arrows in Fig. 3c), and a part of the deep sand bodies which is overlapped with Fig. 3c can also be seen in Fig. 3d. It can be found that the overlapping areas correspond to the deep geological bodies of 20 Hz. In other words, 20-Hz deep single-stage sand bodies appear in 40-Hz isofrequency magnitude slices of shallow horizons, since the interference changes the frequency components. Therefore, the isofrequency magnitude slice is still the superposition of different frequency components, which is not easy to interpret relatively close geological bodies. That is to say, the 40-Hz ISD-based isofrequency magnitude slice does not only contain the geological bodies of 40 Hz, but also the geological structures of other frequencies. However, the MISD-based spectral amplitude slice at 125 ms ( Fig. 3e) indicates only three sand bodies in the upper layer, which reflects the real distribution characteristics of the upper layer. Moreover, there are no obvious imprints of 20-Hz deep single-stage sand bodies even at a relatively deep slice, as shown in Fig. 3f. Figure 4 shows the comparisons among the 20-Hz spectral amplitude slices obtained by using different spectral decomposition methods. The CWT-based results (Fig. 4a,  b) show that the two layers interfere together and cannot be identified separately. Although the ISD-based spectral amplitude attribute at a relatively shallow horizon slice does not image any part of 20-Hz deep sand bodies, 40-Hz shallow sand bodies are clearly appear in the shallow spectral amplitude slice, as shown in Fig. 4c. Moreover, the deep ISD-based spectral amplitude slice visually presents geometry of two 20-Hz deep sand bodies well,  Fig. 4d. It is mainly because the overlapping parts interfere with each other to change the frequency components. Nevertheless, as expected, MISD does not image any structures of two-stage sand bodies in the shallow spectral amplitude slice (Fig. 4e), while can clarify the distributions of two deep sand bodies perfectly in the deep spectral amplitude slice (Fig. 4f). The results in both Figs. 3 and 4 suggest that MISD can yield a 4D spectrum with higher time-frequency resolution than CWT and ISD, To make an intuitive insight that MISD has a higher time-frequency resolution, we extract three seismic traces from the 3D non-stationary synthetic data, which are related to only one upper geological body, only one lower geological body, and both upper and lower geological bodies, as shown in Figs. 5a, 6a and 7a. The other subfigures in Figs. 5, 6 and 7 display the time-frequency magnitude spectra obtained by using CWT, ISD and MISD, respectively. It can be seen that the first vertical white line on the left (20 Hz) of CWT (Fig. 5b) passes through the 40-Hz shallow energy cluster. Although ISD is superior to the CWT method (Fig. 5c), the first vertical white line passes through the boundaries of the 40-Hz energy cluster. It can be clearly observed that the left white line of MISD does not pass through the 40-Hz energy cluster, as indicated in Fig. 5d. In Fig. 6, the observed phenomena for the right white lines (40 Hz) are similar to Fig. 5. The right vertical white line of MISD not only does not pass through the 20-Hz deep energy cluster, but also is far away, suggesting that MISD has a higher frequency resolution. Figures 5 and 6 show that both CWT and ISD suffer from strong frequency mixture issue, because of the poor time-frequency focusing property. One can see from Fig. 7b-d that the upper and lower layers of CWT and ISD are interfered together, and CWT is the most serious, while the upper and lower layers of MISD are separate. The upper and lower layers interfered together will affect the frequency components of local signals, especially for the deep. To further explain that MISD has a higher time resolution, we extract 40-Hz frequency components as denoted in the right white lines of Fig. 7b-d, as shown in Fig. 8a-c. It is obvious that the leakage at 40-Hz frequency components obtained via CWT is the most serious. Although the time-frequency resolution of ISD has been improved, there is still leakage, which affects the structure interpretation. It is why a little leak of ISD has an effect on the spectral amplitude slice that cannot be ignored. However, MISD present only 40-Hz shallow sand-layer associated spectral components, meaning that it does not suffer from the frequency leakage. Moreover, the spectral magnitudes of both CWT and ISD are not accurate and focused, as indicated by the black ellipses in Fig. 8a, b, whereas the MISD-based result (Fig. 8c) is as perfect as expected. Similarly, we extract 20-Hz frequency components from Fig. 7b-d, as shown in Fig. 9a-c. The results corresponding to CWT and ISD still have leakage, but MISD can only show the 20-Hz deep sandlayer associated spectral components. Moreover, the width (the red arrow) is the narrowest in contrast to CWT and ISD, indicating that MISD has a highest time resolution. In conclusion, MISD has a highest time-frequency resolution and the best energy concentration, which provides a potential accurate way for seismic data interpreters to analyze seismic signals.

3D field data example
A 3D migrated field dataset from the Northwest China with a size of 594 Inlines × 594 Crosslines × 251 time samples is utilized to test the applicability of MISD. The research region covers an area of approximately 230 km 2 , and the target reservoir is located at an average depth of approximately 5 km. A seismic reflection profile extracted from the processed data volume is displayed in Fig. 10a, where the lateral black line denotes the interpreted horizon for the top layer of the target reservoir. It can be clearly seen that there exist strong beadlike reflections near and below the interpreted horizon. Figure 10b shows an original full-band seismic amplitude slice of the horizon corresponding to the lateral black line in Fig. 10a. Due to the obvious difference between the high-value amplitude and the surrounding lowvalue amplitude, several channel-like structures are visible, but their overall outlines are not clear enough. In addition, there is a large area of high-amplitude anomaly in the South, which is demonstrated to be associated with a karst slope fracture. As the amplitude spectrum of seismic data is shown in Fig. 10c, the band width of data is narrow, the frequency corresponding to a global maximum amplitude is approximately 27 Hz, and the frequency corresponding to a locally maximum amplitude is approximately 17 Hz. Figure 11 displays the corresponding slices of horizon through 17-Hz and 25-Hz spectral amplitude volumes extracted from 4D magnitude volumes, which are obtained by using CWT, ISD and MISD. Figure 11a, b displays channels and karst features, which cannot be separated at two close frequencies. These both CWTbased spectral amplitude slices do not provide each other with additional information. In contrast, there is relatively large difference between 17-and 25-Hz spectral amplitude slices for ISD, as shown in Fig. 11c, d. It can be observed that there are two low-frequency high-amplitude channels on the east side of the middle clearly and a large area of high-frequency high-amplitude anomaly in the South. However, only the 17-Hz spectral amplitude slice of MISD (Fig. 11e) clearly delineates channels in the middle, and especially highlights narrow channels in the East (denoted in a red ellipse), which cannot be easily interpreted in the original amplitude slice and Fig. 11a-d. Moreover, the 25-Hz spectral amplitude slice (Fig. 11f) is obviously different from 17-Hz one (Fig. 11e), and the high-frequency spectrum mainly images the karst (denoted in a black rectangle) and relatively low-amplitude channels in the Southeast. In conclusion, compared with ISD and CWT, MISD has a highest time-frequency resolution and leads to a high spatial resolution to better identify different geological bodies.

Discussion
The inverse SD imposes an additional constraint related to the time-frequency focusing or time-frequency spectrum sparsity on the data misfit derived from the inverse WT, in contrast to CWT. Therefore, the inversion-based time-frequency spectral results are improved greatly, as demonstrated in our examples. In this paper, both CWT and ISD using a typical l 1 norm can be understood as a special case of MISD. When the total iteration number is 1 and the regularization parameter is 0, the MISD-based result is equivalent to the CWT-based one. When p = 1, MISD can be derived to be the same as ISD using a typical l 1 norm. In this way, these three spectral decomposition methods can be unified in the framework of MISD. Due to only one iteration, CWT is the fastest. However, the CWT-based time-frequency resolution is the lowest, since there is no any constraint on the resulting time-frequency spectrum. For the same algorithm flow and parameters except for different p values, MISD with a typical p value of 0.5 can yield the higher time-frequency resolution than ISD, mainly because larger thresholds are adopted to make the updated time-frequency spectrum sparser within a limited number of iterations. Consequently, MISD can obtain the highest time-frequency resolution and the most flexible.
The MISD method can be extended to decompose non-Gaussian noisy seismic data, but noise in this paper should be random and satisfy the Gaussian distribution. Compared with CWT and ISD, MISD has more parameters to determine, mainly including the p value and the regularization parameter. By testing different p values, and considering some theoretical developments (e.g., Xu et al. 2010), it is feasible to take 0.5 for p in our examples. Similar to the p value, the regularization parameter can also affect the sparsity of spectral results. Generally, the larger the regularization value is, the sparser the resulting time-frequency spectrum is, but there is a larger risk of losing the weak amplitude reflection. In this paper, the selection of the regularization parameter is mainly based on the integrity and clarity of geological bodies in the spectral amplitude slice of horizon as a quality-control standard. How to pick a reasonable p and the regularization parameter adaptively (e.g., Li et al. 2019) is also worth studying in the future.
Compared with CWT, MISD is time-consuming, since a large-scale inverse problem is required to solve, but MISD brings a potential overwhelming advantage in identifying relatively close geological anomalies at the frequency along a horizon or those along the time direction. It is noticeable that four main schemes are utilized to alleviate the time-consuming issue of MISD greatly. (1) Only relatively cheap matrix-vector multiplication without any inverse operation of the largescale matrix is used to update the time-frequency spectrum.
(2) All matrix-vector multiplication operators in each iteration are implemented fast by using the Fourier transform. (3) A clever update in each iteration is employed to accelerate convergence by considering a very specific linear combination of two previous updates. (4) Parallel computing is adopted to decompose different seismic traces at the same time. Due to these schemes, MISD is not so expensive to process hundreds of square kilometers of real data, and to further obtain high-resolution time-frequency spectrum at the expense of acceptable time.
We only use the time-frequency amplitude spectrum of MISD in this paper. In fact, we can also readily obtain the corresponding time-frequency phase spectrum by implementing the inverse tangent operation on the inverted optimal time-frequency complex-value spectrum. Compared with the spectral amplitude attribute, time-frequency phase spectrum should be more sensitive to weak amplitude. This may provide some new information. The time-frequency phase spectrum attribute can be combined with time-frequency phase spectrum attribute to reduce the uncertainty of seismic interpretation in the future.

Conclusions
In this paper, we use MISD based on an l p -norm constraint of spectral coefficients to detect different geological bodies, which can reach a higher energy concentration in the time-frequency plane than CWT and ISD. There are less frequency mixing for spectral amplitude computed by using MISD than those using CWT and ISD, which is helpful for obtaining precise spectral amplitude, even when both the time variation of wavelet and the interference of thin beds complicate the time-frequency spectrum. The 4D spectral amplitude obtained by using MISD can be utilized to interpret channels, karst caves, faults and alluvial fans of  Fig. 11 Comparisons among the spectral amplitude slices obtained by using CWT at a 17-Hz and b 25-Hz, ISD at c 17-Hz and d 25-Hz, as well as MISD at e 17-Hz and f 25-Hz. There are significant differences between relatively low-and high-frequency spectral amplitude components of MISD. The MISD-based low-frequency spectral amplitude slice delineates high-amplitude channels, especially for those narrow channels in the East (the red ellipse), which cannot be indicated in the full-band amplitude slice, the CWT-based spectral slice and the ISD-based spectral slice, but only high-frequency one clearly images the high-amplitude karst slope fracture zone (the black rectangle) horizons better than the original broadband data, CWT and ISD, as demonstrated in both 3D synthetic data and 3D field data examples. Furthermore, the MISD-based isofrequency horizon magnitude slice at the chosen dominant frequency presents the less number of geological structures due to the l p -norm sparsity promotion, and thus giving a high spatial resolution to make seismic interpretation easier. MISD and the coherence algorithm will be combined to implement fine seismic interpretation in the future. The extension to quantitative seismic interpretation is also the future plan.