Waveform-based simulated annealing of crosshole transmission data: a semi-global method for

22 Pages·2014·7.38 MB·English

Geophysical Journal International
Geophys.J.Int.(2014)199,1586–1607 doi:10.1093/gji/ggu307
GJISeismology
Waveform-based simulated annealing of crosshole transmission data:
a semi-global method for estimating seismic anisotropy
Michael V. Afanasiev,1,∗ R. Gerhard Pratt,2 Rie Kamei3,∗ and Glenn McDowell4
1InstituteofGeophysics,ETHZu¨rich,Zu¨rich8053,Switzerland.E-mail:[email protected]
2DepartmentofEarthSciences,UniversityofWesternOntario,London,ON,CanadaN6A5B7
3SchoolofEarthandEnvironment,UniversityofWesternAustralia,Crawley,WA6009,Australia
4ValeTechnologyDevelopment(Canada)Ltd.,Sudbury,Ontario,Canada
D
o
Accepted2014August4.Received2014July30;inoriginalform2014January9 w
n
lo
a
d
e
SWUesMuccMessAfuRllyYapplythesemi-globalinversemethodofsimulatedannealingtodeterminethe d from
best-fitting1-Danisotropymodelforuseinacousticfrequencydomainwaveformtomography. http
Ourforwardproblemisbasedonanumericalsolutionofthefrequencydomainacousticwave s
equation,andweminimizewavefieldphaseresidualsthroughrandomperturbationstoa1-D ://a
c
a
verticallyvaryinganisotropyprofile.Bothrealandsyntheticexamplesarepresentedinorder d
e
m
todemonstrateandvalidatetheapproach.
ic
Fortherealdataexample,weprocessedandinvertedacross-boreholedatasetacquiredby .o
u
p
ValeTechnologyDevelopment(Canada)Ltd.intheEasternDeepsdeposit,locatedinVoisey’s .c
o
Bay, Labrador, Canada. The inversion workflow comprises the full suite of acquisition, data m
/g
processing,startingmodelbuildingthroughtraveltimetomography,simulatedannealingand ji/a
finallywaveformtomography. rtic
Waveformtomographyisahighresolutionmethodthatrequiresanaccuratestartingmodel. le
-a
A cycle-skipping issue observed in our initial starting model was hypothesized to be due b
s
to an erroneous anisotropy model from traveltime tomography. This motivated the use of tra
c
simulatedannealingasasemi-globalmethodforanisotropyestimation.Weinitiallytestedthe t/1
9
simulatedannealingapproachonasyntheticdatasetbasedontheVoisey’sBayenvironment; 9/3
thesetestsweresuccessfulandledtotheapplicationofthesimulatedannealingapproachto /1
5
8
therealdataset.Similarbehaviourwasobservedintheanisotropymodelsobtainedthrough 6
/6
traveltime tomography in both the real and synthetic data sets, where simulated annealing 0
9
producedananisotropymodelwhichsolvedthecycle-skippingissue.Intherealdataexample, 07
8
simulatedannealingledtoafinalmodelthatcompareswellwiththevelocitiesindependently b
y
estimated from borehole logs. By comparing the calculated ray paths and wave paths, we g
u
e
attributed the failure of anisotropic traveltime tomography to the breakdown of the ray- s
theoreticalapproximationinthevicinityofstrongvelocitydiscontinuities. t on
0
9
Key words: Inverse theory; Controlled source seismology; Seismic anisotropy; Seismic A
p
tomography;Computationalseismology;Wavepropagation. ril 2
0
1
9
tomographyistoestimateasetofelasticEarthparameters(usually
1 INTRODUCTION
velocities)thatminimizethedifferencesbetweenasetofobserved
Seismictomographyisaversatilesubsurfaceimagingmethodap- andpredictedseismictraveltimes.Althoughtraveltimetomography
propriate for use with transmitted seismic waveforms on a wide has been used widely and successfully, the spatial resolution of
rangeofscales.Dataarecommonlyprocessedusingtraveltimeto- ray-theoreticalmethodsarelimited,andstructuressmallerthanthe
mography, often based on a ray-theoretical approximation of the orderofthewidthofthefirstFresnelzoneareeffectivelyinvisible
elasticoracousticwaveequation(Cˇerveny1972;Cˇerveny&Jech (Williamson1991;Williamson&Worthington1993).
1982;Zelt&Barton1998;Zeltetal.2003).Thegoaloftraveltime Asthepowerofcomputershasincreased,ithasbecomefeasible
tomovebeyondtheasymptotic,high-frequencyrayapproximation,
andinsteadtobasetomographicmethodsonnumericalmodelsof
∗Formerlyat:DepartmentofEarthSciences,UniversityofWesternOntario, seismic wave propagation via the elastic or acoustic wave equa-
London,ON,CanadaN6A5B7. tions.Tomographicmethodsusingthesetechniquesarecommonly
1586 (cid:3)C 2014TheAuthorsPublishedbyOxfordUniversityPressonbehalfoftheRoyalAstronomicalSociety
Estimatinganisotropybysimulatedannealing 1587
referredtoas‘waveformtomography’(Brenders&Pratt2007a,b), thehalf-cyclecriterioninourstartingmodelsforwaveformtomog-
andarebasedonadjointmethodsforwaveforminversionfirstin- raphy.Itisthepurposeofthispapertoshowthatsemi-globalinverse
troduced by Lailly (1983) and Tarantola (1984). As in traveltime methodscanovercomethecycle-skippingphenomena.Todothis,
tomography,inwaveformtomographyweseektoiterativelyupdate wereducethedimensionalityofthemodelspacetoalevelwhich
an earth model with the goal of producing an accurate image of permitsthetime-efficientcomputationofmanyforwardmodels.
subsurfaceelasticparametersthatminimizesthemisfitbetweenob- Since an anisotropy model which fits the data to the degree
served and predicted data. Whereas in the traveltime method the required for waveform tomography could not be found through
dataarereducedtoasetoftraveltime‘picks’,thedatainwaveform manual picking and traveltime tomography, a novel method was
tomographyareeditedversionsoftherecordedseismicwaveforms developedtodiscovera‘best-fitting’1-Danisotropymodelusing
themselves.Furthermore,becausewaveformtomographyreplaces waveformdataratherthantraveltimedata.Wepresentasemi-global
ray methods for data prediction with numerical modelling of the method, waveform-based simulated annealing, as an intermediate
waveequation,thisnaturallynegatestheneedtoapproximatewave stepthatallowsustoextractalayer-by-layerestimateoftheseismic
propagationasoccurringalonganinfinitelythinraypath,andcom- anisotropyofsufficientaccuracytoallowforsuccessfulwaveform
plex,heterogeneousmediaisaccuratelyaccountedfor. tomography. The traveltime tomography method for determining
In implementation, waveform tomography redistributes wave- anisotropywasabandoned,andattentionwasgiventoreducingthe D
o
formresidualsoverspatiallybroad‘wavepaths’(Woodward1992). waveformresidualsthroughperturbationstotheanisotropymodel w
n
Thisisamorecorrectmodelofwavepropagationphysicsthatleads usingsimulatedannealing(whilepreservingthehorizontalvelocity lo
a
to an improved resolution limit on the order of one half of the modelobtainedfromthetraveltimes). d
e
sceyicslmiciaclwnaatvuerleenogftwhsavoefftohremddaatata(Wtheuw&avTeofkosromzt1o9m8o7g).raDpuheytporothbe- relTevhainstpinavpeerrsiisonomrgeatnhiozdeds,ainscflouldlionwgst:hoWseeubseegdininwaintihsoatrroepviicewtravo-f d from
lemishighlynon-linear,andasaconsequenceofthelocalinverse eltimetomographyandinseismicwaveformtomography.Aftera h
methodscommonlyusedthesuccessfulapplicationofthetechnique brief discussion of the geological setting of the Voisey’s Bay de- ttps
requiresastartingmodelwhichpredictsfirstarrivalwaveformsto posit,wedescribethedataacquisition,instrumentationandsurvey ://a
c
within one-half cycle of the lowest frequency present (Sirgue & procedures for the seismic data set. We then present some sam- a
d
Pratt 2004). This is known as the ‘half-cycle criterion’, and it is ple data gathers and comment on data quality, and we illustrate em
a condition that is often difficult to satisfy in real data examples, theeffectsofpre-processing.Followingthis,wediscussinitialve- ic
.o
inwhichlowfrequencydatamaybedifficulttoacquire.Failureto locityandanisotropymodelsfoundthroughanisotropictraveltime u
p
satisfythisconditionleadsto‘cycle-skipping’duringtheinversions. tomography,andthecycle-skipphenomenonisidentifiedinthese .c
o
CrossholeseismictomographyhasbeenusedbyValeTechnol- results.Afterhypothesizingthattheresultsareduetoanerroneous m
/g
ogyDevelopment(Canada)Ltd.(referredtoasValehereafter)to anisotropy model, we present a synthetic study in which results ji/a
provide geophysical images of massive sulphides in the Eastern of a similar nature are reproduced via numerical waveform sim- rtic
Deeps deposit at Voisey’s Bay, Labrador, Canada (Enescu et al. ulations. We demonstrate that traveltime tomography may indeed le
-a
2002;Perozzietal.2012).Thepurposeofthesesurveysistoobtain produceincorrectanisotropymodels,andwethenillustratehowwe b
s
information on ore-body geometries, and to use this information, successfully overcome this issue with the waveform-based simu- tra
c
alongwithanyrelevantgeologicalinformation,toplacespatialcon- latedannealingapproachwehavedeveloped.Thisisfollowedby t/1
straintsonthelocationandextentofeconomicores.Thegoalofthis asuccessfulapplicationofthesimulatedannealingmethodtothe 9
9
particular study is to provide Vale with a high resolution seismic original data set, along with a validation of our results by com- /3
/1
velocity map of the target region via waveform tomography, and parisonsbetweentheinvertedvelocitiesandindependentborehole 5
8
6
to investigate the feasibility of waveform tomography in similar velocitymeasurements,andacomparisonwithcorelogs.Wefinish /6
hardrockenvironments.Thegeologicalsettingishighlyheteroge- withadiscussionofthephysicallimitationsofraytheoryandtrav- 09
0
neous,withstrongvelocitycontrastsandvariablelevelsofseismic eltimetomographyinastronglyheterogeneous,anisotropicmodel 7
8
anisotropy (Enescu et al. 2002; Cosma & Enescu 2003); Perozzi inordertobetterunderstandthebreakdownofthemethodinthis b
y
etal.(2012)performedtraveltimetomographyonthesamedataset example. g
u
presentedinthispaper,howevertheydidnotconsideranisotropyto es
beafactor. t o
n
Seismicvelocityanisotropycanbesignificantinhardrockset- 0
9
tings (Pratt et al. 1993). Cosma & Enescu (2003) estimate the 2 METHODS A
p
anisotropy within the host rock at Voisey’s Bay at approximately ril 2
3percent,whilelaboratorymeasurementshavereportedontrocto- 2.1 Inversemethods 0
1
litesshowinganisotropyashighas13percent(Babuska1968),with 9
We begin with a quick review of inverse techniques, followed by
Iturrinoetal.(1991)reportingvaluesashighas6percent.While
anintroductiontothespecifictechniquesrelevanttotheworkpre-
these laboratory measurements are very accurate, their relevance
sented.Inthemostgeneralsense,Tarantola(2005)definesthesolu-
totrue,insituanisotropyparametersatcrossholefrequenciesare
tionoftheinverseproblemasaprobabilitydistributionP(m∈M)
uncertain,asinsitupressuresmaycloseanyfractureplaneswhich
overthespaceofallmodelsM,withtheform
are present at atmospheric pressure (Tsvankin 2012). Regardless,
(cid:2) (cid:3)
whenapplyinghigh-resolutionseismicimagingtechniquestothis 1
area, the proper quantification of anisotropy at the relevant scale P(m∈M)=exp − E(m) . (1)
2
lengthisessential.
Thisgeologicalsettingpresentsparticularchallenges,including Here, E(m), known as the ‘objective’ or ‘misfit’ function, is an
significantvelocitydiscontinuities.Thisleadstothefailureofray- attempt to objectively measure of the quality of the model m. A
basedtraveltimetomographyforthiscasestudytocorrectlyextract common way to define this function is by the L norm of the
2
a useful anisotropy model, initially preventing us from satisfying differences between observed data (d) and estimated data g(m),
1588 M.V.Afanasievetal.
combined with appropriate a priori information about the function to direct the search. While extremely robust, global and
modelm semi-globalinversemethodshavenotbeenwidelyadoptedintrav-
0
eltimeorwaveformtomography,althoughexamplesdoexist(Sen
E(m)=[g(m)−d]TC−1[g(m)−d]+(m−m )TC−1(m−m ) (2)
d 0 m 0 &Stoffa1991;Varelaetal.1998;Sambridge&Mosegaard2002;
(after Tarantola 2005). The notation g(m) is used in eq. (2) to Kaufl et al. 2013). This is mainly due to the extreme size of the
symbolize forward modelling, which produces an estimated data parameterspace,andthesignificantamountoftimeneededtoper-
setbyapplyingafunctionaloperatorgtoasetofmodelparameters. formtheforwardmodelling.Acomprehensivereviewofglobaland
Wehaveincludedameasureofmodelanddatauncertaintythrough semi-globalinversemethodswithinageophysicalcontextisgiven
theincorporationofmodelanddatacovariancematrices,C and bySambridge&Mosegaard(2002).
m
C ,respectively.
d
The probabilistic interpretation of inverse problems requires a
fulldescriptionofthebehaviouroftheobjectivefunctionoverthe 2.1.3 Simulatedannealing
space of all possible models. A less ambitious but highly useful
Of specific relevance to this paper is the semi-global method of
goalistofindasetofmodelsthataremostlikelygiventhedata,
simulatedannealing(Rothman1985,1986;Kirkpatricketal.1987;
andgiventhepriormodel.Toaccomplishthis,wetrytoexplorethe D
Mosegaard&Vestergaarad1991;Sen&Stoffa1991;Varelaetal. o
maximumofeq.(1),whichcorrespondstoexploringtheminimum w
of E(m). There are a variety of approaches that can be used to 1998).Theconceptsbehindsimulatedannealingareborrowedfrom nlo
minimizeE(m)ineq.(2).Thesemaybeclassifiedinto statisticalmechanics,specificallyfromtheanalysisofasolidcool- ad
e
ing in a heat bath. The method attempts to sample a probability d
(i) Localmethodsthatsearchintheneighbourhoodofastarting distributionoftheform fro
m
mo(diie)lGfolorbnaelarobrysmemini-imglao,baanldmethodsthatattemptamoregeneral P(m)=exp[−E(m)/T], (4) http
s
searchofthemodelspacetoidentifyglobalminima,whichmaybe wheretheprobabilityassignedtoacertainmodelP(m)isdependent ://a
unreachablewithlocalmethods. c
onthecurrentvalueofthemisfitE(m),andacontrolparameterT. a
d
There are several ways to generate samples of the distribution in em
2.1.1 Localinversemethods eq. (4) (Sen & Stoffa 1995). The one we adopt is the Metropolis ic.o
algorithm(Metropolisetal.1953),whichproceedsasfollows:First, u
p
In local inverse methods, updates to an initial model are usually arandomperturbation,(cid:3)mismadetoasingle,randomlychosen .c
o
foundbyusinginformationonthelocalgradientofthemisfitfunc- modelparameterbelongingtomodelm(k),resultinginm(k+1)where m
tuisoinng(Plirnaettaeritzaeld.1sy9s9t8e)m,asnodfmthoedfeolrmupdatesareiterativelyestimated m(k+1)=m(k)+(cid:6)m;(cid:4)m(k),m(k+1)∈M(cid:5). (5) /gji/artic
m(k+1)=m(k)−α(k)∇mE(m(k)), (3) ThemisfitfunctionE(m)isevaluatedform(k)andm(k+1).Depend- le-a
wherekreferstotheiterationnumber,andαisknownasthe‘step ianpgroonbatbhielivtyalAue(smo(kf+E1)(,mT(k()k)+a1)n)dwEh(emre(k+1)),m(k+1)isacceptedwith bstra
lenWgthhi’l,eatqhueraenatirteymthaantymvuasrtiabteioensstiomnaeteqd.(a3s)a,aplalrlotcoaflthinevperrosecemdeutrhe-. A(cid:6)m(k+1), T(k+1)(cid:7) ct/19
9
odsrelyonthesameconcept:thatofsuccessiveminimizationsto (cid:8) (cid:9) (cid:10) (cid:6) (cid:7) (cid:6) (cid:7) /3
the misfit function along the gradient direction, with the require- = exp −E(m(k+T1)(k)+−1E)(m(k)) if E(cid:6)m(k+1)(cid:7)> E(cid:6)m(k)(cid:7) (6) /158
mentthattheinitialmodelisstrictlywithinthe‘zoneofattraction’ 1 if E m(k+1) ≤ E m(k) 6/6
oftheglobalminimumofE(m). 0
9
ThisrecipeforacceptanceisknownastheMetropoliscriterion.We 07
canwritetheprobabilityP(mj),ofacceptinganymodelmj ∈M 8 b
2.1.2 Globalorsemi-globalinversemethods as y g
u
AGnloablatlerinnavteivrseetomleotchaoldisnvienrssteeamdettrhyodtosafirendgltohbealbiensvtefirstemmoedthelodbsy. P(mj, T(k))=G(mj)×A(mj,T(k)) est o
einxfpolromriantgiotnh.eTmhoedeslimsppalceestMexwamithplliettliesraefgerreidncseetaorclho,cawldheicrhivatetisvtes G(mj)= M1 , (7) n 09 A
p
adnisdcrwehteicmhoredteulrsnmsa(kn)einnsMem,bwleheorfemkordeeflesrswthoicthhebietsetrfiattiothnenduamtab(ears, wvehresreelyGo(nmthje)nisumthbeegreonfeerlaetmioenntpsroMbainbiMlity.Tohfemsuj,bsdcerpipetnjdiinsgusiend- ril 20
1
evaluatedbythemisfitfunction).Animportant,pointmustbemade heretorepresentacertainrealizationofmodelparameters,while 9
here.IndefiningP(m)ineq.(1),weinterpretedthesolutionofthe the superscript k retains its role as an iteration counter. Note that
inverseproblemasaprobabilitydistributionoverM.Inthelocal herewearediscretizingandassumingafinitesizeforM,andthat
optimizationroutinesdiscussedaboveweseekasinglesolutionm(k) thegenerationprobabilityisuniform.
whichbestfitsthedata.Thisisbecauseofthedeterministicnature To use simulated annealing with the Metropolis algorithm, we
ofgradient algorithms,which in mostcases sample M much too begin with a high value for T(0), and slowly decrease T(k) over
narrowly to make any rigorous claims about the full distribution the course of the inversion. While T(k) is large, most model per-
P(m).Globaloptimizationmethodsovercomethisproblemasthey turbationswillbeaccepted.AsT(k) decreases,onlyperturbations
attempttosampletheentiremodelspace,andtheirsolutionslead whichdecreaseE(m(k)),orwhichonlyincreaseE(m(k))byasmall
tothenaturalinterpretationofeq.(1)asaprobabilitydistribution. amount, are likely to be accepted. Eventually, once T(k) drops to
Similartoglobalinversemethodsare‘semi-global’inversemeth- a low enough value, it will be overwhelmingly unlikely that any
ods which attempt to search a large portion of the model space, model perturbations that increase the misfit will be accepted, we
but also exploit some information on the structure of the misfit assumeconvergence.
Estimatinganisotropybysimulatedannealing 1589
Revisitingthediscussiononthepitfallsoflocalsearchmethods, (Tsvankin2012).Asingleparameterisusedtoquantifytherela-
notethemainadvantageofthesimulatedannealingmethod:aslong tivevariancesbetweenisotropicandanisotropicparameters.Ifthis
astheinitialtemperatureishighenoughthefinalmodeldoesnotde- parameterissmall,thisallowstheintroductionofanisotropyinto
pendontheinitial(orseed)model.Inthehightemperaturelimit,all the model without significant increases in the objective function
randomlyproposedmodelsareaccepted,andwearetrulyworking E(m). In contrast, if this parameter is large, this is equivalent to
withaglobalinversemethod.However,thealgorithmcomeswith assumingthevarianceoftheanisotropyparametersismuchsmaller
amajorlimitation:ifanyfinite-timecoolingscheduleisused,there than the a priori (isotropic) model variance, and perturbations to
isachancethattheglobalminimumofE(m)willnotbefound.To any anisotropic parameters will thereby be penalized due to their
mitigatethispossibility,weshouldsetthecoolingscheduletobeas significantinfluenceonthevalueofE(m).
slowaspossible.However,practicalrequirements dictatethatthe Transverselyisotropicmediawitheitheratiltedsymmetryaxis
algorithmfinishinareasonableamountoftime.Asacompromise, (TTI) or a symmetry axis aligned with the vertical (VTI) are in-
thecoolingscheduleisusuallysetsothatsimulatedannealingwill vestigatedinthispaper.TheparametersweusefordescribingTTI
completeafteracertainamountofcomputertime,orafteracertain or VTI systems are those proposed by Thomsen (1986): V (the
p
numberofiterations(Sen&Stoffa1995).Bydoingthis,wetrade velocityinthefastdirection),ε(thefractionaldifferencebetween
a guarantee of convergence for a high likelihood of convergence, fastandslowvelocity)andδ(theparameterdescribingtheanellip- D
o
although unless the full model space is searched one cannot say ticity).Notethatalthoughthisapproachallowsforgeneral,2-DTTI w
n
for certain whether the final model obtained is the model which anisotropy,ultimatelyouracousticwaveformtomographyscheme lo
producestheminimummisfit. (describedbelow)islimitedtoVTIsystemswithelliptical(δ=ε) ad
e
adUhoncfoarntudnpartoeblyle,msetdteinpgencdoeonlitn(gSescnh&edSutloefsfain1t9h9i5s;mTaanrannetroilsal2a0rg0e5l)y, aniTsohterodpeygrteheastvoafriferesewdiotmhdaevpatihlaablloenfeo.rtheinversionschemene- d from
althoughNulton&Salamon(1988)andMosegaard&Vestergaarad cessitateaproperstrategyforchoosingvaluesfortherequiredreg- h
(1991) develop implementations which set the cooling schedule ularizationparameters.Pratt&Chapman(1992)suggestedthatthe ttps
accordingtostatisticalinformationgainedduringthemodelsearch. regularizationparametersbechosenaccordingtotheappearanceof ://a
c
AnotherdifficultyistheestimationoftheinitialtemperatureT(0). asuiteoftomograms,andrelaxingtheconstraintsuntilamodelfits a
d
We must ensure that at this temperature the Metropolis criterion the data to within an acceptable margin of error while remaining em
willbepassedforalargevarietyofmodels;otherwise,theinitial geologicallyplausible.Decidingwhatan‘acceptablemarginofer- ic
.o
model will have an impact on the final result. This choice of the ror’iscanbesomewhatambiguous,andPratt&Chapman’s(1992) u
p
initialtemperatureisalsoproblemdependent,althoughastrategy suggestionistotryandfitthedatatowithintheestimatedpicking .c
o
issuggestedinalatersection. error.Asdiscussedbelow,ifthegoaloftraveltimetomographyisthe m
/g
It should be noted that the use of the term ‘simulated anneal- generationofastartingmodelforwaveformtomography,wehave ji/a
ing’does notautomatically implythatmodelsareperturbedonly theadditionalrequirementthattheerrorinthetraveltimeresiduals rtic
oneparameteratatime,orthatthegenerationprobabilityisuni- belessthanone-halfoftheperiodofthelowestfrequencysignal le
-a
formandtemperatureindependent.Indeed,Szu&Hartley(1987) presentinthesurveydata(Sirgue&Pratt2004). b
s
generate perturbations with a temperature dependent probability. tra
c
Thebehaviourofthealgorithmdependsonhowtheperturbations t/1
aretreated.Inourformulation,whereasinglemodelparameteris 2.3 Acousticwaveformtomographyin2-D, 99
perturbedateachiteration,notallpointsinmodelspaceareacces- weaklyanisotropicmedia /3/1
sibleatoneinstant:wemayonly‘crawl’alongoneaxisatatime. 5
8
Asmentionedintheintroduction,thepurposeofwaveformtomog- 6
Ttehrsis,iinswinhcicohntcraassettaomthueltaidpipmroeancshioonfalp‘ejrutumrpb’intghraolulgmhomdeoldpealrsapmacee- raphyistodiscoverasubsurfacemodelthatminimizesthediffer- /609
encebetweenobservedandcalculatedwaveforms.Afulloverview 0
is possible. The choice of a ‘crawling’ or ‘jumping’ formulation 7
islikelyproblemdependent,althoughforproblemswherecertain ofthetheoryforisotropicmediaisgiveninPrattetal.(1998),Shin 8 b
modelparametersarebetterconstrainedthanothersthechoiceof &Cha(2009)andKamei&Pratt(2013). y g
awecreawncloinugntfeorremduilnatoiounrmfiealydsdpaeteadsuept,caonndvethrgeebnecnee.fiTthsisofisatdhoepctiansge forFmoltlhoewriencgotrhdeeddetvimeleo-pdmomenatinofwKavamefieeiledtua(lx.,(2t0)1i2n)t,owtheefiLrsatptlraacnes–- uest o
thecrawlingapproachwillbeexpandeduponinsection5. Fourier domain via a (discrete version) of the complex-valued n 0
LaplaceTransform: 9 A
(cid:11) p
u(x,s)= ∞u(x,t)exp[−st]dt, (8) ril 2
0
2.2 Traveltimetomographyin2-D, 0 19
weaklyanisotropicmedia where
Our anisotropic traveltime tomography approach is based on the s =σ +iω, (9)
work of Chapman & Pratt (1992) and Pratt & Chapman (1992);
thereaderisreferredtothesepapersforadetaileddescriptionof with a real-valued Laplace constant σ and a real-valued angular
themethodology.Herewewillquicklyreviewsomecriticalaspects frequencyω.Eq.(8)canbeexpandedtoyield
(cid:11)
whichareespeciallyrelevanttothispaper. ∞
Pratt&Chapman(1992)advocatedasystemforcontrollingthea u(x,s)= u(x,t)exp[−σt]exp[−iωt]dt. (10)
priorimodelvarianceinwhichthevarianceoftheprior(isotropic) 0
velocitycanbesetindependentlyfromthevarianceoftheparame- Weinterprettheparameterσ astheinverseofacharacteristicdecay
terscontrollinganisotropy.Thistwo-tieredsystemisrequireddueto time
thefactthattheparameterscontrollinganisotropyareoftenmuch 1
more poorly determined than those controlling isotropic velocity τ = σ, (11)
1590 M.V.Afanasievetal.
so that eq. (10) represents the Fourier transform of each time- conditionsthekinematiceffectsoftheanisotropymaybeaccounted
dampedtrace for by implementing a layer-by-layer geometrical stretch propor-
(cid:11) ∞(cid:12) (cid:2) t(cid:3)(cid:13) tionaltotheamountofanisotropyineachlayer,anapproachorigi-
u(x,s)= u(x,t)exp − exp[−iωt]dt. (12) nallysuggestedbyDellinger(1991).Thislimitationrequiresthatwe
τ
0 attempttorepresentan-ellipticanisotropywithanellipticmodel;as
ThepurposeofformulatingwaveformtomographyintheLaplace– suchtheparameterεhastheaddedresponsibilityofapproximating
Fourier domain, instead of just the Fourier domain, is to use the theeffectofδovertherangeofanglesinthesurvey.Theanisotropy
parameterτasadatapre-conditioner(Sirgue&Pratt2004;Brenders isassumedtobefullydescribedinthestartingmodel,andisnot
&Pratt2007a;Kameietal.2013),whichpreferentiallyweightsdata updatedduringthewaveformtomographystages.
correspondingtoearlyarrivals.Wecanfurtherrewriteeq.(12)as
(cid:11)
∞
u(x,s)= u(x,t)exp[−i(cid:9)t]dt (13)
3 ACQUISITION AND INITIAL
−∞
PROCESSING
=u(x,(cid:9)), (14) Themainpurposeoftheexperimentalcrossholestudywastoin- D
o
vestigate the continuity of a mineralized massive sulphide zone w
where(cid:9)=ω−i/τ isacomplex-valuedfrequency.Inthisdevel- locatedinVoisey’sBay,Labrador,Canada.Weusedseismicwave- nlo
opmentweseethatthecomplexLaplaceandFouriertransformsare formdatafromVale’spropertyobtainedfromacrossholeseismic ad
e
e‘Lquapivlaacleen–tF,oaunrdiearsdsoumchaiwnewreavfeerfiteoldth’e(Strhainns&forCmhead2w0a0v9e)fi.eldasthe stiuornv.eTyhdeesdiagtnaesdettoinmclaupdethdethcoensteinisumitiycowfatvheefosurmlphdiadtea,mthineesraoluizrcae- d from
Weemploya‘logarithmicresidual’,or andreceivergeometry,andacousticboreholeloggingresultsfrom h
δdj =ln(cid:14)udj(cid:15) (15) thesourceandreceiverboreholes. ttps://a
j c
a
d
(Shin & Min 2006), where u and d refer to the modelled and e
j j m
3.1 Geologicalandgeophysicalsetting
recorded wavefield, respectively, at receiver j. Expressing the ic
complex-valued data in terms of their amplitude (A) and phase The Voisey’s Bay deposit, located on the East Coast of Canada .ou
components(θ),wecanwrite inLabrador,containssignificanteconomicquantitiesofsulphide- p.c
o
uj = Aujexp[iθuj], (16) h20o0st0e)d. Tnhicekeml,inceoraplpizeerdanddepocosibtsalatroerpesart(Eovfatnhse-LVaomissewyo’sodBaeyt Ianl-. m/gji/a
d = A exp[iθ ], (17) trusiveSuite,whichconsistsprimarilyoflayeredtroctolite-gabbros rtic
j dj dj intrudingintoArcheanquartz-feldspar-biotitegneissesoftheNain le
-a
soeq.(15)becomes province,andProterozoicsulfidicgarnetiferousparagneissesofthe b
s
(cid:16) (cid:17) Churchillprovince(Ryanetal.1995).Therearetwomagmacham- tra
δd =ln Auj +i(cid:6)θ −θ (cid:7). (18) bersassociatedwiththeVoisey’sBaydeposit,knownastheEastern ct/1
j A uj dj DeepsandWesternDeepschambers.Thisprojectdealssolelywith 9
dj the Eastern Deeps magmatic system, consisting of a feeder tube 9/3
This allows fora trivialseparation ofamplitude andphase in the (theEasternDeepsFeeder),whichleadsintothelargeEasterDeeps /15
8
misfitfunctionintherealandimaginarypartsofthedataresiduals. 6
Chamber. The main control on the location of mineralization is /6
Notethatwhenusingthisdevelopment,weassumethatthephases structural,whichimpliesthatdepositioniscontrolledbymorpho- 09
of u and d have been unwrapped to the same degree, meaning 0
j j logical changes in the conduit geometry, and not by gravitational 7
|θuj −θdj|≤π. This is critical, as a perturbation to θuj by 2nπ, settling within the chamber itself (Evans-Lamswood et al. 2000). 8 by
wherenisanyinteger,willleavethemisfitunchanged(duetothe Apartfromtheseareasofmassivemineralization,thechamberis g
u
periodicity of the phase). This leads to the realization that if the dominated by troctolite-gabbros, containing a variable amount of es
predictedandobservedwavefieldsaremorethanπ(onehalf-cycle) sulphidefragments,andgranite-syenite. t o
n
outofphaseintheinitialmodel,anylocaloptimizationmethodisat Thespecificboreholesusedinourcrossholestudyweredrilledin 0
9
riskofconvergencetowardsthewrongcycle.Thisiswhatismeant 2000byBoartLongyear(Inco2000),andtheyweredesignedtoin- A
wcrhiteenriowneirseftheretcoonthdeiticoyncloen-skthipepiinnigtiaplhemnoodmeelntohna;t tihserehqaulfir-ecdyctloe tersectazoneofmassivesulphidemineralizationatdepth,referred pril 2
toasthe‘target’.Amock-upofthetargetareaisschematicallyillus- 0
strictlyavoidthisphenomenon. tratedinFig.1.Asexpectedfromapriorigeologicalinformation, 19
the borehole paths mainly intersected troctolites containing vari-
able amounts of sulphide in addition to successfully intersecting
2.3.1 Weakanisotropy
the target: an approximately 60-m-thick mineralized zone, corre-
In order to make use of the waveform tomographic approach de- spondingtothelocationoftheEasternDeepsFeederatadepthof
scribed above where the media is anisotropic, a fully anisotropic approximately720m.Afterdrilling,theboreholegeometrieswere
forward and inverse approach would be required. Although such surveyedwiththeuseofanorth-seekinggyroscope,whichrecorded
an approach has been proposed (Alkhalifah 2000; Brossier et al. theazimuthalorientationanddipoftheholesat30mintervals.
2010;Plessix&Perkins2010;Warneretal.2013),wefollowhere Borehole sonic logs, also presented in Fig. 1, give a high res-
theapproachdevelopedandappliedbyPrattetal.(2004),inwhich olution measurement of the vertical seismic wave velocity in the
theisotropicapproachismodifiedinasimplewaythatallowsus immediatevicinityoftheborehole(atfrequenciesontheorderof
to incorporate a 1-D background model consisting of VTI, ellip- 20kHz).Weobserveasignificantandrapiddecreaseinsonicveloc-
tical anisotropy that varies on a layer-by-layer basis. Under these ityassociatedwiththemineralizedtargetzone.Thisisconfirmation
Estimatinganisotropybysimulatedannealing 1591
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/g
ji/a
rtic
le
-a
b
s
tra
c
t/1
9
9
/3
/1
5
8
6
/6
0
9
0
7
8
b
y
Figure 1. An idealized mock-up of the geological cross-section illustrating (schematically) the expected geologyat survey depths [derived from Evans- g
u
Lamswoodetal.(2000)].Themajorityofthecountryrockistroctolite-gabbro,whichsurroundsamineralizedmassivesulphidezone.Thetwoboreholesused e
s
inthesurvey,andasubsetoftheseismicraypathsareoverlain.Alsopresentaretheactualsonicvelocitylogsinthereceiver(left-handpanel)andsource t o
(right-handpanel)holes.Thelowvelocitiesseeninthesoniclogsareassociatedwiththemineralizedzone. n 0
9
A
p
thattheseismicvelocityofthemassivesulphidebodyisinmajor receiver locations. The sources were deployed at 1m intervals in ril 2
contrasttothetroctolitesofthehostrock,anobservationsupported thesourcehole,and121shotswererecordedforeachreceiverlay- 01
9
bylaboratorymeasurements(Eatonetal.2010).Sinceseismicto- out.ThesourcewasbasedontheSweptImpactSeismicTechnique
mography attempts to produce a map of subsurface velocity, this (SIST;Parketal.1996).Thisallowsanacceptablesignal-to-noise
contrastmotivatestheuseofthemethodintheVoisey’sBayarea, ratiotobeachievedwithoutusinghighenergyimpulsivesources,
asweexpecttoseeastrongvelocityanomalyassociatedwiththe which can damage, and sometimes collapse, the host boreholes.
targetzoneinourresults. BoththesourceandthehydrophonesweremanufacturedbyVibro-
metricLtd.
Werefertothecombinationofastaticreceiverarraydepthand
the corresponding shots as a ‘panel’; in total three panels were
3.2 Dataacquisition
shot.Astherewassomeoverlapintheshotlocationafteralayout
Anarrayof29hydrophoneswasusedinthereceiverhole,witha change, the total number of unique shot locations was 236. The
hydrophoneintervalof2m.Thisarraywasdeployedatthreesepa- shotsthusextendedforatotalof236m,beginningapproximately
ratedepths,spaced60mapart.Thisresultedinatotalof87unique 100m above the top of the massive sulphide zone, and ending
1592 M.V.Afanasievetal.
D
o
w
n
lo
a
d
e
d
fro
m
Figure2. SelectedrawreceivergathersfromtheVoisey’sBaydata;(a)[608m]and(b)[670m]arerepresentativeexamplesfromthehighvelocitytroctolites h
(asseeninthesoniclogsinFig.1);(c)[714m]showsareceiverlocatedwithinthelowvelocitymassivesulphidezone.Redlinesonthefiguresdepictthefirst ttps
arrivalpicks. ://a
c
a
d
approximately50mbelowthebottom.Duetothedipofthetarget AlsoevidentinFig.2(especiallyinFig.2c)arehighamplitude e
m
zone, the first receiver was also approximately 100m above the arrivals which appear to have a linear, steep slope representing ic
.o
topofthemassivesulphidezone.Unfortunately,duetothelimited wavespropagatinginthesourceboreholewithalowandconstant u
p
depthofthereceiverborehole,thefinalreceiverwaslocatedonly propagation velocity. These are caused by the presence of tube .c
o
afewmetresbelowthetargetzone,limitingangularcoveragefor waves, which are waves propagating vertically within the water- m
/g
the deepest part of the target zone. The hole-to-hole offset at the filled boreholes. These waves are the result of body wave energy ji/a
sauprpvreoyxidmeapttehlyw2a6smaptpowroaxridmsattheelybo3t0tomm,oalfththoeughholeths.is narrowed to bpreoinpgaginatjeecwteidthinatoctohneswtaanttervceololucimtynoofftahpepsrooxuirmceahteollye,1w5h0e0rmetsh−e1y. rticle
Whenthesetubewavesencountergeologicaldiscontinuities,such -ab
s
asthoseattheedgesofthemassivesulphidezone,theirenergyis tra
3.3 Datadescription convertedtobodywaves,whicheventuallyreachthereceiverhole. ct/1
Asourwaveformtomographyroutinedoesnotmodeltubewaves, 9
9
Thesignalqualityfromthecrossholesurveywasingeneralgood, their presence in the data can lead to artifacts in the final image. /3
/1
withclearfirstarrivalwaveformsandalowlevelofnoise.Asmall However,theirremovalischallenging(asdescribedbelow). 5
8
6
number of receiver gathers exhibited corruption, and these were /6
identifiedandremovedfromthedatasetwithoutfurtherprocessing. 09
0
ExamplesofsinglereceivergathersfromtheVoisey’sBaysurvey, 7
eachshowingthereceivedwaveformsfrom121uniquesourcepo- 3.4 Traveltimetomography 8 b
y
sitions,aregiveninFig.2.Forthesedata,thefirstarrivaltimeswere Inordertogenerateanappropriatestartingmodelforwaveformto- g
u
pickedbyhandusingProMAX.WherepresentinFig.2,firstarrival mography,traveltimetomographywasappliedtothedatafollowing es
timesarehighlightedinred.Notalltraceswerepicked;wherepick- traveltimepicking.Asdescribedpreviously,ourtraveltimetomog- t o
n
ingwasdifficultduetonoiselevelsnoarrivaltimewasassigned. raphyroutineisbasedontheworkofChapman&Pratt(1992)and 0
9
Anabsenceofanarrivaltimepickresultedintherejectionofthat Pratt&Chapman(1992).Inadditionto2-Dimagesofthehorizon- A
p
p(ba)rtsihcuowlarrtercaecieveinrssluobcsaetqeduewntitwhianvethfoerhmigthomveolgorcaiptyhytr.oFcitgosli2te(az)oanned. toaflTPh-womavseenv’esloecpistyil,otnhe(εt)raavnedltimdeeltainv(δe)r,siaolnosngincwluitdhee2s-tDimiamteasgoesf ril 20
1
Here,thewaveformswerewellbehaved,andthefirstarrivaltimes theorientationofsymmetryaxesundertheassumptionof2-DTTI 9
wererelativelyeasytopickonthewaveformdata.Fig.2(c)shows media. The initial results of traveltime tomography, consisting of
areceivergatherwithinthelowvelocityzone.Thewaveformshere imagesoftheisotropicvelocityandthemagnitudeofε,aregiven
were significantly more complicated than those in the troctolite in Fig. 3. These results were produced via five outer (nonlinear)
zone,duemainlytotheinfluenceofthestrongvelocitycontrastat iterations, between each of which new raypaths were calculated
thetroctolite-sulphideboundary,andarrivaltimepickingwasquite bytheray-bendingmethodofUm&Thurber(1987).Withineach
difficult. nonlinear iteration,the raypaths were keptconstant, and the trav-
A‘precursor’componentofthewaveformsmaybeobservedin eltime residuals were minimized via the LSQR algorithm (Paige
Fig.2,visibleasoscillationsappearingpriortothefirstarrivalpicks. &Saunders1982).Theuseofatiltedsymmetryaxis,insteadofa
Theseprecursors,whichlikelyarisefromanacausallowcutfilter verticalaxis(asisrequiredbyourwaveformtomographyroutine),
applied either during or shortly after the data acquisition, made wasnecessary.Assumingaverticalaxisofsymmetryledtodegra-
pickingevenmoredifficult.Thefilterappearstohavebeenusedto dations in the plausibility of the model; in particular resulting in
removemostoftheenergyinthewaveformsbelow1000Hz. horizontalvelocitiesattheedgesofthemineralizedzonethatwere
Estimatinganisotropybysimulatedannealing 1593
3.5.1 f-kdomainandattemptstoremovetubewaves
Sincetubewavesappearatanearconstantslope,itiscommonto
removethemfromthedatathroughfrequency-wavenumberfilter-
ing(f–kfiltering)(Chenetal.1990).Unfortunately,thehigh-pass
filteringatacquisitionresultedinthetubewavesbeingaliasedin
theusablefrequencyband(above1000Hz).Duetothisaliasing,the
tubewaveenergycouldnotberemovedcleanlyfromthedataus-
ingsimplef–kfiltering,andthemethodwasnotapplied.Wemight
therefore expect this to result in artefacts in our final waveform
tomographyimageswithoutfurtherprocessing.
3.5.2 Datapre-conditioning:windowing,bandpassfiltering,
D
editingandnormalization o
w
n
The waveforms were pre-conditioned for waveform tomography lo
a
via high-cut filtering to less than 2000Hz (to avoid the necessity d
e
oedfimtinogdeblylinregmhoigvhal-forfeqanuoenmcayloduastas)o,uwrcinedaonwdi/nogrrienctehiveetrimgaethdeormsaainnd, d from
normalizationofthedatasothateachsourcegather,andeachre- h
ceivergathercontainconsistentamplitudes.Thewindowingserved ttps
tomuteanynoisepresentbeforethefirstarrivals,andalsotoex- ://a
c
clude(asmuchaspossible)late-arrivingtubewaves,aswellasany a
d
shear waves in the data (which are not modelled by the acoustic em
waveequation).Thewindowbegan0.5msbeforethefirstarrival ic
.o
picks,andextendedto3msafter. u
Figure 3. Traveltime velocity and anisotropy models. (a) Horizontal Imagesofselectedfrequencydomaindatafilesforseveralfre- p.c
velocity;(b)Thomsen’sε. o
quencies,afterbandpassfiltering,timewindowing,editingandnor- m
/g
malization are given in Fig. 4. In these images, sources increase ji/a
inplaceslessthan2kms−1.Theconstraintappliedtoperturbations f1rotom817toal2o3n6gatlhoengvtehretichaolriazxoinst.aTlahxreise,adnisdtirneccteibvleorcskisncarreeasveisfirbolme, rticle
toparameterscontrollinganisotropywasthesameasthatapplied -a
and each block corresponds to one survey panel. Note that some b
to perturbations to isotropic velocity. A heavier constraint on the sourcesarepresentinmorethanonepanel,illustratingtherepetition stra
parameters controlling anisotropy led to an increase in traveltime c
residuals,alongwithvelocitymodelswhichweregeologicallyun- osofusrocme–erescoeuirvceerppaoisristiwonitsh.tDhueecltoosethsteogfefsoemtseatrpypeoafrtahpepsruorxvimeya,ttehlye t/19
9
realistic.Thermstraveltimeresidualafterfivenon-lineariterations along the main diagonal in the images, and source–receiver pairs /3
was0.049ms.Anisotropyinthehighvelocitytroctolites(Fig.3b; /1
withapproximatelyequaloffsetsappearalongsubdiagonalsparallel 5
ε=4–5percent)coincideswellwiththepreviousstudiescitedin 86
tothemaindiagonal. /6
theintroduction(i.e.Cosma&Enescu2003).Aroundtheedgesof While time-windowing of the data introduces a low frequency 09
thelowvelocityzone,astrangelobe-likestructureisobservedin 0
component due to spectral leakage, it is obvious from Fig. 4 that 7
theanisotropymodel.Webelievethisisaninversionartefact,anda mostoftheselowfrequencydataarecorruptedandnoisy.Whilethe 8 b
discussionofitsoriginwillbereturnedtoinSection5.3.1. waveformsat500Hzshowsomecoherency(Fig.4a),thosebetween y g
u
We should also note that although we assume a transversely 600and1000Hzareveryincoherent(arepresentativeexampleat es
isotropicmodel(withanarbitraryaxisofsymmetry)intraveltime 800Hz is givenin Fig.4b).Thisechoes ourobservation thatthe t o
n
tomography,thisislikelyasimplificationofreality.Theborehole datawerefilteredatacquisitiontoremovedatabelow1000Hz.We 0
9
corelogs(Inco2000)reportonhighlyfracturedrock,withmultiple expect the lack of low frequency data to increase the chances of A
fArascwtuerlel,spyomstmemetprylaacxeems,eonvtesreprprienntetidneonmleotcaamloflropwhiasnmd(sIhnecaor2f0ab0r0i)ciss. cycle-skipping. pril 2
Above1000Hz,theprocesseddataareofgoodqualityandco- 0
knowntobehighlyanisotropic(Kern1993),andmayaddadditional herent.Representativeexamplesat1300and1600Hzcanbeseen 19
anisotropicpropertiesalongsomecracks.Tocompletelycharacter-
inFigs4(c)and(d),respectively.
ize the anisotropy would likely require the use of orthorhombic,
OtherrelevantinformationcanbevisuallyextractedfromFig.4.
monoclinicoreventriclinicsymmetrysystems.
Over most of the survey (i.e. in the high velocity troctolite host
rocks),thecoherentpatternsinthesedisplaysaresubparalleltothe
maindiagonal,andcycleinphaseperpendiculartothemaindiag-
onal.However,forthosesourcesandreceiversinthelowvelocity
3.5 Datapreprocessing
mineralizedzone(sources130–180,receivers58–87),weseethat
Prior to waveform tomography, waveform data must be prepro- thesecoherentpatternsinthedatachangeinslope.Withinthiszone
cessedinanattempttoattenuateundesirablesignals.Thesesignals thecoherencyisnolongerparalleltothethemaindiagonal.Thisis
include records of tube waves, ambient noise, and elastic effects. asignthatthedatainthisregionarespatiallyaliased,whichstrongly
Duringthepreprocessingstage,itwasobservedthatthedatahad suggestsacorrespondingdecreaseinseismicvelocity(shortening
beenhigh-passfilteredtoabove1000Hzatacquisition. thewavelengthstobelowtheNyquistcriterion).
1594 M.V.Afanasievetal.
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/g
ji/a
rtic
le
-a
b
s
tra
c
Figure4. FrequencydomaindisplaysoftheVoisey’sBaydatafollowingbandpassfiltering,timewindowingandamplitudenormalizationoverthesourceand t/1
9
receivergathers.(a)500Hz;(b)800Hz;(c)1300Hz;(d)1600Hz.Weseeampleevidenceofcoherentpatternsparalleltothemaindiagonal,especiallyon 9
/3
panels(c)and(d).Asthedistancesbetweensourcesandreceiversareincreased(i.e.perpendiculartothediagonals)thesignalmaybeobservedtocyclein /1
phase[evidentasachangefromredtobluecausedbythechangeinthepolarityofRe(u)]. 58
6
/6
Discontinuitiesparalleltothereceiveraxis(theverticalaxis)are (ε=δ),andthatthesymmetryaxesareclosetovertical.Examin- 09
0
alsoapparentatseverallocations(forexample,aroundsource75). ingFig.3(b),andnotingthattheimagewascreatedbyallowingfor 7
8
Asthehorizontalaxisrepresentssourcelocation,thesediscontinu- (ε (cid:9)= δ) and allowing a tilted symmetry axis, we see that these b
y
itiesmayberelatedtosource-boreholecouplingissues,aswellasto assumptionsareunrealistic.Nevertheless,arealistictraveltimeto- g
u
possiblesourcemis-locationissues.Finally,whitezonesparallelto mographyvelocityimagecouldonlybefoundbyassumingtilted es
thesourceaxisarealsopresent.Theserepresentdeletedreceivers, transverseisotropy.WethereforeinterpretFig.5asameasureofthe t o
n
whichwereexcludedduetoexcessivenoise.Specifically,the23rd averagemagnitudeofεateachdepthlevel—regardlessofsymme- 0
9
receiver in the array was consistently corrupted, and is therefore trydirection.Thevelocitymodelwasinterpolatedontoa96×500 A
p
absentfromeachpanelonFig.4. cellfinitedifferencegrid,withthecelledgelengthsetat0.5m. ril 2
Timedomainwaveformsweregeneratedthroughthisanisotropic 0
1
velocitymodel.Themainpurposeofperformingtime-domainmod- 9
elling was to determine if the half-cycle criterion was satisfied.
3.6 Assessment—waveformmodellingoftraveltime
Theresultingpredictedreceivergathersfromthetraveltimemodel
tomographyresult
(Fig.6b)showamismatchinthefirstarrivaltimesatlargeroffsets
Inordertoassessthequalityofthetraveltimetomographyresults ofmorethanahalfcycle.Specifically,weseethisforshotsdeeper
forthepurposesofwaveformtomography,wesimulatedwaveform than730m.Atnearoffsets,thefirstarrivalswereconsistentlypre-
data numerically in the traveltime velocity model (Fig. 3a) and dicted to within one half wavelength of the dominant frequency.
comparedthemwiththereal(pre-processed)data.Topreparethe Thepatternofagoodfitatnearoffsets,andofcycle-skippingat
modelforwaveformmodelling,wecreateda1-D,verticallyvary- higheroffsets,wasobservedeverywherewithinornearthelowve-
ing anisotropy profile by averaging the ε model in Fig. 3(b) over locityzone.Sinceanincreaseinoffsetcorrespondstoanincrease
horizontallayers.TheresultantprofilecanbeseeninFig.5.Note in angle away from the horizontal, this strongly suggests that the
thatthisapproachtowaveformtomographyassumesthatthe2-D averaged,elliptical,VTIanisotropymodelfoundthroughtraveltime
anisotropy varies only with depth, that the anisotropy is elliptical tomographyisinsufficientforuseinwaveformtomography.
Estimatinganisotropybysimulatedannealing 1595
ofsimulatedannealingtotheanisotropyproblem,weconsiderthis
sufficient.Theanisotropywasellipticinnature(ε=δ).
The true models for V and ε are given in Fig. 7. The source-
p
receiver geometry was the same as that for the field data, with
boreholesapproximately30mapart,sourcesspaced1mapart,and
receivers spaced 2m apart. The background velocity was set to
be7.5kms−1.Adipping,lowvelocitylayerwasintroduced,sim-
ulatingtheactuallocationandorientationofthemassivesulphide
zone, as estimated from the borehole sonic logs (see Fig. 1). As
well,alayerofslightlyhighervelocitywasplacedjustabovethe
lowvelocitylayer,correspondingtotheslightincreaseinborehole
sonic velocities seen above the massive sulphide zone. Random
heterogeneitieswerethenaddedtothemodelusingthemethodsof
Ikelleetal.(1993)andKameietal.(2005);themagnitudesofthese
heterogeneities were estimated from observed fluctuations in the D
o
boreholelogs.Inthehighvelocityzone,theheterogeneitieshada w
standarddeviationof300ms−1,whileinthelowvelocitylayerthe nlo
heterogeneitieshadastandarddeviationof210ms−1.Averticalto ad
e
Figure5. 1-Danisotropyprofileobtainedbyaveragingthemagnitudeofε hhoetreizroognetanleriatiteios,ftooritnhveecsotirgraetleattihoenilmenpgatchtsooffh1o:4riwzoanstaulseladyfeorrintgheosne d from
inthe2-DmodelshowninFig.3(b)overhorizontallayers. theinvertedanisotropymodels.Thelongaxisoftheheterogeneities h
3.7 Necessityofasemi-globalapproach added to the low velocity zone was rotated to coincide with the ttps
dipofthezoneitself.Densityvaluesforthehighandlowvelocity ://a
Ausableanisotropymodelcouldnotbefoundthroughmanualpick- zoneswerechosentocorrespondtotheaveragedensityoftrocto- ca
d
ing and traveltime tomography. Instead, generating such a model lite(2.7gcc−1)andpyrrhotite(4.6gcc−1),respectively.TheP-to em
required a novel method; such a method would have to respect S-wavevelocityratiowasseteverywhereequalto1.87.Inthehigh ic
thewaveformnatureofthedata.Sincecycle-skippingwaspresent, velocityzone,εwassettoaconstant0.15,whileinthelowvelocity .ou
p
it was clear that a local inverse method for waveform tomogra- layeritdroppedto0.05. .c
o
phywouldnotbesufficient.Oursolutiontothiswastodevelopa Forbothmodels,anexplosive,Ku¨pper-waveletsourcewasmod- m
waveform-basedsimulatedannealingapproachthatwouldfindthe elled,withapeakfrequencyof666Hzandamaximumfrequency /gji/a
apphpy.roTpordiaetteeramniisnoetrthoepyvimaboidlietyl ofofrthseubsismequulaetnetdwanavneefaolirnmgatopmproogarcah- opfonaepnptrooxfitmheatwelayv2efi00el0dHwza,sanexdtrthacetescdaalatrecauchrl-rferceeeiv(perrepsosusirteio)nc.om- rticle
asaglobalmethodforanisotropyinversions,atestofthemethod -a
b
s
wasinitiallyperformedonsyntheticallygeneratedcrossholedata, tra
forwhichthetruemodelswouldbeavailableforcomparisonpur- 4.1 Traveltimetomography ct/1
poses.SyntheticmodelsbasedontheVoisey’sBaygeologicalset- 9
ting were created, and seismic data were calculated numerically Topreparethesyntheticdatafortraveltimetomography,littleeffort 9/3
withinthesemodels,fullyrespectingtheanisotropydistributions. wasrequired.Duetothenoise-freenatureofthesyntheticdata,the /15
Thenextsectionfollowsthedevelopmentofsimulatedannealing majorityofthefirstbreakpickingwasveryeasy,andcouldmostly 86
strategiesandmethodologieswithreferencetothesyntheticexam- becarriedoutsemi-automatically.Thisbeingsaid,therewerestill /60
9
ples,andwewillreturntotherealVoisey’sBaydatainSection. significantdifficultiesattheedgesofthelowvelocityzone,where 07
somemanualinterventionwasrequired. 8 b
Traveltimetomographywasperformedfollowingthesamephi- y g
4 SIMULATED losophygiveninSection3.4.Fivenonlinearinversionswerecarried ue
s
ANNEALING—DEVELOPMENT AND out,duringwhichtheconstraintonmodelroughnesswasgradually t o
APPLICATION TO SYNTHETIC relaxed.Traveltimevelocityandanisotropyresultsfortheelliptical n 0
CROSSHOLE DATA synthetic data are depicted in Fig. 8. The values of ε in the high 9 A
velocityhostrockzoneareoverestimatedbyapproximately5per p
The synthetic velocity and anisotropy models we use in this sec- cent.Lobesofstronganisotropyaroundthelowvelocityzone,sim- ril 2
tiontotesttheeffectivenessofsimulatedannealingweredesignedto 0
ilar to those observed in Fig. 3, are also present. The failures of 1
simulatetheessentialelementsoftheactualVoisey’sBayfielddata. 9
thetraveltimetomographymethodisattributedtothefinescaleof
Inordertodothis,wesimulatedtheforwarddatawithafullelastic,
theheterogeneitiesthatlieoutsidetheresolutioncapabilityofray
anisotropicmethodyieldingrealisticsyntheticdatathatincludeboth
tomography;adiscussionofthiseffectfollowsinSection5.3.
qPandqSwaves.ElastictimedomainVTIfinitedifferencemod-
ellingwasperformedusingthefreelyavailableopensourcesoftware
packageMadagascar(Yan&Sava2011).Madagascarwaschosen
4.2 Waveformtomography
sothatarealistic,anelliptic2-Danisotropymodelcouldbeused,
and aswell toavoid thepitfalls ofperformingan ‘inversecrime’
4.2.1 Frequencydomainwaveformdata
(Colton&Kress1998;Wirgin2004).Theanisotropyweintroduced
matchedthatobservedinthetraveltimeinversionsabovewiththe Sincethesyntheticdatawerenoisefreeandnotcorruptedbytube
realdata,withoneimportantdistinction:thesyntheticanisotropy waves,itwasnotnecessarytoapplytimedampingaswasdonewith
modelsassumedaverticalaxisofsymmetry(VTI).Sincethepri- therealdata,andthusthedampingparameterτ ineq.(10)wasset
marygoaloftheseexperimentswastoinvestigatetheapplicability to infinity. This results in eq. (10) reducing to the classic Fourier