**2. Methodology and Principle**

### *2.1. Basic Theory on LSM*

The forward problem in the classical migration imaging can be expressed as:

$$d = \mathcal{L}m,\tag{1}$$

where, *L* is the linear operator in the linear Born simulation, which expresses the relation between underground reflectivity model m and seismic scattering data *d*. The migration images can be obtained by projecting the adjoint Born operator *LT* into the seismic scattering data *d*:

$$I = L^T d.\tag{2}$$

The adjoint Borning operators *L<sup>T</sup>* can correctly reverse the propagation effects on the travel time and phase. However, due to problems such as limited observation apertures and limited data frequency band, the adjoint operator is not a good approximation for the inverse Born operators, which leads to unpreserved amplitude and low spatial resolution.

The LSM enables to provide high-fidelity images by minimizing the difference between synthetic and observed data, and the misfit function of the least-squares migration can be defined as,

$$\mathbb{C}(m) = \frac{1}{2}||Lm - d\_{obs}||\_2^2. \tag{3}$$

When the objective function reaches the minimum, the reflectivity model *m* satisfies the following equation:

$$I = \left(L^T L\right) = Hm.\tag{4}$$

In this equation, *H* is the Hessian operator, which is the second derivative of the misfit function (3). This equation is always named as the Newton normal equation, which establishes the basic principle for the LSM in the imaging domain and physically reveals that the migration image obtained by the conventional migration method is the Hessianblurring version of the subsurface reflectivity model m. Thus, the seeking reflectivity model can be realized using the generalized inverse under the sense of the norm *L*2.

$$m = \left(L^T L\right)^{-1} I = H^{-1} I. \tag{5}$$

*2.2. LSM in Imaging Domain*

2.2.1. The Global Space-Varying Point Spread Function

According to the acoustic-wave equation, Hessian can be expressed in the frequency domain as:

$$H(\mathbf{x}\_1, \mathbf{x}\_2, \omega) = L(\mathbf{x}\_1, \omega; \mathbf{x}\_{s\prime}, \mathbf{x}\_r)^T L(\mathbf{x}\_2, \omega; \mathbf{x}\_{s\prime}, \mathbf{x}\_r),\tag{6}$$

where *xs*, *xr* are the positions of shots and receivers, and *L* denotes the sensitivity kernel, satisfying

$$L(\mathbf{x}, \omega; \mathbf{x}\_s, \mathbf{x}\_r) = f(\omega)G(\mathbf{x}, \mathbf{x}\_s, \omega)G(\mathbf{x}, \mathbf{x}\_r, \omega),\tag{7}$$

where, *f*(*ω*) is a seismic wave, *G*(*x*1, *xs*, *ω*) is the Green function excited from the shot position *xs*.

The PSF is commonly used as the approximation for Hessian, including a row of elements in the Hessian matrix *x*0.

$$K(\mathbf{x}\_2, \mathbf{x}\_0, \omega) = L(\mathbf{x}\_2, \omega; \mathbf{x}\_s, \mathbf{x}\_r)^T [L(\mathbf{x}\_0, \omega; \mathbf{x}\_s, \mathbf{x}\_r) \delta(\mathbf{x} - \mathbf{x}\_0)].\tag{8}$$

Here, *K*(*x*2, *x*0, *ω*) reflects the PSF of the perturbation at the point *x*0.

Substitute Equation (8) into Equation (6), mark the point spread function at *xi* as *Ki*, and Hessian can be expressed as,

$$H = \begin{bmatrix} \ K\_1 & \cdots & \ K\_i & \cdots & \ K\_m \end{bmatrix}^T. \tag{9}$$

Substitute the above equation into Equation (4) to obtain:

$$I\_i = K\_i m.\tag{10}$$

The reflectivity model can be obtained by

$$m = \frac{I\_i}{\left(\mathbf{K}\_i + \boldsymbol{\kappa}\right)'} \tag{11}$$

where α is the regularization parameter. The Equations (10) and (11) can be also expressed in the space domain as:

$$I\_i = \mathbb{K}\_i \ast \mathfrak{m}\_\prime \tag{12}$$

and:

$$
\mathfrak{m} = I\_{\mathfrak{i}} \otimes \mathbf{K}\_{\mathfrak{i}\mathfrak{i}} \tag{13}
$$

where, ∗ and ⊗ represent the operation of convolution and deconvolution, respectively. It reveals that, the migration image at *xi* is the superposition of the convolution of the global reflectivity model and the point spread function at *xi*; on the contrary, the reflectivity model at *xi* can be computed within a high-dimensional deconvolution of the corresponding point spread functions.
