**1. Introduction**

High-resolution seismic imaging is a critical tool to acquire information on underground structures from observed seismic data [1,2]. Reverse time migration (RTM) has no assumption for the high-frequency approximation in traditional ray methods and shows good performance in lateral velocity varying media [3,4]. However, as the real data may have limited borehole diameter for observation, irregular observation systems, and limited wavelet frequency band, and thus images obtained from RTM have limited resolution and poor amplitude fidelity, providing inaccurate estimates of underground reflectivity [5].

By using the reverse Hessian operator on the conventional migration, the least-squares migration (LSM) allows for higher resolution, fewer migration artifacts, and better fidelity [6]. Since the Hessian matrix is too large in scale and the direct inversion process is not realistic, the LSM method is always recast as a data-driven linear optimization problem. Due to the good performance on improving resolution, it has been investigated in seismic imaging for acoustic-wave single component data [7–18] and elastic-wave multiple component data [19–22]. Besides, the elastic least-squares migration is extended to obtain multiparameter images, such as P- and S-wave velocity and density, to cope with the trade-off effects [23,24]. However, since this data-domain LSM has a high dependency on the accuracy of wavelet and velocity, most achievements for LSM are mainly obtained in theoretical experiments. Furthermore, in 3D data processing in practice, this data-domain LSM always suffers from unstable convergence and excessive computational cost, which make it quite difficult in practical applications.

**Citation:** Li, B.; Sun, M.; Xiang, C.; Bai, Y. Least-Squares Reverse Time Migration in Imaging Domain Based on Global Space-Varying Deconvolution. *Appl. Sci.* **2022**, *12*, 2361. https://doi.org/10.3390/ app12052361

Academic Editor: Amerigo Capria

Received: 29 November 2021 Accepted: 19 February 2022 Published: 24 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In contrast, the imaging-domain LSM aims to find a reasonable approximation for the Hessian instead of performing the expensive iterative data-fitting process. The diagonal Hessian operator is commonly used to enhance the amplitude balance by considering the illumination for observation, providing similar results as the classical amplitude-preserving migration [25]. In horizontally layered media, the Green function for migration can be also used as a deconvolution operator, thus processing the results from migration [26,27]. An unsteady-phase filter can be estimated using the images from the front and rear iterations in data-domain LSM thus replacing the inversion in Hessian matrix to calibrate and filter the migration results [28–31]. Point spread function (PSF) from one scattering point is consistent with a row of elements in the Hessian matrix, which physically describes the impact of single-point-scattered energy on underground media, including local illumination characteristics of the observation geometry, the space variation in velocity model, and the band-limiting effect in seismic wavelet and received data [32,33]. These approximations for Hessian in imaging domain LSM can be also used as preconditioners for the data domain LSM, which, after multiple iterations, can more quickly approximate the real reflection coefficient [34–37].

However, as PSFs are of significant space-varying characters, it is rough to apply a space-invariant PSF from a single scattering point on the whole area. The novelty of our method is that it introduces a regional division strategy in the image-domain least-squares migration based on a global space-varying PSF. As a result, it can divide the global PSF into sub-blocks, use a high-dimensional space-varying deconvolution in each sub-block, and eventually perform the data reduction for all the sub-blocks. Since the continuity of velocity model and illumination is much better in local sub-blocks, the proposed method can provide high-resolution images and have computation efficiency.
