5D seismic data completion and denoising using a novel class of tensor decompositions

14 Pages·2015·2.18 MB·English

5D seismic data completion and denoising using a novel
class of tensor decompositions
The MIT Faculty has made this article openly available. Please share
how this access benefits you. Your story matters.
Citation Ely, Gregory, Shuchin Aeron, Ning Hao, and Misha E. Kilmer.
“5D Seismic Data Completion and Denoising Using a Novel
Class of Tensor Decompositions.” GEOPHYSICS 80, no. 4 (June
8, 2015): V83–V95. © 2015 Society of Exploration Geophysicists
As Published http://dx.doi.org/10.1190/GEO2014-0467.1
Publisher Society of Exploration Geophysicists
Version Final published version
Accessed Wed Jan 23 00:25:01 EST 2019
Citable Link http://hdl.handle.net/1721.1/98498
Terms of Use Article is made available in accordance with the publisher's policy
and may be subject to US copyright law. Please refer to the
publisher's site for terms of use.
Detailed Terms
GEOPHYSICS,VOL.80,NO.4(JULY-AUGUST2015);P.V83–V95,11FIGS.
10.1190/GEO2014-0467.1
g/
r
o
g. 5D seismic data completion and denoising using
e
s
y. a novel class of tensor decompositions
r
a
r
b
p://li
htt
at Gregory Ely1, Shuchin Aeron2, Ning Hao3, and Misha E. Kilmer3
e
s
U
f
o
s
m
er cialandphysicalconstraints.Inparticular,theshotandreceiverpo-
T
e ABSTRACT sitions are planned prior to acquisition but the true sampling
e
s locationsusuallydeviatefromtheoriginallyintendedonesresulting
ht; Wehavedevelopedanovelstrategyforsimultaneousin- in less desirable midpoint and azimuth sampling. Moreover, the
g
ri terpolation and denoising of prestack seismic data. Most midpoint-offset-azimuth domain is also not fully sampled due to
y
p seismic surveys fail to cover all possible source-receiver theintrinsicnatureofseismicacquisitiongeometries.Theseeffects
o
r c combinations,leadingtomissingdataespeciallyinthemid- are more pronounced in land surveys. For methods that rely on
e o point-offset domain. This undersampling can complicate griddeddata,the extentof undersamplingdependsonthe manner
ns certain data processing steps such as amplitude-variation- inwhich the data arebinned. Alsonote thatevenundercomplete
ce with-offset analysis and migration. Data interpolation can acquisitionforagivengeometry,theseismicwavefieldmaystillbe
G li mitigatetheimpactofmissingtraces.Weconsideredthepre- undersampled in the domain in which other processes are carried
E stackdataasa5Dmultidimensionalarrayorotherwisere- out. This spatial undersampling adversely affects later processing
S
o ferredtoasa5Dtensor.Usingsyntheticdatasets,wefirst stepssuchasmigrationandquantitativeinterpretation(QI)studies
ct t foundthatprestackdatacanbewellapproximatedbyalow- such as amplitude versus angle analysis (Sacchi and Liu, 2005;
e ranktensorunderarecentlyproposedframeworkfortensor
bj Huntetal.,2010).Recoveryofprestacktracescanmitigatetheim-
u singular value decomposition (tSVD). Under this low-rank
s pact of incomplete sampling. In this paper, we examine methods
n assumption,weproposedacomplexity-penalizedalgorithm
o based in linear algebra for recovery of missing traces.
uti fortherecoveryofmissingtracesanddatadenoising.Inthis Prestackseismicdatacanbeviewedasafifth-ordertensorconsist-
b algorithm, the complexity regularization wascontrolled by
ri ingofonetimeorfrequencydimensionandfourspatialdimensions
dist ttuesnti.nWgeatseisntgeldetrheeguplearrfiozramtioanncpeaorafmtheeteprruospionsgedaasltgaotirsittihcmal describingthelocationofthesourcesandthereceiversonthesurface,
Re eitherintheoriginalsource-receiverðrx;ry;sx;syÞcoordinateframe
on synthetic and real data to show that missing data can
1.3. be reliably recovered under heavydownsampling. In addi- aoprpinrotaecrmhesso,fsumcihdpasoiinntsTarandd(o2f0fs0e9ts)ðaxn;dyK;hrexi;mhyeÞr.aSnodmSeacocfhthie(2c0u1r2rebn)t,
1. tion,wedemonstratedthatcompressibility,i.e.,approxima-
5 workthef-xdomaininwhichthedatacompletioniscarriedoutsep-
8. tionofthedatabyalow-ranktensor,ofseismicdataunder
o 1 tSVDdependedonthevelocitymodelcomplexityandshot araPterleystoacvkerstheeism4Dicslriecceosncsotrrurecstpioonnditnegchtnoiqeuaecsh ftryepqiuceanllcyy.exploit
15 t acnodrrerleacteeidvewristphatchinegr.eWcoevefruyrthoefrmfoisusnindgthdaattcaobmepcareussseibhiliigtyh the fact that under some assumptions, the fully sampled seismic
4/ wavefieldishighlyredundantandcanbecompressedunderanap-
1 compressibility implied good recovery and vice versa.
9/ propriate representation. Most of these methods fall into two cat-
0
d egories: techniques that explicitly define a transform domain and
e
d methodsthatdonotdefineabasis,insteadrepresentingtheseismic
a
o data as a low-rank matrix or tensor. Numerous methods exploit
nl INTRODUCTION
w compactnessinthefrequency-space(f-x)domain(Duijndametal.,
o
D An ideal seismic survey would cover all possible source and 1999; Hindriks and Duijndam, 2000; Liu and Sacchi, 2004; Xu
receivercombinations.Thesegeometriesrarelyoccurduetofinan- et al., 2005; Zwartjes and Sacchi, 2007; Trad, 2009; Curry,
ManuscriptreceivedbytheEditor1October2014;revisedmanuscriptreceived26March2015;publishedonline8June2015.
1MassachusettsInstituteofTechnology,EarthResourcesLaboratory,Cambridge,Massachusetts,USA.E-mail:[e-mail].
2TuftsUniversity,DepartmentofElectricalEngineering,Medford,Massachusetts,USA.E-mail:[e-mail].
3TuftsUniversity,DepartmentofMathematics,Medford,Massachusetts,USA.E-mail:[e-mail];[e-mail].
©2015SocietyofExplorationGeophysicists.Allrightsreserved.
V83
V84 Ely et al.
2010) and other domains such as the curvelet (Hennenfent and cannotclaimthatourapproachissuperiortotheirs.Extensivecom-
Herrmann, 2006; Herrmann and Hennenfent, 2008; Naghizadeh parison of our approach with those of Kreimer et al. (2013) is a
andSacchi,2010a),seislets(FomelandLiu,2010),andRadon(Ka- subjectoffutureworkandwillbeundertakeninthefuture.InSilva
birandVerschuur,1995).Othermethodsexploitthestatisticalcom- and Herrmann (2013), the authors use the HT decomposition for
pactness or the predictably of a signal in the f-x domain (Spitz, datacompletionwithanovelalgorithmthatexploitsthemanifold
1991;Gulunay,2003;NaghizadehandSacchi,2007,2010b).Inad- structureoftensorswithfixedHTrank.SimilartotheHOSVDcase,
g/ ditiontomethodsthateitherimposesparsityexplicitlyorcompact- theperformanceheavilydependsonthepriorknowledgeofthedi-
or ness throughthresholding inabasis,low-rank matrixmethodsdo mensiontree.Furthermore,itiscomputationallydifficulttosearch
eg. not define a single basis and instead impose compactness on the for the best dimension tree for a given data set. In Trickett et al.
y.s signal’ssubspace.Intheselow-rankmethods,the5Dseismicdata (2013), CANDECOMP (CP)-type factorization is used for data
r
a areembeddedintoa2DHankelmatrixinwhichthetruenoiseless completion.Althoughthisapproachrequiresselectingonlyasingle
r
b
p://li a2n0d10c)o.mIteprlaettievesitghnraelshisoladsisnugmaenddtiombpeosoinfgloswparrasnitky(oTnritchkeesttinegtualla.r, pcoarmampuettaetri,onnaalmlyelcyh,atlhleengCiPngr,aannkd, ifninsdoinmgetchaeseCs,Pthdeecboemstpro-rsaitnioknapis-
htt valuesofthematrixenforcesthereducedrankofthedata,requiring proximations may not even exist (Silva and Lim, 2008).
at numerouscostlysingularvaluedecompositions(SVDs).Othervar- We exploit a novel tensor decomposition, tensor singular value
se iationsenforcethelowrankthroughlessexpensivematrixdecom- decomposition (tSVD), for completion of prestack seismic data.
U
f position, completely avoiding the calculation of the SVD (Gao Kilmer and Martin (2011a) develop tSVD for third-order tensors
s o et al., 2011, 2013; Kumar et al., 2013) or applying the SVD to and Martin et al. (2013) extend it to higher orders. Ely et al.
m a randomized subset of the original matrix (Oropeza and Sacchi, (2013)presenttheinitialresultsofthisapproachforseismicappli-
r
Te 2010, 2011). In contrast to the aforementioned methods, other cations,andZhangetal.(2014)presenttheaccompanyingresults
e methods use a priori velocity information and the wave equation forcomputervisionapplications.For5Dsyntheticseismicdata,we
e
ht; s 2to00i3n;teKrpaoplalaten emtisasl.i,ng20t1ra0c),esth(eRreobnyeni,nc1o9r8p7o;raStitnoglt,ac2t0u0a2l;phFyosmiceall, slohwowintfhoartmthaetiporneasltaccokmspeliesxmitiycdeaxtpareasresecdomasptrheessliobwle;tie.en.s,otrhenyuchlaevaer
g
ri information and principles governing the system. See Stanton norm(TNN)inthetSVDdomain.Thiscompressibilitystemsfrom
y
op and Sacchi (2013) for a detailed survey of current reconstruction theconvolutionalstructureofthet-productbecausethetSVD-based
c
r methods. representationisparticularlyefficientatdescribingshiftandscaling
o
e More recently, methods for seismic completion have exploited operationsbetweentraces.Hence,thedatacanbereliablyrecovered
s
n multilinear algebraic approaches based on tensor rank minimiza- underlimitedsamplingusingtheappropriatecomplexitypenalized
e
G lic ttihoinrd.-Uonrlikheigmheartroicredse,rfoternwsohrisc,hththeerraannkkhiseauvniliyqudeelypednedfsineodn,tfhoer (CPT)healtgSoVriDthdmif.fersfromtheHOSVD,CANDECOMP,andHTfor-
SE choiceoftensorfactorizationandrankpenalizationcanencourage mats because the tSVD retains a tensor’s orientation information
o different structures of the tensor. SVD-like tensor decompositions whileremainingcomputationallyefficient.TheTNNderivedfrom
ct t canbedividedintothreeclasses:canonicaldecomposition(CAN- tSVDdifferssignificantlyfromthesumofnuclearnormsapproach
e
bj DECOMP, CP) (Kolda and Bader, 2009), Tucker or higher-order usedbyGandy etal.(2011)andKreimer etal.(2013),where the
u
s SVD (HOSVD) (Lathauwer et al., 2000), and hierarchical tucker summationisoverthenuclearnormsofallthedifferentmodeun-
n
o (HT) (Grasedyck, 2010). All of these decompositions have been foldingsofthetensor.OurapproachcomputestheTNN,asderived
uti appliedtotheproblemof5Dseismicdatareconstructionwithgood fromthetSVD,byfirstapplyingtheFouriertransformrecursively
b
ri performance on real and synthetic data. Seismic data completion alongthedifferentdimensionsandthensummingthenuclearnorms
st
di using tensor rank minimization under HOSVD is proposed by of the frontal slices of the resulting tensor. In contrast to the ap-
e
R Kreimer and Sacchi (2012b) and Kreimer et al. (2013). Although proachinKreimeretal.(2013),theTNNofatensorascomputed
3. theperformancewasshowntobegood,theproposedalgorithmre- fromthetSVDdependssignificantlyontheorderingofdimensions.
1.1. quirespreselectionofthetruncationranksalongeachdimensionfor ThisisalsoincontrasttotheHOSVDapproachusedbyKreimer
5 effectiveinterpolation.Nosystematicapproachhasbeenproposed and Sacchi (2012a) and the approaches based on HTas used by
8.
1 toselectthesetruncationranks.Onewaytoselectthesetruncation Silva and Herrmann (2013), which are orientation independent.
o
15 t rthaneskespisatroampeerteforsrmfotresthtseornesatsomfathllepdieactae.oAfltthheoduagthatsheitsaanpdptrhoeanchusies Fasorðte;xrxa;mrpy;les,x;twsyoÞatenndsothresoctohnetraiinnidnegxetdheassaðmrxe;rdya;tsax,;soyn;etÞi,nwdeoxueldd
14/ dataadaptive,itisnotcleariftheranksshouldremainconstantover have two different TNNs as computed from tSVD. The approach
9/ theentiredataset.Ifthisisnotthecase,thenitmayrequireground basedonthesumofnuclearnormsofmodeunfoldingsisindepen-
0
d truthandmodelselectionforeachsectionofthedata.Kreimeretal. dent of orientation unless an a priori weighted sum is applied to
e
d (2013) propose seismic data completion based on nuclear norm each of the unfolded matrices. This weighting, when unknown
a
nlo minimization. This approach flattens the data intovarious modes, orsettobeuniform,canimposelow-rankalongfoldingsthatshould
w whichmaybesuboptimalbecauseitdestroysthemultidimensional beofhighrankorviceversa.Similarconsiderationsalsoapplyto
o
D structureofthedata.Recently,Zhangetal.(2014)demonstratethis the cases when using the HOSVD and HT decompositions, in
informationlossfortheproblemof3Dand4Dvideodatacomple- which the choice of truncation ranks is orientation independent.
tionfrommissingpixels.Eventhoughvideodataandseismicdata Weorganizetheremainderofthearticleasfollows:Wefirstout-
areverydifferentdatatypes,weexpectthatretainingthemultidi- line the mathematical problem formulation and review the back-
mensionalstructureofthedatashouldalsohelpwithrecoveryinthe ground of the t-product and tSVD. We first perform a detailed
seismiccase.Nevertheless,itisimportanttonotethatbecausewe casestudy of prestack seismic data compressibility under twoex-
donotcompareourmethodwiththatinKreimeretal.(2013),we tremecases:asetofBornscattersandfinite-differencedatasetsfor
5D seismic tensor completion V85
twodifferentvelocitymodelsandtwosurveygeometries.Wethen Notethatv~⋆u~¼u~⋆v~.Underthismultiplicationoperationandthe
apply the proposed algorithm for a variety of sampling rates and usualvectoraddition,theoperationofatensorX ofsizen ×n ×
1 2
plotthereconstructionperformance.Wethenapplythesimultane- n onM~ ofsizen ×1×n isanotherorientedmatrixZ~ ¼X⋆M~
3 2 3
ousdenoisingandreconstructionalgorithmtoarealdataexample. Pof size n1×1×n3 whose ith tubal element is given by Z~i¼
Finally, we provide conclusions and future research directions. nj¼21Xi;j⋆M~j.ThisisillustratedinFigure1.Onecanextendthis
constructionanddefinitiontodescribethemultiplicationoftwoten-
g/ METHOD sors X andY of sizes n1×n2×n3 and n2×k×n3, respectively,
or resultinginatensorC¼X⋆Y ofsizen1×k×n3.Thismultipli-
g. Notation:Wewilldenotevectorsbyboldfacesmalllettersx,ma- cation between the tensors is referred to as the t-product.
e
y.s trices by boldface capital letters X, and tensors (of order >2) by
ar capitalboldfacescriptlettersX.Throughoutthetext,weuseMAT- Tensor singular value decomposition under the tensor product
r
b LAB notation and indexing to denote the elements and fibers of
at http://li ftdhiebenesoreteivnetthcoteotrhfsier,smtpafargtorein,cteaasln,dsalinXcdeð∶toe;fn1sXÞorw,sX.illFðoid;rejne;ox∶taÞemwtphilleel,fdXiersnðto∶tc;eo∶la;u1mtÞennwsooilrfl ltihneUeatsSrinVogpDethroaeftonaronotivnoe1nr×othfnet2hs×eetnt-o3pfrtooednriuseocntrteaXdndmavsaitesrihwcoeiwnsg,noainn3eDFciagtneunrcseoom2rpaausntdea
e X and so on. given in the following equation:
s
U
of Problem setup X ¼U⋆S⋆V⊤; (4)
s
m
er ThetrueseismicdataMarespatiallyundersampledeitherinthe where U and S are orthogonal tensors of size n1×n1×n3 and
T source-receiverorinthemidpoint-offsetdomain.Themeasurement n ×n ×n , respectively, and S is a diagonal tensor of size
ee YcanberepresentedbyalinearoperatorA,resultinginthesparse n2×n2×n3. The construction of the tSVD has the same form
s 1 2 3
ht; observed data Y under additive noise N as asthematrixsingularvaluedecompositionX¼USVT exceptthat
g t-productandtensor-transposeoperationsreplacethematrixprod-
yri Y ¼AðXÞþN. (1) uctandmatrixtranspose.Atensorcanberegardedasamatrixof
p
o fibersortubesalongthethirddimensionofatensorX.ThetSVDis
c
or Theseismicdatacompletionproblembecomesreliablyestimating analogoustoamatrixSVDifweassumethatthediagonaltensorS
e X fromYunderthesamplingoperatorA.Becausethenumberof consistsofsingular“tubes”or“vectors”onthediagonal,analogous
s
en observedmeasurementsissignificantlylessthanthenumberofel- tosingularvaluesonthediagonalintraditionalSVD.Thecompo-
G lic esomlvenedtsdinirXect,ltyhewiinthveorustepplraocbinlegmciosnssetrvaeirnetlsyoilnlpXos.eAdsanndotceadnninotthbee nIenatndtenVs⊤o⋆rsVU¼anId.VDeofbineiytiothnesofortrhtoegnosnoarltirtyancsopnodsietiðo·nÞ⊤sUan⊤d⋆iUden¼-
E
S introduction,theseconstraintsariseascomplexitymeasuresonthe tity tensor I are given below.Definition 1.
bject to ugnodriethrlmyisngoftrtuheedtyatpaeM.OnecanreliablyrecoverX usingCPal- Xthe⊤Tferinossnotthraeltrsanlin2cse×psonas1ne:d×Ltnhe3etnXterenvbsoeerrsaiontbgetnatshioneredoorfdbesyrizoterfantnr1sap×nosnspi2ons×gedne3af,rcohtnhteoanfl
n su minhðXÞ s:t:kY−AðXÞk≤σn; (2) slices two through n3.Definition 2.
distributio mnwoheiaessereuvrhaeðrsXiathnÞec∈ec.oRNmþopte∪letxfhi0atytgeoiqsfuatihvneaolnternnuetelgydaotaitnvaeeXrceaanalnavdlasloσu2nesdoislmvteahpefoparidnadgcittohivnae-t wfrohBInodteseacenlatfsuiitlrsyisectettsfehrnoeasnroett-rap:zlreosTrldoihcu.eectiisdisetnhdteeitfynin×etednnvsioidarencItiirtc∈yulmRanant×trnci×xon,n3vanoisdluatailoltneo,nthsthoeerr
e strained formulation:
R computation of tSVD can be efficiently calculated using the fast
1.1.3. minkY−AðXÞk s:t:hðXÞ≤c0; (3) FXo,urwieer ftirrasntsafoprpmly(tFhFeTF).FTToalcoonmgptuhtee tthhierdtSdVimDenosfioan3aDndtesntosorer
8.5 for an appropriate choice of the constant c . the Fourier transformed tensor X^. We then compute the standard
1 0 matrix SVD of each frontal slice of the transformed tensor X^
o
5 t The t-product
1
4/
1 Inthissection,wereviewthet-productproposedbyKilmerand
09/ Martin(2011b)andfurtheranalyzedbyBraman(2010)andKilmer
d etal.(2013).Wewillfocuson3Dtensorsforeaseofexpositionand
e
ad interpretation.Webeginbyviewinga3DtensorX ∈Rn1×n2×n3 as
nlo ann1×n2 matrixoftubes(vectorsorientedintothepage),whose
ow i;jthentry,denotedbyXi;j,isgivenbyXi;j≜Xði;j;∶Þ.Onecan
D consideran ×1×n tensorasavectoroftubes.Wewillcallsuch
1 3
tensorsorientedmatricesanddenotethembyM~,withthejthtubal
~
elementdenotedbyMj.Todefinethe3Dtensorasalinearoperator
~
onthe set of oriented matrices M, Braman(2010) defines acom-
mutative multiplication operation between two tubes v~∈R1×1×n3
andu~∈R1×1×n3resultinginanothertubeofsamelength.Thiscom-
mutativeoperationisgivenbycircularconvolutiondenotedby⋆. Figure 1. 3D tensors as operators on oriented matrices.
V86 Ely et al.
andstoreatemporarysliceofthecomponenttensorsU^,S^,andV^. rithm for recovering X from linear measurements by minimizing
Finally,weapplyaninverseFFTtothethirddimensionofthecom- the sum of the elements of the tubal-rank vector of X as
ponenttensorstocomputethefinaltSVDdecomposition.Thiscon-
XN
ceptofthetSVDcanberecursivelyextendedtodefinethet-product CP∶min rankðX^ð∶;∶;iÞÞ s:t:Y ¼AðXÞ. (5)
andtSVDforhigherordertensors(Martinetal.,2013).Forexam-
i¼1
ple,for4Ddata,onecanbeginbyconsideringitasamatrixofma-
g/ trices. Then one defines a commutative operation that maps two However, as in the recovery of undersampled matrices via linear
or matricesofsamesizeandoutputsanothermatrixofthesamesize. operatorsbyminimizingrank,theCPproblemisknowntobecom-
eg. This commutative operation is defined to be the 2D circular con- putationallyinfeasible;seeRechtetal.(2010).Wethereforerelax
s
ry. volution.Operationally,intheextensionoftSVDtohigherdimen- thecomplexPitymeasuretoanormthatwecallastheTNNgivenby
p://libra osaipuoprnlsiae,pdpthlrieeccauFtiroosunivr,ieeplry¼tirna4ndsoifmrorepmn¼siaon5nd. tAhinrlgeveoerrsitthehrmoFuog1uhroiepurt.litNnraeonstsefthoterhmastteafporesr t1keX9n9-k21T)nNaoNlrsm¼okonnoNit¼wh1enkmaX^saðtt∶rhi;ex∶ns;iuinÞcgkleunualcar,rnwvoahrlmeuree.sFkion·rktthnheuecandrogeniusometleeesnsstth(ceWasSaect,showante-,
htt to compute the tSVD for an order-p tensor. Note that in the solve for OPT_TNN(1) and for the noisy case, we solve for
at algorithm, X^ is an order p tensor and the notation X^ð∶;∶;iÞ de- OPT_TNN(2) as outlined below:
Use scribes the ith frontal slice of size n1×n2 of the tensor X^ when
of the tensor is stored in MATLAB format; see Martin et al. (2013) Algorithm 1. tSVD
s for more details.
m
The general t-product framework (including computation of
r
e
e T ttSheVDdi)mreenqsuiiorness ofonretthoeffiixrstananodriesnetcaotinodnionfditchees.teInnsocro,nit.rea.s,tfitxo Input: X ∈Rn1×n2:::×np
ght; se tdheepeonthdeerntt,enthseorfidxeicnogmopfodsiimtioennsmioenthoorddse,rswmhiackheasrteheortiSenVtDatioonrieinn-- XN^ ¼¼nX3n;4%::I:nniptialization
yri tationdependent.Inmanyapplicationswherethetensororder,i.e., for i¼3 to p do
p
o thenumberofdimensions,isnotverylarge,suchastheonecon- X^←fftðX^;½(cid:2);iÞ;
c
or sideredinthispaperwherethetensorsareofordersfourandfive, end for
se onecanfindthebestorientationbasedonsomepriorexperiments for i¼1 to N do
n with synthetic and real data.
ce ½U^;S^;V^(cid:2)¼svdðX^ð∶;∶;iÞÞ
G li Seismic data completion using tensor singular value U^ð∶;∶;iÞ¼U^;S^ð∶;∶;iÞ¼S^
SE decomposition V^ð∶;∶;iÞ¼V^;
o
ect t Based on tSVD, one can define the notion of the tensor rank efonrdifo¼r3 to p do
subj (dKefiilnmeedr ettoalb.,e20a13v).ecTthoer mofulttihraenkraonfkasteonfsotrheusifnrgonttSaVlDslwiceass U←ifftðU^;½(cid:2);iÞ;S←ifftðS^;½(cid:2);iÞ;
on X^ð∶;∶;iÞ;i¼1;2;:::;N, where N¼n3n4:::np (for the nota- V←ifftðV^;½(cid:2);iÞ;
uti tion,seeAlgorithm1).Onecantakethesumoftheranksofthese end for
b
ri frontal slices as acomplexitymeasure. This motivates aCPalgo-
st
di
e
R Figure2. ThetSVDofann ×n ×n tensor.The
3. construction of the tSVD is1sim2ilar to3 the matrix
1. SVDexceptthatthet-productandtensortranspose
51. substitute the equivalent matrix operations. Like n3 n3 n
8. the matrix SVD, the tSVD can also be written as 3
o 1 the sum of outer tensor products. =
15 t n n2
4/ n 1 n
1 1 2
9/
0
d
e n n
ad 2 1 Singular values (tubes)
o
nl n
w 3
o n
D min{n ,n } n 3
= 1 2 2
n
1
j=1
Singular vector (of tubes)
5D seismic tensor completion V87
(cid:2) (cid:3)
XN 1 1
OPTTNNð1Þ∶min kX^ð∶;∶;iÞknuc s:t:Y ¼AðXÞ Zkþ1 ¼arg mZin ρkZkTNNþ2kZ−ðXkþ1þBkÞk2F ;
i¼1
(9)
(6)
and and
org/ XN Bkþ1 ¼BkþðXkþ1−Zkþ1Þ; (10)
g. OPTTNNð2Þ∶minkY−AðXÞk2þλ kX^ð∶;∶;iÞk .
e 2 nuc
rary.s i¼1 (7) fwXhe∶rAe ðρX>Þ¼0,Y1gY¼thAaðtXaÞssiusmethseaivnadliuceatoofr0fufonrctaiollnXofsattihsefyisnegt
b
at http://li dOoPNeTso_ntTeoNtthaNarti(s2ae)l,,thwλoeuisgchotnhisenidrleoearwlits(ehtieesnrmseoifcro)a-rrptaphnlekicsatautkinoeinnosgfthcfoeamcntoopril.seetTelenhseessscp.aroIsne- tthheeTlshianemecaplroliscneogdn-osftoprreamrinattsooarl,nuditsioagnivvtaoelnueeqbuoyafti∞ono8t,hwerhweinseA. correspondsto
e posedoptimizationproblemisconvex(Semercietal.,2014;Zhang
Us etal.,2014)and,therefore,canbesolvedusingexistingtechniques. Xkþ1 ¼ðI−AÞðZk−BkÞþY; (11)
of WewouldliketopointoutagainthattheTNN,asderivedfrom
ms thetSVD,differssignificantlyfromthesumofnuclearnormsap-
Ter proachusedbyGandyetal.(2011)andKreimeretal.(2013).tSVD Algorithm 2. The ADMM algorithm for seismic data
e computes TNN by first applying Fourier transform recursively completion using tSVD
se alongthedifferentdimensionsandthentakingthesumofnuclear
ht; norms of the frontal slices of the resulting tensor, which are N¼
rig n3×n4×n5 in number. In contrast,Kreimer et al. (2013)use the X ¼Z¼B¼0 // Initialize variables.
y
op nuclear norms of matrices formed by mode unfoldings along the N¼n3n4:::np
c five modes of the tensor if operating in the space-time domain,
r while Not Converged do
o or along the four modes if operating in the f-x domain.
e if Constrained then
s AKolmogorovSmirnov(KS)testmethoddeterminesthechoice
cen ofλasdiscussedinthesection“Selectionofregularizationparam- // Projection onto Y¼AðXÞwhen A is the sampling
EG li etitoenr.”prIonbaledmditcioann,btehseonlvoeisdyuasnindgntohiesealletesrsnvateirnsgiodnisreocftitohnemopetthimodizoa-f operaXtor¼ðI−AÞðZ−BÞþY
S
o multipliers (ADMM) methods presented in the following section, else
ct t “Algorithmforsolvingtheoptimizationproblems.”Thesemethods X ¼ðρIþAÞ−1ðρðZ−BÞþYÞ;
bje are also used byKreimer et al. (2013) to solvethe resulting opti- end if
su mizationproblem.ADMMmethodsareideallysuitedtothedown- M¼XþB
n sampling problem because the ADMM class of algorithms
o // Calculate tSVD and shrink in the Fourier domain
buti convergesin a few number ofiterations but requires the typically D¼M
stri ewxaprednosipveeraptroorj.eNctoiotentohfatthteherepsridoujeacltioonntoonthtoenthuellnsuplalcsepoafcethoefftohre- for i¼3 to N do
edi downsampling operator is very inexpensive to compute. If, how- D←fftðD;½(cid:2);iÞ;
R
3. ever,thedatawereacquiredusingagenericsamplingoperator,such end for
1. asinthecaseinwhichthetraceswererandomlysummedorsub- // shrinking via singular value thresholding on each slice
8.51. tcroamctpedutfer.oImnothneesaencoatsheesr,,tfhirestp-orordjeecrtioorniwncoruelmdebnetavlemryeethxopdenss(iBveertto- for i¼1 to N do
o 1 sekas, 2011; Juditsky and Nemirovski, 2012) would be better ½U;S;V(cid:2)¼ðDð∶;∶;iÞÞ;
5 t suited.SeeElyandAeron(2013)foranapplicationofsuchmethods // If constrained ε¼1ρ, If unconstrained ε¼ρλ
1
4/ to the problem of hydraulic fracture monitoring. S¼Shε½S(cid:2);
1
9/ Uð∶;∶;iÞ¼U;Sð∶;∶;iÞ¼S;Vð∶;∶;iÞ¼V;
0
d Algorithm for solving the optimization problems end for
e
ad To solve OPT_TNN(1) and OPT_TNN(2), we use the ADMM // Going back to the original domain
o
nl (Boyd,2011).TheADMMmethod,asadaptedtoourcase,works for i¼3 to p do
w
o bydividingthevariableintwo,withonevariablecapturingtheTNN U←ifftðU;½(cid:2);iÞ;S←ifftðS;½(cid:2);iÞ;V←ifftðV;½(cid:2);iÞ;
D
penalty and the other variable capturing the data-matching con- end for
straints(asanindicatorfunction).Inparticular,theADMMrecur- Z¼U(cid:3)S(cid:3)VT
sion steps (in scaled form) for OPT_TNN(1) are given by
// Dual Update
(cid:2) (cid:3)
ρ B¼BþX−Z
Xkþ1 ¼arg mXin 1Y¼AðXÞþ2kX −ðZk−Bkk2F ; (8) end while
V88 Ely et al.
where I denotes the identity operator. For the unconstrained Thesolutiontoequation9isobtainedbysingularvaluethreshold-
optimization case OPT_TNN(2), the update equation for Xkþ1 ing of the frontal slices in the Fourier domain as shown in Algo-
becomes rithm 2, which shows the overall pseudocode for solving the
problems OPT_TNN(1) and OPT_TNN(2). In Algorithm 2, Shϵ
isanelementwiseshrinkagefunctionthatappliesasoftthresholding
Xkþ1 ¼ðρIþAÞ−1ðYþρðZk−BkÞÞ. (12) according to equation 11 as
g/
r
o
g.
y.se a) Receiver-Y Original
r 1
a
r
b
http://li 1.1
at
e
s
U s1.2
f
o
s
m
er 1.3
T
e
e
s
ht; 1.4
g
yri 5 6 Receiver-X 7 8
p
co b) Receiver-Y Measured
or 1
e
s
n
e
c
G li 1.1
E
S
o
ct t s1.2
e
bj
u
s
n 1.3
o
uti
b
ri
st 1.4
di
Re 5 6 Receiver-X 7 8
3. c) Receiver-Y Reconstruction
1.1. 1
5
8.
1
o
5 t 1.1
1
4/
1
09/ s1.2
d
e
d
a
o
nl 1.3
w
o
D
1.4
5 6 7 8
Receiver-X
Figure3. (a)Thisfigureshowsthefullsyntheticdataforfourdifferentreceiver-sourceslicesaswellas(b)theundersampledmeasureddatafor
thecaseinwhich60%ofthetraceswereremoved.Inaddition,(c)thereconstructionforthe5Dslicesisshownaswell.Allsubfiguresare
shown with the same gain and clip values.
5D seismic tensor completion V89
(
x−ε; if x>ε; thedataarenotsampled)canbeofthesameorlowerrankasthe
Shε½x(cid:2)¼ xþε; if x<ε; (13) originalmatrix.Clearlyinthiscase,thenuclear-normheuristicwill
0; otherwise: fail.Ontheotherhand,suchacasedoesnothappenwithhighprob-
For the constrained noiseless case ε¼ð1∕ρÞ and for the uncon-
strained noisy case ε¼ðλ∕ρÞ where ρ is the step size and λ is 100 Singular value decay
g/ thelow-ranktuningfactorcontrollingtheseverityoftheTNNpen-
or alty. The convergencecriteria in Algorithm 2 are set according to 10−1
eg. relativetolerance.Foralloftheinstancesofthealgorithm,wesetρ
y.s tobeequalto1.Thealgorithmisrobustforthechoiceofthisparam- 10−2
r e
a eter,althoughforsomeverysmallorlargevalues,theperformance d
http://libr of the algorithm degrades. Magnitu1100−−43
at PERFORMANCE EVALUATION: SYNTHETIC 10−5
e DATA
s
U
of Compressibility of seismic data in the tensor singular 10−60 0.5 1 1.5 2 2.5 3 3.5 4 4.5
s value decomposition domain Singular value x 105
m
r
Te Todemonstratethecompressibilityoftheseismicdatawithnon- Figure4. Decayofsingularvaluesofthesyntheticseismicdatathat
e linear reflectors, we generate a synthetic 5D survey in which empirically obey a power law decay.
e
ht; s ssohuortcaensdarnedceirveecreisvpearcsinwgeirnethpelaxc-edanodny-adi1re6c×tio1n6s.gTrhidrewesityhnt5h0e-tmic
yrig Born scatters were placed below the surface. Traces consisting 1 Reconstruction error
op of Nt samples were generated for all the possible source-receiver 0.9 54DD
r c combinations, resulting in a 5D tensor with dimensions 16× 0.8
e o 16×16×16×Nt. Figure 3 shows several receiver gathers from 0.7
ns the 5D synthetic survey for the source located at indices sx¼3
G lice aonfdthseyf¼ou7rwsuitbhpleoatcshcsounbtapilnostsghaothweirnsgffoorutrhreecceoinvsetrangtatXherrse.cEeiavcehr RMSE00..56
E indices(5,6,7,8)andallreceiverYindices.ThetSVDwasapplied N0.4
S
o to th[pe synthetic data. Figure 4 shows the decay of the singular 0.3
ct t values of the matrices X^ð∶;∶;iÞ as computed using the tSVD. 0.2
bje ForthesyntheticdaPtacase,theyobeyanempiricalpowerlawdecay, 0.1
n su swuhreichinimthpeliteSsVthDatdomNi¼a1ink.X^Aðs∶;a∶;reiÞskunltu,cwiseaegxopoedctcoAmlgpolerixthitmym2etao- 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
o Sampling rate
uti recover the data from a limited number of measurements.
b Figure5. Reconstructionerrorasfunctionofthesamplingfraction
ri
st Simple reflectors for the 4D frequency by frequency and full 5D reconstruction.
di Reconstruction frequency by frequency results in nearly identical
e
R Toevaluateouralgorithm,weusethesyntheticdatasetdescribed performance as completion in the original 5D domain.
3. previouslyandremove10%to90%ofthemea-
1.
1. sured traces at 5% intervals, resulting in highly a) 4500
5
8. undersampled data as shown in Figure 3b for
o 1 60% of traces removed. Algorithm 2 was then 4000
5 t applied in two ways: (1) Recovering the data
14/1 bqyuecnocmypalnedtio(2n)orefc4oDvertienngsothrsefernetqiureen5cDysbeyismfreic- m/s) 3500
d 09/ vdoowlunmseamatploinncge.oIpneraadtodri,tiownetoextpeestriinmgenatrawnidthoma b) 9 Y-Sli1c7e Index 24 ocity ( 3000
ade regularlyspacedoperatorinwhicheverysecond Vel
o 2500
nl orthirdshotorreceiverwasremoved.Underthe
w regular downsampling, completion using the
Do TNN proved to be unreliable when only 50% 2000
ofthetraceswereremoved.Thereasonthatregu-
lar downsampling does not work well can be 9 17 24 1500
Y-Slice Index
attributedtothe reasons whyitfails inthe case
ofsimplematrixcompletion.Onereasonbehind
Figure6. Threeverticalslicesofthevelocitymodelusedtogeneratethesyntheticshots.
thisisthatforregulardownsampling,theresult- (a) The simplest velocity model with water and sedimentary layers. (b) Flat layered
ingmatrix(paddedwithzerosatplacesinwhich velocity model with a normal fault.
V90 Ely et al.
abilitywhenusingrandomdownsampling.Werefertheinterested wouldneedtobeselectedforeachfrequency.Thisisduetothefact
reader to Recht et al. (2011) for a more technical discussion, thatthesourcewaveletisbandlimited,resultinginadifferentS/N
implicit in the analysis there. An analysis of such conditions for across frequencies, which in turn necessitates optimizing overthe
our case is a subject of current investigations. The normalized regularization parameter for each frequency.
root-mean-square error (NRMSE) measures the recovery error.
Forthesyntheticdata,Figure5showstheperformanceasfunction Finite-difference synthetic
org/ oufndsearmsapmlipnlgedratdea.tAa,lththoeughreftlheectoNrRsMaSreEsistilsligvniisfiibclaentinforthheighrley- BecauseBornscattersinaccuratelyrepresenttypicaldatasets,we
eg. construction; refer to Figure 3. Therefore, the NRMSE appears generateamorerealisticsetofsyntheticsurveyovertwovelocity
s models of varying complexity: (1) a flat layered model and (2) a
y. to underestimate reconstruction quality. For all sampling rates,
r layered model with a normal fault. Figure 6 shows cross sections
a thecompletionperformanceinthe5Dor4Ddomainisnearlyiden-
br for the two 2000-m cubed velocity models. For each of the two
p://li tFicigaul,rew5it.h reconstruction error varying less than 1% as shown in mgeoodmelest,ritewsowoitrhthtohgeosnaamlesutrovtaelysnuwmebreergoenfesrhaotetsdawndithredciefifveerersn.tSsyhno-t
at htt FoAurlitehroudgohmaminantoy erexcpolonisttrtuhcetisopnartseictyhnoiqfutehsemdautsatcoupbeeraitnetihnetfh-xe theticshotswerefiredovera16×16gridanda32×32-receiver
Use domain,theperformanceofourmethoddoesnotappeartodepend gprlaidc.eFdoorvtheresam1a0l0l0sp×ati1al0s0u0r-vmeyagreeaomceentrtye,reshdoitnsatnhdermecideidvleersowf tehree
f onthetemporal Fouriertransform. Thisis dueto thefact thatthe velocity model, resulting in shot spacing of 62.5 m and receiver
o
s constructionofthetSVDinvolvestherecursiveapplicationofthe spacing of 31.25 m. For the large spatial geometry, the shots
m
r Fouriertransformalongthetemporalaswellasthespatialdimen- andreceiverswerespreadacrosstheentire2000×2000-msurface
e
T sions of tensor. After taking the Fourier transform along the tem- withshotspacingof125mandreceiverspacingof62.5m.Wethen
ee poral dimensions, our method implicitly operates in the spatial used PySit, a seismic inversion and forward-modeling toolkit, to
s
ht; Fourier domain. Given that time and Fourier domain reconstruc- generate the synthetic traces using a finite-difference solver such
g tionsofthetensorarenearlyidentical,itcouldbecomputationally that the resulting data sets had a spatial dimension of 16×16×
yri advantageous to reconstruct the data in a frequency-by-frequency 32×32 in the shot-receiver domain. For each of the twovelocity
p
co mannerbecausethisistrivialtoparallelizeacrossnumerousnodes. modelsandtworeceivergeometries,weapplyAlgorithm1tocal-
or Fora5DdatasetwithNtsamples,thecompletionwouldonlyneed culate the tSVD decomposition of the four data sets. Using these
se to be performed on Nt∕2 4D tensors due to the symmetry of the decompositions,wegeneratetruncatedrank(compressed)versions
en Fourier transform. The number of computations could further be of the data sets and compare the error between the low-rank
c
EG li rNetd∕u2cefrdeqifuetnhceiessigisnnaleeidsebda.nGdivleinmtihteadtthanedtSoVnDlysharifnrkacatgieonopoefratthoer roefcsoinnsgturulacrtivoenctaonrdpoairrisgiunsaeldd.aFtaigausrea7fushnocwtiosnthoefnthoermtoatlaizlendummebaenr
S alreadyactsinthefrequencydomainandtheSVDisalreadycom- squarederror(NMSE)ofthereconstructionsversusthefractionof
o
ect t pduotmedainoninacucrosmnpoleaxdmdiatitorinxa,ltrcaonmsfpourtmatiinognatlhecodsattsa. intotheFourier sairnrigvualla.rFvroecmtorthpeaseirscuforvretsh,ewweavseeefietlhdatartrhiveinsgurvafetyergtehoemdeitrreiecst
ubj Despite these advantages of the Fourier domain, we choose to andcomplexityofthevelocitysignificantlyaffectthedata’scom-
n s workinthetimedomainformostofthesyntheticandrealdatasets pressibility,withsurveygeometrybeingmostdramatic.Withawide
butio tpoarsaimmpetleifrys.thIneappaprrtiocauclhara,nidf ttohererdeuccoensthtreucntuiomnbearlgoofrriethgmulawriazsatiaopn- smhaottioanndtharenctehievesrmsaplalcsipnagti,aelagcehomtreatcrey,cionnwtahinicshmeaocrheturanciqeuceonintafionrs-
stri pliedtonoisydatainthef-xdomain,thenaregularizationparameter more redundant information. For any survey geometry, the com-
di plexity of the model inversely correlates with
e
R compressibility. Given that the more complex
1.1.3. a) 0.2 Reconstruction error b)100 CompressionF acuultr v(semall) vwealvoecfiiteyldmgoedoemlewtroyu,lwdeincvoonlcvleudaemthoartemcoormepslienx-
18.5 00..1168 10−1 LFaayueltr (eldar (gsem)all) gwualvaerfiveelcdt.orFsrowmoutlhdebdeenmeoendsetdratteodrerperlaetsieonntshthipe
o Layered (large)
5 t 0.14 among compressibility, survey geometry, and
14/1 MSE00.1.12 MSE10−2 tfhoerrveeclooncsittyrumctoiodnelo,fwmeisesxipnegcttraacseism.Pilearrfboermhaavnicoer
09/ N0.08 N10−3 could also be improved by applying the recon-
d struction on a small spatiotemporal window to
e 0.06
oad 0.04 10−4 reduce the complexity and nonlinearity of the
wnl 0.02 traFceosr.thefoursurveys,between25%and97.5%
Do 00 20 40 60 80 10−6 10−4 10−2 100 of the source-receiver pairs are randomly re-
Sampling rate % Fraction of SV used moved,andweusethefullnoiseless5Dcomple-
tionalgorithmtoreconstructthemissingtraces.
Figure7. (a)Reconstructionerrorforthefourdatasetsasafunctionofthedownsam- WeapplyAlgorithm2tothedecimatedtracesfor
plingrate.(b)CompressibilityofthedatasetsinthetSVDdomainasafunctionofthe
500iterationsfromthefoursurveysanddifferent
fractionofsingularvectorpairsusedtoreconstructthedataset.Notethatifthedataare
morecompressibleinthetSVDdomain,thenthecorrespondingrecoveryerrorisalso samplingrates.Figure7showstheNMSEforthe
small. wavefield after the direct arrival. Regardless of
5D seismic tensor completion V91
a) Original
7 13 20 Receiver-Y 26
0.5
1
1.5
g/ s
r 2
o
g.
e
s 2.5
y.
r
a
br 3
p://li 0 1R0eceive2r0-X 30 0 10 20 30 0 10 20 30 0 10 20 30
at htt b)0.5 7 13 20% 20 Receiver-Y 26
e
s
U
f 1
o
s
m
1.5
r
e
T s
e 2
e
s
ht; 2.5
g
ri
y
p 3
o 0 10 20 30 0 10 20 30 0 10 20 30 0 10 20 30
c
or Receiver-X
e 10%
ns c) 7 13 20 Receiver-Y 26
e 0.5 5
c
G li
E 1 1
S
o
ct t 1.5 5
bje s 2 2
u
s
n
o 2.5 5
uti
b
ri 3
st 0 10 20 30 0 10 20 30 0 10 20 30 0 10 20 30
di
e Receiver-X
R
5%
3. d) 7 13 20 Receiver-Y 26
1. 0.5
1.
5
8. 1
1
o
5 t 1.5
1
4/ s
1 2
9/
0
ed 2.5
d
a
o
wnl 3 0 10 20 30 0 10 20 30 0 10 20 30 0 10 20 30
Do Receiver-X
Figure8. Thetopsetoftracesinthisfigureshowsfourgathersforashotatthecenterofthesurveyat(x,y)sourceindicesof(8,8).Each
panel shows a constant receiver-y-index. The bottom three panels show the reconstructed traces at three different sampling rates,
20%, 10%, and 5%.