Open Collections

UBC Faculty Research and Publications

Application of stable signal recovery to seismic interpolation Hennenfent, Gilles; Herrmann, Felix J. 2006-03-05

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

Item Metadata


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

Full Text

Application of stable signal recovery to seismic data interpolationGilles Hennenfent* and Felix J. HerrmannEarth & Ocean Sciences Department, University of British ColumbiaSUMMARYWe propose a method for seismic data interpolation based on 1) thereformulation of the problem as a stable signal recovery problem and2) the fact that seismic data is sparsely represented by curvelets. Thismethod does not require information on the seismic velocities. Mostimportantly, this formulation potentially leads to an explicit recoverycondition. We also propose a large-scale problem solver for the lscript1-regularization minimization involved in the recovery and successfullyillustrate the performance of our algorithm on 2D synthetic and realexamples.INTRODUCTIONSeismic data is often irregularly (and sparsely) sampled along spatialcoordinates due to practical and economical constraints. This problemof irregularly sampled and missing data mainly occurs in 3D seismicsettings although it does occasionally happen in 2D. On land, it can bedue to the presence of a lake, dead or severely contaminated traces, etc.In marine surveys, it can be caused, for example, by the proximity ofa 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 transforminto image artifacts and poor spatial resolution in the migrated image.Relatively simple approaches can be used to handle missing traces inotherwise regularly sampled data. For instance, one can linearly inter-polate neighboring traces. In the more general case of missing tracesin irregularly sampled seismic data, one can collect the traces in binsand stack them. Unfortunately, these approaches do not account forthe presence of wavefronts in the data and consequently these typesof 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 mappingmethods 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 interpolationproblem 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 andUlrych, 1996; Duijndam and Schonewille, 1999; Zwartjes, 2005), thelinear Radon transform (see e.g. Thorson and Claerbout, 1985) and theparabolic 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`et al. (2005b), to seismic data interpolation. This transform-basedmethod explores the sparsity of seismic data in the curvelet domainand does not require information on the seismic velocities. Most im-portantly, it potentially leads to an explicit recovery condition.THEORYSeismic data interpolationSeismic data interpolation is an underdetermined inverse problem. In-terpolation approaches typically involve some sort of a minimizationproblem that aims at jointly minimizing a quadratic distance to the ac-quired data penalized by a functional on the unknown (regularizationterm). The general form for these minimization problems is as followsminx 12bardbly-Axbardbl22 +lambdaJ(x). (1)In this expression, the n-vector y represents the acquired data, lambda theregularization parameter balancing the emphasis of data misfit versusa 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 interpolatedseismic data in the transform domain ? to the acquired data y. Thesuccess of data interpolation depends upon choices for 1) the operatorA and the penalty term J(x), and 2) the regularization parameter lambda.High-resolution methods explore for example the sparsity (minimumstructure) of seismic data in a transform domain. Sacchi and Ulrych(1996) followed later by Zwartjes (2005) use for instance the sparsityof seismic data in the Fourier domain. Although these methods aresuccessful, it is unclear how to assess their performance.Stable signal recoveryConsider a vastly underdetermined system of linear equations. Theproblem is then to recover an unknown m-vector x0 from n lessmuchm con-taminated observations y =Ax0 +n, where A is an n?m matrix and nis an n-vector Gaussian error term of variance sigma2. Such systems haveinfinitely many solutions and one may wonder whether it is possible torecover x0 accurately. Recent results from information theory (Cand`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 thesolution ? of the following minimization problemminx bardblxbardbl1 s.t. bardblAx-ybardbl2 <=<epsilon (2)is within the noise level bardbl? -x0bardbl2 <=<C1 ? epsilon. In this expression, epsilon isthe 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` et al., 2005b). Roughly speaking, the condition on A requiresthat there exists a positive constant S <m such that every subset of S orless columns of A approximately behaves like an orthonormal system.In this case, solving (2) recovers any sparse signal x0 with the numberof 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 asecond orthonormal basis M. 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 Hermitiantranspose). In that case, A := RMSH , where R is a restriction matrixwhich extracts n rows from the m?m orthonormal matrix MSH , whilenormalizing the columns to unit norm. Then the solution ? tobracelefttpbraceexbraceleftmidbraceexbraceleftbtminx bardblxbardbl1 s.t. bardbly-Axbardbl2 <=<epsilon?f =SH ?x(3)verifiesbardbl? -x0bardbl2 <=<C2 ? epsilon +C3 ? bardblx0 -x0,Sbardbl1S (4)providedS <=<C4 ? 1?2 ? n (5)ignoring log-like factors. Eq. (5) is called the recovery condition. Inthese expressions, x0,S is the truncated vector corresponding to the Slargest entries (in absolute value) of x0 =Sf0, (Ci)i=2,3,4 positive wellbehaved constants, and ? the mutual coherence given by?(M, S) =radicalm? maxk,l|mk(sl)H | (6)with mk the kth row of M and sl the lth row of S.This result is significant for several reasons. First of all, it essentiallystates that recovering the S largest entries (in absolute value) of thedecomposition x0 in the sparsity basis S of f0 requires only n proportional S ? ?2observations in the measurement basis M. It thus provides a robustrecovery criterion. Secondly, this criterion directly links the requirednumber of measurements n to the sparsity S of the to-be-recoveredfunction and the mutual coherence ? (dependent on the choice of Sand M). Note that the more similar mk and sl, the larger the mutualcoherence and consequently the larger the number of observations re-quired for the recovery. Thirdly, this result gives an upper bound of theapproximation error which shows that the estimate is almost as goodas if we knew the S-largest entries of x0. In the seismic context, thisframework can thus address important questions such as 1) how wellone can expect to recover seismic data given an acquisition geometry,and 2) what the ?optimal? acquisition geometry is in order to recoverseismic 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/spikebasis is chosen as the measurement basis M because it has a low mu-tual coherence with the sparsity basis DCT. Indeed, there is no vectorfrom S that has a large inner product with any vector of M. Experi-mentally, we measure ?(DCT,Dirac) =radical2. In order to recover x0 andconsequently f0 by solving Eq. (3), 5 random observations of f0 aretypically required (Tsaig and Donoho, 2005). The recovered signal ?(Fig. 1(b)) virtually overlays f0.(a) (b)Figure 1: Perfect recovery of a 1024-sample signal from 5 randomobservations in the Dirac basis. (a) The original signal f0 is randomlymeasured 5 times in the Dirac basis; (b) the recovered signal ? virtuallyoverlays f0.Application of SSR to seismic data interpolationIn the case of missing traces in otherwise regularly sampled data, theacquired 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-calledpicking operator) R on the ideal data. The seismic experiment can thusbe written asy =RMf0 +n. (7)M is the identity/Dirac basis I. Each of its row contains a discreteDirac function that singles out a unique time sample among all possi-ble time samples in the vector f0. R extracts those rows from M thatrepresent time samples actually acquired. R and M are thus fixed bythe acquisition.The recovery condition is a tradeoff between sparsity/compressibilityand mutual coherence. For a given number of observations, an increaseof mutual coherence by a factor of 2 requires the solution to be 4 timessparser in order to still recover it for example. Herrmann and Hennen-fent (2006) show that the emphasis should be put on compressibilityrather than mutual coherence and propose to solve for f0 as a sparsesuperposition of curvelet coefficients (i.e. S defined as the curvelettransform analysis operator C). Curvelets achieve indeed a much bet-ter compressibility of seismic data than Fourier (most incoherent withseismic spatial sampling) at the expense of a moderate loss of inco-herency with the seismic acquisition (see e.g. Cand` et al., 2005a;Herrmann and Hennenfent, 2006, for more details about the curvelettransform and its applications to seismic data processing). The inter-polated data ? is then obtained by ? =CH ?, where? =argminx bardblxbardbl1 s.t. bardbly-RICHbracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightAxbardbl2 <=<epsilon (8)with epsilon the size of the noise present in the acquired data y. This mini-mization problem has the significant advantage over Eq. (1) to relievethe user of testing several regularization parameters lambda.Large-scale problem solver for lscript1-regularization minimizationWe propose a parallelizable (see Thomson et al., 2006) algorithm tosolve Eq. (8). The solver is based on cooling method optimizationand an iterative thresholding algorithm (Daubechies et al., 2004). Thecooling method aims at finding the optimal multiplier lambdaasteriskmath for L(x, lambda) :=lambdabardblxbardbl1 +bardblAx-ybardbl22 -epsilon2, the Lagrangian function of (8), such that theresidual r :=bardblAx-ybardbl2 <=<epsilon. The algorithm is as followsx0 := initial guesslambda0 := initial Lagrange multiplierwhile r >epsilonminx L(x, lambdak)lambdak+1 =alphak lambdak with 0 <alphak <1end while.The critical part of this algorithm is the minimization of L(x, lambdak) doneby the iterative thresholding algorithm presented in Daubechies et al.(2004). At each sub-iteration, evaluation ofxi+1 =Slambdak parenleftbigxi +AH (y-Axi)parenrightbig (9)withSlambdak (x) :=sign(x) ? max(|x| -lambdak, 0) (10)yields an approximate estimate for x which converges to the solutionof the subproblem for a large enough number of iterations. In prac-tice, one only needs to approximately solve each subproblem, whichsignificantly accelerates the overall procedure.RESULTSSynthetic examplesFor the first example, we limit ourselves to 150 iterations (30 updatesof the Lagrange multiplier and 5 sub-iterations). Fig. 2(a) shows thesynthetic dataset with two events, one linear and the other hyperbolic,with 40% randomly-distributed missing traces. The events show anamplitude-versus-offset (AVO) effect. Fig. 2(b) shows the recoveredevents (signal-to-noise ratio (SNR) = 39.2 dB) and Fig. 2(c) the differ-ence between the original and recovered dataset plotted using the samecolor 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 eventsshow an AVO effect. (a) input data with gaps; (b) recovered events;(c) difference between original and recovered events plotted using thesame color scale as (a) and (b); (c) Original AVO effect (plain line)and recovered AVO effect (dash line). Shape and amplitude of eventsare 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 stillpractical for algorithms based on fast transforms. Fig. 3(a) showsthe 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 exampleIn this example, we limit ourselves to 400 iterations (80 updates ofthe Lagrange multiplier and 5 sub-iterations). Fig. 4(a) shows a realmarine dataset, Fig. 4(b) the input data with 40% randomly-distributedmissing traces, Fig. 4(c) the recovered seismic signal (SNR = 20.5 dB),and Fig. 4(d) the difference between the original and recovered datasetplotted 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 the150-180 trace range where the local ratio of missing traces is higher,the overall recovery quality is good. Note again that amplitudes andcontinuity are quite well preserved along wavefronts.CONCLUSIONSWe have presented a method for seismic data interpolation which is anextension to the stable signal recovery as introduced in Cand` et al.(2005b). Instead of resorting to an orthonormal basis, we proposed theredundant curvelet transform as the sparsity representation for seismicdata. This fast-transform-based method ? the curvelet transform is ofsame numerical complexity as the FFT ? exploits the multiscale, multi-directional and continuity properties of seismic wavefronts which leadto sparsity of seismic data in the curvelet domain without requiring in-formation on the seismic velocities. The recovery is carried out with alarge-scale lscript1-regularization minimization problem solver and the ex-(a)(b)(c)(d)Figure 3: Stable seismic signal recovery. (a) synthetic data; (b) inputdata with 20% randomly-distributed traces missing; (c) interpolateddata; (d) difference between (a) and (c) plotted using the same colorscale.(a)(b)(c)(d)Figure 4: Stable seismic signal recovery. (a) real marine data; (b) inputdata with 40% randomly-distributed traces missing; (c) interpolateddata; (d) difference between (a) and (c) plotted using the same colorscale.amples clearly illustrate the performance of our algorithm in the 2Dsource-receiver domain on synthetic and real examples. Extension to3D is straightforward and our methodology may ? as long as the con-tinuity along the wavefronts is preserved (the wavefronts may containcaustics) ? also be integrated in so-called data mapping approaches.ACKNOWLEDGMENTSThe authors would like to thank S. Fomel for his new seismic pro-cessing software RSF, the authors of CurveLab, D.J. (Eric) Verschuurfor the synthetic dataset and ExxonMobil for the real dataset. Thiswork was carried out in part when G.H. was hosted at the UniversidadeFederal do Rio Grande do Norte (UFRN - Natal, Brazil) by S.J. Bez-erra and is part of the SINBAD project with financial support, securedthrough the Industry Technology Facilitator (ITF), from the followingorganizations: BG, BP, Chevron, ExxonMobil and SHELL. Additionalfunding came from the NSERC Discovery Grant 22R81254.REFERENCESCand` E. J., L. Demanet, D. L. Donoho, and L. Ying, 2005a, Curve-lab: Fast Discrete Curvelet Transform.Cand` 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 inaccuratemeasurements.Claerbout, J. F., 1992, Earth soundings analysis: Processing versusinversion: 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 ofband-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-constraineddata continuation with frames: Applications to missing traces andaliased signals in 2/3D: Presented at the 75th SEG Conference &Exhibition.Herrmann, F. J. and G. Hennenfent, 2006, Non-parametric seismicdata recovery with curvelet frames. (to be submitted).Sacchi, M. D. and T. J. Ulrych, 1996, Estimation of the discrete Fouriertransform, 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 DMOoperators: 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 seismicprocessing. (submitted for presentation to the 76th SEG Confer-ence & Exhibition).Thorson, J. R. and J. F. Claerbout, 1985, Velocity-stack and slant-stackstochastic 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: PhDthesis, Delft university of technology.


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