Open Collections

UBC Faculty Research and Publications

Application of stable signal recovery to seismic interpolation Hennenfent, Gilles 2008

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

Item Metadata


hennenfent06seg.pdf [ 280.8kB ]
JSON: 1.0107398.json
JSON-LD: 1.0107398+ld.json
RDF/XML (Pretty): 1.0107398.xml
RDF/JSON: 1.0107398+rdf.json
Turtle: 1.0107398+rdf-turtle.txt
N-Triples: 1.0107398+rdf-ntriples.txt

Full Text

Application of stable signal recovery to seismic data interpolation Gilles Hennenfent* and Felix J. Herrmann Earth & Ocean Sciences Department, University of British Columbia SUMMARY We propose a method for seismic data interpolation based on 1) the reformulation of the problem as a stable signal recovery problem and 2) the fact that seismic data is sparsely represented by curvelets. This method does not require information on the seismic velocities. Most importantly, this formulation potentially leads to an explicit recovery condition. We also propose a large-scale problem solver for the `1- regularization minimization involved in the recovery and successfully illustrate the performance of our algorithm on 2D synthetic and real examples. INTRODUCTION Seismic data is often irregularly (and sparsely) sampled along spatial coordinates due to practical and economical constraints. This problem of irregularly sampled and missing data mainly occurs in 3D seismic settings although it does occasionally happen in 2D. On land, it can be due to the presence of a lake, dead or severely contaminated traces, etc. In marine surveys, it can be caused, for example, by the proximity of a platform or by the cable feathering. As most of the multi-trace pro- cessing algorithms do not handle irregular sampled (and aliased) data, interpolation to a regular grid is necessary in order to process and sub- sequently interpret the data. Indeed, data irregularities often transform into image artifacts and poor spatial resolution in the migrated image. Relatively simple approaches can be used to handle missing traces in otherwise regularly sampled data. For instance, one can linearly inter- polate neighboring traces. In the more general case of missing traces in irregularly sampled seismic data, one can collect the traces in bins and stack them. Unfortunately, these approaches do not account for the presence of wavefronts in the data and consequently these types of solutions to the seismic data interpolation problem do not neces- sarily give optimal results. In the practice of seismic exploration, several methods published in the geophysical literature are applied. These methods can be divided in three main groups. Data mapping methods are based on migration-like operators – e.g. offset contin- uation (see e.g. Stovas and Fomel, 1993), shot continuation (see e.g. Schwab, 1993). Claerbout (1992) and Spitz (1991) describe a sec- ond type of methods based on the formulation of the data interpolation problem as an iterative optimization problem with a convolution op- erator. Finally, there are methods based on particular transforms – e.g. the (non)uniform discrete Fourier transform (see e.g. Sacchi and Ulrych, 1996; Duijndam and Schonewille, 1999; Zwartjes, 2005), the linear Radon transform (see e.g. Thorson and Claerbout, 1985) and the parabolic Radon transform (see e.g. Hampson, 1986). In this paper, we limit ourselves to the case of missing traces in oth- erwise regularly sampled data and improve upon Hennenfent and Her- rmann (2005) by extending the idea of stable signal recovery (SSR) from few and noise-contaminated observations, as introduced in Candès et al. (2005b), to seismic data interpolation. This transform-based method explores the sparsity of seismic data in the curvelet domain and does not require information on the seismic velocities. Most im- portantly, it potentially leads to an explicit recovery condition. THEORY Seismic data interpolation Seismic data interpolation is an underdetermined inverse problem. In- terpolation approaches typically involve some sort of a minimization problem that aims at jointly minimizing a quadratic distance to the ac- quired data penalized by a functional on the unknown (regularization term). The general form for these minimization problems is as follows min x 1 2 ‖y−Ax‖22+λJ(x). (1) In this expression, the n-vector y represents the acquired data, λ the regularization parameter balancing the emphasis of data misfit versus a priori information residing in the penalty term J(x), and A the op- erator – e.g. restriction operator combined with a generic transform – that relates the m-vector x – e.g. representation of the interpolated seismic data in the transform domain – to the acquired data y. The success of data interpolation depends upon choices for 1) the operator A and the penalty term J(x), and 2) the regularization parameter λ . High-resolution methods explore for example the sparsity (minimum structure) of seismic data in a transform domain. Sacchi and Ulrych (1996) followed later by Zwartjes (2005) use for instance the sparsity of seismic data in the Fourier domain. Although these methods are successful, it is unclear how to assess their performance. Stable signal recovery Consider a vastly underdetermined system of linear equations. The problem is then to recover an unknown m-vector x0 from n m con- taminated observations y=Ax0+n, where A is an n×mmatrix and n is an n-vector Gaussian error term of variance σ2. Such systems have infinitely many solutions and one may wonder whether it is possible to recover x0 accurately. Recent results from information theory (Candès et al., 2004, 2005b) show that, if A obeys some specific (weak) con- dition and, if x0 has few nonzero entries (i.e. x0 is sparse), then the solution x̃ of the following minimization problem min x ‖x‖1 s.t. ‖Ax−y‖2 ≤ ε (2) is within the noise level ‖x̃− x0‖2 ≤ C1 · ε . In this expression, ε is the size of the error term n and C1 a positive well behaved (i.e. small) constant. Note that, in the noise-free case, x0 can be recovered exactly (Candès et al., 2005b). Roughly speaking, the condition on A requires that there exists a positive constant S<m such that every subset of S or less columns of A approximately behaves like an orthonormal system. In this case, solving (2) recovers any sparse signal x0 with the number of nonzero entries of about S or less. Similar results were obtained in a more realistic setting where the sig- nal of interest f0 is assumed to be approximately sparse or compress- ible – i.e. sorted absolute values of entries decay like a power law – in some orthonormal sparsity basis S and incompletely measured in a second orthonormal basisM. The measurement – also called observa- tion – y of the vector f0 against the vector m is defined as the projec- tion of f0 onto m – i.e y=mH f0 (the symbol H denotes the Hermitian transpose). In that case, A := RMSH , where R is a restriction matrix which extracts n rows from the m×m orthonormal matrixMSH , while normalizing the columns to unit norm. Then the solution x̃ to  minx ‖x‖1 s.t. ‖y−Ax‖2 ≤ ε f̃= SH x̃ (3) verifies ‖x̃−x0‖2 ≤C2 · ε+C3 · ‖x0−x0,S‖1√ S (4) provided S≤C4 · 1µ2 ·n (5) ignoring log-like factors. Eq. (5) is called the recovery condition. In these expressions, x0,S is the truncated vector corresponding to the S largest entries (in absolute value) of x0 = Sf0, (Ci)i=2,3,4 positive well behaved constants, and µ the mutual coherence given by µ(M,S) = √ m ·max k,l |mk(sl)H | (6) with mk the kth row ofM and sl the lth row of S. This result is significant for several reasons. First of all, it essentially states that recovering the S largest entries (in absolute value) of the decomposition x0 in the sparsity basis S of f0 requires only n ∝ S ·µ2 observations in the measurement basis M. It thus provides a robust recovery criterion. Secondly, this criterion directly links the required number of measurements n to the sparsity S of the to-be-recovered function and the mutual coherence µ (dependent on the choice of S and M). Note that the more similar mk and sl , the larger the mutual coherence and consequently the larger the number of observations re- quired for the recovery. Thirdly, this result gives an upper bound of the approximation error which shows that the estimate is almost as good as if we knew the S-largest entries of x0. In the seismic context, this framework can thus address important questions such as 1) how well one can expect to recover seismic data given an acquisition geometry, and 2) what the ’optimal’ acquisition geometry is in order to recover seismic data within a given accuracy. Fig. 1 illustrates this recovery result with a simple example. The origi- nal 1024-sample signal f0 (Fig. 1(a)) is defined as a single nonzero en- try in the Discrete Cosine Transform (DCT) domain. The Dirac/spike basis is chosen as the measurement basis M because it has a low mu- tual coherence with the sparsity basis DCT. Indeed, there is no vector from S that has a large inner product with any vector of M. Experi- mentally, we measure µ(DCT,Dirac) = √ 2. In order to recover x0 and consequently f0 by solving Eq. (3), 5 random observations of f0 are typically required (Tsaig and Donoho, 2005). The recovered signal f̃ (Fig. 1(b)) virtually overlays f0. (a) (b) Figure 1: Perfect recovery of a 1024-sample signal from 5 random observations in the Dirac basis. (a) The original signal f0 is randomly measured 5 times in the Dirac basis; (b) the recovered signal f̃ virtually overlays f0. Application of SSR to seismic data interpolation In the case of missing traces in otherwise regularly sampled data, the acquired data y is a subgroup of traces of the ideal data f0. Mathemat- ically it translates as the action of a restriction operator (or so-called picking operator)R on the ideal data. The seismic experiment can thus be written as y= RMf0+n. (7) M is the identity/Dirac basis I. Each of its row contains a discrete Dirac function that singles out a unique time sample among all possi- ble time samples in the vector f0. R extracts those rows from M that represent time samples actually acquired. R and M are thus fixed by the acquisition. The recovery condition is a tradeoff between sparsity/compressibility and mutual coherence. For a given number of observations, an increase of mutual coherence by a factor of 2 requires the solution to be 4 times sparser in order to still recover it for example. Herrmann and Hennen- fent (2006) show that the emphasis should be put on compressibility rather than mutual coherence and propose to solve for f0 as a sparse superposition of curvelet coefficients (i.e. S defined as the curvelet transform analysis operator C). Curvelets achieve indeed a much bet- ter compressibility of seismic data than Fourier (most incoherent with seismic spatial sampling) at the expense of a moderate loss of inco- herency with the seismic acquisition (see e.g. Candès et al., 2005a; Herrmann and Hennenfent, 2006, for more details about the curvelet transform and its applications to seismic data processing). The inter- polated data f̃ is then obtained by f̃= CH x̃, where x̃= argmin x ‖x‖1 s.t. ‖y−RICH︸ ︷︷ ︸ A x‖2 ≤ ε (8) with ε the size of the noise present in the acquired data y. This mini- mization problem has the significant advantage over Eq. (1) to relieve the user of testing several regularization parameters λ . Large-scale problem solver for `1-regularization minimization We propose a parallelizable (see Thomson et al., 2006) algorithm to solve Eq. (8). The solver is based on cooling method optimization and an iterative thresholding algorithm (Daubechies et al., 2004). The cooling method aims at finding the optimal multiplier λ ∗ forL (x,λ ) := λ‖x‖1+‖Ax−y‖22−ε2, the Lagrangian function of (8), such that the residual r := ‖Ax−y‖2 ≤ ε . The algorithm is as follows x0 := initial guess λ0 := initial Lagrange multiplier while r> ε minxL (x,λk) λk+1 = αk λk with 0< αk < 1 end while. The critical part of this algorithm is the minimization ofL (x,λk) done by the iterative thresholding algorithm presented in Daubechies et al. (2004). At each sub-iteration, evaluation of xi+1 =Sλk ( xi+AH(y−Axi) ) (9) with Sλk (x) := sign(x) ·max(|x|−λk,0) (10) yields an approximate estimate for x which converges to the solution of the subproblem for a large enough number of iterations. In prac- tice, one only needs to approximately solve each subproblem, which significantly accelerates the overall procedure. RESULTS Synthetic examples For the first example, we limit ourselves to 150 iterations (30 updates of the Lagrange multiplier and 5 sub-iterations). Fig. 2(a) shows the synthetic dataset with two events, one linear and the other hyperbolic, with 40% randomly-distributed missing traces. The events show an amplitude-versus-offset (AVO) effect. Fig. 2(b) shows the recovered events (signal-to-noise ratio (SNR) = 39.2 dB) and Fig. 2(c) the differ- ence between the original and recovered dataset plotted using the same color scale as Fig. 2(a) and Fig. 2(b). There is hardly any recovery er- ror beside a light dimming in the largest gap around trace number 310, which is more visible in the recovered linear event AVO effect com- pared to the original one (Fig. 2(d)). Note that this small amplitude (a) (b) (c) (d) Figure 2: Reconstruction of linear and hyperbolic events from syn- thetic data with 40% randomly-distributed missing traces. The events show an AVO effect. (a) input data with gaps; (b) recovered events; (c) difference between original and recovered events plotted using the same color scale as (a) and (b); (c) Original AVO effect (plain line) and recovered AVO effect (dash line). Shape and amplitude of events are recovered with high fidelity. decline can be corrected by doing more Lagrange multiplier updates – i.e. more iterations. For the second example, we limit ourselves to 400 iterations (80 up- dates of the Lagrange multiplier and 5 sub-iterations), which is still practical for algorithms based on fast transforms. Fig. 3(a) shows the synthetic dataset, Fig. 3(b) the input data with 20% randomly- distributed missing traces, Fig. 3(c) the recovered seismic signal (SNR = 24.4 dB), and Fig. 3(d) the difference between the original and re- covered dataset plotted using the same color scale as Fig. 3(a),Fig. 3(b) and Fig. 3(c). Despite small errors mainly concentrated near the apex, the overall recovery quality is satisfactory. Note especially that ampli- tudes and continuity are well preserved along wavefronts. Field data example In this example, we limit ourselves to 400 iterations (80 updates of the Lagrange multiplier and 5 sub-iterations). Fig. 4(a) shows a real marine dataset, Fig. 4(b) the input data with 40% randomly-distributed missing traces, Fig. 4(c) the recovered seismic signal (SNR = 20.5 dB), and Fig. 4(d) the difference between the original and recovered dataset plotted using the same color scale as Fig. 4(a), Fig. 4(b) and Fig. 4(c). Despite small errors mainly concentrated near the zero offset and in the 150-180 trace range where the local ratio of missing traces is higher, the overall recovery quality is good. Note again that amplitudes and continuity are quite well preserved along wavefronts. CONCLUSIONS We have presented a method for seismic data interpolation which is an extension to the stable signal recovery as introduced in Candès et al. (2005b). Instead of resorting to an orthonormal basis, we proposed the redundant curvelet transform as the sparsity representation for seismic data. This fast-transform-based method – the curvelet transform is of same numerical complexity as the FFT – exploits the multiscale, multi- directional and continuity properties of seismic wavefronts which lead to sparsity of seismic data in the curvelet domain without requiring in- formation on the seismic velocities. The recovery is carried out with a large-scale `1-regularization minimization problem solver and the ex- (a) (b) (c) (d) Figure 3: Stable seismic signal recovery. (a) synthetic data; (b) input data with 20% randomly-distributed traces missing; (c) interpolated data; (d) difference between (a) and (c) plotted using the same color scale. (a) (b) (c) (d) Figure 4: Stable seismic signal recovery. (a) real marine data; (b) input data with 40% randomly-distributed traces missing; (c) interpolated data; (d) difference between (a) and (c) plotted using the same color scale. amples clearly illustrate the performance of our algorithm in the 2D source-receiver domain on synthetic and real examples. Extension to 3D is straightforward and our methodology may – as long as the con- tinuity along the wavefronts is preserved (the wavefronts may contain caustics) – also be integrated in so-called data mapping approaches. ACKNOWLEDGMENTS The authors would like to thank S. Fomel for his new seismic pro- cessing software RSF, the authors of CurveLab, D.J. (Eric) Verschuur for the synthetic dataset and ExxonMobil for the real dataset. This work was carried out in part when G.H. was hosted at the Universidade Federal do Rio Grande do Norte (UFRN - Natal, Brazil) by S.J. Bez- erra and is part of the SINBAD project with financial support, secured through the Industry Technology Facilitator (ITF), from the following organizations: BG, BP, Chevron, ExxonMobil and SHELL. Additional funding came from the NSERC Discovery Grant 22R81254. REFERENCES Candès, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2005a, Curve- lab: Fast Discrete Curvelet Transform. Candès, E. J., J. Romberg, and T. Tao, 2004, Robust uncertainty prin- ciples: Exact signal reconstruction from highly incomplete fre- quency information. ——– 2005b, Stable signal recovery from incomplete and inaccurate measurements. Claerbout, J. F., 1992, Earth soundings analysis: Processing versus inversion: Blackwell Scientific publishing. Daubechies, I., M. Defrise, and C. de Mol, 2004, An iterative thresh- olding algorithm for linear inverse problems with a sparsity con- straint: Comm. Pure Appl. Math., 1413–1457. Duijndam, A. J. W. and M. A. Schonewille, 1999, Reconstruction of band-limited signals, irregularly sampled along spatial direction: Geophysics, 64, 524–538. Hampson, D., 1986, Inverse velocity stacking for multiple elimination: J. Can. Soc. of Expl. Geophys., 22, 44–55. Hennenfent, G. and F. J. Herrmann, 2005, Sparseness-constrained data continuation with frames: Applications to missing traces and aliased signals in 2/3D: Presented at the 75th SEG Conference & Exhibition. Herrmann, F. J. and G. Hennenfent, 2006, Non-parametric seismic data recovery with curvelet frames. (to be submitted). Sacchi, M. D. and T. J. Ulrych, 1996, Estimation of the discrete Fourier transform, a linear inversion approach: Geophysics, 61, 1128– 1136. Schwab, M., 1993, Shot gather continuation: Technical Report 77, SEP. Spitz, S., 1991, Seismic trace interpolation in the f-x domain: Geo- physics, 56, 785–794. Stovas, A. M. and S. B. Fomel, 1993, Kinematically equivalent DMO operators: Russian Geology and Geophysics, 37, 102–113. Thomson, D., G. Hennenfent, H. Modzelewski, and F. J. Herrmann, 2006, A parallel windowed curvelet transform applied to seismic processing. (submitted for presentation to the 76th SEG Confer- ence & Exhibition). Thorson, J. R. and J. F. Claerbout, 1985, Velocity-stack and slant-stack stochastic inversion: Geophysics, 50, 2727–2741. Tsaig, Y. and D. L. Donoho, 2005, Extensions of compressed sensing: Technical report, Stanford university. Zwartjes, P., 2005, Fourier reconstruction with sparse inversion: PhD thesis, Delft university of technology.


Citation Scheme:


Usage Statistics

Country Views Downloads
China 1 0
Japan 1 0
City Views Downloads
Beijing 1 0
Tokyo 1 0

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


Share to:


Related Items