## 4 The Multi-Resolution Pre-Stack Migration Imaging Method
### The Finite-Difference Migration Imaging Method  
In a Cartesian coordinate system, the governing equation for the wavefield in the passive region is as follows [33]:  
\[\
abla^{2}\mathbf{u}-\frac{1}{v^{2}}\frac{\partial^{2}\mathbf{u}}{\partial t^{2}}=0, \tag{16}\]  
In this paper, upgoing waves are used for migration imaging, Equation (16) is rewritten as a two-dimensional wave equation:  
\[\frac{\partial^{2}\mathbf{u}}{\partial x^{2}}+\frac{\partial^{2}\mathbf{u}}{\partial z ^{2}}=\frac{1}{v^{2}}\frac{\partial^{2}\mathbf{u}}{\partial t^{2}}. \tag{17}\]  
By applying a coordinate transformation to Equation (17), let \(x^{\prime}=x,z^{\prime}=z\), and \(t^{\prime}=t+z/v\). Since the wave field remains invariant under the change of coordinates, we have \(\mathbf{u}(x,z,t)=\mathbf{u}(x^{\prime},z^{\prime},t^{\prime})\). Based on this, the following expression can be derived:  
\[\left\{ \begin{array}{c}\frac{\partial^{2}\mathbf{u}}{\partial x^{2}}=\frac{ \partial^{2}\mathbf{u}}{\partial x^{2}}\\ \frac{\partial^{2}\mathbf{u}}{\partial z^{2}}=\frac{\partial^{2}\mathbf{u}}{\partial z ^{2}}+\frac{1}{\partial z}\frac{\partial^{2}\mathbf{u}}{\partial z^{2}}+\frac{1}{ v^{2}}\frac{\partial^{2}\mathbf{u}}{\partial t^{2}}+\frac{1}{v^{2}}\frac{\partial^{2} \mathbf{u}}{\partial t^{2}}.\\ \frac{\partial^{2}\mathbf{u}}{\partial t^{2}}=\frac{\partial^{2}\mathbf{u}}{\partial t^ {2}}\end{array} \right. \tag{18}\]  
Substituting Equation (18) into Equation (16) yields the wave equation in the floating coordinate system:  
\[\frac{\partial^{2}\mathbf{u}}{\partial x^{\prime 2}}+\frac{2}{v}\frac{\partial^{2} \mathbf{u}}{\partial z^{\prime}\partial t^{\prime}}+\frac{\partial^{2}\mathbf{u}}{ \partial z^{\prime 2}}=0, \tag{19}\]  
By performing a Fourier transform on \(x^{\prime}\) and \(t^{\prime}\) in Equation (19), the first-order approximation equation can be obtained [33]:  
\[\frac{\partial^{2}\mathbf{u}}{\partial z\partial t}+\frac{v}{4}\frac{\partial^{2 }\mathbf{u}}{\partial x^{2}}=0. \tag{20}\]  
Figure 5: The schematic diagram of the preconditioned precise integration time-sweeping wavefield inverse transformation. In the above equation, the velocity \(v\) can either be a constant or a function of \(x\) and \(z\), denoted as \(v(x,z)\). Assume that Equation (20) exists in the spatial domain \(\Omega(-\infty<x<\infty,0\leq t\leq T,\)\(0\leq z<Z)\), with the following mixed boundary conditions given:  
\[\left\{ \begin{aligned} \mathbf{u}(x,z,T)&=0\\ \frac{\partial(x,z,T)}{\partial t}&=0\\ \mathbf{u}(x,z=0,t)&=\varphi(x,t)\\ \mathbf{u}(x,Z,t)&=0\end{aligned} \right.\text{.} \tag{21}\]  
where \(\mathbf{u}(x,z=0,t)=\mathbf{\varphi}(x,t)\) represents the wave field observed at the surface. By solving the above partial differential equation, the wave field values \(\mathbf{u}(x,z,t)\) of the upgoing waves at any point in the subsurface \((x,z>0)\) can be obtained. In this paper, the Crank-Nicolson symmetric implicit difference scheme is employed, as illustrated in Figure 6.  
Next, the following notation is used to derive the equation:  
\[\left\{ \begin{aligned} \mathbf{u}_{i,j}^{n}&=\mathbf{u}(i \Delta x,n\Delta z,j\Delta t)=\mathbf{u}(x,z,t)\\ \delta_{t}\mathbf{u}_{i,j+1/2}^{n}&=\left(\mathbf{u}_{i,j+1 }^{n}-\mathbf{u}_{i,j}^{n}\right)/\Delta t\\ \delta_{xx}\mathbf{u}_{i,j+1/2}^{n}&=\left(\mathbf{u}_{i+1,j+1/2}^{n}-2\mathbf{u}_{i,j+1/2}^{n}+\mathbf{u}_{i-1,j+1/2}^{n}\right)/\Delta x^{2} \end{aligned} \right., \tag{22}\]  
where \(\Delta x\) and \(\Delta z\) are the spatial steps in the \(x\) and \(z\) directions, respectively, and \(\Delta t\) is the time step. For the central point in Figure 6, the wavefield is denoted as \(\mathbf{u}_{i,j+1/2}^{n+1/2}\).