Seismic signal sparse time-frequency analysis by Lp-quasinorm constraint
Yingpin Chen, Zhenming Peng, Ali Gholami, Jingwen Yan, Shu Li
1 Abstract —Time-frequency analysis has been applied successfully in many fields. However, the traditional methods, like short time Fourier transform and Cohen distribution, suffer from the low resolution or the interference of the cross terms. To solve these issues, we put forward a new sparse time-frequency analysis model by using the Lp-quasinorm constraint, which is capable of fitting the sparsity prior knowledge in the frequency domain. In the proposed model, we regard the short time truncated data as the observation of sparse representation and design a dictionary matrix, which builds up the relationship between the short time measurement and the sparse spectrum. Based on the relationship and the Lp-quasinorm feasible domain, the proposed model is established. The alternating direction method of multipliers (ADMM) is adopted to solve the proposed model. Experiments are then conducted on several theoretical signals and applied to the seismic signal spectrum decomposition, indicating that the proposed method is able to obtain a higher time-frequency distribution than state-of-the-art time-frequency methods. Thus, the proposed method is of great importance to reservoir exploration.
Index Terms —ADMM, sparse representation, Lp-quasinorm constraint, seismic signal decomposition, sparse time-frequency analysis. I. I NTRODUCTION N seismic exploration, when the seismic signals go through the field containing the gas or oil, the amplitude of the high-frequency components would decrease significantly. Therefore, the time-frequency analysis is one of the key technologies to predict the position of the reservoir. At present, the time-frequency analysis technologies used in seismic exploration can be mainly divided into two categories. One is the linear time-frequency analysis technology, including short-time Fourier transform (STFT) [1], wavelet transform [2], S transform [3] and Chirplet transform [4], etc . The linear time-frequency analysis technology is calculated efficiently but suffered from low resolution because of the constraint of the Heisenberg uncertainty principle [1]. The other is the bilinear time-frequency analysis technology, which is summarized as a uniform form called Cohen distribution [5]. Cohen distribution *The corresponding author of this paper is Zhenming Peng. Yingpin Chen is a Ph.D. candidate of the School of Information and Communication Engineering, University of Electronic Science and Technology of China, Chengdu 610054 China (e-mail: [email protected]). Zhenming Peng is now the professor and doctor’s supervisor of University of Electronic Science and Technology of China, Chengdu 610054, China (e-mail: [email protected]). He is also the professor of Center for Information Geoscience, University of Electronic Science and Technology of China, Chengdu, China. He is the corresponding author of this paper. suffered from the interference of cross terms [6]. Considering the drawbacks of the linear and bilinear time-frequency technologies, many efforts are made to improve the traditional time-frequency technologies [7-14]. These efforts mainly focus on obtaining a higher resolution time-frequency distribution (TFD). Recently, sparse representation (SR) [15-17] has become a new branch of technology in signal processing. Time-frequency analysis technologies based on SR have appeared. For example, Liu et al. [16] proposed a short time sparse representation algorithm based on the smooth L0-norm algorithm (STSR-SL0). Stanković et al. [18] regarded the polynomial Fourier transform as a sparse domain and obtained a sparse time-frequency distribution by the sparse representation in a sparse domain based on L1-norm. Flandrin and Borgnat [19] introduced the compressed sensing (CS) [20] into bilinear time-frequency analysis and proposed the sparse Cohen distribution (SCD) based on the under-sampling in the ambiguity domain and the L1-norm constraint. Because of the excellent performance of SCD, Whitelonis and Ling [21] applied the SCD to radar signature analysis. However, in the framework of SCD, a square kernel is adopted to obtain the under-sampling data in the ambiguity domain. As a result, the cross term interference may not be eliminated completely [19, 22]. Jokanovic et al. [23] improved the SCD by using the adaptive optimal kernel method [11]. Gholami [24] proposed a sparse time-frequency decomposition (STFD) by using the split Bregman iterations. In the framework of STFD, the Kronecker product operator is used. Therefore, a large scale matrix is inevitable, resulting in the heavy calculation. Hu et al. [25] added a L2-norm regularization based on the work of Gholami to adjust the sparsity of the spectrum. Wang et al. [26] proposed the sparse S transform by using L1-norm. In summary, the sparse time-frequency analysis methods are mainly based on L1-norm constraint. However, the L1-norm is only the convex relaxation of L0-norm [27]. The th p pow of Lp-norm ( Np pi ip i x x , Ali Gholami is with the Institute of Geophysics, University of Tehran, Tehran 14155/6466, Iran (e-mail: [email protected]). Jingwen Yan is now the professor and doctor’s supervisor of Shantou University, Shantou 515063, China (e-mail: [email protected]). Shu Li is a Ph.D. candidate of the School of Information and Communication Engineering, University of Electronic Science and Technology of China, Chengdu 610054 China (e-mail: [email protected]). Seismic signal sparse time-frequency analysis by Lp-quasinorm constraint
Yingpin Chen, Zhenming Peng*, Member, IEEE, Ali Gholami, Jingwen Yan and Shu Li I p , for simplicity, we name pi p x as Lp-quasinorm) is another relaxation of L0-norm. In fact, L1-norm constraint is a particular case of Lp-quasinorm. Gribonval and Nielsen [28] proved the uniqueness of the solution to Lp- quasinorm based sparse representation (LpSR) theoretically. Recently, Woodworth and Chartrand [29] pointed out that the Lp-quasinorm is better able to approximate the original L0-norm than L1-norm and developed an iterative Lp-quasinorm shrinkage (LpS) solving LpSR. Woodworth and Chartrand also proved the convergence of iterative LpS for solving LpSR. Then, the generalized shrinkage operator was applied in sparse representation. For instance, Zhang et al. [30] applied the LpS operator to computed tomography image reconstruction which provided a better quality of reconstruction than the one by using L1-norm constraint. Zuo et al. [31] adopted the LpS instead of L1-norm based soft-thresholding shrinkage in image blind deconvolution and obtained a satisfactory performance. Since Lp-quasinorm is the approximation of the L0-norm with a degree of freedom, which can obtain a more precise solution than the L1-norm, we put forward a Lp-quasinorm based local time inversion model in the frequency domain. The advantages of Lp-quasinorm regularization are listed as follows. 1) The LpS operator may lead to a precise and convergent solution. 2) Lp-quasinorm is more flexible than L1-norm. This may fit the degree of sparsity with respect to the processed signal. 3) The Lp-quasinorm feasible domain makes the solution robust to noise. The proposed model absorbs the truncating thought of STFT, that is, cutting off the signal by a sliding window. In this way, the size of the involved matrix would be very small, which may decrease the amount of calculation. To reduce the interference of the false frequency introduced by the surrounding data in the proposed model, we adopt Gaussian window as the sliding window. Regarding the short time measurement as a part of signal reconstructed by sparse spectrum, the relationship between the short time measurement and the sparse spectrum is then established. Based on this relationship, the Lp-quasinorm regularization is adopted to the proposed model, fitting the prior sparsity information. To solve the proposed model, the alternating direction method of multipliers (ADMM) [32] is adopted. In order to evaluate the performance of the proposed method, experiments based on several theoretical signal data and real seismic signals are conducted by comparing with some state-of-the-art time-frequency analysis methods. The indicators we adopt are the peak signal to noise ratio (PSNR), Renyi entropy, concentration measurement (CM) , relative error (RE) and cost time [33-35]. It can be found from the experiments that the proposed method has the capacity of obtaining an accurate time-frequency representation. The main contributions are listed as follows. 1) We put forward a new time-frequency inversion model, which is a general framework for short time spectrum inversion problem. Therefore, the other sparse representation methods can be easily adopted in the proposed model. 2) The LpS is adopted to the proposed model, which can obtain a precise TFD. 3) Considering the form of the proposed model, we solve the model by ADMM, which provides a precise solution of the model. The rest of the paper is organized as follows. In Section Ⅱ , we provide some preliminary knowledge with regard to SR. In Section Ⅲ , we introduce SR into STFT and put forward a Lp-quasinorm regularization based time-frequency inversion model and discuss how to solve the model by ADMM in detail. In Section Ⅳ , based on three theoretical signal data, the experiments are conducted to compare with the traditional methods and some state-of-the-art sparse time-frequency methods. In Section Ⅴ , the proposed method is applied to the field of reservoir exploration. Finally, we summarize this paper in Section Ⅵ . II. P RELIMINARY
Sparse representation is a very active branch of signal processing. A diagrammatic sketch of sparse representation is shown in Fig. 1. M y is the observation of SR. ( ) M N
M N is the dictionary of SR. N x is the sparse representation coefficient in which the number of non-zero entries is small. It is shown in Fig. 1 that there are only two non-zero entries. Clearly, to represent y , only the second and fifth columns of the dictionary, which are boxed in red, are selected. The SR framework can be modeled as the P issue shown in (1) : min , . . , P s t x y x (1) where L0-norm x is defined as the number of the non-zero entries of x . Since the issue described in (1) is a non-deterministic polynomial hard (NP-Hard) problem, people often relax the P issue as the P issue shown in (2). : min , . . , P s t x y x (2) where = N ii x x is the L1-norm of x . III. P ROPOSED METHOD A. The sparse time-frequency inversion model
Signals in the real world are non-stationary. However, the Fig. 1. The diagrammatic sketch of the sparse representation. x y N s as the sub-signals ( 1,2, , ) Mi i N s where M is an odd number denoting the length of the window. In order to ensure that the length of every sub signal is M , we need to pad the original signal on the right and left boundaries with ( 1) / 2 M zeros. After padding, the ( 1) / 2 M i -th entry of the padding signal is regarded as the center of i s . Then, we truncate the padding signal by extracting the center of i s and the ( 1) / 2 M entries on both sides of the center from the padding signal. After padding and truncating, i s is obtained. In order to reduce the influence of the data far away from the center of i s and enhance the effect of data near the center, a Gaussian window M g is needed to reweight i s , that is = , i i y g s (3) where the symbol denotes the componentwise multiplication. Suppose Ni x is the sparse spectrum with respect to i y , thus, the corresponding signal in the time domain is -1 1 Ni F x where N N F denotes the inverse Fourier transform matrix. However, the size of observation vector i y is Mi y , thus, a matrix is required to establish the relationship between i y and -1 i F x . -1 N Ni F x is the reconstructed signal in time domain, which is assumed to an approximation of the observed measurement i y . Regarding i y as the first M entries of -1 i F x , consequently, i y and i x can be modeled as -1 , i i y SF x (4) where M N S is the selecting matrix, which is defined as , O S I (5) where
M M I is the identity matrix and ( ) M N M O is the zero matrix whose entries are all zeros. The F in (4) is the Fourier transform matrix defined as NN N NNN N NN N N NN N N
W W WW W WW W W F (6) where N W j N . In a traditional STFT framework, after truncation and reweighting, the length of i y is M . However, the size of the spectrum is N . Therefore, a zero-padding operation is carried out on i y . Therefore, extra frequency components are brought into STFT. To avoid extra frequency components, a reverse idea is created. In fact, the truncation observation in STFT can be regarded as the observation in SR. Comparing the definition of with (4), we can find that the dictionary -1 SF . S plays a role of selecting the first M entries of -1 i F x . To illustrate the differences between L1-norm and Lp-quasinorm, we show the contour line of Lp-quasinorm in Fig. 2. As is shown in Fig. 2, the smaller the parameter p is, the sparser the feasible domain of Lp-quasinorm is. In Fig. 3, we demonstrate the advantages of Lp-quasinorm. Suppose the signal is disturbed by Gaussian noise (the standard deviation is ) as seen in Fig. 3(a). Obviously, the dotted line in Fig. 3(c) would intersect the contour line of Lp-quasinorm near the axes, inducing a sparser solution (which is also robust to noise) than the one of L1-norm and L2-norm. As is mentioned, Lp-quasinorm may be a better choice. Therefore, we establish the model as (7). -1 : min , . . = , pp i i ip P s t SF x y x (7) where i p x is the Lp-norm defined as =( ) N p pip i x x . In the proposed model, the local sparse frequency is the inverse objective, while we only focus on the similarity between the local measurement and the observation interval of the objective in the time domain. Thus, extra frequency components would not be brought in the reconstructed spectrum. (a) (b) (c) (d) Fig. 2. The contour line of Lp-quasinorm. (a) p=2; (b) p=1; (c) p=0.8; (d) p=0.6. B. The solver based on ADMM and LpS operator
The issue shown in (7) is a constrained problem, which is equivalent to the unconstrained form shown in (8), i pi i ip x x y x (8) where is the coefficient, which balances the Lp-quasinorm regularization and the fidelity term i i y x . To solve (8) by ADMM, the splitting variable i z x and the Lagrange coefficient are required. Then we have the augmented Lagrangian function as i p i ipi i J x , z z y xz x z x (9) where is the penalty coefficient. For the convenience of completing the square in subsequent calculations, we set = z , then (9) becomes (10).
22, 22 i p i ipi i J x zz z y xz z x z x (10) With respect to the i x subproblem, the sub-objective function is
22 ( ) ( )2 2 i i k ki i i J x x y x z x z (11) Setting the first-order derivative of i J x with regard to i x zero, that is, i i J x x , then we have ( 1) 1 ( ) ( ) ( ) ( ( )), k H H k ki i x J y z z (12) where N N J is the identity matrix. It should be pointed out that H N N J is a large scale matrix and the computation complexity of calculating ( ) H J directly is ( ) O N . To decrease the computation complexity of ( ) H J , the following theorem is required. Theorem 3.1 : with respect to any matrix N K A , N M B , M K C , the matrix ( ) A BC satisfies the following equation [36], ( ) ( ) ,
A BC A A B I CA B CA (13) where
M M I is the identity matrix. Let = A J , = H B , = C . Obviously, ( ) H I is equal to the expression shown in (14). ( ) = ( ) H H H
J J I (14) Since M is usually much smaller than N , as a result, the computation complexity of ( ) H M M I becomes ( ) O M . Therefore, the computation complexity of ( ) H J becomes ( ) O N M . In this way, we reduce the computation complexity of the whole algorithm remarkably. With respect to the z sub-problem, fixing z and i x , the sub-objective function becomes
2( ) ( 1) ( 1) 2 min{ , }.2 p k k ki ip J z z z z z x z x (15) By using the method of completing square, then we have
2( 1) ( ) ( +1) ( +1) 22( ) ( +1) ( +1) 2( ) 2 ( ) 2 2( +1) ( ) ( ) 222( +1) ( ) 2 arg min{ , }2= arg min{ , 2+ ( ) ( ) }= arg min{ ( ) }2= arg min{ }.2 pk k k ki ipp k k ki ipk kp k k kipp k kip zzzz z z z z x z xz z z x z xz zz z x z zz z x z (16) The solution of (16) is ( 1) ( +1) ( ) ( + , ), k k kp i shrink z x z (17) With regard to real number, Lp shrinkage operator is [29] [37] ( , )= max{ } . ppp shrink - ,0 (18) (a) (b) (c) Fig. 3. The feasible region of Lp-quasinorm. (a) p=2; (b) p=1; (c) 0
STFA-LpS
Input: , , p s . Output: ( 1,2, , ) i i N x . Initialize: (1) (1) (1) , , , , , , , , i M N Max = 0 = 0 = 0 x z z .
1: Truncate the sub-signal ( 1,2, , ) i i N s from the processed signal and weight the sub-signal as shown in (3); Repeat For k Max do ( 1) 1 ( ) ( ) ( ) ( ( )) k H H k ki i x J y z z ; 5: ( 1) ( +1) ( ) ( + , ) k k kp i shrink z x z ; 6: ( +1) ( ) ( 1) ( 1) ( ) k k k ki z z x z ; 7: End for i i ; Until every weighted sub-signal i y is processed; 10: ( 1) ( 1) = ( ) ( 1,2, , ) k ki i fftshift i N x x ; 11: Return ( 1) ( 1,2, , ) ki i N x as ( 1, 2, , ) i i N x . where the symbol Max denotes the maximal iteration number of the iterations. In our experiments, the parameters are set as =0.1 p , =1 , =0.5 , =1 , =11 M , Max . It is worth pointing out that all the sub-signals can be processed in a parallel fashion. IV. N UMERICAL EXPERIMENTS A. The running environment and evaluating indicators
The experiments are all implemented in MATLAB 2016 and executed on a laptop equipped with an Intel Core i5-3210M 2.5GHz CPU and 6G of RAM. In this section, we evaluate the proposed method by general measurements, which are the concentration measurement (CM), Renyi entropy, peak signal to noise ratio (PSNR), relative error (RE) and cost time. ( , ) d d ,( , ) d d
TF t f t fCM TF t f t f (23) where the symbol ( , ) TF t f represents the time-frequency distribution. The larger the CM is, the better the quality of TFD is. The Renyi entropy shown in (24) is adopted to measure the aggregation of TFD. d d1( ( , )) log ( , ) ,1 R TF t f TF t f t f (24) where =3 . The smaller the R is, the better the quality of TFD is. In this paper, we carry out experiments on three ideal model, whose ideal time-frequency distribution is known. Therefore, the PSNR can be used to evaluate the methods. All the TFDs are normalized to the range [0,1] and then compared with the ideal TFD by PSNR. The PSNR is defined as ( )( , )=10lg ,1 ( ) N N ij iji j
MAXPSNR N
YX Y X Y (25) where Y represents the ideal time-frequency image and the X is the time-frequency distribution obtained by certain time-frequency analysis algorithm. The larger the PSNR is, the greater the similarity of X and Y is. The relative error is defined as ( , )= . RE X - YX Y X (26) As is defined in (26), the smaller the relative error is, the better the quality of time-frequency distribution is. B. Signal 1 : mono-component signal Linear frequency modulated (LFM) signal is often used to measure the quality of time-frequency analysis. The closed form of the tested LFM signal is defined as ( )=exp( ), s t jπt (27) where time range is from -8 s to 8 s and the sampling rate is 16 Hz. The instantaneous phase of the LFM signal is ( ) t t , thus, the instantaneous frequency is d tf t tdt (28) Therefore, the ideal TFD should be a straight line. Table I provides the evaluating indicators of TFDs shown in Fig. 4, where the TFDs are obtained by STFT [1], synchrosqueezed wavelet transform (SSWT) [38], Hilbert-Huang Transform (HHT) [13], STFD [24], SCD [20], STSR-SL0 [16] and the proposed method with regard to signal 1. It is obvious that the TFD obtained by the proposed method has the best performance in PSNR, Renyi entropy, and CM. As is shown in Fig. 4 that the result of the proposed method is closest to the ideal TFD shown in Fig. 4(b). The other TFDs have some disadvantages shown in Fig. 4. For example, the TFD by STFT does not have a satisfactory time-frequency resolution. Since the SSWT is based on the wavelet transform, the SSWT can only represent the positive frequency components shown in Fig. 4(d). In addition, the SSWT losses the low-frequency information comparing with other methods. Fig. 4(e) shows that HHT can only represent the positive frequency component because of the use of Hilbert transform. To demonstrate the effectiveness of the LpS, we figure out the relative error curve in Fig. 5 with (a) (b) (c) (d) (e) (f) (g) (h) (i) Fig. 4. The real part of the processed signal and the TFDs. (a) The real part of the processed signal; (b) the ideal time-frequency distribution; (c) the result of STFT; (d) the result of SSWT; (e) the result of HHT; (f) the result of STFD; (g) the result of SCD; (h) the result of STSR-SL0; (i) the result of the proposed method. TABLE
I T HE E VALUATING I NDICATORS OF
TFD
S IN F IG . SSWT 23.04 9.28 4.0E-3 0.17 HHT 22.90 9.13 3.3E-3 0.72 STFD 23.65 10.10 1.0E-3 8.73 SCD 22.12 9.94 1.2E-3 2.67 STSR-SL0 20.55 11.09 9.0E-4 3.36 STFA-LpS t in Fig. 4 (b). As is shown, the smaller the p is, the smaller the RE with respect to the ideal TFD is. C. Signal 2 : parabola frequency modulated signal The second signal we adopt is the parabola frequency modulated signal, whose closed form is shown as follow ( ) exp( ),12 ts t j (29) where the time range is [ 8 ,8 ] t s s and the sampling rate is 16Hz. The instantaneous phase of the second signal is ( ) 12 tt , consequently, the instantaneous frequency is d t tf t dt (30) Thus, the shape of the ideal time-frequency distribution with respect to the second signal is a parabola. Therefore, the instantaneous frequency of the signal changes fast. As a result, the short time truncated signal may include many components. In other words, the spectrum of the truncated signal is not sparse. However, the ideal time-frequency distribution of parabola frequency modulated signal is sparse according to (30). Under this circumstance, the weighted function is required to eliminate the influence of the data near the center of the truncated signal and force the sparse spectrum to be sparse. With regard to signal 2, we provide the TFDs by STFT, SSWT, HHT, STFD, SCD, STSR-SL0 and the proposed method. The performances of the algorithms are recorded in Table Ⅱ . It can be found that the PSNR and CM obtained by our method are the highest while the Renyi entropy by our method is the smallest. Therefore, the proposed algorithm outperforms the other algorithms. Fig. 6 provides the real part of signal 2 and TFDs computed (a) (b) (c) (d) (e) (f) (g) (h) (i) Fig. 6. The real part of the processed signal and the TFDs. (a) The real part of the processed signal; (b) the ideal time-frequency distribution; (c) the result of STFT; (d) the result of SSWT; (e) the result of HHT; (f) the result of STFD; (g) the result of SCD; (h) the result of STSR-SL0; (i) the result of the proposed method. TABLE Ⅱ T HE E VALUATING I NDICATORS OF
TFD
S IN F IG . SSWT 23.71 9.92 3.0E-3 0.30 HHT 22.90 9.13 3.6E-3 0.19 STFD 22.27 10.01 1.2E-3 8.72 SCD 19.05 10.02 1.8E-3 2.58 STSR-SL0 20.83 9.86 9.3E-4 3.22 STFA-LpS by different methods. In Fig. 6, it is observed that the TFD calculated by the proposed method shown in figure 6(i) is the closest TFD with respect to the ideal TFD shown in Fig. 6(b). On the one hand, by using the Gaussian window, the spectrum becomes very sparse and avoids the influence of the surrounding data like STFT shown in Fig. 6(c). On the other hand, it is not interfered by the cross terms like Fig. 6(g) and the false frequency components like HHT shown in Fig. 6(e). It is observed that the result of SSWT shown in Fig. 6 (d) losses some frequency information among the range [ 4 , 4 ] t s s while the result of ours depict the time-frequency trajectory precisely. D. Signal 3 : multi-component signal The aforementioned signals are mono-component signals. However, the real signals are always composed of many components. Without loss of generality, we carry out the comparisons on a multi-component signal, which is ( ) exp( ) exp( ) exp( 8 )12 12exp( (14 cos(2 ))), t ts t j j j tj t t (31) where [ 8 ,8 ] t s s and the sampling rate is 16 Hz. As is defined in (31), the third signal contains the stationary component exp( 8 ) j t , the fast time-varying frequency component exp( (14 cos(2 ))) j t t and the slow time-varying frequency components exp( )12 tj and exp( )12 tj . With regard to signal 3, we provide the TFDs by STFT, SSWT, HHT, STFD, SCD, STSR-SL0 and the proposed method. The performances of the algorithms are recorded in Table Ⅲ and TFDs are listed in Fig. 7. It is observed that the PSNR, Renyi entropy, and CM of the (a) (b) (c) (d) (e) (f) (g) (h) (i) Fig. 7. The real part of the processed signal and the TFDs. (a) The real part of the processed signal; (b) the ideal time-frequency distribution; (c) the result of STFT; (d) the result of SSWT; (e) the result of HHT; (f) the result of STFD; (g) the result of SCD; (h) the result of STSR-SL0; (i) the result of the proposed method. TABLE Ⅲ T HE E VALUATING I NDICATORS OF
TFD
S IN F IG . SSWT 20.53 11.48 8.17e-4 0.17 HHT 20.45 10.35 1.8E-3 0.67 STFD 21.03 9.99 1.9E-3 10.61 SCD 20.54 11.43 8.60E-4 4.42 STSR-SL0 19.96 12.30 3.7E-4 2.75 STFA-LpS
9 proposed method with regard to the third signal are still best shown in Table Ⅲ . As is shown in Fig. 7, the result of STFT shown in Fig. 7(c) does not have a high time-frequency resolution. It is observed in Fig. 7(d)-(e) that the result of SSWT and HHT cannot show the negative frequency component. The result of STFD shown in Fig. 7(f) cannot reflect the fast time-varying frequency component. The TFD calculated by SCD is interfered by cross terms shown in Fig. 7(g). In terms of STSR-SL0 shown in Fig. 7(h), because the sliding window of STSR-SL0 is a rectangular window , some false frequency is involved in the result of STFR-SL0. Thus, the TFD by using the proposed method may be the most satisfactory among the TFDs shown in Fig. 7(i). V. A PPLICATION IN SEISMIC SIGNAL PROCESSING A. Experiment on single seismic signal
The seismic signal shown in Fig. 8(a) comes from an oil field of Sichuan Basin, China. To illustrate the performance of our method, we carry out experiments on STFT, SSWT, SCD, STFD, and the proposed method shown in Fig. 8 (b)-(f). Considering that there is no standard TFD with respect to the seismic signal, we only use Renyi entropy, CM and cost time to measure the performances of TFDs. The performances are recorded in Table Ⅳ . It is worth pointing out that the Renyi entropy by the proposed method is the smallest and CM obtained by our method is the largest, indicating that the proposed method can complete with state-of-the-art methods in the quality of time-frequency distribution. In terms of cost time, the sparse time-frequency methods are slower than the traditional methods while the SSWT runs fastest among the sparse time-frequency methods. As is seen in Fig. 8, the seismic signal is a signal with fast frequency changing. As it is observed in Fig. 8(a), the signal in the range 1850-1950 ms starts with low frequency, then the frequency goes even lower and gets to the lowest frequency point at the time 1910 ms for signal near this area changes slowest. Then in the range 1910-1950 ms, the frequency of signal becomes larger gradually. Thus, the ideal time-frequency trajectory of the seismic signal in the range 1850-1950 ms should start from a low value, then decrease to the lowest position at the time 1910 ms and finally increase to a higher value as the process progressed. To a certain extent, the result of STFT shown in Fig. 8(b) reflects the process. However, the resolution of STFT is unsatisfactory. For the operation of synchrosqueezing, the energy in the TFD of SSWT (a) (b) (c) (d) (e) (f) Fig. 8. The seismic signal and the TFDs. (a) The real seismic signal; (b) the result of STFT; (c) the result of SSWT ; (d) the result of SCD; (e) the result of STFD; (f) the result of STFA-LpS. TABLE Ⅳ T HE E VALUATING I NDICATORS OF
TFD
S IN F IG . SSWT 11.95 1.0E-3 0.58 SCD 9.97 2.2E-3 2.50 STFD 10.50 1.0E-3 3.83 STFA-LpS X B. Application in seismic signal spectral decomposition
In this part, we show the performances of different time-frequency analysis methods based on the two-dimensional seismic signal derived from the Sichuan basin, China. The seismic data is shown in Fig. 9 where the sample rate is 512 Hz. The seismic section from 1814 to 2068 ms is the target stratum for seismic interpretation. In the figure, the time range is 1814–2068 ms, as the ordinate, and the range of common depth point (CDP) is from 51 to 150, as the abscissa. The well X near the CDP 131 is a prolific natural gas well, whose daily capacity is 54.18*104 m3/d. According to the exploration information, the location of the well full of gas is 1932-1948 ms. The area full of gas is marked by the red dotted line shown in Fig. 9. Fig.10 shows the frequency slices calculated by the proposed method. Each slice takes about 30 s. Clearly, the main frequency slice at 30 Hz indicates the location full of gas precisely. It is worth pointing out that the 15 Hz frequency slice shows the low-frequency shadow phenomenon. That is, the amplitude of low-frequency slice below the reservoir is larger than the amplitude on the reservoir. As is shown in Fig. 10 (a), there is a stronger response marked by the white dotted line than the one of the reservoir. The area marked by the white dotted line is below the reservoir. On the contrary, the amplitude of this area becomes smaller than the amplitude of the area full of gas in the 30 Hz slice. The result fits in with the low-frequency shadow phenomenon [39, 40]. In the high-frequency slice shown in Fig. 10(c)-(d), the area containing natural gas has an obvious energy attenuation, which is consistent with high-frequency attenuation in the gas-bearing reservoir. For comparison, we extract the frequency slices by STSR-SL0. The slices are shown in Fig. 11. Each slice obtained by STSR-SL0 takes about 70 s. It is observed that the slices of each frequency obtained by the proposed method have higher resolution than the ones of STSR-SL0. In summary, the frequency slices obtained by the proposed method fit the real frequency characteristic of the reservoir without the interference of cross terms. Therefore, the proposed method is capable of being applied to the reservoir exploration because of the high time-frequency resolution. VI. C ONCLUSION
In this paper, we first reviewed the SR theory, and then, explored the relationship between SR and STFT. Since the short time sampling of STFT can be regarded as the measurement of SR, we build up the local time sparse spectrum construction model. In order to obtain a sparser spectrum, the Lp-quasinorm is used to substitute for the L1 regularization. As the proposed model is established in the framework of STFT, the time-frequency representation by our model is not interfered by the cross terms. To solve the proposed model, the ADMM framework and the LpS is adopted. As a result, the constrained problem is changed to an unconstrained problem while the variables in the augmented Lagrangian function are all decoupled. Consequently, the proposed model can be solved by several simple sub-problems. For the constraint of the sparsity in the
Lp-quasinorm feasible domain, we can easily obtain the time-frequency representation precisely. According to the special. (a) (b) (c) (d)
Fig. 10. The single frequency slices obtained by the proposed method. (a) 15 Hz frequency slice; (b) 30 Hz frequency slice; (c) 45 Hz frequency slice; (d) 55 Hz frequency slice.
11 form of the dictionary in the proposed method, we optimized the matrix inversion which reduces the running time strikingly Based on three theoretical signals, we carried out the experiments comparing the proposed method with the other methods and evaluated the methods by PSNR, Renyi entropy, RE, CM, and the cost time. The results showed that the proposed method is capable of calculating high-resolution time-frequency distribution. Finally, the proposed method is applied in the reservoir exploration and obtained high-resolution spectrum decomposition, which is fitting with the low-frequency shadow phenomenon and the high-frequency attenuation. Although we try our best to optimize the matrix inversion, which can reduce the computational complexity, the proposed method is not as fast as the traditional time-frequency analysis methods. Therefore, we will focus on improving the efficiency of calculating the proposed model in the future. A
CKNOWLEDGMENT
This work is supported by National Natural Science Foundation of China [grant numbers 61571096, 41274127, 61775030], Natural Science Foundation of Fujian Province [grant number 2015J01270], Education And Scientific Research Foundation of Education Department of Fujian Province for Middle-aged and Young Teachers [grant number JAT170352], and the Open Foundation of Digital Signal and Image Processing Key Laboratory of Guangdong Province [grant number 2017GDDSIPL-01]. We appreciate Dr. Xiaobo Qu of Xiamen University for constructive suggestions that greatly improved the manuscript. References [1] D. Gabor, “Theory of communication. Part 1: The analysis of information,”
Journal of the Institution of Electrical Engineers - Part III: Radio and Communication Engineering, vol. 93, no. 26, pp. 429-441, 1946. [2] M. J. Shensa, “The discrete wavelet transform: wedding the a trous and Mallat algorithms,”
IEEE Trans. Signal. Process., vol. 40, no. 10, pp. 2464-2482, 1992. [3] R. G. Stockwell, L. Mansinha, and R. Lowe, “Localization of the complex spectrum: the S transform,”
IEEE Trans. Signal. Process., vol. 44, no. 4, pp. 998-1001, 1996. [4] S. Mann, and S. Haykin, “Adaptive chirplet transform: an adaptive generalization of the wavelet transform,”
OptEn, vol. 31, no. 6, pp. 1243-1256, 1992. [5] L. Cohen, “Generalized phase-space distribution functions, ” J. Math. Phys, vol. 7, no. 5, pp. 781-786, 1966. [6] L. Stankovic, M. Dakovic, and T. Thayaparan,
Time-frequency signal analysis with applications , Boston: Artech house, 2014. [7] Y. Chen, Z. Peng, Z. He, L. Tian, and D. Zhang, “The optimal fractional Gabor transform based on the adaptive window function and its application,”
Appl. Geophys., vol. 10, no. 3, pp. 305-313, 2013. [8] Y. Wang, and Z. Peng, “The optimal fractional S transform of seismic signal based on the normalized second-order central moment,”
J. Appl. Geophys., vol. 129, pp. 8-16, 2016. [9] Y. Yang, Z. K. Peng, W. M. Zhang, and G. Meng, “Frequency-varying group delay estimation using frequency domain polynomial chirplet transform,”
Mech. Syst. Signal. Pr., vol. 46, no. 1, pp. 146-162, 2014. [10] Y. Yang, Z. Peng, X. Dong, W. Zhang, and G. Meng, “Nonlinear time-varying vibration system identification using parametric time–frequency transform with spline kernel,”
Nonlinear. Dynam., vol. 85, no. 3, pp. 1679-1694, 2016. [11] R. G. Baraniuk, and D. L. Jones, “A signal-dependent time-frequency representation: optimal kernel design,”
IEEE Trans. Signal. Process., vol. 41, no. 4, pp. 1589-1602, 1993. [12] Y. Chen, Z. Peng, Z. Cheng, and L. Tian, “Seismic signal time-frequency analysis based on multi-directional window using greedy strategy,”
J. Appl. Geophys., vol. 143, pp. 116-128, 2017. [13] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,”
Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 454, no. 1971, pp. 903-995, 1998. (a) (b) (c) (d)
Fig. 11. The single frequency slices obtained by STSR-SL0. (a) 15 Hz frequency slice; (b) 30 Hz frequency slice; (c) 45 Hz frequency slice; (d) 55 Hz frequency slice. [14] N. E. Huang, and Z. Wu, “ A review on Hilbert ‐ Huang transform: method and its applications to geophysical studies,”
Rev. Geophys., vol. 46, no. 2, pp. 1-23, 2008. [15] G. Mohimani, M. Babaie-Zadeh, and C. Jutten, “Fast sparse representation based on smoothed ℓ 0 norm,” in Independent Component Analysis and Signal Separation, London, UK, 2007, pp. 389-396. [16] Z. Liu, P. You, X. Wei, D. Liao, and X. Li, “High resolution time-frequency distribution based on short-time sparse representation,”
Circ. Syst. Signal. Pr., vol. 33, no. 12, pp. 3949-3965, 2014. [17] J. Wright, A. Y. Yang, A. Ganesh, S. S. Sastry, and Y. Ma, “Robust face recognition via sparse representation,”
IEEE T. Pattern. Anal., vol. 31, no. 2, pp. 210-227, 2009. [18] S. Stanković, I. Orović, and L. Stanković, “Polynomial Fourier domain as a domain of signal sparsity,”
Signal. Process., vol. 130, pp. 243-253, 2016. [19] P. Flandrin, and P. Borgnat, “Time-frequency energy distributions meet compressed sensing,”
IEEE Trans. Signal. Process., vol. 58, no. 6, pp. 2974-2982, 2010. [20] D. L. Donoho, “ Compressed sensing,”
IEEE T. Inform. Theory., vol. 52, no. 4, pp. 1289-1306, 2006. [21] N. Whitelonis, and H. Ling, “Radar signature analysis using a joint time-frequency distribution based on compressed sensing,”
IEEE Trans. Antennas Propag., vol. 62, no. 2, pp. 755-763, 2014. [22] P. Flandrin, N. Pustelnik, and P. Borgnat, “On Wigner-based sparse time-frequency distributions,” in 2015 IEEE 6th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, Saint Martin, French, 2015, pp. 65-68. [23] B. Jokanovic, M. G. Amin, Y. D. Zhang, and F. Ahmad, “Time-frequency kernel design for sparse joint-variable signal representations,” in 2014 Proceedings of the 22nd European Signal Processing Conference (EUSIPCO), Lisbon, Portugal, 2014, pp. 2100-2104. [24] A. Gholami, “Sparse time–frequency decomposition and some applications,”
IEEE Trans. Geosci. Remote., vol. 51, no. 6, pp. 3598-3604, 2013. [25] J. Hu, X. He, W. Li, H. Ai, H. Li, and J. Xi, “Parameter estimation of maneuvering targets in OTHR based on sparse time-frequency representation,”
J. Syst. Eng. Electron., vol. 27, no. 3, pp. 574-580, 2016. [26] Y. Wang, Z. Peng, and Y. He, “Time-frequency representation for seismic data using sparse S transform,” in 2016 2nd IEEE International Conference on Computer and Communications, Chengdu, China, 2016, pp. 1923-1926. [27] D. L. Donoho, and M. Elad, “Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ1 minimization,”
Proceedings of the National Academy of Sciences, vol. 100, no. 5, pp. 2197-2202, 2003. [28] R. Gribonval, and M. Nielsen, “Sparse representations in unions of bases,”
IEEE T. Inform. Theory., vol. 49, no. 12, pp. 3320-3325, 2003. [29] J. Woodworth, and R. Chartrand, “Compressed sensing recovery via nonconvex shrinkage penalties,”
InvPr, vol. 32, no. 7, pp. 75004-75028, 2016. [30] H. Zhang, L. Wang, B. Yan, L. Li, A. Cai, and G. Hu, “Constrained total generalized p-variation minimization for few-view x-ray computed tomography image reconstruction,”
PLoS. One., vol. 11, no. 2, pp. e0149899, 2016. [31] W. Zuo, D. Ren, D. Zhang, S. Gu, and L. Zhang, “Learning iteration-wise generalized shrinkage–thresholding operators for blind deconvolution,”
IEEE Trans. Image. Process., vol. 25, no. 4, pp. 1751-1764, 2016. [32] S. Boyd, “Distributed optimization and statistical learning via the alternating direction method of multipliers,”
Foundations and Trends® in Machine Learning, vol. 3, no. 1, pp. 1-122, 2010. [33] B. Boashash, N. A. Khan, and T. Ben-Jabeur, “Time–frequency features for pattern recognition using high-resolution TFDs: A tutorial review,”
Digit. Signal Process., vol. 40, pp. 1-30, 2015. [34] E. Sejdić, I. Djurović, and J. Jiang, “Time–frequency feature representation using energy concentration: An overview of recent advances,”
Digit. Signal Process., vol. 19, no. 1, pp. 153-183, 2009. [35] L. Stanković, “A measure of some time–frequency distributions concentration,”
Signal. Process., vol. 81, no. 3, pp. 621-631, 2001. [36] K. B. Petersen, and M. S. Pedersen, “The matrix cookbook,”
Technical University of Denmark, vol. 7, no. 15, pp. 1-25, 2008. [37] R. Chartrand, “Shrinkage mappings and their induced penalty functions,” in 2014 IEEE International Conference on Acoustics, Speech and Signal Processing, Florence, Italy, 2014, pp. 1026-1029. [38] I. Daubechies, J. Lu, and H.-T. Wu, “Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool,”
Applied and Computational Harmonic Analysis, vol. 30, no. 2, pp. 243-261, 2011. [39] Z. He, X. Xiong, and L. Bian, “Numerical simulation of seismic low-frequency shadows and its application,”
Appl. Geophys., vol. 5, no. 4, pp. 301-306, 2008. [40] X. H. Chen, Z. H. He, S.-X. Zhu, W. Liu, and W.-L. Zhong, “Seismic low-frequency-based calculation of reservoir fluid mobility and its applications,”
Appl. Geophys., vol. 9, no. 3, pp. 326-332, 2012.
Yingpin Chen received his bachelor's degree from Fuzhou University, Fuzhou, Fujian, China, in 2009 while his major is electronic science and technology. He received his master’s degree in signal and information processing from University of Electronic and Science Technology of China, Chengdu, Sichuan, China, in 2013. He is now a staff of School of Physics and Information Engineering, Minnan Normal University, Zhangzhou, China. He is now a Ph.D. candidate at University of Electronic and Science Technology of China, Chengdu, Sichuan, China.
His research interests include time-frequency analysis, compressed sensing and convex optimization.
Zhenming Peng (M’06) received the Ph.D. degree in geo-detection and information technology from Chengdu University of Technology, Chengdu, China, in 2001. From 2001 to 2003, he was a Postdoctoral Researcher at the Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu, China. He is currently a professor at University of Electronic Science and Technology, Chengdu, China. He is also the professor of Center for Information Geoscience, University of Electronic Science and Technology of China, Chengdu, China. His research interests include image processing, radar signal processing, and target recognition and tracking.
Ali Gholami received the Ph.D. degree in geophysics (exploration seismology) from the University of Tehran, Tehran, Iran, in 2010. He is currently an Assistant Professor of geophysics with the Institute of Geophysics, University of Tehran. His research interests include inverse problems, seismic signal processing, and seismic imaging.
Jingwen Yan received his MS’s ’ degree from Changchun Geography Institute of Chinese Academy of Sciences and Ph.D. degree from Changchun Geography Institute and Changchun Institute of Optics, Chinese Academy of Sciences, Changchun, Jilin, China. Now, he is a full professor and doctor’s supervisor at Shantou University, Shantou, Guangdong, China. His research interests include wavelet analysis and applications, sparse representation of signals, remote sensing, embedded systems, and FPGA/DSP hardware design. 13 Shu Li received his B.S. degree in Communication Engineering from Hebei University, Baoding, Hebei, China, in 2007.