UBC Faculty Research and Publications

Compressed wavefield extrapolation with curvelets Lin, Tim T. Y.; Herrmann, Felix J. 2007

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata


lin2007seg.pdf [ 177.42kB ]
JSON: 1.0107406.json
JSON-LD: 1.0107406+ld.json
RDF/XML (Pretty): 1.0107406.xml
RDF/JSON: 1.0107406+rdf.json
Turtle: 1.0107406+rdf-turtle.txt
N-Triples: 1.0107406+rdf-ntriples.txt
Original Record: 1.0107406 +original-record.json
Full Text

Full Text

Compressed wavefield extrapolation with curveletsTim T.Y. Lin and Felix J. HerrmannSUMMARYAn explicit algorithm for the extrapolation of one-way wave-fields is proposed which combines recent developments in in-formation theory and theoretical signal processing with thephysics of wave propagation. Because of excessive memoryrequirements, explicit formulations for wave propagation haveproven to be a challenge in 3-D. By using ideas from ?com-pressed sensing?, we are able to formulate the (inverse) wave-field extrapolation problem on small subsets of the data vol-ume, thereby reducing the size of the operators. According tocompressed sensing theory, signals can successfully be recov-ered from an imcomplete set of measurements when the mea-surement basis is incoherent with the representation in whichthe wavefield is sparse. In this new approach, the eigenfunc-tions of the Helmholtz operator are recognized as a basis that isincoherent with curvelets that are known to compress seismicwavefields. By casting the wavefield extrapolation problem inthis framework, wavefields can successfully be extrapolated inthe modal domain via a computationally cheaper operatoion.A proof of principle for the ?compressed sensing? method isgiven for wavefield extrapolation in 2-D. The results show thatour method is stable and produces identical results comparedto the direct application of the full extrapolation operator.INTRODUCTIONThis paper introduces an explicit algorithm for the compressedevaluation of one-way inverse wavefield extrapolation. Themain challenge in this method is the size of the eigenvectorsin the modal domain (Grimbergen et al., 1998) that diagonal-ize the one-way wavefield extrapolation operators. We providea solution strategy that addresses this issues via a ?compres-sion? in the size of monochromatic extrapolation operators thatare solutions of the one-way wave equation in d dimensions(Claerbout, 1971; Berkhout, 1982; Hale et al., 1992). This ismade possible by exploiting the sparsity of seismic wavefieldsin the curvelet domain (Cand` et al., 2006b; Donoho, 2006;Tsaig and Donoho, 2006). As we shall see, this leads to thepossibility of using much smaller linear operators in the ex-trapolation of one-way wavefields.Compressed operatorsThe main result from the relatively new field of ?compressedsensing? (Cand` et al., 2006b; Donoho, 2006; Tsaig and Donoho,2006) states that an arbitrary k non-zero sparse spike train oflength N greatermuchk can exactly be recovered from m incoherent mea-surements with m similar k (similar means proportional to within logNconstants). The term incoherent here refers to a quality be-tween two bases. Qualitatively speaking, a basis is incoherentwith respect to another basis if a sparse signal in one of thebases generally does not have a sparse representation in theother basis (see appendix B for a formal definition of this qual-ity). A classical example of two bases that are incoherent witheach other is the identity basis and the Fourier basis.This result means that the unknown spike train can, for in-stance, exactly be recovered from m random Fourier measure-ments. These measurements correspond to taking inner prod-ucts between random rows, selected from the Fourier matrix,and the vector containing the sparse spike train. This measure-ment could explicitly be stated asy =Ax0 (1)with A := RMST element Cm?N the synthesis matrix, M the mea-surement matrix defined in terms of the Fourier matrix (M :=F with F denoting the discrete Fourier transform matrix),and R a restriction matrix randomly selecting m rows from M.Here, S and ST are the sparsity analysis and synthesis matricesfor a domain that compresses the signal. The restriction matrixis defined such that the columns of A are 2-norm normalizedto unity. The symbol := is used to denote definition. Simplyspeaking, Eq. 1 corresponds to randomly selecting m Fouriercoefficients from the Fourier transform of x0.Spike trains are sparse in the Dirac/identity basis so we setST := I. The rows of I are incoherent with the rows of theFourier measurement matrix that consists of complex expo-nentials. Because x0 has few non zeros, it is sparse and thissparse vector can exactly be recovered by solving the follow-ing nonlinear optimization probleme =argminxbardblxbardbl1 =NXi=1|xi| s.t. Ax =y (2)with the symbol ehereby reserved for quantities obtained bysolving an optimization problem. The argminx stands for theargument of the minimum, i.e., the value of the given argu-ment for which the value of the expression attains its mini-mum value. This recovery is successful when the measure-ment and sparsity representations are incoherent and when mis large enough compared to the number of non-zero entriesin x0. Since m lessmuchN this recovery involves the inversion of anunderdetermined system. As long as the vector x0 is sparseenough, recovery according to Eq. 2 is successful. Typicallyfor Fourier measurements 5 coefficients per non-zero entry aresufficient for full recovery (Cand? 2007).Instead of asking ourselves the question of how to recover x0from incomplete data suppose now that we ask ourselves howto apply an integer shift by tau to an arbitrary but sparse vectorx0, without having to shift each single entry. We all know thatshifts translate to phase rotations in the Fourier domain andthat the Fourier basis functions (rows of the Fourier matrix, F)are incoherent with the Dirac basis, I. More formally, considerthe approximate shift operation, defined in terms of the expo-nentiation of the discrete difference matrix D element RM?M. In thatcase, the shift by tau can be written asu =e-Dtau v =LejOmegaOmegaOmegatau LHv, (3)where the decomposition matrix LH, with the symbol H de-noting the Hermitian transpose, is derived from the eigenvalueproblemD =LOmegaLH. (4)In this expression, Omega is a diagonal matrix with the eigenval-ues omega = diag(Omega) on its diagonal. These eigenvalues corre-spond to the angular frequencies, while the orthonormal (de-)composition matrices LH, L correspond, when applying Neu-mann boundary conditions, to the forward and inverse discretecosine transforms, respectively. The accuracy of this discreteapproximation of the shift operator depends on the type andorder of the finite-difference approximation in D. Becausethe eigenvectors of the above shift operation correspond to therows of Fourier-like (discrete cosine) measurement matrix ofthe previously posed recovery problem, we can define an alter-native ?compressed? procedure for applying the shift by solv-ing the following nonlinear optimization program(yprime =RejOmegaOmegaOmegatau Fv =RMprimeve =argminu bardblubardbl1 s.t. Au =yprime (5)in which we took the liberty to overload the symbol F with thediscrete cosine transform. The input for this nonlinear programis given by the phase-rotated Fourier transform of v, restrictedto a (small) random set of m frequencies. The symbol prime ishereby reserved for phase rotated quantities. The shifted spiketrain is obtained by nonlinear recovery of the phase-rotatedmeasurement vector y. Instead of applying a full matrix-vectormultiplication involving all temporal frequency components asin Eq. 3, the shift according to the above program involves therepeated evaluation of the matrix A element Cm?N and its transpose.In the extreme case of a vector with a single non-zero entryfor v, the matrix A will usually only need to be of size 5?N,leading to a significant reduction for the size of the matrix. Anexample of the above procedure is included in Fig. 1, where5 spikes with random positions and amplitudes in a vector oflength N = 200 are circular shifted by 20 samples. Compari-son of the results of applying the full shift operator (cf. Eq. 3)and the compressed shift operator according to Eq. 5 showsthat these results are identical. Only 15 random Fourier mea-surements were necessary for the recovery of the shifted spiketrain. Instead of applying a full 200 ?200 operator, applica-tion of the compressed operator of size 15 ?200 is sufficient.These results were calculated with the lscript1-solver of Basis Pur-suit (Chen et al., 2001).One-way wavefield extrapolation operatorMotivated by the above reduction of the shifting operator, wepropose that the same technique can be applied to the operatorwhich extrapolates one-way wavefields. After discretizationalong the space and time axes, the multi-frequency downwardwavefield extrapolation of a downgoing time-domain wave-field vector p+|x3 element RM at depthlevel x3 can be written in theform of a matrix-vector product (Grimbergen et al., 1998)u =ejH1deltax3 v =Wv (6)with u = p+|xprime3 , v = p+|x3 for x3 > xprime3. The wavefield vectorscontain the reordered entries of discretized wavefield . Thepropagation over the interval deltax3 = x3 -xprime3 is determined byFigure 1: Example of compressed shifting of length 200 with5 arbitrary spikes. Top: the 5 spikes; Middle: shifted spikesby 20 samples according to Eq. 3. Bottom: the same but ac-cording to the compressed program of Eq. 5. Notice that thereis virtually no difference.the square-root of the block-diagonal multi-frequency M ?MHelmholtz matrix H2 =H1H1 with M =nnu ?nf the size of thediscretized wavefield, nnu the number of total spatial samplesand nf the number of angular frequencies. Similar expressionshold for upward extrapolation of an upgoing wavefield.The evaluation of the matrix W element CM?M is expensive since itinvolves the solution of a sparse eigenvalue problem that di-agonalizes the discretized Helmholtz matrix H2. The solutionof this eigenproblem leads to the following factorization of theextrapolation operatorW = FHLej???1/2deltax3 LHF (7)= MT Mprime (8)with LH the orthonormal modal decomposition matrix. Themodal transform matrix L contains the 2-norm normalized eigen-functions that solveH2 =L?LH (9)with ? a diagonal matrix with the eigenvalues on the diagonal.From this point on the system of eigenvectors and eigenval-ues pertains to the Helmholtz equation. Commensurate withthe language of compressed sensing, we will call the matrixM := LHF the ?measurement? matrix and Mprime := ej???1/2deltax3 Mthe phase-rotated measurement matrix.COMPRESSED WAVEFIELD EXTRAPOLATIONThis section describes how the above one-way wavefield ex-trapolation operator can be formulated in the ?compressed op-erator? framework much like the shift operator above. Ac-cording to the theory of compressed sensing, successful recov-ery from the non-linear recovery problem in Eq. 2 depends onthree main factors:? the compression rate of the to-be-recovered signal e.This compression rate is quantified by a power-law de-cay rate for the magnitude sorted coefficients|uielementI| <=<Cri-r with r greaterequal1, (10)with I the indices such that uI(1) greaterequal uI(2) greaterequal ? ? ? greaterequal uI(N)and r the compression rate with Cr a constant depend-ing on the signal?s energy. The faster the decay thelarger r.? the mutual coherence between the rows of the measure-ment matrix M and the basis of e. In the previous spikeshift example this is respectively the Fourier transformand the Dirac (identity) basis. However, in order to fa-cilitate the above compression rate criteria, we can alsochoose to represent our signal in any arbitrary sparsity-promoting basis. Provided that we can find a synthesismatrix ST which transforms the signal to the desiredbasis, our signal to be recovered is ST e. We can thusgenerally consider the mutual coherence between Mand ST , which is defined by?(M, S) =radicalN max(i, j)element[1???N]?[1???N]|angbracketleftmi, sjangbracketright| (11)with mi and sj the rows of M and S, respectively. Themutual coherence between the Dirac-Fourier pair is min-imal (? =1). The smaller the coherence the fewer ob-servations are required for successful recovery (Cand`et al., 2006b; Tsaig and Donoho, 2006; Hennenfent andHerrmann, 2006a).? the randomness of the restriction that is optimal forcases where R selects rows according to a uniform dis-tribution.Depending on the mutual coherence and the compression rateobtained by the sparsity representation, a generalized versionof the ?compressed operation? program as described by Eq. 5is able to recover e for data yprime with large percentages (up to80% in 2-D) missing. In the case of one-way wavefield extrap-olation, the measurement basis M is LHF, the spectral eigen-basis of the Helmholtz matrix H2. Thus in order to fulfill themutual coherence and the compression rate requirements weneed to find a basis which simultaneously compresses seismicwavefield signals while maintaining a low mutual coherencewith LHF.It has been shown in literature that seismic signals, includ-ing wavefields, are highly compressible in the curvelet basis,which expands the wavefield in terms of localized, multiscaleand multidirectional prototype waveforms that are anisotropi-cally shaped. As such, the curvelet transform provides a near-optimal sparsity matrix ST for the purpose of wavefield com-pression Cand` and Donoho (2000). Furthermore, it can bequalitatively argued that curvelets hold a very low mutual co-herence with LHF, due to its rapidly decaying nature in thephysical domain. Interested readers can refer to (Cand` et al.,2006a; Ying et al., 2005) for a more detailed analysis on thebehavior of curvelets.Choosing the discrete curvelet transform as ST , we proposeto cast the one-way forward wavefield extrapolation into thefollowing nonlinear optimization problemW1 :8>>><>>>:yprime =RMprimevA :=RMCTe =argminx bardblxbardbl1 s.t. Ax =yprimee =CT e,(12)for an appropriately chosen restriction matrix and for the mea-surement matrix M, defined as in Eq. (7), and with its rotationdefined byMprime :=ej???1/2deltax3 M. (13)In Eq. 12, CT refers to curvelet synthesis by the fast inversecurvelet transform (Cand` et al., 2006a; Ying et al., 2005;Hennenfent and Herrmann, 2006b, and Appendix B). Duringthe compressed extrapolation, the wavefield is extrapolated witha compressed forward extrapolation operator b prime =RMprime elementCm?M.The forward extrapolated wavefield is subsequently recoveredby the nonlinear inversion, promoting sparsity in the curveletdomain. This compressed formulation, which we write in shorthand, ase =W1[v], (14)is designed to yield the same results as for the full forwardextrapolation, u =Wv.Restriction strategiesThe restriction matrix in the compressed formulation for theone-way wavefield extrapolation problem allows for a reduc-tion of the computational cost to synthesize the propagationmatrices as well as of their storage requirements. Missed sam-ples of the modal plane can simply be left unconstructed, orignored in the applying the extrapolation operator. While de-vising a strategy for the restriction, the following issues need tobe taken into consideration for the evaluation of the synthesismatrix Aprime? the cost of solving the eigenvalue decomposition foreach temporal frequency. For a wavefield of size M,i.e. u element RM with nf = O(M1/d) angular frequenciesand nnu = O(M(d-1)/d) spatial samples, the computa-tion of the eigenfunction decompositions for the sparseHelmholtz matrix requires O(M(2d-1)/d) operations forthe full problem and O(m2nu ?mf ) for the compressedproblem. The reduced number of spatial samples equalsmnu = pnu ? O(M(d-1)/d) and the reduced number of fre-quencies mf = pf ? nf with pnu, pf <1 the fractions ofthe size of the restriction over the total size of the wave-field?s discretization;? the cost of applying (repeated) matrix-vector products,which is O(M2(d-1)/d) for the full problem and pnu ?pf ? O(M2(d-1)/d) for the compressed problem;? the memory use for each frequency amounting to stor-ing a M(d-1)/d ?M(d-1)/d matrix for the full problemand a mnu ?m with m = mnu ?mf for the compressedproblem;? an additional O(N logN) for the computation of thecurvelet transform;For d =3, the size of the wavefield vector grows cubicly (M =nf ?n1 ?n2), which illustrates the formidable challenge weare faced with when designing explicit (inverse) wavefield ex-trapolation operators in higher dimensions.Depending on the requirements (memory imprint versus num-ber of flops), the restriction can be designed to? select a subset of temporal frequencies that leads toa direct reduction of the number of ?block diagonal?eigenproblems that need be solved and to a reductionof the number of eigenvectors that need to be stored;? select a subset of eigenvalues. In principle, discreteeigenfunctions can be calculated at a random subset ofdiscrete eigenvalues by applying the appropriate shiftstowards the eigenvalues. This random selection alsoleads to a reduction of the number of eigenfunctionsthat need to be stored in memory.? a combination of temporal-frequency and eigenvaluerestrictions;EXTRAPOLATION EXAMPLESIn this section, we apply our proposed method to wavefieldpropagation in a lateral varying medium profile, consisting of asingle Gaussian shaped low-velocity zone as shown in Fig 2(a).The Helmholtz operator is discretized with deltax1 = 4m, deltat =0.004s with nt =nnu =256. During the experiments, a chain ofhorizontally oriented fine-scale curvelets is extrapolated down-ward with the propagation direction over a distance of deltax3 =600m perpendicular to the direction in which the medium varies.(a) (b)Figure 2: (a) Lateral velocity profiles with background veloc-ity 2000ms-1 and a velocity dip of 1200ms-1. Spatial sam-pling interval of the profile is set to 4m with 256 samples,while its width at half maximum amplitude is set to 80m. (b)A chain of horizontally-oriented fine-scale curvelets playingthe role of a ?plane wave?.The results for the forward extrapolation are plotted in Figs 3.This figure include the results with the full modal-domain for-ward extrapolation operator, u = Wv, and with the nonlin-ear compressed operator, e = W1[v] (cf. Eq. (12)). The ini-tial wavefield v is plotted in Fig. 2(b). The discretized wave-fields consist of M = 216 samples. The compressed extrapo-lation results are obtained with a cooling method outlined in(Daubechies et al., 2005; Elad et al., 2005) which proves ef-fective for L1-optimization problems. The examples show thata fraction p = 0.16 is sufficient to recover the forward propa-gated wavefields for the Gaussian low and high (with p denot-ing the total fractional reduction of the measurement matrix,calculated as the product of pf and pnu). The compressed for-mulation leads to a gain of in memory usage and more impor-tantly a gain in the cost for computing the operators.(a) (b)(c) (d)Figure 3: Compressed forward extrapolation according to W1(cf. Eq. (12)) for different restrictions. The velocity model cor-responds is plotted in Fig. 2(a). The initial source wavefieldv is plotted in Fig. 2(b). (a) The full extrapolated wavefieldu = Wv is included for reference; (b) The compressed for-ward propagated wavefield with pf = 0.2 and pnu = 0.2; (c)The same as (b) but with pf = 0.4 and pnu = 0.4; (d) Thesame as (b) but with pf = 0.6 and pnu = 0.4. Observe thatthe forward propagated wavefield is largely recovered for therestriction in (c).CONCLUSIONSIn this paper, we proposed a compressed extrapolation algo-rithm that combines the diagonalization of the Helmholtz op-erator by the modal transform with recent insights from thefield of compressed sensing. In our approach, the eigenfunc-tions spanning the modal transform are recognized as a suit-able measurement basis that is incoherent with the curveletframes that are known to compress seismic wavefields. Thisobservation allows for a compression of operators by castingthe action of the operator into a compressed sampling and sub-sequent recovery problem. The data is compressively sampledin the modal domain, where the operator is diagonalized, andsubsequently recovered with the nonlinear sparsity promotingtechniques from compressed sensing. This procedure leads toa reduction of the size of the operators with an overall compu-tational gain that is contingent upon the cost for the recovery.REFERENCESBerkhout, A. J., 1982, Seismic migration. Imaging of acousticenergy by wave field extrapolation: Elsevier, Amsterdam.Cand? E., 2007, Compressive sensing: Presented at the 06/07IAM Seminars, UBC.Cand` E., L. Demanet, D. Donoho, and L. Ying, 2006a,Fast discrete curvelet transforms: SIAM Multiscale Model.Simul., 5, 861?899.Cand` E., J. Romberg, and T. Tao, 2006b, Stable signalrecovery from incomplete and inaccurate measurements:Comm. Pure Appl. Math., 59, 1207?1223.Cand` E. J. and D. L. Donoho, 2000, Recovering Edges inIll-posed Problems: Optimality of Curvelet Frames: Ann.Statist., 30, 784?842.Chen, S. S., D. L. Donoho, and M. A. Saunders, 2001, Atomicdecomposition by basis pursuit: SIAM Journal on ScientificComputing, 43, 129?159.Claerbout, J. F., 1971, Toward a unified theory of reflectormapping: Geophysics, 36, 467?481.Daubechies, I., M. Defrise, and C. de Mol, 2005, An iterativethresholding algorithm for linear inverse problems with asparsity constrains: Comm. Pure Appl. Math., 1413?1457.Donoho, D. L., 2006, Compressed sensing: IEEE Trans. In-form. Theory, 52, 1289?1306.Elad, M., J. Starck, P. Querre, and D. Donoho, 2005, Simula-taneous Cartoon and Texture Image Inpainting using Mor-phological Component Analysis (MCA): Appl. Comput.Harmon. Anal., 19, 340?358.Grimbergen, J., F. Dessing, and C. Wapenaar, 1998, Modalexpansion of one-way operator on laterally varying media:Geophysics, 63, 995?1005.Hale, D., N. R. Hill, and J. Stefani, 1992, Imaging salt withturning seismic waves: Geophysics, 57, 1453?1462. Dis-cussion and reply by authors in GEO-58-8-1205-1206.Hennenfent, G. and F. J. Herrmann, 2006a, Application of sta-ble signal recovery to seismic interpolation: Presented atthe SEG International Exposition and 76th Annual Meet-ing.???, 2006b, Seismic denoising with non-uniformly sampledcurvelets: Comp. in Sci. and Eng., 8, 16?25.Tsaig, Y. and D. Donoho, 2006, Extensions of compressedsensing: Signal processing, 86, 549?571.Ying, L., L. Demanet, and E. Cand? 2005, 3d discretecurvelet transform: Presented at the Wavelets XI, SPIE.


Citation Scheme:


Usage Statistics

Country Views Downloads
United States 11 0
China 8 8
Japan 6 0
Germany 3 0
Poland 1 0
France 1 0
Canada 1 0
City Views Downloads
Tokyo 6 0
Ashburn 6 0
Unknown 6 0
Beijing 6 0
Sunnyvale 3 0
Shenzhen 2 8
Ottawa 1 0
Cambridge 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items