A comprehensive seismic project including conventional and multi-component surface seismic and borehole seismic was acquired in a producing field in northeast British Columbia, Canada. Delineating productive zones within the fractured shale reservoir was the primary goal of the seismic recorded in the project. The borehole seismic includes a conventional P-wave 3D VSP which was expected to provide useful images, and provide borehole calibration information for accurate processing of the surface seismic data.

This article covers some unique solutions to common problems in land 3D VSP processing. These problems needed to be solved before the basic goal of providing high resolution images could be met. This 3D VSP data set offered us the opportunity to make in-situ measurements of velocity, anisotropy and Q. Some of the most important multiple generators in the field were identified, and a de-multiple technique for VSP images was developed. We then looked at determining fracture orientation information from this same data set. The results from a feasibility study of elastic Full Waveform Inversion (eFWI) using data from this 3D VSP are also shown in this article. eFWI plays an important part in the future of borehole seismic, as it allows us to extract the most important medium properties directly from the unprocessed vector data recordings.

Fig. 01
Figure 1. Project location in northeast British Columbia, Canada. The map shows source locations for the 3D VSP as red dots, receivers are black dots. Map courtesy Google Earth.

3D VSP background

The 3D VSP dataset was acquired for Nexen Energy ULC. in Northeast British Columbia, Canada (Figure 1). The 3D VSP was acquired simultaneously with the conventional surface seismic program. Vibrator sources from the surface seismic survey which were within 2500 m from the wellhead were used for the 3D VSP analysis. The three-component down-hole receiver array had 135 depth levels extending from below the top reservoir to surface (Figure 2).

Fig. 02
Figure 2. Well logs with 3D VSP receiver depths plotted in right track.

Processing and imaging

Prior to processing the downhole data using a VSP processing flow, noise attenuation and amplitude compensation methods used in surface seismic processing were adapted for application to the VSP data. Anomalous Amplitude Attenuation (AAA) is commonly used as a surface seismic noise reduction technique. AAA aims to attenuate random noise by transforming the processing gather into the frequency domain and applying a spatial median filter. Frequency bands that deviate from the median amplitude by a specified threshold are either zeroed, or replaced by good frequency bands interpolated from neighboring traces. By using AAA, the high amplitude noise was suppressed from the raw VSP data with minimal loss of signal.

With land data, due to variability of the near-surface conditions and inconsistent coupling of the vibroseis source, it is typical to observe a rather random scatter in the recorded amplitudes versus offset. Surface Consistent Amplitude Compensation (SCAC) compensates for this amplitude variation from shot to shot. SCAC is commonly used in land surface seismic processing and compensates for shot, detector and offset amplitude variations that are caused by acquisition effects and are not a consequence of the subsurface geology. For this land 3D VSP, SCAC was used by applying the compensation along the offset and shot axis (Figure 3).

Fig. 03
Figure 3. Walkaway line, Z axis, input (A) and output (B) after AAA and SCAC. Variable source to source amplitudes have been compensated. The shallow noisy receivers were deleted before SCAC.

To generate upgoing waves as input to migration, traditional median filtering enhanced and removed downgoing waves from the vertical channel wavefields. To maximize the offsets used in the p-wave migration, shear noise was reduced using a Radon filter prior to deconvolution. An inverse Q filter was also applied prior to deconvolution.

Q is commonly estimated from the downgoing arrivals of the borehole data. This data set allows us to estimate Q, using a spectral ratio technique, from the near surface to the reservoir. A log of effective Q was determined by setting a shallow receiver as reference and calculating Q between the reference and all other trace spectra. Q values from one near shot are plotted in Figure 2. An average Q value of 100 was estimated from several near offset sources, and applied as an inverse filter. The data were deconvolved using standard up-by-down deterministic VSP deconvolution.

Fig. 04
Figure 4. 2D east-west slice through the wellhead from P-wave GRT image volume from single component processing.

The deconvolved upgoing waves were migrated using GRT (Miller et al., 1987) migration. Figure 4 shows a 2D slice extracted from a 3D volume centered on the wellhead. In the same figure, Vs from the dipole shear log, density and Vp from compressional sonic logs are plotted on top of the image. The reservoir shale does appear to generate an interpretable reflection on these Pp images at about 2,100m depth. More work is currently being done on a quantitative interpretation of the images.


This 3D VSP dataset has a geometry ideal for the estimation of anisotropic parameters. We have receivers in situ, from the reservoir level to the near surface. The surface seismic imaging project could benefit from a more accurate estimate of anisotropy, or at least a calibration of the anisotropy being used in the NMO correction and imaging steps. Estimation of VTI anisotropy is an integral part of the imaging model calibration.

For the depth imaging of our 3D VSP, an initial v(z) profile at the well location was generated from the sonic logs. The picked direct 3D VSP travel times were iteratively fit by perturbing interval anisotropy parameters and velocities to determine best-fit anisotropy. The resulting VTI model is shown in Figure 5. In this case, the gamma VTI measurement of shear anisotropy from sonic measurements, together with 1/Vp and a compaction law constraint were used to constrain the magnitude of anisotropy in the inversion. Using these constraints, a depth-variable anisotropy model was created.

Fig. 05
Figure 5. The calibrated VTI model and residuals versus offset (m). Top: Vp, Vs, density in the left 3 tracks. Epsilon (black) and delta (red) in track 4 and shear sonic gamma in the rightmost track. Bottom: Residuals versus offset for the isotropic model (red), and the VTI model (green).

The velocity model must explain both the surface seismic and borehole seismic data. To confirm that the model explains the borehole data, we ray-traced direct arrivals through the VTI model, and compared the modeled arrival times with the observed direct arrival times. The residual time should be near zero for a perfect model. Figure 5 shows the residual times in red for an isotropic version of the model, and green for the VTI model. Residuals are small and centered on zero for the VTI model.


VSP deconvolution is typically used for zero phasing the data and in the attenuation and identification of multiple energy. In this survey, a shot near the wellhead is the equivalent to a standard Zero-offset VSP (ZVSP). Inside and outside corridor stacks of the deconvolved ZVSP data can be used to indentify interbed (or internal) multiples and their generators (Lines and Burton, 1996). The outside corridor stack should be multiple free along the receiver array and the inside corridor stack should contain both interbed multiples and primaries. The blue-shaded area on the deconvolved upgoing from a near offset shot shown in Figure 6 defines the data that went into the inside corridor stack. The pink-shaded area shows the data that went into the outside corridor stack. An interbed multiple can be clearly seen on the inside corridor stack between 1.3s and 1.4s. In depth, this is very close to the top of the reservoir at about 2100m. The multiple is observed to be interfering with the target reflection. The inside corridor can be targeted to examine multiples generated by high reflectivity interfaces.

Fig. 06
Figure 6. Deconvolved upgoing from a near offset shot on left, outside corridor and inside corridor on far right. The blue box indicates the time location of an interbed multiple.

The interbed, or pegleg, multiples were generated from the reflectivity estimated from the compressional sonic log. This is an excellent tool for locating multiples in seismic sections. These synthetics were used confirm the multiples detected by our inside/outside corridor stack methodology. Our multiple of interest was confirmed on the multiple only synthetic, and the primary plus multiple synthetic.

Fig. 07
Figure 7. Migrated P-wave image before and after the VSP de-multiple process. The red arrow indicates the location of the interbed multiple being attenuated.

From the inside and outside corridor stacks, we used adaptive subtraction to estimate an interbed multiple model for this 3D VSP. The multiple model is the right-most stack in Figure 6. After estimating the interbed multiple model, and using a 1D assumption, the multiples were adaptively subtracted from the 3D VSP image cube. The initial results of this methodology demonstrate the promise of the technique (Figure 7) in attenuating the interbed multiple interfering with the top reservoir.

Fracture Orientation from Shear Waves

The downgoing mode-converted shear waves recorded in a walkaround VSP can be informative on fracture or stress orientation. This technique requires a circle of constant offset sources, with the horizontal channels oriented to the radial and transverse directions. Downgoing mode-converted shear waves recorded on a walkaround VSP will split upon encountering a zone of fractured rock. As described by Horne (2003), the presence of energy on the transverse channel as the source moves in a circle around the well is a good indicator of the shear splitting caused by azimuthal anisotropy. There are certain azimuths where the energy on the transverse channel dims as the source moves around the well. These azimuths correspond to vertical symmetry planes aligned with the fracture plane and fracture normal directions

From the full 3D VSP data set, a circle of sources was extracted for the walkaround analysis. An optimal offset distance where mode conversions were consistently generated was chosen. On the radial channel data we observed a strong downgoing shear wave passing through our zone of interest. This shear wave was analyzed on both the radial and transverse channels for evidence of splitting as it propagated past the receivers. Ideally, the shear wave being analyzed should have it’s origin close to the zone of interest, to avoid the possibility of multiple shear splittings.

Fig. 08
Figure 8. Radial and transverse channels for one receiver, flattened at 200ms time. R/T ratio plotted as graph on transverse channel, and as a polar plot below.

The Root Mean Square (RMS) amplitude of the downgoing shear was estimated on both the radial and transverse channels, in a window around the downgoing mode converted shear wave. The R/T ratio should go to a maximum where the transverse channel is dimmest. Figure 8 shows clear dimming and brightening in the amplitudes on the transverse channel, every 90 degrees for the receiver shown. The dimmings and brightenings on the transverse channel show up as maxima and minima on the R/T graph plotted on top of the transverse channel data. The R/T ratio numbers are very large, indicating the effect we are looking at is very weak, but is still measurable. The maxima on the R/T polar plot indicate the symmetry axis orientation directions. Where the horizontal data were good quality, there was a fairly consistent trend of maxima in the R/T ratio, approximately in the east-west (70 to 80 degrees azimuth), and north-south azimuths.

Fig. 09
Figure 9. From left to right: Gamma Ray (GR), fast and slow shear sonic slowness, fast shear azimuth, Vs1/Vs2 and Sh gamma.

Figure 9 shows sonic shear anisotropy results over the bottom section of the well. The fast shear azimuth is seen to be consistent and close to what is seen (N70E) in the converted shear R/T analysis. Also shown in Figure 9 is the sonic “gamma_tiv” curve, which is the horizontally polarized shear (Sh) VTI anisotropy number, and since it is significant it is of interest to consider an orthotropic reservoir model. Gamma is seen to be lower but positive where HTI is large and higher above and below. Thus a VTI/orthorhombic/VTI layering is suggested.

Elastic Full Waveform Inversion

Even more valuable than the high resolution images provided by this 3D VSP would be a direct estimation of the medium parameters Vp, Vs, density, epsilon and delta in a volume around the well. The recent developments in elastic Full Waveform Inversion (eFWI) detailed in Podgornova et al, 2014, make this possible. The method is currently 2.5D, uses the vertical channel and the in-line horizontal channel, and is fully anisotropic. EFWI as described here uses elastic wave equations in anisotropic VTI media to model the wave propagation. Five parameters are inverted: compressional velocity, ratio of the compressional and shear velocities, density, and the Thomsen anisotropic parameters epsilon and delta.

The inversion minimizes the misfit between the modeled data and the real data by updating the model’s medium parameters. The starting model is the borehole calibrated VTI model generated from the well data. The modeled data are saved as synthetic seismograms. The FWI algorithm performs a simultaneous inversion for medium parameters and the source wavelet. Inverting for the source wavelet is a critical component of FWI with land data due to the variability in surface conditions.

Fig. 10
Figure 10. Shot gather at 1050m offset from the well. A and B are the real data, C and D are the synthetic data. Yellow arrow indicates upgoing p-waves, red indicates downgoing mode-converted shear, green is upgoing mode-converted shear. Adapted from Podgornova et al, 2014.

The inversion was run on an extracted 2D line of sources using frequencies up to 40 Hz. The starting model was generated from a smoothed version of the model shown in Figure 5. Synthetic seismograms generated in the inversion are compared to the input data in Figure 10. The input data vertical (A in Figure 10) and horizontal traces (B) for a medium offset shot compare extremely well to the modeled data (C and D) after 150 iterations of inversion. There are complex wavefields in these recordings which are well recovered in the modeled data. Strong upgoing and downgoing mode conversions are evident in the input data, and are accurately reproduced in the synthetic data. This gives us confidence in the medium parameters recovered by the inversion.

Fig. 11
Figure 11. Inverted Vp and Vs sections, receivers are white circles.

The eFWI updated the initial model and recovered reasonable layered Vp and Vs structures (shown in Figure 11) and similar results for epsilon and delta anisotropy parameters (not shown). In Figure 12, all the high-contrast and many small-contrast layers seen in the Vp and Vs sections agree with the sonic data with respect to the frequency content of the data used in the inversion.

Fig. 12
Figure 12. Extracted logs from near the wellbore. Black curves are the sonic data, green are smoothed sonic data used as the starting model, blue curves are extracted from the inverted model at the well location.

The eFWI results from this feasibility study prove the viability of the method in land 3D VSPs (or walkaway VSPs). We were able to estimate the elastic parameters of the medium directly from the unprocessed vector data without the need for more complicated imaging/inversion workflows. While the method is currently essentially 2D, work is underway to progress to 3D eFWI.


To create high resolution images from this land 3D VSP, land-specific algorithms such as AAA and SCAC were used to properly pre-process the data. Q compensation was incorporated in the pre-processing. The imaging model needed to include strong anisotropy, an important learning from this project.

3D VSPs offer more than high resolution images. VTI parameters were estimated for depth imaging of the 3D VSP. These parameters are also useful for surface seismic imaging. HTI parameters were estimated from shear splitting observed on the downgoing mode-converted shear waves and from sonic measurements. Both measurements were consistent with an ENE fast azimuth. To attenuate interbed multiples, the traditional method of inside-outside corridor stacks was used to first identify interbed multiples. Then a newly implemented adaptive methodology was used to attenuate the interbed multiples in the image volumes.

The data from this 3D VSP were used in a very successful feasibility study for eFWI. The velocity and anisotropy structure in a 2D section through the well were accurately recovered. We believe eFWI is going to be a key component to the future of borehole seismic as it allows direct computation of the key medium parameters from the unprocessed vector data. eFWI should be seen as an enabling technology which will in the future extend the usefulness of borehole seismic data into areas such as effective stress estimation and pore pressure prediction.



We gratefully acknowledge Nexen Energy ULC for permission to publish this work. We also acknowledge Marwan Charara and the eFWI developers at Schlumberger Moscow Research Center.


About the Author(s)

Allan Campbell graduated with a B.S in Geophysics from the University of Calgary in 1990. Allan has worked at Schlumberger in various roles in the borehole seismic business since graduation. He spent 7 years in Calgary acquiring and processing VSP data. In 1997, Allan moved to Houston to manage the VSP processing business there. Since 1997, he has been closely involved in the development of many advanced borehole seismic techniques. Allan’s current role is Borehole Seismic Marketing Manager for Schlumberger PTS. His major geophysical interest is the development of new applications for borehole seismic in all its forms. He is a member of the CSEG, SEG, APEGA, and the EAGE.

Scott Leaney is a Geophysics Advisor in Schlumberger, based in Houston. After a Masters in Geophysics at UBC (B.Sc. Manitoba) he joined SLB in 1988 in Engineering in France developing software for 3C VSP processing and AVO modeling. From 1993-1997 he was Southeast Asia Area Geophysicist for Wireline Operations based in Indonesia; 1998-2002 manager of a group dedicated to developing integrated seismic processing for WesternGeco, based in Gatwick, England; since 2003 VSP processing development and microseismic R&D Advisor in Houston focusing on ray trace modeling and inversion in anisotropic media. (leaney@slb.com).

Jitendra Gulati is a Principal Geophysicist at Schlumberger. He has extensive experience in 3D and 4D VSP in several basins of the world. He has published articles in peer reviewed journals and filed patents in downhole seismic. He has also been actively involved in organizing industry workshops including SEG’s first downhole seismic workshop focused on 3D VSP, and the first fiber-based distributed acoustic sensing (DAS) workshop. At Schlumberger, he currently holds the responsibility of business development for downhole microseismic, VSP and crosswell seismic. He has received geophysics degrees from Indian Institute of Technology, University of Calgary and Stanford University.

Olga Podgornova graduated from Novosibirsk State University in Russia. She received her PhD in Keldysh Institute of Applied Mathematics in Moscow and joined Schlumberger Research in 2009. Olga’s interests include numerical methods for wave propagation problems and inverse problems in seismics.

Jennifer M. Leslie-Panek graduated with a Ph.D. from the University of Calgary specializing in Depth Imaging and Anisotropy. Worked in Seismic Data Processing for several years, processing land and marine data in both time and depth, before moving to Nexen in 2006. Started in technical support as a processing specialist and researcher, moved into exploration for the Canadian New Ventures business unit, and then the Shale Gas business unit in 2009. Worked on seismic projects in Colombia, Yemen, Nigeria, Norway and Western Canada. Extensively involved in research into land inter-bed multiple elimination, VSP and microseismic technology. Actively involved with the CSEG in various roles, currently as a Director on the CSEG Foundation Board.

Eric von Lunen was born and raised in western Pennsylvania, USA. He attended Penn State University earning a BS in Geological Science-Geophysics with extensive minor course study in physics and mathematics. He earned a MS in Geology from the University of Texas.

His career experience includes extensive technical work (and travel) in international settings such as Latin America, Middle East, Africa, Eastern Europe and Southeast Asia. He presently lives in Calgary, Alberta.

Eric has done extensive consulting and contracting for geophysical application research for major and large petroleum independents such as Exxon-Mobil, Aramco, Hunt Oil, Marathon and YPF-Maxus among others. In 2006, sold consultancy and worked as chief geophysicist at Kerogen Resources planning and implementing geophysical strategies in tight reservoir and shale gas plays. In 2014 he was working as senior geophysical advisor and team lead on shale gas projects with secondary responsibilities including 3-D/3-C seismic technology and microseismic research and innovation, plus extensive technical mentoring. These technology interests are coupled with a strong interest in geophysical applications for resource and reserve estimation.

Affiliations include SEG, CSEG, EAEG, AAPG, and SPWLA.

Eric looks forward to mentoring graduate students on technical subjects at University of Texas, Colorado School of Mines, Queens, University of Calgary, University of British Columbia, and University of Toronto, among others.


Horne, S., 2003, Fracture Characterization from walkaround VSPs: Geophysical Prospecting, 51, 493-499

Prospecting, 51, 493-499 Lines L., Burton, A., 1996, VSP detection of interbed multiples using inside-outside corridor stacking: Canadian Journal of Exploration Geophysics, 32, 2, pp. 113-120

Miller, D.E., Orstaglio, M. and Beylkin, G. 1987, A new slant on seismic imaging: Migration and integral geometry. Geophysics, 52, pp. 943-964.

Podgornova, O., Leaney, W. S., Charara, M., von Lunen, E., Elastic full waveform Inversion for land walkaway VSP data, CSEG Geoconvention 2014 Abstract

Suggested Readings

Chapman, C.H., 2004, Fundamentals of Seismic Wave propagation, Cambridge University Press.

Leaney, W.S., 2008, Polar Anisotropy from Walkaway VSPs: The Leading Edge, October, 2008

Schoenberg, M. and Sayers, C. 1995, Seismic anisotropy of fractured rock: Geophysics, 60, 1, 204–211


Join the Conversation

Interested in starting, or contributing to a conversation about an article or issue of the RECORDER? Join our CSEG LinkedIn Group.

Share This Article