--- title: Removing density effects in LS-RTM with a low-rank matched filter runninghead: matching density by a low-rank matched filter author: | Mengmeng Yang^1\*^, Rajiv Kumar^1^, Rongrong Wang^2^ and Felix J. Herrmann^1^\ 1 Seismic Laboratory for Imaging and Modeling (SLIM), Georgia Institute of Technology \ 2 Michigan State University bibliography: - myang.bib --- ## Abstract Least-squares reverse-time migration faces difficulties when it inverts the data containing strong components related to density variation with velocity-only Born modeling operator. The strong density perturbation will be inverted as strong dummy velocity perturbations, which influence the amplitudes and phase of the velocity perturbations in the inverted model. The traditional method is to invert the additional density variations by developing Born operator with respect to both density and velocity or modify the image condition. In this work, we develop a matched-filter based LS-RTM for velocity-only Born modeling operator, which removes the artifacts in the imaging created by the strong density variation. This method doesn't call for extra work of finite difference stencil and is more general. In the experiment part, we use a complex discontinuous layered medium with strong density variations at the ocean bottom, and show the efficacy of the propose formulation. ## Introduction Least-squares reverse-time migration (LS-RTM, @guitton2006least; @dai2012multi; @plessix2002amplitude) tries to fit observed reflection data in a least-squares sense to overcome RTM's shortcomings in producing high resolution and high-fidelity amplitudes. In other words, LS-RTM attempts to invert the linearized Born modeling operator iteratively whereas RTM directly treats the adjoint of that operator as its inverse. So far, LS-RTM has demonstrated an ability to produce high-resolution images in combination with an efficient computational framework [@herrmann2015EAGEfom; @yang2016time], overcoming drawbacks of overfitting artifacts [@herrmann2012efficient] caused by minimizing the ``\ell_2`` norm. In addition to the developments listed above, people working on LS-RTM made lots of progress on incorporating multi-parameters (elastic parameters [@duan2016elastic], visco-acoustic parameters [@dutta2014attenuation] and so on) for complex geological structures. The corresponding Born modeling operators are linearized with respect to these elastic parameters, which allow us to mimic the elastic wave-propagation effects during the inversion. While important progress has been made handling elastic effects--- e.g. by grouping subsets of elastic parameters that give rise to different radiation patterns [@operto2013guided]--- working with multiple elastic parameters remains challenging. Among the different parameters that rule the leading order behaviour of wave propagation, we count velocity and density as the most important pair [@beylkin1987multiparameter]. The products of these two, i.e., the seismic impedance, determines the amplitudes of the seismic waves to leading order. Perturbations in density generate reflection events even for a constant velocity model. This means that if we invert data generated by strong density perturbations with a Born modeling that accounts for velocity changes only we can expect strong artifacts degrading the quality of migrated images [@przebindowska2012role, @plessix2013multiparameter]. There are two main reasons for these artifacts: *(i)* the wavefields scattered by the velocity and density parameters exhibit similar behaviours for some scattering angles [@operto2013guided], and *(ii)* if we only invert for velocity perturbations without incorporating the true density perturbations in Born modeling operator then LS-RTM will try to fit the amplitudes and phase of the observed seismic data in terms of velocities only [@bai2014simultaneous]. This can lead to dummy reflection events in the LS-RTM along with incorrect amplitudes and phase distortions of the true reflectivity yielded by the velocity perturbations. In this work, we propose to use a matched-filter approach to remove artifacts caused by strong unmodeled density perturbations in the context of imaging with surface-related multiples. Specifically, we are interested in imaging based on linearized inversions that derive from velocity-only acoustic Born modeling that handles strong water-bottom multiples generated by strong density changes at the ocean bottom. We find that the proposed matched-filter, when organized as a matrix, exhibits low-rank structure. This is due to the fact that the matched-filter tries to approximate the difference between radiation patterns of velocity and a (strong ocean bottom) density contrast, which varies smoothly with offset. Inspired by this observation, we propose to simultaneously estimate velocity perturbations and a low-rank matched-filter, which maps nonlinear (observed) data that contains components related to both density and velocity perturbations into "linearized" data close to data generated by Born modeling for perturbations in velocity only. The proposed method does not require the explicit Born operator for density but does need terms that depend on a smoothly varying background density. Our paper is organized as follows. First, we form the objective function for our extended LS-RTM with a low-rank matrix constraint on the matched filter and conclude by describing a computationally efficient algorithm. Next, we evaluate the performance of the proposed approach on a quasi-layered model with faults where the density varies strongly at the ocean bottom. Finally, we show the benefits of inverting for the matched filter to correctly image both the amplitude and phase of the reflectivity for imaging problems that contain a strong density contrast at the ocean bottom. ## Methodology We start with a brief overview of LS-RTM, and then propose a matched-filter based formulation to handle strong density-related effects in imaging. As we mentioned before, LS-RTM attempts to minimize the ``\ell_2`` norm of the data residual between the observed and synthetic data by solving the following (unconstrained) optimization problem: ```math #OBJ \begin{aligned} \min_{\vector{x}} f(\vector{x}) = \sum_{i=1}^{n_f} \ & \Vert \vector{B}_i - \ctensor{J}_i(\vector{m}_0)(\vector{x}) \Vert_{F} , \end{aligned} ``` where ``\ctensor{J}_i(\vector{m}_0)`` represents the monochromatic Jacobian with respect to velocity for all shots and followed by a matrication putting monochromatic shots in its columns. The vector ``\vector{x}`` stands for the unknown velocity perturbations. Finally, the matrix ``\vector{B}_i`` is the ``i^{th}`` frequency slice of the observed (nonlinear) data in the S-R domain (source-receiver domain). In this work, we think of non-linear data as the difference between the response of the "true earth"---i.e., a "hard" model for velocities and densities and the response of a "smoothed background earth" where both velocity and density vary smoothly. The symbol ``\Vert \cdot \Vert_F`` denotes the Frobenius norm. The above equation entails a "velocity only" linearization, which is accurate in the absence of strong density variations and a good background velocity model with respect to which the linearization is carried out. Under those conditions, equation #OBJ can produce good quality migrated images, which can then be used to perform reservoir characterization. However, if the observed seismic data contains strong density effects generated by a strong ocean bottom, then the linearization undergirding equation #OBJ is no longer valid and as a consequence this may lead to a degradation in quality of migrated images [@przebindowska2012role, @plessix2013multiparameter]. One way to address this issue is to include density into equation #OBJ and re-linearize the Born modeling operator with respect to both velocity and density. While this approach is certainly a viable option, it is challenging and perhaps excessive to form a Born modeling operator that includes density, especially because it is well known that simultaneously inverting for both velocity and density is difficult because the two parameters have similar radiation patterns for certain scattering angles [@operto2013guided]. As a result, LS-RTM runs the risk to map the perturbations in density to the velocity, which results in crosstalk in the imaging [@bai2014simultaneous]. To address these issues, we propose to include a matched filter, which allows us to compensate for certain leading order density effects while inverting for the velocity perturbations only. Under the assumption that we can find such a matched filter ``\vector{M}``, we modify equation #OBJ as follows: ```math \min_{\vector{x},\vector{M}_i} f(\vector{x},\vector{M}_i) = \sum_{i=1}^{n_f} \Vert \vector{B}_i \vector{M}_i - \ctensor{J}_i(\vector{m}_0)\vector{x} \Vert_{F} , ``` where ``\vector{M}_i`` is the matched filter matrix for ``i^{th}`` frequency. As in our earlier work involving on-the-fly source estimation [@tu2013fast,@yang2016time], we use variable projections to solve for ``\vector{M}_i`` while minimizing the above objective with respect to the velocity perturbations collected in the vector ``\vector{x}``. However, contrary to finding a single time signature for the wavelet, the above matched filter involves for each frequency a full wavefield opening the risk of overfitting. To counter this problem, we control the rank ``k_i`` for each frequency. We motivate this choice by the fact that the difference in radiation patterns of velocity and the strong density contrasts at the ocean bottom vary smoothly over offset and we use this to stabilize the inversion. We now solve ```math #OBJ1 \begin{aligned} \min_{\vector{x},\vector{M}_i} f(\vector{x},\vector{M}_i) = \sum_{i=1}^{n_f} \ & \Vert \vector{B}_i \vector{M}_i - \ctensor{J}_i(\vector{m}_0)\vector{x} \Vert_{F} ,\\ & \text{s.t.} \ \text{rank}({\vector{M}_i})=k_i, \end{aligned} ``` For simplicity, we will now focus on solving the above problem for one single frequency and drop the subscript ``i`` accordingly and implicitly sum over frequencies, i.e., ``\sum_{i=1}^{n_f}`` for the remainder of this section. As we mentioned before, we solve the above problem #OBJ1 with variable projections [@golub2003separable, @tu2016source, @yang2016time]. This involves computing gradient steps with respect to ``\vector{x}`` that minimize the ``\Vert \cdot \Vert_F`` norm ---i.e, ```math #xupdate \vector{x}_{k+1} = \vector{x}_k + \mathbf{s} \nabla_{\vector{x}} f(\vector{x},\vector{M})|_{\vector{x}=\vector{x}_k, \vector{M}=\vector{M}_k}, ``` where ``{s}`` is the stepsize. As prescribed by variable projection, we solve for ``\vector{M}`` by minimizing #OBJ1 for ``\vector{x}_k`` fixed. Since rank-minimization problems are NP hard, we use its convex relaxation instead, i.e., we replace the rank constraint by a nuclear-norm constraint [@recht2010guaranteed], following the strategy proposed by @aravkin2014SISCfmd\. ## Numerical Experiments To test the performance of the proposed method, we conduct experiments on a quasi layered model with strong density contrast at the ocean bottom (Figure #f1). Both velocity and density models are ``1\mathrm{km}`` deep and ``2\mathrm{km}`` wide and the underlying grid is discretized to ``10\mathrm{m}``. The background velocity model ``\vector{m}_0`` and background density model ``\vector{\rho}_0`` are smooth version of the true velocity and density models respectively and kinematically correct. We use a Ricker wavelet centered at ``10\mathrm{hz}`` as source wavelet and record the data for 4 seconds. The ``100`` shots and ``100`` receivers are spread over the model with ``10\mathrm{m}`` spacing. To test this method on this model, we compare the image inverted from ideal linearized data related to velocity perturbation and images inverted from nonlinear data with and without the matched-filter approach. #### Figure: {#f1} ![Velocity background](fig/vel_bg.png){#f1a width=30%} ![Velocity perturbation](fig/v_pert.png){#f1b width=30%} ![True velocity](fig/vel.png){#f1c width=30%}\ ![Density background](fig/den_bg.png){#f1d width=30%} ![Density perturbation](fig/den_pert.png){#f1e width=30%} ![True density](fig/den.png){#f1f width=30%}\ : Quasi layered model with faults. (a, b, c) background velocity, the corresponding perturbation and the true veloicty. (c, d, e) Same for density model. Given the true velocity and density models and their backgrounds, we generate "observed" nonlinear data by subtracting the response of the background models for varying velocity and density from the response yielded by the true velocity and density models. We compare this response with linear data obtained by applying linearized velocity only Born scattering operator to the true velocity perturbations plotted in Figure #f1\. The monochromatic frequency slices of the linear data, nonlinear data and the matched-filter matrix ``\vector{M}`` at 5 and 10 Hz are shown in Figure [#f0]. It is clear from the figures that the estimated matched-filter varies smoothly across sources, thus, exhibits low-rank structure. These results also validate our belief that the the difference between radiation patterns of velocity and the strong density contrast exhibits smooth structure over offset. #### Figure:{#f0 .wide} ![Linear data at 5hz](fig/lineardata_5hz.png){#f0a width=30%} ![Nonlinear data at 5hz](fig/nonlineardata_5hz.png){#f0b width=30%} ![Matrix ``\vector{M}`` matching nonlinear data to linear data at 5hz](fig/M_5hz.png){#f0c width=30%}\ ![Linear data at 10hz](fig/lineardata_10hz.png){#f0d width=30%} ![Nonlinear data at 10hz](fig/nonlineardata_10hz.png){#f0e width=30%} ![Matrix ``\vector{M}`` matching nonlinear data to linear data at 10hz](fig/M_10hz.png){#f0f width=30%}\ : Visualization of the linearized data, non-linearized data and the corresponding matched-filter ``\vector{M}`` at 5 and 10 Hz. As expected, the matched-filter ``\vector{M}`` varies smoothly over offset, thus exhibits low-rank structure. Figure [#f2a] shows the idealized LS-RTM results using the linearized data (Figure #f0a). We can see that for the ideal scenario, the layer interfaces are sharp and amplitudes are in the correct range. Also it is clear that the velocity perturbations at shallower interfaces especially at the ocean bottom are much weaker than those at the deeper interfaces. Next, we invert the non-linearized data (Figure #f0b) without any matched-filter approach using the velocity-only Born modeling operator. It is evident from the inverted image (Figure #f2b) that the LS-RTM maps the strong density perturbations to the velocity perturbations, thus, creating the dummy strong reflectors at the ocean bottom. Moreover, the amplitudes and phase of the subsequent deeper reflectors are wrong. Finally, we use the proposed matched-filter approach to perform the LS-RTM (Figure [#f2c]). Using the propose method, we are able to remove the effects of the strong density perturbations at both the shallow and deeper sections. Thus, the estimated matched-filter can handle the strong density perturbation related effects while inverting for the velocity perturbation only, and can successfully remove the crosstalk created by density perturbations. #### Figure: {#f2 .wide} ![Image inverted from linear data](fig/ImageLinear.png){#f2a width=50%} ![Image inverted from nonlinear data](fig/ImageNonlinear.png){#f2b width=50%}\ ![Image inverted from nonlinear data using matched-filter approach](fig/ImageNonlinear_Lowrank_new.png){#f2c width=50%}\ : LS-RTM results using velocity-only Born modeling operator. (a) Idealized linearized data, non-linearized data (b) without, and (c) with matched-filter approach. We can clearly see that the matched-filter approach can remove the strong density effects at the ocean bottom and correct the amplitudes and phase in both shallow and deeper parts of the inverted velocity perturbations. ## Conclusions \vspace{-1em} In this abstract, we propose a matched-filter based least-squares reverse time migration formulation to remove the strong density variation related components in the observed data. In contrast to other methods which invert density as one additional output or reform image condition, our method doesn't require the work related to finite difference. Our modified formulation invert for the matched-filter and velocity perturbation simultaneously that matching the nonlinear observed data to the linear data with respect to velocity-only. In the experiment part, we used a complex discontinuous layered model with strong density variations at ocean bottom to test our method. We showed that the proposed matched-filter approach can remove the artifacts from density perturbation and get artifacts free migrated image related to the velocity perturbation. Future work is to incorporate the surface related multiples into the formulation since strong density variation can leads to strong surface related multiples, which can further enhance the resolution of the inverted images. ##Acknowledgements \vspace{-1em} This work was in part sponsored by the SINBAD project. ```math_def \def\argmin{\mathop{\rm arg\,min}} \def\vec{\mbox{$\mathrm{vec}$}} \def\ivec{\mbox{$\mathrm{vec}^{-1}$}} %\newcommand{\Id}{\mbox{$\tensor{I}\,$}} \def\pector#1{\mathrm{\mathbf{#1}}} \def\cector#1{#1} \def\censor#1{#1} \def\vector#1{\mathbf{#1}} \def\fvector#1{{\widehat{\vector{#1}}}} \def\evector#1{{\widetilde{\vector{#1}}}} \def\pvector#1{{\breve{\vector{#1}}}} \def\pector#1{\mathrm{#1}} \def\ctensor#1{\mathbf{\mathcal{#1}}} %\def\tensorm#1{{#1}} \def\tensorm#1{\bm{#1}} %\def\tensor#1{\bm{#1}}7 \def\tensor#1{\vector{#1}} %\def\cector#1{{\vector{\underline{#1}}}} \def\hensor#1{\tensor{#1}} \def\uensor#1{\underline{\bm{#1}}} %\def\hensor#1{{\boldsymbol{\mathsf{#1}}}} \def\hector#1{\vector{#1}} %\def\hector#1{{\boldsymbol{\mathsf{#1}}}} % \def\hfector#1{\hat{\boldsymbol{\mathsf{#1}}}} % \def\ctensor#1{{\tensor{\underline{\underline{\mathsf{#1}}}}}} \DeclareMathOperator\minim{\hbox{minimize}} \def\ftensor#1{{\widehat{\tensor{#1}}}} \def\calsor#1{{\boldsymbol{\mathcal{#1}}}} \def\optensor#1{{\boldsymbol{\mathcal{#1}}}} \def\hvector#1{\hat{\boldsymbol{\mathbf{#1}}}} \def\uvector#1{{\underline{\vector{#1}}}} \def\utensor#1{{\underline{\tensor{#1}}}} ```