Geophysical Journal InternationalGeophys. J. Int. (2013) doi: 10.1093/gji/ggt305

GJI

Mar

ine

geos

cien

ces

and

appl

ied

geop

hysics

Acoustic full waveform tomography in the presence of attenuation:a sensitivity analysis

A. Kurzmann,1 A. Przebindowska,1 D. Kohn2 and T. Bohlen1

1Karlsruhe Institute of Technology, Geophysical Institute, Hertzstr. 16, 76187 Karlsruhe, Germany. E-mail: andre.kurzmann@kit.edu2Department of Geophysics, Institute of Geosciences, Kiel University, Otto-Hahn-Platz 1, 24118 Kiel, Germany

Accepted 2013 July 26. Received 2013 July 25; in original form 2012 December 18

S U M M A R YFull waveform tomography (FWT) is a powerful velocity building method to exploit thefull richness of seismic waveforms in complex media. Most applications today neglect thefrequency-dependent amplitude decrease and phase velocity dispersion caused by intrinsicattenuation. In this study, we present a numerical investigation of the influence of attenuationon the recovered velocity model. Based on the generalized standard linear solid as rheologicalmodel, we incorporate attenuation into a 2-D time-domain acoustic FWT scheme. Attenuationis considered here as a modelling and not an inversion parameter. We investigate two reflectionseismic experiments: (1) a visco-acoustic 1-D model and (2) a visco-acoustic version ofthe Marmousi model to which realistic quality factors are assigned. In the presence of softrocks with pronounced absorption we observe a poor recovery of the velocity model whenattenuation effects are not taken into account in the modelling. By considering an appropriateattenuation model in the forward modelling of the FWT, the accuracy of the reconstructedvelocity model improves significantly in both cases. Even a homogeneous background qualityfactor model might allow a satisfactory recovery of the velocity model, provided that it is aquite good representation of the shallow structures. Our results suggest to consider attenuationas a smooth background modelling parameter in reflection seismic configurations to improvevelocity model building by a purely acoustic inversion scheme.

Key words: Inverse theory; Numerical approximations and analysis; Seismic attenuation;Seismic tomography; Wave propagation.

1 I N T RO D U C T I O N

The aim of full waveform tomography (FWT) is to find a sub-surface model which explains the recorded seismic data, that is, ititeratively minimizes the difference between observed and syntheticseismograms. In the majority of FWT applications, we face a mul-tiparameter problem. In particular, attenuation might significantlyaffect amplitudes and phases of seismic signals due to frequency-dependent amplitude decrease and phase velocities (Causse et al.1999).

Nowadays, acoustic FWT is usually applied to recover the P-wavevelocity model in transmission (e.g. Wang & Rao 2006; Maurer et al.2009) and reflection seismic configurations (e.g. Hu et al. 2009;Bleibinhaus & Hilberg 2012). Most applications, thus neglect theimpact of attenuation on seismic waveforms. The purpose of thiswork is to quantitatively investigate the validity of an acoustic FWTin the presence of attenuation.

There are two different ways to take attenuation into account.Attenuation may be used as a passive parameter, that is, as mod-elling parameter only. The aim is to improve the accuracy of theP-wave velocity model (e.g. Brenders & Pratt 2007). Alternatively,

a multiparameter inversion, such as a visco-acoustic FWT (Liao &McMechan 1995, 1996; Askan et al. 2007; Hak & Mulder 2008,2011; Kamei & Pratt 2008) or the asymptotic visco-acoustic wave-form inversion (Ribodetti et al. 2004), involves attenuation as anadditional inversion parameter.

Field data applications that take attenuation into account aremainly conducted in the frequency domain. Hicks & Pratt (2001)obtained reliable quality factor subsurface models from shallowseismic data recorded in the North Sea. Takam Takougang & Calvert(2011) reconstructed realistic velocity models from marine reflec-tion data. They combined different inversion strategies, such as amultistage approach with incremental frequency selection (com-pare Bunks et al. 1995; Sirgue & Pratt 2004) and separate inversionof near and far offsets. The vP-only inversion at low frequenciesconnected with the joint inversion for vP and QP at higher frequen-cies improved the recovery of both vP and QP models. Malinowskiet al. (2011) applied a frequency-domain visco-acoustic FWT towide-aperture seismic field data recorded in the Polish basin. Theyfocus on the applicability of a visco-acoustic joint inversion forP-wave velocity vP and QP. They reconstructed satisfactory sub-surface models for both vP and QP coinciding with the expected

C© The Authors 2013. Published by Oxford University Press on behalf of The Royal Astronomical Society. 1

Geophysical Journal International Advance Access published August 28, 2013

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

2 A. Kurzmann et al.

geology. In this context, Mulder & Hak (2009) discussed problemsof a joint inversion with respect to short-aperture data in reflectionseismics. They found that the ill-posedness of the inverse problemcauses a very poor recovery of both phase velocity and attenuation.A brief overview of recent publications in this area is also given byVirieux & Operto (2009).

In this paper, we use a time-domain implementation to studythe role of attenuation for weakly and strongly attenuative mediain marine environments. We investigate the validity of an acous-tic inversion scheme in presence of attenuation by systematicallyquantifying the errors in the inverted velocity model. Attenuation isincorporated as a passive parameter in forward modelling only. Weapply a purely acoustic FWT including acoustic or visco-acousticmodelling to analyse and compare the impact of attenuation on thereconstruction of the velocity models from visco-acoustic reflectiondata sets.

Although this study concentrates on 2-D acoustic FWT in thetime domain, it is also targeted to 3-D applications. While time-domain modelling (e.g. Bohlen 2002) is commonly used in 3-Dapplications, the high efficiency of 2-D frequency-domain mod-elling including a straightforward implementation of attenuationcontrasts with highly demanding 3-D frequency-domain modellingusing direct solvers (e.g. Wang & de Hoop 2011). However, thereare robust and efficient frequency-domain approaches, such as itera-tive solvers (Plessix 2007). Furthermore, the time-domain approachhas some advantageous features: straightforward time-windowingof data and the consideration of broad frequency bands (instead ofsingle frequencies).

This work comprises implementation of finite-difference visco-acoustic modelling into the time-domain acoustic FWT and its ap-plication to visco-acoustic data. We choose two examples: a simple1-D medium with a reflection acquisition geometry providing a highray coverage and the 2-D Marmousi model with a towed streamergeometry. For both experiments we perform the same inversiontests. We investigate the resulting data misfits and model errors todemonstrate the footprint of attenuation on the recovered velocitymodels.

2 M E T H O D O L O G Y

The following section gives a brief overview on the FWT methodused in our implementation. While the inversion scheme is purelyacoustic, attenuation is incorporated into visco-acoustic modellingby using a quality factor model QP. In principle, the implementa-tion of visco-acoustic modelling relies on the generalized standardlinear solid (Liu et al. 1976) and yields the following visco-acousticwave equation (e.g. Emmerich & Korn 1987; Carcione et al. 1988a;Robertsson et al. 1994):

p (x, t) = κr (x) ∇ · w (x, t) [1 + L τP (x)] +L∑

l=1

rl (x, t) ,

w (x, t) = 1

ρ (x)∇ p (x, t) ,

rl (x, t) = − 1

τp,l[κr τP (x) ∇ · w (x, t)+rl (x, t)] , l ={1, . . . , L},

(1)

where p(x, t) and vector w(x, t) denote pressure and particle veloc-ity field. Attenuation is characterized by L relaxation mechanismsrequiring additional variables rl (x, t) (with l = {1, . . . , L})—the so-

called ‘memory variables’. In contrast to the acoustic wave equation,the acoustic bulk modulus is replaced by the relaxed bulk modulusκ r. Further symbols denote the relaxation parameter (τ P) obtainedfrom the QP model, the relaxation times (τ p, l) and density (ρ). Thetheory of visco-acoustic modelling including the implementation ofperfectly matched layers (PML) is described in detail in Appendix.

The solution of the inverse problem is mainly based on the time-domain FWT introduced by Tarantola (1984) and Mora (1987). Itcomprises the adjoint method and the conjugate gradient method us-ing the least-squares misfit definition. In contrast to the generalizedderivations given by these authors, we assume some simplifications.

Within the scope of this work, we focus on the reconstruction ofthe P-wave velocity model in a visco-acoustic medium. To avoidany unwanted side effects, we use homogeneous density, which isnot subject to the inversion and can be neglected in subsequent con-siderations. Furthermore, there is no inversion for the source timefunction due to the usage of the true source signal. Hence, insteadof parametrization m(x) = [κκκ(x) , ρρρ(x) , q(x, t)]T with bulk mod-ulus κκκ(x), density ρρρ(x) and source term q(x), the inverse problemonly considers m(x) = κκκ(x). In the following considerations, therelaxed bulk modulus κκκ r is treated as bulk modulus κκκ . The acous-tic forward problem is solved by the general function for pressure:p = p(κκκ,ρρρ, q). In general, the observed pressure data pobs are notcompletely explained by the given initial model. The aim of the non-linear inversion is to minimize the difference between synthetic datap (m) and observed data pobs. The definition of data residuals �pis based on the least-squares misfit function E (Crase et al. 1990)and calculated as a sum over all data samples Nt, receivers Nr andsources Ns:

E(m) = 1

2‖�p‖2 = ‖p (m) − pobs‖2 . (2)

The gradient of the misfit function involves the transpose of theFrechet derivatives matrix, DT = (∂p / ∂m)T , which can be derivedby linearization of the forward problem, δp = D δm (Tarantola1984; Mora 1987), in the framework of the Born approximation.The partial derivative of the functional E with respect to m yieldsthe steepest ascent gradient vector g:

g = ∂ E

∂m=

(∂p

∂m

)T

�p = DT �p. (3)

If g = 0, the minimization problem is solved. An iterative solutionat iteration k is given by

mk =mk−1−μk Hk

{DT

k [p(mk−1) − pobs]}=mk − μk Hk gk, (4)

with the inverse Hessian matrix Hk and the step length μk of thegradient algorithm. The computation of Hk is highly demanding(Pratt 1998). However, there are matrix-free methods, such as thel-BFGS quasi-Newton method, which computes in a finite-difference sense the product of the inverse of the Hessian withthe gradient from few gradients and few solutions of previous iter-ations (Nocedal & Wright 1999), or the truncated Newton method(Metivier et al. 2012). Using a gradient method, the inverse Hessianmatrix can be approximated by Hk ≈ I. The resulting model updateis given by

mk = mk−1 − μk I gk . (5)

The steepest ascent gradient can be written in terms of model per-turbations δm (Mora 1987). Due to the focus on inversion for bulkmodulus (or P-wave velocity, respectively), we concentrate on thecorresponding gradient gκκκ . It describes the perturbation of bulk

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

Acoustic FWT in presence of attenuation 3

modulus given by a cross-correlation of forward and residual wave-fields summed over all Ns sources (Tarantola 1984):

gκκκ = 1

κ2(x)

∑Ns

∫tdt p

[∑Nr

G ∗ δ p

]= 1

κ2(x)

∑Ns

∫tdt p p′, (6)

with the source wavefield p := p(xs, x, t), the residual wavefieldp′ := p′(xs, xr, x, t) and the Green’s function G := G(x, −t ; xr, 0).While p is forward-propagated from sources xs to receiver locationsxr, p′ involves the superposition of all data residuals δ p over Nr

receivers and is back-propagated from receivers to sources. Apartfrom the parametrization m = κκκ , it is desirable to have the P-wavevelocity vP instead of the bulk modulusκκκ . On the assumption that themedium is parametrized by one parameter, the gradient is definedin terms of the new parameter vP by:

gvP= ∂κκκ

∂vPgκκκ , (7)

involving the Jacobian (Mora 1987):

JvP := ∂κ

∂vP= 2 ρ vP. (8)

At iteration k, the P-wave velocity is updated by

vP|k = vP|k−1 − μk gvP|k . (9)

Furthermore, the preconditioned conjugate gradient method (e.g.Polak & Ribiere 1969; Luenberger 1984; Nocedal & Wright 1999)improves the convergence of the steepest descent gradient algo-rithm. This required preconditioned gradient is obtained by applica-tion of a preconditioning operator Pk to the steepest ascent gradient.In our work, Pk is used to weight gk , that is, it damps source artefactsor excludes the model update at predefined locations. In detail, weapply a wavefield-based preconditioning of the gradient [modifi-cation of the method proposed by Igel et al. (1996) and Fichtneret al. (2009)]. The operator is computed from the source wavefieldsp(x, t) and residual wavefields p′(x, t):

Pk(x) = b(x)

maxx b(x), (10)

with

b(x)= 1

a(x)+Cstab aand a(x)=max

t|p(x, t)|+max

t|p′(x, t)|.

The auxiliary parameter a denotes the spatial average of a(x). Thecoefficient Cstab stabilizes the computation of Pk and thus scalesits strength and dynamic range. However, Cstab has to be estimatedmanually, for example, Cstab = 0.001. The application of Pk is anelement-wise operation at each spatial location. As a similar ap-proach, one can use the diagonal terms of the approximate Hessianmatrix (Shin et al. 2001), which is computed from the source wave-field only and also takes effects caused by geometrical spreadinginto account.

At each iteration the estimation of the optimal step length μk isperformed. It relies on Pica (1990) and is composed of two parts:a user-defined relative factor μrel|k and a factor used to scale thegradient to the maximum range of the model parameter. This allowsa meaningful physical unit and a proper distance of the gradient.Eq. (9) is rewritten as

vP|k = vP|k−1 − μrel|kmax

(vP|k

)max

(∣∣gvP|k∣∣) gvP|k . (11)

The relative factor μrel|k is approximated by a parabolic curve fittingmethod (compare the line-search method with quadratic or cubicinterpolation in Nocedal & Wright 1999). An initial step lengthμrel,ini has to be provided at iteration k = 1. Two additional valuesare computed by applying a constant factor a > 1: μrel,low = μrel,ini

aand μrel,high = a μrel,ini. The minimum of the parabola, fitted to thecorresponding data misfits [E(μrel,ini), E(μrel,low), E(μrel,high)], in-dicates the optimal step length at iteration k. The adaptive estimationreveals a significant and stable reduction of the data misfit function(e.g. Kurzmann et al. 2008).

The following paragraph summarizes the involvement of modelsin our inversion scheme. At the first iteration, the method requiresan initial model κκκ0. It is obtained from the user-defined initial refer-ence model vP,ref |init by applying the model relaxation (A5). On theone hand, the relaxed bulk modulus is required for visco-acousticmodelling and, on the other hand, it is treated as an acoustic parame-ter by the inversion algorithm (eq. 9). Due to the different meaningsof relaxed model (defined at frequency ω → 0) and reference model(defined at reference frequency ω0, see Appendix), the desired finaloutput has to be a reference velocity model to ensure the physi-cal consistency with the initial reference velocity model. First, weobtain the updated velocity model vP|result from eq. (9). Then, wecompute the resulting reference model vP,ref |result by revoking themodel relaxation (A5), which can be rewritten in terms of P-wavevelocity:

vP,ref |result = vP|result

√√√√1 +L∑

l=1

ω20/ω

2r,l

1 + ω20/ω

2r,l

τττ P. (12)

3 S Y N T H E T I C WAV E F O R MT O M O G R A P H Y E X P E R I M E N T S

We investigate the impact of attenuation in applications of acousticFWT to visco-acoustic data with and without a passive QP model.Several acoustic full waveform inversion tests are employed to anal-yse the footprint of attenuation on vP inversion results. We analysethe effects of both visco-acoustic and acoustic modelling in acousticFWT. The tests are applied to both a simple 1-D medium and themore complex Marmousi model using different initial vP and pas-sive QP models. A list of all tests can be found in Table 1. To avoidunwanted side effects, some general restrictions and only necessary

Table 1. List of all FWT tests for both the 1-D and the Marmousi experiment.

FWT Figures of vP results Data Modelling Initial Passivetest 1-D model Marmousi in FWT vP model QP model

Reference 6(a) 10(a) Acoustic Acoustic Smooth –Test 1 6(b) 10(b) Visco-acoustic Acoustic True –Test 2 6(c) 10(c) Visco-acoustic Acoustic Smooth –Test 3 6(d) 10(d) Visco-acoustic Visco-acoustic Smooth TrueTest 4 6(e) 10(e) Visco-acoustic Visco-acoustic Smooth SmoothTest 5 6(f) 10(f) Visco-acoustic Visco-acoustic Smooth Homogeneous

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

4 A. Kurzmann et al.

preconditioning features are applied. The FWT tests comprise thefollowing set-up and constraints:

(1) the P-wave velocity vP is the only inversion parameter,(2) apart from the gradient preconditioning mentioned above,

we suppress the model update within the water layer due to theassumption of known model parameters,

(3) Marmousi experiment: low-pass filtering of observed dataand source wavelet over multiple stages is applied to improve theprogress of the inversion (compare Bunks et al. 1995; Sirgue &Pratt 2004),

(4) stopping criterion for 1-D experiment to obtain most optimalresults: FWT is unable to reduce data misfit (relative thresholdvalue between two successive iterations is 0.0001 per cent), or it isnot possible to compute a meaningful step length,

(5) stopping criterion for shifting within multiple stages in theMarmousi experiment: due to high computational efforts, the rel-ative threshold value of the data misfit between two successiveiterations is 1 per cent.

The reference computation applies acoustic inversion to acousticdata. It aims to show the performance of FWT for a given geologyand geometry of both experiments. Tests 1 and 2 analyse the effectof attenuation on an acoustic inversion with acoustic modelling.Tests 3, 4 and 5 involve visco-acoustic modelling in acoustic FWTand investigate three different passive QP models in conjunctionwith the more realistic initial vP model (compare Table 1).

We calculate data misfits and model errors ε to quantify the per-formance of all tests. They are computed with respect to the truemodel vP|true and to the observed data pobs. Due to the usage of theleast-squares norm in the FWT algorithm, the data misfits are ex-pressed as normalized squared L2 norms, whereas the model errorsare normalized L1 norms (to ensure comparability with deviationfigures in case of the Marmousi experiment):

ε(vP|init) = ‖vP|init − vP|true‖1

‖vP|true‖1(initial model error), (13a)

ε(vP|result) = ‖vP|result − vP|true‖1

‖vP|true‖1(final model error), (13b)

ε(pinit) = ‖pinit − pobs‖22

‖pobs‖22

(initial data misfit), (13c)

ε(presult) = ‖presult − pobs‖22

‖pobs‖22

(final data misfit), (13d)

where pinit and presult denote the synthetic data for initial model vP|init

and final model vP|result, respectively.

3.1 Synthetic experiment: layered 1-D model

The first experiment uses a 2-D model with a 1-D geology (here-inafter referred to as ‘1-D model’) consisting of four layers overa half-space: a water layer on top, followed by highly and weaklyattenuative sedimentary rocks. The corresponding vP and QP mod-els are shown in Fig. 1. Due to the occurrence of a thick layer withhigh attenuation on top of the sediments, this model represents ashallow marine environment. The total model size is 130 × 308 mwith a spatial discretization of �h = 0.5 m. However, due to the ap-plication of a PML (width = 15 m) in finite-difference modelling,all model-related figures are limited to the relevant area (excludingthe PML layer). The acquisition geometry is located at the water

Figure 1. 1-D experiment: Vertical cross-sections of the 1-D models: (a)True and initial velocity model. The range of phase velocity dispersion dueto attenuation is illustrated by shaded areas. The layers are labelled withRoman numbers (water layer and half-space are denoted by ‘I’ and ‘V’,respectively). (b) shows true as well as smooth and homogeneous passiveQP models (‘homogeneous’ with respect to the sub-seafloor area).

Figure 2. 1-D experiment: Approximation of the quality factor. (a) showsthe normalized amplitude spectrum of all observed visco-acoustic data.(b) and (c) illustrate the quality factor approximation and phase velocitydispersion for the second layer using the relaxation frequencies for QP, 0 = 10and the average QP, 0 = 74. The shaded areas represent the bandwidth withrespect to the observed data. The solid and dashed lines show Q(f) andcorresponding phase velocity dispersion for the same relaxation frequenciesfr, l = (1.202, 17.62, 179.4) Hz but true and average QP, respectively.

surface and consists of 24 explosive sources with a spacing of 12 mas well as 278 hydrophones with a spacing of 1 m. For each shotgather all receivers are used. The resulting offsets range from 0.5to 292.5 m. The source signal is a Ricker wavelet with a peak fre-quency fpeak = 80 Hz and the recording time of synthetic seismicdata is 0.21 s.

For visco-acoustic modelling, we determine the relaxation pa-rameters such that we get an optimal representation of constant QP

within the bandwidth of the seismic data (see Fig. 2b). Here, wedefine the bandwidth as the contiguous frequency range, whereinthe decay of the amplitude spectrum with respect to the maximumamplitude is less than 10 dB (shaded areas in Fig. 2). The rangeof phase velocity dispersion of all layers is estimated within this

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

Acoustic FWT in presence of attenuation 5

Figure 3. 1-D experiment: (a) shows the difference between acoustic and visco-acoustic data for the central shot located at x = 160 m. It exhibits trueamplitudes which are clipped to ±4 per cent of the maximum visco-acoustic amplitude. (b) illustrates traces at exemplary offsets. For a better visualization,they are normalized independently for each offset and the direct wave of the zero-offset trace is clipped. Acoustic and visco-acoustic amplitudes are stillcomparable.

bandwidth (see exemplary dispersion for the second layer withQP, 0 = 10 in Fig. 2c) and visualized by shaded areas in verticalsections across the 1-D medium (Fig. 1a). The aim is to analyseif the recovered vP model can be explained by the minimum andmaximum velocity dispersion.

The approximation of relaxation parameters is based on the refer-ence frequency f0 := fpeak and an average QP, 0 = 74 computed fromthe true quality factor model within the sub-seafloor area (Fig. 1b).The acoustic velocity model vP, ref (Fig. 1a) is defined at f0, that is,no dispersion occurs at f0. The bandwidth of the seismic data is lim-ited to � f = [19.6, 141] Hz. This corresponds to a dynamic rangeof 2.8 octaves. Based on the rule of one relaxation mechanism peroctave (Blanch et al. 1995), we use three relaxation mechanisms.The resulting optimal set of relaxation frequencies is fr, l = (1.202,17.62, 179.4) Hz. We found that it is not necessary to obtain relax-ation frequencies for all quality factors of the true model, that is,it is not essential to provide models containing relaxation frequen-cies. Using the relaxation parameters computed from the averageQP, 0 = 74, we obtain quite accurate approximations for all qualityfactors given in the 1-D model. An exemplary quality factor ap-proximation is shown for the second layer with QP = 10 (Fig. 2b).Obviously, for this example the deviation in corresponding phasevelocity dispersion is negligible (compare Fig. 2c).

The L1-based QP approximation error is given with respectto constant QP, 0 and is quantified within the seismic bandwidth.

For all layers, we obtain acceptable approximation errors ofabout 3 per cent. While the approximation at reference frequencyis perfect, the largest errors can be observed at the upper endof the bandwidth. These errors are nearly identical to those ofthe accurate QP approximation using exact QP values of eachlayer.

The smooth initial vP and passive smooth QP models for wave-form tomography are generated by the application of a 2-D Gaussianfilter (size 100 × 100 m, σ = 33) to the sub-seafloor area of thetrue model (see Figs 1a and b). In addition, Fig. 1(a) illustrates theminimum and maximum velocity dispersion with respect to the truevP and QP models.

A comparison of acoustic and visco-acoustic data computed fromthe true model at the central shot location is shown in Fig. 3. Dueto very weak attenuation in water, the direct waves are nearly iden-tical. At near offsets we can observe a quite good match in phasesand amplitudes of the seafloor reflection. However, at larger off-sets and for all later reflection events, there are significant differ-ences between acoustic and visco-acoustic data. The phase mis-fit is explained by the highly dispersive properties of the secondlayer.

In the following two paragraphs, we give a more detailedoverview of the impact of attenuation on the seismic data by inves-tigating both the sensitivity on model perturbations and the cycle-skipping phenomenon.

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

6 A. Kurzmann et al.

3.1.1 Sensitivity analysis

Although attenuation is not subject to the inversion and density isneglected in the synthetic experiments, we investigate the trade-offbetween model perturbations in the visco-acoustic parametrizationm = (vP, QP, ρ). Therefore, we compute the partial derivative pres-sure wavefield with respect to a background medium and a per-turbed medium (see study done by Malinowski et al. 2011). In thissensitivity analysis, the entire subsurface of the 1-D medium is re-placed by the density ρback = 2100 kg

m3 and the parameters of layer IV(vP,back = 2000 m

s , QP, back = 25). These homogeneous backgroundmodels are used to obtain the reference data pback. A point scattererwith a size of 0.5 m is placed beneath the first source location ata single grid point at (x0, z0) = (16, 100) m. At (x0, z0), perturba-tions in vP and QP are applied independently of each other. On theone hand, we choose perturbations vP, perturb = (1, 3, 25 per cent),QP, perturb = 1 per cent as well as ρperturb = (1, 3, 25 per cent). On theother hand, we apply QP, perturb = 296 per cent, which corresponds tothe average value QP = 74 of the sub-seafloor area and QP, perturb =∞representing an acoustic scatterer. The scattering-angle interval forthe given acquisition geometry is [−0.6◦, 70.2◦]. A correspondingvirtual receiver array is used to record pback and the perturbed datapperturb. It is located around the scatterer at a defined distance of50 m to exclude the influence of geometrical spreading. We calcu-late a dimensionless data misfit from the scattered partial derivativewavefield with respect to the model parameters m = (vP, QP, ρ) ofbackground and scatterer:

εperturb (m) =

∥∥∥ pperturb−pback

mperturb−mbackmback

∥∥∥2

‖pback‖2

. (14)

As pointed out by Malinowski et al. (2011), the impact of a velocityperturbation on the partial derivative wavefield is of higher signifi-cance compared to a perturbation of the quality factor (Fig. 4a). Thechoice of an identical relative perturbation vP, perturb and QP, perturb

(i.e. the term (mperturb − mback)/mback is identical and ensures afair assessment of model perturbations) results in a misfit differ-ence of more than one order of magnitude (compare graphs forvP, perturb = QP, perturb = 1 per cent in Fig. 4a). However, QP is notsubject to the inversion in this study and, thus, the direct im-pact of model perturbations on the wavefield—which is not bi-ased by the order of magnitude of model parameters—is nowof interest (see Fig. 4b). While a small QP perturbation, such as

QP, perturb = 1 per cent, leads to insignificant data misfit, contrastingscatterers with QP, perturb = 296 per cent or QP, perturb = ∞ result inmisfits which correspond to very small vP perturbations between 1and 3 per cent. In case of vP, perturb and QP, perturb, the scattered am-plitude does not depend on the scattering angle, which agrees withtheory of P-wave scattering in a purely acoustic medium (Wu 1989).Consequently, a QP perturbation acts like a vP perturbation. Apartfrom the amplitude, the phase misfit is almost negligible, which iscaused by the small size of the scatterer compared to the dominantwavelength of 25 m. While the impact of a QP structure on thewavefield is a second-order effect, the impact of QP along the wavepath is more significant. In particular, this propagation effect is ofinterest in the following FWT experiments.

The consideration of density in this work is a relevant simpli-fication which might lead to insufficient model reconstructions inFWT applications influenced by media with inhomogeneous den-sity. While vP modifies both phases and amplitudes, ρ only affectsamplitudes of waveforms. The sensitivity of identical relative per-turbations vP, perturb and ρperturb ranges within the same order of mag-nitude. In particular, at short apertures and thus small scatteringangles, we can observe the maximum impact of ρperturb on the data(dashed graphs in Fig. 4b). From zero scattering angle to maximumscattering angle, the sensitivity εperturb(ρ) decreases by more than30 per cent and thus describes a function of the scattering angle(according to Wu 1989).

3.1.2 Impact of attenuation on the cycle-skipping problem

Apart from the choice of the initial vP model, attenuation and,thus, dispersion might lead to cycle skipping. We compare bothamplitude and phase misfits between visco-acoustic-observed dataand the initial synthetic data of contrasting inversion tests 3 and1 (see Table 1). On the one hand, the combination of the smoothinitial vP model and the true QP model is used to analyse the impactof vP (see test 3). On the other hand, we choose the true vP modeland neglect all information on attenuation to focus on the impact ofQP (see test 1). Figs 5(a) and (b) (test 3) and Figs 5(c) and (d) (test1) show the results of a time-frequency analysis applied to the datafor the central source location. We applied a moving time windowwith a cosine shape and a length of one dominant cycle to theinitial data, computed amplitude and phase spectra and calculated

Figure 4. 1-D experiment: Impact of scatterers in vP, QP or ρ on the partial differential wavefield (a) and the pure wavefield (b). (a) and (b) show the amplitudemisfits (normalized to the amplitude of the unperturbed data) for different relative model perturbations as a function of the scattering angle and with respect toa homogeneous background medium. The normalized amplitude misfit for QP, perturb = ∞ in (a) amounts to 3.2 × 10−7.

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

Acoustic FWT in presence of attenuation 7

Figure 5. 1-D experiment: Time-frequency analysis of the residuals between observed data and the initial synthetic data for inversion test 3 (a and b) and test1 (c and d). Both amplitude misfits (left column) and phase misfits (right column) are averages over the bandwidth of � f = [19.6, 141] Hz.

the average misfit (normalized to the observed data) with respect tothe bandwidth.

The smooth initial vP model (Figs 1a and b) is a quite goodrepresentation of the shallow structures, that is, the seafloor reflec-tion is characterized by a moderate amplitude misfit (Fig. 5a) and anon-existent phase misfit (Fig. 5b). All later events reveal strongerdeviations in amplitudes and phase misfits up to half a cycle. How-ever, phase-misfit discontinuities along the time direction mightindicate the appearance of cycle skipping in test 3.

In contrast, the choice of the true vP model as initial modeland the assumption of an acoustic medium (e.g. compare infiniteQP and QP = 10 in layer II) is reflected in significant amplitudemisfits rather than phase misfits (Figs 5c and d). In particular, atthe seafloor reflection, we can observe high amplitude misfits andthe appearance of phase misfits, even though they are smaller thanπ/6. Altogether, the phase misfits increase up to half a cycle butdo not show any discontinuous behaviour along the time direction.Consequently, we do not expect cycle skipping in test 1.

3.1.3 Inversion tests

Figs 6(a)–(f) show the results obtained by the reference FWT as wellas from acoustic FWT of visco-acoustic data (compare Table 1),which will be discussed in the following. According to Fig. 1(a),all subfigures contain auxiliary plots of the dispersion range. Thevertical section of all inversion results is computed by lateral av-eraging of a representative model area within the interval x = [80,228] m. This avoids the involvement of unreliably recovered veloc-ities, mainly related to the areas close to lateral model boundaries.Fig. 7 depicts the evolution of the corresponding data misfits andmodel errors.

A reliable interpretation of the effects caused by attenuation canbe done by computing a reference result which comprises acousticinversion of pure acoustic observed data. The nearly optimal confor-mity of the true and the final model (Fig. 6a) ensures the resolvingpower of FWT with respect to the given model and geometry. The

reference FWT is characterized by the strongest reduction of bothdata and model error (see Fig. 7). It stops after 1762 iterations.Due to the computation of too small step lengths and the limitedaccuracy of single precision, the model update stagnated.

Test 1 investigates the effect of neglecting QP information in anacoustic FWT applied to visco-acoustic data (Fig. 6b). In spiteof using the true vP model as initial model, the FWT starts atthe highest initial data misfit (Fig. 7a), which is reduced at theexpense of the accuracy of the velocity model. On the one hand,the interface locations are still recognizable, but on the other hand,both a significant model error (increasing from 0 to 4.5 per cent,see Fig. 7b) and an artificial alteration of the velocity model canbe observed. In particular, layer III is smeared heavily. Only layerII is recovered within the maximum range of velocity dispersion.This indicates that the effects of attenuation might be negligiblein near-surface areas of the given model. Although we did notexpect the occurrence of cycle skipping in this test (Figs 5c and d),the choice of a misfit definition being sensitive to amplitudes andthus the wrong explanation of amplitudes results in a failure of theinversion. The FWT of test 1 stops after 61 iterations due to theinability of computing a meaningful step length.

Test 2 represents a common FWT application (Fig. 6c). In this casewe neglect QP and use the smooth initial vP model (Fig. 1a). Thefinal velocity model is artificially recovered by the attempt of theFWT to explain the observed data. While the final vP model showssome improvements, it is still very similar to the initial model.Both the data and the model error are reduced slightly producing asmooth velocity model. Furthermore, the large-scale structures arecomparable to the result of test 1. The FWT of test 2 stops after 77iterations due to the inability of computing a meaningful step length.Although this inversion test already produced an unsatisfactory re-sult, it still represents an ideal case with respect to the considerationof noise-free data. The appearance of noise might hide the effectsof attenuation, particularly at later (multiple) reflection events. The

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

8 A. Kurzmann et al.

Figure 6. 1-D experiment: (a) shows the FWT result for the acoustic reference computation. (b) to (f) illustrate the results of acoustic FWT applied tovisco-acoustic data (tests 1–5). The shaded areas denote phase velocity dispersion.

Figure 7. 1-D experiment: (a) shows the evolution of the L2-based data misfit with respect to the number of iteration. The data misfit of the reference FWT(‘R’) is normalized to its own initial value. The data misfits of tests 1–5 (‘T1’–‘T5’) are normalized to the initial value of test 3 due to its best comparabilitywith the reference FWT. (b) illustrates the evolution of corresponding L1-based model errors with respect to the true vP model. Colour coding is equivalent toFig. 6.

misinterpretation of noisy data in terms of acoustic propagation ef-fects leads to a further degradation of the vP model. Nevertheless,the primary reflections are crucial in reconstructing the large-scalestructures of the 1-D model.

Test 3 includes the true QP model (Fig. 1b) as a passive modelparameter. It can be clearly seen that the result resembles the acous-tic reference result (Fig. 6d). The velocity of layer II is explainedwithin the range of relevant phase velocity dispersion (see Fig. 2c).

Especially in case of low attenuation (layer III and half-space), wedo not observe this effect. Although the initial residuals indicatedthe appearance of cycle skipping (Figs 5a and b), the usage of thetrue QP model allows the correct explanation of amplitude misfits(in contrast to test 1). Thus, test 3 is characterized by a strong reduc-tion of both data misfit and model error (see green plots in Figs 7aand b). This verifies the methodology of combining visco-acousticmodelling with the acoustic inversion scheme, that is, the gradi-ent computation (6) and model update (9) use the relaxed model

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

Acoustic FWT in presence of attenuation 9

Figure 8. Marmousi experiment: (a) shows the true vP model and (b) depicts the initial vP model. (c) and (d) illustrate the true QP model and the smoothpassive QP model.

parameter but are based on the acoustic equations without any mod-ification. The revocation of relaxation (12) yields an inversion resultcomparable to the reference result. The FWT of test 3 stops after4191 iterations due to the threshold stop criterion mentioned above.

Test 4 shows a more realistic FWT application (Fig. 6e). We usethe smooth passive QP model (Fig. 1b). The upper model areas arerecovered quite well. In contrast to the previous result, we can ob-serve a decreasing resolution of the velocity model with increasingdepth. Although the passive QP model is a quite good approximationof the true model, a possible reason is the wrong fit of waveformamplitudes, that is, QP is overestimated in layers II and IV, whileQP is underestimated in layers III and V resulting in a loss of high-frequency information. However, there is a remarkable reduction ofboth data and model error (see Fig. 7) resulting in a qualitativelygood identification of layers and interfaces. The FWT of test 4 stopsafter 195 iterations due to the inability of computing a meaningfulstep length.

Test 5 deals with the simplest case of implementing attenuation(Fig. 6f). Here, we use a homogeneous passive QP model, that is,this model consists of the water layer over a homogeneous half-space with the average quality factor QP = 74. This test yieldsthe smallest misfit reduction among all tests with visco-acousticmodelling (Fig. 7a) and a poor recovery of the velocity modelwhich still is characterized by a high similarity to the initial model(compare Fig. 7b). The FWT of test 5 stops after 111 iterations dueto the inability of computing a meaningful step length.

The observations for all tests coincide with the evolution of boththe data misfit and the model error with respect to FWT iterations(Fig. 7a). The neglection of attenuation in tests 1 and 2 result ininappropriate velocity models. Surprisingly, the usage of true vP

as initial model in test 2 causes a higher initial data misfit. How-ever, both FWT tests end at the same high misfit level failing inthe attempt to explain visco-acoustic data with acoustic modelling.Furthermore, the homogeneous QP model in test 5 yields a misfitevolution being nearly identical to test 2. In spite of incorporatingattenuation, the result is very similar to the velocity model obtainedin test 2. Obviously, in case of the 1-D experiment, the applicationof a homogenous passive QP = 74—which is an incorrect rep-resentation of the subsurface—is insufficient for a successful vP

reconstruction. In contrast, the usage of the true QP model in test 3results in a continuous reduction of the data misfit. This indicatesthe most optimal convergence to the global minimum of the misfitfunction. An acceptable trade-off is achieved by using the smoothQP model in test 4. In practice, a good QP model is usually un-known. Empirical relations can be used to derive a QP model fromthe initial velocity model. However, in general, they do not accountfor all rock types and physical conditions occurring in the givensubsurface.

As mentioned above, a certain choice of a homogeneous pas-sive QP model might cause unsatisfactory results. In addition, weperformed an inversion using a homogeneous passive model withQP = 10 (representing layer II in Fig. 1b). In contrast to test 5,the inverted vP model shows a much better recovery of the upper

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

10 A. Kurzmann et al.

layers. The reconstruction of layer II is nearly identical to the refer-ence result. Layers III and IV show a small vP overestimation andunderestimation, respectively. The average deviation from the truemodel is less than 3 per cent. However, the velocity of layer V ischaracterized by a large overestimation of about 13 per cent. Conse-quently, the passive QP model should rather be a good representationof the shallow structures of the model.

3.2 Synthetic experiment: Marmousi model

Using a modified section of the Marmousi-II model (Martin et al.2006, shown in Fig. 8), we repeat all previous investigations. Thetotal model size is 3 × 10 km with a spatial discretization of�h = 5 m. To reduce computational efforts velocities are clippedto vP = [1.5, 4] km s−1 (Fig. 8), which affects the deep salt layerextended from x ≈ 6 km to x = 10 km. The acquisition geometryconsists of 32 explosive sources and a maximum number of 300hydrophones per source. It forms a marine streamer geometry at thewater surface moving from the right to the left model boundary. Wechoose a streamer length of 5980 m, a receiver spacing of 20 m anda nearest offset of 45 m. However, due to the existence of the rightmodel boundary, only the receiver arrays for shots 20–32 providethe full streamer length, while shot one is equipped with the shorteststreamer containing 18 receivers only. The source time function isa Ricker wavelet with a peak source frequency fpeak = 9 Hz. Therecording time of seismic data is 5.15 s.

The true QP model (Fig. 8c) is derived from the velocity model(Fig. 8a) by applying an empirical vP–QP relation (Hamilton 1972):

1

QP= αP

vP

π f − α2Pv2

P4π f

. (15)

The intrinsic attenuation αP is assigned to the structures of thevelocity model. We used laboratory αP values for the frequencyrange of the given example in marine sedimentary layers (Attewell& Ramana 1966)—ranging from 10−5 to 10−3 m−1. The qualityfactor ranges from 21 in the shallow sedimentary layers to 707 indeeper high-velocity zones.

The approximation of relaxation parameters is based on the ref-erence frequency f0 := fpeak, L = 3 and an average QP, 0 = 62(harmonic mean) within the area beneath the seafloor. The band-width of the seismic data is limited to � f = [3.3, 16.5] Hz. Thiscorresponds to a dynamic range of 2.3 octaves. Consequently, weuse a sufficiently high number of three relaxation mechanisms. Theresulting optimal set of relaxation frequencies is fr, l = (0.1513,1.925, 18.94) Hz. Based on the true model, both acoustic and visco-acoustic data as well as their residuals are shown for shot 24 locatedat x = 2555 m (Fig. 9). For a better illustration, the residual seis-mogram is clipped to ±4 per cent of maximum residual amplitude(Fig. 9a) and data traces are normalized to the maximum ampli-tude of visco-acoustic-observed data (Fig. 9b). Both amplitudes andphases of the waveforms are significantly affected by attenuation.While the strongest amplitude misfit is observed for the seafloor

Figure 9. Marmousi experiment: (a) shows the residual seismogram with true amplitudes, which are clipped to ±4 per cent of the maximum visco-acousticamplitude. The residuals are computed from acoustic and visco-acoustic data recorded at shot location x = 2555 m. (b) illustrates corresponding observedtraces at exemplary offsets. For a better visualization they are normalized independently for each offset. Acoustic and visco-acoustic amplitudes are stillcomparable.

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

Acoustic FWT in presence of attenuation 11

Figure 10. Marmousi experiment: (a)–(f) show the recovered vP modelsfor acoustic reference FWT of pure acoustic data (a) and acoustic FWT ofvisco-acoustic data for tests 1–5 (b)–(f).

reflection, reflections from shallow structures and the refractedwaves, the phase misfit continually increases with traveltime.

For the Marmousi model, we perform the same inversion tests asfor the 1-D experiment (Table 1). The smooth initial vP model is gen-erated by application of a 2-D Gaussian filter (size 1250 × 1250 m,σ = 51) to the sub-seafloor area of the true model vP, true (see Fig. 8b).While we use relation (15) to compute the smooth QP model (seeFig. 8d) from the initial smoothed vP model, the homogeneous pas-sive QP model consists of the water layer over a half-space with theaverage quality factor QP = 62. All inversion tests are decomposedinto multiple stages to reduce non-linearity of the inverse problemand to mitigate the cycle-skipping phenomenon: By applicationof low-pass filters, the inversion is performed for five frequencyranges, moving from low to high frequencies. The peak frequenciesof low-pass-filtered data are not equidistantly spaced, as suggestedby Sirgue & Pratt (2004): fpeak = (1.7, 2.6, 3.7, 5.0, 9.0) Hz.

The inverted velocity models are shown in Fig. 10 and the cor-responding relative deviations from the true model can be found inFig. 11. For a better visualization, the deviation images are clippedto ±20 per cent. The actual maximum range is up to ±50 per cent.Due to the resolution limit of FWT of about half the wavelength,this is caused by the lack of very high frequencies. However, theincorporation of higher frequency contents and thus smaller wave-lengths results in a finer finite-difference discretization (based onthe grid dispersion criterion of approximately 16 grid points perminimum wavelength for a second-order finite-difference scheme).Consequently, in order to handle a reasonable grid size, very small-scale structures and high-contrast interfaces can not be recovered.Furthermore, Table 2 summarizes the data and model errors com-

Figure 11. Marmousi experiment: (a)–(f) show the relative deviation of theresults in Fig. 10 with respect to the true model in Fig. 8(a). For a bettervisualization, the images are clipped at ±20 per cent of the relative deviation.

puted by eqs (13). For all FWT tests, it compares the change betweeninitial and final errors.

Again, as described for the 1-D experiment, an acoustic inversionof pure acoustic data is used as the reference result (Fig. 10a). Con-sidering the quite low peak frequency of the data, we can observe agood match of the reconstructed model and the true model (comparerelative model deviation in Fig. 11a). Both data and model errorsare reduced significantly (see Table 2).

In case of neglecting QP information (tests 1 and 2), the FWT isunable to recover subsurface structures properly from visco-acousticdata (Figs 10b and c). The resolution decreases dramatically withincreasing depth causing a poor final vP model with high modeland data errors (compare initial and final errors for tests 1 and 2 inTable 2, as well as Figs 11b and c).

In the following, we discuss the inversion progress of test 2 inmore detail. Fig. 12 depicts the intermediate inversion results at theend of every frequency-filtering stage. Test 2 shows a satisfactoryreconstruction of large-scale structures during the inversion of thelow-frequency content (see Figs 12a–c). By including higher fre-quencies, artefacts appear in the upper sedimentary structures, whilethere are no improvements in deeper regions. In contrast, we canobserve a destruction of large-scale structures which have alreadybeen recovered (Figs 12e and f). The minimum model error is ob-tained within the fourth stage with a peak frequency fpeak = 5.0 Hz(Fig. 12d). Throughout the remaining inversion progress, the modelerror is increasing continuously, while the data misfit is decreas-ing. For a further investigation of this phenomenon, we repeat test2 without multistage approach. We invert for the full frequencycontent at once and obtain a similar final result. However, theFWT skips the computation of an accurate intermediate result and

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

12 A. Kurzmann et al.

Table 2. Marmousi experiment: List of errors ε with respect to the initial model vP|init and itscorresponding initial data pinit as well as for the resulting model vP|result and correspondingdata presult. Using eqs (13), the errors are computed with respect to the true model vP|true

and observed data pobs. The arrows indicate the strength of error ratios ε(presult) /ε(pinit)and ε(vP|result)/ε(vP|init).

Data error with respect to Model error with respect to Change ofFWT observed data (in per cent) true model (in per cent) Data Model

ε(pinit) ε(presult) ε(vP|init) ε(vP|result) error error

Reference 53.4 0.0659 8.6 2.9

Test 1 28.8 14.2 0.0 7.4

Test 2 45.2 16.1 8.6 8.4

Test 3 12.4 0.0204 8.6 3.1

Test 4 13.2 0.186 8.6 4.1

Test 5 19.2 2.69 8.6 5.0

Figure 12. Marmousi experiment: (a)–(f) show the evolution of the vP model with respect to the inversion progress for the test 2. (a)–(c), (e) and (f) illustrateintermediate vP models at the end of each stage of low-pass filtering. (d) shows the model with the lowest model error (for comparison: the model error of theinitial model is 8.6 per cent). (f) corresponds to Fig. 10(c).

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

Acoustic FWT in presence of attenuation 13

directly produces an artificial model. Only until the sixth iterationthere is a marginal improvement of the model. Apparently, the inver-sion of low-frequency contents is less sensitive to attenuation. Thisobservation coincides with inversion strategies in visco-acousticfrequency-domain FWT applied by several authors. For example,Malinowski et al. (2011) analyse the sensitivity on attenuationin theoretical considerations and numerical experiments. TakamTakougang & Calvert (2011) use a similar multistage approach ina field data application. At low frequencies they only invert for vP.Later on, a combined inversion for vP and QP is applied.

Obviously, tests 1 and 2 demonstrate that the choice of the initialvP model does not affect significantly the outcome of the acousticinversion of visco-acoustic data within the framework of the Mar-mousi experiment. In both cases, the velocity model is artificially al-tered in order to minimize the data misfit (compare increased modelerrors in Table 2). The misfit between visco-acoustic-observed dataand acoustic synthetic data is mapped to the velocity model.

In contrast to tests 1 and 2, the involvement of a QP model im-proves the vP recovery significantly. Using QP, true as passive parame-ter, as performed in test 3, the reconstructed vP model is comparableto the optimal acoustic reference result (Figs 10d and 11d). Test 4employs the smooth QP model and, thus, is the most realistic case.Considering this imperfect QP information, the FWT produces asatisfactory vP model (see final result in Fig. 10e and the modeldeviation in Fig. 11e). Although the smooth QP model does notallow the reconstruction of a vP model with high resolution, thecomparison of test 4 with test 3 only shows a minor increase of themodel error (Table 2). Even the implementation of a homogeneousQP information within the sub-seafloor region yields a surprisinglygood result (Figs 10f and 11f). Apart from decreasing vP resolutionwith increasing depth, the qualitative Marmousi geology is clearlynotable.

In general, the Marmousi experiment resembles the results ofthe 1-D experiment. In case of test 5, we observe different perfor-mances of the FWT. For both examples, there is a nearly identicalreduction of the data misfit. However, in case of the Marmousi ex-periment, there is a significantly stronger reduction of the modelerror (42 per cent versus 25 per cent for the 1-D example). This ob-servation is explained by the choice of the homogeneous QP model.For example, there is a huge QP discrepancy between the secondlayer of the true 1-D model and homogeneous model (true QP = 10versus passive QP = 74). This causes incorrect data and an artificialmodel reconstruction. In contrast, there is a much better match ofthe homogeneous model and the upper structures of the Marmousimodel (the arithmetic mean value of the true model and homoge-neous QP = 62 are equivalent for depths 470 m ≤ z ≤ 1830 m).Concluding, a good QP model should be a good representation ofthe upper sub-seafloor regions.

4 C O N C LU S I O N S

In this work, we investigate the impact of attenuation on 2-D acousticFWT in the time domain. The acoustic inversion scheme is appliedto two visco-acoustic data sets generated for a 1-D structured modeland the Marmousi model. The neglection of attenuation as a passivemodelling parameter causes an unsuccessful recovery of the vP

model. The attenuation-related data misfit is mapped to the velocitymodel by generating remarkable artefacts. If we use the true velocitymodel as an initial model, then the footprint of attenuation can beclearly observed in the artificially altered vP model. Although theimpact of scattering QP structures on the data is a second-order

effect, the influence of attenuation along the wave path in thesereflection experiments is not negligible.

The application of a passive QP model significantly improves thevP model reconstruction. The usage of a smooth QP model resultsin a sufficient near-surface recovery of vP but the resolution isdecreasing with increasing depth. This is caused by an attenuation-related loss of high-frequency information with increasing depth.Depending on the deviation from the true QP model, the choice ofa homogeneous QP model increases the risk of an unsatisfactory vP

recovery (see 1-D experiment).In case of the appearance of soft sediments, the FWT has to

take attenuation into consideration. The availability of a sufficientlygood passive quality factor model allows the reconstruction of areliable velocity model by applying the acoustic inversion scheme.However, such an appropriate good model does not necessarilyhave to be characterized by a high complexity being close to thetrue model. For example, the passive involvement of QP, which isderived from the initial vP model, might significantly improve theresolution of the vP model, provided that the QP model is at least anappropriate representation of the uppermost subsurface structures.In conclusion, it is not advisable to neglect attenuation or to usepotentially poor attenuation information in FWT applications toreal data recorded in marine environments with soft sediments.

A C K N OW L E D G E M E N T S

This work was performed within the programme Geotechnologienfunded by the German Ministry of Education and Research (BMBF)and the German Research Foundation (DFG)—grant 03G0752. Itis also supported by the Verbundnetz Gas AG and the sponsors ofthe Wave Inversion Technology consortium (WIT). In particular, wewant to thank the reviewer S. Operto and the editor J. Trampert forvery helpful and constructive comments. The FWT computationswere performed on the supercomputers JUROPA at Julich Super-computing Centre and HP XC3000 at the Karlsruhe Institute ofTechnology as well as on the bwGRiD cluster computers at severalBaden-Wurttemberg state universities.

R E F E R E N C E S

Askan, A., Akcelik, V., Bielak, J. & Ghattas, O., 2007. Full waveform inver-sion for seismic velocity and an elastic losses in heterogeneous structures,Bull. seism. Soc. Am., 97(6), 1990–2008.

Attewell, P.B. & Ramana, Y.V., 1966. Wave attenuation and internal frictionas functions of frequency in rocks, Geophysics, 31, 1049–1056.

Berenger, J.P., 1994. A perfectly matched layer for the absorption of elec-tromagnetic waves, J. Comput. Phys., 114, 185–200.

Blanch, J.O., Robertsson, J.O.A. & Symes, W.W., 1995. Modeling of aconstant Q: methodology and algorithm for an efficient and optimallyinexpensive viscoelastic technique, Geophysics, 60(1), 176–184.

Bleibinhaus, F. & Hilberg, S., 2012. Shape and structure of the SalzachValley, Austria, from seismic traveltime tomography and full waveforminversion, Geophys. J. Int., 189, 1701–1716.

Bohlen, T., 2002. Parallel 3D viscoelastic finite difference seismic mod-elling, Comput. Geosci., 28, 887–899.

Brenders, A.J. & Pratt, R.G., 2007. Full waveform tomography for litho-spheric imaging: results from a blind test in a realistic crustal model,Geophys. J. Int., 168, 133–151.

Bunks, C., Saleck, F.M., Zaleski, S. & Chavent, G., 1995. Multiscale seismicwaveform inversion, Geophysics, 60(5), 1457–1473.

Carcione, J.M., Kosloff, D. & Kosloff, R., 1988a. Wave propagation simu-lation in a linear viscoacoustic medium, Geophys. J. R. astr. Soc., 93(2),393–407.

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

14 A. Kurzmann et al.

Carcione, J.M., Kosloff, D. & Kosloff, R., 1988b. Wave propagation sim-ulation in a linear viscoelastic medium, Geophys. J. R. astr. Soc., 95(3),597–611.

Causse, E., Mittet, R. & Ursin, B., 1999. Preconditioning of full-waveforminversion in viscoacoustic media, Geophysics, 64(1), 130–145.

Chew, W.C. & Weedon, W.H., 1994. A 3D perfectly matched medium frommodified Maxwell’s equations with stretched coordinates, Microw. Opt.Technol. Lett., 7(13), 599–604.

Crase, E., Pica, A., Noble, M., McDonald, J. & Tarantola, A., 1990. Robustelastic nonlinear waveform inversion: application to real data, Geophysics,55(5), 527–538.

Emmerich, H. & Korn, M., 1987. Incorporation of attenuation into time-domain computations of seismic wave fields, Geophysics, 52(9), 1252–1264.

Fichtner, A., Kennett, B.L.N., Igel, H. & Bunge, H., 2009. Full seismic wave-form tomography for upper-mantle structure in the Australasian regionusing adjoint methods, Geophys. J. Int., 179, 1703–1725.

Grote, M.J. & Sim, I., 2009. Perfectly matched layer for the second-orderwave equation, in Proceedings of 9th Int. Conf. on Math. and NumericalAspects of Wave Propagation (WAVES 2009), Pau, France, pp. 370–371.

Hak, B. & Mulder, W., 2008. Preconditioning for linearized inversionof attenuation and velocity perturbations, in Extended Abstracts, 70thConference and Technical Exhibition, European Association of Geosci-entists & Engineers, Available at: http://www.earthdoc.org/publication/publicationdetails/?publication=9906 (last accessed 12 August 2013).

Hak, B. & Mulder, W.A., 2011. Seismic attenuation imaging with causality,Geophys. J. Int., 184, 439–451.

Hamilton, E.L., 1972. Compressional-wave attenuation in marine sediments,Geophysics, 37(4), 620–646.

Hicks, G.J. & Pratt, R.G., 2001. Reflection waveform inversion using localdescent methods: estimating attenuation and velocity over a gas-sanddeposit, Geophysics, 66(2), 598–612.

Hu, W., Abubakar, A. & Habashy, T.M., 2009. Simultaneous multifrequencyinversion of full-waveform seismic data, Geophysics, 74, R1–R14.

Igel, H., Djikpesse, H. & Tarantola, A., 1996. Waveform inversion of marinereflection seismograms for P impedance and Poisson’s ratio, Geophys. J.Int., 124, 363–371.

Johnston, D.H., 1981. Attenuation: a state-of-art summary, in: Seismic WaveAttenuation, pp. 123–139, eds Toksoz, M.N. & Johnston, D.H., Geo-physics Reprint Series, no. 2: Society of Exploration Geophysicists.

Kamei, R. & Pratt, R.G., 2008. Waveform tomography strategies for imag-ing attenuation structure with cross-hole data, in Extended Abstracts, 70thConference and Technical Exhibition, European Association of Geosci-entists & Engineers, Available at: http://www.earthdoc.org/publication/publicationdetails/?publication=9827 (last accessed 12 August 2013).

Kurzmann, A., Kohn, D. & Bohlen, T., 2008. Comparison of acoustic fullwaveform tomography in the time- and frequency-domain, in ExtendedAbstracts, 70th Conference and Technical Exhibition, European Associa-tion of Geoscientists & Engineers, Available at: http://www.earthdoc.org/publication/publicationdetails/?publication=10539 (last accessed 12 Au-gust 2013).

Liao, Q. & McMechan, G.A., 1995. 2.5D full-wavefield viscoacoustic in-version, Geophys. Prospect., 43, 1043–1059.

Liao, Q. & McMechan, G.A., 1996. Multifrequency viscoacoustic modelingand inversion, Geophysics, 61, 1371–1378.

Liu, H.P., Anderson, D.L. & Kanamori, H., 1976. Velocity dispersion due toanelasticity: implications for seismology and mantle composition, Geo-phys. J. R. astr. Soc., 47(1), 41–58.

Luenberger, D., 1984. Linear and Nonlinear Programming, 2nd edn,Addison-Wesley.

Malinowski, M., Operto, S. & Ribodetti, A., 2011. High-resolution seismicattenuation imaging from wide-aperture onshore data by visco-acousticfrequency-domain full-waveform inversion, Geophys. J. Int., 186, 1179–1204.

Martin, G. S., Wiley, R. & Marfurt, K.J., 2006. Marmousi2—an elasticupgrade for Marmousi, Leading Edge 25, 156–166.

Martin, R. & Komatitsch, D., 2009. An unsplit convolutional perfectlymatched layer technique improved at grazing incidence for the viscoelas-tic wave equation, Geophys. J. Int., 179, 333–344.

Maurer, H., Greenhalgh, S. & Latzel, S., 2009. Frequency and spatial sam-pling strategies for crosshole seismic waveform spectral inversion exper-iments, Geophysics, 74, WCC79–WCC89.

Metivier, L., Brossier, R., Virieux, J. & Operto, S., 2012. Toward Gauss-Newton and Exact Newton optimization for full waveform inversion, inExtended Abstracts, 74th Conference and Technical Exhibition, Euro-pean Association of Geoscientists & Engineers, Available at: http://www.earthdoc.org/publication/publicationdetails/?publication=59194 (last ac-cessed 12 August 2013).

Mora, P., 1987. Nonlinear two-dimensional elastic inversion of multioffsetseismic data, Geophysics, 52(9), 1211–1228.

Mulder, W.A. & Hak, W., 2009. Simultaneous imaging of velocity andattenuation perturbations from seismic data is nearly impossible, inExtended Abstracts, 71st Conference and Technical Exhibition, Euro-pean Association of Geoscientists & Engineers, Available at: http://www.earthdoc.org/publication/publicationdetails/?publication=23986 (last ac-cessed 12 August 2013).

Nocedal, J. & Wright, S.J., 1999. Numerical Optimization, Springer.Pica, A., Diet, J.P. & Tarantola, A., 1990. Nonlinear inversion of seismic

reflection data in a laterally invariant medium, Geophysics, 55, 284–292.

Plessix, R.-E., 2007, A Helmholtz iterative solver for 3D seismic-imagingproblems, Geophysics, 72(5), SM185–SM194.

Polak, E. & Ribiere, G., 1969. Notes sur la convergence de methodes dedirections conjugees, Revue Fr. Inf. Rech Oper, 3(R1), 35–43.

Pratt, R.G., 1999. Seismic waveform inversion in the frequency domain,part 1: theory and verification in a physical scale model, Geophysics, 64,888–901.

Pratt, R.G., Shin, C. & Hicks, G.J., 1998. Gauss-Newton and Full-Newtonmethods in frequency-domain seismic waveform inversion, Geophys. J.Int., 133, 341–362.

Ribodetti, A., Gaffet, S., Operto, S., Virieux, J. & Saracco, G., 2004. Asymp-totic waveform inversion for unbiased velocity and attenuation measure-ments: numerical tests and application for vesuvius lava sample analysis,Geophys. J. Int., 158, 353–371.

Robertsson, J.O.A., Blanch, J.O. & Symes, W.W., 1994. Viscoelastic finite-difference modeling, Geophysics, 59(9), 1444–1456.

Shin, C., Yoon, K., Marfurt, K.J., Park, K., Yang, D., Lim, H.Y., Chung, S.& Shin, S., 2001. Efficient calculation of a partial derivative wavefieldusing reciprocity for seismic imaging and inversion, Geophysics, 66(6),1856–1863.

Sirgue, L. & Pratt, R.G., 2004. Efficient waveform inversion and imaging: astrategy for selecting temporal frequencies, Geophysics, 69(1), 231–248.

Takam Takougang, E.M. & Calvert, A.J., 2011. Application of waveformtomography to marine seismic reflection data from the Queen CharlotteBasin of western Canada, Geophysics, 76(2), B55–B70.

Tarantola, A., 1984. Inversion of seismic reflection data in the acousticapproximation, Geophysics, 49, 1259–1266.

Virieux, J. & Operto, S., 2009. An overview of full-waveform inversion inexploration geophysics, Geophysics, 74(6), WCC1–WCC26.

Wang, Y. & Rao, Y., 2006. Cross-hole seismic waveform tomography—I.Strategy for real data application, Geophys. J. Int., 166, 1224–1236.

Wang, S. & de Hoop, M.V., 2011. On 3D modeling of seismic wave prop-agation via a structured parallel multifrontal direct Helmholtz solver,Geophys. Prospect., 59, 857–873.

Wu, R.S., 1989. The perturbation method in elastic wave scattering, Pureappl. Geophys., 131, 605–637.

Yuan, X., Borup, D., Wiskin, J., Berggren, M. & Johnson, S.A., 1999. Sim-ulation of acoustic wave propagation in dispersive media with relaxationlosses by using FDTD method with PML absorbing boundary condition,IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 46(1), 14–23.

A P P E N D I X : V I S C O - A C O U S T I CM O D E L L I N G I N T H E T I M E D O M A I N

The implementation of attenuation into acoustic or elastic mod-elling has been described by numerous authors (e.g. Liu et al. 1976;Emmerich & Korn 1987; Carcione et al. 1988a,b; Robertsson et al.

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

Acoustic FWT in presence of attenuation 15

1994; Blanch et al. 1995; Bohlen 2002). In modelling of frequency-domain FWT (e.g. Pratt 1999), the intrinsic attenuation is includedquite easily by introducing complex velocities (Johnston 1981). Incontrast, this direct implementation of attenuation is not possiblein the time domain and must be approximated by a suitable rheol-ogy, which is represented by the generalized standard linear solid(GSLS) (Liu et al. 1976). It consists of a parallel connection ofa Hooke and L Maxwell bodies. The Hooke body represents pureacoustic properties. A Maxwell body consists of a Hooke and aNewton body where the latter one characterizes the viscosity ofthe medium. Thus, attenuation is described by L relaxation mecha-nisms.

While an acoustic model is only characterized by density ρ andbulk modulus κ , the visco-acoustic medium is defined by ρ, therelaxed bulk modulus κ r and additional 2L relaxation parametersτ p, l and τ ε, l for each mechanism of the GSLS with l = {1, . . . ,L}. The relaxation times τ p, l and retardation times τ ε, l are requiredfor an appropriate approximation of the quality factor QP, which isproportional to the reciprocal of attenuation αP. We only considerconstant QP models, that is QP( f ) = const. =: QP,0 within the de-sired frequency range fmin ≤ f ≤ fmax of the seismic waveforms. Thegeneral frequency-dependent relation of QP(ω) (with ω = 2π f) andrelaxation parameters is given by

QP(ω, τp,l , τε,l ) =1 − L + ∑L

l=11+ω2 τε,l τp,l

1+ω2 τ2p,l∑L

l=1ω(τε,l −τp,l)

1+ω2 τ2p,l

. (A1)

This equation can be simplified by defining a new L-independentparameter (Blanch et al. 1995)

τP := τε,l

τp,l− 1, (A2)

which replaces the L-dependent retardation times, yielding

QP(ω, τp,l , τP) =1 + ∑L

l=1

ω2τ2p,l

1+ω2τ2p,l

τP∑Ll=1

ω τp,l

1+ω2τ2p,l

τP

. (A3)

In order to achieve an approximation of QP(ω) = const. =: QP,0

the non-linear equation (A3) has to be minimized by applicationof a least-squares optimization algorithm (e.g. Blanch et al. 1995).Hence, the constant QP model is defined by L relaxation times andQP, 0. Instead of defining relaxation times, a more common choiceof relaxation frequencies is ωr,l = 2π fr,l = 1

τp,l.

The implementation of attenuation requires τ P rather than QP. Incase of low attenuation, relation (A3) can be simplified (Emmerich& Korn 1987; Blanch et al. 1995) and rewritten in terms of τ P(QP)and ω = ω0:

1

τP=

L∑l=1

ω0/ωr,l

1 + ω20/ω

2r,l

QP,0 for QP,0 � 1. (A4)

For ωr,l := ω0 and L = 1, this equation reduces to the approximationτ P ≈ 2/QP, 0. Furthermore, the relaxation parameters are used tocompute relaxed frequency-dependent bulk modulus κ r (or relaxedP-wave phase velocity vP, r, respectively):

κr =ρ v2P,ref

(1+

L∑l=1

ω20/ω

2r,l

1 + ω20/ω

2r,l

τP

)−1

for QP,0 � 1, (A5)

with the angular reference frequency ω0 = 2π f0. At reference fre-quency ω0, the reference velocity vP, ref corresponds to the acoustic

phase velocity. A useful choice of f0 is the peak frequency fpeak ofthe source wavelet or the seismic data.

Apart from QP(ω, τ p, l, τ P), the frequency-dependent dispersionvP|ω := vP(ω, τp,l , τP) is computed from the relaxation parameters.Provided that no dispersion occurs at the reference frequency, it isdefined by

vP|ω =vP,ref

⎛⎝

√√√√1+L∑

l=1

ω2τ 2p,l

1 + ω2τ 2p,l

τP−√√√√1 +

L∑l=1

ω20τ

2p,l

1 + ω20τ

2p,l

τP

⎞⎠ .

(A6)

The expressions for minimum and maximum dispersion are:

min vP|ω = limω→0

vP

(ω, τp,l , τP

) = vP,ref

√√√√√ 1

1 + ∑Ll=1

ω20τ2

p,l

1+ω20τ2

p,lτP

,

(A7a)

max vP|ω = limω→∞

vP

(ω, τp,l , τP

) = vP,ref

√√√√√ 1 + L τP

1 + ∑Ll=1

ω20τ2

p,l

1+ω20τ2

p,lτP

.

(A7b)

In addition to the relaxation of the model parameter, attenuationhas to be implemented into the wave equation. This requires thefirst-order partial differential equations (pressure-velocity formula-tion):

p (x, t) = κ (x) ∇ · w (x, t) ,

w (x, t) = 1

ρ (x)∇ p (x, t) , (A8)

where p(x, t) and vector w(x, t) denote pressure and particle veloc-ity field. The incorporation of attenuation based on the GSLS withL relaxation mechanisms yields the following system of partial dif-ferential equations (e.g. Emmerich & Korn 1987; Carcione et al.1988a; Robertsson et al. 1994):

p (x, t) = κr (x) ∇ · w (x, t)

[1+

L∑l=1

(τε,l

τp,l− 1

)]+

L∑l=1

rl (x, t) ,

w (x, t) = 1

ρ (x)∇ p (x, t) ,

rl (x, t) = − 1

τp,l

[κr

(τε,l

τp,l− 1

)∇ · w (x, t) + rl (x, t)

](A9)

with l = {1, . . . , L}. Substitutions based on relation (A2) result in

p (x, t) = κr (x) ∇ · w (x, t) [1+L τP (x)]+L∑

l=1

rl (x, t) ,

w (x, t) = 1

ρ (x)∇ p (x, t) ,

rl (x, t) = − 1

τp,l[κr τP (x) ∇ · w (x, t)+rl (x, t)] , l ={1, . . . , L}.

(1)

Visco-acoustic modelling in the time domain needs additional Lwavefield variables rl (x, t) and L equations. The so-called ‘memoryvariables’ rl (x, t) characterize the memory of the visco-acoustic

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from

16 A. Kurzmann et al.

medium. The system of partial differential equations requires thefollowing initial conditions:

p (x, t = 0) = p (x, t = 0) = 0,

w (x, t = 0) = w (x, t = 0) = 0,

rl (x, t = 0) = rl (x, t = 0) = 0. (A10)

In order to suppress artificial reflections at model boundaries, thetime-domain modelling includes perfectly matched layers (PML) asboundary condition. Their involvement into the system of first-orderwave equations bases on the application of the so-called complexcoordinate stretching (Berenger 1994; Chew & Weedon 1994). Al-though this implementation uses a similar rheology, it differs fromexisting work, such as Yuan et al. (1999). A similar approach forthe viscoelastic case in the frequency domain has been proposed byMartin & Komatitsch (2009). The 2-D visco-acoustic wave equationwith PML boundary condition is given by:

p = κr (1 + L τP) (∇ · w + e · u)

+ (1 + e · φφφ + ϕ)L∑

l=1

rl − (σx + σy + θ

)p, (A11a)

w = 1

ρ∇ p −

[σx 0

0 σy

]w, (A11b)

rl = − 1

τp,l[κrτP (∇ · w + e · u) + (1 + e · φφφ + ϕ) rl ]

− (σx + σy + θ

)rl , (A11c)

u =[

σy∂x 0

0 σx∂y

]w, (A11d)

θ = σxσy, (A11e)

φφφ =[

σx

σy

], (A11f)

ϕ = ψ, (A11g)

ψ = σxσy (A11h)

with l = {1, . . . , L}, unit vector e and additional auxiliary variablesu, θ , φφφ, ϕ and ψ . The PML coefficients σ x and σ y are defined by

σx = σy = 0 within the interior of the model,

σx > 0, σy > 0 within the PML boundary. (A12)

Furthermore, appropriate PML coefficients σ x and σ y can be com-puted from several mathematical functions fx, y, such as quadratic,exponential or cosine functions (e.g. Grote & Sim 2009). This en-sures a smooth transition from the interior of the model domain tothe boundary and within the PML boundary. The coefficients aredefined as follows:

σx,y = σ0 fx,y with 0 ≤ fx,y ≤ 1 and σ0 = − vP ln (R)

B,

(A13)

where vP is the average P-wave velocity, B is the width and R is therelative reflection of the absorbing frame. The relative reflection islimited to the range 0 < R ≤ 1. Useful relative reflection valueshave been estimated empirically: 10−5 ≤ R ≤ 10−3.

Finally, the 2-D modelling comprises (10 + L) equations:

(1) one equation for the update of pressure p,(2) two equations for the update of particle velocities wx and wy,(3) L equations for the update of memory variables rl,(4) seven equations for additional auxiliary PML variables ux, uy,

θ , φx, φy, ϕ and ψ .

at UB

Kiel on A

ugust 31, 2013http://gji.oxfordjournals.org/

Dow

nloaded from