UBC Faculty Research and Publications

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

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


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

Full Text

Compressed wavefield extrapolation with curvelets Tim T.Y. Lin and Felix J. Herrmann SUMMARY  each other is the identity basis and the Fourier basis.  An explicit algorithm for the extrapolation of one-way wavefields is proposed which combines recent developments in information 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 “compressed sensing”, we are able to formulate the (inverse) wavefield extrapolation problem on small subsets of the data volume, thereby reducing the size of the operators. According to compressed sensing theory, signals can successfully be recovered from an imcomplete set of measurements when the measurement basis is incoherent with the representation in which the wavefield is sparse. In this new approach, the eigenfunctions 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.  This result means that the unknown spike train can, for instance, exactly be recovered from m random Fourier measurements. These measurements correspond to taking inner products between random rows, selected from the Fourier matrix, and the vector containing the sparse spike train. This measurement could explicitly be stated as y = Ax0  (1)  with A := RMST ∈ Cm×N the synthesis matrix, M the measurement 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 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 exponentials. Because x0 has few non zeros, it is sparse and this sparse vector can exactly be recovered by solving the following nonlinear optimization problem  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 diagonalize the one-way wavefield extrapolation operators. We provide a solution strategy that addresses this issues via a “compression” 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`es 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 extrapolation of one-way wavefields. Compressed operators The main result from the relatively new field of “compressed sensing” (Cand`es et al., 2006b; Donoho, 2006; Tsaig and Donoho, 2006) states that an arbitrary k non-zero sparse spike train of length N k can exactly be recovered from m incoherent measurements with m ∼ k (∼ means proportional to within log N constants). The term incoherent here refers to a quality between 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 quality). A classical example of two bases that are incoherent with  e x = arg min x x  1=  N X  |xi | s.t. Ax = y  (2)  i=1  with the symbol e hereby reserved for quantities obtained by solving an optimization problem. The arg minx stands for the argument of the minimum, i.e., the value of the given argument for which the value of the expression attains its minimum value. This recovery is successful when the measurement 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´es, 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 exponentiation of the discrete difference matrix D ∈ RM×M . In that case, the shift by τ can be written as Ωτ H u = e−Dτ v = Le jΩ L v,  (3)  where the decomposition matrix LH , with the symbol H denoting the Hermitian transpose, is derived from the eigenvalue problem ΩLH . D = LΩ (4) In this expression, Ω is a diagonal matrix with the eigenvalΩ) on its diagonal. These eigenvalues correues ω = diag(Ω spond to the angular frequencies, while the orthonormal (de)composition matrices LH , L correspond, when applying Neumann 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 alternative “compressed” procedure for applying the shift by solving the following nonlinear optimization program ( Ωτ F v = RM v y = Re jΩ (5) e u = arg minu u 1 s.t. Au = y in which we took the liberty to overload the symbol F 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. Comparison 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 measurements were necessary for the recovery of the shifted spike train. Instead of applying a full 200 × 200 operator, application of the compressed operator of size 15 × 200 is sufficient. These results were calculated with the 1 -solver of Basis Pursuit (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 wavefield 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 ∆x3 v = Wv p+ |x3 ,  p+ |x3  (6)  with u = v= for x3 > x3 . The wavefield vectors contain the reordered entries of discretized wavefield . The propagation over the interval ∆x3 = x3 − x3 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 according 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 matrix H2 = H1 H1 with M = 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 matrix W ∈ CM×M is expensive since it involves the solution of a sparse eigenvalue problem that diagonalizes the discretized Helmholtz matrix H2 . The solution of this eigenproblem leads to the following factorization of the extrapolation operator W  = =  1/2  Λ F H Le jΛ  ∆x3  LH F  T  M M  (7) (8)  with LH the orthonormal modal decomposition matrix. The modal transform matrix L contains the 2-norm normalized eigenfunctions that solve ΛLH H2 = LΛ (9) with Λ a diagonal matrix with the eigenvalues on the diagonal. From this point on the system of eigenvectors and eigenvalues pertains to the Helmholtz equation. Commensurate with the language of compressed sensing, we will call the matrix Λ1/2 ∆x3 M M := LH F the “measurement” matrix and M := e jΛ the phase-rotated measurement matrix.  COMPRESSED WAVEFIELD EXTRAPOLATION This section describes how the above one-way wavefield extrapolation operator can be formulated in the “compressed operator” framework much like the shift operator above. According to the theory of compressed sensing, successful recovery from the non-linear recovery problem in Eq. 2 depends on three main factors:  • the compression rate of the to-be-recovered signal e u. This compression rate is quantified by a power-law decay rate for the magnitude sorted coefficients |ui∈I | ≤ Cr i−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 depending on the signal’s energy. The faster the decay the larger r. • the mutual coherence between the rows of the measurement matrix M and the basis of e u. In the previous spike shift example this is respectively the Fourier transform and the Dirac (identity) basis. However, in order to facilitate the above compression rate criteria, we can also choose to represent our signal in any arbitrary sparsitypromoting basis. Provided that we can find a synthesis matrix ST which transforms the signal to the desired basis, our signal to be recovered is ST e u. We can thus generally consider the mutual coherence between M and ST , which is defined by √ µ(M, S) = N max | mi , s j | (11) (i, j)∈[1···N]×[1···N]  with mi and s j the rows of M and S, respectively. The mutual coherence between the Dirac-Fourier pair is minimal (µ = 1). The smaller the coherence the fewer observations are required for successful recovery (Cand`es 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 distribution. 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 e u for data y with large percentages (up to 80 % in 2-D) missing. In the case of one-way wavefield extrapolation, the measurement basis M is LH F , the spectral eigenbasis 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 LH F . It has been shown in literature that seismic signals, including wavefields, are highly compressible in the curvelet basis, which expands the wavefield in terms of localized, multiscale and multidirectional prototype waveforms that are anisotropically shaped. As such, the curvelet transform provides a nearoptimal sparsity matrix ST for the purpose of wavefield compression Cand`es and Donoho (2000). Furthermore, it can be qualitatively argued that curvelets hold a very low mutual coherence with LH F , due to its rapidly decaying nature in the physical domain. Interested readers can refer to (Cand`es et al., 2006a; Ying et al., 2005) for a more detailed analysis on the behavior of curvelets. ST ,  Choosing the discrete curvelet transform as we propose to cast the one-way forward wavefield extrapolation into the  following nonlinear optimization problem 8 y = RM v > > > <A := RMCT W1 : > e x = arg minx x 1 s.t. Ax = y > > : e u = CT e x,  (12)  for an appropriately chosen restriction matrix and for the measurement matrix M, defined as in Eq. (7), and with its rotation defined by Λ1/2 ∆x3 M := e jΛ M. (13) In Eq. 12, CT refers to curvelet synthesis by the fast inverse curvelet transform (Cand`es et al., 2006a; Ying et al., 2005; Hennenfent and Herrmann, 2006b, and Appendix B). During the compressed extrapolation, the wavefield is extrapolated with b = RM ∈ Cm×M . a compressed forward extrapolation operator W 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 e u = 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 reduction of the computational cost to synthesize the propagation matrices as well as of their storage requirements. Missed samples of the modal plane can simply be left unconstructed, or ignored in the applying the extrapolation operator. While devising 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(M 1/d ) angular frequencies and nν = O(M (d−1)/d ) spatial samples, the computation of the eigenfunction decompositions for the sparse Helmholtz matrix requires O(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 frequencies 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 wavefield’s discretization; • the cost of applying (repeated) matrix-vector products, which is O(M 2(d−1)/d ) for the full problem and pν · p f · O(M 2(d−1)/d ) for the compressed problem; • the memory use for each frequency amounting to storing 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 log N) 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 extrapolation operators in higher dimensions.  ing the total fractional reduction of the measurement matrix, calculated as the product of p f and pν ). The compressed formulation leads to a gain of in memory usage and more importantly a gain in the cost for computing the operators.  Depending on the requirements (memory imprint versus number 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;  (a)  (b)  (c)  (d)  • 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 = 4 m, ∆t = 0.004 s with nt = nν = 256. During the experiments, a chain of horizontally oriented fine-scale curvelets is extrapolated downward with the propagation direction over a distance of ∆x3 = 600 m perpendicular to the direction in which the medium varies.  Figure 3: Compressed forward extrapolation according to W1 (cf. Eq. (12)) for different restrictions. The velocity model corresponds 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 forward 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  (a)  (b)  Figure 2: (a) Lateral velocity profiles with background velocity 2000 ms−1 and a velocity dip of 1200 ms−1 . Spatial sampling interval of the profile is set to 4 m with 256 samples, while its width at half maximum amplitude is set to 80 m. (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 forward extrapolation operator, u = Wv, and with the nonlinear compressed operator, e u = W1 [v] (cf. Eq. (12)). The initial wavefield v is plotted in Fig. 2(b). The discretized wavefields consist of M = 216 samples. The compressed extrapolation results are obtained with a cooling method outlined in (Daubechies et al., 2005; Elad et al., 2005) which proves effective for L1-optimization problems. The examples show that a fraction p = 0.16 is sufficient to recover the forward propagated wavefields for the Gaussian low and high (with p denot-  In this paper, we proposed a compressed extrapolation algorithm that combines the diagonalization of the Helmholtz operator by the modal transform with recent insights from the field of compressed sensing. In our approach, the eigenfunctions spanning the modal transform are recognized as a suitable 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 subsequent 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 computational 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´es, E., 2007, Compressive sensing: Presented at the 06/07 IAM Seminars, UBC. Cand`es, E., L. Demanet, D. Donoho, and L. Ying, 2006a, Fast discrete curvelet transforms: SIAM Multiscale Model. Simul., 5, 861–899. Cand`es, E., J. Romberg, and T. Tao, 2006b, Stable signal recovery from incomplete and inaccurate measurements: Comm. Pure Appl. Math., 59, 1207–1223. Cand`es, 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. Inform. Theory, 52, 1289–1306. Elad, M., J. Starck, P. Querre, and D. Donoho, 2005, Simulataneous Cartoon and Texture Image Inpainting using Morphological 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. Discussion and reply by authors in GEO-58-8-1205-1206. Hennenfent, G. and F. J. Herrmann, 2006a, Application of stable signal recovery to seismic interpolation: Presented at the SEG International Exposition and 76th Annual Meeting. ——–, 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´es, 2005, 3d discrete curvelet transform: Presented at the Wavelets XI, SPIE.  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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