Geophysical Prospecting, 2012, 60, 138–152 doi: 10.1111/j.1365-2478.2011.00978.x

Application of 3D vertical seismic profile multi-component datato tight gas sands‡

Yousheng Yan1∗, Zengkui Xu1, Mingli Yi2 and Xin Wei21Information Technology Center, BGP, China National Petroleum Corporation (CNPC), 1501, Gehua Building, Dongcheng District, Beijing100007, China, and 2Geophysical Research Institute, BGP, China National Petroleum Corporation (CNPC), Jiaxiu Road, Zhuozhou, Hebei072751, China

Received March 2010, revision accepted April 2011

ABSTRACTDue to the complicated geophysical character of tight gas sands in the Sulige gasfield ofChina, conventional surface seismic has faced great challenges in reservoir delineation.In order to improve this situation, a large-scale 3D-3C vertical seismic profiling (VSP)survey (more than 15 000 shots) was conducted simultaneously with 3D-3C surfaceseismic data acquisition in this area in 2005. This paper presents a case study on thedelineation of tight gas sands by use of multi-component 3D VSP technology. Twoimaging volumes (PP compressional wave; PSv converted wave) were generated with3D-3C VSP data processing. By comparison, the dominant frequencies of the 3D VSPimages were 10–15 Hz higher than that of surface seismic images. Delineation ofthe tight gas sands is achieved by using the multi-component information in the VSPdata leading to reduce uncertainties in data interpretation. We performed a routinedata interpretation on these images and developed a new attribute titled ‘CentroidFrequency Ratio of PSv and PP Waves’ for indication of the tight gas sands. Theresults demonstrated that the new attribute was sensitive to this type of reservoir.By combining geologic, drilling and log data, a comprehensive evaluation based onthe 3D VSP data was conducted and a new well location for drilling was proposed.The major results in this paper tell us that successful application of 3D-3C VSPtechnologies are only accomplished through a synthesis of many disciplines. Weneed detailed analysis to evaluate each step in planning, acquisition, processing andinterpretation to achieve our objectives. High resolution, successful processing ofmulti-component information, combination of PP and PSv volumes to extract usefulattributes, receiver depth information and offset/ azimuth-dependent anisotropy inthe 3D VSP data are the major accomplishments derived from our attention to detailin the above steps.

Key words: 3D, Centroid frequency, Multi-component, Tight gas sand, Verticalseismic profiling.

‡This paper is based on expanded abstract T036 presented at the 71st

EAGE Conference & Exhibition Incorporating SPE EUROPEC 2009,8–11 June 2009 in Amsterdam, the Netherlands.∗E-mail: yanyousheng@cnpc.com.cn

INTRODUCTION

The Sulige gasfield is located in the Erdos basin in centralChina. In this area, the target formations are Permian in ageand quite thin (generally less than 20 m), with poor poros-ity and low permeability deposited in a deltaic sedimentary

138 C© 2011 BGP/China national Petroleum Corporation

3D vertical seismic profile multi-component data 139

environment. In general, the P-wave impedance of dry sandsis higher than that of cap beds (shale). However, if thesand is saturated with gas, its P-wave impedance decreasessuch that the contrast of P-wave impedance between thesand and shale decreases. Conventional P-wave explorationfaces challenges in detecting and characterizing this type ofreservoir.

2D surface seismic technology is the primary tool employedfor gas exploration in the Erdos basin. Surface conditions arevery complicated because of the loess plateau, this makes 3Dsurface seismic acquisition logistically and technically diffi-cult. The local geophysicists and geologists have made achieve-ments in the delineation of these tight gas sands using con-ventional 2D seismic technology. However, the limitations ofconventional 2D surface seismic become more apparent asthe exploration continues. Methods such as amplitude varia-tion with offset (AVO) analysis and prestack inversion havebeen applied in recent years with limited success, the poorsurface and subsurface conditions have reduced their effec-tiveness. This is because the reliability of seismic attributes,particularly those related to shear modulus, as estimated byAVO analysis or prestack inversion is dominated by the lackof fidelity in the raw data recorded by a vertical sensor inconventional seismic surveys.

In recent years, many applications of 3D VSP technologyhave been presented. However, only a few of them showedresults based on the integration of PP- and PSv-waves. Someof the typical case studies are 3D VSP PP- and PS-imagingused for structural interpretation onshore (Gomes andRonholt 2005), 3D VSP PP- and PS-imaging for carbonatereservoirs (Yan et al. 2007), etc. Most of the authors triedto extract as much shear-wave information as possible fromthe 3D-3C VSP data to combine with P-waves for reservoiranalysis but ultimately had to give up because of a series ofdifficulties relating to data processing and interpretation. As aresult, the multi-component attributes present in 3D-3C VSPdata have not been fully applied to reservoir analysis.

To solve some of the problems in tight gas sand reservoir ex-ploration in this area, a multi-component exploration projectwas proposed and sponsored by the China National PetroleumCorporation. In this paper, we will only discuss the applica-tion of 3D-3C VSP data to the tight gas sands. Our workincluded optimizing geometry design in data acquisition, prin-cipal methods in data processing and new multi-componentattributes applied in data interpretation. By combining the3D-3C VSP data with geologic, drilling and log data, we com-pleted a comprehensive evaluation on the tight gas sands andproposed a new well location for drilling.

CHALLENGES A ND OBJECTIVES

Geological setting

The Erdos basin in central China is a large cratonic basin withseveral sedimentary cycles. The gas source rocks are in theTaiyuan and Shanxi group of Permian age and contain alter-nating continental and marine sedimentary sequences, whichinclude sands, shales and coal deposits.

The Sulige gasfield is located in the middle of the Erdosbasin and is the largest onshore gasfield in China, its locationis shown in (Fig. 1). The major gas-bearing sands are in TP andTC formations respectively in the Permian and Carboniferous(TP: seismic reflection horizons in the Permian; TC: seismicreflection horizons in the Carboniferous). The TP8 forma-tion in this area consists of shale and channel sands depositedin a braided drainage pattern in the delta plain. The aver-age porosity and permeability in a single layer are 12∼15%and 0.15∼2 mD, respectively. The Shan1 formation in TCis generally deposited in a meandering sedimentary environ-ment with thin layers (generally less than 20 m). Its averageporosity and permeability in a single layer are 7∼10% and0.13∼0.16 mD, respectively. The alluvial and deltaic depo-sitional environments are responsible for the development ofthe stratigraphic traps. All of the favourable configurationsof the hydrocarbon source rock, reservoir, cap bed, migrationpath, trap and preservation form the gas reservoirs.

Local logging and drilling results indicate that the targetformations are generally thin (20–30 m), with poor porosity(less than 10%), low permeability (less than 1 mD) and stronglateral variability. Figure 2 shows the distribution of channelsands tied with several wells that is derived from interpretation

Figure 1 Location of the 3D VSP project.

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

140 Yousheng Yan et al.

Figure 2 Well-tie interpretation from log-ging data. Gamma ray logs are the red curveand acoustic logs are in blue.

Figure 3 Attenuation estimated from zero-offset VSP data.

of log and drilling data. In this figure, gamma ray logs arethe red curve and acoustic logs are blue. From this figure it isobvious that the connectivity of the sands is poor at best. Thesereservoir features, as well as their geophysical characteristicssuch as velocity, impedance and stratification indicated by thelogging are complicated not only laterally but also vertically.

Geophysical characteristics and challenges

The Sulige gasfield was found by seismic exploration in 1999.Surface seismic technology has been widely used in the explo-ration and exploitation of this reservoir. As the explorationproceeds more challenges present themselves. The surface con-ditions are complex because of sand dunes, alkali flats, alkalilakes and hard rock outcrops. The depths of the major reser-voirs are generally greater than 3000 m with thicknesses ofaround 10 m. Average porosities and permeabilities are less

than 10% and 1 mD respectively typical of tight gas sands.The work and experience in the last ten years tell us that themajor challenges are as follows:1 Seismic waves are strongly attenuated in the shallow sub-

surface. The major seismic attributes influenced by the at-tenuation, include amplitude, frequency and phase. Thephase variation cannot be visualized easily so we will con-centrate our analysis on the frequency and amplitude por-tion of the wavefield. Figure 3 shows the attenuation offrequency and amplitude versus depth, which was derivedfrom the first breaks of zero-offset VSP data recorded ata depth range of 480–3230 m with a receiver spacing of10 m. The vertical axis stands for depth and horizontalaxes represent the frequency and normalized amplitude.The horizontal dashed line indicates a depth of 1200 m. Thetwo vertical dashed lines around 36 Hz and 0.13 (normal-ized amplitude) are baselines in frequency and amplituderespectively. The following discussion on the attenuationwill focus only on the absorption and ignores geometricspreading attenuation.It can be easily found in the figure that the dominant fre-quency of first breaks in the P-wave data rapidly decreasesfrom 46 Hz at 480 m to 36 Hz at 1200 m in depth (at-tenuation of 11.6 Hz/1000 m). The variation below thedepth of 1200 m trends to a constant (attenuation of0.5 Hz/1000 m). This means the main absorption in thefrequency is from the layers above 1200 m in depth. If weuse a linear fit to the frequency variation in the depth inter-val of 0–1200 m, the dominant frequency at 0 m is around50 Hz and based on this assumption, about 30% of thehigh-frequency components of the seismic waves propa-gated in these layers are absorbed. What we would like toemphasize is that the attenuation shown in this figure is

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

3D vertical seismic profile multi-component data 141

derived from one way propagation instead of two waypropagation, because the calculation is performed on thefirst breaks of a zero-offset VSP data set. Therefore, forsurface seismic data acquisition, the attenuation should bedouble or more. This is probably why the dominant fre-quency of surface seismic data in the target zones is gen-erally less than 30 Hz. Even though some of the high-frequency components attenuated in the shallow subsur-face can be compensated and enhanced by inverse Q fil-tering, deconvolution methods etc. in data processing, onthe whole this compensation and enhancement is still un-satisfactory. Hence, methodologies to acquire the surfaceseismic data with a wider frequency band are very impor-tant to future exploration in this area.Figure 3 also shows the variation of the amplitude versusdepth. Most of the seismic energy emitted on the surface ishighly attenuated by the shallow subsurface above 1200 m.The amplitude attenuation of the one-way propagation isapproximately –27 dB over the depth interval between0–1200 m. Because the depth of the target is far below1200 m, the signal-to-noise (S/N) ratio, frequency and am-plitude decay are important considerations to be addressedin the design of the 3D-3C VSP to be acquired.

2 Conventional seismic attributes, such as P-wave amplitudeand impedance etc., are not sensitive to the tight gas sandsin this area. For gas exploration, the amplitude anomaliespresent in seismic data have been one of the most impor-tant attributes used in data interpretation. In the Suligegasfield, however, the application of amplitude attributesis of limited use in the exploration of the tight gas sands.The difficulty is that the amplitude anomalies from the thingas charged reservoirs are too weak to detect. Poor surfacecondition leads to many problems in statics, attenuationand velocity, which influence the accuracy of the amplitudedisplayed in the resultant seismic sections. As a result, allof these factors increase the uncertainty in the data analysisand interpretation. Figure 4 shows a cross-plot of P-wavevelocity and S-wave velocity that are derived from loggingand drilling data, the horizontal axis stands for P-wave ve-locity Vp, the vertical axis stands for S-wave velocity Vs, themagenta triangles for dry sand, blue squares for shale andred circles for gas-bearing sand. If we divide the figure intofour quadrants as shown by the cyan dashed lines, the datapoints are primarily distributed in the 1st, 2nd and 3rd quad-rants. The majority of the dry sands (magenta triangles) arelocated in the 1st quadrant, which corresponds to high P-and high shear (S-) wave velocity. Most of the shale’s (bluesquares) are located in the 3rd quadrant, which correspond

Figure 4 Crossplot of P- and S-wave velocity from logging data.

to low P-wave velocity and low shear-wave velocity. Thegas-bearing sands are mainly located in the 2nd quadrant,which correspond to low P-wave velocity (Vp) and highshear-wave (Vs) velocity. As a result the P-wave impedanceof dry sand is greater than that of cap bed (shale) if weassume their density is a constant. When the sand is satu-rated with gas, its P-wave impedance decreases. Thus, thecontrast of P-wave impedances between gas-bearing sandsand shale becomes smaller. Obviously, it is not a typicalbright-spot gas reservoir that could be easily delineated byP-wave post-stack seismic attributes. Fortunately, as can beseen in (Fig. 4), the gas-bearing sand and shale can be sep-arated in the shear-wave domain due to big differences intheir shear-wave velocity and impedance (here we assumethe density is constant). This is why conventional methodssuch as instantaneous attributes and post-stack inversion ofP-wave data have been replaced gradually with AVO andprestack inversion in data processing and interpretation.Integration of compressional and shear-wave attributes inAVO analysis and prestack inversion makes the delineationof the gas reservoirs more reliable. However attributes suchas shear-wave impedance and shear modulus, calculated byAVO analysis and prestack inversion, are derived from thez-component data using Zoeppritz equations, or its sim-plified formulas. We believe that these estimated resultsare a pseudo shear-wave attribute, which displays differ-ent characteristics than the shear-waves recorded by multi-component sensors. In some areas, the fidelity of raw datais seriously influenced by poor surface and subsurface con-ditions. As a result, AVO analysis or prestack inversion ofthe data for reservoir analysis contains many uncertainties.So far we can summarize the challenge as two prob-

lems to be solved. The first problem is how to acquire the

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

142 Yousheng Yan et al.

multi-component data with high resolution and high S/N ra-tio. The second is how to apply appropriate multi-componentseismic attributes to delineate the tight gas sands.

Geological objectives

The major geological objectives in the multi-component ex-ploration project are:1 stratigraphic mapping.2 delineation of tight gas sands in TP and TC formations,

particularly in TP7 and TP8.3 proposal of new well locations for drilling.

Because the lateral imaging coverage of VSP data is smallerthan that of surface seismic survey, the objectives for the VSPdata will focus on the area around the geophone deploymentwell.

METHODS

In terms of the major problems described above, the acquisi-tion of multi-component seismic data with higher resolutionand the integration of different types of seismic attributes forgas reservoir analysis are the most important issues to study.We arrived at the solution to these problems in three steps.First, the key point to acquire broadband seismic data wasto reduce the attenuation of seismic waves propagated in theshallow layers. The range of attenuative shallow layers, asshown in (Fig. 3), is in the depth interval between 0–1200 m.Obviously, if the shots and/or receivers can be put belowthe depth of 1200 m, higher-frequency components of thereflection signal can be maintained and recorded. Second,3-component geophones are necessary to record the multi-component information from the target zones. Third is thatwe need to combine the P- and S-wave image volumes to per-form a comprehensive analysis for the gas reservoir interval.

In order to take advantage of the different elastic propertiesof the formations described above, we chose the applicationof 3D-3C VSP technologies to image this type of reservoir.Before we introduce what we did, we will review the benefitsof 3D-3C VSP technologies.1 The attenuation of seismic reflection signals caused by

near-surface and shallow layers only affects the downgoingwavefield (one way in propagation through the media) withappropriate placement of the geophone array at depth.With deployment of the geophone borehole seismic arraybelow 1200 m, the reflection signals recorded at the seis-mic sensors will not propagate through the near-surfaceformations. In this way we could preserve more of the high-frequency components in the recorded data.

2 Multi-component data are highly desirable for detailedreservoir analysis. In a VSP survey when 3C geophonesare used for recording we gain the benefit of acquiring theentire wavefield. This enables independent interpretation ofthe P- and S-wave sections, which can highlight impedancedifferences inherent in the different wave modes, gas cloudvisualization, bright spot confirmation, etc. Secondly, theattributes calculated from these wave modes can be inter-preted stand-alone and/or simultaneously analysed for fur-ther analysis in detail.

3 Accurate depth information from the downhole deploy-ment of receivers in a VSP is very useful in accurate seismicvelocity determination. This gives us a big advantage in ty-ing the surface seismic data to known geology, when com-pared to surface seismic measurements on their own. Theseaccurate velocity measurements also lay a good foundationfor the VSP data processing, particularly for the imaging ofthe multi-component shear-wave mode data.

4 First breaks from different offsets and azimuths can be di-rectly used for anisotropy analysis. In a 3D VSP configu-ration receivers in the downhole array record the signalsfrom many offsets and azimuths. Under the assumptionof anisotropic theory, such as vertical transverse isotropy(VTI) media, anisotropic parameters can be estimated bycombining receiver depths, offsets and traveltime mea-surements of the first arrivals. Geophysicists generally ap-ply anisotropy analysis in two regards. First is that theanisotropic parameters can be used to improve imagingaccuracy and clarity. Second is to use the anisotropy mea-surements in reservoir analysis to determine rock propertiesthat are related to fracture density and pattern such as hor-izontal, vertical fractures and their orientations.In late 2004 the China National Petroleum Corporation

(CNPC) completed feasibility studies of new explorationmethods. At this point CNPC made the decision to spon-sor a multi-component exploration project in this area, whichincluded a surface 3D-3C seismic survey and a 3D-3C VSPsurvey. In this paper, we only discuss the application of 3D-3C VSP technologies. Most of our work centred on acquisitionplanning, target-oriented data processing, data interpretation,multi-component attributes analysis and integrated reservoirevaluation.

Geometry design and planning optimization

Per the design specification for the multi-component seis-mic exploration project, the 3D-3C VSP survey was not astandalone project and would be combined with a 3D-3C

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

3D vertical seismic profile multi-component data 143

surface seismic data acquisition project. Therefore in the 3D-3C VSP geometry design, the main shooting parameters in-cluding shot numbers, shot positions, hole depths and explo-sive sizes are shared with those in the surface seismic survey.Only the borehole seismic array parameters including obser-vation depth interval and receiver spacing were yet to be de-termined. Thus, the VSP geometry design will only focus onthe configuration of the downhole receiver array. The geom-etry design, independent of shooting parameters was focusedtowards acquisition of data with higher resolution (frequency)and a higher S/N ratio. The geometry design is divided intothree steps.1 Determination of the optimal depth range for the borehole

seismic receiver array. As shown in (Fig. 3), most of thehigh-frequency components and amplitude of the seismicwaves are attenuated at the depth interval between 0–1200m. This means that one useful method to reduce the atten-uation in the data acquisition is to put the receiver arraybelow the depth of 1200 m. In this way, the recorded re-flection signals from the target zone consist of two parts.The first is the downgoing wavefield from the source to thetarget zone, which will be attenuated by the shallow layers.The second part is the upgoing wavefield from the targetzone to the receivers, this portion of the travelpath will befree of the attenuation caused by the shallow layers withthe above configuration. Thus, more high-frequency com-ponents can be recorded in the raw data, which will resultin higher resolution imaging. In addition, the receiver arrayis generally needed to deploy to the target depth as closeas possible to better record the reflection signals from thetarget zones. However, for a given offset in a VSP survey,the relationship between receiver depth and the imagingcoverage looks like focusing in a camera, the shallower thereceiver depth, the bigger the imaging coverage. If the re-ceiver array is deployed in a deep depth far from the surface,the imaging coverage will be greatly reduced. Thereby weneed to optimize the receiver depth range between 1200 mand the target depth through forward modelling to strike abalance between the maximum imaging coverage and betterrecording.

2 Initial geometry designs. The longest downhole receiver ar-ray we had at this time consisted of an 8-level array of3-component geophones. This was a highly limiting fac-tor in the acquisition design but allowed for many possibledeployment depths at or above the target zones. Based onthe geologic objectives and the constraint that the obser-vation depth should be greater than 1200 m, several 3D-3C VSP survey geometries were modelled by incorporating

Table 1 Acquisition parameters designed for 3D-3C VSP survey

Acquisition parameters Contents

Total of shots 12730 (explosives)Shooting line spacing 280 mShooting point spacing 40 mReceiving interval 1500–1640 mReceiver spacing 20 mReceiver array 8-level 3C geophones (MaxiWave)

well data, near-surface structure, velocity, forward mod-elling (ray trace and wave equation) and charge size testshooting. As a result of this acquisition modelling manyparameters, including imaging coverage, reflection fold, in-cident angles, illumination maps and response of the shearwave on the target zone were calculated within this process.

3 Planning optimization. Planning optimization is based onthe initial geometry designs. In fact, planning optimizationis a procedure to find an optimal geometry among the initialgeometry designs that are technically and economically fea-sible. All of the calculated parameters are used for compari-son and evaluation. Taking the receiver depths for example,in general, the shallower the depths, the bigger the cover-age of imaging and the lower the reflection fold is at thetarget zone, sometimes we need to make tradeoffs amongthem in this planning. The costs of field operations shouldbe considered, in this case as it was a simultaneous acqui-sition with the surface 3D-3C survey we had little impacton the acquisition cost. Another important step in planningoptimization is to analyse the field test data to optimizethe hole depths, explosive sizes and possible downhole re-ceiver conditions. On the basis of this evaluation, the finalacquisition geometry can be implemented.Thus, the final design parameters for the 3D-3C VSP survey

are described in Table 1.

Target-oriented data processing

Design and optimization of the 3D-3C VSP geometry was justthe first step to realize our aims. After the 3D-3C data acquisi-tion, turning the multi-component data into a useful productfor the reservoir analysis was the major consideration thatwe needed to address in data processing. We performed thefollowing analysis to try and outline the crucial steps neededin data processing, some of which will be introduced here indetail.

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

144 Yousheng Yan et al.

Table 2 Acquisition parameters in the 3D-3C VSP survey

Acquisition parameters Contents

Total of shots 15 294 (12 730) / explosivesShooting line spacing 280 (280) m; 140 m (the area

marked by a white dashedrectangle)

Shooting point spacing 40 (40) m; 200 m (the area markedby a white dashed rectangle)

Receiving interval 1500–1640 (1500–1640) mReceiver spacing 20 (20) mReceiver array 6 (8)-level 3C geophones

(GeoWave R©)

The parameters in parenthesis are from the geometry design listed in Table 1.

Based on the acquisition design parameters, there is littledoubt that the 3D-3C VSP survey is massive, with more than12 730 shots and shot coverage of 173 km2. Because thenumber and location of the shots are from the design of thesurface seismic survey, some of the shots that are far fromthe wellhead are useless to the VSP image volume. It is well-known that the larger offset shot locations in a VSP config-uration will probably only generate refracted events, insteadof useful reflection data. The maximum source to the surfacewellhead offset that was recorded in the 3D VSP survey wasapproximately 10 km. Because of this some of the longer off-set shot data were not included in the VSP data processingsequence. Forward modelling was used to make the determi-nation of the useful shot to wellhead distance to include inprocessing. Finally, the 3C raw data within an offset rangeof 0–8000 m (from wellhead to sources) were selected forprocessing.

Application of multi-component information to gas reser-voir analysis is an innovative and systematic process. First,the full wavefield is recorded by 3C geophones. This leads toproblems of effective wavefield separation of the primary re-flection wavefields, mainly upgoing PP- and PSv-waves from3C data. Secondly, how to transform these signals to PP andPSv volume images that can be easily integrated for interpre-tation? Finally, how to integrate the two imaging volumes inreservoir analysis? What kind of multi-component attributesare sensitive to the tight gas sand and how to extract them?In fact, the answers to the above questions imply that wave-field separation, velocity model building, imaging methods,multi-component attribute calculation and application in dataprocessing are the key to success. Each key point, however,presents unconventional problems that need to be solved inthe implementation. For the target-oriented data processing,we designed a detailed full 3C data processing flow, which

was linked by these key points. In the following text, we willdiscuss some of them in detail.

Wavefield separation

In contrast to surface seismic data, the wavefields recordedin typical VSP data sets are much more complicated becausethe receivers are located downhole instead of on surface. Ba-sically, for a P-wave source, the major wave modes presentin the VSP data include the following five wave modes, theyare downgoing P-wave, downgoing converted shear wave, up-going P-wave (PP), upgoing converted shear wave (PSv) andeven downgoing shear wave (SS) (because the source is not aperfect P-wave source in practice). The wavefield separationdiscussed in this paper is under the assumption of wave prop-agation in isotopic media. In most cases, the upgoing waves,which propagate in the same direction as its polarization ori-entation (P-wave), are commonly used for imaging. But forour project, we have to image the upgoing converted waves(PSv) that propagate in the direction perpendicular to its po-larization orientation, as well as the upgoing P-wave. Theupgoing PP- and PSv-waves are often not only mixed witheach other but are also heavily contaminated by the downgo-ing wavefield. Because the different wavefields are distributedin the 3C data, the complexities of wavefield separation areobvious.

The main target for wavefield separation of 3D-3C VSPdata is to extract the primary reflection signals for imaging,those are, upgoing PP- and PSv-waves. A common wavefieldseparation for 3D-3C VSP data is usually implemented in twosteps. First, the two primary reflection signals are isolated bypolarization reorientation of the multi-component data. Sec-ond, filtering methods are applied to remove other residualwaves from the two primary reflection signals, respectively.The major difficulties are at the first step. As we know, thehodogram analysis is a routine method for polarization re-orientation and widely used in VSP data processing. But itcannot be fully applied to 3D-3C VSP data. Hodogram analy-sis is based on the amplitudes in the short time windows givenby the first breaks of 3C data. Its basic principle is to use theamplitude cross-plot of first breaks in an orthogonal coordi-nate system to find the wavefield’s polarization angle betweentwo component data. The problem is that the estimated an-gle is a linear fitting result based on the cross-plot. In mostcases, this method is generally used in wavefield separationof offset VSP data. Due to a constant offset in this geometry,the variations of polarization orientation for different waves,particularly for the primary reflection signals, are limited to

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

3D vertical seismic profile multi-component data 145

a small range. So, some of the wavefields in offset VSP datacan be reoriented with these angles correctly. But for a 3D-3CVSP or with walkaway VSP data, their offsets change fromnear to far. At this point, it is impossible for us to reorientthe data very well because the angles estimated by first breaksare of great deviation with the polarization orientation of pri-mary reflection signals. Another problem occurs when usingmultichannel filtering on the data. This occurs because weare equipment constrained in the acquisition of this 3D VSPsurvey, we used a small downhole receiver array with veryfew levels. Thus, most of the routine filtering methods do notwork on the common shot gathered data because we have fewtraces in the transform. This means that we have to performthe filtering in a common receiver gather domain in which thenon-linear distribution of wavefields increases the difficultiesin wavefield separation. Even though the shot static correc-tion is also a problem in the filtering in a common receivergather domain, the shot statics derived from surface seismicdata acquisition and processing can be used to improve thissituation.

In order to improve the effectiveness of wavefield separa-tion, we developed a new approach called time-variant vec-tor analysis. The method is based on the calculation of thepropagation direction and polarization orientation for variouswaves at different traveltimes. This means the calculations arein the time-spatial domain instead of the frequency domain(f -x or f -k domain) and other domains. Obviously, the timeand spatial variant propagation and polarization informationcalculated in this way are more accurate and helpful to isolatethe various wavefields recorded in the time-spatial domaineffectively. We implement this method in three steps. First,the raw 3C data (Z, H1, H2) are transformed to (Z, R, T),where Z is a vertical component, H1 and H2 are commonlycalled horizontal components and R and T are respectivelythe radial and transverse components. Second, the verticaland radial components (Z, R) are generally used for furtherprocessing because most of the upgoing PP- and PSv-wavesare in these two components. The calculations of the propa-gation direction and polarization orientation versus time forthe major wavefields are done with forward modelling. Withthe assumption of linear polarization of the wavefields, the up-going PP- and PSv-waves in the Z and radial components arereoriented on their own polarization planes in terms of or-thogonalities of displacement orientation for compressionaland shear waves. The key point in this method is to set uptwo polarization planes, one for the upgoing P-wave and an-other for the upgoing converted wave. In the procedure oftime-variant vector analysis, the other wavefields in the two

components, such as downgoing waves will be residual (mi-nor constituents) in the two planes. We used a high resolutionRadon transform to remove the other residual wavefields onthese two data subsets, respectively.

Velocity model building

Can we simply apply velocity analysis methods used in sur-face seismic data processing to 3D VSP data processing? Theanswer is negative because it is impossible to sort the raw VSPdata into common-midpoint (CMP) gathers in the same wayas we do with surface seismic data. The fact that a downholereceiving array for the 3D VSP survey does not move or movesvertically leads to this impossibility. Therefore, utilization oftraveltimes and receiver depth information in VSP data to es-timate velocity is naturally a well-known means. But velocityestimation for us is not a routine procedure. The procedureof velocity model building is related tightly with complexitiesof surface and subsurface conditions. We try to simplify theprocedure as much as possible. In our method, inversion algo-rithms iteratively carry out the velocity estimation. We buildthe velocity model in one of three ways.1 One-dimensional velocity model building. The essence of

traveltime calculation in the 3D case can be expressed asfollows

T =∫

s (x, y, z) dl, (1)

where T is traveltime, x and y are horizontal coordinates,z is the vertical coordinate, dl is the differential distancealong the seismic ray and s(x, y, z) is the slowness(reciprocalvelocity) at the point (x, y, z). For a 1D velocity model,we assume that the subsurface consists of horizontal layerswith homogenous media. The slowness in the subsurfaceis only related to z. The corresponding inversion problem(linear inversion assumption) is

�T = D�s, (2)

where �T is a vector whose components are the differencebetween the traveltime calculated for the model and theobserved traveltime, �s is a vector whose components arethe differences in slowness between the model and the truesolution and D is a Jacobian matrix whose elements arethe derivatives of T with respect to slowness. The aboveequations are just general forms. The inversion of velocityis carried out by solving equation (2) iteratively. Some cor-rections and constraints are needed to get more accurateresults.

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

146 Yousheng Yan et al.

For estimation of velocity in a one-dimensional model, thefirst breaks of P- and S-waves (if any) on zero-offset VSPdata are picked and used for the calculation of intervalvelocity by using equation (2). The velocity model built inthis method is accurate vertically near the observation well.For a small imaging radius (several square km) and sim-ple subsurface conditions, this type of model is sometimessufficient.

2 Anisotropy analysis. The anisotropy analysis we performedin this paper was primarily used for he evaluation of itsinfluences on imaging. The analysis begins with the com-parison of traveltimes. First, we calculate the differencesbetween the first breaks computed for a 1D velocity modeland the observed first breaks in the same VSP geometry.For an area whose geologic structures are simple and for-mations are relatively flat, the small differences indicatethat the media are approximately isotropic and big differ-ences remind us to pay more attention to the anisotropy insubsurface formations. Further consideration to anisotropyis to estimate the anisotropic parameters. Thomsen’s theory(1986) demonstrated the relation among traveltime, phasevelocity and anisotropic parameters such as ε, δ and γ un-der assumption of VTI or horizontal transverse isotropy(HTI) media. With the same basic principles as described inequations (1) and (2), we can estimate the anisotropic pa-rameters with an inversion algorithm. The main differencesin the inversion procedure are that the slowness expressionin equation (1) and the elements in the Jacobian matrixexpressed in equation (2) are different. For estimation ofanisotropic parameters, the elements in matrix D are com-posed of the derivatives of T with respect to anisotropicparameters instead of slowness. In this case, we can use1D velocity and anisotropic parameters to construct thevelocity model.

3 Complicated velocity model building. For a complicatedsubsurface condition, velocity model building faces greatchallenges. Our strategy is to take full advantage of trav-eltime information in the VSP data. Not only first breaksbut also reflection traveltimes are needed. An inversion al-gorithm also performs this velocity estimation. Of coursetomography is a popular method. Usually, we apply oneof two methods to build the velocity model. First is anadaptive tomography algorithm that does not need to builda geologic model before inversion. The other is a type ofmodel-driven method that needs to input a guess at the ve-locity model either derived from the surface seismic velocityestimation or a geologic model provided by the interpreters.Shear-wave velocity model building is much more compli-

cated than that used for P-waves and should be done aftercompletion of P-wave velocity model building (Yan et al.

2004). In general, the traveltime data of the upgoing PSv-wave on the major horizons shall be picked. The calculationof the Jacobian matrix in the inversion is just related to theraypath of the upgoing PSv-wave because the downgoingwave propagates in the P-wave velocity model. A very im-portant point we shall pay more attention to in the inversionis the ambiguity between depth and velocity. In this case,taking advantage of the depth information in the VSP datacan reduce the uncertainty.

Images of PP- and PSv-waves

For integration of the multi-component information in reser-voir analysis, we performed the imaging of the PP- and PSv-wave data in the depth domain, instead of in the time domain.Its advantages are that there are no scaling problems on thematching of PP- and PSv-images in the depth domain, if thevelocity models are reasonable and accurate. Even thoughVSP common-depth-point (CDP) is a general method for VSPdata imaging, Kirchhoff or wave equation migration is recom-mended. For the areas where the formations display stronganisotropy, anisotropic migration should be considered. Inthis paper, we performed the depth migration on 3D VSP datawith a 1D velocity model and VTI anisotropic parameters.

Calculation of centroid frequency ratio

In terms of seismic theory, different media such as sand andfluid influence the propagations of seismic waves differently.In this area, (Fig. 4) shows the response of compressionaland shear waves to dry sands, shale’s and gas-bearing sands.It can be seen that the shear-wave velocity is not sensitiveto the gas in sands and the compressional wave velocity issignificantly influenced by the presence of gas in the sand.Therefore, the ratio of the shear-wave attribute to the com-pressional wave attribute gives prominence to the response ofgas in sands. Usually P-wave amplitudes (or impedance) arecombined with S-wave amplitudes (or impedance) to identifytight gas sands. However, how to preserve their true (abso-lute or relative) amplitudes in data acquisition and processingand how to calculate their impedances (particularly for shear-wave impedance) by inversion are very tough. So, we made anattempt to use an attribute related to frequency for analysis.

The method of centroid frequency calculation was proposedby Quan and Harris (1997) and used to estimate the at-tenuation on crosswell seismic data (Yan et al. 2001). We

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

3D vertical seismic profile multi-component data 147

used the centroid frequency of various wavefields to gen-erate a new attribute for the delineation of the gas-bearingsands.

The centroid frequency calculation is given by

centroid frequency =∫

R ( f ) f d f∫R ( f ) df

, (3)

where R(f ) stands for the amplitude spectrum of seismic dataand f for frequency. The centroid frequency is calculated in thefrequency domain. Formula (3) indicates that a centroid fre-quency value is independent of amplitudes because the numer-ator is normalized by the integration of the amplitude spec-trum along the frequency. This means its values are primarilydominated by the variation of frequency and less influencedby the amplitude. This character can simplify the complexitiesof amplitude consistency and preservation in data processingand imaging. The previous work (Quan and Harris 1997)demonstrated that the calculation of the centroid frequencyis stable and its result is more reliable than amplitude in theattenuation analysis.

For multi-component attributes the shear modulus is theo-retically independent of a fluid saturation level. This propertyof shear modulus can be derived from other shear-wave at-tributes such as velocity and impedance. In our project, wepropose a new attribute titled the centroid frequency ratioand try to apply it to indicate the presence of gas-bearingsands. The centroid frequency ratio is defined as follows

RCF (x, y, z) = centroid frequencyPSv (x, y, z)centroid frequencyPP (x, y, z)

, (4)

where centroid frequencyPP is the centroid frequency of thePP-wave, centroid frequency PSv is the centroid frequency ofthe PSv-wave and RCF is their ratio. As the saturation levelof gas in the sands increases, the attenuation of frequencycomponents for compressional and shear waves is different.

The centroid frequencyPP becomes smaller because the high-frequency components are greatly attenuated and the centroidfrequencyPSv changes little. This results in an increase in theRCF ratio, we believe this response is indicative of the presenceof a tight gas sand anomaly.

R E S U L T S

The 3D-3C VSP data acquisition was conducted simultane-ously with a 3D-3C surface seismic survey in 2005. At thattime, we designed the VSP geometry with an 8-level (maxi-mum) 3C downhole geophone array. But we were told in thefield test that 2-levels of the geophone array did not work. Tocompensate for the acquisition with the smaller receiving ar-ray, about 2500 shots were added to the data acquisition. Thenumber of the additional shots was determined by forwardmodelling. The basic principle is that the additional shots canapproximately compensate the reduction of reflection folds onthe target zone by using 6-level geophones instead of 8-levelgeophones.

Thus, the final acquisition parameters for this project werechanged as follows.

Figure 5 shows the map of the shots where dots stand for lo-cations of the shots, the colours indicate the surface elevation,red colour stands for high altitude and blue for low, the redcircle is the location of the deployed well and the red line tiedwith the well is about 10 km long. Note that the added shots(about 2500) are located in the area marked by a white dashedrectangle. The maximum difference in elevation is greater than200 m. To aid the data processing and interpretation in thisproject, two zero-offset VSP surveys were conducted in thesame well with P- and S-sources, respectively.

The data processing started with an analysis of raw data.Common receiver gathers instead of common source gathers

Figure 5 Map of shots in the 3D VSP sur-vey. The white dashed line indicates the areawhere 2500 shots were added. The red lineis about 10 km long and the red circle is thewell where the receiver array is deployed.The colours represent the surface elevation.

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

148 Yousheng Yan et al.

were chosen for analysis due to only having data six channelsin the latter. The two zero-offset VSP data sets were analysedand processed at the same time. By doing so, both vertical andlateral propagation variations of the seismic waves were out-lined. The vertical attenuation is mainly attributed to the shal-low layers as shown in (Fig. 3). The contours of first breaks ina common receiver gather are composed of concentric circles,which mean that the velocity in this area is without stronglateral variation. But at an offset of 7000 m the difference

Figure 6 Interval velocity calculated from zero VSP and VTIanisotropic parameters estimated from 3D VSP data. An azimuthwalkaway data in the 3D VSP data set is used for estimation of theanisotropic parameters.

between the first breaks of the P-wave arrival calculated for a1D velocity model and the observed first breaks of the P-waverecorded at the depth of 1600 m is up to 160 ms, which isto say the anisotropy from the overburden sediments has tobe considered in the imaging because the type of traps in thisarea is generally stratigraphy instead of structure.

It took about two months to complete the data processing.In the first phase, a series of data processing steps includingsurface consistent amplitude correction, shot static correction,attenuation compensation and 3C polarization reorientationwere implemented. Subsequently, using the time-variant vec-tor analysis and high resolution Radon transform we effec-tively separated the upgoing PP- and PSv-waves. Then, wefound that the overburden formations are deposited stably,well stratified and flat but had a highly VTI anisotropic char-acter. At this point anisotropic parameters (ε, δ) were esti-mated by an inversion algorithm and integrated with the in-terval velocity derived from the zero-offset VSP data to buildthe velocity model for imaging. Figure 6 shows the isotropicinterval velocity and VTI anisotropic parameters (ε, δ) abovea depth of 1600 m (the deepest receiver in the 3D VSP sur-vey). Finally we applied Kirchhoff migration on this wavefieldseparated data using the above derived velocity models. Wethen produced two image volumes (PP- and PSv-waves), witha bin size of 20∗20 m and depth sampling of 2 m. The imagingareas of PP- and PSv-volumes at the target depth cover around18 km2 and 6 km2, respectively. Figure 7 shows the two imag-ing volumes (left: PP-wave imaging; right: PSv-imaging) wherethe PP-wave imaging volume was muted to the same coverageas PSv-imaging.

Figure 7 3D VSP PP (left) and PSv (right) images. The green dash polygrons highlight the target area.

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

3D vertical seismic profile multi-component data 149

Figure 8 Comparison of VSP (right) and surface seismic (left) images. Top: PP-wave. Bottom: PSv-wave. The red rectangles highlight the targetarea.

Figure 9 Comparison of the amplitude spectrum of surface seismic imaging and VSP imaging. The black line stands for surface seismic data andthe red line for VSP data. Left: PP-wave. Right: PSv-wave.

In order to compare VSP data with surface seismic data,the imaging volumes were converted from the depth domainback into the time domain. The green polygons in the figureindicate the zones of interest at a time interval of1650–1900 ms. The bottom line of the polygons is locatedin the trough of a strong reflection event that corresponds toa coal bed. Obviously, both of the images display good sig-nal quality. Figure 8 shows the comparisons of surface seismicimages and the 3D VSP images. The top display is for compari-son of P-wave images and the bottom one is for comparison ofPSv-wave images, with the target zones marked by red rectan-gles. Two points can be summarized in the comparison. One isthat there is a high degree of similarity between these images,

which demonstrates the reliability in imaging. The other isthat the VSP images are of higher resolution. Quantitatively,the calculations of amplitude spectra on these images indi-cate that the dominant frequency of VSP images is 10–15 Hzhigher than that of surface seismic images as shown in(Fig. 9) where black and red lines represent amplitude spectraof surface seismic data and VSP data, respectively.

The two image volumes were interpreted independently.The main target horizons in this area are in the Permianand include the TP and TC formations, which are in a depthrange of 3100–3500 m (a time range of 1650–1900 ms). TheTP7 and TP8 formations are the major pay zones. Data in-terpretation started with horizon identification via tying the

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

150 Yousheng Yan et al.

Figure 10 Target horizons in the 3D VSP imaging are identified byzero-offset VSP data. The red arrows indicate the major interestinghorizons such as TP7, TP8, TP9 and TC3.

Figure 11 Time difference between TP8 and TP7 horizons is basedon PP imaging of 3D VSP data. The red rectangle indicates the areaof PSv imaging.

zero-offset VSP data integrated with log and drilling dataand was then used to identify several horizons of interest.Figure 10 shows the horizon identification for the PP-imaging.We performed the same interpretation on the PSv-imaging byusing the zero-offset VSP data acquired with a shear-wavesource. Because most of the reservoirs in this area are dom-inated by stratigraphy and not by structures, the interpretedhorizons are mainly used to indicate spatial configuration ofthe target zones and additionally used to calculate the thick-ness between two horizons such as sand and shale formations.Figure 11 shows the map of the PP-imaging time differencesbetween the TP8 and TP7 horizons, where the red rectangleindicates the coverage area of PSv-imaging. The colours in thefigure stand for the thickness in time, ranging from 15–30 ms.

The centroid frequencies of the PP- and PSv-images werecalculated with formula (3). For calculation of the centroidfrequency ratio, the PP-imaging volume was muted to thesame coverage as the PSv-imaging volume. Thus, two cen-troid frequency volumes were generated over a coincidentarea. Simply, the centroid frequency volume of the PSv-wavewas divided by the centroid frequency volume of the PP-waveto generate the centroid frequency ratio volume. Figure 12shows the centroid frequency ratio images in a chair (left)and a horizontal slice (right), at a depth of 3340 m (1810 msin time), where the red colour represents large values for thecentroid frequency ratio and the blue colour represents lowcentroid frequency ratio values. It can be found that most ofthe area is filled with a blue colour representing a low cen-troid frequency ratio, with two high centroid frequency ratioanomaly zones, indicated by the red colour, which are obvi-ous in the north-south direction. Drilling results in this areaindicate that the TP8 formation lies in the depth interval be-tween 3287–3345 m. This places the depth of the anomalyzone near the bottom of the TP8 horizon.

The local geologic setting tells us that the source rocks arelocated in the coal beds and black shale of the Taiyuan andShanxi formations of Permian age. The TP8 formation in thisarea is generally a sand rich formation, deposited in a largedeltaic system that is 10–20 km wide in an east-west extentand about 200 km long in the north-south direction. Thechannel sands are primarily distributed in the north-southdirection. Figures 11 and 12 show the consistency of the seis-mic response within this geologic constraint. That is, the twoanomaly zones correspond to the prospect geometry expectedwithin this geologic setting. Therefore, the interpretation ofhorizons on the 3D VSP data was used to indicate where thesands are and their relative thickness. These calculations andanalysis of multi-component attributes were used to delineatethe potential reservoir targets. For the interpretation of thecentroid frequency ratio presented in this paper, we believethe centroid frequencyPP values reflect the frequency attenu-ation from sands and gas (if any) and centroid frequencyPSv

is only affected by the presence of sand, their ratio minimizesthe frequency variations caused by sands and amplifies thefrequency variations attributed to the presence of gas. Manyresults from conventional seismic data also demonstrated thatthe dominant frequencies of PP-wave data are strongly attenu-ated by gas in the formations. Based on the analysis describedabove, the two anomaly zones shown in (Fig. 12) were in-terpreted as gas-bearing zones based on the integration oflogs, drilling and geologic data. In addition, a series of con-ventional methods for attribute estimation were applied. All

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

3D vertical seismic profile multi-component data 151

Figure 12 Centroid frequency ratio volume (left) and a horizon slice of the centroid frequency ratio (right). Two anomaly zones delineated bythe centroid frequency ratio were interpreted as gas-bearing sands.

of them were used in the final comprehensive evaluation. Asa result, one well location for drilling was proposed in thenorth-west portion of (Fig. 12). A positive feedback on thisproject was received.

CONCLUSIONS

This paper indicated the feasibility and effectiveness of theapplication of 3D-3C VSP data to delineation of the tightgas sands in the Sulige gasfield. The concept of using target-oriented 3D-3C VSP exploration in reservoir developmentdrove the implementation process that included data acqui-sition, processing and interpretation. The major conclusionspresented in this paper are as follows.1 Optimizing geometry design greatly improves the qual-

ity of raw data recorded in the 3D-3C VSP survey. Thedownhole receiver array was configured below a depth of1200 m for reduction of the attenuation from shallow lay-ers and did record the upgoing reflected waves with thehigher-frequency components (only one-way attenuation)preserved. The fact that the resolution of VSP reflectionimages is higher than that of surface seismic images makesa more detailed data interpretation possible.

2 Target-oriented data processing is essential in generatinghigh quality PP- and PSv-images, thus enabling the applica-tion of multi-component information to reservoir analysis.Accurate velocity model building and effective wavefieldseparation methodologies are the most important aspectsof target-oriented data processing. High precision images

of PP- and PSv-waves in the depth domain avoid the scal-ing problems encountered with time-domain volumes. Theintegration of PP- and PSv-image volumes in data interpre-tation provides more information for reservoir analysis andreduces uncertainty.

3 A new attribute called the centroid frequency ratio of PP-and PSv-waves is proposed in this paper. This attributecalculation minimizes the influences from a possible incon-sistency of seismic amplitudes and makes the results morereliable. The attribute combines multi-component informa-tion and is expressed by a ratio that can target on the seis-mic responses from pore fluids. So, the new attribute canbe used as a direct indicator for gas-bearing sands.In conclusion, we believe 3D-3C VSP technology can be

successfully used for the delineation of tight gas sands in thisarea. This does not mean the methods presented in this paperwill work everywhere. For different objectives, what we needto do is to investigate the problems, come up with a strategy,perform detailed planning, improve and develop the methodsand see if we arrive at the expected solution.

ACKNOWLEDGEMENTS

We would like to thank BGP for permission to publish thispaper and also Mr Yabin Guo, who provided some figures ondata interpretation. We would also like to thank Marc Sterlingfrom HiPoint Reservoir Imaging, LLC and APEX HiPointLLC in Littleton, Colorado, USA for his help in editing thispaper for final publication.

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152

152 Yousheng Yan et al.

REFERENCES

Gomes F. and Ronholt G. 2005. Multi-well 3D VSP PP and PS imag-ing used for structural interpretation in the onshore, CAM-field,Brazil. 75th SEG meeting, Houston, Texas, USA, Expanded Ab-stracts, 2645–2648.

Quan Y. and Harris J. 1997. Seismic attenuation tomography usingthe frequency shift method. Geophysics 62, 895–905.

Thomsen L. 1986. Weak elastic anisotropy. Geophysics 51,1954–1966.

Yan Y., Xu Z., Yi M. and Wei X. 2007. 3D VSP PP and PSv

imaging for carbonate reservoirs. 69th EAGE meeting, London,UK, Expanded Abstracts, H016.

Yan Y., Xu Z., Yi M. and Wei X. 2009. Gas reservoir-oriented3D-3C VSP data processing. 71st EAGE meeting, Amsterdam, theNetherlands, Expanded Abstracts, T036.

Yan Y., Yi M., Wei X. and Wan W. 2001. Joint tomographic imag-ing for cross-hole seismic velocity and Q value. Oil GeophysicalProspecting 36, 9–17.

Yan Y., Yi M., Wei X. and Xu Z. 2004. C wave processing of 3DVSP data. 74th SEG meeting, Denver, Colorado, USA, ExpandedAbstracts, 2509–2512.

C© 2011 BGP/China national Petroleum Corporation, Geophysical Prospecting, 60, 138–152