Wave-equation reflection tomography: annihilators and sensitivity kernels

21 Pages·2006·0.86 MB·English

Geophys.J.Int. (2006)167,1332–1352 doi:10.1111/j.1365-246X.2006.03132.x
Wave-equation reflection tomography: annihilators and sensitivity
kernels
Maarten V. de Hoop,1 Robert D. van der Hilst2 and Peng Shen3
1CenterforComputationalandAppliedMathematics,PurdueUniversity,150N.UniversityStreet,WestLafayetteIN47907,USA
2DepartmentofEarth,AtmosphericandPlanetarySciences,MassachusettsInstituteofTechnology,Rm54-517A,Cambridge,MA02139,USA.
E-mail:[email protected]
3TotalE&PUSA,800Gessner,Suite700,HoustonTX7024,USA
Accepted2006July6.Received2006July2;inoriginalform2005December13
SUMMARY
Inseismictomography,thefinitefrequencycontentofbroad-banddataleadstointerference
effectsintheprocessofmediumreconstruction,whichareignoredintraditionalraytheoretical
implementations. Various ways of looking at these effects in the framework of transmission
tomographycanbefoundintheliterature.Here,weconsiderinversescatteringofbodywavesto
developamethodofwave-equationreflectiontomographywithbroad-bandwaveformdata—
whichinexplorationseismicsisidentifiedasamethodofwave-equationmigrationvelocity
analysis.Inthetransitionfromtransmissiontoreflectiontomographytheusualcrosscorrelation
betweenmodelledandobservedwaveformsofaparticularphasearrivalisreplacedbytheaction
ofoperators(annihilators)totheobservedbroad-bandwavefields.Usingthegeneralizedscreen
expansionforone-waywavepropagation,wedeveloptheFre´chet(orsensitivity)kernel,and
showhowitcanbeevaluatedwithanadjointstatemethod.Wecastthereflectiontomography
intoanoptimizationprocedure;thekernelappearsinthegradientofthisprocedure.Weinclude
y
g anumericalexampleofevaluatingthekernelinamodifiedMarmousimodel,whichillustrates
o
ol thecomplexdependencyofthekernelonfrequencybandand,hence,scale.Inheterogeneous
m mediathekernelsreflectproperwavedynamicsanddonotrevealaself-similardependence
s
ei on frequency: low-frequency wave components sample preferentially the smoother parts of
S
themodel,whereasthehigh-frequencydataare—asexpected—moresensitivetothestronger
I
J heterogeneity.Wedeveloptheconceptforacousticwavesbuttherearenoinherentlimitations
G
fortheextensiontothefullyelasticcase.
Keywords: migrationvelocityanalysis,reflectiontomography,sensitivitykernels.
1 INTRODUCTION
Theevergrowingvolumesofdenselysampledbroad-banddata,bothinindustryapplicationsandinglobal(array)seismology,areproviding
new opportunities and challenges for the development of techniques for subsurface imaging that exploit efficiently the rich information
containedinseismicwaveforms.Here,weanalyseamethodofwave-equationreflectiontomographywithfinite-frequencydata(derivedfrom
DeHoop&VanderHilst2003);inexplorationseismicsonereferstosuchmethodsingeneralaswave-equationmigrationvelocityanalysis
(MVA,Biondi&Sava1999).Itmakesuseoftheredundancyinthescatteredwavefielddata;thescatteringisassumedtooccursomewhere
inthesubsurface(Fig.1,bottom),butthelocationsofthescatterersneednotbeknown.
Thegeneralobjectiveoflinearizedseismictomographyistofindbyinversionamodel,oraclassofmodels,ofelasticpropertiesthat
throughforwardmodellingpredictsthedata.Thedatafitisusuallyevaluatedthroughaparticularcriterion,oftenexpressedasacostorpenalty
functionalthatistobeoptimized.Intheprocessofoptimization,thekeyquantitytobefoundisthesocalledFre´chetorsensitivitykernel.
Thecriterion—and,therefore,thecharacterofthesensitivitykernel—dependsonthetypeofobservationbeinginterpreted.
Intraveltimetomographytheobserveddataarethemeasuredarrivaltimesofknownseismicphases,andtheinversiontypicallyinvolves
a backprojection algorithm for estimating medium wave speeds relative to some reference model. In the idealized case of infinitely high
frequencywavepropagationthetraveltimeandkernel(raypaths)arewelldefined.Withfinitefrequencydata,however,thedirectdetection
ofthetraveltimeisoftenobscuredbywaveinterferenceanddiffractionsothatalternativewaysofcomparingmodelledandobserveddata
mustbesought.Intheframeworkoftransmittedbodywavesvariouswaysofaccountingforfinitefrequencyeffectsintheprocessofmedium
estimationhavebeenproposedandanalysed(Luo&Schuster1991;Woodward1992;Dahlenetal.2000;DeHoop&VanderHilst2005).
Insuchapproachestofinite-frequencytomographyonecandistinguishbetweenmethodsbasedonexplicitbackprojectionoverFresnel-like
1332 (cid:2)C 2006TheAuthors
Journalcompilation(cid:2)C 2006RAS
Onwave-equationreflectiontomography 1333
S S
410
410
SS
660
S S
660
CMB
ScS
SS
S S
410
S S
660
S s
ScS 410
S s
660
S s
660
S s
410
Figure1. Top:Examplesofreflectedwavefieldsthatcanbeconsideredinthewaveequationreflectiontomographydevelopedinthispaper.Themedium
reconstruction(thatis,wavespeedestimation)takesplacewhereupanddowngoingwavefieldsinteract:thehigherthedensityofscatterers(e.g.reflectors),the
moreconstrainedthesolution.Bottom:Cartoonofsourceandreceiverarrayconceptexploitedbythereflectiontomographyapproach.
volumes(insteadofbackprojectionalonginfinitesimallynarrowrays),see,forinstance,Dahlenetal.(2000)andDeHoop&VanderHilst
(2005),andthosecombiningthewaveequationwiththeBornapproximationusinganadjointstatemethod(Vascoetal.1995;Trompetal.
2005).Inwave-equationtomographythewavefieldsinthekernelcanbecomputeddirectlyfromthetime-domainwaveequation,usingGreen’s
functionscalculated,forinstance,byspectralelements,normalmodesummation(Zhao&Jordan1998),orfrequency-domainone-waywave
propagation.
Herewefocusonamethodofwaveequationreflectiontomography,nottodetectandinvestigatereflectorsbuttoestimatesubsurface
wavespeedvariations.Weuseafrequencydomainformulationandassumesinglescattering.Incontrasttotomographicapproachesthat
interpretspecificphasearrivals—measured,forinstance,byphasepickingorwaveformcrosscorrelation—weconsidertheentire,broad-band
wavefield(althoughsometimewindowingorotherprocessingmayneedtobeappliedtomutethepartsofthedatainfluencedbymultiple
scattering).Insteadofmeasuringsingletime-shifts,itinvolvestheapplicationofannihilatorstotheobserveddata.Theseannihilatorsare
operatorswhoseactiononthewavefieldvanishesifthewavespeedmodelpredictstheobservationsused:thelevelofannihilationisthus
a measure for the accuracy of estimated wave speed variations. The annihilation criterion has a significant advantage over conventional
traveltimemismatchcriteriainthatitexplicitlyaccountsfortheeffectsofscatteringinEarth’sinterior,whileitexploitstheredundancyin
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS
1334 M.V.deHoop,R.D.vanderHilstandP.Shen
the(scattered)wavefield.Suchredundancycanarisefromobservationsatmultipleoffsets(inexplorationseismology)oratmultipleangular
epicentraldistancesandazimuths(inglobalseismology).Forexample,inthestudiesofsedimentarybasinsonecanusewavesscatteredfrom
themanysubsurfacereflectorsandfaultstoelucidatetheinterveningseismicvelocities(e.g.Zeltetal.2003),andonaglobalscaleonecould
usethewavefieldsformedbyreflectionoffthecoremantleboundaryortheuppermantlediscontinuities(seeFig.1,top)tostudymantle
heterogeneity(insteadofmantlelayering,whichhasbeendonebeforeSipkin&Jordan1976;Revenaugh&Jordan1991).
Themethodpresentedhereinvolvesseveralsteps,whichalsodefinethestructureofthispaper.First,wedecomposethebroad-band
wavefieldintoupanddowngoingwaves(Section2).Thisdirectionaldecomposition(e.g.Fishman&McCoy1984a,b)leadstoacoupled
systemofone-waywaveequations.Second,weusethissystemtoobtainaso-calleddouble-square-root(DSR)equationforthedownward
continuationofthedatatoselecteddepthsinthesubsurface(Clayton1978;Claerbout1985),simultaneouslyforthesourceandthereceiver
side(Section3).Third,afterdownwardcontinuation,thedecomposedwavefieldsaresubjectedtoa(wave-equation)angletransform,which
enablestheexploitationofredundancyassociatedwithdifferentscatteringanglesandazimuthsatthesedepths(Section4.1).Thistransform,
whichisrelatedtobeamforming(Scherbaumetal.1997),isusedtogeneratecommon-image-pointgathers(DeBruinetal.1990);inSava
etal.(1999)andSava&Fomel(2003)adifferentbutrelatedtransformimplementedasapost-processingoperatorisused.Thesegatherscan
bere-arrangedintomultipleimagesofthesamepartofEarth’sreflectingstructurefromwavesscatteredoverarangeofangles(Stolk&De
Hoop2006).Fourth,inSection4.2wediscusshowtheannihilationofthedownwardcontinueddata(again,atspecificdepthsinthemedium)
formsthecriterionforselectingacceptablewavespeedmodels.InSection5wederivethesensitivitykernelsunderlyingtheoptimizationof
thiswavefieldannihilation.
InSection6weillustratethemultiscaleaspectsofthefinite-frequencysensitivitykernel,thecharacterofthecostfunctionalusedfor
optimization,andtheeffectivenessoftheangletransform.Weshowthatinheterogeneousmediathedependenceofthekernelsonfrequency
iscomplexandnotself-similar.Low-frequencydatasampleheterogeneousmodelsinwaysthatareverydifferentfromhigh-frequencydata.
Itisthisfinitefrequencybehaviourthatallowsonetorelatescalesinthedata(frequencies)toscalesinthemodel(spatialscales)andinthe
physicalprocessesthatcausethem.
Theimagegathersproducedbytheangletransformarewithoutartefacts(thatis,falseevents)ifthemediumwavespeedsarecorrect,
eveninthepresenceofcaustics.(ThiswasnotthecaseintheapproachfollowedbyBrandsberg-Dahletal.2003,whichistheprecursorof
thedevelopmentpresentedhere).Forincorrectwavespeedmodels,however,theimagegatherswillrevealdependencyonscatteringangle
andazimuth;thatis,therewillbearesidualmoveout.Thismoveoutistheunderlyingdiagnosticexploitedhereandisevaluatedusingdata
annihilators(zeromoveoutanduniformamplitudeswillproduceperfectannihilation).Theannihilationoperatorsdependonthewavespeed
modelanddeterminehowwellobservationsarematchedbymodelpredictions(Stolk&DeHoop2002);awavespeedmodelisacceptableif
thedataareintherangeof,thatis,whentheycanbepredictedbyourmodellingoperator.1Thetomographicinversioncanthenbeformulated
as the problem of minimizing the action of these annihilators on the data. An important and unique property of annihilators is that they
smoothlydependonthewavespeed(Stolk&Symes2003).
WedevelopthetheoryuptotheexplicitexpressionfortheFre´chetkernelunderlyingtheannihilatorapproach(Section5).Theevaluation
ofthiskernelisformulatedasanadjointstatemethod.Suchmethodshavebeenexploitedinseismicimagingandinversescattering(see,
forexample,Tarantola1987;Stolketal.2006)andhavebeguntobeusedinwaveequationtransmissiontomography(Trompetal.2005).
Oncethekernelisavailable,standardoptimizationmethods(suchasconjugategradients)canbeinvokedtoselectaparticularwavespeed
model.Weillustratethemethodanditseffectivenesswithasyntheticdataexample.Forsimplicityofformulationweconsiderhereacoustic
(compressional) waves, possibly in transversely isotropic media with vertical symmetry axis, described by a scalar wave equation (see,
e.g.Alkhalifah1998;Schoenberg&DeHoop2000),andweusetheflattenedEarthassumption.Therearenoobstructions,however,for
extendingtheformalismtothefullelasticcase.
Theannihilator-basedapproachtowave-equationreflectiontomography(ortoMVA)isakintothenotionofdifferentialsemblancein
explorationseismics.Systematicmethodsforupdatingwavespeedmodelsbyoptimizationusingthedifferentialsemblancecriterionhave
beenintroducedbySymes&Carazzone(1991),andfurtherdevelopedandappliedbyChauris&Noble(2001),Mulder&TenKroode(2002),
andothers.Wave-equationMVAbasedonfocussinginsubsurfaceoffset(relatedtoannihilatorsofimagegathers,asexplainedinStolk&
DeHoop2001,2006)andshortrecordmigrationratherthandownwardcontinuationmigrationwasmorerecentlydevelopedbyoneofthe
authors(Peng2004).Theworkpresentedherebuildsonvariouspreviousdevelopments:ThegeneralizedBremmercouplingseriestomodel
seismicreflectiondata(DeHoop1996),inversescatteringbaseduponthefirst-ordertermofthisseries,thewave-equationangletransform
andannihilators(Stolk&DeHoop2001,2006),therelationbetweenwavefieldreciprocityandoptimization(DeHoop&DeHoop2000),
andthegeneralizedscreenexpansionforone-waywavepropagation(DeHoopetal.2000;LeRousseau&DeHoop2001).
2 DIRECTIONAL WAVEFIELD DECOMPOSITION:
SINGLE SQUARE-ROOT OPERATOR
In our method for wave-equation reflection tomography the measure for selecting acceptable models is the successful annihilation of the
observedwavefield.Theannihilatingactionis,here,carriedout—atselecteddepthsinthemodel—onimagegathersthatareconstructedfrom
1TheprincipleofcharacterizingtherangeofmodellingoperatorscanbefoundinGuillemin&Uhlmann(1981)andhasbeendevelopedforseismicbody-wave
scatteringbyDe Hoop&Uhlmann(2005).
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS
Onwave-equationreflectiontomography 1335
downwardcontinueddata(Section4).Thedownwardcontinuationoftheobservedwavefieldinvolvesdecompositionintoupanddowngoing
constituents.Forthispurposeweintroduceheretheone-waywaveequationsandtheirassociatedwavefieldpropagators.Initially,wecarry
outthedirectionaldecompositionforwavepropagationinasmoothbackgroundmodel(forthegeneralformulation,seeFishman&McCoy
1984a,b);inlatersectionsweaccountforscatteringofwavesoffreflectorssuperimposedonthisbackground.
2.1 Thewaveequation,evolutionindepth
Wemakedepthz—insteadoftimet—theevolutionparameterforwavepropagationanddatacontinuation.Thehorizontalcoordinatesare
collectedinx.ThepartialderivativesaredenotedbyD =−i∂ , D =−i∂ , D =−i∂ sothattheirFourierdomaincounterparts
x1,...,n−1 x1,...,n−1 z z t t
becomemultiplicationsbyξ1,...,n−1(thenegativeofhorizontalwavevector),ζ (thenegativeofverticalwavevector)andω(frequency).Here,
n=2or3for2-Dor3-Dseismics,respectively.Analternativenotationforthewavevectorusedintheseismicliteratureis(k ,k ,k )for
x y z
(−ξ ,−ξ ,−ζ)(n=3).
1 2 (cid:2)
Tobringouttheuseofzastheevolutionparameter,thescalarwaveequationforcompressionalwaves,(cid:6)u−c (x,z)−2u¨ ≡( n−1∂2 +
∂2−c (x,z)−2∂2)u= f,iswrittenasthefirst-ordersystem 0 i=1 xi
z (cid:3) 0(cid:4) (cid:3) t (cid:4)(cid:5) (cid:6) (cid:3) (cid:4)
∂ u 0 1 u 0
= + , (1)
∂z ∂u −L(x,z,D ,D)0 ∂u f
∂z x t ∂z
(cid:2)
whereuisascalarwavefieldquantitysuchaspressure,thepartialdifferentialoperator L(x,z,D ,D) = c (x,z)−2D2− n−1D2 (see,
x t 0 t i=1 xi
e.g.Aki&Richards2002,p.25),c isthewavespeedinthe(background)medium,andf isasourceterm.(N.B.herewetakeforf avolume
0
injectionsource,butageneralbodyforcecanbetreatedinasimilarmanner.)Withexp[i(ξx +ωt)]representingaFourierfactorinaplane
wave,L(x,z,ξ,ω)isdefinedas
L(x,z,ξ,ω)=exp[−i(ξx+ωt)]L(x,z,D ,D)exp[i(ξx+ωt)]=c (x,z)−2ω2−(cid:4)ξ(cid:4)2. (2)
x t 0
Forspacedimensionn>2,ξxisavectordotproduct.Wenotethat,forgivenz,(x,t,ξ,ω)arecoordinatesonphasespace,andeq.(2)defines
afunctiononthisspace.Forwavepropagationinheterogeneousmedia,carryingouttheanalysisinphasespacehasimportantadvantages
overanalysisinmoreconventional(n +1)-Dspace–time.Heuristically,itallowsadescriptionofwavesintermsoflocallyplanewaves.
Theexplicitdefinitionandconstructionoftheone-waywaveequationsbelow,theircorrespondingfundamentalsolutions(orpropagators)
in the form of path integrals, and the marching computational algorithms are all rooted in phase space (Fishman et al. 1987; Fishman
2004).
2.2 Thesystemofone-waywaveequations
Inphasespace,therootsofL(x,z,ξ,ω)correspond,locally,totransitionsfrompropagating(L>0)toevanescent(L<0)waveconstituents.
Awayfromtheseroots,inparticularinthepropagatingregime,system(1)canbetransformedintodiagonalform(e.g.Stolk&DeHoop
2005).WithDandUdenotingthedown-andupwardpropagatingconstituentsofthewavefield,respectively,operatormatricesQ(z)=Q(x,
z,D ,D )canbeconstructedsuchthat
(cid:3) (cid:4)x t (cid:3) (cid:4) (cid:3) (cid:4) (cid:3) (cid:4)
u u f 0
D =Q(z) , D =Q(z) ,
u ∂u f f
U ∂z U
effectivelysatisfytheone-waywaveequations
(cid:3) (cid:4) (cid:3) (cid:4)
∂ ∂
+iB (x,z,D ,D) u = f , −iB (x,z,D ,D) u = f . (3)
∂z D x t D D ∂z U x t U U
AnycouplingbetweenDandUconstituentsisoftheformQ(z)∂ Q(z)−1.
∂z
Inthepropagatingwaveregime,theoperatorB(=B orB )admitsanintegralrepresentationofthetype
(cid:7) (cid:7) (cid:7) (cid:7) U D
(Bu)(x,t,z)=(2π)−n b(x(cid:5),z,ξ,ω)exp[iξ(x−x(cid:5))]exp[iω(t−t(cid:5))]u(x(cid:5),t(cid:5),z)dx(cid:5)dt(cid:5)dξdω. (4)
Here b(x(cid:5), z, ξ, ω) is a smooth function on phase space2 and is, hence, better suited for manipulations (such a(cid:8)s taking the wave speed
derivativesneededbelow)thantheoperatorBitself.Forsystem(1),forhighfrequencies,3 b = b(x,z,ξ,ω) = ω 1 −ω−2(cid:4)ξ(cid:4)2,and
(cid:9) (cid:9) c0(x,z)2
forω>0,bcoincideswith L(x,z,ξ,ω)whileforω<0,bcoincideswith− L(x,z,ξ,ω).Forthisreason,Bisoftenreferredtoasthe
‘single-square-rootoperator’.
2Inthemathematicsliteratureb(x(cid:5),z,ξ,ω)wouldbereferredtoasthe(antistandard)symbolofoperatorB.Likewise,L(x,z,ξ,ω)ineq.(2)isthesymbolof
L(x,z,Dx,Dt).Wedeviatefromthecalculusofstandardsymbolstofacilitatethedevelopmentofourpreferredrepresentationforthesensitivitykernel.
3ForhighfrequenciesthesymbolofBreducestowhatiscalledtheprincipalsymbol.Forsimplicity,wewilldenotetheprincipalsymbolofBbybfromnow
on.Weemphasize,however,thattheprincipalsymbolexpressionforbisnotexact.ThefrequencydependenceofoperatorBwasdiscussedbyFishman&
McCoy(1984a).
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS
1336 M.V.deHoop,R.D.vanderHilstandP.Shen
Withbrepresentingtheverticalwavenumber,ω−1bhastheappearanceofaverticalwaveslownessatthepoint(x,z)whereasω−1ξ isa
horizontalslowness.For(x,z,ξ,ω)suchthatbisreal,theone-waywaveequationsareofhyperbolictype,describingpropagatingwaves.4
(Weremarkthateventhoughtheseparationbetweenlocallypropagatingandevanescentwavesis,here,basedupontheanalysisofprincipal
symbols,thebasicpictureoflocallypropagatingandevanescentphasespaceregimesremainsintactforamorecompletedescriptionofthe
symbols.)
Inthefrequencydomain,weintroducetheoperatornotationBˆ inaccordancewithBu= F−1 BˆF uwhereFdenotestheFouriertransform.5
ω→t t→ω
WenotethatBˆ = Bˆ(x,z,D ,ω).
x
WechooseanormalizationforoperatorQ(z)suchthateq.(3)isself-adjoint;thus
(cid:5) (cid:6)
(Q∗)−1 −HQ
Q= 1 D D , (5)
2 (Q∗)−1 HQ
U U
where ∗ denotes the adjoint, H denotes the Hilbert transform in time, and QD,U = QD,U(z) = QD,U(x, z, Dx, Dt) are operators with
pverirnticciapla-lacsoyumsbtiocl-spo(wc0(eωxr2,zfl)2u−x((cid:4)Dξe(cid:4)2H)−o1o/p4.1T9h9e6p).hWysiitchalthmiseannoirnmgaolifztahtiiosnc,htohiececooufpQlinDg,UQis(zth)a∂tQth(ez)d−o1wbnetawnedenupDgoainndgUfiecldosnsatrietuneonrtmsbaleiczoemdeins
∂z
asymptoticallyoflowerorderand,therefore,isneglectedhere;weget
f =−1HQ f
u= Q∗u +Q∗u , D 2 D
D D U U fU= 12HQUf, (6)
withQ∗u representingtheup-andQ∗u representingthedowngoingconstituentofwavefieldu.
U U D D
2.3 One-waywavepropagators
Ifsingularitiespropagatenowherehorizontallybetweenzandz ,withz<z ,thereisawell-definedsolutionoperatorG (z,z )oftheinitial
0 0 U 0
valueproblemforu givenbyeq.(3)with f =0(Stolk2004);thisoperatordescribespropagationinthe(upward)directionfromz toz
U U 0
(Stolk&DeHoop2005,Section2).(Analternativedescriptionofone-waypropagators,inwhichistheissueofnowherehorizontalraysis
suppressed,canbefoundinFishmanetal.1997;DeHoop&Gautesen2000).IntheIntermezzoinSection(4.1)webrieflyindicatehowthe
full-wavesolutioncanbeusedaswell.Asolutionfortheinhomogeneouseq.(3)isthengivenby
(cid:7) (cid:7) (cid:7)
∞
u (x,t,z)= (G (z,z ))(x,t−t ,x )f (x ,t ,z )dx dt dz , (7)
U U 0 0 0 U 0 0 0 0 0 0
z
or,inoperatorformby
(cid:7)
∞
u (.,z)= G (z,z )f (.,z )dz . (8)
U U 0 U 0 0
z
Here,theGreen’sfunction(G (z,z )),relatingthewavefieldat(x,t)anddepthztotheforceat(x ,t )atdepthz ,isthekernelofone-way
U 0 0 0 0
propagatorG (z,z ).6 Theadjoint G (z,z )∗ describesdownwardpropagationfromztoz ofeq.(3)orupwardfromz tozinreversed
U 0 U 0 0 0
time.7
2.4 Modelrepresentationandperturbation
Fortheoptimizationunderlyingourreflectiontomographyweneedthe(Fre´chet)derivativesoftheabovementionedsquare-rootoperators
withrespecttothebackgroundwavespeedc (x,z).Suchderivatives,whicharederivedinAppendixA,aretypicallyexpressedintermsof
0
aperturbationδc andcanbeinterpretedasfollows.Weconsiderafinite-dimensionalsubspaceof(smooth)wavespeedmodelsc ,spanned
0 0
byafinitesetofbasisfunctions{φ }:
k
(cid:10)M
E :(c )→c (x,z)= c φ (x,z)= E (c )(x,z).
c k 0 k k c k
k=1
Thissubspacecouldbegenerated,forexample,byexpandingthewavespeedintocubicB-splines.Thederivativesineq.(A4)canthenbe
(cid:2)
expressedintermsofderivativeswithrespecttothecoordinates{c }throughδc (x,z) = M δc φ (x,z).Theadjoint E∗ of E ,which
k 0 k=1 k k c c
satisfies
(cid:11) (cid:12)
(cid:9)Ec(ck),c0(cid:10)(x,z)=(cid:9)(ck), Ec∗c0 (cid:10)IRM,
4Intheregimewherebisreal,theprincipalsymbolsoftheone-waywaveoperators—ζ±b(x,z,ξ,ω)uptoafactori—canbeidentifiedasHamiltoniansfor
downandupgoingwaves,andcanbeusedtotracerayswithzasevolutionparameter.Thenextordersymbolneedstobeaccountedfortoobtainthetransport
equationleadingtothecorrectasymptoticsolutions.
(cid:13) (cid:13)
5 Ft(cid:5)→ωu= exp(−iωt(cid:5))u(t(cid:5))dt(cid:5), Fω−→1tuˆ =(2π)−1 exp(iωt)uˆ(ω)dω, andsimilarlyforFx−(cid:5)→1ξand Fξ−→1x indimensionn−1.
6Theupgoingconstituent‘U’ofthefull-waveGreen’sfunctionisgeneratedby 12HQ∗U(z)GU(z,z0)QU(z0)(cf.6).
7Wenotethat,eventhoughtheprincipalsymbolapproximationisahigh-frequencyapproximation,itcanstillbeusedtoobtainanapproximationtothe
propagatorrevealingawavebehaviour.
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS
Onwave-equationreflectiontomography 1337
istheprojection(interpolation)
(cid:3)(cid:7) (cid:7) (cid:4)
(cid:11) (cid:12)
c → c (x,z)φ (x,z)dxdz = E∗c .
0 0 k c 0 k
Thisadjointappearsintheapplicationoftheadjointstatemethodtoevaluatingthesensitivitykernel.
3 DATA CONTINUATION IN DEPTH: DOUBLE-SQUARE-ROOT OPERATOR
Thebackgroundmediumc (withsmoothperturbationsδc )consideredintheprevioussectioncanproducecausticsbutnoscattering(such
0 0
asreflections).Wenowdescribehow,withtheBornapproximation,single-scatteredphasesaremodelledintheframeworkoftheone-way
wavetheorydevelopedabove.Makinguseofreciprocity,wemodelthedatapurelybyupcomingpropagators.Whereasupanddowngoing
wavesaremodelledwithasinglesquare-rootoperatorB(seeSection2.2),theinteractionofthesecontributions(throughatimeconvolution)
yieldsaDSRoperator.
3.1 Upwardcontinuationandmodellingofreflectiondata
Weassumewavefieldscatteringoffacontrastδcthatcontainssingularvariationsinmediumwavespeed(forinstance,reflectors).Thetotal
mediumwavespeedfollowsthedecomposition(withδc−2(cid:11)−2c−3δc)
(cid:14) (cid:15) 0
c−2(x,z)=γ−2(z)+ c−2(x,z)−γ−2(z) −2c−3(x,z)δc(x,z).
0 0 0 0
Toobtainadownward/upwardcontinuationrepresentationofthesinglescatteredwaves,weintroducetheextendedmediumcontrast(Stolk
&DeHoop2005),
(cid:3) (cid:4)
1 (cid:11) (cid:12) r(cid:5)+s(cid:5)
R(s(cid:5),r(cid:5),t(cid:5),z(cid:5))= δ(t(cid:5))δ(r(cid:5)−s(cid:5)) c−3δc ,z(cid:5) . (9)
2 0 2
OnecanviewRasadata-likerepresentationofthecontrast.Usingthefootnotebeloweq.(8)concerningtherelationbetweenone-wayGreen’s
functionsandtheGreen’sfunctionoftheoriginalwaveequationandusingthatH2=−I,thesingularpartofthedata(asrecordedatEarth’s
surface,thatis,atz=0)canbemodelledas
(cid:7)
∞
d(s,r,t)= Q∗U,s(0)Q∗U,r(0) (H(0,z(cid:5))QU,s(z(cid:5))QU,r(z(cid:5))Dt2R(.,z(cid:5)))(s,r,t)dz(cid:5), (10)
0
whereQU,s(z)isshortforQU(cid:7)(s,z,Ds,Dt)andQU,r(z)isshortforQU(r,z,Dr,Dt),andthekernelofH(z,z(cid:5))isgivenbythetimeconvolution
(H(z,z(cid:5)))(s,r,t,s(cid:5),r(cid:5),t(cid:5))= (G (z,z(cid:5)))(s,t−t(cid:5)−t¯,s(cid:5),0)(G (z,z(cid:5)))(r,t¯,r(cid:5),0)dt¯, forz>z(cid:5). (11)
U U
IR
ThroughH asineq.(11),werecognizeineq.(10)thetimeconvolution(integrationovert¯)ofanincomingGreen’sfunctionthatconnects
thesourceats(atzerodepth)withthescatteringpointatr(cid:5)=s(cid:5)(atdepthz(cid:5))withtheoutgoingGreen’sfunctionthatconnectsthisscattering
pointwiththereceiveratr(atzerodepth).ThisprocessgivesrisetoaDSRpropagator;thespatialnotationisillustratedinFig.2(seealso
Stolk&DeHoop2005).
Thedata,ineq.(10),canberegardedasthesolutiontoa(Cauchy)initialvalueproblem:atthesurface(z =0)thesinglescattered
wavefieldurepresentsthedatad,thatis,d(s,r,t)=Q∗ (0)Q∗ (0)u(s,r,t,z=0).Forz>0,u=u(s,r,t,z)representsthe(decomposed)
U,s U,r
wavefield—treatedas‘virtual’data—asafunctionoftime,generatedatdepthzbelowthesurfacebyfictitioussourcessandobservedby
fictitiousreceiversr:
(cid:7)
∞
u(s,r,t,z)= (H(z,z(cid:5))QU,s(z(cid:5))QU,r(z(cid:5))Dt2R(.,z(cid:5)))(s,r,t)dz(cid:5).
z
Figure2. DSRpropagator(eq.11).Dataatdepthz(cid:5)arepropagatedtowardsthesurfaceatz=0;s(cid:5)andr(cid:5)indicatefictitioussourceandreceiverlocations
whereassandrindicateactualsourceandreceiverlocations.Thesubsurfaceoffsethandmidpointxoccurringintheangletransformareindicatedaswell.
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS
1338 M.V.deHoop,R.D.vanderHilstandP.Shen
Then,usolvesaninhomogeneousDSRequation
(cid:3) (cid:4)
∂
∂z −iBU(s,z,Ds,Dt)−iBU(r,z,Dr,Dt) u=g, g= QU,s(z)QU,r(z)Dt2R(.,z). (12)
Thisequationissolvedinthe(upward)directionofdecreasingz(withvanishinginitialconditionforsomelargez,say,atthebottomofthe
model).
TheDSRequationplaysarolelaterintheevaluationofthesensitivitykernelforreflectiontomography.Infact,H isthepropagator
associatedwitheq.(12),andishencereferredtoastheDSRpropagator.Forlateruse,weintroducetheDSRoperator
C(s,r,z,D ,D ,D)= B (s,z,D ,D)+B (r,z,D ,D), (13)
s r t U s t U r t
sothateq.(12)attainstheform(∂ −iC(s,r,z,D ,D ,D))u=g.TheDSRoperatorCdependsonthebackgroundmodel,andapertubation
∂z s r t
δc inc yieldsaperturbationδCinC,whichisderivedfromtheresultsinAppendixA,SubsectionA2.8
0 0
3.2 Downwardcontinuationandimagingofreflectiondata
TheadjointoperatorH(0,z)∗isusedtopropagatethedatabackwardtodepthz.Weconsider
ψd (cid:12)→D= H(0,z)∗Q∗ (0)−1Q∗ (0)−1ψd, (14)
U,s U,r
whereψ isa(source,receiver,andtime)taperthatsuppressesthepartsofthedatathatonedoesnotwanttoconsiderintheimagingprocess
(suchasdataassociatedwithturningrays).Ineq.(14),Disafunctionof(s,r,t,z);fort >0,D(s,r,t,z)canbeidentifiedwithu(s,r,t,z)
ineq.(12),andissometimesreferredtoasa‘sunkensurvey’.Indeed,Dcanbeobtainedbysolvingtheevolutionequation
(cid:3) (cid:4)
∂
−iC(s,r,z,D ,D ,D)∗ u=0, (15)
∂z s r t
inthedirectionofincreasingz(downward),subjecttotheinitialconditionu(s,r,t,z=0)=−Q∗ (0)−1Q∗ (0)−1ψd.Inthepropagating
U,s U,r
regime,C(s,r,z,D ,D ,D )∗=C(s,r,z,D ,D ,D ).
s r t s r t
FollowingClaerbout(1985),animageofδc—thatisthesingularwavespeedvariations(includingreflectors)—canbeobtainedfromD
byapplyingimagingconditions:
(cid:11) (cid:12)
I(x,z)=D s =x− h2,r =x+ h2,t,z |h=0,t=0, (16)
wherewehaveintroducedsubsurfacehorizontalmidpoint-offsetcoordinates,(x,h).Ourtomographyaimstoestimatewavespeedvariations
inthebackground,butweuseimagegatherstoevaluateassociatedmodelupdates.Theseimagegathersareobtainedbyreplacingtheimaging
conditionbyanangletransform(seethenextsection).
4 WAVE-EQUATION ANGLE TRANSFORM AND DATA ANNIHILATORS
BuildingontheworkbyStolk&DeHoop(2001)andBrandsberg-Dahletal.(2003)wediscussherehowredundancyinthedatacanbe
exploitedforthepurposeoftomography.Weassumethatwehavedatafrommultiplesourcesandmultiplereceivers(seeFig.1,bottom).In
theimagingprocess,thereflectionpointbecomestheimagepoint,andtheredundancyinthedatabecomesmanifestinthemultipleimages
thataregenerated(atthatimagepoint)fromthedifferentscatteringanglesandazimuths.Thisisaccomplishedbytheangletransform.In
thistransform,weencounteravariablepthatisessentiallyhalfthedifferencebetweenthehorizontalslownessvectorsassociatedwithsource
andreceiverraysatthedepthwherethetransformisapplied.TherelationbetweenpandscatteringangleandazimuthisdescribedbyDe
Hoopetal.(2003,eqs88–90);byvirtueofthisrelation,inimaging,theangletransformcanbethoughtofasareplacementofthegeneralized
Radontransform.
Foranacceptablewavespeedmodeltheimagesareindependentofscatteringangleandazimuth.Equivalently,thewavespeedmodelis
acceptableifthedataarepredictableby(inotherwords,theyareintherangeof)themodellingoperatorineq.(10).Inthissection,wederive
theannihilatorsusedtoevaluatethiscriterion.(Weemphasizethatthenotionthatamodelisacceptabledoesnotimplythatitisuniquely
determined;infact,theretypicallyisaclassofacceptablemodels.Itcanbeshown,however,thatthehigherthescattererdensity,thesmaller
thecollectionofacceptablemodels).
8Forcompletenesssake,inAppendixBweshowthatgeneralizedscreenexpansionsforQU,s(z)andQU,r(z)leadtowavespeedderivativesinafashionsimilar
totheonesforBU.Thesederivativesaresuppressedinthetreatmentpresentedhere,however.
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS
Onwave-equationreflectiontomography 1339
4.1 Angletransformandcommonimage-pointgathers
With,asbefore,u=u(s,r,t,z),andhthehorizontaloffsetbetweenrands(Fig.2),weintroduce
(cid:7) (cid:3) (cid:4)
h h
(Ru)(x,z,p)= u x− ,x+ ,ph,z χ(x,z,h)dh
IRn−1 (cid:7) (cid:7) 2 (cid:3) 2 (cid:4)
1 h h
= uˆ x− ,x+ ,ω,z exp(−iωph)χ(x,z,h)dhdω, (17)
2π IRn−1 2 2
(notethezerointercepttime:t =0+ ph,withphamultiplicationforn =2andadotproductforn ≥3;thistransformwasoriginally
introducedbyDeBruinetal.1990).Here,χ(x,z,h)isataperinh,suchthatχ(x,z,0)=1.WeapplyRtothedownwardcontinueddata,
D(s,r,t,z)(cf.eq.14),toobtainthewave-equationangletransform,A,
Ad = RH(0,z)∗Q∗ (0)−1Q∗ (0)−1ψd. (18)
U,s U,r
Theresultofthistransformappliedtothedatais
¯I(x,z,p)=(Ad)(x,z,p)=(RD)(x,z,p). (19)
Foreachx,¯I(x,z,p)isacommonimage-pointgatherin(z,p).9Withanappropriatechoiceofχthecommonimage-pointgathersareartefact
free(thatis,freeoffalseevents),eveninthepresenceofcaustics(Stolk&DeHoop2001,2006).Asmentionedabove,thedependenceon
preflectstheredundancyinthedata(adimensioncountconfirmsthis:(s,r,t)∈IR2n−1 while(x,z)∈IRn and2n−1−n=n−1isthe
dimensionofp).Foranacceptablemodelc ,Amapsdatadtoap-familyofreconstructionsofδc,andforeachpthesameimage(¯I(x,z,p))
0
ofδc(x,z)isobtained.
4.2 Intermezzo:AngletransformforglobalEarthapplications
InglobalEarthapplications,inparticularinviewofturningrays,theone-waywavepropagatorsintroducedinSection2maybereplacedby
full-wavepropagators.Accordingtoeq.(12),the(fictitious)sourcessandreceiversrareatthesamedepth(z).Forapplicationtoearthquake
data,however,weconsider(clusters)ofearthquakeswithhypocentresatz=z andstationsatz=0(seeFig.1,bottom).Then,wecanreplace
s
thekernel(H(0,z(cid:5)))(s,r,t,s(cid:5),r(cid:5),t(cid:5))ofH(0,z(cid:5))by
(cid:7)
(H(0,z(cid:5)))(s,r,t,s(cid:5),r(cid:5),t(cid:5)):= G(s,z ,t−t(cid:5)−t¯,s(cid:5),z(cid:5))G(r,0,t¯,r(cid:5),z(cid:5))dt¯. (20)
s
IR
Inviewofreciprocitythekerneloftheadjointisthesame,sothattheangletransformbecomes
(cid:3) (cid:4) (cid:7) (cid:3)(cid:7) (cid:3) (cid:4) (cid:3) (cid:4) (cid:4)
h h h h
D x− ,x+ ,t(cid:5),z = G s,z ,−(t(cid:5)−t+t¯),x− ,z G r,0,t¯,x+ ,z dt¯ (ψd)(s,r,t)dsdrdt, (21)
2 2 s 2 2
IR
followedbyeq.(17);here,z >z foralls.WecanuseanynumericalmethodtocomputetheGreen’sfunctions,forexamplenormalmode
s
summation,leadingthentoa‘normalmodeangletransform’.Eq.(21)canberewritten,uponachangeofvariableofintegrationt¯¯=−(t¯−t),
intheform
(cid:3) (cid:4) (cid:7) (cid:16) (cid:3) (cid:4)(cid:16)(cid:7)
h h h
D x− ,x+ ,t(cid:5),z = G s,z ,t¯¯−t(cid:5),x− ,z
2 2 s 2
(cid:3) (cid:4) (cid:17) (cid:17)
h
G r,0,−(t¯¯−t),x+ ,z (ψd)(s,r,t)drdt dt¯¯ ds; (22)
2
att(cid:5)=0andh=0,thisformrevealsthestructureofshotrecordmigration,whileonerecognizesthenotionofdoublefocussing(Berkhout
1997;Thorbecke1997).
Essentially, eq. (17) is a special case of beamforming with the downward continued data (cf. 18) in sources and receivers (see also
Scherbaumetal.1997).
4.3 Annihilatorsofthedownwardcontinueddata
Sincetheoutcome(Ad)(x,z,p)shouldbeindependentofpwecandefineannihilatorsW,whoseactiononthedatad istoyieldzero,as
follows:fori =1,...,n−1,W =(∂)−1(cid:9)A−1(cid:10) ∂ A(where(cid:9)A−1(cid:10)indicatesaregularizedinverseofA),whichindeedyieldsWd=0.We
considertheannihilatorsnotintiheda∂ttabutinth∂epiimagedomain.Therefore,removingthemapping(cid:9)A−1(cid:10)fromimagegatherstiodata,we
considerthecompanionoperators ∂ A(∂)−1(notethat(∂)−1hastoactinthedatadomain),which,whenappliedtoeq.(17),bringsouta
∂pi ∂t ∂t
(cid:13)
9Savaetal.(1999)usearelatedbutdifferenttransform,viz.,¯¯I(x,Z,p)= IRn−1D(x− h2,x+ h2,t,z)χ(x,z,h)|t=0,z=Z+phdh.UnlikeA,thistransform
cannotbecast,withappropriatelychosenweights,intoanestimationofthereflectioncoefficientinducedbyc0,δc.
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS
1340 M.V.deHoop,R.D.vanderHilstandP.Shen
multiplicationbyfactorh :
(cid:7) (cid:3)i (cid:4)
h h
(R(cid:5)u)(x,z,p)= u x− ,x+ ,ph,z h χ(x,z,h)dh
i IRn−1 (cid:7) (cid:7) 2 (cid:3) 2 i (cid:4)
1 h h
= uˆ x− ,x+ ,ω,z exp(−iωph)h χ(x,z,h)dhdω. (23)
2π IRn−1 2 2 i
Theannihilationofthedataisthusreplacedbyanannihilationofthesetofsubsurfaceimagegathers,R(cid:5)D,parametrizedbyp.
i
Forthepurposeoftomographywealsoneedtheadjoint(R(cid:5) )∗ofR(cid:5) .LetIdenoteatrialimageasafunctionof(x,z,p),then
i i
(cid:7) (cid:3) (cid:7) (cid:7)
1
(cid:9)Ri(cid:5)u,I(cid:10)(x,p) = I(cid:3)R2n−2 2π IRn−1 (cid:4) (cid:4)
h h
uˆ x− ,x+ ,ω,z exp(−iωph)h χ(x,z,h)dhdω I(x,z,p)dxdp
2 2 i
(cid:7) (cid:7)
1
= uˆ(s,r,ω,z)
2π IR2n−2
(cid:3)(cid:7) (cid:4)
(cid:11) (cid:12)
(r−s) exp[iωp(r−s)]I 1(s+r),z,p dp χ(1(s+r),z,r−s)dsdrdω
IRn−1 i 2 2
=(cid:9)uˆ,(Rˆi(cid:5))∗I(cid:10)(s,r,ω) forgivenz, (24)
sothat
(cid:7)
(cid:11) (cid:12)
(Rˆ(cid:5))∗I(s,r,ω,z)= (r−s) exp[iωp(r−s)]I 1(s+r),z,p dp. (25)
i IRn−1 i 2
Here, indicatescomplexconjugation.Wenotethat(Rˆ(cid:5))∗,ineq.(24),yieldsanextensionofa(differentiated)imagegathertofictitious
i
data,inthenatureofRdefinedineq.(9).Ineq.(24)wechangedvariablesofintegrationfrom(x,h)to(s,r).Byremovingthefactorh inthe
i
integrand,weimmediatelyobtainanexpressionforRˆ∗.ByParseval’stheorem,wealsofindR∗,
(cid:7)
(R∗I)(s,r,t,z)= δ(t− p(r−s))I(1(s+r),z,p)dp
IRn−1 2
withthepropertythat(cid:9)Ru,I(cid:10)(x,p) = (cid:9)u,R∗I(cid:10)(s,r,t).NoticethatneitherRnor Ri(cid:5) dependsonc0 and,thus,thattheyareinsensitivetosmooth
perturbationsδc inc .
0 0
5 AN OPTIMIZATION PROCEDURE FOR REFLECTION TOMOGRAPHY
Withthereflectiontomographydevelopedhereweaimtoestimatethebackgroundmedium(c +δc ),withasameasureofsuccessthe
0 0
annihilationofimagegathers.Indeed,ifannihilationwithoperator R(cid:5) appliedtothedownwardcontinueddataineq.(23)occurs,thenthe
i
backgroundmediumisconsideredacceptable.
5.1 Costfunctional
Withthisnotion,themodelestimationiscastintotheminimizationofthefunctional
(cid:7) (cid:7)
(cid:10) (cid:18) (cid:18)
J[c ]= 1 (cid:18)(R(cid:5)H(0,z)∗Q∗ (0)−1Q∗ (0)−1ψd)(x,z,p)(cid:18)2dxdzdp, (26)
0 2 i U,s U,r
i
whichcorrespondswiththeeffectivedataannihilationintegratedoverallimage(scattering)points(x,z)and‘angles’p.Here,H(0,z)∗(and
ψ)dependonc andaresensitivetoperturbationsδc ;however,sinceweassumedthatδc vanishesnearz=0,weneednotconsiderwave
0 0 0
speedderivativesofQ∗ (0)−1Q∗ (0)−1.Comparedwithconventionaltomography,themeasureoftraveltimemismatch(orwaveformcross
U,s U,r
correlation)hasthusbeenreplacedbytheoveralleffectoftheannihilationoperatorR(cid:5).
i
5.2 Gradientofcostfunctional;sensitivitykernel
AstandardmethodforoptimizationcanbeinvokedtocarryouttheminimizationofJ.Here,wediscussamethodforevaluatingtheFre´chet
kernelorgradientofJ.Themethodforevaluatingthegradientofafunctionalderivedfromthesolutionofapartial(orpseudo-)differential
equationisknownastheadjointstatemethod(Wunsch1996).IthasbeenintroducedinseismologybyTarantola(1987)andusedbymany
others.Thegradientisrequiredfortheoptimization,whereastheFre´chetorsensitivitykernelprovidesinformationabouttheresolutionof
ourapproachtowave-equationreflectiontomography.
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS
Onwave-equationreflectiontomography 1341
TheperturbationofthefunctionalJ,underasmoothperturbationδc ofc ,isderivedfrom
(cid:19) (cid:5) (cid:6) 0 0
(cid:7) (cid:7)
(cid:10)
δJ[c ]= (Rˆ(cid:5))∗R(cid:5)H(0,z)∗Q∗ (0)−1Q∗ (0)−1ψd (s,r,ω,z)
0 i i U,s U,r
i
(cid:20)
(cid:11) (cid:12)
(δH 0,z)∗Q∗ (0)−1Q∗ (0)−1ψd ˆ(s,r,ω,z)dsdrdω dz, (27)
U,s U,r
and evaluated in Appendix C. The evaluation procedure, which is common to many different imaging schemes, consists of solving two
evolutionequationsfollowedbyanimagingoperation.Weusetheobservationmadebeloweq.(14)andrelateDtosolutions,u,oftheDSR
equation.Thefirstevolutionequationiseq.(15),whichissolvedinthedirectionofincreasingz(downward);inthefrequencydomainthis
equationreads:
(cid:3) (cid:4)
∂ (cid:11) (cid:12)
−iCˆ(s,r,z,D ,D ,ω)∗ uˆ =0, uˆ(s,r,ω,z=0)=− Q∗ (0)−1Q∗ (0)−1ψd ˆ(s,r,ω). (28)
∂z s r U,s U,r
Thesecondistheadjointfieldequation,
(cid:3) (cid:4)
∂ (cid:10)
−iCˆ(s,r,z,D ,D ,ω) vˆ = (Rˆ(cid:5))∗R(cid:5)u, (29)
∂z s r i i
i
whichissolvedinthe(upward)directionofdecreasingz(withvanishinginitialconditionforsomelargez,say,atthebottomofthemodel).
Theright-handsidequantifiesthefailureofannihilationandrepresentsamism(cid:13)atchsourcedistribution.
Thesolutionsuˆ andvˆ combinedformthekernel(or‘image’)KinδJ = K(x,z)δc (x,z)d(x,z),andisgivenby(cf.C12andA7)
⎡ 0
(cid:7) (cid:7)
1 ∞ (cid:10)N ⎢
K(x,z)∼ Re dω(−i) ⎣ uˆ(x,r,ω,z)S(cid:5)(x,z)F−1 A (σ,ω,z)F vˆ(.,r,ω,z)dr
π (cid:24)j σ→x(cid:25)(cid:26)j s→σ(cid:27)
0 j=0
generalizedscreen ⎤
(cid:7)
⎥
+ uˆ(s,x,ω,z)S(cid:5)(x,z)F−1 A (ρ,ω,z)F vˆ(s,.,ω,z)ds⎦. (30)
(cid:24)j ρ→x (cid:25)(cid:26)j r→ρ(cid:27)
generalizedscreen
Here,S andA appearinthegeneralized-screenexpansionofthesingle-square-rootoperator,whichisexplainedindetailinAppendixA(cf.
j j
A1andA4).Weusedthesymmetryinfrequencytorestricttheevaluationtopositivevalues.InvokingthemodelrepresentationinSection2.4
yieldsδJ =(cid:9)K,δc0(cid:10)(x,z)=(cid:9)(δck),(Ec∗K)(cid:10)IRM.
Eq.(30)isoftheformofatimecrosscorrelationofthedownwardcontinueddataandtheadjointfieldexcitedbyamismatchforce
nestedinageneralizedscreenoperation.Itdiffersfromastandardimagingconditioninparticularthroughtheintegrationsoversandr.The
proceduretoevaluateeq.(30)isillustratedinFig.3.Itconsistsofthefollowingsteps:
(i) startingatthetop,downwardcontinuestep-by-stepthedata(inthefrequencydomain)allthewaytothebottom(cf.14)whilestoring
theresultsatallintermediatedepths;
(ii) startingatthebottom,evaluatethesuccessofannihilation,thatis,themismatchsource,withthedownwardcontinueddata,
Figure3. Theevaluationofthekernel.Thesingleanddoublearrowsrefertothetwotermsinbetweenbracketsineq.(30);thegreyarrowsindicatethe
generalizedscreenoperations.Thefieldvisupwardcontinuedwithsourcesdeterminedbythefailuretoannihilatethedataintheimagedomain;urepresents
thedownwardcontinueddata.Atanydepth,z,ascreenactionisappliedtov(indicatedbygreyarrows)inparallelinthesourcevariable(left)andinthereceiver
variable(right).Then,inparallel,theinnerproducts(includingtimecorrelation)aretakeninreceivers(left)andsources(right).Finally,thetwocontributions
areaddedtogether.
(cid:2)C 2006TheAuthors,GJI,167,1332–1352
Journalcompilation(cid:2)C 2006RAS