Least-squares Reverse Time Migration Research Articles

Inversion-based imaging: from LSRTM to FWI imaging

Summary Least Squares RTM (LSRTM) is a powerful inversion-based imaging algorithm which minimizes the data misfit between observed seismic recordings and forward-modelled synthetic data. The algorithm, which can be implemented in either data or image domains, carries a fundamental limitation because it is based on a linear inversion theory which cannot accommodate velocity refinement as part of its model update process. Successful application of LSRTM therefore requires highly accurate velocity information, and if the velocity model is in significant error, modelled events will not be aligned kinematically with the observed data, and the algorithm will tend to produce unsatisfactory results. FWI is another inversion-based algorithm that enjoys widespread industry use. Unlike LSRTM, FWI poses its inverse problem within a non-linear framework whereby it updates the velocity model and associated wave paths throughout its iterative process, gradually aligning modelled events with observed events. With the recent convergence of FWI and LSRTM methodologies, FWI is not only being used as a velocity update tool, but also as a direct imaging tool, thereby achieving two key imaging goals, namely refining the velocity model and deriving a better-quality seismic image. The latter process, which is known as ‘FWI imaging’, has recently been gaining a lot of industry attention as it offers the possibility of high-quality imaging along with workflow simplification. In this article we will compare and contrast LSRTM and FWI. We conclude that the process of generating the FWI-imaging essentially amounts to nonlinear, data-domain inversion. This recognition facilitates a ready comparison against the data-domain form of LSRTM, the latter being a linear, data-domain inversion.

Least-Squares Reverse Time Migration Using the Gradient Preconditioning Based on Transmitted Wave Energy

A gradient preconditioning approach based on transmitted wave energy for least-squares reverse time migration (LSRTM) is proposed in this study. The gradient is preconditioned by using the energy of “approximate transmission wavefield,” which is calculated based on the non-reflecting acoustic equation. The proposed method can effectively avoid a huge amount of calculation and storage required by the Hessian matrix or approximated Hessian matrix and can overcome the influence of reflected waves, multiples, and other wavefields on the gradient in gradient preconditioning based on seismic wave energy (GPSWE). The numerical experiments, compared with that using GPSWE, show that LSRTM using the gradient preconditioning based on transmitted wave energy (GPTWE) can significantly improve the imaging accuracy of deep target and accelerate the convergence rate without trivial increased calculation.

Imaging Complex Subsurface Structures for Geothermal Exploration at Pirouette Mountain and Eleven-Mile Canyon in Nevada

Accurate imaging of subsurface complex structures with faults is crucial for geothermal exploration because faults are generally the primary conduit of hydrothermal flow. It is very challenging to image geothermal exploration areas because of complex geologic structures with various faults and noisy surface seismic data with strong and coherent ground-roll noise. In addition, fracture zones and most geologic formations behave as anisotropic media for seismic-wave propagation. Properly suppressing ground-roll noise and accounting for subsurface anisotropic properties are essential for high-resolution imaging of subsurface structures and faults for geothermal exploration. We develop a novel wavenumber-adaptive bandpass filter to suppress the ground-roll noise without affecting useful seismic signals. This filter adaptively exploits both characteristics of the lower frequency and the smaller velocity of the ground-roll noise than those of the signals. Consequently, this filter can effectively differentiate the ground-roll noise from the signal. We use our novel filter to attenuate the ground-roll noise in seismic data along five survey lines acquired by the U.S. Navy Geothermal Program Office at Pirouette Mountain and Eleven-Mile Canyon in Nevada, United States. We then apply our novel anisotropic least-squares reverse-time migration algorithm to the resulting data for imaging subsurface structures at the Pirouette Mountain and Eleven-Mile Canyon geothermal exploration areas. The migration method employs an efficient implicit wavefield-separation scheme to reduce image artifacts and improve the image quality. Our results demonstrate that our wavenumber-adaptive bandpass filtering method successfully suppresses the strong and coherent ground-roll noise in the land seismic data, and our anisotropic least-squares reverse-time migration produces high-resolution subsurface images of Pirouette Mountain and Eleven-Mile Canyon, facilitating accurate fault interpretation for geothermal exploration.

Efficient iterative prestack least-squares reverse time migration through surface-offset gather compression

Extended least-squares inversion is superior to stack-based least-squares inversion for imaging the subsurface because it can better account for amplitude-variation-with-offset (AVO) effects as well as residual moveout (RMO) effects induced by erroneous velocity models. Surface-offset extensions have proven to be a robust alternative to angle gathers as well as subsurface extensions when applied to narrow-azimuth data acquisitions, especially when using erroneous velocity models. As such, least-squares reverse time migration (LSRTM) applied to surface-offset gathers (SOGs) obtains accurate surface-offset-dependent estimates of the reflectivity with better AVO behavior, while respecting the curvatures of the events in the gathers. Nevertheless, the computational expense incurred by SOG demigration generally renders this process unfeasible in many practical situations. We have exploited a compression scheme for SOGs that captures AVO and some RMO effects to improve efficiency of extended LSRTM. The decompression operator commutes with the demigration operator; therefore, gathers compressed in the model domain may be decompressed in the data domain. This obviates the need to demigrate all SOGs, requiring only the demigration of a few compressed gathers. We determine the accuracy of this compression in the model and data domains with a synthetic 2D data set. We then use our model-compression/data-decompression scheme to SOG-extended iterative LSRTM for two field data examples from offshore Brazil. These examples demonstrate that our compression can capture most AVO and some RMO information accurately, while greatly improving efficiency in many practical scenarios.

Least-squares reverse time migration using convolutional neural networks

Least-squares reverse time migration (LSRTM) has the potential to reconstruct a high-resolution image of subsurface reflectivity. However, the current data-domain LSRTM approach, which iteratively updates the subsurface reflectivity by minimizing the data residuals, is a computationally expensive task. To alleviate this problem and improve imaging quality, we have developed an LSRTM approach using convolutional neural networks (CNNs), which is referred to as CNN-LSRTM. Specifically, the LSRTM problem can be implemented via a gradient-like iterative scheme, in which the updating component in each iteration is learned via a CNN model. To make the most of the observation data and migration velocity model at hand, we use the common-source reverse time migration (RTM) image, the stacked RTM image, and the migration velocity model rather than only the stacked RTM image as the input data of CNN. We have successfully trained the constructed CNN model on the training data sets with a total of 5000 randomly layered and fault models. Based on the well-trained CNN model, we have proved that our approach can efficiently recover the high-resolution reflection image for the layered, fault, and overthrust models. Through a marine field data experiment, we can determine the benefit of our constructed CNN model in terms of its computational efficiency. In addition, we analyze the influence of input data of the constructed CNN model on the reconstruction quality of the reflection image.

From acoustic to elastic inverse extended Born modeling: A first insight in the marine environment

Elastic least-squares reverse time migration is a state-of-the-art linear imaging technique used to retrieve high-resolution quantitative subsurface images. Its successful application requires many migration/modeling cycles. To accelerate the convergence rate, various pseudoinverse Born operators have been proposed, providing quantitative results within a single iteration, while having roughly the same computational cost as reverse time migration. However, these are based on the acoustic approximation, leading to possible inaccurate amplitude predictions as well as the ignorance of S-wave effects. To solve this problem, we have extended the pseudoinverse Born operator from acoustic to elastic media to account for the elastic amplitudes of PP reflections and provide an estimate of physical density P- and S-wave impedance models. We restrict the extension to marine environments, with the recording of pressure waves at the receiver positions. First, we replace the acoustic Green’s functions by their elastic version, without modifying the structure of the original pseudoinverse Born operator. We then apply a Radon transform to the results of the first step to calculate the angle-dependent response. Finally, we simultaneously invert for the physical parameters using a weighted least-squares method. Through numerical experiments, we first examine the consequences of the acoustic approximation on elastic data, leading to inaccurate parameter inversion as well as to artificial reflector inclusion. Then, we determine that our method can simultaneously invert for elastic parameters in the presence of complex uncorrelated structures, inaccurate background models, and Gaussian noisy data.

The importance of including density in multiparameter asymptotic linearized direct waveform inversion: a case study from the Eastern Nankai Trough

SUMMARY Iterative least-squares reverse time migration (LSRTM) is the state-of-the-art linearized waveform inversion method to obtain quantitative subsurface parameters. The main drawback of such an iterative imaging scheme is the significant computational expense of many modelling/adjoint cycles through iterations. In the context of the extended domain, an interesting alternative to LSRTM is the asymptotic linearized direct waveform inversion, providing quantitative results with only a single iteration. This approach was first proposed for constant-density acoustics and recently extended to the variable-density case. The former is based on the application of the asymptotic inverse Born operator, whereas the latter has two more extra steps: building an angle-dependent response of the asymptotic inverse Born operator and then solving a weighted least-squares approach for simultaneous inversion of two acoustic parameters. To examine the importance of accounting for density variations, we compare constant- and variable-density linearized direct waveform inversion techniques applied to a marine real data set from the Eastern Nankai Trough, offshore Japan. The inversion results confirm the efficiency of the asymptotic linearized direct waveform inversion in estimating quantitative parameters within a single iteration. The variable-density direct inversion yields subsurface images that (1) exhibit a superior resolution and (2) better reconstruct the field data than does the constant-density approach, even if the data set does not contain large enough surface offset to fully decompose velocity and density perturbations.

Prestack correlative elastic least-squares reverse time migration based on wavefield decomposition

Standard elastic least-squares reverse time migration (ELSRTM) can iteratively improve the quality of stacked reflection images of P- and S-wave impedance by minimizing the L2 norm of waveform residual. However, the resulting images are sensitive to the velocity errors and amplitude errors in the observed data. Moreover, the crosstalk artifacts created by the different wave modes may degrade the imaging quality. To mitigate these problems, a correlative ELSRTM approach based on wavefield decomposition and prestack implementation is developed. This approach includes two key points. The first is that the proposed approach aims at finding the optimal reflection images for each shot recording via minimizing the normalized zero-lag cross-correlation objective function. The final reflection images can be produced by summing the optimal reflection images for all shot recordings. The second is that we develop a wavefield decomposition scheme, which is compatible with the correlative ELSRTM approach, as a gradient precondition to reduce crosstalk artifacts and accelerate the convergence. Numerical examples using several 2D experiments have determined that the proposed ELSRTM approach can provide higher quality images with fewer crosstalk artifacts, compared to the correlative ELSRTM approach without wavefield decomposition. Furthermore, it is less sensitive to velocity errors than the poststack correlative ELSRTM approach.

Elastic least-squares reverse time migration based on decoupled wave equations

Elastic reverse time migration (ERTM) is developed for better characterization of complex structures by imaging multicomponent seismic data. However, conventional ERTM is subject to limitations such as finite recording aperture, limited bandwidth, and imperfect illumination. Elastic least-squares reverse time migration (ELSRTM) can improve imaging accuracy gradually with iterations by minimizing the residuals between the observed and calculated multicomponent data. Conventional ELSRTM suffers from crosstalk artifacts caused by coupled elastic wavefields with different wave modes. Decomposing the coupled elastic wavefields into pure P- and S-waves is an effective method to suppress these crosstalk artifacts. Considering the trade-off between calculation accuracy and efficiency, we have developed a new ELSRTM scheme for isotropic media based on decoupled wave equations to suppress these wave mode-related crosstalk artifacts in the images of conventional ELSRTM. Pure wavefields are obtained by solving the decoupled wave equations using the finite-difference method in our new ELSRTM method. We also derive new decoupled adjoint-state wave equations that are suitable for the elastic velocity-stress equations in isotropic media. Furthermore, we use the gradient equations based on pure wavefields to update the reflectivity images. Synthetic examples demonstrate that our new ELSRTM method can generate images that better represent the subsurface when compared with conventional ERTM and conventional ELSRTM.

Sparse Constrained Least-Squares Reverse Time Migration Based on Kirchhoff Approximation

Least-squares reverse time migration (LSRTM) is powerful for imaging complex geological structures. Most researches are based on Born modeling operator with the assumption of small perturbation. However, studies have shown that LSRTM based on Kirchhoff approximation performs better; in particular, it generates a more explicit reflected subsurface and fits large offset data well. Moreover, minimizing the difference between predicted and observed data in a least-squares sense leads to an average solution with relatively low quality. This study applies L1-norm regularization to LSRTM (L1-LSRTM) based on Kirchhoff approximation to compensate for the shortcomings of conventional LSRTM, which obtains a better reflectivity image and gets the residual and resolution in balance. Several numerical examples demonstrate that our method can effectively mitigate the deficiencies of conventional LSRTM and provide a higher resolution image profile.

Hessian-based reflection waveform inversion: Theory and applications to synthetic and field data

Wave-equation reflection waveform inversion (RWI) is a promising method to reconstruct the background velocity model with reflection data. But it is difficult to precondition this highly nonlinear inverse problem for efficient convergence and reliable model updating. In the context of full-waveform inversion (FWI), the second-order derivative of the objective function with respect to the model parameters (i.e., the Hessian operator) can be used with the second-order optimization to overcome similar difficulties. To our knowledge, however, there are few such studies for RWI. This motivates us to derive and investigate the Hessian operator using the Born approximation in the context of RWI. For the same experimental setup with a toy model, we observe that the Hessian matrix of RWI exhibits a more compact structure than its counterpart in the context of FWI because they are constructed with the sensitivity kernels of very different spatial distribution features. Accordingly, we have developed a truncated Gauss-Newton RWI method consisting of two nested loops: An outer loop updates the background model in the manner of nonlinear waveform inversion, and an inner loop aims to iteratively estimate the background model update in a matrix-free manner by solving the Gauss-Newton equation. The workflow consists of alternating estimation of the perturbation and the background velocity model components, of which the former is done with a linear least-squares reverse time migration using the conjugate-gradient method. A synthetic example on the overthrust model and a field data example demonstrate that the Hessian-based RWI can significantly improve the estimation of intermediate-wavenumber components and promote velocity model building and seismic imaging with precritical reflection data.

Full wavefield least-squares reverse time migration

Waveform inversion based on least-squares reverse time migration (LSRTM) usually involves Born modeling, which models the primary-only data. As a result, the inversion process handles only primaries and corresponding multiple elimination preprocessing of the input data is required prior to imaging and inversion. Otherwise, multiples left in the input data are mapped as false reflectors, also known as crosstalk, in the final image. At the same time, the developed full-wavefield migration (FWM) methodology can handle internal multiples in an inversion-based imaging process. However, because it is based on the framework of the one-way wave equation, it cannot image dips close to and beyond 90°. Therefore, we have aimed at upgrading the LSRTM framework by bringing in the functionality of FWM to handle internal multiples. We inject the secondary source term, as used in the original formulation of FWM, to define a wavefield relationship that allows modeling multiple scattering via reflectivity. The secondary source term is based on the estimated reflectivity and can be injected into the pressure component when simulating the two-way wave equation using finite-difference modeling. We use this modified forward model for estimating the reflectivity model in an FWM-type manner and validate the method on synthetic and field data containing visible internal multiples.

Reflection angle/azimuth-dependent least-squares reverse time migration

Seismic images under complex overburdens such as salt are strongly affected by illumination variations due to overburden velocity variations and imperfect acquisition geometries, making it difficult to obtain reliable image amplitudes. Least-squares reverse time migration (LSRTM) addresses these issues by formulating full wave-equation imaging as a linear inverse problem and solving for a reflectivity model that explains the recorded seismic data. Because subsurface reflection coefficients depend on the incident angle, and possibly on the azimuth, quantitative interpretation under complex overburdens requires LSRTM with output in terms of image gathers, e.g., as a function of the reflection angle or angle and azimuth. We have developed a reflection angle- or angle/azimuth-dependent LSRTM method aimed at obtaining physically meaningful image amplitudes interpretable in terms of angle- or angle/azimuth-dependent reflection coefficients. The method is formulated as a linear inverse problem solved iteratively with the conjugate gradient method. It requires an adjoint pair of linear operators for reflection angle/azimuth-dependent migration and demigration based on full wave-equation propagation. We implement these operators in an efficient way by using a mapping approach between migrated shot gathers and subsurface reflection angle/azimuth gathers. To accelerate convergence of the iterative inversion, we apply image-domain preconditioning operators computed from a single de-remigration step. An angle continuity constraint and a structural dip constraint, implemented via shaping regularization, are used to stabilize the solution in the presence of limited illumination and to control the effects of coherent noise. We examine the method on a synthetic data example and on a wide-azimuth streamer data set from the Gulf of Mexico, where we find that angle/azimuth-dependent LSRTM can achieve significant uplift in subsalt image quality, with overburden- and acquisition-related illumination variation effects on angle/azimuth-dependent image amplitudes largely removed.

Angle-domain common-image gathers from plane-wave least-squares reverse time migration

Angle-domain common-image gathers (ADCIGs) that can be used for migration velocity analysis and amplitude-versus-angle analysis are important for seismic exploration. However, because of the limited acquisition geometry and seismic frequency band, ADCIGs extracted by reverse time migration (RTM) suffer from illumination gaps, migration artifacts, and low resolution. We have developed a reflection angle-domain pseudoextended plane-wave least-squares RTM method for obtaining high-quality ADCIGs. We build the mapping relations between the ADCIGs and the plane-wave sections using an angle-domain pseudoextended Born modeling operator and an adjoint operator, based on which we formulate the extraction of ADCIGs as an inverse problem. The inverse problem is iteratively solved by a preconditioned stochastic conjugate-gradient method, allowing for reduction in computational cost by migrating only a subset instead of the whole data set and improving the image quality thanks to preconditioners. Numerical tests on synthetic and field data verify that our method can compensate for illumination gaps, suppress migration artifacts, and improve resolution of the ADCIGs and the stacked images. Therefore, compared to RTM, our method provides a more reliable input for migration velocity analysis and amplitude-versus-angle analysis. Moreover, it also provides much better stacked images for seismic interpretation.

Q least-squares reverse time migration based on the first-order viscoacoustic quasidifferential equations

Seismic wave attenuation caused by subsurface viscoelasticity reduces the quality of migration and the reliability of interpretation. A variety of Q-compensated migration methods have been developed based on the second-order viscoacoustic quasidifferential equations. However, these second-order wave-equation-based methods are difficult to handle with density perturbation and surface topography. In addition, the staggered grid scheme, which has an advantage over the collocated grid scheme because of its reduced numerical dispersion and enhanced stability, works in first-order wave-equation-based methods. We have developed a Q least-squares reverse time migration method based on the first-order viscoacoustic quasidifferential equations by deriving Q-compensated forward-propagated operators, Q-compensated adjoint operators, and Q-attenuated Born modeling operators. In addition, our method using curvilinear grids is available even when the attenuating medium has surface topography and can conduct Q-compensated migration with density perturbation. The results of numerical tests on two synthetic and one field data set indicate that our method improves the imaging quality with iterations and produces better imaging results with clearer structures, higher signal-to-noise ratio, higher resolution, and more balanced amplitude by correcting the energy loss and phase distortion caused by Q attenuation. It also suppresses the scattering and diffracted noise caused by the surface topography.

Shear reflectivity compensation in full-waveform inversion using least-squares reverse-time migration

SUMMARY The computational cost of elastic-waveform inversion is too high for inverting PP reflections, while using acoustic full-waveform inversion (FWI) is inaccurate because it does not depend on the shear modulus/velocity/impedance that affects elastic PP wavefield amplitudes. To solve this problem, we develop a waveform inversion method that uses acoustic least-squares reverse-time migration (LSRTM) to compensate the shear reflectivity for acoustic FWI. Our method is based on the quasi-elastic-wave equation developed by Chapman et al. (2014). The quasi-elastic-wave equation uses a linearized acoustic-wave equation with shear modulus μ as a virtual source to correct the acoustic PP wavefield amplitudes toward elastic ones. Our waveform inversion method inverts for elastic parameters by minimizing the L2 norm of the difference between recorded and predicted PP reflections modelled using the quasi-elastic-wave equation. Numerical tests on synthetic and field data show that our method can properly handle the amplitudes of elastic PP reflections and provides an accurate estimate of the P- and S-wave velocities/impedances and, in some cases, the density. The method does not need the computationally expensive numerical solution to the elastic-wave equation. It also gives a better estimate of elastic parameters than a pure LSRTM method for elastic PP reflections.

Least-squares reverse-time migration with sparsity constraints

Abstract Least-squares reverse-time migration (RTM) works with an inverse operation, rather than an adjoint operation in a conventional RTM, and thus produces an image with a higher resolution and more balanced amplitude than the conventional RTM image. However, least-squares RTM introduces two side effects: sidelobes around reflectors and high-wavenumber migration artifacts. These side effects are caused mainly by the limited bandwidth of seismic data, the limited coverage of receiver arrays and the inaccuracy of the modeling kernel. To mitigate these side effects and to further boost resolution, we employed two sparsity constraints in the least-squares inversion operation, namely the Cauchy and L1-norm constraints. For solving the Cauchy-constrained least-squares RTM, we used a preconditioned nonlinear conjugate-gradient method. For solving the L1-norm constrained least-squares RTM, we modified the iterative soft thresholding method. While adopting these two solution methods, the Cauchy-constrained least-squares RTM converged faster than the L1-norm constrained least-squares RTM. Application examples with synthetic data and laboratory modeling data demonstrated that the constrained least-squares RTM methods can mitigate the side effects and promote image resolution.

Shenzi OBN: An imaging step change

The story of seismic imaging over BHP's Shenzi Gulf of Mexico production field follows the history of offshore seismic imaging, from 2D to 3D narrow-azimuth streamer acquisition and to its leading the wide-azimuth movement with the Shenzi rich-azimuth (RAZ) survey. Each RAZ reprocessing project over the last 15 years applied the latest processing technology, culminating in hundreds of scenario tests to refine the salt model, but eventually the RAZ data reached a technical limit. A new ocean-bottom-node (OBN) survey acquired in 2020 has produced a step-change improvement over the legacy RAZ image. The uplift can be attributed to several factors. First, an OBN feasibility and survey design study demonstrated that a core of dense nodes combined with sparse nodes would improve the accuracy and resolution of the full-waveform inversion (FWI) solution. Second, the OBN data acquired following the survey design and employing FWI as the main model-building tool realized the predicted improvement. The result was a substantial change to the complex salt model, verified by a salt proximity survey as well as other salt markers, and improvement in imaging over the entire field. In addition to the improvement arising from a more accurate FWI velocity model, the steep-dip imaging also benefited from the new full-azimuth and long-offset data. However, the best steep-dip and fault imaging comes from the FWI image, a direct estimation of reflectivity from the FWI velocity. As the maximum frequency used by FWI moves toward the maximum frequency of the final reverse time migration (RTM), the FWI image approaches the resolution necessary to compete as the primary interpretation volume. Its subsalt illumination surpassed that of the RTM and even the least-squares RTM volumes. These imaging improvements are providing a new understanding of the faults and stratigraphic relationships of the field.

Sparsity-promoting multiparameter pseudoinverse Born inversion in acoustic media

Least-squares reverse time migration (RTM) has become the method of choice for quantitative seismic imaging. The main drawback of such a scheme is that it requires many migration/modeling cycles. The convergence of least-squares RTM can be accelerated by using a suitable preconditioner. In the context of an extended domain in variable-density acoustic media, the pseudoinverse Born operator is the recommended preconditioner, providing quantitative results within a single iteration. This method consists of two steps: application of the pseudoinverse Born operator and inversion of two parameters using an efficient weighted least-squares approach based on the Radon transform. As expected, crosstalk artifacts are generated in the second step due to limited acquisition. We have developed a variable-density pseudoinverse Born operator constrained with the [Formula: see text]-norm for each model parameter to suppress the artifacts. The fast iterative shrinkage-thresholding algorithm is used to carry out the optimization problem. In classic iterative least-squares migration, the [Formula: see text]-norm constraints would affect the whole imaging process. Because the imaging method is split into two steps, only the Radon transform part is modified, where no wave-based operators are involved. Through numerical experiments, we verify the robustness of our method against different migration artifacts including the parameter crosstalk, interfaces with abrupt truncations, sparse shot acquisition geometry, noisy data, and high-contrast complex structures.

Acoustic-elastic coupled least-squares reverse time migration in marine environment with rugged seabed interface

ABSTRACT The severe rugged seabed interfaces bring great difficulties to seismic imaging in the marine environment. To accurately image submarine structures under the rugged seabed interface, an acoustic-elastic coupled curvilinear-coordinated least-squares reverse time migration (AE-CLSRTM) is proposed. This method is based on the coupled equation method, which uses the acoustic wave equations in seawater and the elastic wave equations in the underlying elastic medium. The pressure in the acoustic wave equations and the stresses in the elastic wave equations are transmitted steadily and continuously by using acoustic-elastic control equations at the seabed interface. To overcome the influence of the rugged seabed interface, the acoustic-elastic model is meshed into non-uniform curvilinear grids, and the corresponding mapping technique is used to transform the model with the rugged seabed interface to a horizontal one in the curvilinear coordinate system through the coordinate transformation. Based on the acoustic-elastic coupled equations in the curvilinear coordinate system, the acoustic-elastic coupled LSRTM algorithm in the rugged seafloor structure is realised. The numerical examples on a simple model and an actual area model show that the proposed LSRTM method can obtain the accurate imaging results of submarine structures in both P- and S-velocity components.