Seismic Full-Waveform Inversion Using Decomposed P-Wavefield

JOURNAL OF SEISMIC EXPLORATION 29, 201-224 (2020) 201 SEISMIC FULL-WAVEFORM INVERSION USING DECOMPOSED P-WAVEFIELD SOOYOON KIM1, WOOKEEN CHUNG2, SUNGRYUL SHIN2 and DAWOON LEE1 1 Department of Ocean Energy and Resources Engineering, Korea Maritime and Ocean University, Busan, South Korea. [email protected] 2 Department of Energy and Resources Engineering, Korea Maritime and Ocean University, Busan, South Korea. (Received May 21, 2019; revised version accepted February 25, 2020) ABSTRACT Kim, S.Y., Chung, W.K., Shin, S.R. and Lee, D.W., 2020. Seismic full-waveform inversion using decomposed P-wavefield. Journal of Seismic Exploration, 29: 201-224. Here we describe the development of a seismic full-waveform inversion method which employs P-wavefield decomposition to obtain accurate velocity information. Briefly, P-wavefield decomposition for multi-component data was performed with Helmholtz decomposition in elastic media and an objective function. To achieve efficient inversion, application of a back-propagation technique is essential. Therefore, a stress tensor was used for P-wavefield decomposition to allow application of a backpropagation technique. Our proposed inversion algorithm was validated with synthetic data obtained from the Marmousi2 velocity model which simulated an ocean bottom, multi-component survey. The subsurface information obtained with our inversion method was more accurate in regard to velocity and structure compared with a conventional elastic inversion method. In addition, the application of our inversion method to synthetic data simulating an ocean bottom seismometer survey which uses a small number of receivers also obtained better results in a numerical test. KEY WORDS: full-waveform inversion, P-wavefield decomposition, multi-component, elastic. INTRODUCTION Seismic full-waveform inversion is the process of obtaining a subsurface property (i.e., velocity) of an exploration area. This process involves an objective function which calculates the difference between field 0963-0651//20/$5.00 © 2020 Geophysical Press Ltd. 202 data and modeling data. This difference is minimized when subsurface properties are updated. The final model is then assumed to represent an exploration area of interest. It has been proposed that calculation of gradient direction can improve the objective function of an inversion process, with zero-lag cross-correlation of partial derivative wavefield and residual of data performed. However, this approach requires a large number of calculations. Lailly (1983) and Tarantola (1984) both proposed efficient calculations for gradient direction with use of a back-propagation technique which uses wave equation-based migration. These approaches have reduced calculation time and the resources needed for waveform inversion work, thereby facilitating inversion studies. Initially, a subsurface was assumed to represent an acoustic media. Accordingly, Kolb (1986) and Gauthier et al. (1986) performed waveform inversion for acoustic media in the time domain. However, representing the subsurface as an acoustic media does not account for possible effects of complex waves such as S-waves, surface waves, and converted waves. Therefore, Mora (1987) expanded waveform inversion to elastic media, and this was subsequently studied in various domains. For example, Crase et al. (1990), Shipp and Singh (2002), and Sheen et al. (2006) studied waveform inversion for elastic media in the time domain, while Choi et al. (2008), Lee et al. (2010), and Jeong et al. (2012) studied it in the frequency domain. Data processing of an elastic media can provide various types of information, including P-wave and S-wave data. However, calculations of gradient direction for updating models are based on determinations of PP, PS, SP, and SS reflectivity at each grid point in a data processing approach which uses cross-correlation with parameters such as migration. Furthermore, when performing full-waveform inversion in elastic media, crosstalk from the scattering of each parameter can lead to contamination (Operto et al., 2013; Operto et al., 2018; Pan et al., 2019). To reduce this crosstalk, wavefield decomposition has been studied by Dankbaar (1987), Dellinger and Etgen (1990), Amundsen et al. (1998), and Sun (1999). In addition to the application of wavefield decomposition to the processing of seismic data, wavefield decomposition has also been applied to migration by Yan and Sava (2008, 2009), Yan (2010), Zhang and McMechan (2011), Chung et al. (2012), and Ha et al. (2015). Among these studies, Chung et al. (2012) defined a new migration technique involving P-wavefield decomposition, a procedure which was originally proposed by Dellinger and Etgen (1990). This techique obtained clear subsurface information in numerical tests. In the present study, P-wavefield decomposition was applied to waveform inversion to obtain accurate subsurface property information. For this, we set a new objective function for P-wavefield decomposition. 203 However, when Dellinger and Etgen (1990) defined waveform inversion theory with P-wavefield decomposition, it was demonstrated that derivation with a back-propagation algorithm is very difficult due to divergence of the operator in the objective function. Ha et al. (2015) proposed elastic reverse time migration using a stress tensor. When we applied this decomposition method to inversion and a back-propagation algorithm for efficient inversion, we were able to develop a full-waveform inversion theory using a decomposed P-wavefield procedure. To demonstrate that the proposed inversion method can be appropriately applied to multi-component data, inversion results generated with our proposed method were compared with those obtained with a conventional inversion method in elastic media in numerical tests. (Hereafter, conventional inversion refers to multi-parameter elastic inversion.) Numerical tests were performed with synthetic data and the Marmousi2 model, and modeling and inversion were subsequently performed to simulate ocean bottom cable (OBC) and ocean bottom seismometer (OBS) multi-component survey data. Inversion results were also compared with true velocity models and depth profiles. THEORY The objective function in seismic waveform inversion is composed of a l2-norm between field data and modeling data for a single source and can be defined as follows (Claerbout and Muir, 1973): !"#$ !!"# ( 𝒖! 𝑡 − 𝒅! 𝑡 )! , 𝐸= (1) !!! !!! where 𝐸 is an objective function, 𝑁𝑟𝑒𝑐 is the number of receivers, 𝒖 indicates modeling data, 𝒅 indicates field data, and 𝑇!"# is the total record time. Gradient direction minimizes the objective function and is the first derivative form for the model parameter, 𝑝! , as shown in eq. (2). However, gradient direction must be calculated. Alternatively, it is possible to perform an indirect calculation [see eqs. (2)-(5)] by using a backpropagation technique. !" !!! = !"#$ !!! !!"# !𝒖! (!) !!! { !! ! ! ⊗ (𝒖! (𝑡) − 𝒅! (𝑡))} , (2) 204 𝜕𝐸 = 𝜕𝑝! 𝜕𝐸 = 𝜕𝑝! 𝜕𝐸 = 𝜕𝑝! where !"#$ !!"# ! { 𝑽⋆ ∗ 𝑮! ⊗ 𝒓! (𝑡)}, 𝒓! (t) = [𝒖! (𝑡) − 𝒅! (𝑡)] (3) !!! !!! !"#$ !!"# ! { 𝑽⋆ ∗ 𝑮! ∗ 𝒓! (𝑇!"# − 𝑡)}, (4) !!! !!! !!"# ! { 𝑽⋆ ∗ (𝑩(𝑡))}, (5) !!! !" !!! is the gradient direction with 𝑝! as a model parameter obtained from eq. (1), !𝒖(!) !!! is a partial derivative wavefield at the k-th point, ⊗ is a cross correlation operator, 𝑮 is Green’s function, 𝑽⋆ is the virtual source, 𝒓 is residual, ∗ is a convolution operator, and 𝑩 is a back-propagated wavefield. The latter, represented as 𝑩(𝑡) in eq. (5), is the same as the 𝑮 ∗ (𝒓(𝑇!"# − 𝑡)) term in eq. (4). Furthermore, this back-propagated wavefield is obtained from modeling with residual between modeling data and field data as sources. By using a back-propagated wavefield and a virtual source as shown in eq. (5), the gradient direction is efficiently calculated and this method is referred to as a back-propagation technique. P-wavefield decomposition In this study, two P-wavefield decomposition techniques were applied to perform waveform inversion: the wavefield decomposition techniques proposed by Dellinger and Etgen (1990) and by Ha et al. (2015). The P-wavefield decomposition proposed by Dellinger and Etgen (1990) based on Helmholtz decomposition is described as follows: ∇∙𝒖= 𝜕𝑢! 𝜕𝑢! + , 𝜕𝑥 𝜕𝑧 (6) 205 where ∇ ∙ 𝒖 is the decomposed P-wavefield of modeling data and 𝑢! and 𝑢! represent displacement along the x- and z-axes, respectively. Pwavefield decomposition with the stress tensor proposed by Ha et al. (2015) is described in eq. (7). This method is based on use of a time domain, twodimensional (2D) staggered grid finite difference method for elastic media as described in eqs. (8)-(10) ∇∙𝒖= !!! !" + !!! !" 𝜏!! = 𝜆 + 2𝜇 𝜏!! = 𝜆 + 2𝜇 = !!! !!!! ! !!! (7) , 𝜕𝑢! 𝜕𝑢! +𝜆 , 𝜕𝑥 𝜕𝑧 !!! !" +𝜆 𝜏!! + 𝜏!! = 2 𝜆 + 𝜇 !!! !" !!! !" (8) (9) , + !!! !" , (10) where 𝜏!! and 𝜏!! are stress tensors and 𝜆 and 𝜇 are Lamé constants. To check the accuracy of P-wavefield decomposition in eq. (7) compared to that in eq. (6), we performed a simple numerical test in elastic homogeneous media. The P-wave velocity of the model was set to 2.0 km/s, S-wave velocity was set to 1.0 km/s, and density was set to 1000.0 g/cm3. The positions of the source and receivers are shown in Fig. 1. Snapshots and seismograms at 1.2 s are presented in Figs. 2 and 3. Fig. 1. Elastic homogeneous model for P-wavefield decomposition tests. 206 In Fig. 2, snapshots are provided of horizontal (a) and vertical (b) displacement, as well as a decomposed P-wavefield with divergence operator (c). In panels (d), (e), and (f), the seismogram is at the receiver position. Fig. 2. P-wavefield decomposition with a divergence operator in homogeneous media. The following snapshots are presented with the seismogram at the receiver position: (a, d) horizontal displacement, (b, e) vertical displacement, and (c, f ) decomposed Pwavefield with divergence operator. 207 The stress tensor in elastic media has a different component at horizontal and vertical. In the case of using the stress tensor in elastic media like an eq. (7), we can confirm the decomposed P-wavefield. In Fig. 3, snapshots are provided of the horizontal (a) and vertical (b) stress tensor, as well as a decomposed P-wavefield with a stress tensor (c). Panels (d), (e), and (f ) have the seismogram at the receiver position. Fig. 3. P-wavefield decomposition with a stress tensor in homogeneous media. The following snapshots are provided with the seismogram at the receiver position: (a, d) horizontal displacement, (b, e) vertical displacement, and (c, f ) decomposed Pwavefield with a divergence operator. The decomposed P-wavefields shown in Figs. 2 and 3 were found to correspond to each other in numerical tests. Thus, P-wavefield decomposition with a stress tensor as described in eq. (7) can replace decomposition with a divergence operator as described in eq. (6). 208 Full waveform inversion using P-wavefield decomposition To apply P-wavefield decomposition to waveform inversion, we proposed a new objective function in this study. An objective function for a single source is described in eq. (11), and the gradient direction of the k-th model parameter is described in eq. (12): !"#$ !!"# ( 𝛻 ∙ 𝒖! 𝑡 − 𝛻 ∙ 𝒅! 𝑡 )! , 𝐸= (11) !!! !!! 𝜕𝐸 = 𝜕𝑝! !"#$ !!"# !!! !!! In eq. (12), 𝜕𝛻 • 𝒖! (𝑡)! { ⊗ (𝛻 • 𝒖! (𝑡) − 𝛻 • 𝒅! (𝑡))}. 𝜕𝑝! !"•𝒖(!) !!! (12) is a partial derivative wavefield of decomposed P- wavefield. A back-propagation technique should be used for efficient calculation of gradient direction. However, this is difficult due to the divergence operator present in the partial derivative wavefield shown in eq. (12). Therefore, we applied P-wavefield decomposition with a stress tensor as proposed by Ha et al. (2015) to inversion and a back-propagation algorithm. The developed gradient direction of model parameter, which is applied above decomposition to the partial derivative wavefield in eq. (12), is shown in eq. (13), which is subsequently expanded in eqs. (14) and (15). 𝜕𝐸 = 𝜕𝑝! 𝜕𝐸 = 𝜕𝑝! !"#$ !!"# !!! !!! 𝜕 𝜏!!! + 𝜏!!! ! { ( ) ⊗ (𝛻 ∙ 𝒖! (𝑡) − 𝛻 ∙ 𝒅! (𝑡))}, 𝜕𝑝! 2 𝜆! + 𝜇! !"#$ !!"# ( !!! !!! 1 2 𝜆! + 𝜇! ) 𝜕𝜏!!! ! 𝜕𝑝! ⊗ 𝒓! 𝑡 + 𝜕𝜏!!! ! 𝜕𝑝! (13) ⊗ 𝒓! 𝑡 , (14) 𝒓! (t) = [𝛻 ∙ 𝒖! (𝑡) − 𝛻 ∙ 𝒅! (𝑡)] 209 𝜕𝐸 = 𝜕𝑝! !!"# ! ! { 𝑽!!! ⋆ ∗ 𝑮! ⊗ 𝒓(𝑡) + 𝑽!!! ⋆ ∗ 𝑮! ⊗ 𝒓(𝑡)}. !!! (15) 𝒓 t = 𝒓𝒋 t 2 𝜆! + 𝜇! , [𝑗 = 1,2, ⋯ , 𝑁𝑟𝑒𝑐] If the model parameter is 𝜆, the virtual sources 𝑽!!! ⋆ and 𝑽!!! ⋆ in eq. (15) are composed as follows: ⋆ ⋆ 𝑽!!! = 𝑽!!! = !!! !" + !!! !" (16) . Elastic parameters such as 𝜆 , 𝜇 , and 𝜌 can be updated with a calculation of gradient direction according to this process, with the only difference being the definition of the virtual sources. It should also be noted that calculation of gradient direction can result in large update values proximal to the source position. In the present study, final parameter updates were performed with gradient direction scaled in a pseudo-hessian manner as follows (Shin et al., 2001; Ha et al., 2009): 𝑝! !!! = 𝑝! ! − 𝛼 ∙ NRM[(diag 𝑯! + 𝛾𝐼 )!! 𝜕𝐸 ], 𝜕𝑝! (17) where 𝑝! ! is the k-th model parameter at the i-th update, 𝛼 is the step length and has a constant value, 𝑯𝒑 is the pseudo-hessian parameter defined as the square of the virtual source, 𝛾 is a damping factor, 𝐼 is a diagonal matrix, and NRM[ ] is the normalizing operation. 210 NUMERICAL TESTS We performed numerical tests with our proposed inversion algorithm, and then compared these results with those from tests performed with elastic inversion as described in eq. (5) (which is referred to conventional inversion hereafter). Briefly, numerical tests were performed with synthetic data simulating an ocean bottom multi-component seismic survey, and then were compared with the results obtained from developed inversion with elastic inversion for multi-component data. The Marmousi2 model (Martin et al., 2006) was used for modeling and inversion because this model provides Pwave velocity, S-wave velocity, and density models for elastic inversion. This model also includes various structures for benchmarking. Recently, ocean bottom multi-component surveys have been actively used for elastic inversion studies. Therefore, we simulated various receiver numbers based on available OBC and OBS surveys. Fig. 4. True Marmousi2 model of (a) P-wave and (b) S-wave velocities. 211 Inversion simulating of an ocean bottom survey with use of a large number of receivers We compared the results of conventional inversion with those from our proposed inversion for simulating an OBC survey with the Marmousi2 model. True models of P-wave and S-wave velocities which were used for inversion are shown in Fig. 4. The initial P-wave and S-wave velocity models for inversion (shown in Fig. 5) were generated with the smooth2 function of Seismic Unix. The density model was assumed to be homogeneous. Fig. 5. The smoothed initial model of (a) P-wave and (b) S-wave velocities. 212 To obtain modeling data, pressure sources at the water surface were defined at 52 points which were spaced at 320-m intervals. In addition, receivers simulating the OBC survey recorded horizontal and vertical displacements for 850 points spaced at 20-m intervals. Survey geometry is shown in Fig. 6. Fig. 6. Survey geometry in the marine Marmousi2 model. Fig. 7. Updated P-wave velocity models from (a) conventional elastic inversion and (b) our proposed inversion after 200 iterations. 213 The source waveform was obtained from a first derivative Gaussian function with a 15 Hz cutoff frequency. The maximum recording time was 6 s and the sampling interval was 0.002 s. After performing an inversion with the above parameters and survey geometry, P-wave velocity models after 200 iterations were generated (Fig. 7). S-wave velocity models are shown in Fig. 8. Fig. 8. Updated S-wave velocity models from (a) conventional inversion and (b) our proposed inversion after 200 iterations. The results from the conventional inversion technique and our proposed inversion technique both accurately show a thin, high-velocity layer (Figs. 7 and 8). Both approaches also clearly provide various structures such as a gas layer, fault, and turtle back. However, the images obtained with the proposed inversion method are more accurate. The differences between these two methods are apparent from the depth profiles obtained at three points which include various structures such as gas sand, salt layer, fault, turtle back, oil cap, and strata (Fig. 9). 214 Fig. 9. The three points evaluated in depth profiles of the Marmousi2 model. The depth profile of the inverted velocities at the 3-km point is shown in Fig. 10. The P-wave velocity trends observed in this profiles accurately indicate details in the structures. However, the proposed inversion method provides more exact information in the area within the dotted circles in regard to P-wave (a) and S-wave (b) velocities. Fig. 10. Depth profiles of (a) P-wave and (b) S-wave velocities at the 3-km point. 215 The depth profiles at the 10.4-km point are shown in Fig. 11. Similar to the results obtained at the 3-km point, the proposed inversion method also provides more accurate information for both P-wave (a) and S-wave (b) velocities. In particular, the profile for the turtle-back structure below the high-velocity layer and low-velocity information for the oil trap were more exact (see regions within the dotted circles). Fig. 11. Depth profiles of (a) P-wave and (b) S-wave velocities at the 10.4-km point. Fig. 12 presents depth profiles obtained for the 14-km point. The proposed inversion method provides more accurate velocity information in the salt layer and below. 216 Fig. 12. Depth profiles of (a) P-wave and (b) S-wave velocities at the 14.0-km point. Thus, the results obtained from the proposed inversion method to simulate the ocean bottom and multi-component surveys provides more accurate velocity information at greater depths. In addition, the results obtained are similar to the true velocity recorded in the area below the highvelocity layer. Full waveform inversion simulating OBS survey To test the proposed inversion method, we simulated an OBS survey which uses node-type receivers. In the present study, the number of receivers was small, and the intervals were wider than those used in the OBC survey. Specifically, 22 OBS receivers were simulated at the ocean bottom at 400-m intervals. The corresponding survey geometry is shown in Fig. 13. The other parameters for inversion were the same as described in the previous example above. 217 Fig. 13. Survey geometry in the Marmousi2 model with an OBS receiver. Results from the updated conventional elastic inversion method (a) and the proposed inversion method are presented in Figs. 14 and 15, respectively. Both accuracy and the continuity of the structures are worse than the previous test, although the proposed inversion method provided clear P-wave and S-wave velocity information in the updated model. Fig. 14. Updated P-wave velocity models from (a) conventional elastic inversion and (b) proposed inversion methods after 200 iterations. 218 Fig. 15. Updated S-wave velocity models from (a) conventional elastic inversion and (b) proposed inversion methods after 200 iterations. To identify specific differences between the two methods, depth profiles were generated at three points as described above. The depth profile presented in Fig. 16 for the 3-km point shows that both inversion methods can find trends in P-wave and S-wave velocities with a small number of receivers. However, use of the proposed inversion method did provide more accurate results near the high-velocity layers (indicated in dotted circles). 219 Fig. 16. Depth profiles of (a) P-wave and (b) S-wave velocities at the 3-km point with OBS receivers. Fig. 17 presents the depth profiles obtained for the 10.4-km point. The proposed inversion method provided more accurate information of the turtle back and complex structures present (shown in dotted circles). In particular, low velocity is observed below a continuous high-low velocity area in the circled regions in (a) and (b). 220 Fig. 17. Depth profiles of (a) P-wave and (b) S-wave velocities at the 10.4-km point with OBS receivers. The depth profiles at 14.0 km are shown in Fig. 18. The proposed inversion method provided accurate velocity information for the salt layer in regard to P-wave and S-wave velocities, despite use of a small number of receivers. Moreover, at greater depths below the salt layer, the depth profile closely represented the true model. 221 Fig. 18. Depth profiles of (a) P-wave and (b) S-wave velocities at the 14.0-km point with OBS receivers. Thus, the inversion results simulating an OBS multi-component survey exhibited less continuity than the previous example. However, the proposed inversion method did provide more accurate information at greater depths, particularly in regard to velocity. CONCLUSIONS We developed a full waveform inversion method which employs Pwavefield decomposition by using numerical, multi-component seismic survey data. In elastic inversion, various wave modes such as PP, PS, SP, and SS act as noise. Therefore, we hypothesized that elastic inversion would be improved by applying P-wavefield decomposition. Furthermore, we proposed the application of a back-propagation algorithm for efficient inversion to be possible by employing a stress tensor in the objective function. This proposed inversion method was subsequently evaluated in numerical tests with synthetic data from a Marmousi2 model. In the test- 222 simulating OBC survey, our proposed inversion method accurately predicted a thin, high-velocity layer, strata, and gas layer in updated velocity models. In particular, clear structures and accurate velocity information were predicted at greater depths compared with the elastic inversion method. These improved results are attributed to the use of P-wavefield decomposition, which prevents the corruption of data by S-waves at greater depths. In the numerical test to simulate an OBS survey, both inversion methods provided unclear images with use of a small number of receivers. However, our proposed inversion method did provide more accurate boundary information and velocity of the salt layer. Thus, inversion with Pwavefield decomposition provided better results with a small number of receivers. Finally, P-wavefield decomposition in this study can be considered as a pressure field, whereby pressure data obtained at the ocean bottom can be applied for inversion. As a result, these pressure data can provide information regarding elastic properties with use of our proposed inversion method. ACKNOWLEDGMENTS This research was supported by the Basic Research Project (18-3312) of the Korea Institute of Geoscience and Mineral Resources (KIGAM) which is funded by the Ministry of Science, ICT and Future Planning, and the Korea Institute of Energy Technology Evaluation and Planning (KETEP) with a grant funded by the Ministry of Trade, Industry and Energy, Republic of Korea (No. 20182510102470). REFERENCES Amundsen, L., Ikelle, L. and Martin, J., 1998. Multi attenuation and P/S splitting of OBC data: A heterogeneous sea floor. Expanded Abstr., 68th Ann. Internat. SEG Mtg., New Orleans: 722-725. Claerbout, J.F. and Muir, F., 1973. Robust modeling with erratic data. Geophysics, 38: 826-844. Crase, E., Pica, A., Nobel, M., McDonald, J. and Tarantola, A., 1990. Robust elastic nonlinear waveform inversion: Application to real data. Geophysics, 55: 527538. Choi, Y., Min, D. and Shin, C., 2008. Frequency-domain full waveform inversion using the new pseudo-Hessian matrix: Experience of elastic Marmousi-2 synthetic data. Bull. Seismol. Soc. Am., 98: 2402-2415. Chung, W., Pyun, S., Bae, H., Shin, C. and Marfurt, K.J., 2012. Implementation of elastic reverse-time migration using wavefield separation in the frequency domain. Geophys. J. Internat., 189: 1611-1625. 223 Dankbaar, M., 1987. Vertical seismic profiling - separation of P- and S-waves. Geophys. Prosp., 35: 803-814. Dellinger, J. and Etgen, J., 1990. Wave-field separation in two-dimensional anisotropic media. Geophysics, 55: 914-919. Gauthier, O., Virieux, J. and Tarantola, A., 1986. Two-dimensional nonlinear inversion of seismic waveforms: Numerical results. Geophysics, 51: 1387-1403. Ha, T., Chung, W. and Shin, C., 2009. Waveform inversion using a back-propagation algorithm and a Huber function norm. Geophysics, 74(3): R15-R24. Ha, J., Shin, S., Shin, C. and Chung, W., 2015. Efficient elastic reverse-time migration for the decomposed P-wavefield using stress tensor in the time domain. J. Appl. Geophys., 116: 121-134. Jeong, W., Lee, H. and Min, D., 2012. Full waveform inversion strategy for density in the frequency domain. Geophys. J. Internat., 188: 1221-1242. Kolb, P., Collino, F. and Lailly, P., 1986. Pre-stack inversion of a 1-D medium. Proc. IEEE, 74: 498-508. Lee, H., Koo, J., Min, D., Kwon, B. and Yoo, H., 2010. Frequency-domain elastic full waveform inversion for VTI media. Geophys. J. Internat., 183: 884-904. Lailly, P., 1983. The seismic inverse problem as a sequence of before stack migration. Conference on Inverse Scattering, Theory and Application, Soc. Industr. Appl. Mathemat., 206-220. Martin, G., Wiley, R. and Marfurt, K.J., 2006. Marmousi2: An elastic upgrade for Marmousi. The Leading Edge, 25: 156-166. Mora, P., 1987, Nonlinear two-dimensional elastic inversion of multi-offset seismic data. Geophysics, 52: 1211-1228. Operto, S., Gholami, Y., Prieux, V., Ribodetti, A., Brossier, R.. Metivier, L. and Virieux, J., 2013. A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice. The Leading Edge, 32: 10401054. Operto, S. and Miniussi, A., 2018. On the role of density and attenuation in threedimensional multiparameter viscoacoustic VTI frequency-domain FWI: An OBC case study from the North Sea. Geophys. J. Internat., 213: 2037-2059. Pan, W., Innanen, K.A., Yu, G. and Li, J., 2019. Interparameter trade-off quantification for isotropic-elastic full-waveform inversion with various model parameterizations. Geophysics, 84(2): R185-R206. Pratt, R.G., 1990. Inverse theory applied to multi-source cross-hole tomography. part 2: elastic wave-equation method. Geophys. Prosp., 38: 311-329. Pratt, R.G., Shin, C. and Hicks, G.J., 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophys. J. Internat., 133: 341362. Sheen, D., Tuncay, K., Baag, C.E. and Ortoleva, P.J., 2006. Time domain Gauss-Newton seismic waveform inversion in elastic media. Geophys. J. Internat., 167: 13731384. Shin, C., Yoon, K., Marfurt, K.J., Park, K., Yang, D., Lim, H.Y., Chung, S. and Shin, S., 2001. Efficient calculation of a partial-derivative wavefield using reciprocity for seismic imaging and inversion. Geophysics, 66: 1856-1863. Shipp, R. and Singh, S., 2002. Two-dimensional full waveform inversion of wideaperture marine seismic streamer data. Geophys. J. Internat., 151: 325-344. Sun, R., 1999. Separating P- and S-waves in a prestack 2-dimensional elastic seismogram. Extended Abstr., 61st EAGE Conf., Helsinki: 6-23. Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49: 1259-1266. 224 Yan, J. and Sava, P., 2008. Isotropic angle-domain elastic reverse-time migration. Geophysics, 73(6), S229-S239. Yan, J. and Sava, P., 2009. Elastic wave-mode separation for VTI media. Geophysics, 74(5): WB19-WB32. Yan, J., 2010. Wave-mode Separation for Elastic Imaging in Transversely Isotropic Media. Ph.D. thesis, Colorado School of Mines, Golden, CO. Zhang, Q. and McMechan, G.A., 2011. Common-image gathers in the incident phaseangle domain from reverse time migration in 2D elastic VTI media. Geophysics, 76(6): S197-S206.