Inverse Gaussian beam stacking method for imaging both primary reflections and free-surface multiples in 2D VSP

In vertical seismic profiling (VSP) exploration, it is difficult to produce an accurate image for large-offset reflections only using reflection waves and the image resolution is low in traditional VSP-CDP stacking as is the number of folds of reflection points. To mitigate these problems, we present an inverse Gaussian beam stacking method for imaging both primary reflections and free-surface multiples. We first compute the stacking weighted functions at each trace location by Gaussian beam forward modeling, and then apply an inverse projection for VSP data to produce common shot gathers (CRP). Since inverse Gaussian beam stacking maps the common-shot data along finite-frequency wave-paths instead of single rays as the traditional ray-based stacking method does, it enlarges the reflection-point coverage, increases stacking fold and reduces the requirement for large bin sizes. We incorporate free-surface multiples into the proposed inverse Gaussian beam stacking, which enables us to expand the horizontal imaging aperture and mitigate the low-fold problem of primary reflections in the shallow large-offset regions for VSP surveys. Numerical examples for synthetic and field data demonstrate the feasibility and adaptability of the proposed inverse Gaussian beam stacking method for VSP data.


Introduction
As an important acquisition method in exploration seismology, VSP has been widely used in fine-scale oil/gas reservoir characterisation (Mateeva et al. 2014;Abedi et al. 2019), geothermal resource detection (Reiser et al. 2017), shale fracturing study and CO 2 injection and monitoring. In comparison to surface seismic imaging methods, VSP techniques can produce images with higher spatial resolution, and provide accurate local velocity models and time-to-depth conversion relation for well ties (Sandoval et al. 2019;Shi et al. 2019). Ray-based stacking is the most commonly used method for VSP data. Based on the assumption of 1D velocity model, the concept of VSP-CDP stacking was introduced (Wyatt 1981). This approach does not consider lateral velocity variations and therefore has limited application for complex structures. The common reflection surface stacking method was produced and had higher signal-to-noise ratios than the conventional stacking methods in complex structures. An improved 852 common reflection surface stacking method was presented to solve the problem of same phase-axis intersection on a zerooffset section (Schwarz et al. 2014). The Gaussian beam method was introduced into CDP stacking for VSP data (Yang et al. 2015), and the stacking fold and image quality were improved by decomposing single-trace data into different local plane waves. The VSP-CDP stack technique based on bin iteration was proposed and used in fold calculation (Wei et al. 2017). This method is based on the lateral velocity field instead of the single well velocity model, which further improve the accuracy of VSP imaging. The inverse Gaussian beam stacking method has been used in gas cloud exploration . This method decomposes one trace data in a common-shot point (CSP) gather into multiple points in a common-reflection-point gather, which makes the folds of reflection points more uniform. In traditional surface and VSP imaging, only primary reflections are considered to be effective signals, and the free-surface and interbed multiples are commonly suppressed in pre-processing (Rao & Wang 2016;Ogagarue & Nwankwo 2019). However, the multiples hold valuable information on subsurface structures and lithology. Multiples are mainly divided into free-surface multiples and internal multiples. Commonly, free-surface multiples occur in marine seismic data and internal multiples occur in land data. In VSP data and non-zero-offset VSP data with a large offset, free-surface multiples are well developed. VSP free-surface multiples are downward reflected waves, which can be identified based on the relationship of the slope of a seismic event between free-surface multiples and a direct wave. There are lots of methods in seismic data wave field separation, such as the median filtering method, singular value decomposition method and high precision -p domain seismic wave field separation method. The first-order multiple operator was first introduced into Kirchhoff depth migration for ocean bottom hydrophone data, which migrated both primary and deep-water multiples (Reiter et al. 1991). The multiples were transformed into primaries based on the free-surface multiple elimination (SRME) theory and joint migration for primaries, and converted multiples were realised using the traditional imaging operator (Berkhout & Verschuur 2006). A controlled-order LSRTM was developed for migrating surface multiples through selecting different reflection orders (Liu et al. 2018). This approach reduces the crosstalk artifacts and improves image quality. In the face of VSP exploration with a large offset in a deviated well, the area above the well trajectory cannot be imaged when using primary reflections. However, the free-surface multiples can be used for fine imaging in this area, and make up the shortcomings of insufficient illumination from VSP primary reflections in the shallow area of a large offset.
Although VSP migration for surface multiples has been studied by more and more people, in the actual production of VSP, the stacking image method is still the main imaging tool at present, and the VSP stacking image methods for free-surface multiples are currently found in few studies. In this study, we present an inverse Gaussian beam stacking method for imaging both primaries and free-surface multiples in VSP. First, we compute the weighted functions for each beam at receiver locations. Then, we transform the common-shot-gather data including upgoing wave and firstorder free-surface multiples into a common-reflection-point (CRP) domain within a certain range. Finally, we select an appropriate bin size for stacking.

VSP Gaussian beam forward simulation
In Gaussian beam forward modeling, the seismogram at a given receiver is computed by summing the contributions of all beams in their effective vicinity (Červený et al. 1982;Huang et al. 2016). The 2D scalar Gaussian beam has the expression, , where (s, n) is the ray-centered coordinate, v(s) is the velocity, p(s) and q(s) are the dynamic ray tracing parameters, (s) is the travel time of the central rays, n is the distance from the adjacent ray to the center ray, K(s) is the wave front curvature of seismic wave, L s is the effective half width of the seismic wave propagating to the detection point and Φ( ) is the weighted function and is related to the imaging resolution (Sun et al. 2020).
As shown in figure 1, the wave field of primary reflections transmitting from shot point to receiver point by Gaussian beam ray tracing can be expressed as,  and the wave field of free-surface multiples can be expressed as, where GP means primary reflections of a Gaussian beam, and GM denotes free-surface multiples of a Gaussian beam. Based on the Gaussian beam effective neighborhood wave field approximation theory, the wave field at the receiver point is the Gaussian weighted stacking of the ray energy in its effective half-width range, which is consistent with the Gaussian normal probability density function. Based on the physical meaning of the normal distribution function, we approximate the weighted function relationship between adjacent rays and the central ray as, As shown in figure 2 (L s is 50 m), when adjacent rays coincide with the center ray (that is, n = 0), the wave field of re- ceiver point is illustrated by the energy of central ray. When n < L s , the wave field of receiver point is illustrated by the energy of adjacent rays that are in effect half width. When n > L s , the wave field energy of adjacent rays does not contribute to the wave field energy of the receiver point.
It can be seen from figure 2 that as the distance between adjacent rays and the center ray increases, the value of the seismic wave field gradually decreases. If the weighted function is used for VSP-CDP conversion, the wave field of the effective ray calculated by Gaussian normal distribution weight function is just equal to the original. So, the seismic data is amplitude-preserved.

Inverse Gaussian beam stacking
For a horizontal reflector (as shown in figure 3), the travel time of the primary reflections and first-order free-surface multiples in VSP acquisition can be expressed as, where v is the velocity, h is the depth of the layer, X 0 is the horizontal source-well distance, Z R is the depth of the receiver point in well and T denotes the travel time from source to receiver. The symbol '∓' is used to distinguish the primary reflections (−) and free-surface multiples (+). The timedomain VSP stacking projects the common-shot data from the receiver-depth and time domain (Z-T) to the source-well distance and time domain (X P∕M , P∕M ), and then produces a CRP profile. According to the geometric relation illustrated in figure 3, the CRP gather for the primaries and free-surface multiples can be computed as, where the subscript 'P' for the coordinate (X P∕M , P∕M ) denotes the primary reflections, which corresponds to the '−' symbol in equation (7). Similarly, subscript 'M' stands for free-surface multiples and corresponds to the symbol '+' .
P∕M is the two-way traveltime from reflection point to surface.
The theory of the Gaussian beam effective neighborhood states that the wave field size of the receiver point is the weighted stacking of adjacent rays' energy in its effective halfwidth range, and the reflection points corresponding to the receiver point are the set of reflection points related to its effective half-width. Their coordinates can be illustrated as, where L s is the effective half width and is the incidence angle. The symbol '−' is the present primary reflection and the symbol '+' is the present free-surface multiples. The conventional VSP-CDP conversion converts a single sample of the common-shot gather into a single sample of the CRP gather, which can be expressed as, But, in inverse Gaussian beam VSP-CDP conversion, the weighted primary reflections and free-surface multiples at each receiver can be projected to the CRP domain, which can be expressed as, Figure 4 is a schematic diagram of the inverse Gaussian beam stacking for primary reflections and free-surface multiples in the VSP seismic wave fields. The CRP gather data generated by the inverse Gaussian beam VSP-CDP conversion are divided into bins at with proper spatial intervals, and the projected samples in the same bin are summed together to produce one stacked sample, and the final stacking image is obtained. Compared with conventional raybased VSP-CDP stacking method, the number of reflection points in the same bin is increased due to the finite-frequency property of Gaussian beams. In addition, the inverse Gaussian beam VSP-CDP conversion is implemented for both primaries and free-surface multiples, which improves the VSP stacking accuracy and enhances the VSP imaging coverage.

VSP stacking image workflow
The proposed inverse Gaussian beam stacking method for primaries and free-surface multiples in VSP is summarised in the following steps (as shown in figure 5 points, the weight functions of the effective rays for VSP primary reflections and the free-surface multiples using forward Gaussian beam simulation. 2) The primary reflections and the free-surface multiples are separated. A reverse-time Gaussian beam wavefield projection is applied for separated CSP data to produce CRP gathers, and a suitable stack bin size is chosen for final stacking.

Numerical experiments
Based on the VSP inverse Gaussian beam stacking image method, this section will test the correctness and robustness of the method through the horizontal stratigraphic model. The parameters of the geological model and the parameters of the geometry are shown in Tables 1 and 2, respectively. The geological model and reflection points (primary reflection and free-surface multiples) are shown in figure 6. We can see that the illumination range is expanded by the freesurface multiples. Figure 7 shows the wave field of primary reflection and free-surface multiples in VSP. The VSP Gaussian beam ray tracing method is used to calculate the weighted functions of each ray to the energy of the receiver point and the corresponding reflection point coordinates, and one trace 856  of the primary reflection and free-surface multiples wave record is decomposed into multiple reflection points records according to the Gaussian beam effective neighborhood wave field approximation theory. Then, a super commonreflection-point gather is formed and 10 m is chosen for stacking bin in imaging. Figure 8a is a time-domain imaging result obtained by a traditional VSP-CDP stacking method based on a two-point ray tracing method and figure 8 parts b and c are imaging results by the inverse Gaussian beam stacking method. Figure 8b is the imaging result of primary reflections and figure 8c is an imaging result produced using both primary reflections and free-surface multiples in VSP.
As shown in figure 8a and b, because the inverse Gaussian beam stacking method effectively produces more intensive reflection point folds, the resolution of the inverse Gaussian beam stacking imaging result is higher than the traditional VSP-CDP stacking method. From the red box in the figure 8c, we can see that it can effectively increase the VSP imaging range and solve the problem of large-offset imaging by using both primary reflections and the free-surface multiples.

Practical application
To test the performance of the proposed inverse Gaussian beam stacking method, we apply it to a VSP field data for a Bi434 well in NanYang oil field in China. The velocity model (as shown in figure 9) was built combined with surface seismic interpretation results, logging velocity and velocity analysis of zero-offset VSP. The Bi434 well is located in the middle of the chosen profile. There are three shots in VSP exploration. The offsets of every shot are −393.19, +525.33 and +863.03 m, which correspond to shot 1, shot 2 and shot 3, respectively. The total number of traces is 416 and trace interval is 5 m. The depth of the first trace is 435 m, and the depth of the last trace is 2510 m. The recording length is 4 s 857 Figure 8. The VSP imaging result. (a) Traditional VSP-CDP stacking using primary reflections, (b) inverse Gaussian beam stacking using primary reflections and (c) inverse Gaussian beam stacking using both primary reflections and free-surface multiples.
Figure 9. Velocity model. and the sample interval is 0.5 ms. For the common-shot VSP data, we first separate the upgoing and downgoing waves. Taking the third shot as an example, the primary reflections (as shown in figure 10a) are computed by extracting the upgoing waves and the free-surface multiples (as shown in figure 10b) are computed by choosing the downgoing waves. The inverse Gaussian beam stacking results using only primary reflections are present in figure 11a, and the corresponding results using both primary reflections and free-surface multiples are shown in figure 11b. The imaging result is expanded by free-surface multiples in the red box. Compared with the imaging results from surface seismology (as shown in figure 11c), the characteristics of seismic wave groups in the VSP stacking section are consistent with the surface seismic imaging results, verifying the feasibility of the proposed method in this study. The characteristics of spectrum for VSP and surface imaging results show that the frequency band of the VSP result is broader than the surface seismic imaging result, indicating a higher spatial resolution.

Conclusion
Based on Gaussian beam forward modeling, we present an effective inverse Gaussian beam stacking method for imaging both primary reflections and free-surface multiples in VSP. Compared to a traditional two-point ray-based stacking method, this method uses a Gaussian beam to compute the weighted functions of effective rays from source to CRP locations and to receivers, which allow us to consider wave fields reflected from a region on the subsurface rather than only one point. Therefore, each sample in the common-shot VSP data can be inversely projected to common reflection regions along Gaussian beam paths, and stacking all the projected samples produces a CRP gather with increasing folds compared to the traditional method. In addition, incorporating both primaries and free-surface multiples helps us to enlarge the CRP range, especially for shallow large-offset areas. This mitigates the insufficient illumination problem of VSP stacking only using primary reflections. However, the calculation accuracy of this method is related to the velocity model, so the accuracy of the velocity model affects the final imaging quality during actual data processing. Synthetic and field data examples show that the proposed inverse Gaussian beam stacking method can effectively image VSP reflection with good illumination for both shallow and deep target layers, and produce sections with higher spatial resolutions than surface seismology for fine-scale reservoir interpretation and characterisation.