## 3 Numerical Simulation
### Simple Three-Layer Model  
Figure 4 shows the simple three-layer model used to test the performance of CR-MRTM. The model consisted of 2000 grid points in the horizontal direction and 600 grid points in the vertical direction, with both horizontal and vertical grid sizes of 5 m. In addition, a PML with a grid size of 60 was added around the model to eliminate boundary reflections during seismic wave propagation. The region outside the yellow dashed box in Figure 4 represents the PML. The setup of the VSP observation system is as follows. The system included 960 sources with a fixed distance of 10 m between adjacent ones. And all sources were positioned at a depth of 200 m. The positions of the left- and rightmost source points were 200 m and 9800 m, respectively. The vertical observation well was located at a horizontal position of 1500 m, with 300 receiver points spaced 5 m apart. The upper- and lowermost receiver points were located at 500 m and 2000 m, respectively. The positions of the source and receiver points are marked with red and blue solid lines, respectively, in Figure 4. In numerical experiments, a Ricker wavelet with a central frequency of 20 Hz was used as the excitation source, with a maximum recording time of 5 s and a time sampling interval of 0.5 ms. These parameters were selected to ensure that the recorded seismic data accurately reflect the subsurface medium while minimizing computational resource consumption. Furthermore, the parameters must avoid numerical dispersion and ensure stability when simulating seismic wave propagation using the finite difference method. Based on the numerical experiments, we conclude that the chosen parameters meet the required criteria. For typical scenarios and the model selected, the grid size, time sampling interval, and central frequency of the seismic wavelet must satisfy the following equations [32]:  
\[\Delta t\max(v)\sqrt{\frac{1}{\Delta x^{2}}+\frac{1}{\Delta z^{2}}}\leq C \tag{12}\]  
\[\Delta t=\frac{\min(v)}{\max(f)\cdot\Delta x} \tag{13}\]  
Here, \(\Delta t\) represents the time sampling interval, \(\Delta x\) is the grid size, f is the seismic wave frequency, v is the wave velocity, and the value of C can be obtained from the tables. In the numerical experiments presented in this study, C is 0.7774. Figure 5a shows the seismic data corresponding to the source located at \(x=500\) m. To highlight the experimental results, original seismic data were separated in order to ensure that only the compressional (P-wave) component was included as input. The portion of the seismic data within the red dashed box is enlarged in Figure 5b. In addition to the direct wave and the first reflection wave, significant multiple reflections can be observed in the data. As the model is surrounded by PML layers, these multiples only include interbed multiples, excluding surface multiples. In conventional RTM imaging methods, these multiples manifest as artifacts in the imaging results, significantly affecting the quality of the images. Figure 6 shows the image obtained using the conventional RTM method. Given the setup of the velocity model and the position of the observation system in our numerical experiments (the velocity model contains a large-scale anticline, and the observation well covers a limited depth range), it can be inferred that most of the subsurface structural information cannot reach the receiver positions through primary reflections alone. As indicated by the dashed box in Figure 6, the imaging results from the conventional RTM method show significant interference from artifacts for our velocity model. Additionally, the magnified section within the blue solid line box, displayed at the top of Figure 6, highlights the limited illumination area provided via the conventional RTM method. Figure 4: P-wave velocity model of three-layer model; the red and blue solid lines represent source and receiver points, and the dashed yellow box represents the boundary of PML.