- Email: [email protected]
Seismic imaging by 3D partial CDS method in complex media
Author’s Accepted Manuscript Seismic imaging by 3D partial CDS method in complex media M. Soleimani
www.elsevier.com/locate/petrol
PII: DOI: Reference:
S0920-4105(16)30062-6 http://dx.doi.org/10.1016/j.petrol.2016.02.019 PETROL3362
To appear in: Journal of Petroleum Science and Engineering Received date: 19 April 2015 Revised date: 2 December 2015 Accepted date: 22 February 2016 Cite this article as: M. Soleimani, Seismic imaging by 3D partial CDS method in complex media, Journal of Petroleum Science and Engineering, http://dx.doi.org/10.1016/j.petrol.2016.02.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Seismic imaging by 3D partial CDS method in complex media M. Soleimani* Faculty of Mining, Petroleum and Geophysics, University of Shahrood, Shahrood, Iran.
[email protected] *
Tel-Fax: +98 (23) 3239 5509
Abstract
Seismic imaging in complex geological structures is state of art to deliver a high quality section for further geological interpretation. Conventional processing and imaging methods will fail in case of facing complex media, strong heterogeneity and poor quality data. Thus, conventional seismic imaging methods have to be improved inevitably or new approaches have to be extended to produce valid images. The partial common diffraction surface method is based on improving the common diffraction surface operator to handle lateral heterogeneities in complex media. In this study, the partial diffraction surface method was improved to 3D. This method enhances quality of the 3D pre-stack data to be used for further depth imaging. It can also resolve some of the possible ambiguities in geological interpretation from seismic images. This also will enhance the quality of the seismic images, where it suffers from conflicting dips problem. A time variant linear function defines the offset range for each zero offset sample in the stacking operator. This function controls the offset increment for the partial form of the operator in the time-midpoint-offset domain. It should be noted that the velocity model can be considered as smooth and simple as possible in this method for further depth imaging. The introduced method was applied on a 3D synthetic and a real land data set. Results from both data example show the ability of the new method in enhancing the quality of the seismic section in the presence of faults and lateral velocity heterogeneity. This enhancement achieved by a simple velocity model in a strong complex media. Consequently, depth migrated result of the real land data underwent a time conversion to create a time underground contour map, for comparison with the maps obtained with conventional methods. More details of the target formation were observed on the map used result of the proposed method for mapping.
Keywords: Seismic imaging, common diffraction surface, lateral velocity heterogeneity, Complex media, Zagros overthrust. 1.
Introduction
Seismic imaging methods that are based on ray tracing techniques face with problem in imaging the complex media data. Geologically boundaries of complex structures such as salt dome, overthrust zone and faulted - folding system are illuminated by different ray paths. Tracing the paths of these rays cannot be properly handled by conventional imaging techniques (Yilmaz, 2001). Significant improvement in quality of imaging seismic data from such media could be achieved both in performing an accurate velocity model building procedure and/or applying advanced migration algorithm. In imaging complex structures with strong heterogeneity, the conventional migration methods based on the Kirchhoff integral, have
1
theoretical and practical shortcomings. Thus combination of accurate velocity model building and appropriate migration algorithm is the key to successful imaging (Biondi, 2006). Among the problems of imaging in complex structures, the conflicting dips situation has been investigated by different methods. The problem of conflicting dips was also observed in application of the new introduced common reflection surface (CRS) method. Concept of the CRS stack method was developed by Jäger et al. (2001) as a reasonable extension of the traveltime equation to fully exploit the information contained in the pre-stack data. Different applications were derived for the CRS stack method ranging from multiple suppression (Dümmong and Gajewski, 2008), pre-stack data enhancement and regularization (Baykulov et al., 2011) and velocity model building for depth migration (Duveneck, 2004). Dell and Gajewski (2011) had introduced a CRS based workflow for diffraction imaging that was followed by Bakhtiari Rad et al. (2014), which used common reflection surface, based on pre-stack diffraction separation method only to image diffractions rather than reflection events. On the other hand, some of the new introduced data driven imaging techniques mainly aimed to simulate zero offset (ZO) sections from multi-coverage seismic reflection data. Among of these new methods, several velocity independent methods have been introduced. Druzhinin et al. (1999) developed a scpecial prestack depth migration technique which avoids the necessity of constructing a macro velocity model. Gelchinsky et al. (1999) introduced the multi-focusing method for application in multicoverage data. Müller (1999) also formulated the CRS method as a macro velocity model independent technique. Garabito et al. (2001) later on introduced a new parameter search strategy by global optimization. Hertweck et al. (2004) also introduced the workflow of obtaining the kinematic wavefield attributes from the CRS equation. These methods were data driven in the way that they: (a) use a multi-parameter moveout formulas, where the moveout parameters are derived based on the coherency analysis, and (b) do not need explicit knowledge of the macro velocity model (Soleimani et al., 2010).The data drive CRS method provide result with high signal to noise ratio with increased vertical resolution.The CRS method imposes some approximations on the wavefield parameter estimation, which indeed less information of velocity model is required (Prüssmann et al., 2007). However, the CRS stack method sufferes from the problem of conflicting dips. Several strategies were introduced to resolve this problem in the CRS method. Höcht et al. (2009) introduced two techniques, called as target-oriented (TO) scheme and an alternative operator-oriented (OO) scheme to remove this problem in 2D ZO CRS. They have mentioned that a ‘proper’ handling of conflicting dip situations requires to identify the different events in terms of kinematics and waveform. Soleimani et al. (2011) introduced a modification of the CRS method that was called the common diffraction surface (CDS) or data based CDS stack method. Shahsavani et al. (2011) introduced the model based CDS stack method that decreased the large computation time of the CDS stack method. Garabito et al. (2011) introduced the surface operator of the CDS method that was a special case of the data based CDS stack method when only the dominant surface is considered as the operator. Yang et al. (2012) called the phenomenon caused by conflicting dips problem as ‘dip discrimination phenomenon’. This method was called the CRS stack with the output imaging 2
scheme (CRS-OIS). However, Soleimani (2015) showed that application of the 2D CDS method could produce enhanced migrated image of the semi-complex geological structures. 2.
The Conventional 3D CRS imaging
The CRS stack is a multi-parameter stacking method that uses any contributions along any realization of its operator. These contributions are tested by coherence analysis for each ZO sample, and the set of attributes, which yield the highest coherency, are accepted as the parameters of the optimum operator to perform the actual stack (Bergler, 2001). The CRS equation is described by three parameters α, RN and RNIP, as (Müller, 1998): 2
2sinα(xm x0 ) 2t 0cos 2α (xm x0 )2 h2 2 thyperbolic (xm ,h) = t0 + + + v0 v0 RN RNIP
(1)
where h is the half offset, Δx=(xm-x0) is the displacement with respect to the considered common mid-point (CMP) position, t0 corresponds to the ZO two-way traveltime, α is the emergence angle of the ZO ray, RN and RNIP are the radius of curvature of the normal (N) wave and the normal incidence point (NIP) wave, respectively, and finally v0 is the near surface velocity. The CRS stack method assigns merely one optimum stacking operator for each ZO sample to be simulated (Heilmann, 2007). However, in situation where different events in a seismic section intersect with each other and/or themselves, a single stacking operator for that ZO sample would be no longer appropriate. Thus, Mann (2001) proposed to allow small discrete number of stacking operators for a particular ZO sample. Therefore in conflicting dip situations, more than one operator for each ZO sample in ZO stack simulation would be used. Müller (2003) formlulated the 2D CRS operator into 3D. The 3D CRS operator is based on a second order approximation of traveltimes in midpoint and half-offset coordinates, (Cristini et al., 2003). The stacking surface in the time space domain, used to compute the CRS volume, depends on eight parameters that describe the shape and direction of normal ray and wavefront curvatures for Normal (N) and NIP wave. The method does not require the exact knowledge of the velocity model, a priori knowledge of the near surface constant velocity being sufficient, (Marchetti et al., 2006). The 3D CRS processing is done for each sample by summing up all events located on a traveltime surface in the five dimensional pre-stack data hyper volume (Müller, 2003): 2
2t 2 2 T T thyp t0 wm 0 mT Qzyz N N Qxyz m hT Pzyz N NIP Pzyz h v v
(2)
Where thyp describes the 3D hyperbolic traveltime approximation related to the local x, y system, centered at point X0 x0, y0. The parameter t0, is the two-way ZO traveltime at point P0x0, y0, t0, v is the surface velocity, w is the projection of the unit vector of the emerging normal ray at point X0 onto the ground surface. It contains the parameters , i.e. the angles that determine the direction of propagation of the two hypothetical wavefronts, NIP wave and N wave, along the emerging central (normal) ray. Curvature of the NIP and N-waves are defined by the NNIP and the NN matrices (Müller, 2003):
3
1 R N i i ,max 0
0 1 Ri ,min
(3)
The angle parameter , are defined by auxillary parameters P, Q and D as follows: Pxyz Dz ( ) Dy ( ) Dz ( NIP ),
(4a)
cos sin QDxyzz D z ( ) Dy ( ,) Dz (N,),NIP , N sin cos cos D y 0
0 1
cos D z sin
(4b) (5)
sin , cos
, NIP , N
cos eight 0 parameters at X0 are, final DThe y 0 1 are obtained with only eight parameters
NIP N, RNIP-max, RNIP-min, RN-max and RN-min. The two searches (one three-parametric and one fiveparametric search). The simplified flowchart of obtaining these parameters is shown in figure 1. pre-stack data Step I
VNMO,min VNMO,max
Automatic 3D CMP Stack
v
one three-parametric search for VNMO,min ,
CMP Coherency CMP Stack
VNMO,max and v.
Step II
,
3D ZO Stack
N
one five-parametric search for ,
RN,min and
N
RN,min, RN,max ZO Coherency ZO Stack
RN,max
Step III Determination of RNIP parameters
3D CRS Stack Stack using the eight paramenters within the projectet Fresnel Zones using the complete CRS operator
RNIP,min RNIP,max
NIP
From VNMO,min , VNMO,max and v (Step I) and ,(Step II), determination of
RNIP,min, RNIP,max and
NIP
3D CRS Stack
Fig. 1 Flowchart of the search strategy of the 3D CRS stack (Müller, 2003).
3.
Impelemntation of the 3D partial CDS
Soleimani and Mann (2008) introduced the CDS stack method for solving the conflicting dips problem in the 2D CRS stack. In this method, for each sample in time domain, a stacking surface is performed in a range of angles. Therefore, there exist many stacking surfaces for one sample 4
in (t, xm, h) domain that makes a volume of stacking operators and all of them will contribute into stacking for that sample. Thus, the conflicting dips will be treated well by the CDS operators. Solving the problem of conflicting dips in this way will enhance the usually weak diffraction events in the stacked section (Soleimani et al., 2009a). For true diffraction events, the radius of the NIP wavefront and the normal wavefront will coincide, RN = RNIP. Thus the only attribute to be searched for in a fixed emergence angle, is a combination of RN and RNIP that is called RCDS, (Soleimani et al., 2009b): 2
2 2sinα xm x0 + 2t 0cos α xm x0 2 + h 2 t 2 xm h = t 0 + v0 v0 RCDS
(6)
By implicit knowledge of the RCDS, the shape of the operator could be defined. Baykulov and Gajewski (2009) calculate a stacking surface around a specified ZO point defined by its offset and traveltime coordinates in a chosen CMP location and perform the summation of data along that surface (Figure 2a). This strategy is called the partial CRS stack method. Baykulov (2009) had mentioned that the partial CRS stack method needs further investigations on the conflicting dips problem. However, prestack migration of noisy and low quality data produces migrated section of comparably lower quality than the post stack migration of the CRS stack. Baykulov (2009) showed that the CRS traveltime formula, where the dip of the reflector element is incorporated, is used to compute new partially-stacked CRS super-gathers, where each trace is a result of summation of data along the CRS stacking surface. The number and location of traces in the produced super-gathers can be defined. To completely remove the problem of conflicting dips and produce high quality super-gathers for further prestack migration, the idea of the partial CRS stack method was used here to modify the CDS stack equation. The new method differs from the partial CRS stack. This method, that is called the partial CDS, uses this idea for operator volume stack rather than only a surface stack. The partial CDS stack method, calculates collection of stacking surfaces around a specified point. This point is defined by its CMP number, zero offset traveltime, t0, and offset boundary related to that t0 in a zero offset section, and performs the summation of data along that surface. The result of summation is assigned to that sample. The operator equation in the partial CDS is same as the CDS operator with performing limitation in offset range according to one offset banding function. The partial CDS stack operator is shown as yellow surface related to the sample t0 (red point) in figure 2b. To have advantage of partial 2D CDS stack, the conventional 3D CRS stack method was improved here not only to 3D CDS stack, but into a partial form. Hence, the radius of normal wave and NIP wave should coincide. Therefore, eight parameters are reduce to five parameters while: NIP = N, RNIP-max = RN-max and RNIP-min = RN-min. Finally these five parameters are: CDS RCDS-max, and RCDS-min. The partial 3D CDS stack equation would be:
5
2
t
2 hyp
2 2t T T t0 wm 0 mT Qzyz NCDSQxyz m hT Pzyz NCDS Pzyz h v v
(7)
Figure 3 shows the proposed search startegy of the partial 3D CDS equation. By knowing the fact that the CDS stacking parameter varies smoothly for a smooth velocity model, the rays can be calculated on a relatively coarse emergence angle grid. In contrast, the semblance and the stack itself are quite sensitive to variations of the emergence angle rather than the initial stacking velocity. Thus, stacking and semblance value perform on a finer emergence angle grid using interpolated stacking parameters. This shows that result of the CDS stack is more sensitive to emergence angle increment rather than the accuracy of the near surface velocity model.
Fig. 2 a) The partial CRS stack performs summation of data around a specified point on a CMP traveltime curve (magenta line) and assigns result to the same point in a newly generated CRS supergather. b) The partial CDS stacking surface shown with a yellow color coincides locally with the common offset travel time surface, but may be limited in size. (Figure (a) from Baykulov (2009) and figure (b) from Soleimani and Mann (2008)). pre-stack data VNMO,min VNMO,max
Step I Automatic 3D CMP Stack
v Offset banding CMP Coherency CMP Stack
one three-parametric search for VNMO,min ,
VNMO,max and v.
Step III
Step II
3D Partial CDS Stack
3D CDS Stack
Stack using the five paramenters within the projectet Fresnel Zones using the complete CDS operator
one five-parametric search for , CDS
RCDS,min and
RCDS,max
Enhanced prestack data by 3D Partial CDS, (ready for depth imaging)
` 6
Fig. 3 Flowchart of the search strategy of the 3D partial CDS method.
4.
Application of the partial CDS on synthetic data
To verify the improved search strategy of the CDS stack, a 3D synthetic data set was generted for processing. The synthetic model data obtain by forward modeling on a five layered dome shape model. Synthetic data was contaminated by random noise for better evaluation. The results of the CRS stack process were three optimized kinematic wavefield attributes, one coherency section and one optimized stacked section. However, the results of the CDS stack process are just one optimized stacked section and the section that shows the number of used traces in aperture. Therefore only the stacked sections are comparable. Figure 4a shows the velocity model used for synthetic data generation. Figure 4b shows the coherency result of application the CDS method on the synthetic data. Results of the CRS and the CDS method on the data are shown in Figure 4c and 4d. As it could be seen, the stacked section obtained by the CDS method could clearly better image truncation of the reflector, conflicting events and continuity of the reflectors. Result proved that the proposed method could be applied on real data.
(a)
(b)
(c)
(d) Fig. 4- (a) Velocity model used for synthetic data generation by ray tracing. (b) Coherency section of the CDS result. (c) The CRS stacked result and (d) The CDS stacked result on the synthetic data.
5.
Application of the 3D partial CDS on real land data
To see how the new introduced method could overcome some of the problem of imaging of complex structures, it was applied on a complex 3D land data. The related seismic data was 7
obtained on a mountainous topography with complex geology that makes strong velocity heterogeneity in the media. Figure 5a shows location of the study area, in SW of the Zagros overthrust belt. The study structure is an extensively faulted anticline with 74 km long and 6 to 8 km wide, (Figure 5b). The northern 60 km of the southwest flank is defined by a reverse fault with up to 2100 m of throw. Normal faults occur on both flanks and the crest, but are most common on the southwest flank, between the crest and the reverse fault (Figure 5b). Figure 6 shows the stratigraphic column of the study area. The major problematic lithology here is the anhydrite - salt layer called the Gachsaran Formation that is responsible for velocity contrast.
(a) (b) Fig. 5 (a) Location of the Zagros overthrust belt in SW of Iran and location of the study area shown by a red rectangle. (b) Structural map of the target anticline and area of 3D seismic data. Locations of the seismic lines are shown by bold black lines. (Figure (a) from Sherkati et al. (2005))
Fig. 6 Stratigraphic column of the study area that shows thick layer of anhydrate - salt in the area (Sherkati et al., 2004)
8
The target formation here is the Asmari formation, which is the major reservoir formation in SW Iran, with an anticlinal shape beneath its cap rock, the Gachsaran formation. The enhanced pre-stack data obtained by the 3D partial CDS method underwent depth imaging by the Gaussian beam migration (GBM) method, (Alkhaliffah, 1995; Hill, 2001). To see the results of the processing, two in-lines and one cross line of the seismic cube were extracted for visualization. Figure 7a shows result of the conventional prestack depth migration (PSDM) and figure 7b shows depth migrated result of the partially enhanced pre-stack data. As it can be seen, migration of the partially enhanced data, images more steep dips in the section. The anticline is dominated in both sides by faults, which was imaged better by migration of partially enhanced data. Image section in figure 7b also shows more continuity in reflectors below the anhydrite and salt layer. Truncations of layers are also better identified by enhancing data with 3D partial CDS method, which helps for easier structural interpretation. Thick layer of anhydrite is also imaged with more continuity by new method. The next in-line for presenting the result selected from a part of the anticline that is dominated by faults and difficult to define drilling path. Figure 8a shows result of the conventional PSDM and figure 8b shows depth migrated result of partially enhanced prestack data. The other major problem besides imaging those faults is to obtain higher quality image for imaging small scale faults in the left part and inside the anticline. Increasing continuity of the reflectors is obvious in migration result of the partially enhanced data, and two large scale faults on both sides of the anticline are also better imaged in figure 8b. Enhanced parts of the section are highlighted by the black rectangles. Faults inside the anticline, that has important effect on defining drilling path, are better imaged in the result of the new method. As it can be seen, some small faults are only imaged in figure 8b, while in figure 8a, only a blurred picture of small scale faults inside the anticline is obtained. In both sections, the quality of the seismic data appears to be affected, at least in some parts, by the variation of thickness of the Gachsaran (anhydrite and salt) formation. To see better result of application the 3D partial CDS method on data, a cross line section was shown that does not suffers from the problem mentioned above. Figure 9a shows a cross section result of conventional PSDM and figure 9b shows a cross section depth migrated result of partially enhanced pre-stack data. The final migrated result shown in figure 9b has better quality with more continuity in the reflectors, especially in the right down side of the section. The pinch out imaged below the CDP number 1700 in depth of 1600 m is well imaged in figure 9b, while it is not imaged in figure 9a. For better understanding of improvement on results, two depth slices from the seismic cube are also shown in figures 10 and 11. Enhanced parts of the sections are highlighted by the black rectangles. Figures 10a and 11a show depth slices result of conventional PSDM and figures 10b and 11b show the related depth slices result of depth migration of partially enhanced pre-stack data. The quality increasing in slices obtained by the migration of the partially enhanced data is obvious. Tight lines are well imaged in figures 10b and 11b. Almost all of the reflectors on both side of the anticline are imaged by the new method. Thus, it could be concluded that the pre-stack depth migration of the pre-stack data enhanced by 3D partially CDS method can be applied to complex structures to resolve some of the ambiguities of imaging in such media. In all the examples shown, horizons
9
are more continuous in the migration result of partially enhanced data and the increasing signal to noise ratio is obvious.
(a) (b) Fig. 7 Migrated section of in-line number 1800. (a) Result of conventional PSDM with Gaussian beam algorithm and (b) result of depth migration of partially enhanced pre-stack data.
(a) (b) Fig 8 Migrated section of in-line number 3100. (a) Result of conventional PSDM with Gaussian beam algorithm and (b) result of depth migration of partially enhanced pre-stack data.
10
(a) (b) Figure 9- Cross section number 700. (a) PSDM with GBM and (b) migrated result of partially enhanced pre-stack data, migrated by pre-stack Gaussian beam depth migration.
(a) (b) Figure 10- Depth slice of depth 1200m. (a) PSDM (with GBM) and (b) migration of partially enhanced pre-stack data, migrated by pre-stack Gaussian beam depth migration.
(a) (b) Figure 11- Depth slice of depth 1700m. (a) PSDM (with GBM) and (b) depth slices result of depth migration of partially enhanced pre-stack data, migrated by pre-stack Gaussian beam depth migration.
6.
Structural interpretation and velocity model considerations
The Gachsaran formation governs the velocity of the whole swelling part. In the center of the Gachsaran Formation, salt is reported as the main deposit. In areas where the amount of anhydrite is high, the expected total interval velocity would also increase. Therefore, the Gachsaran ridges show anomalous seismic velocity behavior (Sherkati et al., 2006). Usually conventional time-migrated seismic sections are distorted and obscure owing to the presence of inflated Gachsaran bodies (salt thickening) and related lateral velocity variations. In such conditions, strong lateral velocity variations, related to the lithology contrasts between steeply dipping layers bend the seismic rays like an optical lens and distort the sub-surface image. The 11
new method introduced here is considered to be the appropriate method for imaging targets in the presence of heterogeneous overburden, especially in the presence of a salt body. The only way to remove effects of salt in seismic image is to define the velocity variations by building a velocity model and performing depth migration, which compensates for ray-bending propagation effects. The velocity model used here was obtained by the method introduced by Adler et al. (2008). The introduced method was a nonlinear 3D tomographic least-squares inversion approach of residual moveout in PSDM common-image gathers. The velocity model used here is shown in figure 12. As it can be seen, the velocity along the crest of the structure is about 3400 to 3600 m/s, increasing to 3700 - 3800 m/s in the vicinity of the structural closure. In the next step, synthetic seismograms were constructed from sonic logs of 21 wells, (not shown here). A total of 36 well ties were made based on the 21 synthetic seismograms. A regression curve was fitted to the data and later used to aid the initial picking of seismic events at wells with no sonic logs. This time–depth curve was used for preliminary, approximate interpretations subject to adjustments necessitated by further seismic character evaluation and line intersection ties. The interpreted horizons on the seismic sections are shown in figures 13 and 14. Faults, anticline and different horizons are better identified in the up-thrown side of the faults, while in the down-thrown side; they could be traced by the similarity of the layers from the other side of the faults. To remove structural errors inherent in time migration, it is necessary to convert timemigrated images into the depth domain either by migrating the original data with a pre-stack depth migration algorithm or by depth migrating post-stack data after time de-migration, (Cameron et al., 2008). To better compare the final result of interpretation, which is depth map of the target, the industry standard method should be obtained for comparison. This standard method is performing conventional pre-stack time migration on the data, obtain a time map of the target horizon and then convert it to depth by a depth velocity model (Cameron et al., 2008).
(a)
(b) 12
Fig. 12 Velocity model used for depth imaging of seismic data processed by the introduced method. (a) A cube of the velocity model which shows trend of the anticline by red color. (b) A cross section that shows consistency of the velocity model in shape with an overthrust model.
In order to use result of the 3D CDS migrated result into the standard method, depth migrated result should be time converted. Many data processing and interpretation procedures are performed on time domain data. If such procedures are to be applied after a depth imaging project has been conducted, it is necessary to convert the data to the time domain from the depth domain. Although there will be a velocity-depth model available that was used for the depth migration, however, it is inappropriate to use the depth migration velocity model for depth-totime conversion if additional processing is to be performed on the pre-stack data in the time domain (Jones, 2009). The situation is less clear-cut if the objective is to compare the timeconverted depth image to check-shot or interval times. It should be noted that time conversion of depth migrated data is subtly different from depth conversion of time migrated data. For time migration, the velocity field is inherently smooth, this is not the case for depth migrated data (Jones, 2009). Hence we need to introduce a new, separate velocity field for the purpose of depth-to-time conversion. Thus a detail velocity model was obtained by the method introduced by Alaei (2006). This is an integrated procedure for migration velocity analysis in complex structures of thrust belts. Finally, the velocity model of the target formation and the related time contour map is shown in figure 15. Afterwards, this time map was converted to depth by the algorithm proposed by Cameron et al. (2008) for 2D and 3D data (figure 16a). Depth map obtained by conventional method is also shown in figure 16b. As it could be seen in figure 16a, there is high elevation difference in both sides of the over-thrust. There are also frequent small fault in top of the anticline that makes it difficult for a suitable seismic imaging, in case of introducing strong velocity changes. The apparent detail of the fault traces, correlations and contours imply a resolution greater than is warranted by the seismic coverage and quality, which proves that ability of applying 3D partial common diffraction surface method to provide enhanced seismic data for final depth imaging.
(a) (b) Fig. 13 Interpreting seismic sections obtained by the 3D partial common diffraction surface method, (a) In line number 480 and (b) In line number 720. 13
(a) (b) Fig. 14 Interpreting seismic sections obtained by the 3D partial common diffraction surface method, (a) In line number 1800 and (b) In line number 3100.
(a) (b) Fig 15 (a) The time map of the target formation. (b) Velocity model of the target formation used for time to depth conversion.
(a) (b) Fig. 16 (a) Depth map of the target formation obtained from the data processed by the 3D partial common diffraction surface method. (b) The same map obtained by the conventional method. 14
7.
Conclusion
The 3D partial common diffraction surface method was introduced here to overcome some of the problems of seismic imaging in complex geological structures. The result of the partial common diffraction surface followed by depth migration showed that it can serve as a method of enhancing pre-stack data for imaging complex structures. This mainly relates to enhancing weak events compared to other methods that such weak events may be suppressed in process of selecting the most coherent and/or dominant reflections/diffractions. These weak events are mostly diffractions that are related to structures like faults, overthrust, and salt bodies. The partial common diffraction surface method was applied on 3D synthetic and real land data set. After conventional processing, real data set were processed by the conventional PSDM and partial common diffraction surface methods. The velocity model for depth imaging was obtained by 3D tomographic least-squares inversion. The migrated sections showed that the partial common diffraction surface method can clearly define large and small scale faults beneath a thick layer of anhydrite and salt. Truncation of layers in the overthrust zone, extension of the anticline and continuity of reflectors were imaged better by the partial common diffraction surface method. Depth migrated data were converted to time to make a time map and underground contour map of the target formation. Time map were obtained by horizon picking in the enhanced result. Then this time map was converted to depth by hybrid velocity model. The final depth map was an enhanced high resolution image that shows much detail of faults and structures. Therefore, it concludes that the 3D partial common diffraction surface method provides suitable input for migration at least in semi-complex media.
References Adler, F., Baina, R., Soudani, M.A., Cardon, P., and Richard, J.B., 2008. Nonlinear 3D tomographic least-squares inversion of residual moveout in Kirchhoff prestack-depth-migration common-image gathers. Geophysics, 73(5), 13-23. DOI: 10.1190/1.2956427. Alaei, B., 2006. An integrated procedure for migration velocity analysis in complex structures of thrust belts. Journal of Applied Geophysics, 59, 89-105. DOI: 10.1016/j.jappgeo.2005.08.004. Alkhallifah, T., 1995. Gaussian beam depth migration for anisotropic media. Geophysics, 60, 1474-1484. DOI: 10.1190/1.1443881. Baykulov M., 2009. Seismic imaging in complex media with the Common Reflection Surface stack. Ph.D. Thesis, Hamburg University. Baykulov, M., and Gajewski, D., 2009. Prestack seismic data enhancement with partial common-reflection-surface (CRS) stack. Geophysics 74, 49-58. DOI: 10.1190/1.3106182. Baykulov, M., Dummong, S., and Gajewski, D., 2011. From time to depth with CRS attributes. Geophysics, 76, 151-155. DOI: 10.1190/1.3580607.
15
Bakhtiari Rad, P., Schwarz, B., Vanelle, C., and Gajewski, D., 2014. Common reflection surface (CRS) based pre-stack diffraction separation. SEG conference, Denver, 4208-4212. DOI: 10.3997/2214-4609.201413007. Bergler, S., 2001. The Common reflection surface stack for common offset - Theory and application. Master’s thesis, University of Karlsruhe. Biondi B., 2006. 3D Seismic imaging, investigations in geophysics, 14, SEG Publishing, Tulsa. Cameron, M., Fomel, S. and Sethian, J., 2008. Time-to-depth conversion and seismic velocity estimation using time-migration velocity. Geophysics, 73, 205-210. DOI: 10.1190/1.2967501. Cristini, A., Cardone, G., Marchetti, P., Zambonini, R., Hubral, P. and Mann, J., 2003. 3D zero offset CRS stack for narrow azimuth data: formulation and examples: EAGE/SEG Research Workshop, Trieste, Italy. de Bazelaire E., 1988. Normal moveout revisited - inhomogeneous media and curved interfaces. Geophysics, 53, 143-157. DOI: 10.1190/1.1442449. Dell, S., and Gajewski, D., 2011. Common-reflection-surface-based workflow for diffraction imaging: Geophysics, 76, (5), 187–195. DOI: 10.1190/geo2010-0229.1. Doherty, S.M., 1975. Structure independent seismic velocity estimation. Ph.D. Thesis, Stranford University. Duveneck, E., 2004. Velocity model estimation with data derived wavefront attributes. Geophysics, 69, 61-72. DOI: 10.1190/1.1649394. Dümmong, S. and Gajewski, D., 2008. A multiple suppression method via CRS attributes. SEG meeting, Expanded abstracts. Las Vegas, ID: 2008-2531. Druzhinin, A., MacBeth, C. and Hitchen, K., 1999. Prestack depth imaging via modelindependent stacking. Journal of Applied Geophysics, 42, 157–167. DOI:10.1016/S09269851(99)00036-1. Garabito, G., Cruz, J. C., Hubral, P., and Costa, J., 2001. Common reflection surface stack: A new parameter search strategy by global optimization. In Expanded abstracts, 71 st SEG meeting, 2009–2012. Garabito, G., Oliva, P.C. and Cruz, J.C.R., 2011. Numerical analysis of the finite offset common reflection surface traveltime approximations. Journal of Applied Geophysics, 74, 89– 99. DOI:10.1016/j.jappgeo.2011.03.005. Gelchinsky, B., Berkovitch, A., and Keydar, S., (1999), Multifocussing homeomorphic imaging Part 1. Basic concepts and formulas. Journal of Applied Geophysics, 42, 229–242. DOI:10.1016/S0926-9851(99)00038-5. Heilmann, Z., 2007. CRS-stack-based seismic reflection imaging for land data in time and depth domain. Logos Verlag, Berlin. Hertweck, T., Jäger, C., Mann, J., Duveneck, E., and Heilmann, Z., 2004. A seismic reflection imaging workflow based on the Common-Reflection-Surface (CRS) stack: theoretical background and case study. In Expanded abstracts, 74th SEG meeting, ID: SP 43.
16
Hertweck, T., Schleicher, J., and Mann, J., 2007. Data stacking beyond CMP. The Leading Edge, 26 (7), 818–827. DOI: 10.1190/1.2756859. Hill, N.R., 2001. Prestack Gaussian-beam depth migration. Geophysics, 66, 1240-1250. DOI: 10.1190/1.1487071. Höcht, G., Ricarte, P., Bergler, S. and Landa, E., 2009. Operator-oriented CRS interpolation. Geophysical Prospecting, 57, 957–979. DOI: 10.1111/j.1365-2478.2009.00789.x Höcht, G., de Bazelaire, E., Majer, P., and Hubral, P., 1999. Seismic and optics: hyperbolae and curvatures. Journal of Applied Geophysics, 42, 261–281. DOI: 10.3997/22144609.201407918 Jäger, R., Mann, J., Höcht, G., and Hubral, P., 2001. Common reflection surface stack: Image and attributes. Geophysics, 66, 97–109. DOI: 10.1190/1.1444927. Jones, I.,F., 2009. Time conversion of depth migrated data. First break, 27, 51-55. DOI: 10.3997/1365-2397.2009013 Mann, J., 2001. Common-Reflection-Surface stack and conflicting dips. In Extended abstracts, 63rd EAGE Conference Session P077. DOI: 10.1190/1.1815811. Marchetti, P., Follino, P., Borrini, D., and Zamboni, E., 2006. 3D CRS processing: A better use of pre-stack data. 6th Exploration Conference and Exposition on Petroleum Geophysics, Kolkata. Müller, T., 1998. Common reflection surface stack versus NMO/stack and NMO /DMO/stack. th 60 EAGE Conference,, Extended abstracts. DOI: 10.3997/2214-4609.201408166. Müller, T., 1999. The common reflection surface stack method—seismic imaging without explicit knowledge of the velocity model. Ph.D. Thesis, University of Karlsruhe, Germany. Müller. A., 2003. The 3D common-surface-reflection stack-theory and application. Master thesis, University of Karlsruhe. Prüssmann, P., Mann, J., and Buske, S., 2007. High resolution 3D CRS imaging for seismic assessment and monitoring of subsurface CO2 storage sites. First French-German Symposium on Geological Storage of CO2. Shahsavani, H., Mann. J., Piruz. I. and Hubral. P., 2011. A Model based approach to the Common diffraction surface Stack. 73rd EAGE Conference, SPE EUROPEC ,Vienna, P074. Sherkati, S., and Letouzey, J., 2004. Variation of structural style and basin evolution in the central Zagros (Izeh zone and Dezful Embayment), Iran, Marine Petroleum Geology, 21, 535 – 554. DOI: 10.1016/j.marpetgeo.2004.01.007. Sherkati, S., Molinaro, M., de Lamotte, D., and Letouzey, J., 2005., Detachment folding in the central and eastern Zagros fold-belt (Iran), Structural Geology, 27, 1680 – 1696. DOI: 10.1016/j.jsg.2005.05.010. Sherkati, S. Letouzey, J., and de Lamotte, D., 2006. Central Zagros fold-thrust belt (Iran): New insights from seismic data, field observation, and sandbox modeling. Tectonics, 25(4). DOI: 10.1029/2004TC001766. Sherwood, J.W.C., Schultz, P.S., and Judson, D.R., 1978. Equalizing the stacking velocities of dipping events via devilish. 48th SEG meeting.
17
Soleimani, M., 2015. Seismic imaging of mud volcano boundary in the east of Caspian Sea by common diffraction surface stack method. Arabian journal of geoscience, 8(6), 3943-3958. DOI: 10.1007/s12517-014-1497-5. Soleimani M. and Mann J., 2008. Merging aspects of DMO correction and CRS stack to account for conflicting dip situations. Annual WIT Report, 159–166. Soleimani, M., Piruz, I., Mann, J. and Hubral, P., 2009a. Solving the problem of conflicting dips in common reflection surface (CRS) stack. 1st International Petroleum Conference and Exhibition, Shiraz, Iran, A39. Soleimani, M., Piruz, I., Mann, J., and Hubral, P., 2009b. Common reflection surface stack; accounting for conflicting dip situations by considering all possible dips. Journal of seismic exploration, 18, 271–288. Soleimani, M., Mann, J., Adibi, E., Shahsavani, M., and Piruz, I., 2010, Applying the CRS stack method to solve the problem of imaging of complex structures in the Zagros overthrust, south west Iran. 72nd, EAGE conference. Barcelona, P556. Soleimani, M., Adibi E., and Mann, J., 2011. Imaging in complex structures by post-stack time migration and CRS stack. 12th International Congress of the Brazilian Geophysical Society, SBGf. DOI: 10.1190/sbgf2011-212. Yang, K., Bao-shu Chen, B.S., Wang, X.J., Yang, X.J., and Liu, J.R., 2012. Handling dip discrimination phenomenon in common-reflection-surface stack via combination of outputimaging-scheme and migration/demigration. Geophysical Prospecting, 60, 255–269. DOI: 10.1111/j.1365-2478.2011.00981.x. Yilmaz, Ö. and Claerbout, J. F., 1980. Prestack partial migration. Geophysics, 45(12), 1753– 1779. DOI: 10.1190/1.1441064. Yilmaz, Ö., 2001. Seismic data analysis, SEG publishing, Tulsa. Highlights
This research paper introduces and application of a developed method for imaging of 3D seismic data in reservoir located in complex geological media. The proposed method produces an enhanced prestack data which require smooth and simple velocity model for imaging. The paper improves the previously introduced 2D CDS method to 3D CDS. The paper uses advantages of the partial form in the CRS method to the CDS method.
18
www.elsevier.com/locate/petrol
PII: DOI: Reference:
S0920-4105(16)30062-6 http://dx.doi.org/10.1016/j.petrol.2016.02.019 PETROL3362
To appear in: Journal of Petroleum Science and Engineering Received date: 19 April 2015 Revised date: 2 December 2015 Accepted date: 22 February 2016 Cite this article as: M. Soleimani, Seismic imaging by 3D partial CDS method in complex media, Journal of Petroleum Science and Engineering, http://dx.doi.org/10.1016/j.petrol.2016.02.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Seismic imaging by 3D partial CDS method in complex media M. Soleimani* Faculty of Mining, Petroleum and Geophysics, University of Shahrood, Shahrood, Iran.
[email protected] *
Tel-Fax: +98 (23) 3239 5509
Abstract
Seismic imaging in complex geological structures is state of art to deliver a high quality section for further geological interpretation. Conventional processing and imaging methods will fail in case of facing complex media, strong heterogeneity and poor quality data. Thus, conventional seismic imaging methods have to be improved inevitably or new approaches have to be extended to produce valid images. The partial common diffraction surface method is based on improving the common diffraction surface operator to handle lateral heterogeneities in complex media. In this study, the partial diffraction surface method was improved to 3D. This method enhances quality of the 3D pre-stack data to be used for further depth imaging. It can also resolve some of the possible ambiguities in geological interpretation from seismic images. This also will enhance the quality of the seismic images, where it suffers from conflicting dips problem. A time variant linear function defines the offset range for each zero offset sample in the stacking operator. This function controls the offset increment for the partial form of the operator in the time-midpoint-offset domain. It should be noted that the velocity model can be considered as smooth and simple as possible in this method for further depth imaging. The introduced method was applied on a 3D synthetic and a real land data set. Results from both data example show the ability of the new method in enhancing the quality of the seismic section in the presence of faults and lateral velocity heterogeneity. This enhancement achieved by a simple velocity model in a strong complex media. Consequently, depth migrated result of the real land data underwent a time conversion to create a time underground contour map, for comparison with the maps obtained with conventional methods. More details of the target formation were observed on the map used result of the proposed method for mapping.
Keywords: Seismic imaging, common diffraction surface, lateral velocity heterogeneity, Complex media, Zagros overthrust. 1.
Introduction
Seismic imaging methods that are based on ray tracing techniques face with problem in imaging the complex media data. Geologically boundaries of complex structures such as salt dome, overthrust zone and faulted - folding system are illuminated by different ray paths. Tracing the paths of these rays cannot be properly handled by conventional imaging techniques (Yilmaz, 2001). Significant improvement in quality of imaging seismic data from such media could be achieved both in performing an accurate velocity model building procedure and/or applying advanced migration algorithm. In imaging complex structures with strong heterogeneity, the conventional migration methods based on the Kirchhoff integral, have
1
theoretical and practical shortcomings. Thus combination of accurate velocity model building and appropriate migration algorithm is the key to successful imaging (Biondi, 2006). Among the problems of imaging in complex structures, the conflicting dips situation has been investigated by different methods. The problem of conflicting dips was also observed in application of the new introduced common reflection surface (CRS) method. Concept of the CRS stack method was developed by Jäger et al. (2001) as a reasonable extension of the traveltime equation to fully exploit the information contained in the pre-stack data. Different applications were derived for the CRS stack method ranging from multiple suppression (Dümmong and Gajewski, 2008), pre-stack data enhancement and regularization (Baykulov et al., 2011) and velocity model building for depth migration (Duveneck, 2004). Dell and Gajewski (2011) had introduced a CRS based workflow for diffraction imaging that was followed by Bakhtiari Rad et al. (2014), which used common reflection surface, based on pre-stack diffraction separation method only to image diffractions rather than reflection events. On the other hand, some of the new introduced data driven imaging techniques mainly aimed to simulate zero offset (ZO) sections from multi-coverage seismic reflection data. Among of these new methods, several velocity independent methods have been introduced. Druzhinin et al. (1999) developed a scpecial prestack depth migration technique which avoids the necessity of constructing a macro velocity model. Gelchinsky et al. (1999) introduced the multi-focusing method for application in multicoverage data. Müller (1999) also formulated the CRS method as a macro velocity model independent technique. Garabito et al. (2001) later on introduced a new parameter search strategy by global optimization. Hertweck et al. (2004) also introduced the workflow of obtaining the kinematic wavefield attributes from the CRS equation. These methods were data driven in the way that they: (a) use a multi-parameter moveout formulas, where the moveout parameters are derived based on the coherency analysis, and (b) do not need explicit knowledge of the macro velocity model (Soleimani et al., 2010).The data drive CRS method provide result with high signal to noise ratio with increased vertical resolution.The CRS method imposes some approximations on the wavefield parameter estimation, which indeed less information of velocity model is required (Prüssmann et al., 2007). However, the CRS stack method sufferes from the problem of conflicting dips. Several strategies were introduced to resolve this problem in the CRS method. Höcht et al. (2009) introduced two techniques, called as target-oriented (TO) scheme and an alternative operator-oriented (OO) scheme to remove this problem in 2D ZO CRS. They have mentioned that a ‘proper’ handling of conflicting dip situations requires to identify the different events in terms of kinematics and waveform. Soleimani et al. (2011) introduced a modification of the CRS method that was called the common diffraction surface (CDS) or data based CDS stack method. Shahsavani et al. (2011) introduced the model based CDS stack method that decreased the large computation time of the CDS stack method. Garabito et al. (2011) introduced the surface operator of the CDS method that was a special case of the data based CDS stack method when only the dominant surface is considered as the operator. Yang et al. (2012) called the phenomenon caused by conflicting dips problem as ‘dip discrimination phenomenon’. This method was called the CRS stack with the output imaging 2
scheme (CRS-OIS). However, Soleimani (2015) showed that application of the 2D CDS method could produce enhanced migrated image of the semi-complex geological structures. 2.
The Conventional 3D CRS imaging
The CRS stack is a multi-parameter stacking method that uses any contributions along any realization of its operator. These contributions are tested by coherence analysis for each ZO sample, and the set of attributes, which yield the highest coherency, are accepted as the parameters of the optimum operator to perform the actual stack (Bergler, 2001). The CRS equation is described by three parameters α, RN and RNIP, as (Müller, 1998): 2
2sinα(xm x0 ) 2t 0cos 2α (xm x0 )2 h2 2 thyperbolic (xm ,h) = t0 + + + v0 v0 RN RNIP
(1)
where h is the half offset, Δx=(xm-x0) is the displacement with respect to the considered common mid-point (CMP) position, t0 corresponds to the ZO two-way traveltime, α is the emergence angle of the ZO ray, RN and RNIP are the radius of curvature of the normal (N) wave and the normal incidence point (NIP) wave, respectively, and finally v0 is the near surface velocity. The CRS stack method assigns merely one optimum stacking operator for each ZO sample to be simulated (Heilmann, 2007). However, in situation where different events in a seismic section intersect with each other and/or themselves, a single stacking operator for that ZO sample would be no longer appropriate. Thus, Mann (2001) proposed to allow small discrete number of stacking operators for a particular ZO sample. Therefore in conflicting dip situations, more than one operator for each ZO sample in ZO stack simulation would be used. Müller (2003) formlulated the 2D CRS operator into 3D. The 3D CRS operator is based on a second order approximation of traveltimes in midpoint and half-offset coordinates, (Cristini et al., 2003). The stacking surface in the time space domain, used to compute the CRS volume, depends on eight parameters that describe the shape and direction of normal ray and wavefront curvatures for Normal (N) and NIP wave. The method does not require the exact knowledge of the velocity model, a priori knowledge of the near surface constant velocity being sufficient, (Marchetti et al., 2006). The 3D CRS processing is done for each sample by summing up all events located on a traveltime surface in the five dimensional pre-stack data hyper volume (Müller, 2003): 2
2t 2 2 T T thyp t0 wm 0 mT Qzyz N N Qxyz m hT Pzyz N NIP Pzyz h v v
(2)
Where thyp describes the 3D hyperbolic traveltime approximation related to the local x, y system, centered at point X0 x0, y0. The parameter t0, is the two-way ZO traveltime at point P0x0, y0, t0, v is the surface velocity, w is the projection of the unit vector of the emerging normal ray at point X0 onto the ground surface. It contains the parameters , i.e. the angles that determine the direction of propagation of the two hypothetical wavefronts, NIP wave and N wave, along the emerging central (normal) ray. Curvature of the NIP and N-waves are defined by the NNIP and the NN matrices (Müller, 2003):
3
1 R N i i ,max 0
0 1 Ri ,min
(3)
The angle parameter , are defined by auxillary parameters P, Q and D as follows: Pxyz Dz ( ) Dy ( ) Dz ( NIP ),
(4a)
cos sin QDxyzz D z ( ) Dy ( ,) Dz (N,),NIP , N sin cos cos D y 0
0 1
cos D z sin
(4b) (5)
sin , cos
, NIP , N
cos eight 0 parameters at X0 are, final DThe y 0 1 are obtained with only eight parameters
NIP N, RNIP-max, RNIP-min, RN-max and RN-min. The two searches (one three-parametric and one fiveparametric search). The simplified flowchart of obtaining these parameters is shown in figure 1. pre-stack data Step I
VNMO,min VNMO,max
Automatic 3D CMP Stack
v
one three-parametric search for VNMO,min ,
CMP Coherency CMP Stack
VNMO,max and v.
Step II
,
3D ZO Stack
N
one five-parametric search for ,
RN,min and
N
RN,min, RN,max ZO Coherency ZO Stack
RN,max
Step III Determination of RNIP parameters
3D CRS Stack Stack using the eight paramenters within the projectet Fresnel Zones using the complete CRS operator
RNIP,min RNIP,max
NIP
From VNMO,min , VNMO,max and v (Step I) and ,(Step II), determination of
RNIP,min, RNIP,max and
NIP
3D CRS Stack
Fig. 1 Flowchart of the search strategy of the 3D CRS stack (Müller, 2003).
3.
Impelemntation of the 3D partial CDS
Soleimani and Mann (2008) introduced the CDS stack method for solving the conflicting dips problem in the 2D CRS stack. In this method, for each sample in time domain, a stacking surface is performed in a range of angles. Therefore, there exist many stacking surfaces for one sample 4
in (t, xm, h) domain that makes a volume of stacking operators and all of them will contribute into stacking for that sample. Thus, the conflicting dips will be treated well by the CDS operators. Solving the problem of conflicting dips in this way will enhance the usually weak diffraction events in the stacked section (Soleimani et al., 2009a). For true diffraction events, the radius of the NIP wavefront and the normal wavefront will coincide, RN = RNIP. Thus the only attribute to be searched for in a fixed emergence angle, is a combination of RN and RNIP that is called RCDS, (Soleimani et al., 2009b): 2
2 2sinα xm x0 + 2t 0cos α xm x0 2 + h 2 t 2 xm h = t 0 + v0 v0 RCDS
(6)
By implicit knowledge of the RCDS, the shape of the operator could be defined. Baykulov and Gajewski (2009) calculate a stacking surface around a specified ZO point defined by its offset and traveltime coordinates in a chosen CMP location and perform the summation of data along that surface (Figure 2a). This strategy is called the partial CRS stack method. Baykulov (2009) had mentioned that the partial CRS stack method needs further investigations on the conflicting dips problem. However, prestack migration of noisy and low quality data produces migrated section of comparably lower quality than the post stack migration of the CRS stack. Baykulov (2009) showed that the CRS traveltime formula, where the dip of the reflector element is incorporated, is used to compute new partially-stacked CRS super-gathers, where each trace is a result of summation of data along the CRS stacking surface. The number and location of traces in the produced super-gathers can be defined. To completely remove the problem of conflicting dips and produce high quality super-gathers for further prestack migration, the idea of the partial CRS stack method was used here to modify the CDS stack equation. The new method differs from the partial CRS stack. This method, that is called the partial CDS, uses this idea for operator volume stack rather than only a surface stack. The partial CDS stack method, calculates collection of stacking surfaces around a specified point. This point is defined by its CMP number, zero offset traveltime, t0, and offset boundary related to that t0 in a zero offset section, and performs the summation of data along that surface. The result of summation is assigned to that sample. The operator equation in the partial CDS is same as the CDS operator with performing limitation in offset range according to one offset banding function. The partial CDS stack operator is shown as yellow surface related to the sample t0 (red point) in figure 2b. To have advantage of partial 2D CDS stack, the conventional 3D CRS stack method was improved here not only to 3D CDS stack, but into a partial form. Hence, the radius of normal wave and NIP wave should coincide. Therefore, eight parameters are reduce to five parameters while: NIP = N, RNIP-max = RN-max and RNIP-min = RN-min. Finally these five parameters are: CDS RCDS-max, and RCDS-min. The partial 3D CDS stack equation would be:
5
2
t
2 hyp
2 2t T T t0 wm 0 mT Qzyz NCDSQxyz m hT Pzyz NCDS Pzyz h v v
(7)
Figure 3 shows the proposed search startegy of the partial 3D CDS equation. By knowing the fact that the CDS stacking parameter varies smoothly for a smooth velocity model, the rays can be calculated on a relatively coarse emergence angle grid. In contrast, the semblance and the stack itself are quite sensitive to variations of the emergence angle rather than the initial stacking velocity. Thus, stacking and semblance value perform on a finer emergence angle grid using interpolated stacking parameters. This shows that result of the CDS stack is more sensitive to emergence angle increment rather than the accuracy of the near surface velocity model.
Fig. 2 a) The partial CRS stack performs summation of data around a specified point on a CMP traveltime curve (magenta line) and assigns result to the same point in a newly generated CRS supergather. b) The partial CDS stacking surface shown with a yellow color coincides locally with the common offset travel time surface, but may be limited in size. (Figure (a) from Baykulov (2009) and figure (b) from Soleimani and Mann (2008)). pre-stack data VNMO,min VNMO,max
Step I Automatic 3D CMP Stack
v Offset banding CMP Coherency CMP Stack
one three-parametric search for VNMO,min ,
VNMO,max and v.
Step III
Step II
3D Partial CDS Stack
3D CDS Stack
Stack using the five paramenters within the projectet Fresnel Zones using the complete CDS operator
one five-parametric search for , CDS
RCDS,min and
RCDS,max
Enhanced prestack data by 3D Partial CDS, (ready for depth imaging)
` 6
Fig. 3 Flowchart of the search strategy of the 3D partial CDS method.
4.
Application of the partial CDS on synthetic data
To verify the improved search strategy of the CDS stack, a 3D synthetic data set was generted for processing. The synthetic model data obtain by forward modeling on a five layered dome shape model. Synthetic data was contaminated by random noise for better evaluation. The results of the CRS stack process were three optimized kinematic wavefield attributes, one coherency section and one optimized stacked section. However, the results of the CDS stack process are just one optimized stacked section and the section that shows the number of used traces in aperture. Therefore only the stacked sections are comparable. Figure 4a shows the velocity model used for synthetic data generation. Figure 4b shows the coherency result of application the CDS method on the synthetic data. Results of the CRS and the CDS method on the data are shown in Figure 4c and 4d. As it could be seen, the stacked section obtained by the CDS method could clearly better image truncation of the reflector, conflicting events and continuity of the reflectors. Result proved that the proposed method could be applied on real data.
(a)
(b)
(c)
(d) Fig. 4- (a) Velocity model used for synthetic data generation by ray tracing. (b) Coherency section of the CDS result. (c) The CRS stacked result and (d) The CDS stacked result on the synthetic data.
5.
Application of the 3D partial CDS on real land data
To see how the new introduced method could overcome some of the problem of imaging of complex structures, it was applied on a complex 3D land data. The related seismic data was 7
obtained on a mountainous topography with complex geology that makes strong velocity heterogeneity in the media. Figure 5a shows location of the study area, in SW of the Zagros overthrust belt. The study structure is an extensively faulted anticline with 74 km long and 6 to 8 km wide, (Figure 5b). The northern 60 km of the southwest flank is defined by a reverse fault with up to 2100 m of throw. Normal faults occur on both flanks and the crest, but are most common on the southwest flank, between the crest and the reverse fault (Figure 5b). Figure 6 shows the stratigraphic column of the study area. The major problematic lithology here is the anhydrite - salt layer called the Gachsaran Formation that is responsible for velocity contrast.
(a) (b) Fig. 5 (a) Location of the Zagros overthrust belt in SW of Iran and location of the study area shown by a red rectangle. (b) Structural map of the target anticline and area of 3D seismic data. Locations of the seismic lines are shown by bold black lines. (Figure (a) from Sherkati et al. (2005))
Fig. 6 Stratigraphic column of the study area that shows thick layer of anhydrate - salt in the area (Sherkati et al., 2004)
8
The target formation here is the Asmari formation, which is the major reservoir formation in SW Iran, with an anticlinal shape beneath its cap rock, the Gachsaran formation. The enhanced pre-stack data obtained by the 3D partial CDS method underwent depth imaging by the Gaussian beam migration (GBM) method, (Alkhaliffah, 1995; Hill, 2001). To see the results of the processing, two in-lines and one cross line of the seismic cube were extracted for visualization. Figure 7a shows result of the conventional prestack depth migration (PSDM) and figure 7b shows depth migrated result of the partially enhanced pre-stack data. As it can be seen, migration of the partially enhanced data, images more steep dips in the section. The anticline is dominated in both sides by faults, which was imaged better by migration of partially enhanced data. Image section in figure 7b also shows more continuity in reflectors below the anhydrite and salt layer. Truncations of layers are also better identified by enhancing data with 3D partial CDS method, which helps for easier structural interpretation. Thick layer of anhydrite is also imaged with more continuity by new method. The next in-line for presenting the result selected from a part of the anticline that is dominated by faults and difficult to define drilling path. Figure 8a shows result of the conventional PSDM and figure 8b shows depth migrated result of partially enhanced prestack data. The other major problem besides imaging those faults is to obtain higher quality image for imaging small scale faults in the left part and inside the anticline. Increasing continuity of the reflectors is obvious in migration result of the partially enhanced data, and two large scale faults on both sides of the anticline are also better imaged in figure 8b. Enhanced parts of the section are highlighted by the black rectangles. Faults inside the anticline, that has important effect on defining drilling path, are better imaged in the result of the new method. As it can be seen, some small faults are only imaged in figure 8b, while in figure 8a, only a blurred picture of small scale faults inside the anticline is obtained. In both sections, the quality of the seismic data appears to be affected, at least in some parts, by the variation of thickness of the Gachsaran (anhydrite and salt) formation. To see better result of application the 3D partial CDS method on data, a cross line section was shown that does not suffers from the problem mentioned above. Figure 9a shows a cross section result of conventional PSDM and figure 9b shows a cross section depth migrated result of partially enhanced pre-stack data. The final migrated result shown in figure 9b has better quality with more continuity in the reflectors, especially in the right down side of the section. The pinch out imaged below the CDP number 1700 in depth of 1600 m is well imaged in figure 9b, while it is not imaged in figure 9a. For better understanding of improvement on results, two depth slices from the seismic cube are also shown in figures 10 and 11. Enhanced parts of the sections are highlighted by the black rectangles. Figures 10a and 11a show depth slices result of conventional PSDM and figures 10b and 11b show the related depth slices result of depth migration of partially enhanced pre-stack data. The quality increasing in slices obtained by the migration of the partially enhanced data is obvious. Tight lines are well imaged in figures 10b and 11b. Almost all of the reflectors on both side of the anticline are imaged by the new method. Thus, it could be concluded that the pre-stack depth migration of the pre-stack data enhanced by 3D partially CDS method can be applied to complex structures to resolve some of the ambiguities of imaging in such media. In all the examples shown, horizons
9
are more continuous in the migration result of partially enhanced data and the increasing signal to noise ratio is obvious.
(a) (b) Fig. 7 Migrated section of in-line number 1800. (a) Result of conventional PSDM with Gaussian beam algorithm and (b) result of depth migration of partially enhanced pre-stack data.
(a) (b) Fig 8 Migrated section of in-line number 3100. (a) Result of conventional PSDM with Gaussian beam algorithm and (b) result of depth migration of partially enhanced pre-stack data.
10
(a) (b) Figure 9- Cross section number 700. (a) PSDM with GBM and (b) migrated result of partially enhanced pre-stack data, migrated by pre-stack Gaussian beam depth migration.
(a) (b) Figure 10- Depth slice of depth 1200m. (a) PSDM (with GBM) and (b) migration of partially enhanced pre-stack data, migrated by pre-stack Gaussian beam depth migration.
(a) (b) Figure 11- Depth slice of depth 1700m. (a) PSDM (with GBM) and (b) depth slices result of depth migration of partially enhanced pre-stack data, migrated by pre-stack Gaussian beam depth migration.
6.
Structural interpretation and velocity model considerations
The Gachsaran formation governs the velocity of the whole swelling part. In the center of the Gachsaran Formation, salt is reported as the main deposit. In areas where the amount of anhydrite is high, the expected total interval velocity would also increase. Therefore, the Gachsaran ridges show anomalous seismic velocity behavior (Sherkati et al., 2006). Usually conventional time-migrated seismic sections are distorted and obscure owing to the presence of inflated Gachsaran bodies (salt thickening) and related lateral velocity variations. In such conditions, strong lateral velocity variations, related to the lithology contrasts between steeply dipping layers bend the seismic rays like an optical lens and distort the sub-surface image. The 11
new method introduced here is considered to be the appropriate method for imaging targets in the presence of heterogeneous overburden, especially in the presence of a salt body. The only way to remove effects of salt in seismic image is to define the velocity variations by building a velocity model and performing depth migration, which compensates for ray-bending propagation effects. The velocity model used here was obtained by the method introduced by Adler et al. (2008). The introduced method was a nonlinear 3D tomographic least-squares inversion approach of residual moveout in PSDM common-image gathers. The velocity model used here is shown in figure 12. As it can be seen, the velocity along the crest of the structure is about 3400 to 3600 m/s, increasing to 3700 - 3800 m/s in the vicinity of the structural closure. In the next step, synthetic seismograms were constructed from sonic logs of 21 wells, (not shown here). A total of 36 well ties were made based on the 21 synthetic seismograms. A regression curve was fitted to the data and later used to aid the initial picking of seismic events at wells with no sonic logs. This time–depth curve was used for preliminary, approximate interpretations subject to adjustments necessitated by further seismic character evaluation and line intersection ties. The interpreted horizons on the seismic sections are shown in figures 13 and 14. Faults, anticline and different horizons are better identified in the up-thrown side of the faults, while in the down-thrown side; they could be traced by the similarity of the layers from the other side of the faults. To remove structural errors inherent in time migration, it is necessary to convert timemigrated images into the depth domain either by migrating the original data with a pre-stack depth migration algorithm or by depth migrating post-stack data after time de-migration, (Cameron et al., 2008). To better compare the final result of interpretation, which is depth map of the target, the industry standard method should be obtained for comparison. This standard method is performing conventional pre-stack time migration on the data, obtain a time map of the target horizon and then convert it to depth by a depth velocity model (Cameron et al., 2008).
(a)
(b) 12
Fig. 12 Velocity model used for depth imaging of seismic data processed by the introduced method. (a) A cube of the velocity model which shows trend of the anticline by red color. (b) A cross section that shows consistency of the velocity model in shape with an overthrust model.
In order to use result of the 3D CDS migrated result into the standard method, depth migrated result should be time converted. Many data processing and interpretation procedures are performed on time domain data. If such procedures are to be applied after a depth imaging project has been conducted, it is necessary to convert the data to the time domain from the depth domain. Although there will be a velocity-depth model available that was used for the depth migration, however, it is inappropriate to use the depth migration velocity model for depth-totime conversion if additional processing is to be performed on the pre-stack data in the time domain (Jones, 2009). The situation is less clear-cut if the objective is to compare the timeconverted depth image to check-shot or interval times. It should be noted that time conversion of depth migrated data is subtly different from depth conversion of time migrated data. For time migration, the velocity field is inherently smooth, this is not the case for depth migrated data (Jones, 2009). Hence we need to introduce a new, separate velocity field for the purpose of depth-to-time conversion. Thus a detail velocity model was obtained by the method introduced by Alaei (2006). This is an integrated procedure for migration velocity analysis in complex structures of thrust belts. Finally, the velocity model of the target formation and the related time contour map is shown in figure 15. Afterwards, this time map was converted to depth by the algorithm proposed by Cameron et al. (2008) for 2D and 3D data (figure 16a). Depth map obtained by conventional method is also shown in figure 16b. As it could be seen in figure 16a, there is high elevation difference in both sides of the over-thrust. There are also frequent small fault in top of the anticline that makes it difficult for a suitable seismic imaging, in case of introducing strong velocity changes. The apparent detail of the fault traces, correlations and contours imply a resolution greater than is warranted by the seismic coverage and quality, which proves that ability of applying 3D partial common diffraction surface method to provide enhanced seismic data for final depth imaging.
(a) (b) Fig. 13 Interpreting seismic sections obtained by the 3D partial common diffraction surface method, (a) In line number 480 and (b) In line number 720. 13
(a) (b) Fig. 14 Interpreting seismic sections obtained by the 3D partial common diffraction surface method, (a) In line number 1800 and (b) In line number 3100.
(a) (b) Fig 15 (a) The time map of the target formation. (b) Velocity model of the target formation used for time to depth conversion.
(a) (b) Fig. 16 (a) Depth map of the target formation obtained from the data processed by the 3D partial common diffraction surface method. (b) The same map obtained by the conventional method. 14
7.
Conclusion
The 3D partial common diffraction surface method was introduced here to overcome some of the problems of seismic imaging in complex geological structures. The result of the partial common diffraction surface followed by depth migration showed that it can serve as a method of enhancing pre-stack data for imaging complex structures. This mainly relates to enhancing weak events compared to other methods that such weak events may be suppressed in process of selecting the most coherent and/or dominant reflections/diffractions. These weak events are mostly diffractions that are related to structures like faults, overthrust, and salt bodies. The partial common diffraction surface method was applied on 3D synthetic and real land data set. After conventional processing, real data set were processed by the conventional PSDM and partial common diffraction surface methods. The velocity model for depth imaging was obtained by 3D tomographic least-squares inversion. The migrated sections showed that the partial common diffraction surface method can clearly define large and small scale faults beneath a thick layer of anhydrite and salt. Truncation of layers in the overthrust zone, extension of the anticline and continuity of reflectors were imaged better by the partial common diffraction surface method. Depth migrated data were converted to time to make a time map and underground contour map of the target formation. Time map were obtained by horizon picking in the enhanced result. Then this time map was converted to depth by hybrid velocity model. The final depth map was an enhanced high resolution image that shows much detail of faults and structures. Therefore, it concludes that the 3D partial common diffraction surface method provides suitable input for migration at least in semi-complex media.
References Adler, F., Baina, R., Soudani, M.A., Cardon, P., and Richard, J.B., 2008. Nonlinear 3D tomographic least-squares inversion of residual moveout in Kirchhoff prestack-depth-migration common-image gathers. Geophysics, 73(5), 13-23. DOI: 10.1190/1.2956427. Alaei, B., 2006. An integrated procedure for migration velocity analysis in complex structures of thrust belts. Journal of Applied Geophysics, 59, 89-105. DOI: 10.1016/j.jappgeo.2005.08.004. Alkhallifah, T., 1995. Gaussian beam depth migration for anisotropic media. Geophysics, 60, 1474-1484. DOI: 10.1190/1.1443881. Baykulov M., 2009. Seismic imaging in complex media with the Common Reflection Surface stack. Ph.D. Thesis, Hamburg University. Baykulov, M., and Gajewski, D., 2009. Prestack seismic data enhancement with partial common-reflection-surface (CRS) stack. Geophysics 74, 49-58. DOI: 10.1190/1.3106182. Baykulov, M., Dummong, S., and Gajewski, D., 2011. From time to depth with CRS attributes. Geophysics, 76, 151-155. DOI: 10.1190/1.3580607.
15
Bakhtiari Rad, P., Schwarz, B., Vanelle, C., and Gajewski, D., 2014. Common reflection surface (CRS) based pre-stack diffraction separation. SEG conference, Denver, 4208-4212. DOI: 10.3997/2214-4609.201413007. Bergler, S., 2001. The Common reflection surface stack for common offset - Theory and application. Master’s thesis, University of Karlsruhe. Biondi B., 2006. 3D Seismic imaging, investigations in geophysics, 14, SEG Publishing, Tulsa. Cameron, M., Fomel, S. and Sethian, J., 2008. Time-to-depth conversion and seismic velocity estimation using time-migration velocity. Geophysics, 73, 205-210. DOI: 10.1190/1.2967501. Cristini, A., Cardone, G., Marchetti, P., Zambonini, R., Hubral, P. and Mann, J., 2003. 3D zero offset CRS stack for narrow azimuth data: formulation and examples: EAGE/SEG Research Workshop, Trieste, Italy. de Bazelaire E., 1988. Normal moveout revisited - inhomogeneous media and curved interfaces. Geophysics, 53, 143-157. DOI: 10.1190/1.1442449. Dell, S., and Gajewski, D., 2011. Common-reflection-surface-based workflow for diffraction imaging: Geophysics, 76, (5), 187–195. DOI: 10.1190/geo2010-0229.1. Doherty, S.M., 1975. Structure independent seismic velocity estimation. Ph.D. Thesis, Stranford University. Duveneck, E., 2004. Velocity model estimation with data derived wavefront attributes. Geophysics, 69, 61-72. DOI: 10.1190/1.1649394. Dümmong, S. and Gajewski, D., 2008. A multiple suppression method via CRS attributes. SEG meeting, Expanded abstracts. Las Vegas, ID: 2008-2531. Druzhinin, A., MacBeth, C. and Hitchen, K., 1999. Prestack depth imaging via modelindependent stacking. Journal of Applied Geophysics, 42, 157–167. DOI:10.1016/S09269851(99)00036-1. Garabito, G., Cruz, J. C., Hubral, P., and Costa, J., 2001. Common reflection surface stack: A new parameter search strategy by global optimization. In Expanded abstracts, 71 st SEG meeting, 2009–2012. Garabito, G., Oliva, P.C. and Cruz, J.C.R., 2011. Numerical analysis of the finite offset common reflection surface traveltime approximations. Journal of Applied Geophysics, 74, 89– 99. DOI:10.1016/j.jappgeo.2011.03.005. Gelchinsky, B., Berkovitch, A., and Keydar, S., (1999), Multifocussing homeomorphic imaging Part 1. Basic concepts and formulas. Journal of Applied Geophysics, 42, 229–242. DOI:10.1016/S0926-9851(99)00038-5. Heilmann, Z., 2007. CRS-stack-based seismic reflection imaging for land data in time and depth domain. Logos Verlag, Berlin. Hertweck, T., Jäger, C., Mann, J., Duveneck, E., and Heilmann, Z., 2004. A seismic reflection imaging workflow based on the Common-Reflection-Surface (CRS) stack: theoretical background and case study. In Expanded abstracts, 74th SEG meeting, ID: SP 43.
16
Hertweck, T., Schleicher, J., and Mann, J., 2007. Data stacking beyond CMP. The Leading Edge, 26 (7), 818–827. DOI: 10.1190/1.2756859. Hill, N.R., 2001. Prestack Gaussian-beam depth migration. Geophysics, 66, 1240-1250. DOI: 10.1190/1.1487071. Höcht, G., Ricarte, P., Bergler, S. and Landa, E., 2009. Operator-oriented CRS interpolation. Geophysical Prospecting, 57, 957–979. DOI: 10.1111/j.1365-2478.2009.00789.x Höcht, G., de Bazelaire, E., Majer, P., and Hubral, P., 1999. Seismic and optics: hyperbolae and curvatures. Journal of Applied Geophysics, 42, 261–281. DOI: 10.3997/22144609.201407918 Jäger, R., Mann, J., Höcht, G., and Hubral, P., 2001. Common reflection surface stack: Image and attributes. Geophysics, 66, 97–109. DOI: 10.1190/1.1444927. Jones, I.,F., 2009. Time conversion of depth migrated data. First break, 27, 51-55. DOI: 10.3997/1365-2397.2009013 Mann, J., 2001. Common-Reflection-Surface stack and conflicting dips. In Extended abstracts, 63rd EAGE Conference Session P077. DOI: 10.1190/1.1815811. Marchetti, P., Follino, P., Borrini, D., and Zamboni, E., 2006. 3D CRS processing: A better use of pre-stack data. 6th Exploration Conference and Exposition on Petroleum Geophysics, Kolkata. Müller, T., 1998. Common reflection surface stack versus NMO/stack and NMO /DMO/stack. th 60 EAGE Conference,, Extended abstracts. DOI: 10.3997/2214-4609.201408166. Müller, T., 1999. The common reflection surface stack method—seismic imaging without explicit knowledge of the velocity model. Ph.D. Thesis, University of Karlsruhe, Germany. Müller. A., 2003. The 3D common-surface-reflection stack-theory and application. Master thesis, University of Karlsruhe. Prüssmann, P., Mann, J., and Buske, S., 2007. High resolution 3D CRS imaging for seismic assessment and monitoring of subsurface CO2 storage sites. First French-German Symposium on Geological Storage of CO2. Shahsavani, H., Mann. J., Piruz. I. and Hubral. P., 2011. A Model based approach to the Common diffraction surface Stack. 73rd EAGE Conference, SPE EUROPEC ,Vienna, P074. Sherkati, S., and Letouzey, J., 2004. Variation of structural style and basin evolution in the central Zagros (Izeh zone and Dezful Embayment), Iran, Marine Petroleum Geology, 21, 535 – 554. DOI: 10.1016/j.marpetgeo.2004.01.007. Sherkati, S., Molinaro, M., de Lamotte, D., and Letouzey, J., 2005., Detachment folding in the central and eastern Zagros fold-belt (Iran), Structural Geology, 27, 1680 – 1696. DOI: 10.1016/j.jsg.2005.05.010. Sherkati, S. Letouzey, J., and de Lamotte, D., 2006. Central Zagros fold-thrust belt (Iran): New insights from seismic data, field observation, and sandbox modeling. Tectonics, 25(4). DOI: 10.1029/2004TC001766. Sherwood, J.W.C., Schultz, P.S., and Judson, D.R., 1978. Equalizing the stacking velocities of dipping events via devilish. 48th SEG meeting.
17
Soleimani, M., 2015. Seismic imaging of mud volcano boundary in the east of Caspian Sea by common diffraction surface stack method. Arabian journal of geoscience, 8(6), 3943-3958. DOI: 10.1007/s12517-014-1497-5. Soleimani M. and Mann J., 2008. Merging aspects of DMO correction and CRS stack to account for conflicting dip situations. Annual WIT Report, 159–166. Soleimani, M., Piruz, I., Mann, J. and Hubral, P., 2009a. Solving the problem of conflicting dips in common reflection surface (CRS) stack. 1st International Petroleum Conference and Exhibition, Shiraz, Iran, A39. Soleimani, M., Piruz, I., Mann, J., and Hubral, P., 2009b. Common reflection surface stack; accounting for conflicting dip situations by considering all possible dips. Journal of seismic exploration, 18, 271–288. Soleimani, M., Mann, J., Adibi, E., Shahsavani, M., and Piruz, I., 2010, Applying the CRS stack method to solve the problem of imaging of complex structures in the Zagros overthrust, south west Iran. 72nd, EAGE conference. Barcelona, P556. Soleimani, M., Adibi E., and Mann, J., 2011. Imaging in complex structures by post-stack time migration and CRS stack. 12th International Congress of the Brazilian Geophysical Society, SBGf. DOI: 10.1190/sbgf2011-212. Yang, K., Bao-shu Chen, B.S., Wang, X.J., Yang, X.J., and Liu, J.R., 2012. Handling dip discrimination phenomenon in common-reflection-surface stack via combination of outputimaging-scheme and migration/demigration. Geophysical Prospecting, 60, 255–269. DOI: 10.1111/j.1365-2478.2011.00981.x. Yilmaz, Ö. and Claerbout, J. F., 1980. Prestack partial migration. Geophysics, 45(12), 1753– 1779. DOI: 10.1190/1.1441064. Yilmaz, Ö., 2001. Seismic data analysis, SEG publishing, Tulsa. Highlights
This research paper introduces and application of a developed method for imaging of 3D seismic data in reservoir located in complex geological media. The proposed method produces an enhanced prestack data which require smooth and simple velocity model for imaging. The paper improves the previously introduced 2D CDS method to 3D CDS. The paper uses advantages of the partial form in the CRS method to the CDS method.
18