Accelerated plane-wave destruction

Accelerated plane-wave destructionZhonghuan Chen 1,Sergey Fomel 2,and Wenkai Lu 1INTRODUCTIONLocal slope fields have been widely used in geophysical applications,such as wave-field separation and denoising (Harlan et al.,1984;Fomel et al.,2007),antialiased seismic interpolation (Bardan,1987),seislet transform (Fomel and Liu,2010),velocity-independent NMO correction and imaging (Fomel,2007b ;Cooke et al.,2009),predictive painting (Fomel,2010),seismic attribute analysis (Marfurt et al.,1999),etc.Several tools exist for local slope estimation:local slant stack (Ottolini,1983;Harlan et al.,1984),complex trace analysis (Barnes,1996),multiwindow dip search (Marfurt,2006),local structure tensor (Fehmers and Höcker,2003;Hale,2007),and plane-wave destruction (Claerbout,1992;Fomel,2002).Plane-wave destruction (PWD)approximates the local wavefield by a local plane wave,and models it using a linear differential equation (Claerbout,1992).When plane-wave destruction is applied on discrete sampled seis-mic signals,the corresponding differential equation needs to be dis-cretized by finite differences.Claerbout (1992)used explicit finite differences.In this method,plane-wave approximation of the wave-field can be seen as applying a linear finite impulse response (FIR)filter to the wavefield.Slope estimation is equivalent to estimating a parameter of the FIR filter.A least-squares estimator of the local slope can be obtained by minimizing the prediction error of the fil-ter.To improve estimation performance of the explicit finite differ-ence scheme,Schleicher et al.(2009)proposed total least-squares estimation.The implicit finite difference scheme was applied to the differ-ential equation by Fomel (2002).Using an infinite impulse response (IIR)filter,known as the “Thiran allpass filter ”(Thiran,1971),to approximate the phase-shift operator,the plane-wave destruction equation becomes a nonlinear equation of the slope.An iterative algorithm was designed to estimate the slope.To improve stability in the iterative algorithm,a smoothing regularization (Fomel,2007a )of the increment can be applied at each iteration.Iterations of regular-ization can be time-consuming,however,particularly in the 3D case.In this paper,we prove that the plane-wave destruction equation is a polynomial equation of an unknown slope.In the case of a three-point approximation of Thiran ’s filter,the convergence results of the iterative algorithm can be analytically analyzed.In this case,we obtain an analytical estimator of the local slope and show that the smoothing regularization can be applied on the final estimator only once.This approach reduces the computational time signifi-cantly.We present 2D and 3D examples,which demonstrate that the proposed algorithm can obtain a slope-estimation result faster than the iterative algorithm,with a similar or even better accuracy.THEORYReview of PWDThe local plane wave can be represented by the following differ-ential equation (Claerbout,1992)Manuscript received by the Editor 23April 2012;revised manuscript received 7August 2012;published online 25January 2013.1Tsinghua University,Department of Automation,State Key Laboratory of Intelligent Technology and Systems,Tsinghua National Laboratory for Informa-tion Science and Technology,Beijing,China.E-mail:zhonghuanchen@;lwkmf@.2The University of Texas at Austin,Bureau of Economic Geology,Jackson School of Geosciences,Austin,Texas,USA.E-mail:sergey.fomel@beg.utexas .edu.©2013Society of Exploration Geophysicists.All rights reserved.V1GEOPHYSICS,VOL.78,NO.1(JANUARY-FEBRUARY 2013);P.V1–V9,9FIGS.,1TABLE.10.1190/GEO2012-0142.1D o w n l o a d e d 09/22/13 t o 221.226.175.186. R e d i s t r i b u t i o n s u b j e c t t o SE G l i c e n s e o r c o p y r i g h t ; s e e T e r m s o f U s e a t h t t p ://l i b r a r y .s e g .o r g /∂u ∂x þσ∂u ∂t¼0;(1)where σis the local slope in continuous space,with dimension of time/length.The wavefields observed at the two positions x 1,x 2have a time delay which is proportional to their distance,σj x 1−x 2j .In the sampled system with space and time intervals Δx and Δt ,we define the discrete space slope in the unit of Δt ∕Δx as p ¼σΔx ∕Δt .Because p is independent of the sampling interval,it can be directly used in irregular data set (in this case,the unit of the slopes becomes space-variant).The time delay between two adjacent positions is then the slope p Δtu ðx;t Þ¼u ðx þΔx;t þp Δt Þ:(2)With the Z transform applied along time and space directions,the above equation becomesð1−Z x Z p t ÞU ðZ x ;Z t Þ¼0;(3)where Z t is the unit time-shift operator,Z x denotes the unit space-shift operator,and U ðZ x ;Z t Þis the Z transform of u ðx;t Þ.Theoperator 1−Z x Z p t is the plane-wave ing Thiran ’sfractional delay filter H ðZ t Þ¼B ð1∕Z t ÞB ðZ t Þ(Thiran,1971)to approximate the time-shift operator Z p t ¼e j ωp ,where ωis the circular frequency,the plane-wave destructor can be expressed as (Fomel,2002),C ðp Þ¼B ðZ t Þ−Z x B 1Z t;(4)whereB ðZ t Þ¼X N k ¼−Nb k ðp ÞZ −k t ;(5)N is the order of the noncausal temporal filter,and b k ðp Þare func-tions of the local slope p .Equation 4is a 2D filter.Applying the filter at an arbitrary point in the wavefield,the plane-wave destruction equation 3becomes a nonlinear equation for the local slope pC ðp;Z x ;Z t ÞU ðZ x ;Z t Þ≈0:(6)An iterative method,such as Newton ’s method,can be applied to find the slope.In practice,wavefields are polluted by noise and the plane wave assumption may not hold true where faults and conflict-ing boundaries exist.To obtain a stable slope estimation,an addi-tional smoothing regularization process (Fomel,2007a )is needed at each step.The total computational cost of slope estimation by plane-wave destruction becomes O ðN d N f N l N n Þ,where N d is the size of the data,N f ¼2N þ1is the size of the filter,N l is the num-ber of linear iterations for regularization,and N n is the number of nonlinear iterations for solving equation 6.Typical values are N f ¼3,5,N l ¼10–50,and N n ¼5–10.Accelerated PWDGauss-Newton ’s iteration searches the solutions for nonlinearequation 6as follows:Let p k be the estimated slope at step k ,with estimating error (or destructive error)e k ¼C ðp k ÞU ðZ x ;Z t Þ.To findthe correct solution p k þ1that minimizes e k þ1,we need to find theincrement Δp k from the local linearizatione k þ1¼C ðp k ÞU ðZ x ;Z t ÞþC 0ðp k ÞU ðZ x ;Z t ÞΔp k ≈0;(7)where C 0ðp k Þis the derivative of C ðp Þat p k with respect to p .The iterative algorithm stops when a stationary point or a root of C ðp ÞU is reached.They are:•Points where C 0ðp k ÞU ¼0:When p k satisfies C 0ðp ÞU ¼0,then e k þ1¼e k ,and the Δp k dependency in equation 7is removed,stopping further iterations on p k .•Points where C ðp k ÞU ¼0and C 0ðp ÞU ≠0:In this case,Δp k ¼0,thus p k þ1¼p k ,eliminating the need for further improvements on p k .The iterative algorithm for equation 6may converge at different points,depending on the initial point that we chose;p 0¼0is a common practical choice for the initial solution.In this case,the iterative algorithm may converge to the least absolute root,which denotes the event with smallest dip angle.To analyze the convergence results,the maximally flat fractional delay filter (Thiran,1971;Zhang,2009)is designed with polyno-mial coefficientsb k ðp Þ¼ð2N Þ!ð2N Þ!ð4N Þ!ðN þk Þ!ðN −k Þ!Y N −1−km ¼0ðm −2N þp Þ×Y N −1þk m ¼0ðm −2N −p Þ:(8)Details on how to design the filter can be found in Appendix A .Because b k ðp Þis a polynomial of p ,expanding it,we get b k ðp Þ¼P 2N i ¼0c ki p iandB ðZ t ;p Þ¼X N k ¼−N X 2N i ¼0c ki Z −k t p i:(9)From equation 8,it is obvious that b k ðp Þ¼b −k ð−p Þ,thereforec ki ¼ð−1Þi c −k;i andB ðZ t ;p Þ¼B1Z t;−p :(10)Substituting the above two equations,the nonlinear equation 6becomes a 2N th°polynomial equation for pX 2N i ¼0a i p i ¼0;(11)and the coefficients of the polynomial plane-wave destruction canbe expressed asa i ¼½1−ð−1ÞiZ xX N k ¼−Nc ki Z −kt U;(12)which says that the coefficients of the polynomial PWD can be ob-tained by applying a 2D filter on the wavefield u .Moreover,the 2D filter can be decoupled into the cascade of two 1D directional filters;the temporal filter P N k ¼−N c ki Z −k t and the spatial filter 1−ð−1Þi Z x .V2Chen et al.D o w n l o a d e d 09/22/13 t o 221.226.175.186. R e d i s t r i b u t i o n s u b j e c t t o SE G l i c e n s e o r c o p y r i g h t ; s e e T e r m s o f U s e a t h t t p ://l i b r a r y .s e g .o r g /In the special case of N ¼1,we get a three-point approximation of B ðZ t Þ.It takes the following form (Fomel,2002)B ðZ t Þ¼ð1þp Þð2þp Þ12Z −1t þð2þp Þð2−p Þ6þð1−p Þð2−p Þ12Z t :(13)The plane-wave destruction equation 11is a quadratic equation.The coefficients a i ði ¼0;1;2Þcan be solved for and expressed asa i ¼112½1−ð−1Þi Z x v i ;(14)where v i are outputs of the following three-point temporal filtersv 0¼2ðZ −1t þ4þZ t ÞU;(15)v 1¼3ðZ t −Z −1t ÞU;(16)andv 2¼ðZ −1t −2þZ t ÞU:(17)In this case,the quadratic plane-wave destruction equation 11has one stationary point and two roots,which can be analytically ex-pressed as:f −a 12a 2;−a 1ÆffiffiffiDp 2a 2g ,where D ¼a 21−4a 0a 2.The plots in Figure 1show the convergence process of the itera-tive algorithm when we choose p 0¼0as the starting value.Geo-metrically when D ≤0,the iteration converges to the stationarypoint −a 12,as shown in Figure 1a .When D >0,it converges to the least absolute solution of equation 11.Figure 1b and 1c shows the convergence process to the least absolute solution in different cases.We can summarize the convergence result of the iterative algorithm as followsp ¼8>><>>:−a12a 2D ≤0−2a0a 1−ffiffiffiD p D >0;a 1<0−2a 0a 1þffiffiffiDpD >0;a 1>0:(18)As −2a 0a 1ffiffiffiD p ¼−a 1ÆffiffiffiDp 2a 2inthe above equation,we use−2a 0a 1ffiffiffiDp instead of −a 1ÆffiffiffiDp 2a 2to obtainbetter numerical stability.When the data is polluted by noise,to obtain a robust slope es-timation,we can combine the equations in a local window into the following equation setFp ≈g ;(19)where F is a normalized diagonal matrix and g is a vector.Their elements are denominators and numerators of equation 18,respec-tively.When we are solving the above equation set by least squares,we can use Tikhonov ’s regularization (Fomel,2002)or the shaping regularization (Fomel,2007a ,equation 13)to obtain a smooth solu-tion as followsp ¼H ½I þH T ðF T F −I ÞH −1H T F T g ;(20)where H is an appropriate smoothing operator.In this case,the reg-ularization runs only once;therefore,the computational cost isreduced to O ðN d N f N l Þ.In 3D applications,there are two polynomial PWD equations for inline and crossline slopes separately.Note that,using the decou-pling,inline,and crossline slope estimations can share the temporal filtering results in equations 15–17.We can obtain the coefficients of the crossline plane-wave destruction equation asa i ¼112½1−ð−1Þi Z y v i .(21)The five-point or longer approximations of B ðZ t Þcan achieve higher accuracy.Equation 6in this case becomes a higher-order polynomial equation (see details in Appendix A ),which can be solved numerically.However,there are multiple stationary points,and it is difficult to determine the right one analytically.For appli-cations that need five-point or higher accuracy,we suggest obtain-ing an initial slope estimation by the proposed three-point method and using it to make the iterative algorithm converge faster (to de-crease N n ).EXAMPLESSynthetic examplesTo test the performance of the proposed slope-estimation method,we generated a harmonic wave field with constant slope p 0¼0.3,shown in Figure 2.We added different scales of additive white Gaussian noise (AWGN)to the wave field,and estimated the slope by the proposed method.To compare with the iterative algorithm by Fomel (2002),the mean square error (MSE)is used as the criterionMSE ðp Þ¼E fðp −p 0Þ2g :(22)We use five iterations in the iterative method,N n ¼5.For the constant slope model,using a large smoothing window,theFigure 1.The iterative results of Newton ’s algorithm.The initial point is p 0¼0(solid circle),and the convergence result is marked by an empty circle:(a)when D ≤0,(b)when D >0and a 1a 2>0,(c)when D >0,and a 1a 2<0.Accelerated plane-wave destructionV3D o w n l o a d e d 09/22/13 t o 221.226.175.186. R e d i s t r i b u t i o n s u b j e c t t o SE G l i c e n s e o r c o p y r i g h t ; s e e T e r m s o f U s e a t h t t p ://l i b r a r y .s e g .o r g /smoothing regularization can converge faster (with less N l )and we can obtain a better estimation accuracy.For each method,we try the smoothing windows from 2to 200and show the mean square errors in Figure pared with the iterative method,the proposed method has better accuracy at the left upper (low SNR and small smoothing window)and worse accuracy at right bottom corner (high SNR and large smoothing window).We show the total runtime of all the noise scale data in Figure 4.For all smoothing windows in the regularization,the proposed method (solid line)only uses about one-fifth of the run time of the three-point (N ¼1)iterative method (dash line).A more complicated model from Claerbout (1999)and Fomel (2002)is shown in Figure 5a .It has variable slopes in synclines,anticlines,and faults.The slope estimated by the proposed method is shown in Figure 5b .The estimation takes about 20ms.The three-point (N ¼1)iterative method can obtain a similar estimation (shown in Figure 5c )by five iterations,which takes about 130ms.Bothmethods use a four-point smoothing window in the regularization,but the proposed method obtains a slightly smoother estimation.In Figure 5d ,we show the faults detected by the residuals of the pro-posed plane-wave destruction.ApplicationThe 2D seislet transform (Fomel and Liu,2010)uses local slopes to predict and update even and odd traces in the wavelet lifting scheme.The seislet transform itself is fast,but the slope estimation step is comparatively slow.The 3D seislet transform can be con-structed in the same way by using 3D slopes and cascading 2D transforms in inline and crossline directions.So that an efficient transform can be built,the proposed accelerated plane-wave de-struction is applied to slope estimation.Figure 6shows a part of the Teapot Dome image from Wyoming.The 3D slope estimation yields two slope fields along twospaceFigure 2.Harmonic waves with constant slope p 0¼0.3.a)b)Figure 3.Mean square errors of the slope estimations by the proposed method (a)and the three-point (N ¼1)iterative method(b).Figure 4.Run time of the proposed method (solid line)and the three-point (N ¼1)iterative method (dashed line).V4Chen et al.D o w n l o a d e d 09/22/13 t o 221.226.175.186. R e d i s t r i b u t i o n s u b j e c t t o SE G l i c e n s e o r c o p y r i g h t ; s e e T e r m s o f U s e a t h t t p ://l i b r a r y .s e g .o r g /directions.Figure7shows the local inline and crossline slope fields estimated by the proposed algorithm.The two slopes are used by the lifting scheme in the seislet trans-form to obtain the seislet transform coefficients.The coefficients of the2D seislet transform are concentrated at the planes near the zero inline plane,as shown in Figure8a,whereas in Figure8b,the3D transform coefficients are concentrated in the corner region near the origin.To illustrate the compressive performance of the3D seislet trans-form,Figure9shows the reconstruction results at different percen-tages of the compression.Most of the data can be reconstructed using only1%of the seislet coefficients(1∶100compression ratio). To compare with the iterative methods quantitatively,we define the following normalized crosscorrelation(NCC)NCCðx;yÞ¼x T yk x k2k y k2:(23)a)c)b)d)Figure5.2D slope estimation example:(a)synthetic data,(b)slope field estimated by the proposed algorithm,(c)slope field estimated by the iterative algorithmafter five iterations,(d)faults detection by plane-wave destruction with the estimated slope.Figure6.A portion of the3D data from the Teapot data set.Accelerated plane-wave destruction V5 Downloaded9/22/13to221.226.175.186.RedistributionsubjecttoSEGlicenseorcopyright;seeTermsofUseathttp://library.seg.org/Figure7.3D slopes estimated by the proposed algorithm:(a)inline slope(b)crossline slope.Figure8.Coefficients of:(a)the2D seislet transform along inline direction only,(b)the3D seislettransform.Figure9.Reconstruction of3D-seislet-compressed data using:(a)5%of the seislet coefficients,(b)1%of the seislet coefficients.V6Chen et al.Downloaded9/22/13to221.226.175.186.RedistributionsubjecttoSEGlicenseorcopyright;seeTermsofUseathttp://library.seg.org/In the seislet transform,the more accurate dip we use,the better compression we can obtain.The NCC between the reconstructed data and the original data can be used to quantify the compressive performance.In Table 1,the proposed method is compared with three-point (N ¼1)and five-point (N ¼2)iterative methods.In all three methods,we use 10-point smoothing windows in in-line and crossline dimensions,and use five-point smoothing win-dow in time direction.To obtain a similar slope estimation as the proposed method,the three-point (N ¼1)and five-point (N ¼2)iterative methods need about six and five iterations,respectively.The five-point method can obtain a better compression ratio than the three-point method because of its better accuracy in slope es-timation.However,although the iterative methods have smaller PWD residuals,neither of them achieves a better NCC than the proposed method.In this case,using the noniterative estimation as the initial in five-point estimation can save at least 250s.That is to say,the computational time cost is reduced by a factor of about six.CONCLUSIONSIn this paper,we derived an analytical estimator of the local slope in three-point implicit plane-wave destruction.On the basis of this result,we built an accelerated slope-estimation algorithm.Exam-ples show that the proposed method can produce a result that is similar to that of the iterative algorithm,at a reduced computation time.Two or more conflicting slopes can be estimated simultaneously by the iterative algorithm.In this case,the PWD equation is a multi-variable polynomial,and the convergence analysis becomes com-plicated.The proposed noniterative method is not yet suitable for multislope estimation.However,we believe that it can find many applications in situations where one dominant slope is suf-ficient.ACKNOWLEDGMENTSWe thank Alexander Klokov for his helpful technical sugges-tions.We also thank the reviewers and the editors for their carefulreview of the paper.The Teapot Dome data is provided by the Rocky Mountain Oilfield Testing Center,sponsored by the U.S.Department of Energy.Zhonghuan Chen ’s joint research is sponsored by the China Scholarship Council and China State Key Science and Technology Projects (2011ZX05023-005-007).APPENDIX APOLYNOMIAL FORM OF PWDIf all the coefficients of B ðZ t Þare polynomials of p ,equation 4is also a polynomial of p ,and the plane-wave destruction equation becomes in turn a polynomial equation of p .The problem is to de-sign a 2N þ1-points filter B ðZ t Þwith polynomial coefficients suchthat the all pass system H ðZ t Þ¼B ð1∕Z t ÞB ðZ t Þcan approximate the phase-shift operator Z p t ¼e j ωp .Denoting the phase response of the system as θðωÞ,that is H ðe j ωÞ¼e j θðωÞ,the group delay of the system isτðωÞ¼∂θðωÞ∂ω:(A-1)The maximally flat criteria designs a filter with a smoothest phase response.There are 2N unknown coefficients in H ðZ t Þ,so we can add 2N flat constraints for the first 2N -th order deviratives of the phase response.It becomes (Zhang,2009,equation 7)τðωÞ¼p ∂n τðωÞ∂ωn ¼0n ¼1;2;:::;2N;(A-2)which is equivalent to the following linear maximally flat condi-tions (Thiran,1971)X N k ¼−Nðd −k Þ2n þ1b k ¼0;(A-3)where n ¼0;1;:::;2N −1and d ¼p ∕2is the fractional delay of B ð1∕Z t Þor 1∕B ðZ t Þ.To solve b k from the above equations,Thiran (1971)used an ad-ditional condition b 0¼1,which leads b k ðk ≠0Þto be a fractional function of p .Differently from that,we use the following condition,X N k ¼−Nb k ¼1;(A-4)where b k can be proved to be polynomials of p .Let vector b ¼½b 0;b N ;:::;b 1;b −1;:::;b −N T .Combining equations A-3and A-4,we rewrite them into the following matrix formTable 1.Performance of different dip estimation methods in 3D seislet transform:Runtime is the runtime of the slope estimation process;RES-inline is the inline residual;RES-xline is the crossline residual;NCC-1is the normalized crosscorrelation between the original data and the reconstruction using 1%of the seislet coefficients;NCC-5uses 5%of the coefficients.The run time for 2D and 3D seislet transform are about 22.45s and 45.0s,respectively.Method N IterationsRuntime (s)RES-inline RES-xline NCC-1NCC-5Noniterative 1054.630.23140.22420.88580.9618Iterative 16334.50.22670.21720.88140.961Iterative25301.60.21940.20750.88380.9615Accelerated plane-wave destruction V7D o w n l o a d e d 09/22/13 t o 221.226.175.186. R e d i s t r i b u t i o n s u b j e c t t o SE G l i c e n s e o r c o p y r i g h t ; s e e T e r m s o f U s e a t h t t p ://l i b r a r y .s e g .o r g /266666411:::11:::1d d −N :::d −1d þ1:::d þN d 3ðd −N Þ3:::ðd −1Þ3ðd þ1Þ3:::ðd þN Þ3......::::::::::::...d 2N −1ðd −N Þ4N −1:::ðd −1Þ4N −1ðd þ1Þ4N −1:::ðd þN Þ2N −13777775×b ¼2666664100 03777775:The matrix on the left side,denoted as V ,can be split into four blocksA B C Das shown above.Following the lemma of matrix inversion,V −1¼ðA −BD −1C Þ−1−ðA −BD −1C Þ−1D −1−D −1C ðA −BD −1C Þ−1D −1þD −1ðA −BD −1C Þ−1BD −1;(A-5)therefore,the coefficientsb ¼V −1½1;0;:::;0 T¼ ðA −BD −1C Þ−1−D −1C ðA −BD −1C Þ−1:(A-6)Let subindex i ¼−N;−N þ1;:::;−1;1;2;:::;N andx i ¼d þi .Submatrix D can be expressed asD ¼EX¼26666641:::11:::1ðd −N Þ2:::ðd −1Þ2ðd þ1Þ2:::ðd −N Þ2ðd −N Þ4:::ðd −1Þ4ðd þ1Þ4:::ðd þN Þ4...:::......:::...ðd −N Þ4N −2:::ðd −1Þ4N −2ðd þ1Þ4N −2:::ðd þN Þ4N −23777775×diag 2666666664x −N ...x −1x 1...x N3777777775;so D −1¼X −1E −1.Denoting U ¼E −1with elements u ij ,j ¼1;2;:::;2N ,as E is a Vandermonde matrix,u ij and Lagrange intepolating polynomials have the following relationshipX 2N j ¼1u ij x 2j −2¼l i ðx Þ;(A-7)where i ¼−N;:::;−1;1;:::;N ,and l i ðx Þis the Lagrange poly-nomial related to the basis d þi ,l i ðx Þ¼Ym ≠i;m ≠0−N ≤m ≤Nx 2−ðd þm Þ2ðd þi Þ2−ðd þm Þ2:(A-8)Substituting the above equation,u ij and x into equation A-7,we can prove equation A-7.It follows that½E −1C i ¼d l i ðd Þ;(A-9)½D −1C i ¼½X −1E −1C i ¼dd þi l iðd Þ;(A-10)withl i ðd Þ¼ð−1Þi þ1N !N !ðN þi Þ!ðN −i Þ!d þi d Y Nm ¼−N 2d þm2d þm þi:(A-11)Thus,A −BD −1C ¼XN i ¼−Nð−1ÞiN !N !ðN þi Þ!ðN −i Þ!Y N m ¼−N p þmp þm þi ¼ð4N Þ!N !N !1Q m ¼N þ1(A-12)and½A −BD −1C −1¼ð2N Þ!ð2N Þ!ð4N Þ!N !N !Y 2Nm ¼N þ1ðm 2−p 2Þ:(A-13)It is the coefficient b 0,a 2N th°polynomial of p .Substituting it into equation A-6,the coefficients at k ¼Æ1;Æ2;:::ÆN are expressed asb k ¼−½D −1C k ½A −BD −1C −1¼ð2N Þ!ð2N Þ!ð4N Þ!ðN þk Þ!ðN −k Þ!Y N −1−km ¼0ðm −2N þp Þ×Y N −1þk m ¼0ðm −2N −p Þ:(A-14)With the additional condition A-4in 2N þ1points approxima-tion,all the coefficients are polynomials of p of 2N th°.Thus,theplane-wave destruction equation 6is proved to be a polynomial equation of 2N th°.REFERENCESBardan,V .,1987,Trace interpolation in seismic data processing:Geo-physical Prospecting,35,343–358,doi:10.1111/j.1365-2478.1987.tb00822.x .Barnes,A.E.,1996,Theory of 2-D complex seismic trace analysis:Geo-physics,61,264–272,doi:10.1190/1.1443947.Claerbout,J.F.,1992,Earth soundings analysis:Processing versus inver-sion:Blackwell Scientific Publications,/sep/prof/index.html ,accessed 10October 2012.V8Chen et al.D o w n l o a d e d 09/22/13 t o 221.226.175.186. R e d i s t r i b u t i o n s u b j e c t t o SE G l i c e n s e o r c o p y r i g h t ; s e e T e r m s o f U s e a t h t t p ://l i b r a r y .s e g .o r g /。