## 3 Numerical Simulation
### Marmousi Model  
We further demonstrate the effectiveness of CR-MRTM in complex scenarios using the Marmousi model, as well as the differences in the imaging quality between CR-MRTM and conventional RTM methods. The velocity model is illustrated in Figure 9. The model consisted of 2720 grids in the x direction and 560 grids in the z direction, with a grid spacing of 5 m in both directions. A 500-meter-wide PML layer was added around the model to absorb boundary reflections. The region outside the dashed yellow box in the figure represents the PML. The VSP observation system was set up as follows. The observation well is located at \(x=3000\) m, with 500 receiver points uniformly spaced every 5 m. The depth of the uppermost receiver point was 300 m. The source depth was uniformly set to 250 m, with 1335 sources uniformly distributed every 10 m. The positions of the sources and receivers are marked by the red and blue solid lines, respectively, in the figure. We used a Ricker wavelet with a central frequency of 40 Hz to generate the VSP seismic data. The maximum recording time was 6 s, and the sampling interval was 0.5 ms. Only P-waves were retained in the seismic data generated. We first applied the conventional RTM method to image the seismic data corresponding to the 500 sources on the left side. The results are shown in Figure 10. The solid red line indicates the distribution of the sources involved in the imaging. When the offset distance is small, the energy in the seismic data is primarily composed of the first reflection wave. In this case, the reflective layers within the illuminated region can be effectively imaged using the conventional RTM method. As the number of sources involved in the imaging continues to increase, the imaging results using the conventional RTM method are shown in Figure 11. Figure 11a presents  
Figure 10: Migration result using the conventional W-VSP RTM method for seismic data from 500 sources on the left side. The solid red line represents the distribution of the sources used in the imaging. Figure 9: P-wave velocity model of the Marmousi model; the red and blue solid lines represent source and receiver points, and the dashed yellow box represents the boundary of PML. the imaging result obtained by stacking the data from the 900 sources on the left side, while Figure 11b shows the result from stacking all sources. As the offset distance increases, seismic data from sources farther from the observation well exhibit a significant increase in multiples. The conventional RTM method is not effective at handling these multiples. Therefore, compared with the results in Figure 10, the imaging range in Figure 10 does not show significant improvement, and the imaging results contain more artifacts, resulting in a lower imaging quality. The results obtained using CR-MRTM are shown in Figure 12. Compared with the conventional RTM method, our approach provides a larger illumination range and higher imaging energy. In the region indicated by the blue solid line box, our method clearly highlights the steep structural features present in the original velocity model, with the magnified results displayed on the right side of Figure 12. Additionally, in the area near the observation well, our method also provides a higher imaging quality. Figure 11: Migration results using the conventional W-VSP RTM method with an increasing number of sources: (**a**) The imaging result obtained by stacking data from 900 sources on the left side; (**b**) The imaging result obtained by stacking data from all sources. Furthermore, to demonstrate that the CR-MRTM is also effective for surface multiples, we removed the PML layer from the upper part of the velocity model and regenerated the VSP seismic data for imaging. Regenerated VSP seismic data contain surface multiples, with an energy significantly greater than that of the interbed multiples.