UBC Faculty Research and Publications

Compressed wavefield extrapolation with curvelets Lin, Tim T. Y. 2008

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

Full Text

Compressed wavefield extrapolation with curvelets Tim T.Y. Lin and Felix J. Herrmann SUMMARY An 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 the physics of wave propagation. Because of excessive memory requirements, explicit formulations for wave propagation have proven 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 to compressed 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 which the wavefield is sparse. In this new approach, the eigenfunc- tions of the Helmholtz operator are recognized as a basis that is incoherent with curvelets that are known to compress seismic wavefields. By casting the wavefield extrapolation problem in this framework, wavefields can successfully be extrapolated in the modal domain via a computationally cheaper operatoion. A proof of principle for the “compressed sensing” method is given for wavefield extrapolation in 2-D. The results show that our method is stable and produces identical results compared to the direct application of the full extrapolation operator. INTRODUCTION This paper introduces an explicit algorithm for the compressed evaluation of one-way inverse wavefield extrapolation. The main challenge in this method is the size of the eigenvectors in the modal domain (Grimbergen et al., 1998) that diagonal- ize the one-way wavefield extrapolation operators. We provide a solution strategy that addresses this issues via a “compres- sion” in the size of monochromatic extrapolation operators that are solutions of the one-way wave equation in d dimensions (Claerbout, 1971; Berkhout, 1982; Hale et al., 1992). This is made possible by exploiting the sparsity of seismic wavefields in the curvelet domain (Candès et al., 2006b; Donoho, 2006; Tsaig and Donoho, 2006). As we shall see, this leads to the possibility of using much smaller linear operators in the ex- trapolation of one-way wavefields. Compressed operators The main result from the relatively new field of “compressed sensing” (Candès et al., 2006b; Donoho, 2006; Tsaig and Donoho, 2006) states that an arbitrary k non-zero sparse spike train of lengthN k can exactly be recovered fromm incoherent mea- surements with m ∼ k (∼ means proportional to within logN constants). The term incoherent here refers to a quality be- tween two bases. Qualitatively speaking, a basis is incoherent with respect to another basis if a sparse signal in one of the bases generally does not have a sparse representation in the other basis (see appendix B for a formal definition of this qual- ity). A classical example of two bases that are incoherent with each 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 as y= Ax0 (1) with A := RMST ∈ 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 fromM. Here, S and ST are the sparsity analysis and synthesis matrices for a domain that compresses the signal. The restriction matrix is defined such that the columns of A are 2-norm normalized to unity. The symbol := is used to denote definition. Simply speaking, Eq. 1 corresponds to randomly selecting m Fourier coefficients from the Fourier transform of x0. Spike trains are sparse in the Dirac/identity basis so we set ST := I . The rows of I are incoherent with the rows of the Fourier measurement matrix that consists of complex expo- nentials. Because x0 has few non zeros, it is sparse and this sparse vector can exactly be recovered by solving the follow- ing nonlinear optimization problem ex= argmin x ‖x‖1 = NX i=1 |xi| s.t. Ax= y (2) with the symbole hereby reserved for quantities obtained by solving an optimization problem. The argminx stands for the argument 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 m is large enough compared to the number of non-zero entries in x0. Since m N this recovery involves the inversion of an underdetermined system. As long as the vector x0 is sparse enough, recovery according to Eq. 2 is successful. Typically for Fourier measurements 5 coefficients per non-zero entry are sufficient for full recovery (Candés, 2007). Instead of asking ourselves the question of how to recover x0 from incomplete data suppose now that we ask ourselves how to apply an integer shift by τ to an arbitrary but sparse vector x0, without having to shift each single entry. We all know that shifts translate to phase rotations in the Fourier domain and that the Fourier basis functions (rows of the Fourier matrix,F ) are incoherent with the Dirac basis, I . More formally, consider the approximate shift operation, defined in terms of the expo- nentiation of the discrete difference matrix D ∈RM×M . In that case, the shift by τ can be written as u= e−Dτv= Le jΩτLHv, (3) where the decomposition matrix LH , with the symbol H de- noting the Hermitian transpose, is derived from the eigenvalue problem D= LΩLH . (4) In this expression, Ω is a diagonal matrix with the eigenval- ues ω = diag(Ω) 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 discrete cosine transforms, respectively. The accuracy of this discrete approximation of the shift operator depends on the type and order of the finite-difference approximation in D. Because the eigenvectors of the above shift operation correspond to the rows of Fourier-like (discrete cosine) measurement matrix of the previously posed recovery problem, we can define an alter- native “compressed” procedure for applying the shift by solv- ing the following nonlinear optimization program( y′ = Re jΩτFv= RM′veu= argminu ‖u‖1 s.t. Au= y′ (5) in which we took the liberty to overload the symbolF with the discrete cosine transform. The input for this nonlinear program is given by the phase-rotated Fourier transform of v, restricted to a (small) random set of m frequencies. The symbol ′ is hereby reserved for phase rotated quantities. The shifted spike train is obtained by nonlinear recovery of the phase-rotated measurement vector y. Instead of applying a full matrix-vector multiplication involving all temporal frequency components as in Eq. 3, the shift according to the above program involves the repeated evaluation of the matrix A ∈Cm×N and its transpose. In the extreme case of a vector with a single non-zero entry for 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. An example of the above procedure is included in Fig. 1, where 5 spikes with random positions and amplitudes in a vector of length 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 shows that these results are identical. Only 15 random Fourier mea- surements were necessary for the recovery of the shifted spike train. 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 `1-solver of Basis Pur- suit (Chen et al., 2001). One-way wavefield extrapolation operator Motivated by the above reduction of the shifting operator, we propose that the same technique can be applied to the operator which extrapolates one-way wavefields. After discretization along the space and time axes, the multi-frequency downward wavefield extrapolation of a downgoing time-domain wave- field vector p+|x3 ∈ RM at depthlevel x3 can be written in the form of a matrix-vector product (Grimbergen et al., 1998) u= e jH1∆x3v=Wv (6) with u = p+|x′3 , v = p+|x3 for x3 > x′3. The wavefield vectors contain the reordered entries of discretized wavefield . The propagation over the interval ∆x3 = x3− x′3 is determined by Figure 1: Example of compressed shifting of length 200 with 5 arbitrary spikes. Top: the 5 spikes; Middle: shifted spikes by 20 samples according to Eq. 3. Bottom: the same but ac- cording to the compressed program of Eq. 5. Notice that there is virtually no difference. the square-root of the block-diagonal multi-frequency M×M Helmholtz matrixH2 =H1H1 withM= nν×n f the size of the discretized wavefield, nν the number of total spatial samples and n f the number of angular frequencies. Similar expressions hold for upward extrapolation of an upgoing wavefield. The evaluation of the matrixW ∈ CM×M is expensive since it involves the solution of a sparse eigenvalue problem that di- agonalizes the discretized Helmholtz matrix H2. The solution of this eigenproblem leads to the following factorization of the extrapolation operator W = FHLe jΛ 1/2∆x3LHF (7) = MTM′ (8) with LH the orthonormal modal decomposition matrix. The modal transformmatrixL contains the 2-norm normalized eigen- functions that solve H2 = 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 with the language of compressed sensing, we will call the matrix M := LHF the “measurement” matrix and M′ := e jΛ 1/2∆x3M the phase-rotated measurement matrix. COMPRESSEDWAVEFIELD EXTRAPOLATION This 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 on three main factors: • the compression rate of the to-be-recovered signal eu. This compression rate is quantified by a power-law de- cay rate for the magnitude sorted coefficients |ui∈I | ≤Cri−r with r ≥ 1, (10) with I the indices such that uI(1) ≥ uI(2) ≥ ·· · ≥ uI(N) and r the compression rate with Cr a constant depend- ing on the signal’s energy. The faster the decay the larger r. • themutual coherence between the rows of the measure- ment matrixM and the basis of eu. In the previous spike shift example this is respectively the Fourier transform and the Dirac (identity) basis. However, in order to fa- cilitate the above compression rate criteria, we can also choose to represent our signal in any arbitrary sparsity- promoting basis. Provided that we can find a synthesis matrix ST which transforms the signal to the desired basis, our signal to be recovered is STeu. We can thus generally consider the mutual coherence between M and ST , which is defined by µ(M,S) = √ N max (i, j)∈[1···N]×[1···N] |〈mi, s j〉| (11) with mi and s j the rows ofM and S, respectively. The mutual coherence between the Dirac-Fourier pair is min- imal (µ = 1). The smaller the coherence the fewer ob- servations are required for successful recovery (Candès et al., 2006b; Tsaig and Donoho, 2006; Hennenfent and Herrmann, 2006a). • the randomness of the restriction that is optimal for cases where R selects rows according to a uniform dis- tribution. Depending on the mutual coherence and the compression rate obtained by the sparsity representation, a generalized version of the “compressed operation” program as described by Eq. 5 is able to recover eu for data y′ with large percentages (up to 80% in 2-D) missing. In the case of one-way wavefield extrap- olation, the measurement basisM is LHF , the spectral eigen- basis of the Helmholtz matrix H2. Thus in order to fulfill the mutual coherence and the compression rate requirements we need to find a basis which simultaneously compresses seismic wavefield signals while maintaining a low mutual coherence with 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, multiscale and 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ès and Donoho (2000). Furthermore, it can be qualitatively argued that curvelets hold a very low mutual co- herence with LHF , due to its rapidly decaying nature in the physical domain. Interested readers can refer to (Candès et al., 2006a; Ying et al., 2005) for a more detailed analysis on the behavior of curvelets. Choosing the discrete curvelet transform as ST , we propose to cast the one-way forward wavefield extrapolation into the following nonlinear optimization problem W1 : 8>><>>: y′ = RM′v A := RMCTex= argminx ‖x‖1 s.t. Ax= y′eu= CTex, (12) for an appropriately chosen restriction matrix and for the mea- surement matrixM, defined as in Eq. (7), and with its rotation defined by M′ := e jΛ 1/2∆x3M. (13) In Eq. 12, CT refers to curvelet synthesis by the fast inverse curvelet transform (Candès et al., 2006a; Ying et al., 2005; Hennenfent and Herrmann, 2006b, and Appendix B). During the compressed extrapolation, the wavefield is extrapolated with a compressed forward extrapolation operator bW′=RM′ ∈Cm×M . The forward extrapolated wavefield is subsequently recovered by the nonlinear inversion, promoting sparsity in the curvelet domain. This compressed formulation, which we write in short hand, as eu=W1[v], (14) is designed to yield the same results as for the full forward extrapolation, u=Wv. Restriction strategies The restriction matrix in the compressed formulation for the one-way wavefield extrapolation problem allows for a reduc- tion of the computational cost to synthesize the propagation matrices as well as of their storage requirements. Missed sam- ples of the modal plane can simply be left unconstructed, or ignored in the applying the extrapolation operator. While de- vising a strategy for the restriction, the following issues need to be taken into consideration for the evaluation of the synthesis matrix A′ • the cost of solving the eigenvalue decomposition for each temporal frequency. For a wavefield of size M, i.e. u ∈ RM with n f = O(M1/d) angular frequencies and nν = O(M(d−1)/d) spatial samples, the computa- tion of the eigenfunction decompositions for the sparse Helmholtz matrix requiresO(M(2d−1)/d) operations for the full problem and O(m2ν ×m f ) for the compressed problem. The reduced number of spatial samples equals mν = pν ·O(M(d−1)/d) and the reduced number of fre- quencies m f = p f ·n f with pν , p f < 1 the fractions of the 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 pν · p f ·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 problem and a mν ×m with m = mν ×m f for the compressed problem; • an additional O(N logN) for the computation of the curvelet transform; For d = 3, the size of the wavefield vector grows cubicly (M = n f × n1 × n2), which illustrates the formidable challenge we are 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 to a direct reduction of the number of “block diagonal” eigenproblems that need be solved and to a reduction of the number of eigenvectors that need to be stored; • select a subset of eigenvalues. In principle, discrete eigenfunctions can be calculated at a random subset of discrete eigenvalues by applying the appropriate shifts towards the eigenvalues. This random selection also leads to a reduction of the number of eigenfunctions that need to be stored in memory. • a combination of temporal-frequency and eigenvalue restrictions; EXTRAPOLATION EXAMPLES In this section, we apply our proposed method to wavefield propagation in a lateral varying medium profile, consisting of a single Gaussian shaped low-velocity zone as shown in Fig 2(a). The Helmholtz operator is discretized with ∆x1 = 4m, ∆t = 0.004s with nt = nν = 256. During the experiments, a chain of horizontally oriented fine-scale curvelets is extrapolated down- ward with the propagation direction over a distance of ∆x3 = 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 playing the 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, eu = 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 that a 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 p f and pν ). 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 toW1 (cf. Eq. (12)) for different restrictions. The velocity model cor- responds is plotted in Fig. 2(a). The initial source wavefield v is plotted in Fig. 2(b). (a) The full extrapolated wavefield u = Wv is included for reference; (b) The compressed for- ward propagated wavefield with p f = 0.2 and pν = 0.2; (c) The same as (b) but with p f = 0.4 and pν = 0.4; (d) The same as (b) but with p f = 0.6 and pν = 0.4. Observe that the forward propagated wavefield is largely recovered for the restriction in (c). CONCLUSIONS In 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 the field 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 curvelet frames that are known to compress seismic wavefields. This observation allows for a compression of operators by casting the action of the operator into a compressed sampling and sub- sequent recovery problem. The data is compressively sampled in the modal domain, where the operator is diagonalized, and subsequently recovered with the nonlinear sparsity promoting techniques from compressed sensing. This procedure leads to a reduction of the size of the operators with an overall compu- tational gain that is contingent upon the cost for the recovery. REFERENCES Berkhout, A. J., 1982, Seismic migration. Imaging of acoustic energy by wave field extrapolation: Elsevier, Amsterdam. Candés, E., 2007, Compressive sensing: Presented at the 06/07 IAM Seminars, UBC. Candès, E., L. Demanet, D. Donoho, and L. Ying, 2006a, Fast discrete curvelet transforms: SIAMMultiscale Model. Simul., 5, 861–899. Candès, E., J. Romberg, and T. Tao, 2006b, Stable signal recovery from incomplete and inaccurate measurements: Comm. Pure Appl. Math., 59, 1207–1223. Candès, E. J. and D. L. Donoho, 2000, Recovering Edges in Ill-posed Problems: Optimality of Curvelet Frames: Ann. Statist., 30, 784–842. Chen, S. S., D. L. Donoho, and M. A. Saunders, 2001, Atomic decomposition by basis pursuit: SIAM Journal on Scientific Computing, 43, 129–159. Claerbout, J. F., 1971, Toward a unified theory of reflector mapping: Geophysics, 36, 467–481. Daubechies, I., M. Defrise, and C. de Mol, 2005, An iterative thresholding algorithm for linear inverse problems with a sparsity 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, Modal expansion of one-way operator on laterally varying media: Geophysics, 63, 995–1005. Hale, D., N. R. Hill, and J. Stefani, 1992, Imaging salt with turning 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 at the SEG International Exposition and 76th Annual Meet- ing. ——–, 2006b, Seismic denoising with non-uniformly sampled curvelets: Comp. in Sci. and Eng., 8, 16–25. Tsaig, Y. and D. Donoho, 2006, Extensions of compressed sensing: Signal processing, 86, 549–571. Ying, L., L. Demanet, and E. Candés, 2005, 3d discrete curvelet transform: Presented at the Wavelets XI, SPIE.


Citation Scheme:


Usage Statistics

Country Views Downloads
China 6 7
Japan 6 0
United States 3 0
City Views Downloads
Beijing 6 0
Tokyo 6 0
Sunnyvale 2 0
Unknown 1 0

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


Share to:


Related Items