Anti-Aliasing Add-On for Deep Prior Seismic Data Interpolation
Francesco Picetti, Vincenzo Lipari, Paolo Bestagini, Stefano Tubaro
AANTI-ALIASING ADD-ON FOR DEEP PRIOR SEISMIC DATA INTERPOLATION
Francesco Picetti Vincenzo Lipari Paolo Bestagini Stefano Tubaro
Image and Sound Processing Lab (ISPL),Department of Electronics, Information and Bioengineering (DEIB),Politecnico di Milano, Milan, Italy
ABSTRACT
Data interpolation is a fundamental step in any seismic pro-cessing workflow. Among machine learning techniques re-cently proposed to solve data interpolation as an inverse prob-lem, Deep Prior paradigm aims at employing a convolutionalneural network to capture priors on the data in order to reg-ularize the inversion. However, this technique lacks of re-construction precision when interpolating highly decimateddata due to the presence of aliasing. In this work, we pro-pose to improve Deep Prior inversion by adding a directionalLaplacian as regularization term to the problem. This reg-ularizer drives the optimization towards solutions that honorthe slopes estimated from the interpolated data low frequen-cies. We provide some numerical examples to showcase themethodology devised in this manuscript, showing that our re-sults are less prone to aliasing also in presence of noisy andcorrupted data.
Index Terms — spatial aliasing, seismic data, interpola-tion, deep prior, convolutional neural network
1. INTRODUCTION
Economic or environmental constraints, cable feathering inmarine acquisitions and the presence of dead or damagedtraces result in most seismic datasets being poorly and oftenirregularly sampled in space. As a result, the vast majority ofseismic processing and imaging workflows [1, 2, 3] typicallyinclude a preliminary interpolation step. The importanceof the issue is confirmed by the vast number of methodsproposed and the vitality of the relevant scientific literature.Besides traditional processing methods such as transformrepresentations [4, 5, 6], low-rank approximation [7], mul-tichannel singular spectrum analysis (MSSA) [8] and its in-terpolated version I-MSSA [9], many seismic interpolationmethods based on Convolutional Neural Networks (CNNs)have been proposed [10, 11, 12, 13]. However, the vast ma-jority of CNN-based methods work according to a trainingparadigm. In other words, the used CNN learns how to re-construct missing traces during a training step, in which sev-eral pairs of corresponding corrupted and uncorrupted dataare fed to the CNN itself. On one hand, this approach needs a huge amount of training data. On the other hand, CNNs per-formances strongly depend on training examples seen duringtraining, which may prevent the CNN to generalize and workon different kinds of data.A different approach has been proposed interpreting theCNN architecture as a Deep Prior in the framework of inverseproblems to address tasks such as interpolation, denoising orsuper-resolution [14]. Through this approach, the CNN learnsthe inner structure of a 2D image from the corrupted data it-self, without pre-training: this prevents any over-fitting issueas well as the need for training data. This strategy has beenrecently explored also for the task of 2D seismic data interpo-lation using common U-Net as Deep Prior. In particular, DeepPriors have shown effective reconstruction performance whentackling both irregular and regular sampling [15]. However,Deep Prior techniques are not optimally performant when fac-ing highly aliased data.Inspired by our previous studies [16, 17, 18], we proposeto regularize the Deep Prior objective function with a direc-tional Laplacian. Specifically, data orientations are estimatedfrom the low-end spectrum of the CNN output. Then, webuild the directional Laplacian that drives the optimizationtowards solutions that minimizes the energy of events not co-herent with estimated slopes.The rest of the paper is structured as follows. Section 2provides the theoretical formulation of the interpolation prob-lem and the details of the proposed method. Section 3 con-tains the results achieved through our numerical simulationsof three different scenarios. Finally, Section 4 concludes thepaper. The complete code to run our analysis and the datasetsused in this manuscript are accessible to the reader in a repos-itory available at https://bit.ly/3a2ef1t .
2. METHODOLOGY
Seismic data interpolation is an ill-posed inverse problem. Inorder to obtain a reasonable solution it is necessary to con-strain the inversion by adding some kind of a priori infor-mation. A standard way consists in formulating the inverseproblem as finding the solution x ∗ that minimizes an objec- a r X i v : . [ ee ss . SP ] J a n ive function J ( x ) = E ( S x − y obs ) + R ( x ) , (1)where S is a linear sampling operator, y obs are the acquiredseismic data, E ( · ) is the data fidelity term, and R ( · ) is a regu-larizer function formalizing the distance of the solution fromthe space of solutions honoring the desired prior. The choiceof R ( · ) is critical and usually derives from insights of humanexperts.In Deep Prior paradigms, feasible solutions are modeledas the application of a CNN to a random noise realization z .The CNN acts as a parametric non linear function f θ ( · ) andthe problem is recast as finding the set of parameters θ ∗ whichminimizes the objective function J ( θ ) = E ( S f θ ( z ) − y obs ) . (2)In this work, we propose a modification to the objectivecost function (2), adding a regularization term given by theLaplacian operator oriented according to the events slopes.Since data is aliased due to sub-sampling, we need to esti-mate the slopes in a robust way. To this end, we first estimatea low-pass version of the data by applying Deep Prior inter-polation considering the cost function J ( θ ) = E ( S f θ ( z ) − F y obs ) , (3) F being a 1D low-pass filter that acts trace-by-trace, thus re-moving alias effect that could badly impact on standard DeepPrior inversion.By minimizing Eq. (3), we estimate a low-pass versionof the interpolated data. Then, we employ the structure ten-sor algorithm [19] to estimate the slope angles φ and to buildthe directional Laplacian operator −→∇ φ [20]. The latter is em-ployed to regularize the inversion of the broadband seismicdata, allowing us to recast the interpolation problem as theminimization of the cost function J ( θ ) = (cid:107) S f θ ( z ) − y obs (cid:107) + ε w (cid:107) −→∇ φ f θ ( z ) (cid:107) , (4)where w is the vector that collects the anisotropy of the es-timated gradient square tensor for each data sample [19]. Inthe context of seismic data, w can be interpreted as a confi-dence measure of the estimated slopes. Thus, it is added to thecost function to locally weight the directional Laplacian −→∇ φ .Notice that, as the optimization goes on, the interpolated datacan be used to produce a better estimate of the slopes, updat-ing the Laplacian directions at runtime.After the minimization of Eq. (4) has been performed, theinterpolated solution is obtained as x ∗ = f θ ∗ ( z ) . (5) (a) (b) (c) (d) Fig. 1 : Results on synthetic linear events: input decimateddata (a), reference full data (b), results of the proposedmethod (c) and standard Deep Prior optimization (d). (a) (b) (c) (d)
Fig. 2 : FK spectra of synthetic linear events example: inputdecimated data (a), reference full data (b), results of the pro-posed method (c) and standard Deep Prior optimization (d).
3. NUMERICAL RESULTS
In this section, we show numerical results obtained with ourmethodology on three different use cases. The CNN archi-tecture is a standard U-Net design with skip connections.The input noise z is sampled from a Gaussian distribu-tion N (0 , . . Optimization is performed through Adamalgorithm with learning rate . . Low-pass filtering isperformed through a second order Butterworth design. As areconstruction metrics, we compute the total SNR of the inter-polated data with respect to the reference data. We underlinethe fact that the latter is not included in the optimization pro-cess: the network is not fed with the reference full data at anystep. First, we employ a synthetic examples made of 4 linear eventswith different dips. In particular, three events are mildly steep,while the fourth event has been designed to create aliasingeffects. This example is simple yet particularly challengingdue to the slope of the events. The data is made of timesamples with sampling frequency kHz and traces eq-uispaced by m. Data are depicted in Fig. 1b. We build aregular binary mask for simulating the decimation by keeping1 trace over 3, as depicted in Figure 1a. The cutoff frequency f c is Hz; the network is optimized for the first it- a) (b) (c) (d) (e) (f)
Fig. 3 : Results on post-stack field data: input decimated data (a), reference full data (b), results of the proposed method (c) andstandard Deep Prior optimization (d), along with the reconstruction errors (e, f). (a) (b) (c) (d)
Fig. 4 : FK spectra of post-stack field data: input deci-mated data (a), reference full data (b), results of the proposedmethod (c) and standard Deep Prior optimization (d). − − − − (a) . . . . . . (b) Fig. 5 : Slopes φ in degree (a) and related confidence w (b)estimated on the interpolated F3 data.erations to fit the low-passed data, as in Eq. (3). Then, theslopes are estimated and the optimization is performed on thesub-sampled data for iterations more with regularizationweight ε = 5 .Figures 1c and 1d show the results of Deep Prior opti-mization when minimizing the proposed objective (4) andthe standard data fitting, respectively. The achieved SNR is . dB and . dB, respectively. We can notice that thestandard Deep Prior optimization could not resolve the alias-ing of the steeper event, while the proposed methodology ,
000 2 ,
000 3 ,
000 4 ,
000 5 ,
000 6 , − − − − Iterations (cid:107) S f θ ( z ) − y obs (cid:107) (cid:107) −→∇ φ f θ ( z ) (cid:107) Fig. 6 : Objective function components of post-stack field datainterpolation.shows a more homogeneous reconstruction. This effect isparticularly visible also in the spectra depicted in Figure 2d,where we can notice a residual aliasing pattern; on the otherhand, in Figure 2c it has been removed.
Encouraged by this result, we processed a post-stack 2D sec-tion of the F3 Netherlands field data [21]. This section con-sists of time samples, with sampling period ms, and traces with m spacing, as depicted in Figure 3b. Thisexample is quite complex, with different slopes and noiseresiduals coming from previous data processing. The regu-lar binary mask keeps 1 trace over 5, as depicted in Figure 3a.The cutoff frequency f c is Hz; the network is optimizedfor the first iterations to fit the low-passed data, as inEq. (3). Then, the slopes are estimated and the optimizationis performed on the sub-sampled data for iterations morewith regularization weight ε = 0 . . The behavior of the costfunction terms at each iteration can be observed in Figure 6.The big discontinuity in the data fidelity term (in blue) arisesbecause the optimization switches from Eq. (3) to Eq. (4).The orange curve depicts the Laplacian regularizer. For the a) (b) (c) (d) (e) (f) Fig. 7 : Results on pre-stack field data: input decimated data (a), reference full data (b), results of the proposed method (c) andstandard Deep Prior optimization (d), along with the reconstruction errors (e, f). (a) (b) (c) (d)
Fig. 8 : FK spectra of pre-stack field data: input decimateddata (a), reference full data (b), results of the proposedmethod (c) and standard Deep Prior optimization (d).first 1000 iterations it has not been taken into account; then, itis included in the optimization. The discontinuities visible atiteration 2000, 3000 etc. are due to the fact that the slopes φ and confidence w are refined every 1000 iterations from thelow-pass interpolated data.Figures 3c and 3d show the results of our approach anda standard Deep Prior optimization, respectively. Both tech-niques produce a good estimate of the flat zones. On the otherhand, where the interpolated data in Figure 3d is aliased, thedirectional Laplacian regularizer has improved the continu-ity of the steeper reflections. The total SNR of the proposedapproach is . dB, against . dB obtained without regu-larization. Moreover, the Deep Prior paradigm produces esti-mates of the data with less noise content; this is clearly visiblein the spectrum depicted in Figure 4c.Figure 5 depicts the slopes φ and its confidence w es-timated from the low-pass filtered interpolated data. Bysmoothing the gradient with σ = 5 , we are able to discrimi-nate local dips with a good resolution. Finally, we show the limitations of the devised methodologyon pre-stack seismic data taken from the Viking Graeben Line12. This dataset consists of time samples, with sampling period ms, and traces with m spacing. Moreover, thedataset is affected by noise and aliasing, visible in Figure 7band its spectrum in Figure 8b. In order to equalize the ampli-tude of the late events, we apply a linear gain on each trace.The decimation is , as depicted in Figure 7a. The cutofffrequency f c is Hz; the network is optimized for the first iterations to fit the low-passed data, as in Eq. (3). Then,the slopes are estimated and the optimization is performed onthe sub-sampled data for iterations more with regular-ization weight ε = 0 . .This example shows clearly the limitation of the stan-dard Deep Prior inversion. Almost every event is poorly re-constructed, as depicted in Figure 7d. The achieved SNR is . dB. Adopting the proposed technique, data in Figure 7care better interpolated; the total SNR is . dB. The resid-ual depicted in Figure 7e is concentrated almost exclusivelyon the first events. However, a little aliasing effect is stillpresent, as clearly visible in the spectra of Figure 8.
4. CONCLUSIONS
In this manuscript, we propose a Deep Prior method for in-terpolating seismic data strongly affected by aliasing due tohigh regular sub-sampling. First, the interpolation problem iscast in the Deep Prior paradigm. This means we make useof a CNN as an optimizer to minimize a cost function, ratherthan a tool that needs to be trained beforehand. Then, we adda regularization term that relies on directional Laplacian. Thelatter is built upon data slopes estimated from a low-pass ver-sion of the interpolated data produced in a first step.Numerical simulations have been run on synthetic, post-stack and pre-stack data. Results have proven the proposedmethodology to be effective and worth of future efforts, whichwill be devoted to improving the robustness of the inversionand extending this paradigm to 3D data. . REFERENCES [1] W. Chang and G.A. Mcmechan, “3d acoustic prestackreverse-time migration,”
Geophysical Prospecting , vol.38, no. 7, pp. 737–755, 2006.[2] J. Virieux and S. Operto, “An overview of full-waveforminversion in exploration geophysics,”
Geophysics , vol.74, no. 6, pp. WCC1–WCC26, 2009.[3] D.J. Verschuur, A.J. Berkhout, and C.P.A. Wapenaar,“Adaptive surface-related multiple elimination,”
Geo-physics , vol. 57, no. 9, pp. 1166–1177, 1992.[4] S. Spitz, “Seismic trace interpolation in the fx domain,”
Geophysics , vol. 56, no. 6, pp. 785–794, 1991.[5] W. Huang, R. Wu, and R. Wang, “Damped dreamletrepresentation for exploration seismic data interpolationand denoising,”
IEEE Transactions on Geoscience andRemote Sensing (TGRS) , vol. 56, no. 6, pp. 3159–3172,2018.[6] R. Kumar, O. L´opez, D. Davis, A. Y. Aravkin, and F. J.Herrmann, “Beating level-set methods for 5-d seismicdata interpolation: A primal-dual alternating approach,”
IEEE Transactions on Computational Imaging , vol. 3,no. 2, pp. 264–274, 2017.[7] O. L´opez, R. Kumar, ¨O. Yılmaz, and F. J. Herrmann,“Off-the-grid low-rank matrix recovery and seismic datareconstruction,”
IEEE Journal of Selected Topics in Sig-nal Processing (JSTSP) , vol. 10, no. 4, pp. 658–671,2016.[8] M. A. Nazari Siahsar, S. Gholtashi, E. Olyaei Torshizi,W. Chen, and Y. Chen, “Simultaneous denoising and in-terpolation of 3-d seismic data via damped data-drivenoptimal singular value shrinkage,”
IEEE Geoscienceand Remote Sensing Letters , vol. 14, no. 7, pp. 1086–1090, 2017.[9] Fernanda Carozzi and Mauricio D. Sacchi, “Interpolatedmultichannel singular spectrum analysis: A reconstruc-tion method that honors true trace coordinates,”
Geo-physics , vol. 86, no. 1, pp. A1–V89, 2021.[10] Sara Mandelli, Federico Borra, Vincenzo Lipari, PaoloBestagini, Augusto Sarti, and Stefano Tubaro, “Seis-mic data interpolation through convolutional autoen-coder,” in
SEG Technical Program Expanded Abstracts ,pp. 4101–4105. 2018. [11] Ali Siahkoohi, Rajiv Kumar, and F Herrmann, “Seis-mic data reconstruction with generative adversarial net-works,” in , 2018.[12] D.A.B. Oliveira, R.S. Ferreira, R. Silva, and E.V. Brazil,“Interpolating seismic data with conditional generativeadversarial networks,”
IEEE Geoscience and RemoteSensing Letters , vol. 15, no. 12, pp. 1952–1956, 2018.[13] D. Yoon, Z. Yeeh, and J. Byun, “Seismic data recon-struction using deep bidirectional long short-term mem-ory with skip connections,”
IEEE Geoscience and Re-mote Sensing Letters (GRSL) , pp. 1–5, 2020.[14] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Deep imageprior,” in
The IEEE Conference on Computer Vision andPattern Recognition (CVPR) , June 2018.[15] Q. Liu, L. Fu, and M. Zhang, “Deep-seismic-prior-based reconstruction of seismic data using convolutionalneural networks,”
Geophysics , 2019.[16] Fantong. Kong, V. Lipari, P. Bestagini, and S. Tubaro,“A deep prior convolutional autoencoder for seismicdata interpolation,” in , 2020.[17] Fantong Kong, Francesco Picetti, Vincenzo Lipari,Paolo Bestagini, and Stefano Tubaro, “Deep prior-basedseismic data interpolation via multi-res u-net,” in
Soci-ety of Exploration Geophysicists (SEG) Annual Meet-ing , 2020.[18] Fantong Kong, Francesco Picetti, Vincenzo Lipari,Paolo Bestagini, Xiaoming Tang, and Stefano Tubaro,“Deep prior-based unsupervised reconstruction of irreg-ularly sampled seismic data,” 2020.[19] Lucas J. van Vliet and Piet W. Verbeek, “Estimatorsfor orientation and anisotropy in digitized images,” in
Advanced School for Computing and Imaging AnnualConference (ASCI) , 1995.[20] Dave Hale, “Local dip filtering with directional lapla-cians,”