Elastic SH- and Love-wave Full-Waveform Inversion for shallow shear wave velocity with a preconditioned technique

Journal Pre-proof Elastic SH- and Love-wave Full-Waveform Inversion for shallow shear wave velocity with a preconditioned technique

Yingwei Yan, Zhejiang Wang, Jing Li, Nan Huai, Yuntao Liang, Shuanglin Song, Jianmin Zhang, Ling Zhang PII:

S0926-9851(18)30997-2

DOI:

https://doi.org/10.1016/j.jappgeo.2020.103947

Reference:

APPGEO 103947

To appear in:

Journal of Applied Geophysics

Received date:

23 November 2018

Revised date:

20 October 2019

Accepted date:

13 January 2020

Please cite this article as: Y. Yan, Z. Wang, J. Li, et al., Elastic SH- and Love-wave FullWaveform Inversion for shallow shear wave velocity with a preconditioned technique, Journal of Applied Geophysics(2020), https://doi.org/10.1016/j.jappgeo.2020.103947

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier.

Journal Pre-proof

Elastic SH- and Love-wave Full-Waveform Inversion for shallow shear wave velocity with a preconditioned technique Yingwei Yana,c, Zhejiang Wanga,* , Jing Lia , Nan Huaia , Yuntao Liangb, Shuanglin Songb, Jianmin Zhanga , Ling Zhanga College of Geo-Exploration Science and Technology, Jilin University, Changchun, China

b c

State Key laboratory of Coal Mine Safety Technology, CCTEG Shenyang Research Institute, Fushun, China

Academy of Opto-Electronics, China Electronic Technology Group Corporation (AOE CETC), Tianjin, China

*

Corresponding author.

of

a

ro

ABSTRACT

In contrast to the Rayleigh waves, the system of SH- and Love-wave Full-Waveform Inversion

-p

(FWI) is simpler (Only S-wave velocity and density is inverted). However, the gradient is not scaled with the increasing depth in gradient-based FWI. The main reasons for the phenomenon

re

are the limited frequency band of the wavelet, non- uniform coverage of the excitation-receiving

lP

points and the double scattering. The illumination can be enhanced by adding a preconditioned operator to the inversion. To compensate the illumination, we presented a preconditioned

na

technique abiding by the theory of inverse scattering in SH- and Love-wave FWI. The synthetic tests demonstrate that the preconditioned technique can expedite the convergence rate and

ur

enhance the illumination intensity to the deeper medium compared with the conventional gradient-based FWI. The results also indicate that the preconditioned conjugate gradient (PCG)

Jo

method has the best parameter reconstruction performance among the local optimization algorithms. Subsequently, we adopt the PCG method to inverse the field data acquired in the village of Čachtice, northwestern Slovakia. In addition, a multiscale strategy is employed and the Wiechert-Herglotz method is followed to estimate the initial S-wave velocity model in inversion. Eventually, this paper illustrates the credibility of the inverted result from the perspectives of waveform fitting, common-offset profile, signal similarity, relative error of channel energy and normalized misfit function values. This study primarily seeks to delineate the foundations of historical buildings under the survey line. According to the range of S-wave velocity, we assert that there are no foundation remains in the subsurface. Keywords: Full-Waveform Inversion; SH-wave; preconditioned technique; Seismic imaging;

Journal Pre-proof

1. Introduction S-wave velocity is counted as one of the most significant seismic parameters in shallow subsurface, and it is also deemed as a crucial parameter for evaluating the dynamic characteristics of shallow interface (Yilmaz et al., 2006). Furthermore, Imai and Tonouchi (1982) demonstrated that the stratum hardness index N (A parameter of the soil mechanism and foundation engineering) is also correlated with S-wave velocity. In addition, the porosity is also related to the S-wave velocity (Uyanık, 2019). Given the large amplitude characteristics of

of

surface waves and its susceptibility to the S-wave velocity, the dispersion-based methods denoted by multichannel analysis of surface waves (MASW) and multichannel analysis of Love

ro

waves (MALW) have already become one of the most popular tools for ascertaining S-wave

-p

velocity (e.g., Xia et al., 1999, 2012).

It has become a convenient approach to obtain 1D S-wave velocity profile by inverting

re

the basic mode dispersion curve of surface waves (e.g., Hering et al., 1995; Misiek, 1996; Misiek et al., 1997; Xia et al., 1999). A multi- mode joint inversion has also been developed (e.g., Xia et

lP

al., 2000, 2003; Luo et al., 2007) to improve the accuracy of surface waves to depict S-wave velocity. However, the method based on the dispersive property of surface waves can only

na

achieve a description of 1D vertical average structure for the subsurface, which limits its horizontal resolution. To overcome the limitation, a 1.5D S-wave velocity estimation method,

ur

called local slant stack, was proposed (e.g., Bohlen et al., 2004; Kugler et al., 2007; Willken and Rabbel, 2012), whereas it only improve the data utilization. It cannot essentially implement a

Jo

description of true 2D or 3D structure for the subsurface. Compared with the dispersion analysis of surface waves, the Full-Waveform Inversion (FWI) based on matching the observed waveforms can describe the true 2D or 3D structure of the media (e.g., Chen et al., 2017; He and Han, 2018; Jeong et al., 2018; He et al., 2019). Previous researchers on FWI laid great attention on body waves (e.g., Mora, 1987; Crase et al., 1990; Vigh et al., 2014). Only fitting the body waves is enough for deep structural imaging, however, for the shallow surface analysis, the energy of surface waves occupies the most of the wavefield (Xia, 2015). Therefore, the surface waves components cannot be ignored in shallow FWI. One has to pay attention to several issues before starting FWI. Including establishing the misfit function, attaining the gradient, constructing the inverse Hessian operator or

Journal Pre-proof preconditioned operator to reckon with the gradient. The source wavelet and initial model estimation are also an indispensable part of the inversion (Zhang et al., 2017). The misfit is generally guided by the norm, such as L1 (Tarantola, 2005) and L2 norm. The misfit guided by L2 norm has become the most common misfit (e.g., Virieux and Operto, 2009; Métivier et al., 2014). Choi and Alkhalifah (2012) recommended to use the global correlation norm (GCN) to construct the misfit due to the different offset range of the predicted and measured marine streamer data. Due to the different geometric diffusion mechanisms and energy scales between the predicted data and the field data, we also use the GCN to construct the misfit when inverting the field data,

of

and the L2 norm is adopted in the synthesized test. The gradient is considered as the cross correlation between the forward and adjoint wavefield and obtained by the first-order adjoint

ro

state method (Plessix, 2006; Castellanos et al., 2011). Pratt (1998) emphasized the significance

-p

of the inverse Hessian matrix and pointed out it is a diagonal dominant matrix for high- frequency approximation. It is also referred to as the preconditioned operator that reckons with the gradient.

re

The huge dimension of the operator makes it difficult to compute, explicitly. Merely some

lP

implicit methods are adopted to estimate the diagonal elements. Shin (2001) estimated it in prestack depth migration abiding by the theory of inverse scattering, and the diagonal elements were successfully estimated by Operto (2006) through the Gauss-Newton method. Another

na

strategy is to estimate the inverse Hessian matrix through differentiating the gradient of the misfit from previous iterations, and the l-BFGS algorithm is the best known (Byrd et al., 1995;

ur

Nocedal and Wright, 2006; Brossier et al., 2009;) among this series.

Jo

The source wavelet is regarded as one of the most crucial factors in FWI, it directly determines the phase matching degree between the predicted and the observed waveforms. There is a basic assumption in seismology, that is, the signal of the seismic trace is understood as the convolution of the source function and the Green's function of source to the receiving point. The assumption can be utilized to estimate the wavelet, which is the source correction filter (Groos, 2013; Groos et al., 2014). In addition to the synthesized tests and theoretical development, there are only a few 2D FWI cases of field data. Tran and Mcvay (2012) have developed a 2D FWI of Rayleigh waves and employed it to detect sinkholes successfully (Tran et al., 2013). Pan et al. (2016) inverted S-wave velocity of the subsurface medium through calculating the gradient of the misfit explicitly, and gave the subsurface S-wave velocity resolution of Love waves. Li et al. (2017) regard the Rayleigh waves as the skeleton of P-SV wavefield to realize the detection of

Journal Pre-proof the Quaternary faults by wave-equation dispersion inversion (WD). Then, Li et al. (2019) applied this method to the Love-wave inversion and indicated that it has a better convergence property due to the completely insensitive to the P-wave. The Love-wave is more susceptible to S-wave velocity than Rayleigh waves (Dokter et al., 2017). We develop a new preconditioned technique in 2D FWI of SH- and Love-wave based on the theory of inverse scattering in this paper. Given that the density is less contributive to the wavefield, merely the S-wave velocity is updated, and the density remains invariable during the inversion. This reduces the complexity of the inversion system and avoids the multi-parameter

of

(v p and v s) crosstalk problem in Rayleigh waves FWI. The synthesized data of the checkerboard and fault model is inverted by the steepest descent (SD), conjugate gradient (CG),

ro

preconditioned steepest descent (PSD) and preconditioned conjugate gradient (PCG) method,

-p

respectively. The synthesized results show that the PCG method has the best performance of parameter reconstruction. Subsequently, the PCG method is adopted to inverse the field data. To

re

reduce the non-linearity of inversion, a multiscale filtering technique in time domain is added to

lP

the system (Bunks et al., 1995). Finally, we provide the five evidences (waveform fitting, common-offset profile, signal similarity, relative error of trace-energy and normalized misfit

na

function values) to support the credibility of the result.

Jo

2.1. Forward problem

ur

2. Theory of 2D SH- and Love-FWI

In 2D isotropic media, the SH-wave elastodynamic equations that satisfy Hooke’s law (Virieux, 1984) are defined as   x, z  1

  x, z  1

  x, z 

v y  x, z , t  t



 xy  x, z , t  t  zy  x, z , t  t

 xy  x, z , t   

x v y  x, z , t  x



 zy  x, z , t  z

   x, z , t 

(1)

v y  x, z , t  z

where v y (x, z, t) represents the particle velocity in the time domain, ρ (x, z) denotes the density, μ (x, z) indicates the shear modulus, σxy (x, z, t) and σzy (x, z, t) refer to the shear stress in the time domain, φ (x, z, t) is the source wavelet. Forward equation (1) can be solved by a high-order

Journal Pre-proof staggered- grid finite-difference method (Virieux., 1984). The free-surface boundary condition in the top part of the model is achieved by the stress-imaging technique (Graves, 1996), the absorption boundary condition in the left, right and bottom of the model are simulated through the two-axis perfectly matched layer technique (Collino and Tsogka, 2001). 2.2. Framework of the local optimization algorithms In FWI, the most classical misfit is defined by the L2 norm of the data residuals:

1 Lu  p   d 2

(2)

2

of

min f  p  

where L denotes a sampling operator to the receiver locations; u(p) is the numerical solution of

ro

the forward equation (1); p represents the model parameter, d indicates the observed data; the

-p

symbol ||·||2 is the L2 norm.

The huge computation scale of FWI makes it difficult to solve by a global optimization or

re

quasi- global optimization algorithm (Métivier et al., 2014; Guo et al., 2017). From a practical

lP

point of view, only some local optimization algorithms can be adopted in FWI. The core of the local optimization algorithm is starting from an initial model p0 and obtaining the final inverted model through minimizing the misfit (Anagaw and Sacchi, 2018). The algorithm framework is

na

presented as (Métivier et al., 2014)

(3)

ur

pk 1  pk   k pk

where p represents the model parameter, which is the S-wave velocity v s in our paper; Δpk is the

Jo

descent direction; γk is the step size, which can be calculated through a line search and threepoint parabolic method (Bonnans et al., 2006; Nocedal and Wright, 2006; Métivier et al., 2014). The local optimization algorithms for solving FWI is divided into gradient-based and preconditioned algorithms in this paper. The preconditioned operator is described by the inverse Hessian matrix of the misfit function. It is regarded as a scaled factor of the gradient, which can achieve a complementary estimation to the deeper media parameters. To shed light on the role played by the preconditioned technique in FWI, we test the imaging ability of the gradient-based and preconditioned algorithms for the checkerboard and fault model, respectively. The steepest descent, conjugate gradient, preconditioned steepest descent and preconditioned conjugate gradient algorithm are described in detail this below. Steepest descent algorithm:

Journal Pre-proof pkSD  f  pk 

(4)

where the descent direction at the kth iteration is calculated as the opposite of the gradient of the misfit. Conjugate gradient method: CG  p0  f  p0   CG CG  pk  f  pk    k pk 1

(5)

Preconditioned steepest descent algorithm:

pkPSD   Pk f  pk 

of

(6)

Preconditioned conjugate gradient method:

(7)

-p

ro

p0PCG   P0f  p0    PCG PCG  pk   Pk f  pk    k pk 1

where Pk is the preconditioned operator; βk is a scalar parameter that calculated from the adjacent

re

two descent directions (Nocedal J and Wright S, 2006). There are different formulas to calculate

lP

βk, and each defines a unique conjugate gradient method. The βk of conjugate gradient method is defined by the following formula in this paper, which is the famous Polak-Ribiere-Polyak (PRP)

na

formula (Polak and Ribiere, 1969).

f  pk  f  pk   f  pk 1   k  T f  pk 1  f  pk 1  T

ur

(8)

In this paper, the expression of scalar parameter βk in preconditioned conjugate gradient

Jo

algorithm is defined by a series of numerical tests.

 Pk f  pk    Pk f  pk   Pk 1f  pk 1   k   T  Pk 1f  pk 1    Pk 1f  pk 1   T

(9)

To improve the convergence rate, the scalar parameter βk can only take non-negative value (Gilbert and Nocedal, 1992).

k  max k ,0

(10)

2.3. Gradient and preconditioned operator In the above algorithmic framework, the key point is to calculate the gradient and preconditioned operator. In FWI, the gradient is usually determined by the first-order adjoint

Journal Pre-proof state method (Plessix, 2006; Castellanos et al., 2011). The gradient is the ze ro-lag crosscorrelation between the forward wavefield and the adjoint wavefield. The adjoint wavefield is solved by the adjoint equation of the forward equation (1). The forward and adjoint equation can also be expressed as a form in which a partial differential operator is multiplied by a wavefield.

 S  p  u    H H  S  p    L  d  Lu  p  

(11)

In the above formula, S(p) is the partial differential operator related to the forward equation (1); u is the numerical solution of the wavefield; φ represents the source wavelet; the superscript H

of

indicates a conjugate transpose operation on an operator, the adjoint equation is also the

ro

conjugate transpose of the forward-propagation operator S(p); λ is the numerical solution of the adjoint wavefield. It should be clearly pointed out here that the source of the adjoint equation is

-p

the time-reversed sequence of the difference between the observed and the predicted data. For a

re

common-shot seismic record, the gradient of the misfit with respect to the model parameter p is represented as:

lP

f  p   R   p S  p  u  p  ,  

w

(12)

Here, ∂ indicates a partial differential operation, R denotes the real part of the cross-correlation

na

between the forward-propagation and the adjoint wavefield. Where (., .)w is the scalar product in w (inner product space) associated with the choice of the norm in formula (2). Merely the S-

ur

wave velocity is inverted, of which the gradient is given as:

Jo

 n  n f 2 n nt   xy * n   zy * n    - 3     zy   vs  vs n1  t xy t 

(13)

where σxy and σzy are the forward-propagation wavefield; σxy* and σzy* are the adjoint wavefield; nt is the total sampling points in the time domain. The gradient is not scaled with the increasing depth in the gradient-based FWI. The main reasons for the gradient being unscaled are the limited frequency band of the wavelet, the insufficient illumination to the subsurface parameters, non-uniform coverage of the excitationreceiving points and the double scattering. This lack of illumination can be clearly predicted by the Hessian matrix. The inverse Hessian matrix is regarded as a deconvolution operator that can achieve the gradient preconditioning. In real case, the time of explicitly computing the Hessian matrix is almost unacceptable. However, we can obtain the pseudo-Hessian matrix by applying

Journal Pre-proof the reciprocity theorem. We refer to the solution theory of pseudo-Hessian matrix in prestack depth migration (Shin et al., 2001) to give the expression of pseudo-Hessian operator in SH- and Love-wave FWI. The theory of suppressing artifacts generated by the double scattering in gradient vector is called inverse scattering theory. Abiding by the theory, H is expressed as:





H ij  p    pi S  p  u  p   p j S  p  u  p  , i, j = 1, …, M   T

(14)

Here, the superscript T represents the transpose, M is the number of parameters to be reconstructed. H in the case of the S-wave velocity is defined as:

  xy n   xy n   zy n   zy n   j j i i   , i, j = 1, …, M   t t t  n 1  t  

n  nt

(15)

ro

of

4 H ij  vs   2 6  vs

The inverse Hessian matrix has the diagonal dominant property for high- frequency

-p

approximation, therefore only the diagonal terms of the Hessian matrix need to be estimated. The first derivative of the stress versus time has been computed in the process of solving the gradient,

re

so the diagonal terms of the Hessian matrix can be obtained by simple summation. However,

lP

because of the fast decrease of the wavefield, very small values might appear on the diagonal elements of the H. Therefore, the inverse cannot be applied directly due to the numerical

na

instabilities, for this, the inversion is stabilized by Marquart-Levenberg regularization. The preconditioned operator Pk at the kth iteration is adjusted as:

ur

 1 Pk  diag   H ii  vs    Ck k 

  

(16)

Jo

where Ck=max(Hii(v s)k); θ denotes the empirical factor to stabilize the inversion and is determined by trials, which is generally adopted as 10-1 to 10-6 . 2.4. Source correction filter Estimating the source wavelet is an essential step in FWI. We utilize the basic assumption that the signal of each channel is the convolution of the source wavelet and the Green's function, and regard the wavelet estimation as a least square optimization problem (Pratt, 1999). The least square optimization problem which is described by minimizing the objective function:

Journal Pre-proof N 1 M

F  cl    dlk  cl glk

2

(17)

l  0 k 1

therein : dlk  sl Gk glk  sl Gk opt

In the above formula, dlk is the complex Fourier expansion coefficient at angular frequency ωl=lΔω for the field record in the time domain at offset rk to the source; glk is the corresponding complex Fourier expansion coefficient of the simulated record; cl represents the filter coefficient; slopt is the optimized source wavelet in the frequency domain that need to be estimated; Gk denotes the Green’s function of the source to the receiving point at offset rk; sl is the wavelet of

of

the simulated record in the frequency domain; N represents the total points of each channel and the corresponding Fourier transform is also N-point transform; M indicates the total number of

ro

seismic traces that involved in the source wavelet estimation. Due to the influence of geometrical

-p

diffusion, the amplitude of each trace decreases with the increasing offset. To ensure that the weights of the seismic traces participating in the seismic wavelet estimation are consistent, the

re

geometrical diffusion compensation is added to the above formula. To prevent a singular value appearing at a certain frequency of the estimated source wavelet, it is necessary to add a damping

lP

term to the objective function. The new objective function with a geometrical diffusion compensation and damping term is expressed as:

na

N 1 M

F  cl ;     f k 2 dlk  cl glk  ME 2 cl 2

2

(18)

l 0 k 1

ur

Where f k = (rk/1 m)k is the compensation term for geometrical diffusion. Using k to compensate

Jo

the power of each trace by analyzing the attenuation law with offset r. For near-source traces, k is generally less than 0.5. Since the global correlation norm is applied in the inversion of field data, the seismic traces need to be normalized with respect to its maximum absolute amplitude. The operation ensures that the weights of the seismic traces participating in the wavelet estimation are consistent, so the compensation term f k is ignored in extracting the source wavelet. A finite ε is used to keep the least-square solution regular by introducing a water level to the average power. Let the first derivative of formula (18) with respect to cl being zero, and the expression of cl is obtained. For M receivers on the free-surface, the filter coefficients cl is calculated via:

Journal Pre-proof M

cl 

f k 1

2 k

glk dlk

M

ME   f k glk 2

2

(19) 2

k 1

Here, glk* is the complex conjugate of glk. The average power of glk is:

E

1 N 1 M 2  fk glk MN l 0 k 1

2

(20)

The wavelet is estimated by applying the inverse Fourier transform to the optimized source wavelet slopt = slcl. This method achieves the estimation of the wavelet by using the filter

of

coefficient cl to filter the wavelet sl of the simulated record, which is also called the source

ro

correction filter (Groos et al., 2014). If the initial model is consistent with the true model, we can estimate the full phase information of the real source wavelet. In general, the method can

-p

estimate more than 90% of the wavelet phase information. So, the method is applied in sub-

re

bands to achieve multi-scale FWI.

lP

3. Synthesized tests

na

3.1. Acquisition geometry and parameters

All operations in the synthesized tests are done serially with a single-core 4.0-GHz CPU.

ur

To ascertain the effectiveness of our method, synthesized data tests are performed on a checkerboard and fault velocity model, respectively. The size of the model is 49.5 m×24.5 m.

Jo

The calculation domain consists of 100×50 grids with dimensions of 0.5 m×0.5 m. The sampling interval is 0.25 ms in the time domain. The source is a 20 Hz Ricker wavelet with a 30 ms delay. The time-order is established at 2 and space-order is given as 10. The layer grid points being perfectly matched are set as 40. The number of shots is 20, the source points are located on the free-surface, the initial source location is set as 1 m with the interval is 2.5 m in x-direction (the positions of source are represented by the red pentagrams in Fig. 1a and Fig. 5a). The wavefield v y is recorded on the free-surface with the total traces are 100 of each common-shot gather. The source wavelet and density are known during the inversion with a density of 2000 kg / m3 . The maximum iterations are set as 30 and the iteration is automatically terminated once the descent value of two adjacent misfit is less than 0.0001. The empirical factors θ for the checkerboard model and fault model are 10-5 and 10-2 , respectively.

Journal Pre-proof To evaluate the inversion results quantitatively, the comparative error of model is introduced, which is defined as:

errp 

ptrue  pinv ptrue

(21)

2

2

where ptrue is the true model parameters, pinv is the inverted model parameters. 3.2. Reconstructed tests of checkerboard model The checkerboard is 10 m × 5 m (as shown in Fig. 1a) in size with the S-wave velocity of

of

each checkerboard as 400 or 600 m/s determined by its position. The S-wave velocities between

ro

400 and 600 m/s are common in near-surface media. The initial S-wave velocity model is a halfspace with the velocity is 500 m/s. Seismic waves are transmitted in the checkerboard model,

-p

which shall trigger the scattering dramatically (Fig. 1b). Prominent scattering occurs when

re

seismic waves propagate to the boundary of the checkerboard cells, as evidenced by the normalized seismograms. The travel-time of the scattered wave is longer than the first arrival

lP

Love-wave (the red lower cases “a” and “b” point out the two components in Fig. 1b). These two groups of waves occupy most of the energy of the wavefield (Fig. 1b). Scattering is a common

na

wave phenomenon in actual seismic wave propagation, such as there is an abnormal body with obvious wave impedance difference under the ground. Although the checkerboard model cannot b)

Jo

a)

ur

exist in the real world, it can evaluate the effectiveness of the inversion algorithm.

Fig. 1 The checkerboard model and the response of SH- and Love-wave. (a) Checkerboard S-wave velocity model, and the red pentagrams indicate the source position. (b) Normalized seismograms of shot 10 with the receiver spacing is 1 m, and the red lower cases “a” and “b” point out the first arrival Love-wave and scattering wave group.

Journal Pre-proof Arising from the lack of illumination, the amplitude of the normalized descent direction at the first iteration (k=0) of gradient-based algorithms in the deeper medium is far less than that of the shallower medium (Fig. 2a). There are some isolated large artifacts that occurred near the free-surface (Fig. 2a), which greatly affects the numerical stability of the inversion. It would cause that the parameters of the deep media are not significantly updated during the iteration. The amplitude of the shallow part in the normalized descent direction of the preconditioned algorithms has the identical order of magnitude compared with that in the deep part (Fig. 2b). The new formed descent direction contains the local curvature information of the misfit due to

of

the preconditioning. As a result, the deep illumination is enhanced by applying the technique. It accordingly bespeaks that the preconditioned technique can compensate the illumination and

ro

attenuate the artifacts generated by the double scattering and coupling between the source,

-p

receivers and free-surface. The gradient-based algorithms can reconstruct the checkerboard shape of the shallow part, while the shape in the deeper part is not reconstruted (Fig. 3a and 3b).

re

The preconditioned algorithms can reconstruct the checkerboard shape in the shallow and deep

lP

(Fig. 3c and 3d), clearly. The minimum of the S-wave velocity inversion results is less than that of the real S-wave velocity of 400 m/s, which is referred to as the parameter insufficient estimation. The maximum of the S-wave velocity is more than 600 m/s, which is known as the

na

parameter over estimation in this paper. We think the phenomenon is caused by boundary effects and excessive illumination at the free-surface. The reconstructed accuracy might be improved by

ur

adding a priori model constraints to the misfit. In the previous iterations, the steepest descent and

Jo

conjugate gradient methods have the identical value of misfit, and the preconditioned algorithms also take on such features, which shows that convergence sequence does not have conjugacy (β=0) (Fig. 4a). The convergence rate of the preconditioned algorithms is faster than that of the gradient-based algorithms, and the misfit function values in the iterative sequence are always lower than that of the gradient-based algorithms (Fig. 4a). The comparative error in the iterative sequence of the preconditioned algorithms is always less than that in the gradient-based algorithms (Fig. 4b). There are different degrees of insufficient parameter estimation among the four local algorithms (Fig. 4c). The predicted waveforms computed from the inversion results of the four optimization algorithms are matched the primary events that occupy the huge energy of the single trace signal well (Fig. 4d). Yet there are still some differences between the inversion results and the real model, which stems from the short cycle skipping problem. Besides, the

Journal Pre-proof initial model is not within the neighborhood of the global minimum. The preconditioned conjugate gradient method has the smallest comparative error value, which is 0.0930 (Table 1). Although the iterative time of the preconditioned conjugate gradient method is more than 1/4 the other three local algorithms, the required iterations of the preconditioned conjugate algorithm are only 26 (Table 1). Accordingly, factoring in the trade-off between the computation time and inversion accuracy, the preconditioned conjugate gradient method is recommended for the inversion in this paper. b)

lP

re

-p

ro

of

a)

Jo

ur

na

Fig. 2 Comparison of the normalized descent direction of the preconditioned and gradient-based algorithm in the checkerboard model test. (a) The normalized descent direction in the first iteration (k=0) of gradient-based algorithms. (b) The normalized descent direction in the first iteration (k=0) of preconditioned algorithms. a) b)

c)

d)

Journal Pre-proof

d)

Jo

ur

c)

na

lP

re

-p

ro

of

Fig. 3 Comparison of inverted results by local optimization algorithms in the checkerboard model test. (a) The result of the steepest descent algorithm. (b) The result of the conjugate gradient method. (c) The result of the preconditioned steepest descent algorithm. (d) The result of the preconditioned conjugate gradient method. a) b)

Fig. 4 The evaluated curves for checkerboard model test. SD, CG, PSD and PCG represent the steepest descent, the conjugate gradient, the preconditioned steepest descent and the preconditioned conjugate gradient, respectively. (a) The normalized misfit function values plotted against the iteration number. (b) The comparative error curve during the iteration. (c) Vertical profile at x = 45 m for the true and initial model in comparison with inversion results whereby SD, CG, PSD and PCG. (d) Single trace waveform fit at x=4.5 m of shot 10 for the observed and initial waveform in comparison with that of SD, CG, PSD and PCG.

Journal Pre-proof Table 1 Comparative results in performance of four local optimization algorithms in checkerboard model test. Local Optimization Algorithm SD CG PSD PCG

Computation time/iteration (s)

Overall iteration

253.4 262.4 262.7 319.5

25 22 30 26

Normalized misfit function final value 0.0228 0.0207 0.0070 0.0071

Comparative error final value 0.1199 0.1182 0.0945 0.0930

of

3.3. Reconstructed tests of fault model The velocities above and below the fault line are 400 m/s and 600 m/s, respectively (Fig.

ro

5a). The initial model is determined as a velocity- increasing layered model with a priori. The

-p

initial model is divided into 25 layers with the velocity of the first layer is 400 m/s. It can be determined through searching the velocity gradient between layers to minimize the misfit (Bunks

re

et al., 1995). Through searching, the velocity gradient is identified as 8 m/s. We can observe the distinct reflected wave (the red lowercase “a” points out the event)

lP

from the left fault interface, the reflected event generated by the right interface is mixed with the Love-wave event (Fig. 5b) due to the too small traveltime difference. The parameters of deep

na

media have been compensated of the preconditioned algorithms compared with the gradientbased algorithms (Fig. 6a, 6b, 6c and 6d). The fault interface of the inverted results is clearer due

ur

to the addition of the preconditioned technique (F ig. 6a, 6b, 6c and 6d). The reconstructed results

Jo

ascertain the preconditioned technique is significant in the inversion system, again. The final misfit in the preconditioned conjugate gradient method reaches 0.0029, which is less than 1/10 of that in the conjugate gradient method (Table 2 and Fig. 7b). The convergence rate of the preconditioned algorithms is also significantly faster than the above tests (Fig. 4a and 7a). In short, the advantages of using the preconditioned technique here are more evident than the foregoing synthesized test, which could be caused by the more precise initial model in the reconstructed tests. This also proves the significance of a high-precision initial model in FWI. a)

b)

Journal Pre-proof

d)

Jo

ur

c)

na

lP

re

-p

ro

of

Fig. 5 The fault model and the response of SH- and Love-wave. (a) Fault S-wave velocity model, and the source positions are shown by the red pentagrams. (b) Normalized seismograms of shot 10 with the receiver spacing is 1 m, and the red lower case “a” indicates the reflected wave event from the left fault interface. a) b)

Fig. 6 Comparison of inverted results by local optimization algorithms in the fault model test. (a) The result of the steepest descent algorithm. (b) The result of the conjugate gradient method. (c) The result of the preconditioned steepest descent algorithm. (d) The result of the preconditioned conjugate gradient method. a) b)

Journal Pre-proof

iteration number. (b) The comparative error during the iteration.

of

Fig. 7 The evaluated curves for fault model test. (a) The normalized misfit function values plotted against

ro

Table 2 Comparative results in performance of four local optimization algorithms in fault model Overall iteration

171.8 170.9 172.9 183.6

29 27 30 30

lP

4. A field data case

-p

Computation time/iteration (s)

Normalized misfit function final value 0.0437 0.0373 0.0042 0.0029

Comparative error final value 0.0684 0.0652 0.0403 0.0315

na

Local Optimization Algorithm SD CG PSD PCG

re

test.

ur

4.1. SH- and Love-wave survey

Jo

The backyard of a winery in the village of Čachtice, northwestern Slovakia is studied in this paper (Dokter et al., 2017). This survey primarily aims to delineate the foundations and ruins of the historical buildings under the ground. There are 42 geophones placed along the line in this measurement. The first receiver position is (0 m, 0 m) in coordinates with the receiver spacing is 0.5 m. The recording time is 1.2 s with the interval is 0.125 ms. The source is excited in the midway of each receiver pair with the first source deployed on (0.25 m, 0 m) in coordinates. The source spacing is 0.5 m and there are 42 sources excited in the survey. The overall length of the profile is 20.75 m. The SH- and Love-wave is excited by a lateral hammer source. For each record, the receiver receives the signal of opposite polarity (There are two sensors at angle of 180° inside the receiver), then stacks them in the opposite sign to remove the P-wave reflection from the seismograms. The direct wave event is not a straight line, which is the first indication of

Journal Pre-proof the heterogeneity of the subsurface (Fig. 8a and 8b). Dispersive Love-waves and some boundary reflections from the surrounding walls of the survey area are observed from the seismograms (Fig. 8a and 8b). The Love-wave energy is mainly distributed in 5~80 Hz with the dominant band is 10~40 Hz (Fig. 8c and 8d). We can also find faint reflections and diffractions due to the heterogeneity of the subsurface medium (Fig. 8a and 8b). b)

re

-p

ro

of

a)

d)

Jo

ur

na

lP

c)

Fig. 8 Normalized seismograms and average amplitude spectrum of the SH- and Love-wave survey. a – direct wave, b – refracted wave, c – Love-wave, d – diffracted/reflected waves, e – Love-wave reflection in the normalized seismograms. (a) Normalized seismogram of shot 1. (b) Normalized seismogram of shot 42. (c) Average amplitude spectrum of shot 1. (d) Average amplitude spectrum of shot 42.

4.2. Data pre-processing for 2D SH- and Love-wave FWI By virtue of the high signal- to-noise ratio (SNR) of the field data, merely a few simple data pre-processing steps are adopted before the 2D SH- and Love-wave FWI. The data collected in the real world is excited by a point source in 3D space and the source of the 2D forward program used for inversion is equivalent to a line source in 3D space. There is π/4 phase difference between the both types of seismic records due to the different excitation mode of the

Journal Pre-proof source (Forbriger et al., 2014; Schäfer, 2014). In this regard, the observed data need to be transformed to match the 2D forward program. Therefore, we adopt an offset-dependent method of data transformation from 3D to 2D space and spreading correction to achieve the data preprocessing (Forbriger et al., 2014; Schäfer, 2014). The field data pre-processing process is elucidated below: 

To reduce the computational cost, each channel is resampled with the interval 0.25 ms. The recording time is 1.2 s of each channel, of which 1 s is used during the inversion. To remove the low frequency noise from deep formations in the record, a high-pass

of



filter with a corner frequency of 5 Hz is applied to the field data. Convolve each channel data with the convolution factor



For near source channels (| r | ≤ r1) multiply each channel by

t 1

.

-p

ro



phase velocity.

, where v ph is the

2 r v ph

For far source channels (| r | ≥r2) multiply each channel by



Delay each channel for 0.12 s to avoid non-causal parts in the estimated wavelet.



To use the GCN in the misfit, the amplitude is normalized with respect to its

r

2 t 1

.

lP

re



na

maximum absolute amplitude of each seismic trace. For this field data, we choose r1=0.1 m, r2=6 m and v ph =164 m/s in this paper.

ur

4.3. Two-dimensional SH- and Love-wave FWI

Jo

To avoid the local minima in the inversion, a multiscale strategy (Bunks et al., 1995) is added to the inversion. And the preconditioned conjugate gradient method is adopted to invert the field data. The main energy band of the data is 5 to 80 Hz (Fig. 8c and 8d). In this regard, the entire inversion is split into 4 scales, which are 0 to 15Hz, 0 to 20 Hz, 0 to 40Hz and 0 to 80 Hz. As mentioned above, the high-precision initial model can improve the imaging accuracy. For this, the Wiechert-Herglotz (Dokter et al., 2017) method is used to estimate the initial layered S-wave velocity model. As the angle of the acquisition method is limited, the maximum depth of the initial S-wave velocity model established by the method is merely 8 m. Even so, we want to remove the influence of perfectly matched layers at the bottom of the model on the Love-wave dispersion. So, the inversion depth is set as 10 m and the initial S-wave velocity model is properly extrapolated to a depth of 10 m (Fig. 9a) through cubic spline interpolation (the

Journal Pre-proof interpolation method avoids the generation of a velocity discontinuity zone.). Below, we only give the inverted results down to a depth of 8 m considering the angle of receiving the wavefield. For the initial density model, a 1D positive gradient model (Schön, 1996) is built empirically (Fig. 9b). b)

-p

ro

of

a)

Fig. 9. Initial model in the 2D SH- and Love-wave FWI. (a) Initial 1D S-wave velocity model. (b) Initial 1D

re

density model.

lP

The density remains invariable during the inversion. The interval of spatial grid is set as 0.25 m and the sampling interval is 0.25 ms of the inversion. From a series trials, the empirical factor θ is determined for different sub-bands, which are 5.0×10-3 (0-15 Hz), 10-3 (0-20 Hz) and

na

10-4 (0-40 Hz, 0-80 Hz). The valid calculation domain consists of 84 × 41 grids with 20

ur

absorption points in the left, right and bottom of the model. For each common-shot data, the source wavelet is estimated according to the source correction filter reviewed in section 2.4. To

Jo

ensure that the fitting of Love-wave dispersion information between the predicted data and field data is caused by the update of the model, rather than the effect of the far-offset traces, we only utilize the near source seismic traces (offset ≤ 5.0 m) to estimate the source wavelet. To more clearly compare the phase of the source wavelets between the common-shot gathers, only the estimated results in 0.5 seconds are shown in Fig. 10a and 10b. The phase of the wavelets is almost identical, which also indicates that the data has a high SNR (Fig. 10a and 10b). It is difficult to ensure that the wavelets between the shot gathers excited by the hammer source have a high phase consistency. It is elaborated in the discussion section whether the inconsistency of the wavelets could reduce the seismic illumination intensity. a)

b)

Journal Pre-proof

of

Fig. 10. The normalized seismic wavelets estimated by the source correction filter, we give the results of five shot gathers which are shot 1, shot 2, shot 3, shot 4 and shot 5. (a) The estimated seismic wavelets of 0-15 Hz. (b) The estimated seismic wavelets of 0-80 Hz.

ro

The low-frequency information in the normalized seismograms is matched accurately,

except for the noisy trace identified by the black arrow (Fig. 11a and 11b). The noisy trace in the

-p

low- frequency range is opposite in polarity to the predicted data (Fig. 11a and 11b), and the polarity becomes the same after the high frequency information is added (Fig. 11c and 11d). It is

re

speculated that the signal received by the receiver at the position of the noisy trace is affected by

lP

a fixed low frequency noise. And the low frequency noise is highlighted during the application of the multiscale strategy. There are very few seismic traces with significant noise in the entire field

na

data, so these traces are maintained during the inversion without any excessive art ificial operations. The normalized seismograms of the right-side of the source is more accurately fitted

ur

compared with that of the left-side (Fig 11.c and 11.d). The reason of the phenomenon is that the signals in the normalized seismograms of the left-side of the source contain the information from

Jo

the foundation walls at the front of the survey line (the black arrows indicate the components in Figure 11c and 11d), and this part cannot be predicted by the model update. In addition to the waveform fitting, the credibility of the inverted result from the perspective of the common-offset profile, signal similarity, relative error of trace-energy and normalized misfit function values is also analyzed. The measurement of the signal similarity is determined by the global-correlation norm. The formula is:



a(t ), b(t ) a(t ), a(t ) b(t ), b(t )

therein : a(t ), b(t )   dta(t )(b(t )

(22)

Journal Pre-proof where a(t) and b(t) are the observed and predicted signal, and is the dot product. The formula implements a description of the similarity of the two signals with the coefficient between 1 and -1. The difference of energy levels between two seismic traces is described by the relative error of trace energy, it is defined as: errE 

Eobs  E pre Eobs

N

therein : E   vi 2

(23)

i 1

In the above formula, the v i represents an individual sample in the seismic trace; Eobs and Epre are the energy of the observed and predicted trace, respectively; Eobs and Epre follow the same

of

definition as E. The results of analyzing the credibility of the inversion results from these perspectives are given below. Hereinafter, we give the conclusion from these perspectives.

ro

The common-offset gather is extracted from all common-shot data with the offset is 2.75

-p

m (Fig 12.a). The event of common-offset profile is generally an indication of the stratigraphic interface, and the dashed black line depicts the interface in the observed and predicted common-

re

offset profile (Fig 12.a and 12.b). The fluctuations of the observed and predicted interface are

lP

relatively consistent, especially in the transition zone (the area indicated by the black arrows in Fig 12.a and 12.b). The observed and predicted seismic traces of shot 22 have strong signal similarity, except for several traces on the left-side of the source (Fig. 12c). As above, the signals

na

contain the reflected wave from the foundation walls. The relative error of trace-energy for shot 22 has a magnitude of less than 0.5 (Fig 12.d) and the final misfit function value drops to about

ur

1/4 of the initial value due to the application of the multiscale strategy (Fig 12.e). Considering

Jo

the influence of the inherent attenuation of the stratum on the phase of the observed waveform, we believe that such a difference between the observed and predicted data is acceptable. a)

b)

Journal Pre-proof c)

d)

Jo

c)

ur

na

lP

re

-p

ro

of

Fig. 11. Normalized seismograms of the observed and predicted waveforms computed from the inverted results of corresponding frequency band. The shot number is 22 with the receiver spacing is 1 m in the image. The black arrow in (a) and (b) indicates the noisy trace near the source. The black arrows in (c) and (d) point out the influence of the foundation walls on the waveform. (a) Normalized waveforms of 0 to 15 Hz. (b) Normalized waveforms of 0 to 20 Hz. (c) Normalized waveforms of 0 to 40 Hz, (d) Normalized waveforms of 0 to 80 Hz. a) b)

e)

d)

Journal Pre-proof

ro

of

Fig. 12. The evidence of supporting the credibility of inverted result. (a) Common-offset gather of the measured data with the offset is 2.75 m. The dashed black line depicts the maximum absolute amplitudes estimated from the traces in the observed common-offset gather. The black arrows indicate the transition parts of the dashed lines. (b) Predicted common-offset gather. (c) Plots of the similarity coefficient between the observed and predicted traces for shot 22. (d) The relative error of the trace-energy for shot 22 plots against the receiver’s position. (e) The normalized misfit function values for each inverted frequency-band plotted versus iteration number.

-p

The inverted result is split into three velocity zones, and the red dashed lines depict the

possible fluctuations of the stratigraphic interface (Fig 13). The identical transition zone (the

re

black arrow indicates the area in Fig 13) as in the observed common-offset profile can be found

lP

in the inverted result, which also shows the credibility of the result. Compared with the initial model, the parameters of the model have also been significantly updated, which also

na

demonstrates the role of the preconditioned technique in parameter complement estimation. The first layer is a low-velocity layer being approximately 3~4.5 m thick in different horizontal

ur

position, and the velocity ranges basically from 80 to 150 m/s. The depth of bottom interface for the second layer is about 5~7 m determined by its horizontal position and the velocity range is

Jo

approximately 180 to 250 m/s. The strong horizontal heterogeneity of the subsurface is revealed from the fluctuant interfaces between the stratums.

Fig. 13. The inverted S-wave velocity result of the SH- and Love-wave survey. The red dashed lines depict the

Journal Pre-proof possible interfaces between the stratums, and the black arrow indicates the transition zone of the interface that is consistent with the observed common-offset profile.

This survey primarily seeks to delineate the foundations of historic buildings under the ground. Geological interpretation of geophysical imaging results is a task that must be treated carefully. We assert the first layer as a low-velocity weathering layer encompassed by loose soil, and the second and third layer are deemed as the sedimentary layer with slight increase in velocity arising from the compaction. This is also consistent with Dokter's previous results. Stone is the most important material for construction of foundations. The S-wave velocity of

of

sandstone is about 2430 m/s (Carmichael, 2017). Considering the velocity range of the inverted results, we assert that no foundations existed under the survey line. This is also consistent with

ro

Dokter's analysis.

-p

5. Discussion

re

During the acquisition of seismic records, it is difficult to guarantee the consistency of

lP

the seismic wavelets between the shot gathers. Although the field data processed in this paper has a high SNR and strong consistency of wavelets, we try to analyze whether the inconsistency of wavelets would reduce the intensity of seismic illumination. For comparison, we perform a

na

synthesized test on the fault model as described above. The only difference from the synthesized test of the fault model above is that the wavelets of each shot are not consistent. Considering

ur

such an extreme situation, the wavelet of the odd-numbered shot is set as a 20 Hz Ricker wavelet

Jo

with a delay time of 20 ms, and the wavelet of the even-numbered shot is a 30 Hz Gaussian first derivative wavelet with a delay time of 50 ms (Fig. 14a). The fault model is reconstructed by the preconditioned conjugate gradient method. As above, we give the normalized descent direction at the first iteration (k=0), but unfortunately the normalized descent direction does not describe the fault interface (Fig. 14b). And it indicates the existence of a low velocity interlayer at the depth 5-10 m, which is seriously inconsistent with the real model. Since the descent direction is too poor, the reconstructed result is not given here. In the inversion process, compared with the above, the misfit drops slowly and the inversion is terminated automatically at the 15th iteration, but the comparative error does not fall (Fig. 14c and 14d). It is indeed a difficult problem to analyze the physical mechanism that forms this phenomenon. The explanation given by us is that the arrival times of the stratum interface are different due to the different wavelets among all

Journal Pre-proof common-shot gathers, causing the gradient not be enhanced in identical phase, and the seismic illumination is not improved constructively. Here, we try to give an empirical conclusion to guide the actual shallow FWI. The wavelets used in the new synthesized test are negatively correlated with a correlation coefficient of -0.23. Through many synthetic tests, we believe that negatively correlated wavelets would reduce the overall seismic illumination intensity, and shot gathers with a wavelet correlation coefficient greater than 0.55 can be added to shallow FWI. This shows that the acquisition of high SNR data and the stability of the source output is extremely important in shallow surface earthquake engineering. Seismic wavelet shaping and b)

na

lP

re

-p

ro

a)

of

FWI based on interference theory may be a kind idea to solve the inconsistency problem.

d)

Jo

ur

c)

Fig. 14. The evaluated results for the fault model test with different seismic wavelet. (a) “Ricker” represents 20 Hz Ricker wavelet with a delay time of 20 ms, “Gauss1diff” represents 30 Hz Gaussian first derivative wavelet with a delay time of 50 ms. (b) The normalized descent direction of the first iteration (k=0). (c) The normalized misfit function values plotted against the iteration number. (d) The comparative error during the iteration.

6. Conclusion

Journal Pre-proof Abiding by the theory of inverse scattering, we developed a novel preconditioned technique for 2D SH- and Love-wave FWI in time domain. The synthetic tests demonstrate that the preconditioned technique can expedite the convergence rate of the misfit and achieve the complemental estimation for deep medium. Then, the preconditioned conjugate gradient method is adopted to invert the field data collected in the village of Čachtice in the northwestern Slovakia. We also add the multi-scale strategy to the inversion system to avoid the local minimum. Eventually, the credibility of the inverted result is analyzed from the perspective of waveform fitting, common-offset profile, signal similarity, relative error of trace-energy and

of

normalized misfit function values. The strong horizontal heterogeneity of the subsurface is revealed by our presented method. The S-wave velocity of the inverted results are too small

ro

relative to that of sandstone. In this regard, we assert that there are no historical building

-p

foundations under the survey line. Although the preconditioned technique is only applied to the elastic and viscoelastic Rayleigh waves.

lP

Acknowledgements

re

FWI of elastic SH- and Love-wave, a similar method can be realized in viscoelastic Love waves,

na

This work is supported by the National Key R&D Program of China (2018YFC0807900, Topic 2), the National Natural Science Foundation of China (51574148, 41874134) and the

ur

Excellent Youth Fund of Jilin Province (20190103142JH). The authors thank Dr. Eva from the Grant Institute, School of Geosciences, University of Edinburgh, who share the measured data

Jo

collected in the village of Čachtice, northwestern Slovakia, and the initial S-wave velocity model established by the Wiechert-Herglotz method with us. We also thank the two reviewers for their constructive comments.

References Anagaw, Y.A. and M.D. Sacchi, 2018, Model parametrization strategies for Newton-based acoustic full waveform inversion: Journal of Applied Geophysics, 157, 23-36, doi: 10.1016/j.jappgeo.2018.06.004. Behura, J., and R. Snieder, 2013, Virtual real source: Source signature estimation using seismic interferometry: Geophysics, 78, no. 5, Q57–Q68, doi: 10.1190/GEO2013-0069.1. Bohlen, T., S. Kugler, G. Klein, and F. Theilen, 2004, 1.5D inversion of lateral variation of Scholtewave dispersion: Geophysics, 69, no. 2, 330-344, doi: 10.1190/1.1707052.

Journal Pre-proof Bonnans, J. F., J. C. Gilbert, C. Lemaréchal, and C. A. Sagastizábal, 2006, Numerical Optimization: Theoretical and Practical Aspects: Springer. Brossier, R., S. Operto, and, J. Virieux, 2009, Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion: Geophysics, 74, no. 6, WCC63-WCC76, doi: 10.1190/1.3215771. Brossier, R., S. Operto, and, J. Virieux, 2010, Which data residual norm for robust elastic frequencydomain full waveform inversion?: Geophysics, 75, no. 3, R37–R46, doi: 10.1190/1.3379323. Bunks, C., F. Saleck, S. Zaleski, and G. Chavent, 1995, Multiscale seismic waveform inversion: Geophysics, 60, no. 5, 1457–1473, doi: 10.1190/1.1443880.

of

Byrd, R.H., P. Lu, and J. Nocedal, 1995, A limited memory algorithm for bound constrained optimization: SIAM Journal on Scientific Computing, 16, 1190–1208, doi: 10.1137/0916069.

ro

Carmichael, R.S., 2017, Handbook of Physical Properties of Rocks, Vol. II: CRC Press, Inc, doi: 10.1201/9780203712085.

-p

Castellanos, C., V. Etienne, G. Hu, S. Operto, R. Brossier, and J. Virieux, 2011, Algorithmic and

Meeting, SEG, Expanded Abstracts, 2793–2798.

re

methodological developments towards full waveform inversion in 3D elastic media: 81st Annual International

lP

Chen J., C. A. Zelt, and P. Jaiswal, 2017, Detecting a known near-surface target through application of frequency-dependent traveltime tomography and full-waveform inversion to P- and SH-wave seismic refraction data: Geophysics, 82, no. 1, R1-R17, doi: 10.1190/geo2016-0085.1.

na

Choi Y., and T. Alk, 2012, Application of multi-source waveform inversion to marine streamer data using the global correlation norm: Geophysical Prospecting, 60, 748-758, doi: 10.1111/j.1365-2478.2012.

ur

01079.x

Collino, F., and C. Tsogka, 2001, Application of the perfectly matched absorbing layer model to the

10.1190/1.1444908.

Jo

linear elastodynamic problem in anisotropic heterogeneous media: Geophysics, 66, no. 1, 294–307, doi:

Crase, E., A. Pica, M. Noble, J. McDonald, and A. Tarantola, 1990, Robust elastic nonlinear waveform inversion: Application to real data: Geophysics, 55, no. 5, 527–538, doi: 10.1190/1.1442864. Dokter, E., D. Köhn, D. Wilken, D. De Nil, and W. Rabbel, 2017, Full-waveform inversion of SH- and Love-wave data in near-surface prospecting: Geophys ical Prospecting, 65, 216-236, doi: 10.1111/13652478.12549. Forbriger, T., L. Groos, and M. Schäfer, 2014, Line-source simulation for shallow-seismic data. Part 1: Theoretical background: Geophysical Journal International, 198, 1387–1404, doi: 10.1093/gji/ggu199. Gilbert, J. C., and J. Nocedal, 1992, Global convergence properties of conjugate gradient methods for optimization: SIAM Journal on Optimization, 2, 21–42. Graves, R. W., 1996, Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences: Bulletin of the Seismological Society of America, 86, no. 4, 1091–1106, doi: 10.1016/0009-

Journal Pre-proof 2541(96)00016-2. Groos, L., 2013, 2D full waveform inversion of shallow seismic Rayleigh waves: Ph.D. thesis, Karlsruhe Institute of Technology. Groos, L., M. Schäfer, T. Forbriger, and T. Bohlen, 2014, The role of attenuation in 2D full-waveform inversion of shallow-seismic body and Rayleigh waves: Geophysics 79, no. 6, R247–R261, doi: 10.1190/geo2013-0462.1. Guo, Q. and T. Alkhalifah, 2017, Elastic reflection-based waveform inversion with a nonlinear approach: Geophysics, 82, no. 6, R309-R321, doi: 10.1190/geo2016-0407.1. He, B., Y. Liu, H. Lu, H. Shi and J. Yi, 2019, Reflection waveform inversion with variable density:

of

Journal of Applied Geophysics, 170, 103827, doi: 10.1016/j.jappgeo.2019.103-827. He, Q. and B. Han, 2018, Accelerating full waveform inversion using HSS solver and limited memory

ro

conjugate gradient method: Journal of Applied Geophysics, 159, 83-92, doi: 10.1016/j.jappgeo.2018.07.022. Hering, A., R. Misiek, A. Gyulai, T. Ormos, M. Dobroka, and L. Dresen, 1995, A joint inversion

43, 135–156, doi: 10.1111/j.1365-2478.1995.tb00128.x.

-p

algorithm to process geoelectric and surface wave seismic data. Part I: Bas ic ideas: Geophysical Prospecting,

re

Imai, T. and K. Tonouchi, 1982, Correlation of N-value with S-wave velocity and shear modulus:

lP

Proceedings of the 2nd European Symposium on Penetration Testing, Amsterdam. A.A. Balkema Publishers, Netherlands, pp, 57-72.

Jeong, W., C. Tsingas and Y. S. Kim, 2018, Full waveform inversion with angle-dependent gradient

10.1016/j.jappgeo.2018.07.016.

na

preconditioning using wavefield decomposition: Journal of Applied Geophysics, 159, 23-31, doi:

ur

Kugler S., T. Bohlen, T. Forbriger, S. Bussat, and G. Klein, 2007, Scholte-wave tomography for shallow-water marine sediments: Geophysical Journal International, 168, 551–570, doi: 10.1111/j.1365-

Jo

246X.2006.03233.x.

Li, J., Z. Feng and G. Schuster, 2017, Wave-equation dispersion inversion: Geophys ical Journal International, 208, no. 3, 1567-1578, doi: 10.1093/gji/ggw465. Li, J., S. Hanafy, Z. Liu and G. Schuster, 2019, Wave-equation dispersion inversion of Love waves, Geophysics, 84, no. 5, R693-R705, doi: 10.1190/GEO2018-0039.1. Luo, Y., J. Xia, J. Liu, Q. Liu, and S. Xu, 2007, Joint inversion of high frequency surface waves with fundamental

and

higher

modes:

Journal

of

Applied

Geophysics,

62,

no.

4,

375–384,

doi:

10.1016/j.jappgeo.2007.02.004. Métivier, L., F. Bretaudeau, R. Brossier, J. Virieux, and S. Operto, 2014, Full waveform inversion and the truncated Newton method: Quantitative imaging of complex subsurface structures: Geophysical Prospecting, 62, 1353–1375, doi: 10.1111/1365-2478.12136. Mis iek, R., 1996, Surface waves: Application to lithostructural interpretation of near-surface layers in the meter and decameter range: Ph.D. thesis, Ruhr-Universität Bochum.

Journal Pre-proof Mis iek R., A. Liebig, A. Gyulai, T. Ormos, M. Dobroka, and L. Dresen, 1997, A joint inversion algorithm to process geoelectric and surface wave seismic data. Part II: applications: Geophys ical Prospecting, 45, no. 1, 65–85, doi: 10.1046/j.1365-2478.1997.3190241.x. Mora, P., 1987, Nonlinear two-dimensional elastic inversion of multi-offset seismic data: Geophysics, 52, no. 9, 1211–1228, doi: 10.1190/1.1442384. Nocedal, J., and S. J. Wright, 2006, Numerical Optimization: Springer. Operto, S., J. Virieux, J. X. Dessa, and G. Pascal, 2006, Crustal imaging from multifold ocean bottom seismometers data by frequency-domain full-waveform tomography: application to the eastern Nankai trough: Journal of Geophysical Research, 111, B09306, doi: 10.1029/2005jb003835.

of

Pan, Y., J. Xia, Y. Xu, L. Gao, and Z. Xu, 2016, Love-wave waveform inversion for shallow shearwave velocity: Geophysics, 81, no. 1, R1–R14, doi: 10.1190/geo2014-0225.1.

ro

Plessix, R.E., 2006, A review of the adjoint-state method for computing the gradient of a functional with geophysical applications: Geophysical Journal International, 167, 495–503, doi: 10.1111/j.1365-

-p

246X.2006.02978.x.

Polak, E. & Ribiére, G. 1969. Note sur la convergence de methodés de directions conjuguées, Revue

re

Francaise Informat. Reserche Opérationnelle, 3e Année, 16, 35-43.

lP

Pratt, R. G., C. Shin, and G. J. Hicks, 1998, Gauss-Newton and full Newton methods in frequencyspace seismic waveform inversion: Geophysical Journal International, 133, 341–362. Pratt, R. G., 1999, Seismic waveform inversion in the frequency domain, Part I: Theory and

na

verification in a physical scale model: Geophysics, 64, no. 3, 888–901, doi: 10.1190/1.1444597. Schäfer, M., 2014, Application of full-waveform inversion to shallow-seismic Rayleigh waves on 2D

ur

structures: Ph.D. thesis, Karlsruhe Institute of Technology. Schön, J., 1996, Physical Properties of Rocks: Fundamentals and Principles of Petrophysics:

Jo

Pergamon.

Shin, C., S. Jang, and D. J. Min, 2001, Improved amplitude preservation for prestack depth migration by inverse scattering theory: Geophysical Prospecting, 49, 592–606, doi: 10.1046/j.1365-2478.2001.00279.x. Tarantola, A., 2005, Inverse problem theory and methods for model parameter estimation: SIAM, doi: 10.1137/1.9780898717921. Tran, K. T., and M. McVay, 2012, Site characterization using Gauss-Newton inversion of 2D full seismic waveform in the time domain: Soil Dynamics and Earthquake Engineering, 43, 16–24, doi: 10.1016/j.soildyn.2012.07.004. Tran, K. T., M. McVay, M. Faraone, and D. Horhota, 2013, Sinkhole detection using 2D full seismic waveform tomography: Geophysics, 78, no. 5, R175–R183, doi: 10.1190/geo2013-0063.1. Uyanık, O., 2019, Estimation of the porosity of clay soils using seismic P- and S-wave velocities: Journal of Applied Geophysics, 170, 103832, doi: 10.1016/j.jappgeo.2019.103832. Vigh, D., K. Jiao, D. Watts, and D. Sun, 2014, Elastic full-waveform inversion application using

Journal Pre-proof multicomponent measurements of seismic data collection: Geophysics, 79, no. 2, R63–R77, doi: 10.1190/geo2013-0055.1. Virieux, J., 1984, SH wave propagation in heterogeneous media: Velocity stress finite-difference method: Geophysics, 49, 1933–1942. Virieux, J., and S. Operto, 2009, An overview of full-waveform inversion in exploration geophysics: Geophysics, 74, no. 6, WCC127–WCC152, doi: 10.1190/1.3238367. Wilken, D., and W. Rabbel, 2012, On the application of particle swarm optimization strategies on Scholte-wave inversion: Geophysical Journal International, 190, no. 1, 580–594, doi: 10.1111/j.1365246x.2012.05500.x.

of

Xia, J., R. D. Miller, and C. B. Park, 1999, Estimation of near-surface shear wave velocity by inversion of Rayleigh wave: Geophysics, 64, 691–700, doi: 10.1190/1.1444578.

ro

Xia, J., R. D. Miller, and C. B. Park, 2000, Advantages of calculating shear-wave velocity from surface waves with higher modes: 70th Annual International Meeting, SEG, Expanded Abstracts, 1295–1298.

-p

Xia, J., R. D. Miller, C. B. Park, and G. Tian, 2003, Inversion of high frequency surface waves with fundamental and higher modes: Journal of Applied Geophysics, 52, no. 1, 45–57, doi: 10.1016/S0926-

re

9851(02)00239-2.

lP

Xia, J., Y. Xu, Y. Luo, R. D. Miller, C. Cakir, and C. Zeng, 2012, Advantages of using multichannel analys is of Love waves (MALW) to estimate nearsurface shear-wave velocity: Surveys in Geophysics, 33, 841–860, doi: 10.1007/s10712-012-9174-2.

Geosciences.

na

Xia, J., 2015, High-frequency Surface Wave Method (in Chinese): Press of China University of

ur

Yilmaz, Ö., M. Eser, and M. Berilgen, 2006, A case study of seismic zonation in municipal areas : The Leading Edge, 25, no. 3, 319–330, doi: 10.1190/1.2184100.

Jo

Zhang, P., L. Han, Z. Xu, F. Zhang and Y. Wei, 2017, Sparse blind deconvolution based low -frequency seismic data reconstruction for multiscale full waveform inversion: Journal of Applied Geophysics, 139, 91108, doi: 10.1016/j.jappgeo.2017.02.021.

Journal Pre-proof

Jo

ur

na

lP

re

-p

ro

of

Conflict of interest statement We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.

Journal Pre-proof Highlights We develop a novel preconditioned technique for elastic SH- and Love-wave full waveform inversion based on the inverse scattering theory in prestack depth migration. The role of preconditioned operator played in full waveform inversion is stated in detail. We also analyze whether the inconsistency of wavelets among all common-shot gathers would reduce the

Jo

ur

na

lP

re

-p

ro

of

intensity of seismic illumination.