The Application of Wave-equation Datuming to 3D VSP Processing

This paper discusses a novel approach to transform the recorded upgoing P-waves to a flat datum on the surface based on the wave-equation. The transformation results in wavefields which appear as if they are recorded at the surface. The result of the transformation enables the application of CMP-based processing techniques, such as velocity analysis, Normal Moveout (NMO), and static correction to 3D VSP data. A good transformation process requires an accurate measurement of downgoing P wave arrival time at each geophone. The 3D VSP data used in this study was acquired for reservoir monitoring purpose. The geophone arrays are positioned very close to the surface, ranging from 15-950.75 feet depth. A finite-difference 2D elastic modelling is conducted to understand the effect of shallow geophone placement to downgoing and upgoing P waves. This study demonstrates that the downgoing P waves are starting to refract at offset 3300 feet therefore the accuracy of downgoing P waves arrival time at offset farther than 3300 feet decrease. This study also concludes that the final 3D VSP image from the transformation process at the reservoir level is adequately correlatable with the surface seismic image. However, the 3D VSP image has lower resolution since the high frequency content of the VSP data is severely attenuated by the low-velocity zone at the near surface.

Abstract: This paper discusses a novel approach to transform the recorded upgoing P-waves to a flat datum on the surface based on the wave-equation. The transformation results in wavefields which appear as if they are recorded at the surface. The result of the transformation enables the application of CMP-based processing techniques, such as velocity analysis, Normal Moveout (NMO), and static correction to 3D VSP data. A good transformation process requires an accurate measurement of downgoing P wave arrival time at each geophone. The 3D VSP data used in this study was acquired for reservoir monitoring purpose. The geophone arrays are positioned very close to the surface, ranging from 15-950.75 feet depth. A finite-difference 2D elastic modelling is conducted to understand the effect of shallow geophone placement to downgoing and upgoing P waves. This study demonstrates that the downgoing P waves are starting to refract at offset 3300 feet therefore the accuracy of downgoing P waves arrival time at offset farther than 3300 feet decrease. This study also concludes that the final 3D VSP image from the transformation process at the reservoir level is adequately correlatable with the surface seismic image. However, the 3D VSP image has lower resolution since the high frequency content of the VSP data is severely attenuated by the low-velocity zone at the near surface. Keywords: 3D VSP, datum transformation, wave-equation, upgoing waves, modelling Abstrak: Makalah ini membahas sebuah pendekatan baru untuk mentransformasikan rekaman gelombang P pantul atau upgoing wavefields ke sebuah datum datar di permukaan berbasis persamaan gelombang. Transformasi tersebut menghasilkan gelombang seismik yang seolah-olah direkam pada pseudo-receiver di permukaan bumi. Hasil transformasi ini memungkinkan untuk diterapkannya pengolahan data berbasis Common Mid Point (CMP), seperti analisis kecepatan, Normal Moveout (NMO), dan koreksi statik pada data VSP 3D. Proses transformasi yang baik membutuhkan pengukuran waktu tiba gelombang langsung P yang akurat pada setiap geophone. Data VSP 3D yang digunakan dalam studi ini diakuisisi untuk keperluan pengamatan reservoir. Rangkaian Geophone VSP ditempatkan sangat dekat dengan permukaan yaitu mulai dari kedalaman 15 hingga 950.75 feet. Sebuah pemodelan elastik 2D berdasarkan metode beda -hingga dilakukan untuk mengetahui efek peletakan geophone di dekat permukaan ini terhadap penjalaran gelombang langsung dan gelombang pan-tul P. Studi ini menunjukkan bahwa gelombang langsung P mulai mengalami refraksi pada offset 3300 sehingga mengurangi keakuratan pengukuran waktu tiba gelombang P pada offset jauh. Selain itu, studi ini juga menyimpulkan bahwa citra VSP 3D yang dihasilkan dari proses transformasi ini berkorelasi cukup baik dengan data seismik di zona reservoir. Namun, resolusi citra VSP 3D yang dihasilkan lebih rendah dari data seismik permukaan akibat atenuasi frekuensi tinggi oleh zona kecepatan rendah di dekat permukaan. Kata kunci: 3D VSP, transformasi datum, persamaan gelombang, gelombang pantul, pemodelan

INTRODUCTION
Vertical Seismic Profiling (VSP) observes downgoing and upgoing seismic waves that propagate through the subsurface by placing the receiver within a drilled wellpad and the source near the surface. This method allows the analysis of not only the downgoing waves but also the reflected upgoing waves and expands the product of borehole seismic analysis from one dimensional time-depth relationship to 2D or 3D image of the subsurface in the vicinity of the wellbore. Typically, VSP image has higher resolution compared to surface seismic one. This is because the seismic waves travel through the attenuative near surface only once.
In the early years of VSP imaging, the 2D and 3D VSP data were imaged using the VSP-CDP mapping technique. Though this technique is proven to be useful in quick analysis of VSP image, it is still not accurate for area in which complex velocity presents. Even though more robust 3D VSP processing technique is vital to handle complex geology, the development of 3D VSP processing algorithm is not as rapid as in surface seismic data. This paper discusses a novel approach to 3D VSP processing based on the concept of wave-equation datuming presented by Berryhill (1984). This approach allows the recorded upgoing waves to be transformed to a new flat datum on the surface. This datum transformation process is known as upward-continuation process. The transformation results in upgoing wavefields which appear as if they are recorded on the surface. Thus, the transformed upgoing wavefields can now undergo established surface seismic processing techniques such as NMO, velocity analysis, and static corrections. The transformation process requires a good first arrival measurement of the downgoing P waves.  The 3D VSP data used in this study has some disadvantages since the geophones are situated very shallow and might be exposed to low-velocity zones. The downgoing waves coming from far offset shots have the risk of being refracted, especially at shallow geophones. The refraction of downgoing P waves decreases the accuracy of first arrival picks. Therefore, an elastic finite difference modelling is performed to understand at what offset the downgoing P waves start to refract.

3D VSP Data Acquisition
The 3D VSP data was acquired simultaneously with the surface seismic data using dynamite as the source. There are 927 shots that finally used for the processing with the closest shot is 350 feet and the farthest one is 8800 feet away from the VSP well. This acquisition is a part of reservoir monitoring of a CO2 -EOR project. The source information used in the acquisition is described in Table 1 below. The receivers used in the 3D VSP survey are GS-One 3C geophone manufactured by OYO Geospace. There are 20 levels of geophone which depth ranges from 15 to 950.75 feet with 49.25 feet interval. The VSP borehole was filled with cement to suppress the tube waves. Tube waves are typically Rayleigh waves that travel down the borehole along the interface between the fluid and the casing in the wellbore.
The interval of the sandstone reservoir target is approximately from 3260 to 3300 feet. The sandstone reservoir truncates to an overlying chalk formation. With the last of the geophone placed at 950.75 feet, there is 2309.25 feet distance between the maximum depth of the geophone array to the top of the reservoir. This large distance might compromise the image quality of the reservoir. An illumination study is conducted to model the quality of the 3D VSP image at the top of the reservoir (Figure 1). This illumination modelling uses the actual source and receiver positions and assumes that the top of the reservoir is flat. Figure 1 shows that the sources are not evenly distributed. The fold map on the top of the reservoir describes the illumination from all sources. The sparse shot distribution on the surface (circled in yellow) creates low fold area at the top of the reservoir (circled in white). A relatively high fold region is situated in the centre of the image area. Farther from the centre, the fold decreases.

Finite-Difference VSP Modelling
To better understand the wavefield recorded by the VSP, an elastic 2D wave propagation modelling is performed. Several 2D synthetic shots with various offset are constructed by applying finite-difference modelling of elastic wave equation in an isotropic media. Solving the elastic wave equation is necessary to simulate the wave propagation in order to understand the effect of shear waves and mode converted P-Sv waves. The algorithm for staggered finite-difference used in this modelling is based on the work of O'Brien following the workflow described by Etgen (1987). The staggered grids approach is more stable and accurate since the newly computed variables are centered between the old variables resulting in reduction of the grid spacing by half (Carcione et al., 2002). The inputs of the model are Vp, Vs, and RHOB logs that contains the sand reservoir interval ranging from 3260 feet to 3300 feet ( Figure 2). Figure 2 shows that the twenty levels geophone arrays are situated within lower velocity zone. To sufficiently image the relatively thin reservoir, a high -bandwidth model is required. The maximum frequency that can be obtained from this model depends on the lowest value of the velocity and model grid (Carcione et al., 2002). Since the value of shear of velocity of the input model is low, the model grid is set to be relatively fine with a five-foot sample of dz and dx. The maximum frequency of this model is defined as: Where F f und is the dominant frequency of the ricker wavelet set to 70 Hz. Giving the dominant frequency of 70 Hz and average P-wave velocity from the log is 9500 ft/s, the vertical resolution of the seismic model will be approximately 34 feet. The VSP is modeled as a walkaway VSP survey with shots buried at 30 feet below the surface and located 0 -8800 feet away from the borehole. The receivers are 20 levels of vertical component geophones ranging from 15 -950 feet of depth with 49 feet interval. The source signature is considered as a zero-phase ricker wavelet from a point source and the boundary condition is designed to reflect the wavefields back to earth instead of absorbing them.
The result of the 2D elastic modelling (Figure 3) shows that the downgoing P-waves is refracted as the offset increases. The figure shows various offset of synthetic shots with receiver depth (ft) in x-axis and time (ms) in y-axis. Typically, in a plot like this, the downgoing waves is expected to be linear with negative dip and upgoing waves with positive dip. However, as the geophones are placed very shallow, shots coming from further offset will start to bend and refracted. Therefore, at further offsets, the downgoing waves will have a positive dip. This phenomenon will cause uncertainty in first break picking. Even though the farthest  offset is 8800 feet, we limit the modeling to 7700 feet offset. This limitation will not affect the understanding that the downgoing P-waves will start to refract after 3300 feet away from the borehole. Figure 3 also reveals the present of high amplitude and low frequency diving waves which mainly affect the shallower geophone. This diving waves might potentially contaminate the upgoing P-waves coming from the reservoir. The upgoing P-wave reflection from the reservoir is not clear due to the 2309.25 feet distance between the last geophone and the top of the reservoir that reduces the signal-to-noise ratio.

3D VSP Processing
Based on the previous modelling, the shallow 3D VSP acquisition geometry has caused several disadvantages that might affect the final image. The goal of the 3D VSP processing is to separate the upgoing-P wavefields from all the recorded wavefields, enhance, and finally image those signals. The processing will use the recorded wavefields only from the vertical (Z) component. The processing flow is described in Figure 4.

Geometry
The processing flow is initiated by loading the geometry information to the trace headers. At this point, first break picking is also performed to check on the geometry process and as reference to direct P-waves energy for tool orientation.

Tool Orientation
In a three component VSP acquisition, the vertical component is always oriented toward the vertical direction. However, the other two horizontal components (H1 and H2) have different orientation from one geophone to another. Therefore, the H1 components is then mathematically oriented towards a fixed reference by using P-wave polarization plane method. Since the geophone array is positioned shallow, most of the upgoing P-wave energy is already recorded on the vertical component. Therefore, the tool orientation step will not the affect the recorded wavefields in the vertical component. The following processing flow will be implemented on the vertical component only.

First Break Picking
First break picking is an essential step in processing VSP with wave-equation datuming approach. The first break picking is performed on the downgoing P-wave energy. However, as offset increases the downgoing P-wave energy is now refracted and affected the first arrival at all geophones. At the offset of 2500 feet away from the well, these refracted waves have dominated and change the dips of the first arrivals. This phenomenon complicates the first break picking process, especially in the offset further than 2000 feet.

Wavefield Separation
The recorded wavefield shown in Figure 5 confirms the observation from the modelling. The near-offset shots are highly affected by the near surface effect. It is proven by the lower frequency content and highly noisy shot records that contaminate the upgoing P-wave information. The goal of the wavefield separation process is capturing these upgoing Pwaves and separating them for the other wavefields. The process is performed by applying a fan-shaped filter that differentiates wavefields based on dip, velocity, and frequency. The downgoing wavefields are recognized as having negative dips while the upgoing ones having positive dips. At this step, it is necessary to know the arrival time of the upgoing energy from the target. Based on the elastic modelling, the arrival time from the reservoir is in the interval of 750 -900 ms. The result of the wavefield separation process is shown in Figure 6.

True Amplitude Recovery
True amplitude recovery is an amplitude correction for geometrical spreading effect. The amplitude recovery is performed by applying simple power of time correction (x) to the traces (t). In this case, the x is 2.

Surface Consistent (SC) Amplitude Correction (AC) and Deconvolution
The goal of applying SCAC is to correct the amplitude due to the near surface and acquisition effect such as source strength and receiver coupling response. The implementation of SCAC involves three main steps, the first is to calculate the RMS amplitude of the input traces within 1400 ms window, measures the gain by decomposing the RMS amplitude to source and receiver component, then apply the gain to the whole traces. The SC Deconvolution is applied to improve the bandwidth of the seismic data and convert the data to zero phase. The process is similar to SCAC where the seismic data is decomposed into four components, which are source, receiver, offset, and impulse response of the earth. Then, the power spectra from source and receiver component are used to design the deconvolution operator. In this case, a 140 ms operator length is applied.

Upward-Continuation and Survey Re-Gridding
Berryhill (1979) defines wave-equation datuming as the name given to upward and downward continuation of seismic wavefields from one reference surface to another. The datum here defines as the surface where the source and receiver are located. The purpose of wave-equation datuming is to simulate the data as if it would have been recorded at the output datum. (Berryhill, 1979) initially applied the concept of wave-equation datuming to post-stack data then expands the application to common-source and commonreceiver gathers.
The application of the upward continuation concept in 2D and 3D VSP processing has been implemented by Fuller et al. (2008). This algorithm enables the upgoing wavefields to be upward-continued to a new datum in the surface where pseudo-receivers are located. They explain that the extrapolated seismic signal at the pseudo-receiver location (P (t)) is the summation of seismic signals recorded at each downhole geophone (Ri(t)) which is delayed by the travel time from point P to the respected downhole geophone (fi). Figure 7 illustrates equation (3) where Ri is the recorded arrival time of upgoing P-wave at each downhole geophone and fi is equal to the first arrival time of the downgoing wave to each downhole geophone. The implementation of equation (3) to 2D or 3D VSP data transforms the upgoing P-wavefields to the extrapolated data that can be processed as surface seismic data (Figure 8). Figure 8 shows the extrapolated upgoing wavefields for different offsets. The deeper reflections, especially in near offset, are heavily influenced by low frequency noise. These extrapolated wavefields will be re-gridded to match the 3D surface seismic survey. The grid configuration is 82.5 feet inline spacing, 41.25 feet crossline spacing, azimuth of the crossline is 198.01 degree with flat datum and replacement velocity 6000 ft/s.

Velocity Analysis, NMO, and Static Correction
After re-gridding, the 3D VSP survey is now ready to be processed with well-established CMP-based processing technique such as velocity analysis, NMO, and static corrections. The initial velocity analysis and NMO are conducted using the RMS velocity from surface seismic. Then, residual stat- ics is computed using stack-power maximization process as described by Yilmaz (2001).
The residual statics for source and receiver position are calculated by determining time-shift, which is based on the cross correlation between traces in CMP gathers and flattened pilot trace in a time-gate 900 -1000 ms. This time shift correction is then decomposed into source and receiver component and the applied to the original data. The velocity analysis, NMO, and static corrections are applied iteratively.

Trim Statics
Trim statics applies a single time-shift to each individual trace by cross -correlating the input trace with a model trace. The trim statics model is obtained from the stack after five iterations of velocity analyses and static corrections. Figure 9 illustrates the improvement of the stack before and after trim statics.

Pre-Stack Time migration, Stack, and Post-stack Enhancement
The 3D VSP data is migrated using Kirchhoff Pre-Stack Time Migration algorithm using RMS velocity from final iteration of NMO. The Kirchhoff Pre-Stack Time Migration works in offset domain and uses the actual ray path from every source to receiver and constructs a diffraction surface based on these ray paths. The Kirchhoff migration sums the contribution of all data within the migration aperture of 80•. After stacking the migrated gathers, F-XY deconvolution filter is applied to suppress random noise. F-XY deconvolution applies a Fourier transform from time and distance to the frequency and distance domain. Typically, 3D VSP image will be resulted in a cone-shaped image ( Figure 10(lower)). However, due to the heavy contamination of near surface noise, the 3D VSP image is extensively muted. The goal of splicing the VSP image into the surface seismic is to match the reflection character of VSP image into the surface seismic.

RESULTS AND DISCUSSION
Based on the power spectra observation, the surface seismic shows higher frequency content than the 3D VSP image ( Figure 10). Consequently, in this study the surface seismic will have better vertical resolution than the 3D VSP image. The lower frequency content of the 3D VSP image is an indication of severe attenuation of the high frequency content of the VSP data since most of the geophone are situated very shallow and heavily affected by the low velocity zone in the near surface as depicted in Figure 2. The splice between the 3D VSP image and the surface seismic shows moderate to poor correlation especially from deep reflectors (Figure 11). This is due to the low signal-tonoise ratio of the upgoing energy coming from the deeper reflectors. The waves need to travel long distance upward to reach the last geophone. Ideally, the last geophone will be positioned very close to the target reflectors.
Despite all the limitations, the splice is still useful to identify image of the target reservoir in the vicinity of the borehole. The reflection from the distinctive top chalk formation and the reservoir is adequately correctable even though the dip of the events is not perfectly match. The low signal-to-noise ratio at the reservoir interval complicates the velocity analysis at the reservoir interval. This different velocity contributes to the wrong dip at target level.

CONCLUSIONS
A 3D VSP data with shallow geophone array is processed with wave-equation datuming approach. In addition to that, a two-dimensional elastic wave equation and illumination modelling are performed to better understand the effect of shallow geophone placement and acquisition geometry of the recorded wavefields. Based on those observations, we have reach to the following conclusions: (i) The wave-equation datuming that was performed by upward-continuing the recorded upgoing P-waves energy to the pseudo-receiver locations on the flat surface has successfully transformed the VSP data to extrapolated wavefields. These extrapolated wavefields are then processed with CMP-based processing technique such as velocity analysis, NMO, and statics correction commonly known for processing surface seismic data. (ii) This VSP processing approach requires good first arrival measurements in order to transform the recorded upgoing wavefields to the surface. However, the result of 2D elastic modelling based on staggered-grid finite difference approach suggests refracted downgoing wavefields coming from far offset shots. The same phenomenon is also observed in the real data that complicates the first arrival measurement. (iii) The splice of 3D VSP image into surface seismic shows an adequately correlatable event in the reservoir interval. The VSP image is situated up-dip of the wedge. (iv) Shallow geophone position contributes to the inferior quality of the final 3D VSP image since the high frequency content of the recorded upgoing wavefields is severely attenuated.