Seismic signal sparse time-frequency analysis by Lp-quasinorm constraint
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 is…
Authors: Yingpin Chen, Zhenming Peng, Ali Gholami
1 Abstract — Time-frequency analysis has been applied successfully in m any fields. H ow ever, th e traditional m ethods, like short time Fourier tran sform a nd Cohen distribution, suffer from the low resolution or the interference of the cross ter ms. To solve the se issues, we put forward a n ew sparse time-frequency analysis model by using t he Lp -quasinorm constraint, w hich is capable of fitting the sparsity prior knowledge in the frequency do ma i n. In the propose d model, w e regard the short ti me tr uncated data as the observation of sparse repre sentation a nd design a dictionary matrix, w hich builds up the relationship betw een the s hort time measurement and the sparse s pectrum . B ased o n t he relationship and the Lp- quasinorm feasible domain, t he proposed m odel is established. The alternating direction method of multipliers (ADMM) is adop ted to solve th e proposed model. E x periments are then conducted on several theo retical signals and applied to the seismic signal spectrum deco mposition, indicating that the proposed method is able to obtain a higher time-frequency distribution than stat e- of -the-art ti me-frequency methods. Thus, the proposed method is of great importance to reservoir exploration. Index Ter m s — ADMM, spar se representation , Lp -quasino r m constraint, seismic signal decomposition, sparse ti me-frequency analysis. I. I NTRODUCTI ON N seismic explor ation, when t he seis mic signals go through the field containing the gas or oil, the amplitude of the hig h- frequency components would decrease s ignificantly . Therefore, the time-frequency anal ysis i s one of the key tec hnologies to predict the position of the reservoir. At present, t he tim e - 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 transfor m [ 2], S transfor m [ 3] and Chirplet transform [4], etc . T he linear time-frequenc y analysis tec hnology is calc ulated ef ficiently but s uffered fro m low resolution because of the constraint of the Heisenbe rg uncertainty principle [1]. The o ther is the bilinear time - frequency anal ysis technology, which is summarized a s a uniform for m called Cohe n distribution [5]. Cohe n distribution *The corre sponding author of this pape r is Zhenming Peng. Yingpin Chen is a Ph.D. candidate of the School of Information and Communication E ngineering, University of El ectronic S cience and Technolo gy of China, Cheng du 610054 China (e -mail: 1105006 17@163.com). Zhenm in g Peng is now t he p rofesso r and doctor’s supervisor of University of Electronic Science and Technol ogy of China, Chengdu 610054, China (e- mail: zmpeng@ue stc.edu.cn). He is also the profe ssor of Center for Info rmation Geoscience, University o f Electronic Science and Technology of China , Chengdu, China. He is the correspond ing author of this pape r. suffered from the interfere nce of cr oss ter ms [6] . Consideri ng the dra wbacks of the linear and b ilinear ti me-frequen cy technologies, man y e fforts are made to improve t he trad itional time-frequency technologies [7- 14 ]. These efforts m ainl y focus on obtaining a h igher resolution tim e -freque ncy distributio n (TFD). Recently, sparse r epresentation (SR) [ 15 - 17 ] has beco me a new branch of technolog y in signal p rocessing. T ime-frequenc y analysis technologies based on SR have appea r ed . For example, Liu et al. [ 16 ] proposed a short time sparse represe ntation algorithm based on the smooth L0-nor m algorithm (STSR -SL0) . Stanković et al. [ 18 ] regarded the polyno mial Fourier transform as a sparse domain and obtained a sparse time -frequency distribution b y the spar se rep resentation in a sparse d omain based on L1 -nor m . Flandrin and Borgnat [ 19 ] introduced the compressed sensing (CS) [ 20 ] into bilinear time-freque ncy analysis and p roposed the sparse Cohen d istribution (SCD) based on the under-sampling in the ambig uity do main and the L1 -nor m constraint. Because of the excellent p erformance of SCD , W hitelonis a nd Ling [ 21 ] applied the SCD to ra dar signature analysis. However, in the framework of S CD, a square kernel i s adopted to obtain t he under-sa mpling data in the ambiguity do main. As a r esult, the cro ss term interference may not be eli minated co mpletely [ 19 , 22 ]. Jokanovic et al. [ 23 ] improved the SCD b y using the adaptive opti mal kernel method [ 11 ]. Gholami [ 24 ] proposed a sparse ti me-frequenc y decomposition (ST FD) by usi ng the split Bregman iteratio ns. In the fra mework of ST FD, t he K ronecker pro duct oper ator is used . Therefore, a large scale matrix is i nevitable , res ulting in t he heavy calculatio n. Hu et al. [ 25 ] add ed a L2-norm regularization based on the work of Gholami to adjust the sparsity of the sp ectrum. Wan g et al. [ 26 ] proposed the sp arse S transfor m b y using L1 -norm. In s ummary, th e sparse time- frequency anal ysis methods are mainly based on L1 -norm constraint. However, th e L1 -no rm is only the convex relaxation of L0 -norm [ 27 ]. T he th p pow of Lp -norm ( 1 N pp ii p i xx , Ali Gholami is w ith t he I nstitute of Geophysics, University of Tehran, Tehran 14155/6 466, Iran (e-mail: agho lami@ut.ac.ir). Jingwen Y an is now the professor and doctor’s supervisor of Shanto u University, Sha ntou 515063, China (e-mail: jwyan@stu.edu.cn ). Shu Li is a Ph.D. candidate of th e School of Infor mation and Communication Engineering, University of Ele ctronic S cience and Technolo gy of China, Cheng du 610054 China (e-mail: l ishuvip@126.com ). Seismic signal sparse time-frequency analysis by Lp -quasinorm constraint Yingpin Chen, Z henming Peng*, Mem ber, IEEE, Ali Gholami, J ingwen Yan and Shu Li I 2 01 p , f or simplicity, w e name p i p x as Lp-quasinorm) is another relaxa tion o f L0- norm . I n fact, L1 -nor m constrai nt is a particular case of Lp -q uasinorm. Gribonval and Nielsen [ 28 ] proved the uniqueness of the solution to Lp- quasinor m b ased sparse representation (LpSR) theoretica lly . Recent ly , Woodworth a nd Chartra nd [ 29 ] pointed out that the Lp - quasinorm is better able to approximate the o riginal L0 - norm than L1-nor m and d eveloped an iterati ve Lp -quasinor m shrinkage (LpS) solving L pSR . Woodworth and Chartrand also proved th e convergence o f it erative L pS f or solving LpSR . Then , the generalized s hrinkage op erator was applied in spar se representation. For instance, Zhang et al. [ 30 ] applied the LpS operator to co mputed to mography image rec onstruction which provided a better quality o f reconstruction than th e o ne b y using L1 -nor m con straint. Zuo et al. [ 31 ] adopted the LpS instead of L1 -nor m b ased soft-threshol ding shrinkage in image bli nd deconvolution and o btained a satisfactory per formance. Si nce Lp -q uasinorm is the ap proximation o f the L0 - norm with a degree of freedom, which can obtain a more prec ise solution than the L1 -norm, we put forward a Lp -q uasinorm based lo cal time inversion model in the freq uency domain. The ad vantages of Lp -quasinorm regularization ar e listed as follows. 1) The LpS oper ator may lead to a precise and convergent solution . 2 ) Lp -q uasinorm is more flexible tha n L1 -norm. Thi s may fit the degree of spar sity with re spect to the p rocessed signal. 3) The Lp -q uasinorm feasible dom ain makes the solu tion ro bust to noise. The proposed model absorbs the truncating thought o f STFT , that is, cutting off t he signal by a sliding windo w . In this way , the size of the involved matrix would be very small, whic h may decrea se the amount of calculation. T o r educe the interference of th e false frequency introduc ed by the surrounding data in the proposed m odel , w e adopt Gaussia n w indo w as the sliding window . Re garding the short time measurement as a par t of signal reco nstructed by sparse spectru m, t he relatio nship between the short time measurement an d the sparse spectrum is then establi shed. Based on this relatio nship, the Lp -quasinor m regularization is adopted to the proposed model, fitting the pr ior sparsity in formation. T o solve t he prop osed model , t he alternating direction method of multipliers (ADM M) [ 32 ] is adopted . In order to evaluate the performance of the p roposed method, experiments based on several theoretical signal data and real seismic signals are conducted b y comparing with some state- of - the-art time-frequenc y analysis methods. The indicators we adopt are the pe ak si gnal to noise r atio ( PSNR), Renyi entrop y , concentration meas urement (CM) ， relative e rror (RE) and cost time [ 33 - 35 ] . It ca n be found fro m the experiments that the prop osed method has the capacit y o f obtaining a n accurate ti me-frequency representation. T he main contributions are listed as follows. 1 ) W e p ut forward a new time-frequency inversio n model, which is a general fra mework for sh ort time spectrum inv er sion problem . T herefore, the oth er sparse representation methods can b e easily adopted in t he proposed model. 2) The LpS i s adopted to the propo sed model, which can obtain a precise T FD . 3) C onsideri ng th e form of th e proposed model, we solve the model by ADMM, which provides a prec ise solution of the model. The rest of the pap er is or ganized as follo ws. In Sectio n Ⅱ , we provide so me preliminar y knowledge with re gard to SR . In Section Ⅲ , we introduce SR into ST FT and put for ward a Lp - quasinorm regularization based time-frequency invers ion model and discuss how to solve the m odel by ADMM in detail. In Section Ⅳ , b ased on three theoretical signal data, the experiments are conducted to co mpare with the traditional methods and some state - of -the-art sparse time -frequenc y methods. In Sectio n Ⅴ , the pro posed method is applied to the field of reservoir exploratio n. Finally, we summarize this paper in Section Ⅵ . II. P RELI MINARY Sparse representation is a very active branch o f signal processing. A diagra mmatic sketch of sparse rep resentation is shown in Fig. 1. 1 M y is the observation of SR. ( ) MN M N is the dictionar y of SR. 1 N x is the sparse r epresentation coefficient in which t he number of no n- zero en tries is small. It is shown in Fig. 1 that there are only two non-zero entries. Clearly , to represent y , o nly t he seco nd and fifth col umns of the d ictionar y , which are boxed in red, are selected. The SR framework can be modeled as the 0 P issue shown in (1) 0 0 : m in , . . , P s t x y x (1) where L0 -norm 0 x is defined as the number of t he non -zero entries of x . Since the issue describ ed in (1) is a non-deter ministic polynomial hard (NP -Hard) p roblem, p eople o ften relax the 0 P issue as the 1 P issue shown in ( 2 ). 1 1 : m in , . . , P s t x y x (2) where 1 1 = N i i xx is the L1-nor m of x . III. P ROPOSED METHOD A. The sparse time-frequ ency inversion model Signals in the real world are non-stationar y. Ho wever, the Fig. 1. The diagramm atic sketch of the s p arse repre sentati on . x y 3 truncated signal s are assu med to b e stationary. Therefore, the truncated signal weight ed by a Ga ussian windo w ca n be regarded as approximatively spar se in the frequenc y d omain. Our motivation is to establish a r elationship b etween SR and STFT and obtain a h igh-resolution time-frequenc y distribution. In add ition , as we establish the framework base d on the short time measurement, there is n o cross-term interfere nce in t he proposed time-frequency inversion model. To introduce SR i nto ST FT , we truncate the original si gnal 1 N s as the sub-signals 1 ( 1 , 2 , , ) M i iN s where M is an od d number d enoting the length of the window. I n order to ensure that the le ngth of ever y s ub si gnal is M , we need to pad the o riginal signal on the right and left boundaries with ( 1 ) / 2 M zeros. After padding, the ( 1 ) / 2 Mi - th entry o f t he p adding signal is regarded as the ce nter of i s . Then, we truncate the padd ing signal by ex tr acting the center of i s and the ( 1 ) / 2 M entries on both sides of the center from the padding signal. After padding an d truncating , i s is o btained . In order to reduce the influence o f the data far a way from the center of i s and enhance the effect o f data near the center , a Gaussian windo w 1 M g is needed to reweight i s , that is = , i i y g s (3) where the symbol denotes the componentwise multiplication. Suppose 1 N i x is the sparse spectr um with respect to i y , thus, the co rresponding signal in the ti me do main is - 1 1 N i Fx where 1 NN F denotes the inver se F ourier transform matrix. However, t he size o f observation vec tor i y is 1 M i y , thus, a matrix is r equired to establish the relationship bet ween i y and -1 i Fx . -1 NN i Fx is the reconstructed sig nal in ti me do main, which is assumed to an approximation of th e observed m easurement i y . Reg ardi ng i y as the firs t M entries of -1 i Fx , consequentl y, i y and i x can be modeled as -1 , ii y S F x (4) where MN S is the selecting matrix, whic h is defined as , O SI (5) where MM I is the identity matrix and () M N M O is the zero matrix whose e ntries are all zero s. The F in (4) is the Fourier transform matrix defined as 1 2 1 2 4 2 ( 1) 1 2( 1) ( 1) ( 1) 1 1 1 1 1 = 1 , 1 N N N N N N N N N N N N N N N W W W W W W W W W F (6) where 2 exp( ) N Wj N . In a traditional ST FT frame work, after tru ncation and reweighting, the length of i y is M . Ho wever, the size o f the spectrum is N . Therefore , a zero-padding operatio n i s ca rried out o n i y . Therefore, extra frequency components are br ought into STFT . To avoid extra frequenc y compo nents, a reverse idea is created. In fact, the truncation observatio n in ST FT can be regarded as the observation in SR . Comparing the definitio n of with (4), we can fi nd t hat the dictionar y -1 SF . S plays a role of selecti ng the first M entries of -1 i Fx . To illustrate the differences between L1 - norm and Lp - quasinorm, w e show th e co ntour line of Lp -quasinorm i n Fig. 2. As is shown in Fig. 2, the smaller the parameter p is, the sparser the feasible domain of Lp -q uasinorm is . In Fig. 3, we demonstrate t he advantages of Lp -quasinor m. Suppo se the signal is dist urbed by Gaussia n noise (the standard deviation is ) as seen in Fig. 3( a). Obviousl y, the dotted line in Fig. 3(c) would intersect the contour line of Lp -qua sinorm near the axes , inducing a sparser solution ( which is also robust to noise) than the one o f L1-nor m and L2-norm . As is mentioned, Lp - quasinorm may be a better choice. Therefore, we establi sh the model as (7). -1 : min , . . = , p p i i i p P s t SF x y x (7) where i p x is the Lp-nor m defined as 1/ 1 =( ) N p p i p i xx . In the proposed m odel, the local sparse frequency is the inverse objective, while we only focus on t he si milarity between the lo cal measurement and the observation interval of the objec tive in the time domain. Thus, extra freq uency components would not be brought i n the reconstruc ted spectrum. (a) (b) (c) (d) Fig. 2. The contour l ine of Lp -quasinorm . (a) p=2; (b) p=1; (c) p=0.8; (d) p=0.6. 4 B. The solver based on ADMM and LpS operato r The issue sho wn in (7) is a constrained pro blem, which is equivalent to the unconstrained for m shown in (8), 2 2 1 min { }, 2 i p i i i p x x y x (8) where is the coefficient, whic h balances the Lp -quasinor m regularization and t he fidelity ter m 2 2 1 2 ii yx . To solve (8) b y ADMM, the sp litting variab le i zx and the Lagrange coefficie nt are required. Then we have the augmented Lagrangian function as 2 2 2 2 1 max min { 2 , }, 2 i p ii p ii J x , z z y x z x z x (9) where 0 is the penalty coe fficient. For the co nvenience of completing the square in subseq uent calculations, we set = z , then (9) beco mes ( 10 ). 2 2 , 2 2 1 max min { 2 , }. 2 i p ii p i i J xz z z y x z z x z x ( 10 ) Wi th respect to the i x subpro blem , the s ub-objective function is 2 2 ( ) ( ) 2 2 1 m in { + } . 22 i i k k i i i J x x y x z x z ( 11 ) Setting the first-order derivative o f i J x with regard to i x zero , that is, 0 i i J x x , then we have ( 1) 1 ( ) ( ) ( ) ( ( ) ) , k H H k k ii x J y z z ( 12 ) where NN J is the identity matrix. It should be p oint ed out t hat H N N J is a large scale matrix and t he co mputation complexity o f calcu lating 1 () H J directly is 3 () ON . To decrease the computation co mplexity of 1 () H J , the following theorem is required . Theorem 3.1 ： with respect to any matrix NK A , NM B , MK C , the matrix 1 () A B C satisfies the following equation [ 36 ], 1 1 1 1 1 1 ( ) ( ) , A B C A A B I C A B C A ( 13 ) where MM I is the identity matrix. Let = AJ , = H B , = C . Obviously, 1 1 () H I is equal to the expression shown in ( 14 ). 1 1 2 1 1 ( ) = ( ) H H H J J I ( 14 ) Since M is usually much smaller than N , as a result, the computation complexity of 11 () H M M I becomes 3 () OM . Therefore , the co mputation co mplexity of 1 () H J becomes 2 () O N M . In this way, we reduce the co mputation co mplexity of th e whole algorithm remarkably. With respect to the z sub-problem , fixing z and i x , the sub-objective functio n becomes 2 ( ) ( 1) ( 1) 2 min { , } . 2 p k k k i i p J z z z z z x z x ( 15 ) By using the method o f completing square, then we have 2 ( 1 ) ( ) ( +1) ( +1) 2 2 ( ) ( +1) ( +1) 2 ( ) 2 ( ) 2 2 ( +1) ( ) ( ) 2 2 2 ( +1) ( ) 2 arg min { , } 2 = arg min { , 2 + ( ) ( ) } = arg min { ( ) } 2 = arg min { } . 2 p k k k k i i p p k k k i i p kk p k k k i p p k k i p z z z z z z z z x z x z z z x z x zz z z x z z z z x z ( 16 ) The solution of ( 16 ) is ( 1) ( +1 ) ( ) ( + , ) , k k k pi sh ri nk z x z ( 17 ) With r egard to real number, Lp shrinkage operato r is [ 29 ] [ 37 ] 1 2 ( , )= ma x { } . p p p sh rink - ， 0 ( 18 ) (a) (b) (c ) Fig. 3. The fe asible region of Lp -quasinorm. (a) p=2; ( b) p=1; (c) 0
Original Paper
Loading high-quality paper...
Comments & Academic Discussion
Loading comments...
Leave a Comment