2D multi-parameter waveform inversion of land reflection seismic data obtained from the particle-motion response from the vertical geophone

Onshore seismic exploration analyzes seismic wave propagation in elastic media, which includes the conversion between P- and S-waves. The development of multi-wave and multi-component seismic exploration methods provides data that enable onshore elastic wave full-waveform inversion. However, most data sets of onshore exploration are single component obtained from the particle-motion response from the vertical geophone. When the aiming area has a low-velocity zone, the ray path of reflected wave that propagates to the detector is nearly perpendicular to the ground surface, so that we call it P-wave data. In this paper, we focus on multi-parameter waveform inversion using P-wave reflection seismic data. Although only P-wave data are received, it still contains the converted P-wave information, and the converted P-wave energy gradually increases as the offset increases. As seismic acquisition technology, observation systems and science develop, the folds and acquisition offset increase significantly, and the seismic data contain important converted P-wave information. In this paper, the first-order elastic velocity–stress equation is decomposed to obtain the scalar-P-wave equation from which the S-wave velocity is included firstly. Then we present the theoretical framework for onshore multi-parameter full-waveform inversion using P-wave data. In order to explore the inversion potential of the P-wave data (extracting the S-wave velocity from the converted P-wave information) and accuracy and stability of the P- and S-wave velocity inverted by our method, we carry out numerical tests via different inversion strategies, by using the P-wave data regarded as containing converted P-wave information, and get successful results.


Introduction
Full-waveform inversion (FWI) is a quantitative data-fitting procedure which minimizes the residuals between the simulated and the observed seismic data to obtain high-resolution subsurface physical parameters such as P-wave velocity, density, impedance, or anisotropic parameters ). Lailly (1983) and Tarantola (1984) first proposed the least-squares misfit functional FWI framework and built the gradient of the misfit function by cross-correlating the incident wavefield emitted from the source and the back-propagated residual wavefields. Pratt and Worthington (1990) extended FWI to the frequency domain and pointed out that multi-scale inversion can improve the effectiveness and stability of the inversion. At present, FWI has been successfully applied for both simulated and real seismic data (Shipp and Singh 2002;Sears et al. 2008;Sirgue et al. 2010).
Since the 1990s, the emergence of the ocean bottom cable (OBC) system has caused a dramatic change in marine multiwave seismic exploration technology. Since 2000, onshore multi-wave multi-component exploration has gradually developed. For multi-component seismic data, elastic wave fullwaveform inversion (EFWI) is more suitable for high-precision reconstruction of subsurface elastic parameters than acoustic full-wave inversion (AFWI) (Brossier et al. 2009). In some 1 3 complex areas, AFWI cannot describe elastic wave propagation processes (Barnes and Charara 2009), such as the P/S converted wave and amplitude variation with offset (AVO) effects. EFWI can provide more accurate P-and S-wave velocities for seismic data processing and interpretation, and also provide high-quality parameter information for reservoir description and for oil and gas field development (Tatham and Stoffa 1976).
When multiple parameters are simultaneously processed, the nonlinearity of FWI increases because certain parameters are coupled to each other (Wu and Aki 1985). Many scholars (e.g., Tarantola 1986;Mora 1988;Crase et al. 1990) discussed the elastic wave waveform inversion theory and applied it to simulated and actual data. They developed and improve the forward modeling and inversion theory, iterative optimization method, misfit function, and inversion strategy. Pratt (Pratt et al. 1998(Pratt et al. , 1999a also proposed the frequency-domain EFWI theory while presenting the AFWI. Crase et al. (1990) used the first-order velocity-stress equation to invert the P-wave impedance, based on various minimization criteria.
Compared with the P-and S-wave velocity parameters, the density variation mainly affects the amplitude information of the seismic wavefield, and has little effect on wave diffusion parameters such as travel time and phase. Therefore, density is the most difficult parameter to be inverted in multi-parameter inversions. Tarantola (1984) proposed a velocity-density inversion strategy based on the acoustic wave equation with variable density. Mora (1987) pointed out that near-offset seismic data are suitable for impedance parameter inversion, while the far-offset data are suitable for velocity parameter inversion, but neither method can adequately invert the density model. Köhn et al. (2012) showed that the use of velocity-density parameterization in multi-parameter EFWI is better than Lame constant-density parameterization. Jeong et al. (2012) proposed a multi-parameter EFWI strategy in the frequency domain, which first inverts the elastic modulus, then simultaneously inverts the velocity and density.
Multi-wave and multi-component seismic exploration technology are still developing, and many problems still need to be solved. At present, most onshore exploration is based on single component, which is collected from obtained from the particle-motion response from the vertical geophone. In this situation, there is a classical assumption that the subsurface is regarded as a fluid (Raknes et al. 2015). When the area has a low-velocity zone, which causes the ray path of reflected wave that propagates to detector is nearly perpendicular by Snell's law working, we call the field data set the P-wave data. Acoustic wave equation is generally used to match the P-wave data, and AFWI is a well-known method in the application of current onshore full-waveform inversion, but it has no possiblity to invert S-wave velocity. However, the seismic data does not contain the S-wave information directly, while it includes a large amount of converted P-wave information, such as P-S-P-wave. It would therefore be possible to establish a scalar-P-wave multiparameter equation in the elastic wavefield to synthesize record to match the real seismograms.
In this paper, we decompose the elastic first-order velocity-stress equation and obtain the scalar-P-wave equation that includes the S-wave velocity. We apply the conventional acoustic equation, the scalar-P-wave equation and the elastic wave equation to a simple layered model, then analyze the converted wave characteristics of the three equations by comparing the different forward results, and evaluate the S-wave velocity inversion. Then we propose the onshore FWI algorithm based on P-wave data and carry out different inversion strategies using the Marmousi-2 model. Our results confirm the accuracy and stability of the inversion of the P-and S-wave velocities using P-wave data with converted P-wave information.

Why the acoustic approximation can work on land seismic
In exploration geophysics, FWI has come into broad use for applications (Sirgue and Pratt 2004). Generally, scholars make an assumption that subsurface can be considered as a fluid, as we called acoustic approximation (Tarantola 1984). Not considering the limitations in available computer resources, there are two main reasons for this assumption. One is that seismic usually acquire measurements of P-wave, whether current marine data obtained from the pressure response of hydrophone (Vigh et al. 2014) or land data obtained from the particle-motion response from the vertical geophone (details in the next paragraph) correspond to compressional wavefield (Gaiser et al. 2001). The other is that conventional processing focuses on the kinematics of P-waves, and it is well modeled by acoustic wave equation (Raknes et al. 2015).
Different from marine data that obtain pressure response directly, land seismic data usually acquire measurements of the vertical particle velocity (Simmons and Backus 2003). As we all know, the numerical solution of elastic wave equation is more reliable than that of acoustic wave equation to approximate field data ). Actually, land seismic data contains converted S-wave, which could have a very strong impact on the acoustic inversion results. Fortunately, the near surface is a low-velocity zone, which causes the ray path of reflected wave that propagates to detector is nearly perpendicular to the ground surface (Snell's law works), as shown in Fig. 1. Thus, land seismic data obtained from the vertical geophone can be considered as P-wave, and it can be assumed that the P-and S-wave are naturally decoupled in land seismic data as shown in Fig. 1. However, this assumption is not applicable in all cases, such as strong anisotropy of near surface and exposed bedrock environment, which are not considered in this paper.
In acoustic FWI, the subsurface medium is parameterized by density (ρ) and P-wave velocity ( v P ), and well inverted ). In the field, the land seismic data are obtained only from vertical particle velocity from the vertical geophone, while in the real application, the measurements of seismograms is pressure but not velocity, because the most choice of acoustic equation is second-order system (Ravaut et al. 2004;Operto et al. 2006), in which the vertical particle velocity is implicit. In fact, almost all research use simulated pressure data to match vertical velocity data observed, directly. Then we analyze the relationship between the vertical velocity field and the pressure field and explain that why the seismic data can be modeled by acoustic wave equation in the field.
For the convenience of research, we give the first-order system acoustic wave equation as, where P is pressure wavefield, K is bulk modulus, is density, f is the source, = v x , v z is particle velocity vector, and t is the time.
The pressure and velocity record modeled by acoustic wave equation are shown in Fig. 2a-c, and the travel time and phase information of reflected waveform in vertical velocity record are the same as that in pressure record. For the medium with strong contrasts, the pressure date and velocity data have significant amplitude mismatch at far offsets, which could be solved by data regularization (Hobro et al. 2014). This explains that why the pressure field is used directly to match observed data in conventional land seismic FWI. In this paper, we assume the medium is full elastic, and propose a new multi-parameter waveform inversion method to invert the S-wave velocity ( v S ). In our method, the P-wave data are considered to contain rich convert P-wave information to match observed vertical velocity data, which details in the next section.
In addition, the land FWI is also affected by near-surface heterogeneity, anisotropic, attenuation and ground roll, which have a very strong impact on the inversion results and are ignored here.

Scalar-P-wave equation based on elastic decomposition
The application of FWI in industry typically invert twoparameter models, as v P and ρ. In this work, in additional to aforementioned models, we also try to invert S-wave velocity ( v S ). Because of the natural decomposition of P-and S-wave in land seismic data, converted S-wave is not collected directly in vertical geophone, and the S-wave velocity can only be indicated by converted P-wave.
In order to match Methods of decomposition for elastic field are generally divided into three categories. The first is Helmholtz decomposition (Dellinger and Etgen 1990;Sun et al. 2007). The second is vector-based wave equation decomposition (Ma and Zhu 2003;Xiao and Leaney 2010). The third is wavenumber-domain field vector decomposition (Zhang and McMechan 2010). This paper focuses on the elastic  information of P-wavefield based on acoustic approximation and its potential ability to invert v S . Thus by decomposing constitutive equations into the swell-shrinking (pressure) and shear-rotating forms to decompose the elastic wavefield, we obtain the scalar-P-wave equation (details in "Appendix 1"), where Č ijkl is new stiffness tensor, and ij is Kronecker's symbol. All indices change from 1 to 3, and rule in Einstein's summation convention over repeated indices.
In isotropic medium, it can be written as: If we set = 0 , the is equal to K, and scalar-P-wave equation can be rewritten as, Equation (5) is conventional acoustic wave equation, in which the pure pressure P is different from scalar-P-wave in Eq. (2), which can be regarded as a pseudo-pressure in elastic medium. Like acoustic FWI, we must first study the relationship of the vertical velocity field and the pseudopressure field, before the seismic data are modeled by scalar-P-wave equation.
We construct a three-layer elastic model with model size of 2400 × 1500 m, the elastic parameters of first layer are: P-wave velocity of 2500 m/s and S-wave velocity of 1785 m/s; the second layer: P-wave velocity of 3000 m/s and S-wave velocity of 2140 m/s; the third layer: P-wave velocity of 3500 m/s and S-wave velocity of 2500 m/s, while all the density parameter is 1420 kg/m 3 . The source, set at (500 m, 10 m), is Ricker wavelet with peak frequency of 20 Hz. Figure 3a, b shows vertical velocity ( v z ) and vertical P-wave ( v Pz ) velocity elastic wavefield snapshots. Compared with acoustic wavefield (Fig. 3c), elastic v z field has a lot of converted P-and S-wave, which are decoupled as v Pz and v Sz in land seismic data (see in Sect. 2.1). As shown in Fig. 3b, for the vertical P-wave component, the converted waves generated at the reflection points of two interfaces are visible (Tang and McMechan 2017). Moreover, a converted P-S-Pwave is generated at the second interface, which exists in pseudo-pressure field modeled by scalar-P-wave equation (Fig. 3d) in the same elastic case, but not exist in pure pressure field modeled by acoustic wave equation in acoustic medium (setting v S = 0 ). Thus v Pz of elastic equation or pseudo-pressure of scalar-P-wave equation can be directly used to match land seismic data in elastic FWI, because of decomposition of P-and S-wave, whose ray path of seismic wave perpendicular to the surface. However, pseudopressure modeled by scalar-P-wave equation is a better choice to match observed data. There are two reasons: (1) in order to get v Pz , the computation of elastic modeling by P/S-wave decomposed method(Ren and Liu 2016) requires higher computational cost than scalar-P-wave modeling; (2) conventional processing most focuses on the kinematics of P-waves, and it is well modeled by scalar-P-wave data. Comparison of v Pz record and pseudo-pressure record as shown in Fig. 4b, d illustrates that converted waves generated of two fields in elastic case are consistent, and this proves the possibility of inversion of S-wave velocity under pseudo-acoustic approximation, which cannot works on acoustic record (Fig. 4c) under pure acoustic approximation. Tarantola (1984) defines a local optimization problem of least-squares minimization of the misfit function E between an observed record P obs and simulated data P cal , where m i denotes the elastic parameters: v P ,v S , and ρ, x r and x s denote receivers and sources.

2D scalar-P-wave equation FWI
The gradient expression of 2D scalar-P-wave equation FWI can be derived by applying the first-order Born approximation (details in "Appendix 2"), Since the parameterization method using the Lamé constant and density is not a reliable method (Tarantola 1986;Köhn et al. 2012), we chose the velocity-density model parameterization. The gradient of the misfit function with respect to the other parameters can be obtained using the chain rule (Mora 1987): Thus, the velocity-density gradient of the elastic multiparameter inversion can be expressed as, There are various inversion algorithms used in FWI, such as Gauss-Newton, conjugate-gradient, and Newton method. This paper we apply the conjugate-gradient method (Fletcher and Reeves 1964) to update elastic parameters model with the following steps, where n denotes the current iteration round, t denotes vector transpose operator, represents the iteration step length, calculated by a one-dimensional line search, and ∇ m i E represents the gradient of the misfit function with respect to the elastic parameters, k n is update direction.

Marmousi model example
In this section, we test the scalar-P-wave FWI with the wellknown 2D Marmousi model. The model shown in Fig. 6 has a grid size of 551 × 251, a sampling interval of 10 m, and a time interval of 0.8 ms. The source excited 8 m under the (10) g n = ∇ m i E = g t n g n g t n−1 g n−1 ground. The source wavelet is a Ricker wavelet with peak frequency of 7 Hz. Fifty-six shots totally are exploded with the 10 m distance between adjacent sources. The geophones were placed on the surface. We generate elastic pseudo-pressure data by Eq. (2), and Fig. 5 shows the real and initial velocity models. We now carry out two tests: (1) one step: invert P-and S-wave velocities simultaneously by scalar-P-wave FWI; (2) two steps: firstly invert P-wave velocity singly, then invert S-wave velocity. The density is regarded as a constant of 1500 kg/m 3 in all tests.
We apply the one-step inversion method to the elastic multi-parameter waveform inversion: the P-and S-wave velocity are inverted simultaneously. The inversion results are shown in Fig. 6. As we can see, the P-wave velocity inverted model is more reliable, but the S-wave velocity looks worse because of the crosstalk effects between various parameters.
To overcome the drawbacks of the one-step inversion method, we applied the two-step inversion method to the elastic multi-parameter waveform inversion: first, the S-wave velocity (Fig. 5d) is not updated and the P-wave velocity is inverted; then, the S-wave velocity is inverted by using the inverted result from the first step as the initial P-wave velocity model. The results are shown in Fig. 7, and the S-wave velocity inverted model is much better than that inverted simultaneously. Figures 8 and 9 show the depth profiles of P-and S-wave from the true, initial and inverted (singly and simultaneously) models at x-axis locations of 1, 2, 3, 4, and 5 km. The trace contrast curve (Fig. 8) extracted from the P-wave velocity models shows that the resolution of two inverted models computed singly and simultaneously are similar. The depth profiles of S-wave velocity model derived from the singly inversion results (Fig. 9) are greatly improved compared with the simultaneous inversion method, and the S-wave velocity is well recovered. Although the singly inversion produces good results, the cross talk effects between the various parameters still exist and cause some deviation from the true value in the deep layer, which could be mitigated effectively by resorting to Hessian-based method. In conclusion, the two-step inversion strategy is a better option for multi-parameter waveform inversion.

Field data example
In this section, we apply a real 2D seismic data set obtained from the particle-motion response from the vertical geophone to verify the scalar-P-wave equation FWI method. This area has a low-velocity zone in near surface, so that the ray path of reflected wave is nearly perpendicular to the ground surface, thus the data set can be regarded as P-wave data, which can be implemented with our scalar-P-wave FWI method.
We use a Fourier transform method to make data interpolation. The sources and receivers are set at the surface. There are 40 shots from x = 1 to 5 km and the shot space is 100 m, and the largest offset is 4 km. The receiver spacing is 10 m. The time spacing is 0.8 ms, and maximum time is 2.88 s. The source wavelet is obtained from seismograms analysis. Figure 10 shows the initial P-wave velocity models computed by migration velocity analysis (MVA), while S-wave velocity model and density are not given. Thus we must find a perfect match empirical relationship between elastic parameters. In     (Gardner et al. 1974), and S-wave velocity model is given by relationship: v s = v p ∕1.41 . Figure 11 displays the P-and S-wave velocity inversion results simultaneously after 23 iterations. As we can see, the reflections in P-wave velocity model is more continuous. However, there are some discontinuous events and noise in low depth area in S-wave velocity models, and there are two reasons: (1) the trade-off between P-and S-wave; (2) the data set mainly contains P-wave information. We synthesize one shot record by finite difference forward modeling method with inverted P-and S-velocity model, and we compare them with the real seismogram in the same picture, which is illustrated in Fig. 12. Except for surface wave and direction cut in inversion, the main reflections match well.

Conclusions
In this paper, we propose a multi-parameter full-waveform inversion method using the P-wave data with high-quality converted P-wave information. The following conclusions were drawn from model testing.
1. The scalar-P-wave equation obtained by decoupling the first-order velocity-stress equation contains only P-wave information, unlike the conventional acoustic and elastic wave equations. Compared with the acoustic equation, scalar-P-wave equation can simulate the converted P-waves at the elastic interface, hence it can indicate the variability of the S-wave velocity. This provides information similar to that of the elastic wave equation. The P-wave information received from actual seismic data obtained from the particle-motion response from the vertical geophone is also similar, when the aiming area has a low-velocity zone. 2. The scalar-P-wave equation proposed in this paper has the ability to simulate the converted P-waves, and can invert the P-and S-wave velocity and density parameters.  Because of the strong crosstalk effects between the density and the other elastic parameters, density can be generally regarded as a constant. We compared the different inversion methods and their effect on the elastic multi-parameter inversion results, and the singly inversion method performed better in the Marmousi-2 model inversion test. 3. The scalar-P-wave multi-parameter inversion method is the only method that evolves from P-wave exploration to multi-wave and multi-component exploration. As P-wave data still dominates onshore seismic surveys, the proposed method has great significance for inverting S-wave velocity parameters based on P-wave.