VSP Imaging Using Free-Surface Multiples With Wavefield Decomposition: Synthetic and Field Data Examples

Vertical seismic profiling (VSP) is an effective technique to provide high-resolution seismic images of the reservoir area. However, the quality of the images is limited by the poor illumination of primary reflection wave. In conventional VSP imaging, only the upgoing primaries are used. Adding free-surface–related multiples into the imaging process can significantly improve the coverage of the illuminated area. Conventional migration methods using multiples need the complex process of multiple prediction. Data-to-data migration (DDM) is an effective imaging technique for multiples in which the recorded data is migrated directly. To improve the imaging quality of DDM in VSP imaging, we propose separating the wavefield into downgoing and upgoing components using Hilbert transform when reverse-time migration (RTM) is implemented in DDM, and the inverse-scattering imaging condition is further applied to the decomposed wavefields. The proposed method eliminates low-frequency noises and false images generated from the conventional cross-correlation imaging condition, and further enhance the illumination in the VSP imaging. Synthetic examples and application to a walkaway field data demonstrate that it can attenuate the noise and improve the imaging resolution effectively. By using DDM with inverse scattering imaging condition and wavefield decomposition based on Hilbert transform, VSP imaging using free-surface–related multiples becomes a practical complement for conventional VSP imaging.


INTRODUCTION
Vertical seismic profiling (VSP) surveys differ from surface seismic surveys or crosswell surveys in that the surface sources and the borehole receivers are used to record both upgoing and downgoing wavefields (Stewart et al., 1984;Hardage, 1985;Chang and McMechan, 1986;Hinds et al., 1996). The receiver well is placed near the target area to obtain sufficient reflection waves generated from the reservoir (Burch et al., 2010). The configuration of VSP gives the benefits to understand corresponding geologic logs and provide additional seismic interpretation insights. However, the results of VSP imaging are restricted by the illuminated area of the primary reflections (O'Brien, et al., 2013). To greatly extend the subsurface illumination, free-surface related multiples recorded in the VSP surveys are also used in the migration of VSP data (Yu and Schuster, 2001). Jiang et al. (2005) use the mirror imaging condition to migration the first-order multiples in the VSP data. But the method needs to calculate the traveltimes of the raypath. To avoid the picking of traveltimes, seismic interferometry theory (Wapenaar and Fokkema, 2006;Schuster, 2009) is employed to transform different orders of freesurface related multiples into virtual primaries and then applied in the migration process (Yu and Schuster, 2002;Jiang et al., 2007). He et al. (2007) demonstrated the wave-equation interferometric migration generates an image volume with wide coverage for 3D VSP data. Soni and Verschuur (2015) used full-wavefield migration to enhance the illumination for VSP measurements. Recently Marchenko imaging also emerges as an alternative tool to analyze the response of multiple reflections (Singh et al., 2015;Wapenaar et al., 2017;Lomas et al., 2018;Zhang and Slob, 2019) and include the contribution for VSP imaging.
Several methods have been developed to image free-surface related multiples directly in surface seismic surveys. Instead of taking as coherent noise, the multiples are used as areal sources in the migration process (Guitton, 2002;Shan 2003;Verschuur and Berkhout, 2005;Artman, 2006;Muijs et al., 2007). Recently, migration of multiples has shown the significant advantages for enhancing areal illumination to image subsurface structures (Lu et al., 2015;Li et al., 2017;Liu et al., 2018;Zhang et al., 2020). However, most methods involve the process to separate the surface-related multiples from the original data, which is complex and prone to error for real data applications. To avoid the separation of the primaries and free-surface related multiples, Wang et al. (2014a) propose an approach that can simultaneously migrate the primaries and freesurface related multiples. The recorded data containing primaries and free-surface related multiples are backward-propagated as the receiver wavefield, and the recorded data, together with a synthetic wavelet, are forward-propagated as the source wavefield. Wang et al. (2014b) isolate the contribution of multiples and name it as data-to-data migration (DDM), in which the recorded data containing primaries and free-surface related multiples are forward and backward propagated simultaneously. The algorithm is designed for surface seismic profile and can also be applied in the VSP data processing. Using the source-receiver reciprocity, the common receiver gathers is similar to the common shot gathers in the surface seismic surveys except that the virtual source is located in the borehole.
Kirchhoff migration (Keho and Beydoun, 1988;Gray and May 1994;Bevc, 1997;Hua and McMechan, 2003) and wave-equation migration (Gazdag, 1978;Stoffa et al., 1990;Ristow and Rühl, 1994;Sava and Vasconcelos, 2011) are the most common migration algorithms in the migration of free-surface related multiples for VSP (Jiang et al., 2005;He et al., 2007). For the surface seismic surveys, reverse-time migration (RTM, Baysal et al., 1983;McMechan, 1983) has shown its superiority in handling complex geologic structures and potential for true-amplitude, high-resolution migration. The main problem for RTM is that low-frequency noise and false images are generated when the source and receiver wavefields are cross-correlated near the strong velocity gradients (Liu et al., 2011;Fei et al., 2015). This is due to the two-way waveequation is used in the wavefield simulation. To eliminate the high-amplitude, low-frequency noise along the wave paths, Fletcher et al. (2005) propose to apply the directional damping factor. Yoon and Marfurt (2006) use the Poynting vector to calculate the direction of the wavefields and cross-corelate the desired component. The most practical method is the Laplace filter proposed by Zhang and Sun (2009), which has been widely applied in the applications of RTM. Liu et al. (2011) use Hilbert transform to separate the upgoing and downgoing components and avoid the storage of the entire wavefields. Fei et al. (2015) point out that only the crosscorrelation of the downgoing source wavefield and upgoing receiver wavefield are the correct imaging condition when strong velocity contrasts exist in the velocity model. Wang et al. (2017) show the wavefield decomposition method based on Hilbert transform can eliminate the noises and false images in the DDM for surface seismic survey. Zheng et al. (2018) use the similar approach to separate the upgoing and downgoing components in the 3D forward modeling and 3D RTM. The method can be also applied in the migration of VSP data. With the decomposed wavefields, it is possible to apply some advanced imaging condition to obtain better estimation of subsurface reflectivity, such as the deconvolution imaging condition (Valenciano and Biondi, 2003;Guitton et al., 2007;Lu et al., 2015) and inverse-scattering imaging condition (Whitmore and Crawley, 2012;Suh and Wang, 2013). The inverse-scattering imaging condition is derived from the inverse theory and high-frequency approximation, which can generate subsurface images with preserved amplitudes and high resolution (Pestana et al., 2014;Duprat et al., 2015). The wavefield separation used in the imaging can be achieved with high accuracy with Hilbert transform instead of Poynting vectors (Yoon and Marfurt, 2006).
In the following sections, we first introduce the theory of DDM for VSP data, and then we illustrate how to use inverse-scattering imaging condition and wavefield decomposition with Hilbert FIGURE 1 | Schematic diagram of DDM method for VSP data. In conventional migration of walkaway VSP data, the upgoing component recorded at R is backward propagated and cross-correlated with the forwardpropagated source wavelet excited at S 2 to image X 1 ; In migration of multiples. In migration of multiples, the receivers are considered to be virtual sources in the borehole and the sources located at the surface are taken as receivers for imaging multiples. The fist-order multiples are backward extrapolated from S 2 and cross-correlated with the forward-propagated direct wave from S 1 to image X 2 ; Similarly, the second-order multiples are backward extrapolated and cross-correlated with the forward-propagated first-order multiple to image X 3 .
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 730184 transform to improve the results. Synthetic examples are used to validate the effectiveness of our approach. Then the method is applied to walkaway field data which is collected to monitor the injection process of CO 2 . The final part is the conclusion of our work.

METHODOLOGY DDM for VSP Data
The DDM method for surface seismic data has been demonstrated by Wang et al. (2014b). Using the source-receiver reciprocity theory, if the VSP data are resorted into common receiver gathers and the receiver in the borehole is taken as a source. They are similar to the common shot gathers of surface seismic survey except that the source in reciprocal domain is located in the well. Figure 1 illustrates the wavepath in the VSP surveys with free surface. In the gathers recorded at R, the data contains the direct wave excited at S 1 , the primary excited at S 2 , first-order free-surface related multiple excited at S 2 and second-order free-surface related multiple excited at S 3 . Only the primary is the upgoing components and can be separated from the original data by f-k filtering. In the conventional migration of VSP data, the primary is backward propagated and crosscorrelated with the forward-propagated source wavelet excited at R to image X 1 . Here the source-receiver reciprocity is used. In migration of free-surface related multiples, the first-order freesurface related multiple excited at S 2 is backward propagated and cross-correalted with the direct wave excited at S 1 to image the reflector X 2 . And the second-order free-surface related multiple excited at S 3 is backward propagated and cross-correalted with the first-order free-surface related multiple excited at S 2 to image the reflector X 3 . From the comparison of the wavepaths in the diagram illustration, the imaging results of free-surface related multiples has wider coverage than the conventional migration. Moreover, it can clarify the shallow reflectors, which are usually not imaged when using primaries only. The imaging condition of DDM is (Wang et al., 2014b) For VSP data, D represents the downgoing components in the common receiver gathers. The subscript F denotes forward propagated, and B denotes backward propagated. The same data are forward and backward propagated and crosscorrelated to form the subsurface image. The results of DDM contain some artifacts related with the undesired crosscorrelations, such as the cross-correlation of direct waves from different shots or the cross-correlation of the direct wave and second-order free-surface related multiple. The first type mainly exists at the surface and can be muted easily. The second has longer wavepath and much weaker energy.

DDM With Inverse-Scattering Imaging Condition
There are different algorithms to implement DDM, such as oneway wave-equation  or two-way waveequation (Wang et al., 2014b). Now RTM based on two-way wave-equation has shown its benefits in offering high-resolution images and handling complex subsurface structures. The original imaging condition in RTM is where x represents the space location and t represents the time. In the DDM for VSP data, S(x, t) and R(x, t) are the forward and backward propagated wavefields of the common receiver gathers from the shot locations.
As the two-way wave-equation is used in the wavefield simulation of S(x, t) and R(x, t), the upgoing and downgoing components both exist in the wavefields. The subscripts d and u are used to represent the downgoing and upgoing components respectively. Then the image in RTM can be divided into four parts The last two terms, , generate the high-amplitude, low-frequency noise in the RTM results. Liu et al. (2011) use Hilbert transform to eliminate such kind of noise. Fei et al. (2015) show that the second term, S u (x, t)R d (x, t), generates false images near the velocity interface or strong velocity contrasts. Thus the de-primary imaging condition is proposed (Fei et al., 2015;Wang et al., 2016;Wang et al., 2017): It is the cross-correlation imaging with wavefield decomposition. True-amplitude imaging is an attractive topic for RTM. Based on the high-frequency asymptotic and the imaging/inversion theory, it is possible to obtain the estimation of slowness perturbations in wave-equation migration (Kiyashchenko et al., 2007). In RTM, the inversescattering imaging condition is proposed and shown the benefit of better amplitude recovery and higher resolution (Whitmore and Crawley, 2012;Pestana et al., 2014;Duprat et al., 2015). The imaging formula can be expressed as: where v(x) represents the velocity. z t and ∇ are the time derivative and spatial gradient operator, respectively. Compared to the conventional cross-correlation imaging condition, Eq. 5 can preserve the amplitudes and improve the resolution. In the imaging process, it is necessary to separate the wavefields to attenuate the backscattered noise introduced by two-way waveequation. Pestana et al. (2014) uses the Poynting vector to obtain the separated components. However, their results show that the method may damage some reflections and has poor performance on receiver field. In this work, we followed the approach based on Hilbert transform to decompose the wavefield effectively.

Wavefield Decomposition Based on Hilbert Transform
According to the definition of Hilbert transform, the Hilbert transform of a signal f(t) has the following quality: where H t and F t represent the Hilbert transform and Fourier transform along the time axis, respectively. i is the imaginary unit, ω is the frequency, and sgn(ω) is the sign function. With Hilbert transform, the downgoing and upgoing wavefields can be computed by (Zheng et al., 2018): where H z and H t represent the Hilbert transform in depth and time, respectively. As the wavefield modeling and Hilbert transform are both linear operators, the H t (S(x, t)) and H t (R(x, t)) in Eq. 7 is calculated by the forward modeling of the Hilbert transformed wavelet and data to avoid the wavefield storage. At each time step, Hilbert transform in depth is applied to the two wavefields, H t (S(x, t)) and H t (R(x, t)) to get the final image.
The workflow for VSP imaging using free-surface related multiples with inverse scattering imaging condition and wavefield decomposition consists of the following steps: 1) separate the upgoing and downgoing components by f-k filter and resort the shotgathers to common receiver gathers; 2) forward propagate the downgoing data from the source locations and store the boundary values; 3) apply Hilbert transform to the downgoing data and forward propagate it to construct H t (S(x, t)); 4) use the boundary value to reconstruct the wavefield S(x, t) and H t (S(x, t)); 5) back propagate the downgoing data from the source locations to construct R(x, t); 6) apply Hilbert transform to the downgoing data and back propagate it to construct H t (R(x, t)); 7) apply imaging condition in Eq. 5 and stack all the images at all time steps.

SYNTHETIC EXAMPLES
To better illustrate the advantages of VSP imaging using the proposed method, we applied it on a part of the Sigsbee 2A  The generated data without free-surface-related multiples, which contains downgoing direct waves and upgoing primaries. (B) The gathers generated with free surface, which contains direct waves, primaries and free-surface related multiples.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 730184 model. Figure 2 shows the velocity model. 500 shots are deployed on the surface and 80 geophones are placed evenly in the observation well located at the center of the model. The time sampling interval is 1 ms and the grid spacing of the model is 10 m. A Ricker wavelet with peak frequency of 15 Hz is used as the source wavelet. Figure 3 shows the synthetic data generated without and with free-surface related multiples. Figure 4 shows the data for VSP imaging after f-k filtering. The upgoing components are used for conventional imaging with primaries only while the downgoing components are used for VSP imaging using free-surface related multiples. Figure 5 shows the input common receiver gathers for VSP imaging. Figures 6A,B shows the conventional cross-correlation imaging (Eq. 2) results using primaries and free-surface related multiples, respectively. The comparison clearly demonstrates the benefits of VSP imaging using multiples. Figure 6B has much wider imaging area and the reflectors above the geophones are also imaged. The results using cross-correlation imaging condition with wavefield decomposition (Eq. 4) are shown in Figure 7. Figure 7A is the result using primaries only and Figure 7B is the result using multiples. The noised are suppressed but the amplitudes decrease with the depth. We apply the inverse scattering imaging condition with wavefield decomposition (Eq. 5) to the data and the results are shown in Figure 8. Figures 8A,B are the results using primaries only and multiples, respectively. Compared to Figure 7, the proposed  Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 730184 approach yields a better estimation of the true subsurface reflectivity with improved resolution and amplitude recovery. And there are no obvious artifacts shown in Figure 8B, as the true images has stronger contribution to the final image.

FIELD DATA APPLICATION
To verify the adaptability of the proposed method to field data, we use a walkaway field data for further test. The walkaway VSP data is collected to monitoring the reservoir changes during the CO2 injection in northwestern China. Figure 9A shows the migration velocity and seismic geometry in the walkaway VSP survey. The data are recorded by 40 receivers equally spaced from the depth of 390-1,170 m. The time sampling interval is 1 ms and the grid spacing of the model is 10 m. Figure 9B shows a common receiver gather for imaging using multiples. Several shots near the wellbore are missing, which can lead to the lack of multiples with small reflection angles in the generated data. Figures 10A,B shows the images obtained from conventional cross-correlation imaging condition with wavefield decomposition using primaries only and free surface related multiples, respectively. Figures 11A,B shows the results of the proposed method. The inverse-scattering imaging condition with wavefield decomposition improves the resolution and balances the amplitudes.
As we can see from the results, migration with multiples can effectively image the reflector in the shallow zone and enlarge the image range. But the illumination near the borehole is influenced by the reflection angle and missing near-offset traces in the data. Note that it also contains the crosstalk generated from undesired cross-correlation of different seismic events. Figure 12 illustrates the migration result that combines the contribution of primaries and multiples. Compared to the conventional images, the joint migration enhances the illumination and achieves the highquality seismic images for monitoring the CO 2 injection.   Overall, we see that VSP imaging using inverse-scattering imaging condition with wavefield decomposition is capable of handling real datasets.

DISCUSSION
In the proposed workflow, the downgoing components of VSP data are included in the imaging process. Combined with the conventional VSP imaging using upgoing waves, we can have improved imaging results. An important advantage of DDM is that it can image the subsurface without any knowledge of the source information. Thus, the reflector in DDM images can be used to calibrate the conventional images when the source wavelet is inaccurate. The matching filter used in the multiple subtraction (Verschuur et al., 1992;Wang, 2003;Fomel, 2009) can be modified to find the correct combination of the two kinds of imaging results. In the application of DDM to VSP walkway data, the imaging quality is contaminated by the large-amplitude, low frequency noises around the well. By using wavefield decomposition based on Hilbert transform, this kind of noises can be suppressed effectively.
In the wavefield separation using Poynting vectors, the reflection angle between the source wavefield and receiver wavefield is calculated, and the components with the opening angle more than certain values are excluded to obtain the final imaging results. Compared to the approach using Hilbert transform, the method using Poynting vectors can not remove the false images, which are generated by the second term on the right side of Eq. 3. But it is able to generate angle domain common image gathers. The proposed approach is only FIGURE 9 | (A) The migration velocity and walkaway VSP acquisition used in field application. (B) The common receiver gathers obtained from the data containing multiples. Note that several near-offset traces are missing.
FIGURE 10 | Imaging results of the walkaway VSP field data using cross-correlation imaging condition with wavefield decomposition: (A) The image using primaries only. (B) The image using DDM. Compared to (A, B) has enhanced illumination away from the well and images the shallow reflectors. Due to the lack of near-offset traces, the zone near the well is not accurately imaged in (B).
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 730184 designed to separate upgoing and downgoing components in the propagated wavefields, and to provide imaging results without low-frequency noise and false images. If the finite-difference time-domain method and staggered grid are used in RTM, the extra computation cost is negligible for the wavefield separation using Poynting vectors. The computation cost in the proposed approach is doubled compared to traditional RTM, thus some high-performance computation techniques, such as MPI or GPU, are necessary to increase the efficiency of the method. The images generated from conventional cross-correlation imaging condition may suffer from poor resolution and incorrect amplitude responses while inverse-scattering imaging condition partly overcomes the drawbacks. To make the results obtained from inverse-scattering imaging condition become suitable estimation of the slowness perturbation, wavefield decomposition is necessary to avoid the noises introduced by the two-way wave equation when strong velocity gradients exist. It is an approximated solution for true amplitude RTM. The final workflow can generate subsurface images with balanced amplitude and high resolution, which is important for amplitude variation with offset analysis.
As we have demonstrated with numerical examples, the final images obtained from DDM can be superior to conventional imaging results. However, it relies on the designed geometry. The sparse shot array and shallow receiver locations may lead to limited free-surface related multiples in the recorded data. In FIGURE 11 | Imaging results of the walkaway VSP field data using inverse-scattering imaging condition with wavefield decomposition: (A) The image using primaries only. (B) The image using DDM. Compared to Panel 10, the reflection amplitudes are well preserved, and the resolution is improved. Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 730184 8 practice, the results of DDM can be used as import complement of conventional imaging, as shown in the field example.
In the theory of DDM, because we have to cross-correlate all orders of multiples to avoid multiple prediction, the results of DDM contains undesirable artifacts. Assume that the amplitude of the true image is one, then the strongest artifact has the strength of the reflectivity. The final imaging result of DDM is still acceptable in most cases. The noise that leaks into the images can be further suppressed in image domain to improve the imaging quality. Leastsquares migration (Nemeth et al., 1999;Dai et al., 2011) can be applied to remove the artifacts (Zhang and Schuster, 2014;Liu et al., 2018;Zheng et al., 2019;Li et al., 2020). Another option is highresolution parabolic Radon filtering in angle domain common image gathers (Wang et al., 2014b;Zheng et al., 2016) or using 3D wideazimuth acquisition (VerWest and Lin, 2007).

CONCLUSION
We present an effective method to imaging free-surface related multiples in VSP data. The downgoing components, which are muted in conventional VSP imaging, are used to image the subsurface. The separated downgoing components are resorted into common receiver gathers and then DDM method is used to image the subsurface. The inverse-scattering imaging condition with wavefield decomposition is applied to generate a better approximation of the subsurface reflectivity. We have illustrated the explicit workflow to apply the proposed approach with Hilbert transform in RTM. Results with synthetic VSP data and field data validate the effectiveness of the proposed approach. Compared to VSP imaging using primaries only, the results have significantly enhanced coverage and shallow reflectors with improved resolution. The algorithm handles free-surface related multiples without the need of multiple prediction in preprocessing. The method has the potential to improve the results using primaries only and offer better estimation of the geological structures.

DATA AVAILABILITY STATEMENT
The field data analyzed in this study was obtained from National Institute of Clean-and-Low-Carbon Energy and used for research purposes only. Requests to access these datasets should be directed to NICE (nice_contact@chnenergy.com.cn).

AUTHOR CONTRIBUTIONS
YZ and YW developed the theory and performed the applications. XC supervised the findings of this work. All authors discussed the results and contributed to the final manuscript.