Multiparameter Elastic Full Waveform Inversion of Ocean Bottom Seismic Four-Component Data Based on A Modified Acoustic-Elastic Coupled Equation

Minao Sun 1,2 and Shuanggen Jin 2,3* 1 Key Laboratory of Meteorological Disaster, Ministry of Education (KLME), Joint International Research Laborotory of Climate and Environment Change (ILCEC), Colloborative Innovation Center on Forecast and Evaluation of Meteorological Disasters (CIC-FEMD), Nanjing University of Information Science and Technology, Nanjing 210044, China; minaosun@nuist.edu.cn 2 School of Remote Sensing and Geomatics Engineering, Nanjing University of Information Science and Technology, Nanjing 210044, China 3 Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China * Correspondence: sgjin@shao.ac.cn


Introduction
Ocean bottom seismic survey is a modern platform for exploring the Earth's interior, locating seismometers at the seabed for all-weather, long-term, continuous, real-time observations. Unlike the conventional towed-streamer acquisition, OBS can record 1C pressure and 3C displacement data [1] using four-component (4C) detectors. The observed multicomponent data contain plenty of elastic properties of subsurface media, which can be used to deduct the lithology, fluid content, and pore pressure of rocks [2].
In multicomponent data processing, elastic full-waveform inversion (EFWI) plays an increasingly important role [3]. In the manner of classical FWI [4], EFWI computes parameter gradients by cross-correlating forward-and back-propagated wavefields and updates models to minimize the data misfit function. As governed by the elastic wave equation, EFWI can interpret multiple elastic wave phenomena, i.e., wave-mode conversion and AVO effects [5], and provide quantitative estimations for subsurface parameter distributions. Although it costs a large number of computing resources in wavefield simulations, its excellent performance still makes it more and more attractive [6][7][8][9][10].
However, the standard elastic wave equation commonly used in the conventional EFWI approaches cannot directly extract pressure components from elastic wavefields. By solving the acoustic and elastic wave equations in different computing areas, a fluid-solid coupled EFWI approach has been proposed [11][12][13]. It can generate pressure in the water immediately above the seabed and elastic components on the solid seabed. However, it requires to explicitly implement the correct boundary conditions, which is challenging for irregular surfaces [14]. Alternatively, Yu et al. [15] proposed an acoustic-elastic coupled (AEC) equation in elastic imaging of OBS 4C data. It introduces the physical relation between pressure and normal stress into the elastic wave equation. Thus, it can compute the pressure wavefield and avoid applying the boundary conditions. The developed 4C elastic reverse-time migration (ERTM) shows to suppress non-physical artifacts in the back-propagated wavefield and provide better-resolved subsurface images. With consideration of propagating direction and anisotropic property, this equation has been extended with an elastic vector imaging for transverse isotropy media [16,17]. However, these 4C ERTM methods aim to retrieve the subsurface structures but fail to provide quantitative parameter reconstructions.
In this study, we propose a new EFWI method based on a modified AEC equation, which can reconstruct multiple elastic parameters from OBS 4C data. Our method defines a new weighted misfit function; thus, it can adjust the weight of pressure and displacement components, and eliminate the differences in the order of magnitude. As more parameter classes and data components are involved, the blurring effects and parameter couplings [18] in the Hessian operator are prone to be more serious. To better consider the inverse Hessian operator, we reformulate the truncated Gauss-Newton-based (TGN) algorithm [19,20] in the framework of this modified AEC equation. Compared with the preconditioned conjugate gradient (PCG) algorithm, TGN can estimate a more accurate inverse Hessian and provide better parameter update directions. TGN has been widely used in multiparameter inversion for acoustic, elastic, and anisotropic media [21][22][23] and elastic least-squares RTM [24,25]. Besides, a pseudo-diagonal Hessian [26,27] is used as a precondition operator to remove the influences of limited observation apertures, geometry spreading, and frequency-limited wavelet.
The paper is organized as follows. In Section 2, we first review the general formulas of FWI, and then introduce the theory of the AEC-EFWI method and the implementation of the preconditioned TGN algorithm. In Section 3, we numerically analyze multiparameter sensitivity kernels of pressure and displacement components. In Section 4, we use two numerical examples to validate the effectiveness of the proposed method. Before conclusions, we discuss whether the AEC-EFWI method can invert elastic parameters using only 1C pressure data for OBS and towed-streamer acquisitions.

General FWI Formulation
Seismic wave equation can be expressed as where w denotes the subsurface wavefields, f indicates the source wavelet, and S is the parameter derivative matrix. By taking the partial derivative of Equation (1) with respect to parameter, the sensitivity kernel L can be acquired, The gradient of misfit function g can be obtained by applying the adjoint of sensitivity kernel to the data residuals between observed and synthetic data, Based on the Newton optimization, the parameter perturbation can be estimated by solving the Newton equation where H denotes the Hessian operator, which is the second-order partial derivatives of the misfit function with respect to parameter. Thus, the (k + 1)th iteration model m (k+1) can be updated by summing the (k)th iteration model m (k) and the (k + 1)th iteration model perturbation δm (k) scaled with a suitable step-length r,

Acoustic-Elastic Coupled EFWI Method
Compared with the standard elastic wave equation, the original AEC equation requires one more formula to compute the pressure component from elastic wavefields. In this study, we have made some modifications to reduce the number of equations and variables, and thus provide a modified AEC equation (details referred to Appendix A), given by where p, u x , and u z denote the pressure, horizontal , and vertical particle displacement wavefields, respectively. τ s n and τ s s are the S-wave-related normal and deviatoric stress components, respectively. λ and µ are the Lamé constants, and ρ is density. f p indicates the source function applied to the p-component. Compared with the original AEC equation, this modified one can generate OBS 4C records with same accuracy in wavefield simulation but costs less computing resources.
In this study, we define a weighted misfit function to account for both pressure and displacement components, given by where δd x , δd z , and δd p denote the residuals of the horizontal, vertical displacement, and pressure components, respectively. Here, ε is a weighting coefficient, satisfying ε ∈ [0, 1]. A scale factor ζ is used to eliminate the differences in magnitude between pressure and displacement components.
According to the adjoint-state theory [28], the adjoint AEC equation can be given by (see Appendix B for details) where (û x ,û z ,p,τ s n ,τ s s ) is the adjoint wavefields of (u x , u z , p, τ s n , τ s s ), and f x , f z , f p is the multicomponent adjoint source, satisfying The gradients of the Lamé constants g λ , g µ and density g ρ,Lame can be computed by performing zero-lag cross-correlations of the adjoint wavefields (Equation (8)) and the forward wavefields (Equation (6)), given by (see Appendix B for details) Compared with the Lamé constants, the parameterization of seismic velocities is a better choice in multiparameter EFWI [29][30][31][32]. According to the chain rule, the gradients of P-(α) and S-wave velocities (β) and density can be obtained, The models of P-and S-wave velocities and density can be updated as follows,

Preconditioned Truncated Gauss-Newton Algorithm
The inverse Hessian operator is estimated by a preconditioned truncated Gauss-Newton (PTGN) algorithm, as shown in Algorithm 1. In each iteration, we should perform demigration (Lp k ) and migration (L T (Lp k )) processes to update the parameter perturbations. The migration has been illustrated in Equation (8), and the demigration of OBS 4C data can be computed through a first-order Born modeling operator, given by

Algorithm 1 Preconditioned Truncated Gauss-Newton algorithm
Input: Gradient g, Hessian precondition operator H p ; Set x (0) = 0 and r (0) = L T L x (0) − g; Solve H p · y (0) = r (0) for y (0) ; Set p (0) = −r (0) and k = 0; Output: Parameter perturbation x 1: while r (k) > do 2: 3: A diagonal pseudo-Hessian with source-side illumination is used as a precondition operator, given by According to the modified AEC equation, the diagonal blocks of the Hessian for Lamé constants and density satisfy and the ones for P-and S-wave velocities and density can be given by

Sensitivity Analysis
In regional and global seismic explorations, sensitivity kernels are always used to portray subsurface wavepaths [33,34]. For the elastic case, the kernels of displacement components with different parameter classes have been studied on the standard elastic equation [23,35]. In this study, we use the modified AEC equation to simulate elastic wavefields, allowing for the kernels of pressure and displacement components. The experiment is performed on the elastic Marmousi model (see Figure 1), including a shallow water layer below the sea surface. Only one shot is excited to generate a pure P-wave source at the place of (5 km, 0.03 km), and 601 4C receivers are evenly located at the seabed. The peak frequency of the source function is 8 Hz.
Observed pressure and displacement records are back-propagated from the receivers, respectively. The obtained sensitivity kernels of pressure (p) and displacement (u x and u z ) components are presented in Figure 2. In the pressure kernels, most energy is distributed in the shallow part of the model and attenuates rapidly with the increase of depth. In contrast, the displacement kernels contain the S-wave reflection paths, thus it can enhance the illumination for the deep model. The corresponding 2D wavenumber spectrum of these kernels are shown in Figure 3. It is clear that the pressure kernels have higher resolution than the displacement ones in both vertical and horizontal directions. It is because the pressure component is computed by the spatial partial derivatives of the displacement wavefield, and these operators physically increase the frequency (or wavenumber) content of the data. Consequently, the simultaneous utilization of pressure and displacement components can provide a better characterization of subsurface structures.  The kernels for different parameter classes also have remarkable differences. The V p kernels (Figure 2a,d) are mainly formed by the diving waves and reflections associated with P-wave, which show relatively isotropic distributions in the wavenumber spectrum. The Vs kernels contain more information of PS reflections (Figure 2b,e). Besides, the ρ kernels are of the highest wavenumber components and behave as migration images of subsurface interfaces. That is the reason why we can easily obtain the short-wavelength structures of density model but fail to reconstruct the background model from seismic data.

Results
We use two numerical experiments on (1) the Overthrust model and (2) the Marmousi model to validate the proposed method. An O(2, 8) time-space-domain finite-difference staggered-grid solution [36] of the modified AEC equation is used to generate both forward-and back-propagated 4C wavefields. A convolution perfectly matched layer absorbing boundary [37,38] is used around the calculation area without consideration of the sea surface. In these experiments, a Ricker wavelet is adopted to generate pure P-wave sources with a peak frequency of 8 Hz (the bandwidth in [2 Hz, 20 Hz]).

Overthrust Model Test
We first use the Overthrust model to demonstrate the effectiveness of the proposed PTGN algorithm for OBC 4C data. The true V p and Vs models are shown in Figure 4. The density model is set to be 1000 kg/m 3 in the water layer and 2000 kg/m 3 below the seabed. The model is sampled as 801 × 166 grids with intervals of 12.5 m in both horizontal and vertical directions. The initial parameters are generated from the true ones with a smoothing window of 250 m. The acquisition geometry includes 101 shots with an interval of 100 m below the sea surface and 801 OBS receivers at the seabed for each shot. The total recording time is 4.0 s, and the temporal sampling rate is 1.0 ms. Observed seismic data are shown in Figure 5, including horizontal and vertical displacements and pressure components. As a comparison, the inversion is also performed using a preconditioned conjugate gradient (PCG). The maximum number of the loop for parameter update is 21, and that of the inner loop in the PTGN algorithm is 10.    Figure 6 displays the multiparameter gradients of the 1th iteration. We observe that the PCG gradients have insufficient illuminations for the deep part of the model, while the PTGN gradients are much improved and behave as amplitude-preserving subsurface images. The final V p and Vs models are displayed in Figure 7. The inverted V p and Vs using PTGN give better descriptions of structural boundaries with a higher interface continuity and fewer vertical artifacts. The vertical profiles (Figure 8) and the root mean square (RMS) errors (Table 1) can further demonstrate that PTGN can provide more accurate multiparameter inversions. Multicomponent data residuals between observed and simulated are shown in Figure 9. The residuals of the PTGN method are much weaken than those of PCG, which demonstrates that the PTGN can better interpret the observed multicomponent data. The convergence curves ( Figure 10) shows that PTGN has a higher decreasing rate and eventually converges to a lower misfit value.

Marmousi Model Test
Next, we use the elastic Marmousi model to demonstrate the effectiveness of the AEC-EFWI method. Smoothed versions of true V p, Vs, and ρ (see Figure 1) with a window of 300 m are taken as the initial models (see Figure 11). The dimension of the model is 601 × 201, and the intervals are 15 m in the horizontal and vertical directions. We have 61 shots with an interval of 150 m and 601 OBS receivers for each shot. The total recording time is 4.8 s, and the temporal sampling rate is 1.2 ms. A Ricker wavelet with a peak frequency of 8 Hz is adopted to generate a pure P-wave source. Observed multicomponent seismic data are simulated using the modified AEC equation, as shown in Figure 12. As a comparison, an EFWI method for horizontal and vertical displacement components are performed. A maximum of 15 iterations is used for the PTGN loop for both AEC-EFWI and EFWI.    Figure 13 displays the inverted V p and Vs models using the two methods. In the V p results (Figure 13a,c), AEC-EFWI provides better-resolved structures, i.e., anticlines, faults, lithologic interfaces, and high-speed bodies. For Vs (Figure 13b,d), however, the results using the two methods are comparable. With incorporation of the RMS errors in Table 2, we can find that considering the pressure data may not take as great effects on Vs as V p. It may be because the pressure data are more sensitive to the V p perturbations. The vertical profiles extracted at the horizontal distances of 3.0, 4.5, and 6.0 km are displayed in Figure 14. The AEC-EFWI results (marked by red lines) can precisely illustrate the deep reflectors with narrower sidelobes, and they are very close to the true models (marked by black lines). In contrast, as displayed in the corresponding wavenumber spectrum (Figure 15), EFWI underestimates the perturbations, especially for the interfaces with sharp parameter contrasts (see green lines in Figure 14). The data residuals simulated by the inverted results ( Figure 16) and the convergence profiles of the misfit function ( Figure 17) prove that the AEC-EFWI can better match the observed data and have a higher convergence rate.     Figure 13. The scale is consistent with the shot gathers in Figure 12. Panels (a-c) denote the residuals of EFWI, and panels (d-f) indicate those of AEC-EFWI.

Discussion
We have validated the effectiveness of the proposed AEC-EFWI method in processing OBS 4C data and improving inversion accuracy. This study mainly focuses on reconstructing the high-wavenumber components of elastic parameters from good initial models. Of course, this AEC-EFWI method also suffers from the notorious cycle-skipping and other practical issues. To alleviate this problem, this work should be further considered to combine with the reflection waveform inversion (RWI) [39][40][41][42][43] or the migration velocity analysis (MVA) [44][45][46]. In these approaches, the most critical step is to compute the PP and PS reflection paths. It can be easily accomplished by the AEC equation instead of performing a complete P/S decomposition on forward/back-propagated wavefields.
In those experiments, the weighting coefficient ε is set to be 50%. In fact, this value is determined by the difference of observed pressure and displacement components. Supposing that it reduces to zero, we wonder whether the AEC-EFWI method can still provide reasonable inversions for elastic parameters? Figure 18a shows a simple cartoon of the wave propagation process for this case. The incident P-wave excited from the source location generates both PP and PS transmissions at the seabed, and these transmissions are reflected at the interface of Layer 1. Because the water layer is assumed to be precisely known in advance, the observed PPP and PSP waves can be treated as "pseudo-first-order" reflections excited by virtual mixed sources from the seabed. Note that, this PSP wave path carries more information of subsurface Vs distribution, which makes a great contribution to Vs update.  We test the method on the Marmousi model (Figure 1), starting from the same initial models ( Figure 11) with 1C pressure data. The inverted results are displayed in Figure 19. The extracted vertical profiles and the corresponding wavenumber spectrum are shown in Figures 20 and 21, respectively. Although the inversion accuracy and spatial resolution decrease to some extent, the 1C results can still provide acceptable multiparameter inversions and have a certain consistency with the true ones. It demonstrates that this AEC-EFWI method is feasible to recover elastic parameters using OBS 1C pressure data. Besides, we can find that the high-wavenumber components in the results using 1C pressure data are better reconstructed than the low-wavenumber components (see differences between the green and red lines in Figure 21), which highlights the contribution of pressure component on high-wavenumber reconstruction.  Similarly, this method may be also applied for marine towed-streamer pressure data. As displayed in Figure 18b, the PP-PP and PS-SP wave paths help to reveal subsurface V p and Vs distributions. Besides, other converted waves, i.e., PP-SP and PS-PP, can further enhance the illumination of Vs model. Of course, it should be further tested with field streamer data. Although some practical problems, i.e., data preprocessing, low-frequency loss, and initial model building, have not been fully considered, the cartoons and the preliminary results can at least inspire us to use 1C pressure data for elastic parameter inversion and eventually provide a flexible method and idea for increasingly complex marine data processing.

Conclusions
In this study, we have proposed an elastic full-waveform inversion method based on a modified acoustic-elastic coupled equation. This method uses the modified AEC equation to simultaneously compute subsurface pressure and displacement wavefields. It adopts a weighted misfit function to quantify the contributions of pressure and displacement records. With the adjoint operator, it can simultaneously make use of OBS 4C data to reconstruct multiple elastic parameters. The preconditioned TGN algorithm with a multiparameter diagonal Hessian operator is developed to cope with unbalanced illumination and coupling effects. The sensitivity analysis and numerical experiments have validated that the AEC-EFWI method can yield a higher spatial resolution, provide more accurate elastic parameter inversions, and have a higher convergence rate. The discussion reveals the potential of this AEC-EFWI method in inverting elastic parameters using 1C pressure data for OBS and marine towed-streamer cases.  Acknowledgments: We thank Pengfei Yu from Hohai University for providing the code for the original acoustic-elastic coupled equation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A. Derivation of the Modified Acoustic-Elastic Coupled Equation
We start from the standard 2D time-domain displacement-stress elastic wave equation [47], given by where u x and u z denote the horizontal-and vertical-particle displacement components, respectively; τ xx and τ zz are the normal stress components; and τ xz is the shear stress component. λ and µ are the Lamé constants, and ρ is density. f indicates the source function that is implemented in τ xx and τ zz to generate a pure P-wave.
In tensor analysis, the stress tensor T can be decomposed into the isotropic pressure −pI and deviatoric τ s parts, For 2D cases, the pressure wavefield satisfies and the deviatoric stress components become These equations formulate the original acoustic-elastic coupled equation; we refer to the work in [15].
Note that, it is redundant to simultaneously compute the pressure and two normal stress components. Thus, we redefine where τ s n and τ s s are the redefined normal and deviatoric stress components. The simplified equation becomes Compared with the original one, this modified equation can provide subsurface 4C elastic wavefields with same accuracy but less computing resources.
For completeness purpose, we also derive the modified AEC equation in 3D. The 3D original AEC equation is given by Because τ s xx + τ s yy + τ s zz = 0, we can similarly provide a 3D modified AEC equation by replacing τ s yy as −(τ s xx + τ s zz ) to reduce the number of formulas.

Appendix B. Derivation of Gradient Computation for Aec-Efwi Method
The modified AEC equation can be rewritten as where w = (u x , u z , p, τ s n , τ s s ) T denotes the elastic wavefields, f = 0, 0, f p , 0, 0 T denotes the source wavelet, and S denotes the parameter derivative matrix, Substituting Equation (2) into Equation (3), the gradient satisfies With consideration of the self-adjoint assumption S −1 T = S T −1 , it can be rewritten by Here, S T is the adjoint operator of the modified AEC-equation, given by According to the adjoint-state method, we defineŵ = S −1 T f as the solution of the adjoint equation, S Tŵ = f (A14) whereŵ = (û x ,û z ,p,τ s n ,τ s s ) T is the adjoint variables as used in Equation (8). For the parameterization of Lamé constants and density m = [λ, µ, ρ] T , the partial derivative matrices can be given by (A18)