TY - JOUR
AU - Bohlen,, Thomas
AB - SUMMARY Full-waveform inversion (FWI) of surface waves can provide a valuable contribution to near-surface investigations, since they are mainly sensitive to the S-wave velocity and hold a high signal-to-noise ratio. We investigate the performance of the individual 2-D elastic FWI of Rayleigh and Love waves as well as the feasibility of a simultaneous joint FWI of both wave types. We apply these three methods to synthetic data as well as to field data. In synthetic reconstruction tests we compare the performance of the individual wave type inversions and explore the benefits of a simultaneous joint inversion. In these tests both individual wave type inversions perform similarly well, given that the initial P-wave velocity model is accurate enough. In the case of an inaccurate initial P-wave velocity model, we observe artefacts in the reconstructed S-wave velocity models of the Rayleigh wave FWI and the joint FWI. If the P-wave velocity is sufficiently known the joint FWI can further improve the model of the S-wave velocity. In the field data application the Love wave FWI is superior to the Rayleigh wave FWI, possibly due to the insufficient initial P-wave velocity model, while the joint FWI further improves the inversion result also in this case. Full-waveform inversion, Surface wave inversion, Joint inversion, Near-surface seismics 1 INTRODUCTION The analysis of shallow-seismic surface waves provides a valuable contribution to near-surface investigations. Their propagation is mainly influenced by the shear-wave velocity, which is an important geotechnical parameter. The acquisition of shallow-seismic surface waves is simple and cost-efficient, since they can be easily excited and recorded. Furthermore, they show a high signal-to-noise ratio in shallow-seismic field data recordings, which makes them even more attractive to a broad spectrum of near-surface studies. Conventional methods for analysing shallow-seismic surface waves inversion consist of dispersion curves (Dziewonski et al.1969) or Fourier-Bessel expansion coefficients (Forbriger 2003a,b). These approaches transform the seismograms into the velocity/slowness–frequency/wavenumber domain and used by 1-D inversion methods to obtain 1-D subsurface models. However, these approaches are limited to lateral homogeneous or smooth heterogeneous subsurfaces, where in the latter case averaged material properties are obtained. Socco et al. (2010) present an overview of several techniques to overcome the limitation to 1-D subsurface models, like the analysis of data subsets along the profile, where local 1-D depth-dependent models are calculated and subsequently combined to a 2-D subsurface model (Bohlen et al.2004). Nevertheless, all of those methods have a limited lateral resolution and are not applicable in highly heterogeneous media. Full-waveform inversion (FWI) of seismic recordings, as proposed by Lailly (1983) and Tarantola (1984) can reveal 2-D as well as 3-D subsurface structures without limitations regarding the subsurface heterogeneity. FWI makes use of the whole information content included in seismic waves, such as the amplitude and the phase, which allows to achieve a resolution below the size of a wavelength. The main drawback of FWI compared to the conventional methods is the requirement of large computational facilities, which are required for the numerical simulation of wave propagation. While this requirement prevented its application in the past decades, today’s high-performance computing (HPC) systems provide enough computation power to make FWI feasible. Over the last decade, FWI has been successfully applied to a wide range of scales, such as seismology (e.g. Bleibinhaus et al.2007; Fichtner et al.2009), seismic exploration (e.g. Operto et al.2004; Brossier et al.2009) or near-surface investigations (e.g. Gélis et al.2007; Romdhane et al.2011). However, the application of FWI to field data is still challenging, in particular the application to shallow-seismic surface waves is ambitious, since their propagation is highly non-linear in complex earth media. So far, there are only a few publications which present successful 2-D FWI field data applications using shallow-seismic surface waves to obtain shear-wave velocity models. Most of the past publications are dealing with Rayleigh waves (e.g. Groos 2013; Tran et al.2013; Schaefer 2014; Athanasopoulos & Bohlen 2017b), whereas Love waves were only rarely used (Dokter et al.2014; Pan et al.2016). Nevertheless, Xia et al. (2012), who investigated the dispersion curve inversion of Rayleigh and Love waves, observed three main advantages of the Love wave compared to the Rayleigh wave inversion: (1) The inversion of Love wave data is more stable, since they are independent of the P-wave velocity, (2) Love wave dispersion curves are simpler than those of Rayleigh waves, for the same reason as (1), and (3) the dispersion curves of Love waves show a higher signal-to-noise ratio than those of Rayleigh waves. Safani et al. (2005) made similar observations and concluded that Love waves exhibit a higher sensitivity as well as inversion stability and show an improved signal-to-noise ratio in dispersion spectra compared to Rayleigh waves, respectively. Recently, Dokter et al. (2017) and Pan et al. (2018), used the Love waves in a FWI scheme and discussed the benefits of using love waves to reconstruct the shear-wave velocity of the subsurface. However, to the best of our knowledge, there are no studies which compare the performance of the Love wave FWI and the Rayleigh wave FWI in real data, while also investigate the coupling of both inversions in order to perform a joint FWI. 1.1 Objectives and outline The main objectives of this work are (1) to compare the performance of the individual full-waveform inversion of Rayleigh and Love waves and (2) to explore the benefits of a simultaneous joint FWI of both types of surface waves. This paper is organized as follows. First, we briefly review the inverse problem of FWI and provide the adjoint-state gradients for elastic P-SV and SH 2-D FWI. We present our approach to perform a joint FWI of both independent FWI schemes. We then investigate the performance of both individual wave type inversions as well as of the joint inversion in synthetic reconstruction experiments at first and then in a field data application, where we verify the synthetic results by a realistic example. We use a similar test setting during the synthetic and the field data study, in order to make a straight forward comparison. To draw conclusions on the quality and reliability of the field data FWI results, we perform a qualitative comparison to the result of a ground-penetrating radar (GPR) measurement. 2 THEORY The basic concept of FWI is to find a model of the subsurface that describes the observed seismic data in the most accurate way. A model is considered as the best, if it minimizes the misfit between the synthetic data and the observed data. However, due to the non-uniqueness of FWI and the fact that different models correspond to the same data careful interpretation of the results from FWI and coupling with other geophysical methods are necessary. As its name suggests, the FWI uses the full seismic waveforms, hence, the whole information content (e.g. wave amplitude, phase) is taken into account to find an appropriate model. For that reason, the FWI can achieve a resolution below the size of a wavelength. 2.1 Inverse problem To formulate the inversion problem of the FWI, we have to parametrize the model space, m = (m1, ..., mN)T. In seismics two parametrizations are common: a parametrization in terms of density and seismic velocities, which yields |${\bf m}=(\pmb {\rho },{\bf v}_\text{S},{\bf v}_\text{P})^T$| or in terms of density and the Lamé parameters, which yields |${\bf m}=(\pmb {\rho },\pmb {\mu },\pmb {\lambda })^T$|⁠. To obtain synthetic data, dsyn(m), based on a certain model, m, we can use the non-linear forward operator, f: \begin{eqnarray*} {\bf d}_\text{syn}({\bf m})=f({\bf m}). \end{eqnarray*} (1) The data residuals, Δd = (Δd1, ..., ΔdM)T, between the synthetic data, dsyn(m), and the observed data, dobs, are defined as: \begin{eqnarray*} \Delta {\bf d}= {\bf d}_\text{syn}({\bf m}) -{\bf d}_\text{obs}. \end{eqnarray*} (2) To measure the fit of the synthetic data to the observed data, we use the least-squares L2-norm of the data residuals: \begin{eqnarray*} E({\bf m})=\frac{1}{2}\cdot \Delta {\bf d}^T \cdot \Delta {\bf d}, \end{eqnarray*} (3) where E(m) is called misfit or objective function. Thereby, the objective function refers to a summation of the data residuals over all time samples and all source–receiver pairs. Moreover, this definition has a special physical meaning, since it describes the residual energy which cannot be described by the current synthetic model. The aim of the inversion process is to minimize the objective function iteratively and thereby find a model of the subsurface that explains the observed data. In the following, we assume only weak non-linearity of eq. (3) in order to use the Born approximation. In local optimization methods a local minimum of the objective function is searched in the vicinity of an initial model, m0. Therefore, we add a model perturbation, Δm, to the initial model to obtain an updated model: \begin{eqnarray*} {\bf m}={\bf m}_0+\Delta {\bf m}. \end{eqnarray*} (4) We now consider the objective function for this updated model which we expand in a Taylor series up to the second-order accuracy: \begin{eqnarray*} E({\bf m})&=& E({\bf m}_0)+\Delta {\bf m} \left(\displaystyle\frac{\partial E({\bf m}_0)}{\partial {\bf m}} \right)\nonumber \\ && \, +\frac{1}{2}\Delta {\bf m} \left(\displaystyle\frac{\partial ^{2}E({\bf m}_0)}{\partial {\bf m}^2} \right)\Delta {\bf m}^{T} + {\mathcal {O}}({||\Delta {\bf m}||^3}). \end{eqnarray*} (5) To find a minimum of this objective function in the vicinity of the initial model, m0, its derivative with respect to m is required to vanish: \begin{eqnarray*} \frac{\partial E({\bf m})}{\partial {\bf m}}=\frac{\partial E({\bf m}_0)}{\partial {\bf m}}+\Delta {\bf m} \left( \frac{\partial ^{2}E({\bf m}_0)}{\partial {\bf m}^2}\right)\overset{}{=} 0. \end{eqnarray*} (6) Rearranging to the model correction, Δm, gives the desired model update: \begin{eqnarray*} \Delta {\bf m}=-\left( \frac{\displaystyle\partial ^{2}E({\bf m}_0)}{\partial {\bf m}^2}\right)^{-1}\frac{\partial E({\bf m}_0)}{\partial {\bf m}}=-{\bf H}^{-1} \cdot \bigtriangledown _{\bf m} E({\bf m}_0), \end{eqnarray*} (7) where ▽mE(m0) is the gradient of the objective function with respect to the N model parameters mi. The second-order derivative with respect to the model parameters contains the curvature information of the objective function and is called Hessian, H. With eq. (4) and (7) we obtain the model update for iteration K by: \begin{eqnarray*} {\bf m}_{K+1}={\bf m}_{K}+\Delta {\bf m}_{K}={\bf m}_{K}-{\bf H}^{-1}_K \cdot \bigtriangledown _{\bf m} E({\bf m}_K). \end{eqnarray*} (8) This second-order accurate model update is called Newton-method, since it considers the curvature information of the objective function (Nocedal & Wright 2006). This means that a model update with the Newton-method requires the gradient and the Hessian of the objective function have to be calculated. The gradient can be calculated efficiently by the adjoint state method, which we use to derive the gradients shown in Section 2.1.1. However, the Hessian, H, is a dense N × N matrix, where second-order derivatives of the objective function have to be calculated. The calculation of such second-order derivatives can be complicated and computational too expensive for large-scale optimization problems like in the FWI, even on modern HPC systems. To overcome this limitation, we use a quasi-Newton limited-memory BFGS (L-BFGS) method, which avoids an explicit calculation of the Hessian. Instead, the L-BFGS method calculates an approximation to the inverse Hessian implicitly at every iteration by measuring the changes in gradients and models from the n most recent iterations (Nocedal & Wright 2006). We thereby use a Wolfe conditions based step length search in order to ensure the stability of the L-BFGS algorithm. To perform a multiparameter L-BFGS update we follow the approach of Brossier (2011), who proposes a dimensionless multiparameter L-BFGS method, where normalized parameter classes are used within the L-BFGS algorithm. This approach allows to calculate updates for parameters with different units and magnitudes by a single L-BFGS algorithm. In this work, we use the arithmetic mean value for the normalization of the individual parameter classes. 2.1.1 Adjoint-state gradients The calculation of the gradient of the objective function by an actual derivative of the objective function with respect to every single model parameter would need as much forward calculations as model parameters. To overcome this, Tarantola (1984) proposed the usage of the adjoint-state method. This method requires only two forward calculations in order to obtain the descent direction of the objective function. Based on the isotropic and elastic SH and P-SV 2-D wave equations (Virieux 1986, 1984; Lay & Wallace 1995) and by using the adjoint-state approach according to Mora (1987), we derive the gradients for Love (SH) and Rayleigh (P-SV) waves for the density ρ and the Lamé-parameters μ and λ in stress–velocity formulation for one source. We hereby refer to the thesis of Wittkamp (2016) for a detailed derivation of these gradients. For simplicity we omit the temporal and spatial dependencies. The gradients for the SH waves are: \begin{eqnarray*} \frac{\partial E^\text{SH}}{\partial \rho }&=-\int \text{d}t \cdot v_y^{\text{F}}\cdot v_y^{\text{B}}, \end{eqnarray*} (9) \begin{eqnarray*} \frac{\partial E^\text{SH}}{\partial \mu }&=-\int \text{d}t \cdot \frac{\displaystyle\left(\sigma ^{\text{F}}_{xy} \cdot \sigma ^{\text{B}}_{xy}+\sigma ^{\text{F}}_{zy} \cdot \sigma ^{\text{B}}_{zy} \right)}{\mu \cdot \mu }, \end{eqnarray*} (10) and for the P–SV waves: \begin{eqnarray*} \frac{\partial E^\text{P-SV}}{\partial \rho }=&-\int \text{d}t\cdot \left[ v_x^{\text{F}}\cdot v_x^{\text{B}} +v_z^{\text{F}}\cdot v_z^{\text{B}}\right], \end{eqnarray*} (11) \begin{eqnarray*} \frac{\partial E^\text{P-SV}}{\partial \mu }&=&-\int \text{d}t \cdot \Bigg [\displaystyle\frac{\left(\sigma ^{\text{F}}_{xz} \cdot \sigma ^{\text{B}}_{xz}\right)}{\mu \cdot \mu }+\frac{(\sigma ^{\text{F}}_{xx}+\sigma ^{\text{F}}_{zz})(\sigma ^{\text{B}}_{xx}+\sigma ^{\text{B}}_{zz})}{4 \cdot (\lambda +\mu )^2}\nonumber \\ &&\, +\displaystyle\frac{(\sigma ^{\text{F}}_{xx}-\sigma ^{\text{F}}_{zz})(\sigma ^{\text{B}}_{xx}-\sigma ^{\text{B}}_{zz})}{4 \cdot \mu \cdot \mu } \Bigg ], \end{eqnarray*} (12) \begin{eqnarray*} \frac{\partial E^\text{P-SV}}{\partial \lambda }=&-\int \text{d}t \cdot \displaystyle\frac{(\sigma ^{\text{F}}_{xx}+\sigma ^{\text{F}}_{zz})(\sigma ^{\text{B}}_{xx}+\sigma ^{\text{B}}_{zz})}{4(\lambda +\mu )^2}, \end{eqnarray*} (13) where v describes the velocity and σ the stress. The upper indices correspond to the forward (F) propagated incident wavefield or to the backward (B) propagated adjoint (residual) wavefield, respectively. These gradients can be interpreted as a zero-lag cross-correlation between the incident and adjoint wavefield.XXXX For the derivation of the gradients we have chosen a model parametrization in terms of the Lamé-parameters and density, since this is proposed by Mora (1987) to be the simplest one for the derivation of the adjoint-state gradients. As a consequence, we also derived the gradients for these parameters. However, other parametrizations can be resolved with less ambiguities between the individual parameter classes. Especially a parametrization by seismic velocities and density shows less ambiguities (Tarantola 1986; Köhn et al.2012). Hence, we will parametrize the actual FWI experiments by seismic velocities and density |${\bf m}=(\pmb {\rho },{\bf v}_\text{S},{\bf v}_\text{P})^T$|⁠. We derive the gradients for the seismic velocities and the density by using the chain rule of derivatives. The transformation of the parametrization reads (Köhn et al.2012): \begin{eqnarray*} \frac{\partial E}{\partial v_\text{P}}=&\, 2\cdot \rho \cdot v_\text{p} \cdot \displaystyle\frac{\partial E}{\partial \lambda }, \end{eqnarray*} (14) \begin{eqnarray*} \frac{\partial E}{\partial v_\text{S}}=&\, -4\cdot \rho \cdot v_\text{S} \cdot \displaystyle\frac{\partial E}{\partial \lambda } + 2 \cdot \rho \cdot v_\text{S} \cdot \frac{\partial E}{\partial \mu }, \end{eqnarray*} (15) \begin{eqnarray*} \frac{ \partial E}{\partial \rho ^{\prime }}=&\, (v_\text{P}^2-2\cdot v_\text{S}^2)\cdot \displaystyle\frac{\partial E}{\partial \lambda } +v_\text{S}^2 \cdot \frac{\partial E}{\partial \mu } +\frac{\partial E}{\partial \rho }. \end{eqnarray*} (16) Note that due to the change of parametrization, also the density gradient has changed. 2.2 Simultaneous joint inversion In two dimensions the propagation of the P-SV and the SH waves is described by two independent wave equations, thus, the forward as well as the inverse problem of both wave types is decoupled. As a consequence, a manual coupling has to be applied to both individual inversions, which we refer to as joint approach, in order to carry out a joint inversion of both wave types. We call the joint inversion a simultaneous joint inversion, since we will invert the information content of both wave types at the same time in a single inversion. The aim of the joint inversion is to improve the final inversion result and to decrease the vulnerability to local minima and to ambiguities by making use of more information. The joint inversion allows to consider the full information content exploited in a 2-D seismic measurement, since the full 2-D three-component seismic data set can be inverted simultaneously. The first step to couple both individual wave type inversions is to bound both to one single parameter model, because the joint inversion should reveal a single parameter model that accounts for both data sets. Second, we have to merge both individual objective functions, in order to measure the total fit of the synthetic data to the observed data. Moreover, to combine the model update of both individual inversions, we have to apply a joint approach to the gradients of both inversions, since the P-SV as well as the SH inversion are sensitive to the S-wave velocity and to the density. The SH waves are not sensitive to the P-wave velocity, we therefore do not apply a joint approach to the update of the P-wave velocity. In the following, we introduce the joint approach to combine the objective functions and the gradients of both individual wave type inversions. 2.3 Joint objective function To obtain a single measure of the fit between the synthetic and the observed data for the P-SV as well as for the SH waves, we have to combine both individual objective functions. The objective function introduced in eq. (3) describes the residual energy that cannot be explained by the current synthetic model. However, this definition is not practical for the combination of both inversions, since both could contain a different amount of energy, which does not necessarily correspondent to the information quantity included or to the reliability of the specific data set. We therefore weight in the objective function the residual energy with the energy of the observed data: \begin{eqnarray*} E_w({\bf m})=\dfrac{1}{2}\cdot \sum \nolimits _{r,s} \frac{ \Delta {\bf d}^T \cdot \Delta {\bf d}}{{\bf d}_\text{obs}^T \cdot {\bf d}_\text{obs}}. \end{eqnarray*} (17) This weighted objective function is defined as ratio between the residual energy and the energy of the observed data set. A ratio of one would indicate that the residual energy is as big as the energy in the observed data set. Since both wave types should be equally balanced, assuming that the medium is isotropic, we use a simple addition of both weighted objective functions to calculate the joint objective function: \begin{eqnarray*} E^\text{JOINT}({\bf m})=E_w^\text{P-SV}({\bf m})+E_w^\text{SH}({\bf m}). \end{eqnarray*} (18) In the joint inversion this joint objective function is used for the steering through the parameter space. Furthermore, the quasi-Newton L-BFGS method will approximate the Hessian implicitly based on this objective function. 2.4 Joint gradients Both individual wave type inversions are sensitive to the S-wave velocity as well as to the density. Consequently, both return gradients for these two parameter classes, which have to be combined for the sake of a joint inversion. This combination is not as simple as the coupling of the objective functions, due to the lack of an intuitive normalization. The amplitude of a gradient depends on the slope of the objective function, since a gradient is the derivative of an objective function. This means in our case that both gradients would only have a similar amplitude, if both objective functions hold a similar slope. However, this is not necessarily fulfilled since both inversions have their own objective function. We therefore propose a normalized addition of both gradients that is followed by a scaling with the sum of the used normalization factors. We choose the maximum absolute gradient amplitude as normalization factor, respectively. The joint gradient reads: \begin{eqnarray*} \delta \hat{{{\bf g}}}^\text{JOINT}&=&\left( \displaystyle\frac{\delta \hat{{{\bf g}}}^{\text{P-SV}}}{\text{max}(|\delta \hat{{{\bf g}}}^{\text{P-SV}}|)} + \frac{\delta \hat{{{\bf g}}}^{\text{SH}}}{\text{max}(|\delta \hat{{{\bf g}}}^{\text{SH}}|)} \right) \nonumber \\ &&\, \cdot \left(\text{max}\left(|\delta \hat{{{\bf g}}}^{\text{P-SV}}|\right)+\text{max}\left(|\delta \hat{{{\bf g}}}^{\text{SH}}|\right) \right), \end{eqnarray*} (19) where |$\delta \hat{{\bf g}}=$||$\lbrace \frac{\partial E}{\partial \pmb {\rho }}$| or |$\frac{\partial E}{\partial {\bf v}_\text{S}}\rbrace$|⁠. Since the SH waves are not sensitive to the P-wave velocity, we do not calculate a joint gradient for this parameter. This joint gradient approach weights both gradients equally and preserves their amplitude information. The latter is important to provide the amplitude information to the L-BFGS algorithm, which relies on the evaluation of the gradient differences. 3 SYNTHETIC FWI EXPERIMENTS We perform synthetic reconstruction tests to explore the properties of the individual and the joint 2-D elastic full-waveform inversion of Rayleigh and Love waves. For these reconstruction tests we assume a true subsurface model which we then use to generate pseudo-observed seismograms. The knowledge of this model allows us to directly study the reconstruction ability of FWI by comparing the true model to the reconstructed models. In the following, we introduce the subsurface model, the acquisition geometry and the FWI workflow. We choose a similar test setting as during the field measurements, in order to obtain relevant information which will be required in Section 4, where we will present the application of our methodology to a near-surface field data set. 3.1 True and initial models In preparation of the field measurements we choose the synthetic model to be close to the subsurface model at the desired test site of the field measurement. We obtain the subsurface model for this location from previous studies (Groos 2013; Groos et al.2014; Binnig 2015). Their inversion results suggest a predominantly depth dependent 1-D background model, which is superimposed by a shallow small-scale low-velocity trench. This trench represents a low-velocity anomaly which we reconstruct in our field data application. It is a relict of a man made defensive line called 'Ettlinger Linie'. The trench has a triangular shape and a length of 10 m at the surface. The lower edge of the trench lies in a depth of 3.5 m. The maximum anomaly of the trench with respect to the background model is 55 per cent in vS, 28 per cent in vP and 12 per cent in ρ. In order to draw conclusions about the resolution in each individual parameter class as well as to explore trade-off and cross-talk effects between them, we shift the horizontal location of the synthetic trench in each parameter class. In vS, we place the trench in the middle of the model space. In vP, we shift the trench by 5 m to the right and in ρ by 5 m to the left, relative to the vS trench, respectively. We show the assumed true subsurface model in Fig. 3.6, where we focus on an area around the trench. The actual model has a size of 60 m in the horizontal and of 16 m in the vertical direction. In a depth of 6.3 m, we assume a discontinuity corresponding to the groundwater table where all elastic parameters contain a sharp contrast. Below the water table, all three elastic parameters are homogeneous. Above the water table, the true background model consists of a gradient model for vS and of a homogeneous layer for vP and ρ. The gradient in vS is steep in the uppermost part and becomes weaker in a depth of 1 m. The initial models of vS and ρ consist of linear gradient models up to a depth of 9 m. Below 9 m the initial models of this two parameters are identical to the true models. For vS we use a gradient of 20.4 |$\frac{\text{m/s}}{\text{m}}$| and for ρ we use a gradient of 154.2 |$\frac{\text{kg/m}^3}{\text{m}}$|⁠. In contrast to the initial models of vS and ρ, we use a high amount of a priori information for the initial vP model, where we use the true background model as initial model. We thereby assume, that a simple two-layer vP model like in this synthetic test could be obtained in a field measurement in a similar quality, for instance by a P-wave traveltime analysis or by inverting the wavefield sequencially by isolating the contribution of each wave type separetely as shown by Athanasopoulos & Bohlen (2016). We show depth profiles of the initial model and the true model in Fig. 1. Figure 1. View largeDownload slide Depth profiles of the elastic parameters vS, vP and ρ for the true model and for the initial model. At depths, where graphs are not visible, they are identical to the yellow graph. The purple vP profile corresponds to the synthetic study presented in Section 3.6. Figure 1. View largeDownload slide Depth profiles of the elastic parameters vS, vP and ρ for the true model and for the initial model. At depths, where graphs are not visible, they are identical to the yellow graph. The purple vP profile corresponds to the synthetic study presented in Section 3.6. 3.2 Acquisition geometry For the seismic acquisition we use 48 multicomponent receivers, where in case of the P-SV simulations the receivers record the vertical velocity component, vz, and in case of the SH simulations they record the horizontal crossline velocity component, vy. We choose not to include the horizontal inline component in order to be consistent with the field data application which we could not use in the FWI due to low signal-to-noise ratio. However, from synthetic tests and previous studies (Manukyan et al.2012; Maurer et al.2017) we know that this extra component could significally improve the inversion results. For the P-SV simulations the source is a vertical force and for the SH simulations the source is a horizontal force in the crossline direction. These source types correspond to vertical and horizontal hammer blows in field measurements. We distribute both the sources and the receivers along the surface with equidistant spacing of 1  and 2 m, respectively. As a source signal, s(t), we choose a cubed sine: \begin{eqnarray*} s(t)= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\sin ^{3}(f_\text{d} \cdot \pi \cdot t), &\mbox{if } 0 \lt t\lt 1/{f_\text{d}} \\ 0 & \mbox{otherwise} \end{array}\right. , \end{eqnarray*} (20) where fd is the dominant frequency. We set the dominant frequency to 30 Hz, which generates frequencies between 0 and 60 Hz. 3.3 Pseudo-observed data For the simulation, we use a time-domain finite-difference forward solver on a staggered-grid (Bohlen 2002). We discretize the model on an equidistant grid with a spacing of 0.125 m. On the top of the model space we apply a free surface condition based on the mirroring technique of Levander (1988) and at the bottom as well as on the lateral boundaries we assume a C-PML boundary (Komatitsch & Martin 2007). We set the temporal sampling to 2 × 10−5 s and the total recording time to 0.8 s. The temporal derivations are calculated with second-order accuracy and the spatial derivations with a sixth-order FD-stencil. In Fig. 2 we show an exemplary shot-gather for the P-SV velocity component, vz, as well as for the SH component, vy. The corresponding source position for the shown shot-gather lies at the profile coordinate of 6.5 m. We normalized the seismograms trace by trace and apply a bandpass filter between 4 and 60 Hz. The P-SV seismogram, vz, consists mainly of the Rayleigh surface wave and its higher modes. The Rayleigh wave carries most of the seismic energy and is by far the most dominant signal in this recording. P-wave onsets are visible, but they have much smaller amplitude compared to the Rayleigh wave. The SH seismogram, vy, is dominated by surface waves as well, which we identify as Love waves due to the polarization used regarding source/receiver pairs. At the profile coordinate of about 30 m the influence of the low-velocity trench is visible in both recordings. In this area the wave fields are scattered and seismic energy is reflected backwards. Figure 2. View largeDownload slide Trace-wise normalized shot-gather of the pseudo-observed seismic data. The seismograms are bandpass filtered between 4 and 60 Hz, representing the frequency range used in the FWI. The corresponding source is located at the profile coordinate of 6.5 m and is indicated by a yellow star. The centre of the trench lies at the coordinate of 30 m. Figure 2. View largeDownload slide Trace-wise normalized shot-gather of the pseudo-observed seismic data. The seismograms are bandpass filtered between 4 and 60 Hz, representing the frequency range used in the FWI. The corresponding source is located at the profile coordinate of 6.5 m and is indicated by a yellow star. The centre of the trench lies at the coordinate of 30 m. 3.4 FWI setup For the synthetic reconstruction tests we further developed our existing 2-D P-SV FWI code (Köhn 2011) by a SH FWI scheme as well as by the joint approach. For the gradient calculation we use the presented adjoint-state gradients (eqs 9–13), which we transformed to a parametrization in terms of seismic velocities and density (eqs 14–16). We obtain the actual model update by a normalized multiparameter L-BFGS method, where the model and gradient differences of the last 20 iterations (exept if the frequency content is changed) are used to estimate the inverse Hessian. As objective function, we use the weighted L2-error (eq. 17) between the normalized synthetic seismograms and the normalized pseudo-observed seismograms. We apply the normalization trace-wise (using the maximum absolute value), in oder to be consistent with the field data inversions, where we use this to be able to account for uncertainties and differences in the coupling of individual geophones to the ground (Maurer et al.2012). In the P-SV case, we calculate the objective function for the vertical displacement component and in the SH case for the horizontal crossline displacement component. We obtain the displacement seismograms from the recorded velocity seismograms by numerical integration. To precondition the shot-wise gradients we apply circular source tapers, which decay within 5 gridpoints from a value of one to zero at the actual source position, where zero refers to absolute damping. Moreover, we use an approximation of the diagonal elements of the Hessian as a preconditioner for the gradient. This approximation is based on the sum of the amplitudes, ui, of the forward modeled incident wavefield at each gridpoint. The influence of the receivers is included by a geometrical estimation of the receivers Green’s functions. The approximation for a single shot reads (Plessix & Mulder 2004; Wehner et al.2015): \begin{eqnarray*} {\bf H}_a^{-1}({\bf x}_s,{\bf x})&=&\Bigg [\epsilon +\int \text{d}t |u_i({\bf x}_s,{\bf x},t)|^2 \cdot \Bigg (\text{asinh}\left(\frac{x_r^\text{max}-x}{z}\right) \nonumber \\ && \, -\text{asinh}\left(\frac{x_r^\text{min}-x}{z}\Bigg ) \right)\Bigg ]^{-1}, \end{eqnarray*} (21) where |$x_r^\text{max}$| and |$x_r^\text{min}$| correspond to the maximum and minimum horizontal receiver position for the source location, xs, respectively. To stabilize the inversion of the expression above, a water level, ε, is added. We choose this water level empirically to 5 × 10−3 for the SH waves and to 5 × 10−2 for the P-SV waves. We calculate this preconditioner shot-wise and apply it normalized to the shot-wise gradients. In addition, the gradients are smoothed with a 2-D median filter, which has a size of 50 cm (4 gridpoints) in both spatial directions. To ensure stability of the forward solver as well as to obtain a physical meaningful elastic parameter model, we force a minimum vP/vS ratio of 1.5 during the inversion by increasing the P-wave velocity, if necessary. We control the multistage inversion process by an automatic workflow implementation. We use this workflow implementation to increase the frequency content of the data set gradually from 4 to 60 Hz, by increasing the corner frequency, LC, of a fourth-order Butterworth low-pass filter. Moreover, we apply a sequential update strategy to the multiparameter inversion by the workflow implementation. In the first iterations we only allow updates of vS, until an automatic abort criterion, AC, is reached. Afterwards, in case of the individual P-SV as well as of the joint inversion we update vS and vP simultaneously, again until an automatic abort criterion is reached and finally we use a full multiparameter inversion. We decided to use this sequential update of the parameter classes to account for the different sensitivities of the surface waves to the individual elastic parameters. Their propagation is mainly influenced by the S-wave velocity. We divide the workflow in eight separate stages, which we describe in Table 1. As abort criterion we use the reduction of the objective function from the second last to the current iteration: \begin{eqnarray*} AC_K(\text{in per cent}) \gt \left(\frac{E_{K-2}-E_K}{E_{K-2}}\right). \end{eqnarray*} (22) Table 1. Workflow used in the synthetic FWI experiments. Each stage is applied to the inversion until the automatic abort criterion AC is reached (eq. 22). The update columns indicate which of the specific elastic parameter is updated (yes=1) or not (no=0). LP represents the corner frequency of the low-pass frequency filter. Stage Update LP in Hz AC in  per cent vS vP ρ 1 1 0 0 10 20 2 1 1 0 10 10 3 1 1 1 10 1 4 1 1 1 20 1 5 1 1 1 30 1 6 1 1 1 40 1 7 1 1 1 50 1 8 1 1 1 60 0.5 Stage Update LP in Hz AC in  per cent vS vP ρ 1 1 0 0 10 20 2 1 1 0 10 10 3 1 1 1 10 1 4 1 1 1 20 1 5 1 1 1 30 1 6 1 1 1 40 1 7 1 1 1 50 1 8 1 1 1 60 0.5 View Large Table 1. Workflow used in the synthetic FWI experiments. Each stage is applied to the inversion until the automatic abort criterion AC is reached (eq. 22). The update columns indicate which of the specific elastic parameter is updated (yes=1) or not (no=0). LP represents the corner frequency of the low-pass frequency filter. Stage Update LP in Hz AC in  per cent vS vP ρ 1 1 0 0 10 20 2 1 1 0 10 10 3 1 1 1 10 1 4 1 1 1 20 1 5 1 1 1 30 1 6 1 1 1 40 1 7 1 1 1 50 1 8 1 1 1 60 0.5 Stage Update LP in Hz AC in  per cent vS vP ρ 1 1 0 0 10 20 2 1 1 0 10 10 3 1 1 1 10 1 4 1 1 1 20 1 5 1 1 1 30 1 6 1 1 1 40 1 7 1 1 1 50 1 8 1 1 1 60 0.5 View Large 3.5 Results In the following, we discuss the result of each reconstruction experiment individually. In Fig. 3 we show a direct comparison of the final reconstructed models as well as the true model, where we focus on an area around the trench. Figs 4(a)-(c) shows a comparison of seismograms for the initial model, for the final reconstructed model and the pseudo-observed data set for each inversion scheme, respectively. Figure 3. View largeDownload slide Inversions results of the synthetic reconstruction experiments. Comparison of the true subsurface model (first column) with the reconstructed model of individual Love wave FWI (second column), individual Rayleigh wave FWI (third column) and simultaneous joint FWI (fourth column). The elastic parameters vS, vP and ρ are shown row-wise from top to bottom, respectively. Figure 3. View largeDownload slide Inversions results of the synthetic reconstruction experiments. Comparison of the true subsurface model (first column) with the reconstructed model of individual Love wave FWI (second column), individual Rayleigh wave FWI (third column) and simultaneous joint FWI (fourth column). The elastic parameters vS, vP and ρ are shown row-wise from top to bottom, respectively. Figure 4. View largeDownload slide Data fit of the synthetic reconstruction experiments. Comparison of the seismograms for the initial model, for the final reconstructed model and the pseudo-observed data set for all three inversions (a–c). Evolution of the objective function over the iterations (d). Figure 4. View largeDownload slide Data fit of the synthetic reconstruction experiments. Comparison of the seismograms for the initial model, for the final reconstructed model and the pseudo-observed data set for all three inversions (a–c). Evolution of the objective function over the iterations (d). 3.5.1 Individual Love wave FWI The individual Love wave FWI reconstructed the vS model successfully by recovering the low-velocity trench sharply and in full extension, as presented in Fig. 3 (second column). At the depth of 6.3 m, the water table is smoothly visible. The SH waves are not sensitive to the P-wave velocity, thus, the vP model is not updated. The reconstruction of the ρ model is surprisingly well, especially if we consider that we use trace-wise normalized seismograms to calculate the objective function and the fact that the impact of density is mainly to the absolute wave amplitude as a function of offset. Nonetheless, the Love wave FWI reconstructed the trench in the ρ model satisfactorily in both its size and density value, but the contour of the trench is less sharp compared to the vS model. The water table is not sharply restored in the ρ model. As a result of a cross-talk by vS, the reconstructed ρ model shows a footprint of the vS anomaly. Especially the outline of the vS anomaly is clearly visible within the ρ model. The anomaly in the ρ model is not visible in the reconstructed vS model, which indicates that the individual Love wave FWI could restore the vS model with less ambiguity than the ρ model, due to a higher sensitivity of the Love wave to vS than to ρ. The fit of the final synthetic seismograms (red) to the the observed seismograms (black) is very well at all offsets and times (Fig. 4b). A remaining residual between the synthetic and observed seismograms is only hardly visible. The decrease in the objective function is smooth and in total the inversion reduced the misfit by four orders of magnitude, as shown in Fig. 4(d). At each new workflow stage the algorithm increases the frequency content and the objective function increases as well, due to the fact that more information is considered for the residual calculation. 3.5.2 Individual Rayleigh wave FWI The individual Rayleigh wave FWI revealed the trench in the vS model accurately, where both its velocity value and shape are correct (see Fig. 3, third column). In a depth of 6.3 m, the inversion imaged the water table sharply. Thereby it is likely that the Rayleigh wave FWI benefits from the sharp water table included in the initial vP model. Nevertheless, the vS model suffers from small-scale artefacts that are present inside the trench and especially to the right side of it. The reconstruction of the trench in the vP model is satisfactorily. The contour is not clearly visible and vertically orientated artefacts are observed inside the trench. These artefacts could be caused by wrong P-wave velocities in combination with receiver related artefacts, since they occur directly underneath the receiver positions at areas with wrong P-wave velocities and not where the initial model contains the true velocity. We observed similar artefacts in the case study on the influence of an inaccurate initial vP model, which we will present in Section 3.6. The artefacts within the vP trench correspond to the artefacts in the vS model, where the vS model suffers from a cross-talk by vP. The vP model itself suffers from a slight footprint of the vS trench. In general, we expect the resolution in the vP model to be lower compared to the vS model, due to the longer wavelengths of the Pwaves than of the Swaves. The reconstructed ρ model matches the true model satisfactorily. The inversion recovered the shape of the trench in the ρ model sufficiently, but reproduced the density value slightly higher than its actual value. Moreover, small-scale artefacts are present in the ρ model at the position of the vS and vP trench, most likely caused by a cross-talk. The water table is visible as a sharp contrast, where the Rayleigh wave FWI could again benefit from the initial vP model. The fit of the synthetic seismograms to the pseudo-observed seismograms is very well, since the residuals are nearly vanished. The inversion decreased the misfit smoothly in each frequency stage and reduced the objective function in total by four orders of magnitude, as shown in Fig. 4(d). 3.5.3 Simultaneous joint FWI The simultaneous joint FWI reconstructed the trench in the vS model very well in terms of shape and velocity values (see Fig. 3, fourth column). The inversion imaged the water level sharply, where it could benefit from the sharp contrast included in the initial vP model. There are no cross-talk effects of the vP or the ρ model visible in the vS model. The reconstruction of the trench in the vP model is acceptable, although the outline of the trench is not as sharp as in the case of the vS model. The final vP model as well as the ρ model suffer from a light footprint of the vS trench. Besides this cross-talk effect, the inversion recovered the ρ model successfully in its shape and value. The reconstructed trench is sharp and the density values matches the true value. The water table is sharply visible in the ρ model, where a positive influence from the initial vP model is most likely. The joint FWI fitted the P-SV seismograms as well as the SH seismograms to the pseudo-observed seismograms without significant residuals. The joint objective function could be reduced by four orders of magnitude, as shown in Fig. 4(d). Compared to the individual wave type inversions the joint FWI was able to reduce the individual objective functions of the Love and the Rayleigh wave inversion even further. More precisely, in the last iteration the absolute misfit value of the P-SV and the SH waves is lower in the case of the joint FWI than in the case of the individual wave type inversions. However, since we calculated the joint objective function as a sum of both individual misfits (Fig. 4d, dotted lines), the total joint objective function is higher than both individual ones. 3.6 Influence of the initial P-wave velocity In this part, we discuss the influence of the initial P-wave velocity model. In shallow seismic applications (i.e. upper 30 m), the amplitudes of the shear waves are much higher than the ones of the compressional Pwaves. In these particular cases the FWI fails to properly incorporate the contribution of the P waves, as the misfit calculation is exclusively dominated by the high amplitude shear waves. Nevertheless, as the full-waveforms contain the P-wave onsets and the fact that the Pwave is mainly influenced by the vP model, the FWI theoretically could reconstruct the vP model. Athanasopoulos & Bohlen (2016) showed how an inaccurate vP model could affect also the vS and suggested a sequential FWI strategy in which initially the low-amplitude refracted P waves are inverted and then subsequently the full wavefield including the Rayleigh waves. Another strategy to improve the vP model would be to apply a sequential approach. First obtain a vS model from either the joint or the SH FWI and then use this as an initial model to invert the vP model from the P-SV waves. Finally, an approach that one could use in order to improve the vP considering that little information can be extracted from the data (e.g. very shallow seismics), is to constrain the model updates through the structural similarity between model parameters. Manukyan et al. (2018) recently shown that such constraints, as for example between vP, vS and density can improve the reconstructed from FWI models. In our case, an inaccurate initial vP model does not influence the individual Love wave FWI, as the P-wave velocity does not affect SH waves. Since the individual Rayleigh wave FWI is highly influenced by an inaccurate vP model (Fig. 1), this by extend will influence the joint FWI (Wittkamp 2016) as shown in Fig. 5. The reconstructed vS model suffers from circular low-velocity artefacts around the receiver positions. The shape of the trench in the reconstructed vS model deviates slightly from the true model, but the velocity inside the trench is rough. The vP model is reconstructed inaccurately and does not fit the true model. These artefacts are also foot-printed to the vS model. The updated ρ model contains systematically too high density values. Figure 5. View largeDownload slide Inversions results of the synthetic reconstruction experiments for an inaccurate and accurate initial vP model of the simultaneous joint FWI (first and second column, respectively). The elastic parameter vS, vP and ρ are shown row-wise from top to bottom, respectively. Figure 5. View largeDownload slide Inversions results of the synthetic reconstruction experiments for an inaccurate and accurate initial vP model of the simultaneous joint FWI (first and second column, respectively). The elastic parameter vS, vP and ρ are shown row-wise from top to bottom, respectively. 4 FIELD DATA APPLICATION In this section, we present the application of the individual and the joint 2-D elastic full-waveform inversion of Rayleigh and Love waves to a near-surface field data set. We recorded a field data set on an airfield near Karlsruhe (Germany). Previous Rayleigh wave FWI studies (Groos et al.2017) took place on the same test site and proved the suitability of it for 2-D FWI. They propose a predominantly depth dependent subsurface that is superimposed by a shallow small-scale low-velocity trench. We assumed such a 2-D subsurface model in the synthetic experiments, where we verified the theoretical applicability of both the individual wave type inversions as well as the simultaneous joint FWI. In this experiment, we investigate the applicability of all three inversions to the recorded near-surface field data set and evaluate their performance. 4.1 Test site The location of the test site is on a glider airfield in Rheinstetten near Karlsruhe (Germany). Fig. 6 shows an overall map of the area. This test site exhibits a planar surface and does not suffer from humanistic noise in the direct surrounding. The geological map by Hüttner et al. (1968) states that the subsurface consists of layered fluviatile sediments of the late Pleistocene. Several shallow-seismic studies were carried out on this area, which provide further information on the subsurface. Groos (2013) and Schaefer (2014) performed a dispersion curve inversion and a FWI of Rayleigh waves on the northwest part of the airfield and propose a predominantly depth dependent 1-D subsurface. Lüttschwager (2014) investigated the northeast corner of the runway and discovered a shallow small-scale low-velocity anomaly (trench) that proceeds straight from the northwest to the southeast. Binnig (2015) confirmed this hypothesis by a 2-D Rayleigh wave FWI. His results suggest that this low-velocity trench locally superimposes the 1-D subsurface proposed by the previous studies. According to historic recordings this trench can be identified as the 'Ettlinger Linie'. It was originally excavated to serve as a line of defense and was refilled several decades ago (Lang 1907). Outside the borders of the airfield the 'Ettlinger Linie' is still uncovered and traceable, which allows to easily interpolate the course of the refilled trench within the airfield. Such a subsurface structure suites well for 2-D FWI experiments. The low-velocity trench proceeds straight for about 70 m and superimposes the lateral homogeneous background subsurface locally. Hence, the assumption of a 2-D subsurface is valid in the case that the acquisition profile crosses the trench vertically, which is important since the 2-D FWI cannot account for signals, such as reflections, caused by anomalies located off the 2-D profile. Figure 6. View largeDownload slide Overview map of the test site. The red line corresponds to the interpolated course of the 'Ettlinger Linie' and the white line denotes the acquisition profile. Source: Google Earth (AeroWest, GeoBasis-DE/BKG). Figure 6. View largeDownload slide Overview map of the test site. The red line corresponds to the interpolated course of the 'Ettlinger Linie' and the white line denotes the acquisition profile. Source: Google Earth (AeroWest, GeoBasis-DE/BKG). 4.2 Acquisition geometry To image the cross-section of the 'Ettlinger Linie' with the 2-D FWI, we placed the acquisition profile to cross the interpolated course of the trench vertically, as illustrated in Fig. 6. We shifted the profile to contain the trench in its centre, in order to obtain a high ray coverage within the low-velocity trench. The orientation of the profile is from northeast (marker one) to southwest (marker two). For the seismic recording we used 48 three-component geophones with 4.5 Hz eigenfrequency. We set the geophone spacing equidistantly to 1 m and adjusted the local orientation of the geophones to the profile, in order to ensure an accurate recording of the horizontal component. The total length of the receiver line was 47 m. For the P-SV data set we recorded the vertical particle velocity and for the SH data set we recorded the horizontal crossline component. We set the spacing between the 24 sources to 2 m, where the first source was located between the first and the second receiver. All source positions were located within the receiver line. The source–receiver offset ranges from 0.5 to 46.5  m. We used vertical hammer blows on a steel plate to excite the P-SV data set and horizontal hammer blows in the crossline direction on a steel source rack to excite the SH data set. 4.3 Observed data The total recording time is 1.5 s and the temporal sampling is 2.5 × 10−4 s. During the measurement we stacked the data with a fold of five to enhance the signal-to-noise ratio. However, the raw field data set is not appropriate for a direct application of the 2-D FWI, thus, we perform a few preparatory steps (Groos et al.2017). First of all, we shorten the data to 0.5  s, due to the absence of any significant energy at later recording times. We then upsample the data by a spline interpolation to a sampling of 1.4 × 10−5 s, in order to satisfy the stability criterion of the finite-difference forward solver. To suppress signals before the actual wave onset, we apply a muting at the beginning of each trace. In order to avoid non-causal effects in the inverted source–time function (see Section 4.5), we delay the whole data set by 0.02 s, which results in a total time length of 0.52 s. Furthermore, we perform a 3-D to 2-D transformation, which is necessary, due to the fact that the recorded wavefields are excited by hammer blows that act like point-sources, while the 2-D forward solver assumes line-sources. We transform the field data to an equivalent 2-D line-source by a trace-wise convolution with |$\sqrt{t^{-1}}$| followed by a multiplication with |$r\cdot \sqrt{2}\cdot \sqrt{t^{-1}}$|⁠, where t denotes the traveltime and r the offset. This transformation is introduced as direct-wave transformation by Forbriger et al. (2014). An exemplarily shot gather of the preprocessed data set is shown in fig. 7. Figure 7. View largeDownload slide Trace-wise normalized shot-gathers of the preprocessed field dataset for the first source located at the profile coordinate of 10.5 m and indicated by a yellow star. The seismograms are band-pass filtered between 4 and 130 Hz, representing the frequency range used in the FWI. Traces near the source were overdriven and thus muted. Figure 7. View largeDownload slide Trace-wise normalized shot-gathers of the preprocessed field dataset for the first source located at the profile coordinate of 10.5 m and indicated by a yellow star. The seismograms are band-pass filtered between 4 and 130 Hz, representing the frequency range used in the FWI. Traces near the source were overdriven and thus muted. We normalized the seismograms trace-wise and applied a bandpass filter between 4  and 130 Hz. The main frequency content is located between 10  and 100 Hz. There is no significant difference in the frequency content of both wave types. The P-SV seismogram, vz, is dominated by the Rayleigh wave. The Rayleigh wave is visible in the fundamental mode as well as in several higher modes. The direct and the refracted Pwave can be identified as well, however they have much smaller amplitude compared to the Rayleigh wave. The SH seismogram, vy, is dominated by the Love wave, which is present in the fundamental mode. Compared to the Love wave the direct as well as the refracted P wave have smaller amplitude and are only slightly visible. After the preprocessing steps that we carried out the P-SV as well as the SH field data set allow for the 2-D full-waveform inversion. 4.4 Initial model Since we use a local optimization method for the full-waveform inversion, we have to estimate initial models for the parameters vS, vP and ρ. These models have to predict the main wave phases well enough to allow local convergence of the inversion and avoid getting stack in local minima. Moreover, we have to derive initial models for the quality factors QS and QP, in order to describe the attenuation properties of the subsurface. The dissipative properties will not be updated during the inversion, instead we will use them as passive model parameters. In the following, we assume the initial models to vary only with depth, since we expect the background model to be predominantly depth dependent. First, we perform a P-wave traveltime analysis to obtain an initial model for vP. We therefore pick and evaluate the onsets of the direct and the refracted Pwave. The obtained model consists of two layers, where the interface lies in a depth of 6.1 m. For the upper layer we calculate a P-wave velocity of 335  m s–1 and for the lower half space we calculate a velocity of 2284  m s–1. We assume that the sharp contrast in a depth of 6.1 m corresponds to the groundwater table. To obtain an initial estimate for the ρ model from the vP model we use the empirical Gardner’s relation (Gardner et al.1974): \begin{eqnarray*} \rho =0.31 \cdot v_\text{P}^{0.25}, \end{eqnarray*} (23) The obtained ρ model has a density of 1325  kg m–3 in the upper layer and 2142  kg m–3 in the lower half space. We obtain an initial vS model based on the minimum and maximum apparent phase velocities of the fundamental modes of the surface waves, which we derive from a dispersion curve analysis. Based on these two values, we created a simple gradient model. We were able to predict all main phases of the Rayleigh and Love waves for a S-wave velocity of 140 m s–1 at the surface and 340  m s–1 in a depth of 9 m. We show the initial model for the three elastic parameters in Fig. 8 (first column). Furthermore, it is important to estimate an attenuation model, since inelastic damping has a significant influence to shallow-seismic recordings (Groos et al.2014). We use two assumptions to obtain the Q-values from the field data set: (1) We assume that the Q-values for P and S waves are identical and (2) that a constant Q-value is sufficient for the whole model space. Then, we approximate a Q-value by calculating and comparing the misfit between synthetic and observed data for different Q-values using the initial models. Following Groos et al. (2014), we calculate a Q-value of 15 and construct the attenuation model for the viscoelastic forward modelling by a Generalize Standard Linear Solid (GSLS) with three relaxation mechanisms (Bohlen 2002). Figure 8. View largeDownload slide Inversions results of the field data application. Comparison of the initial model (first column) with the final model of the individual Love wave FWI (second column), the individual Rayleigh wave FWI (third column) and the simultaneous joint FWI (fourth column). The elastic parameters vS, vP and ρ are shown row-wise from top to bottom, respectively. Figure 8. View largeDownload slide Inversions results of the field data application. Comparison of the initial model (first column) with the final model of the individual Love wave FWI (second column), the individual Rayleigh wave FWI (third column) and the simultaneous joint FWI (fourth column). The elastic parameters vS, vP and ρ are shown row-wise from top to bottom, respectively. 4.5 FWI setup For the field data application we use a similar configuration as for the synthetic FWI experiments, thus, we discuss it only briefly here. The objective function is again the weighted L2-error (eq. 17) between the synthetic and observed seismograms, after normalizing them trace by trace. The model space has a size of 560 gridpoints in the horizontal direction and 160 gridpoints in the vertical direction, resulting in the actual dimensions of 70 m x 20 m (grid spacing is 0.125 m). The location of the first receiver is at 10 m (Fig. 6) and of the last receiver at 57 m. The total propagation time of 0.52 s is identical to the time length of the observed seismograms. 4.5.1 Source wavelet estimation In our field data measurements the actual excited source wavelets were not directly recorded, therefore they represent additional unknowns of the inverse problem (Pratt 1999). To mitigate the effect of an unknown source wavelet, we perform a separate source time function inversion. For the estimation of a source–time function we calculate a wavelet correction filter by a stabilized deconvolution of the observed seismograms with the synthetic seismograms (Virieux & Operto 2009; Groos et al.2014). By then applying a convolution of the synthetic wavelet with the estimated wavelet correction filter we estimate the source wavelet, which is updated only at the first iteration of each stage. However, the improved wavelet does not necessarily represent the actual source wavelet excited in the field measurement, instead it is the wavelet that minimizes the residuals between synthetic and observed seismograms. Therefore, the estimated wavelet might suffer from a trade-off as it could account for residuals caused by an inaccurate parameter model. As initial guess we use a cubed sine wavelet with a dominant frequency of 100 Hz. 4.5.2 Preconditioning To precondition the shot-wise gradients we apply once more circular tapers around the sources and the approximation to the diagonal elements of the Hessian (eq. 21) as preconditioner for the gradients, as in the synthetic experiments. In order to smooth the inversion result, we apply a 2-D median filter to the gradients, where the filter has a size of 1 m (8 gridpoints). We force again a minimum vP/vS ratio of 1.2 by increasing the P-wave velocity, if necessary. We use a normalized multiparameter L-BFGS method, however for frequencies above 50 Hz the L-BFGS method did not converge, since no appropriate step length could be found. The reason could be that the initial approximation of the Hessian identity matrix (Nocedal & Wright 2006) is not accurate enough to allow further convergence. This problem might be solved by providing an external first guess of the size of the Hessian to the L-BFGS algorithm, for example by the second-order adjoint method (Fichtner & Trampert 2011). To overcome this issue in the field data inversions, we switch the optimization method to a conjugate gradient method for frequency stages above 50 Hz. In this case we use an inaccurate step length search in combination with a parabolic fitting to estimate an appropriate step length (Nocedal & Wright 2006), which provided stable updates also for frequencies above 50 Hz. 4.5.3 Multistage approach We increase the frequency content of the dataset gradually from 4 to 130 Hz by lifting up the corner frequency of a low-pass filter in steps of 5 Hz. In addition, we again apply a sequential update strategy to the multiparameter inversion. In the first iterations we only allow updates of vS and vP (for Love FWI only vS), until the automatic abort criterion is reached and finally we use a full multiparameter inversion. In total, the workflow is divided in 26 separate stages, which we present in Table 2. Table 2. Workflow used in the field data FWI. Each stage is applied to the inversion until the automatic abort criterion AC (eq. 22) is reached. The update column indicates which of the specific elastic parameter is updated (yes=1) or not (no=0). The parameter LP represents the corner frequency of the low-pass frequency filter. The method column indicates whether the L-BFGS or the conjugate gradient method is used for optimization. Stage Update LP in Hz AC in  per cent Method vS vP ρ 1 1 1 0 10 10 L-BFGS 2 1 1 1 10 1 L-BFGS 3–9 1 1 1 Increment of 5 1 L-BFGS 10 1 1 1 50 1 L-BFGS 11 1 1 1 55 1 CG 12–25 1 1 1 Increment of 5 1 CG 26 1 1 1 130 0.5 CG Stage Update LP in Hz AC in  per cent Method vS vP ρ 1 1 1 0 10 10 L-BFGS 2 1 1 1 10 1 L-BFGS 3–9 1 1 1 Increment of 5 1 L-BFGS 10 1 1 1 50 1 L-BFGS 11 1 1 1 55 1 CG 12–25 1 1 1 Increment of 5 1 CG 26 1 1 1 130 0.5 CG View Large Table 2. Workflow used in the field data FWI. Each stage is applied to the inversion until the automatic abort criterion AC (eq. 22) is reached. The update column indicates which of the specific elastic parameter is updated (yes=1) or not (no=0). The parameter LP represents the corner frequency of the low-pass frequency filter. The method column indicates whether the L-BFGS or the conjugate gradient method is used for optimization. Stage Update LP in Hz AC in  per cent Method vS vP ρ 1 1 1 0 10 10 L-BFGS 2 1 1 1 10 1 L-BFGS 3–9 1 1 1 Increment of 5 1 L-BFGS 10 1 1 1 50 1 L-BFGS 11 1 1 1 55 1 CG 12–25 1 1 1 Increment of 5 1 CG 26 1 1 1 130 0.5 CG Stage Update LP in Hz AC in  per cent Method vS vP ρ 1 1 1 0 10 10 L-BFGS 2 1 1 1 10 1 L-BFGS 3–9 1 1 1 Increment of 5 1 L-BFGS 10 1 1 1 50 1 L-BFGS 11 1 1 1 55 1 CG 12–25 1 1 1 Increment of 5 1 CG 26 1 1 1 130 0.5 CG View Large 4.6 Results In the following, we discuss the result of each inversion individually. In Fig. 8 we show a direct comparison of the final inversion results as well as the initial model, where we focus on an area around the trench. Figs 9(a)–(c) shows a comparison of seismograms of the inversion result, of the observed data set as well as for the initial model. Figure 9. View largeDownload slide Data fit of the field data inversions. Comparison of the seismograms for the initial model, for the final model and of the observed data set for all three inversions (a–c). Evolution of the objective function over the iterations (d). Figure 9. View largeDownload slide Data fit of the field data inversions. Comparison of the seismograms for the initial model, for the final model and of the observed data set for all three inversions (a–c). Evolution of the objective function over the iterations (d). 4.6.1 Individual Love wave FWI The final vS model is still predominantly depth dependent, but contains 2-D lateral variations. At the expected position of the refilled trench the inversion revealed a low-velocity anomaly. The anomaly of the trench has a smooth triangular shape and exhibits a length of 8 m at the surface and a depth of about 2.6 m. To the left of the trench the vS model contains a second low-velocity anomaly, which is elongated and only present in the shallow part. This anomaly could be related to an increased saturation of the shallow soil, which coincides with the observed soil conditions in this area during the measurements. The SH waves are not sensitive to the P-wave velocity, thus, the vP model is not updated. The 2-D variations observed in the vS model are not visible in the ρ model. The inversion increased the density values in both layers and added strong small-scale variations. The inversion improved the fit of the synthetic seismograms, however, some residual signal still remains. The fit of the near offset traces is better than that of the far offset traces. The waves that arrive at the receivers in the far offset travelled farther and deeper, hence, more subsurface heterogeneities influence those recordings than the recordings at near offset receivers. The inversion decreased the misfit relative to the initial misfit up to the 30 Hz frequency stage, as shown in Fig. 9(d). For higher frequency stages the inversion could not again decrease the misfit below the misfit level of the previous frequency stage, most likely due to an increase of the noise level at higher frequencies. In Fig. 10(a) we show the estimated source wavelets for the final model at 130 Hz. The wavelets are quite similar for all source positions, but the source wavelets for the higher source numbers show a ringing effect. This effect could be related to the saturated soil at the end of the profile, which allowed the steel plate used as a source to slightly oscillate. A delay in the arrivals of certain source wavelets could potentially be related to the acquision set-up and possibly an increasing deviation (in distance) of the source locations from the receiver line over increasing offset. Figure 10. View largeDownload slide Estimated source–time functions for the P-SV sources and the SH sources for the final model at 130 Hz of the individual wave type inversions (first column) and of the simultaneous joint FWI (second column). The SH source wavelets, |$v^s_y$|⁠, are shown in the (a) and the P-SV source wavelets, |$v^s_z$|⁠, in (b). Figure 10. View largeDownload slide Estimated source–time functions for the P-SV sources and the SH sources for the final model at 130 Hz of the individual wave type inversions (first column) and of the simultaneous joint FWI (second column). The SH source wavelets, |$v^s_y$|⁠, are shown in the (a) and the P-SV source wavelets, |$v^s_z$|⁠, in (b). 4.6.2 Individual Rayleigh wave FWI The individual Rayleigh wave FWI revealed several 2-D structures in the vS model, which superimpose the mainly lateral homogeneous background model. In the middle of the profile the vS model contains a square shaped low-velocity anomaly, which corresponds to the refilled trench. The lower edge of this anomaly lies in a depth of approximately 2.2 m. However, the contours of the anomaly are blurred, especially at the surface, where we cannot estimate the horizontal length of the anomaly. Again to the left of the trench, the low-velocity anomaly is reconstructed, like in the case of the individual Love wave FWI. In general, the vS model suffers from slight vertically orientated artefacts underneath some source positions, in particular at the positions of the low-velocity anomalies. We observed similar artefacts in the synthetic example at positions with inaccurate P-wave velocities. However, we expect potential artefacts to be less dominant in the case of the field data application than in the case of the synthetic example, since we set the size of the median smoothing filter to 1 m for the field data FWI, which is twice as big as in the synthetic example. The overall variations in the final vP model are small compared to the vS model. At the position of the trench as well as at the position of the second anomaly in the vS model we observe slightly reduced P-wave velocities. As demonstrated in the synthetic example, the vP model could suffer from a cross-talk by vS, hence, the light anomalies in the vP model might be a result of cross-talk. The small amplitudes of the refracted waves were not fitted by the FWI in the presence of high amplitude surface waves. As shown by Athanasopoulos & Bohlen (2017a), further measures need to be taken in order to properly acount for the refracted waves such is the sequential inversion of each wave-type separetely prior to inverting the total wavefield. The final ρ model does not contain any of the anomalies that are present in the vS model. The inversion increased the density values in the first and second layer and added smooth small-scale vertically orientated artefacts. The final seismogram fit is not as accurate as in the case of the individual Love wave FWI, however, the P-SV wavefield is more complex than the SH wavefield. The inversion fitted the fundamental Rayleigh mode quite well at all offsets. Generally, the fit of the phases is better than of the amplitudes. The evolution of the objective function is shown in Fig. 9(d). The inversion reduced the misfit relative to the initial misfit up to the frequency stage of 30 Hz. From 30 Hz on the inversion only reduced the misfit within each frequency stage. The source time inversion revealed similar wavelets for all source positions, as indicated in Fig. 10(b). 4.6.3 Simultaneous Joint FWI The final vS model of the simultaneous joint FWI contains a low-velocity anomaly at the expected position of the refilled trench. The anomaly has an identical size as in the case of the individual Love wave FWI. The shape of the anomaly is triangular and the contour is sharp. The anomaly of the trench holds higher velocities in the shallow part than in the lower part. Moreover, the final vS model contains an elongated shallow second low-velocity anomaly to the left of the trench, which could be related to an increased saturation of the shallow soil within this area. Altogether, the vS model is still predominantly depth dependent. The variations of the vP model are small compared to the initial model. At the positions of the vS anomalies the vP model shows slightly lower values, which could be a result of a cross-talk. We expect the resolution in the vP model to be low, as explored in the synthetic experiments. The ρ model contains higher values than the initial model, especially in the upper layer. Neither of the two vS anomalies can be observed in the density model. However, the ρ model suffers from vertically orientated small-scale lateral heterogeneities. The fit of the synthetic seismograms to the observed seismograms is not as accurate as in the case of the individual wave type inversions, as opposed to the synthetic experiments, where the joint FWI decreased the misfit of both wave types even further. One reason could be slight anisotropic effects that would influence the propagation of both wave types differently, due to the contrasting polarization of the SH and the P-SV waves. Another reason could be a limited accuracy of the initial vP model. Nonetheless, the fit of the SH data is better than the fit of the P-SV data set. Fig. 9(d) illustrates the evolution of the joint objective function. The joint FWI reached the lowest misfit for the 20 Hz frequency stage. From 20 Hz on the joint FWI decreased the misfit only within each frequency stage and could not decrease it below the misfit of the previous frequency stages. In Figs 10(a) and (b) we present the improved source wavelets for the final frequency stage of 130 Hz. For both wave types the source wavelet estimation revealed homogeneous wavelets across the whole profile, despite a slight ringing effect for the sources with higher source numbers. 4.7 Comparison with ground-penetrating radar To evaluate the reliability of our FWI results, we compare them with the result of a ground-penetrating radar (GPR) measurement. The GPR measurement took place at a later date than the seismic measurements and was carried out in the framework of a master’s thesis by Wegscheider (2017). For the zero-offset GPR measurement they used a radar with a 200 MHz antenna. They time-migrated the data set by a constant-velocity Kirchhoff migration, where they chose the migration velocity to 0.1 m ns–1. We were provided with the final migrated image shown in Fig. 11, where the trench is visible by boundary reflections which reveal a triangular form. Inside the trench few reflections are visible, suggesting a homogeneous filling of the trench. Figure 11. View largeDownload slide Qualitative comparison of the GPR result with the field data FWI results. Overlay of the final S-wave velocity models of individual Love wave FWI (top panel), individual Rayleigh wave FWI (middle panel) and simultaneous joint FWI (bottom panel) with the time-migrated image of the GPR measurement. A colour-bar is not shown, because the transparency effect would falsify the colour representation. Figure 11. View largeDownload slide Qualitative comparison of the GPR result with the field data FWI results. Overlay of the final S-wave velocity models of individual Love wave FWI (top panel), individual Rayleigh wave FWI (middle panel) and simultaneous joint FWI (bottom panel) with the time-migrated image of the GPR measurement. A colour-bar is not shown, because the transparency effect would falsify the colour representation. We used the final S-wave velocity models to compare the FWI results to the GPR image, since the FWI could resolve this parameter with higher resolution and lesser ambiguities compared to the other elastic parameters. We make only a qualitative comparison by linking the location of velocity anomalies to reflections in the GPR image. Due to the difficulties relating an accurate depth-migration of the GPR result, we adjusted the depth axis of the GPR result manually in order to fit the FWI results. We used a velocity factor of 0.086 m ns–1 to transfer the GPR result from the time-domain to the depth-domain of the FWI results. The qualitative comparison of the GPR image to the three FWI velocity models is shown in Fig. 11 as an overlay of both results. The result of the individual Love wave FWI matches the GPR image quite well. The velocity model reveals contrasts at positions where strong reflections at the boundaries of the trench are visible. The extension of the low-velocity anomaly is mainly concentrated to the enclosed part of these boundary reflections. Structures within the trench are not visible in the velocity model. The low-velocity anomaly to the left of the trench correlates with near-surface reflectors, which could represent impermeable layers that lead to accumulated water at the surface. The velocity model of the individual Rayleigh wave FWI satisfactorily matches the migrated GPR image. The centre part of the low-velocity anomaly lies within the enclosed area of the boundary reflections of the trench. However, we could not match the horizontal extension of the trench between both images. The lower part of the trench is not visible in the velocity model. The second low-velocity anomaly to the left of the trench correlates again with shallow reflections. The simultaneous joint FWI revealed a low-velocity anomaly that fills the enclosed part of the boundary reflections of the trench accurately. These boundary reflections are visible as sharp contrast within the velocity model. The lower part of the trench matches between both images. The shallow low-velocity anomaly that is present on the left side correlates again with near-surface reflectors. 5 CONCLUSION To conclude, we propose two recommendations for near-surface investigations of the S-wave velocity distribution by the FWI of shallow-seismic surface waves: In the absence of an accurate initial model for the P-wave velocity, we recommend the individual Love wave FWI for three main reasons: (1) Its convergence behavior is smooth and independent of the P-wave velocity, (2) it therefore holds a smaller parameter space, which leads to less cross-talk effects and (3) the SH wave equation is less complex than the P-SV wave equation, which allows a computationally efficient inversion. In the case an accurate initial model for the P-wave velocity is available, we recommend the individual Love wave FWI against the individual Rayleigh wave FWI for the same reasons. However, in this case a simultaneous joint FWI of both wave types has several advantages compared to both individual wave type inversions: (1) It decreases the ambiguities of the inversion result, since more data is evaluated, (2) it reduces the cross-talk between the elastic parameter classes and (3) it further improves the resolution and accuracy of the inversion result. ACKNOWLEDGEMENTS This work is financially supported by the German Research Foundation (DFG) through CRC 1173 and the German Ministry of Education and Research (BMBF) through the project WAVE, grant 01IH15004A. The authors also appreciate an anonymous reviewer and Weiguang He, for their constructive comments, which greatly improve the quality of this paper. The simulations were performed on the computational resource ForHLR Phase I funded by the Ministry of Science, Research and the Arts Baden-Württemberg and DFG. REFERENCES Athanasopoulos N. , Bohlen T. , 2016 . Sequential full-waveform inversion of refracted and rayleigh waves , in Proceedings of the Near Surface Geoscience 2016-22nd European Meeting of Environmental and Engineering Geophysics , EarthDoc . Athanasopoulos N. , Bohlen T. , 2017a . Aquifer characterization using elastic full-waveform inversion , in Proceedings of the Near Surface Geoscience 2016-23nd European Meeting of Environmental and Engineering Geophysics , EarthDoc . Athanasopoulos N. , Bohlen T. , 2017b . Field data application of sequential full-waveform inversion of refracted and rayleigh waves , in Proceedings of the 79th EAGE Conference and Exhibition 2017 , EarthDoc . Binnig M. , 2015 , Full Waveform Inversion of shallow-seismic Rayleigh waves to characterize the ”Ettlinger Linie” , Master’s thesis , Karlsruher Institut für Technologie (KIT) . Bleibinhaus F. , Hole J.A. , Ryberg T. , Fuis G.S. , 2007 . Structure of the California Coast Ranges and San Andreas Fault at SAFOD from seismic waveform inversion and reflection imaging , J. geophys. Res.: Solid Earth , 112 ( B6 ). https://doi.org/10.1029/2006JB004611 Bohlen T. , 2002 . Parallel 3-D viscoelastic finite difference seismic modelling , Comput. Geosci. , 28 ( 8 ), 887 – 899 .. https://doi.org/10.1016/S0098-3004(02)00006-7 Google Scholar Crossref Search ADS Bohlen T. , Kugler S. , Klein G. , Theilen F. , 2004 . 1.5 D inversion of lateral variation of Scholte-wave dispersion , Geophysics , 69 ( 2 ), 330 – 344 .. https://doi.org/10.1190/1.1707052 Google Scholar Crossref Search ADS Brossier R. , 2011 . Two-dimensional frequency-domain visco-elastic full waveform inversion: Parallel algorithms, optimization and performance , Comput. Geosci. , 37 ( 4 ), 444 – 455 .. https://doi.org/10.1016/j.cageo.2010.09.013 Google Scholar Crossref Search ADS Brossier R. , Operto S. , Virieux J. , 2009 . Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion , Geophysics , 74 ( 6 ), WCC105 – WCC118 .. https://doi.org/10.1190/1.3215771 Google Scholar Crossref Search ADS Dokter E. , Köhn D. , Wilken D. , Rabbel W. , 2014 . Application of elastic 2D waveform inversion to a near surface SH-wave dataset , in Proceedings of the 76th EAGE Conference and Exhibition 2014 , EarthDoc . Dokter E. , Khn D. , Wilken D. , Nil D.D. , Rabbel W. , 2017 . Full waveform inversion of sh and lovewave data in nearsurface prospecting , Geophys. Prospect. , 65 ( S1 ), 216 – 236 .. https://doi.org/10.1111/1365-2478.12549 Google Scholar Crossref Search ADS Dziewonski A. , Bloch S. , Landisman M. , 1969 . A technique for the analysis of transient seismic signals , Bull. seism. Soc. Am. , 59 ( 1 ), 427 Fichtner A. , Trampert J. , 2011 . Hessian kernels of seismic data functionals based upon adjoint techniques , Geophys. J. Int. , 185 , 775 – 798 .. https://doi.org/10.1111/j.1365-246X.2011.04966.x Google Scholar Crossref Search ADS Fichtner A. , Kennett B.L. , Igel H. , Bunge H.-P. , 2009 . Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods , Geophys. J. Int. , 179 ( 3 ), 1703 – 1725 .. https://doi.org/10.1111/j.1365-246X.2009.04368.x Google Scholar Crossref Search ADS Forbriger T. , 2003a . Inversion of shallow-seismic wavefields: I. Wavefield transformation , Geophys. J. Int. , 153 ( 3 ), 719 – 734 .. https://doi.org/10.1046/j.1365-246X.2003.01929.x Google Scholar Crossref Search ADS Forbriger T. , 2003b . Inversion of shallow-seismic wavefields: II. Inferring subsurface properties from wavefield transforms , Geophys. J. Int. , 153 ( 3 ), 735 – 752 .. https://doi.org/10.1046/j.1365-246X.2003.01985.x Google Scholar Crossref Search ADS Forbriger T. , Groos L. , Schäfer M. , 2014 . Line-source simulation for shallow-seismic data. Part 1: theoretical background , Geophys. J. Int. , 198 ( 3 ), 1387 – 1404 .. https://doi.org/10.1093/gji/ggu199 Google Scholar Crossref Search ADS Gardner G. , Gardner L. , Gregory A. , 1974 . Formation velocity and density-the diagnostic basics for stratigraphic traps , Geophysics , 39 ( 6 ), 770 – 780 .. https://doi.org/10.1190/1.1440465 Google Scholar Crossref Search ADS Gélis C. , Virieux J. , Grandjean G. , 2007 . Two-dimensional elastic full waveform inversion using Born and Rytov formulations in the frequency domain , Geophys. J. Int. , 168 ( 2 ), 605 – 633 .. https://doi.org/10.1111/j.1365-246X.2006.03135.x Google Scholar Crossref Search ADS Groos L. , 2013 . 2D full waveform inversion of shallow-seismic Rayleigh waves , Ph.D. thesis , Karlsruhe, Karlsruher Institut für Technologie (KIT) . Groos L. , Schäfer M. , Forbriger T. , Bohlen T. , 2014 . The role of attenuation in 2D full-waveform inversion of shallow-seismic body and Rayleigh waves , Geophysics , 79 ( 6 ), R247 – R261 .. https://doi.org/10.1190/geo2013-0462.1 Google Scholar Crossref Search ADS Groos L. , Schäfer M. , Forbriger T. , Thomas B. , 2017 . Application of a complete workflow for 2D elastic full-waveform inversion to recorded shallow-seismic Rayleigh waves , Geophysics , 82 ( 2 ), R109 – R117 .. https://doi.org/10.1190/geo2016-0284.1 Google Scholar Crossref Search ADS Hüttner R. , Konrad H.-J. , Zitzmann A. , 1968 . Geologische Übersichtskarte 1:200000, Blatt CC7110 Mannheim , Bundesanstalt für Geowissenschaften und Rohstoffe in Zusammenarbeit mit den Geologischen Landesämtern der Bundesrepublik Deutschland und benachbarter Staaten . Köhn D. , 2011 . Time domain 2D elastic full waveform tomography , Ph.D. thesis , Christian-Albrechts-Universität zu Kiel . Köhn D. , De Nil D. , Kurzmann A. , Przebindowska A. , Bohlen T. , 2012 . On the influence of model parametrization in elastic full waveform tomography , Geophys. J. Int. , 191 ( 1 ), 325 – 345 .. https://doi.org/10.1111/j.1365-246X.2012.05633.x Google Scholar Crossref Search ADS Komatitsch D. , Martin R. , 2007 . An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation , Geophysics , 72 ( 5 ), SM155 – SM167 .. https://doi.org/10.1190/1.2757586 Google Scholar Crossref Search ADS Lailly P. , 1983 . The seismic inverse problem as a sequence of before stack migrations , in Proceedings of the Conference on Inverse Scattering, Theory and Application , Society for Industrial and Applied Mathematics , pp. 206 – 220 . Lang K. , 1907 . Die Ettlinger Linien und ihre Geschichte , Veröffentlichungen des Karlsruher Altertumsvereins . Lay T. , Wallace T.C. , 1995 . Modern Global Seismology , Vol. 58 , Academic Press . Levander A.R. , 1988 . Fourth-order finite-difference P-SV seismograms , Geophysics , 53 ( 11 ), 1425 – 1436 .. https://doi.org/10.1190/1.1442422 Google Scholar Crossref Search ADS Lüttschwager G. , 2014 . Simulation and surveying of the near field radiation of seismic vibration sources , Master’s thesis , Karlsruher Institut für Technologie (KIT) . Manukyan E. , Latzel S. , Maurer H. , Marelli S. , Greenhalgh S.A. , 2012 . Exploitation of data-information content in elastic-waveform inversions , Geophysics , 77 ( 2 ), R105 – R115 . Google Scholar Crossref Search ADS Manukyan E. , Maurer H. , Nuber A. , 2018 . Improvements to elastic full-waveform inversion using cross-gradient constraints , Geophysics , 83 ( 2 ), R105 – R115 .. https://doi.org/10.1190/geo2017-0266.1 Google Scholar Crossref Search ADS Maurer H. , Greenhalgh S.A. , Manukyan E. , Marelli S. , Green A.G. , 2012 . Receiver-coupling effects in seismic waveform inversions , Geophysics , 77 ( 1 ), R57 – R63 .. https://doi.org/10.1190/geo2010-0402.1 Google Scholar Crossref Search ADS Maurer H. , Nuber A. , Martiartu N.K. , Reiser F. , Boehm C. , Manukyan E. , Schmelzbach C. , Fichtner A. , 2017 . Chapter one - optimized experimental design in the context of seismic full waveform inversion and seismic waveform imaging , Adv. Geophys. , 58 , 1 – 45 . Google Scholar Crossref Search ADS Mora P. , 1987 . Nonlinear two-dimensional elastic inversion of multioffset seismic data , Geophysics , 52 ( 9 ), 1211 – 1228 .. https://doi.org/10.1190/1.1442384 Google Scholar Crossref Search ADS Nocedal J. , Wright S. , 2006 . Numerical Optimization , Springer Science & Business Media . Operto S. , Ravaut C. , Improta L. , Virieux J. , Herrero A. , Dell’Aversana P. , 2004 . Quantitative imaging of complex structures from dense wide-aperture seismic data by multiscale traveltime and waveform inversions: a case study , Geophys. Prospect. , 52 ( 6 ), 625 – 651 .. https://doi.org/10.1111/j.1365-2478.2004.00452.x Google Scholar Crossref Search ADS Pan Y. , Xia J. , Xu Y. , Gao L. , Xu Z. , 2016 . Love-wave waveform inversion in time domain for shallow shear-wave velocity , Geophysics , 81 ( 1 ), R1 – R14 . Google Scholar Crossref Search ADS Pan Y. , Gao L. , Bohlen T. , 2018 . Time-domain full-waveform inversion of Rayleigh and Love waves in presence of free-surface topography , J. appl. Geophys. , 152 , 77 – 85 .. https://doi.org/10.1016/j.jappgeo.2018.03.006 Google Scholar Crossref Search ADS Plessix R.-E. , Mulder W. , 2004 . Frequency-domain finite-difference amplitude-preserving migration , Geophys. J. Int. , 157 ( 3 ), 975 – 987 .. https://doi.org/10.1111/j.1365-246X.2004.02282.x Google Scholar Crossref Search ADS Pratt R.G. , 1999 . Seismic waveform inversion in the frequency domain, part 1: theory and verification in a physical scale model , Geophysics , 64 ( 3 ), 888 – 901 .. https://doi.org/10.1190/1.1444597 Google Scholar Crossref Search ADS Romdhane A. , Grandjean G. , Brossier R. , Rejiba F. , Operto S. , Virieux J. , 2011 . Shallow-structure characterization by 2D elastic full-waveform inversion , Geophysics , 76 ( 3 ), R81 – R93 .. https://doi.org/10.1190/1.3569798 Google Scholar Crossref Search ADS Safani J. , O’Neill A. , Matsuoka T. , Sanada Y. , 2005 . Applications of Love wave dispersion for improved shear-wave velocity imaging , J. Environ. Eng. Geophys. , 10 ( 2 ), 135 – 150 .. https://doi.org/10.2113/JEEG10.2.135 Google Scholar Crossref Search ADS Schaefer M. , 2014 . Application of full-waveform inversion to shallow-seismic Rayleigh waves on 2D structures , Ph.D. thesis , Karlsruhe, Karlsruher Institut für Technologie (KIT) . Socco L.V. , Foti S. , Boiero D. , 2010 . Surface-wave analysis for building near-surface velocity models—established approaches and new perspectives , Geophysics , 75 ( 5 ), 75A83 – 75A102 .. https://doi.org/10.1190/1.3479491 Google Scholar Crossref Search ADS Tarantola A. , 1984 . Inversion of seismic reflection data in the acoustic approximation , Geophysics , 49 ( 8 ), 1259 – 1266 .. https://doi.org/10.1190/1.1441754 Google Scholar Crossref Search ADS Tarantola A. , 1986 . A strategy for nonlinear elastic inversion of seismic reflection data , Geophysics , 51 ( 10 ), 1893 – 1903 .. https://doi.org/10.1190/1.1442046 Google Scholar Crossref Search ADS Tran K.T. , McVay M. , Faraone M. , Horhota D. , 2013 . Sinkhole detection using 2D full seismic waveform tomography , Geophysics , 78 ( 5 ), R175 – R183 .. https://doi.org/10.1190/geo2013-0063.1 Google Scholar Crossref Search ADS Virieux J. , 1984 . SH-wave propagation in heterogeneous media: Velocity-stress finite-difference method , Geophysics , 49 ( 11 ), 1933 – 1942 .. https://doi.org/10.1190/1.1441605 Google Scholar Crossref Search ADS Virieux J. , 1986 . P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method , Geophysics , 51 ( 4 ), 889 – 901 .. https://doi.org/10.1190/1.1442147 Google Scholar Crossref Search ADS Virieux J. , Operto S. , 2009 . An overview of full-waveform inversion in exploration geophysics , Geophysics , 74 ( 6 ), WCC1 – WCC26 .. https://doi.org/10.1190/1.3238367 Google Scholar Crossref Search ADS Wegscheider S. , 2017 , Geophysikalische Untersuchung der Ettlinger Linie auf dem Segelflugplatz in Rheinstetten , Master’s thesis , Karlsruher Institut für Technologie (KIT) . Wehner D. , Köhn D. , De Nil D. , Schmidt S. , al Hagrey S. , Rabbel W. , 2015 . A combined elastic waveform and gravity inversion for improved density model resolution applied to the Marmousi-II model , in Proceedings of the 77th EAGE Conference and Exhibition 2015 , EarthDoc . Wittkamp F. , 2016 , Individual and joint 2-D elastic full-waveform inversion of Rayleigh and Love waves , Master’s thesis , Karlsruhe Institute of Technology(KIT) , http://dx.doi.org/10.5445/IR/1000058955. Xia J. , Xu Y. , Luo Y. , Miller R.D. , Cakir R. , Zeng C. , 2012 . Advantages of using multichannel analysis of Love waves (MALW) to estimate near-surface shear-wave velocity , Surv. Geophys. , 33 ( 5 ), 841 – 860 .. https://doi.org/10.1007/s10712-012-9174-2 Google Scholar Crossref Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Individual and joint 2-D elastic full-waveform inversion of Rayleigh and Love waves
JF - Geophysical Journal International
DO - 10.1093/gji/ggy432
DA - 2019-01-01
UR - https://www.deepdyve.com/lp/oxford-university-press/individual-and-joint-2-d-elastic-full-waveform-inversion-of-rayleigh-N600rZYT6t
SP - 350
VL - 216
IS - 1
DP - DeepDyve
ER -