## 2 Theory and Methodology
### Common-Receiver-Based Multiple-Wave Reverse Time Migration  
From the illumination analysis shown in Figure 1, it can be observed that using only primary reflections as effective seismic data results in a limited imaging range because of the unique characteristics of the VSP observation system. Therefore, it is necessary to utilize multiples in the VSP seismic data to enhance the illumination range of VSP imaging. In conventional RTM methods, imaging multiples directly often introduces numerous artifacts into the imaging results, thereby degrading the imaging quality. However, CRRTM can effectively suppress artifacts in imaging results. Based on CRRTM, we propose a common-receiver-based multiple-wave reverse time migration (CR-MRTM) imaging method. By exciting both primary reflections and multiples at the source location and simulating the forward propagation of the seismic wavelet at the original receiver location, this method achieves joint imaging of primary reflections and multiples. We demonstrate the effectiveness of CR-MRTM under the cross-correlation imaging condition. The cross-correlation imaging condition was first proposed by [PERSON] in 1971 [30]. The imaging results were obtained through the correlation of the source and receiver wavefields at each time step, followed by a summation of the results. Equation (1) has already presented the source-normalized cross-correlation imaging condition. For the convenience of the subsequent theoretical derivation, the cross-correlation imaging condition without source normalization can be expressed as shown in Equation (7):  
\[I_{(x,z)}=\sum_{t=0}^{T_{\max}}S(x,z,t)\cdot R(x,z,T_{\max}-t) \tag{7}\]  
The receiver wavefield \(R(x,z,T_{\max}-t)\) was obtained by back-propagating seismic data through the subsurface medium. As the original seismic data include both primary and multiple reflections, Equation (7) can be expressed as follows:  
\[I_{(x,z)}=\sum_{t=0}^{T_{\max}}S(x,z,t)\cdot[M_{0}(x,z,T_{\max}-t)+M_{1}(x,z,T _{\max}-t)+M_{2}(x,z,T_{\max}-t)+\cdots] \tag{8}\]  
Here, \(I_{(x,z)}\) represents the imaging value at the point \((x,z)\), \(S(x,z,t)\) denotes the source wavefield generated via the seismic wavelet initiated at the source location, propagating forward through the subsurface medium, \(M_{0}(x,z,T_{\max}-t)\) represents the wavefieldformed through the back-propagation of the primary reflection from the seismic data, and \(M_{1}(x,z,T_{\max}-t)\) and \(M_{2}(x,z,T_{\max}-t)\) represent the wavefields formed via the back-propagation of the first- and second-order multiples, respectively. Furthermore, because seismic data undergo transmission and reflection during their propagation through the subsurface medium, Equation (8) can be further simplified as follows:  
\[I_{(x,z)} =\sum_{t=0}^{T_{\max}}\left[S(x,z,t)\cdot(M_{00}(x,z,T_{\max}-t)+M _{11}(x,z,T_{\max}-t)+M_{22}(x,z,T_{\max}-t)+\cdots)\right.\] \[+S(x,z,t)\cdot(M_{01}(x,z,T_{\max}-t)+M_{02}(x,z,T_{\max}-t)+M_{03 }(x,z,T_{\max}-t)+\cdots)\] \[+S(x,z,t)\cdot(M_{10}(x,z,T_{\max}-t)+M_{12}(x,z,T_{\max}-t)+M_{13 }(x,z,T_{\max}-t)+\cdots)\] \[\left.+\cdots\right] \tag{9}\]  
Here, \(T_{\max}\) represents the maximum recording time of the seismic data, while \(M_{00}\) denotes the wavefield obtained through the reverse propagation of the primary wave along its original path. Similarly, \(M_{11}\) and \(M_{22}\) correspond to the wavefields generated via the reverse propagation of the first-order and second-order multiple waves along their respective propagation paths. In contrast, \(M_{01}\) and \(M_{02}\) represent the wavefields resulting from the reflection and transmission of the primary waves along the propagation paths of the first-order and second-order multiples within the subsurface medium, which are erroneously back-propagated. The first term on the right side of Equation (9) includes the result of the cross-correlation between the seismic data propagated along its original path to the reflection point and the source wavefield.