Application of common reflection surface (CRS) to velocity variation with azimuth (VVAz) inversion of the relatively narrow azimuth 3D seismic land data

Lucky Kriski Muhtar; Wahyu Triyoso; Fatkhan Fatkhan

doi:10.1515/phys-2021-0051

Article Open Access

Application of common reflection surface (CRS) to velocity variation with azimuth (VVAz) inversion of the relatively narrow azimuth 3D seismic land data

Lucky Kriski Muhtar , Wahyu Triyoso and Fatkhan Fatkhan

Published/Copyright: July 29, 2021

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 19 Issue 1

Abstract

The fracture direction and its intensity are critical properties related to hydrocarbon characterization and identification. Both these properties have an essential role in identifying the direction of hydrocarbon migration, determining the sweet spot area, and optimizing the drilling design. The velocity variation with azimuth (VVAz) is a well-known method to estimate the fracture direction and its intensity. This method is of widespread interest because it predicts the properties based on seismic data without any practical constraints. Despite this interest, the technique requires rich azimuth 3D seismic data in our case, which is rare. This study aims to apply regularization and interpolation by including the wave front attributes based on the Common Reflection Surface (CRS) method before the VVAz inversion. The motivation of using the CRS method is to enrich the current azimuth of the 3D seismic data and improve the S/N ratio. The synthetic and the real 3D seismic data are evaluated to examine the interpolation scheme of the proposed CRS method’s performance. Based on the evaluation of the 3D seismic data after regularization, the amplitude versus offset (AVO) phenomena, and the VVAz inversion results are relatively consistent (or matched) with the model. A similar result is found for the case of real 3D seismic data. A significant positive correlation between the fracture intensity of FMI and the real seismic data of about 0.9 is obtained. Therefore, CRS can be used as a regularization and interpolation method before the VVAz inversion of the relatively narrow azimuth 3D seismic data.

Keywords: fracture direction and intensity; CRS; VVAz inversion

1 Introduction

The fracture direction and its intensity are critical properties related to hydrocarbon characterization and identification. These properties have an essential role in identifying the direction of hydrocarbon migration, determining the sweet spot area, and optimizing the drilling design. The use of 3D seismic data to identify fractures in different reservoirs was carried out in the previous study [1,2,3]. Based on ref. [4], the fracture detection based on seismic data could be identified by correlating offset and azimuth. According to ref. [5], there are two approaches to be used in fracture detection – the inversion of amplitude as a function of offset and azimuth (AVAz) and the inversion of velocity as a function of offset and azimuth (VVAz). Based on the previous study [5], AVAz shows better spatial resolution than VVAz. However, AVAz is more unstable in the inversion process to maintain the amplitude preservation when the data contain more noise than VVAz because of the inversion of travel time data [6]. Therefore, the VVAz inversion was applied in this study.

The advantage of VVAz is its ability to better estimate the fracture direction and intensity without using constraints well in the inversion process; for example, the ability to differentiate sweet spot areas in shale gas fields [7,8] and tight sandstone reservoirs [9,10]. Previous researches’ success stories were based on the rich azimuth 3D seismic land data, but, in our case, oil companies having such kind of data are rare.

In 3D prestack seismic data processing, there is a method to enrich the azimuth by the frequency-wavenumber (f-k) interpolation method [11]. Ref. [12] extends the f-k interpolation into a Four-Dimensional (4D) or Five-Dimensional (5D) application. The 4D or 5D f-k interpolation can interpolate between the receiver and shots line to enrich the azimuth of the 3D seismic data. The f-k interpolation method uses the non-alias low-frequency portion of the data to interpolate the alias high frequency. Therefore, it can handle the spatial domain correctly. The multidimensional f-k interpolation of 4D or 5D extracts information from non-aliased spatial axes to interpolate highly oriented axes properly. The synthetic and real data examples need to be obtained to examine the performance of the proposed interpolation scheme since the techniques to reconstruct the missing data in the frequency domain are based on mathematical transforms.

The stacking process in a time domain other than Normal Moveout (NMO) stacking to overcome complex structural conditions was first introduced by ref. [13]. Based on this concept, the Common Reflection Surface (CRS) stacking was developed in which the parameters of the stacking process include the dipping and curvature of the reflector in the subsurface. The CRS method has proven to enhance the S/N ratio of data for complex geological structural conditions [14,15,16]. The successful result of the CRS method in enhancing the real data has been shown by refs. [17,18].

The CRS method can be used as an alternative to the interpolation method that uses wave-front attributes, namely, quantities with specific physical meanings such as the angle of emergence and the curvature of the wave-front. In the CRS attributes, the structural information of subsurface features such as the dip and strike of a reflector is included. It is suggested that including the wave front attributes work on 5-D data space (e.g., common-midpoint coordinates in x and y, offset, azimuth, and time) will lead to a 5-D interpolation technique. Ref. [18] compared CRS and 5D interpolation as a regularization method for narrow azimuth of the 3D seismic data in complex geological conditions and concluded that CRS was a better method.

Thus, this study aims to apply 5D interpolation by including the wave front attributes (the CRS method) before VVAz inversion. The motivation of using the CRS method is to enrich the current azimuth of the 3D seismic data and improve the nominal fold and the S/N ratio. The synthetic and the real 3D seismic data are evaluated to examine the interpolation scheme of the proposed CRS method’s performance.

2 Data and methods

In this study, the synthetic and real data examples are prepared to examine the performance of the proposed interpolation scheme before VVAz inversion is applied. The synthetic data are developed based on the composition of the five formations with different physical properties in each layer, as shown in Figure 1a. In the synthetic model, the fourth layer is anisotropic and the other layers are isotropic. The model contains three classes of elastic properties, isotropy, HTI-1, and HTI-2, as shown in Figure 1b.

Figure 1

The synthetic data are developed based on the composition of the five formations with different physical properties in each layer (a). The fourth layer is anisotropic in the synthetic model, and the other layers are isotropic (b). The model’s acquisition layout and fold distribution are used to generate the synthetic data (c).

Moreover, a Direct Hydrocarbon Indicator (DHI) that is the AVO class 3 is added in the anticline’s crest of the fourth layer. In the synthetic data, which is 9 × 9 km areas of the elastic anisotropy, the 3D seismic data are generated with about 64 maximum fold and 50 m bin size. This template has a 500 m receiver line interval and a 500 m source line interval. There are 80 receivers on each line with a 100 m receiver interval and a 100 m source interval. The illustration of the template can be seen in Figure 1c.

The real 3D seismic data are based on the Teapot dome. It is the land 3d seismic data from Wyoming provided by the U.S. Department of Energy and RMOTC. The data have a maximum fold of about 56 with a bin size of 110 × 110 ft. The acquisition layout is a diagonal template with an angle between the source and receiver of 45°, and the strike from the perpendicular source of the teapot dome structure is N45°E. In contrast, the azimuth of the receiver has an east–west direction, as shown in Figure 2. The vibroseis source is used with the shooting template as shown in the following figure – a receiver interval of 220 ft with a line interval of 880 ft and a source interval of 220 ft with a line interval of 2,200 ft.

Figure 2

The acquisition layout and fold distribution of the Teapot Dome area. The diagonal template with an angle between the source and receiver of 45°, and the strike from the perpendicular source of the teapot dome structure is N45°E. The data have a maximum fold of about 65 with a bin size of 110 × 110 ft.

2.1 Common reflection surface processing (CRS)

The CRS is a method to describe the kinematics of reflections in seismic reflection data utilizing an analytic traveltime approximation centered at zero offsets. The method of the CRS stack is theoretically well-established by refs. [15,16,19]. The method considers layers separated by curved reflectors whose reflection comes from a reflecting segment. The traveltime formula depends on the attributes extracted from the seismic reflection data by fitting the operator to the data for all samples. In 3D, a second-order approximation of transmitted and reflected traveltime rays in a seismic system was developed [20]. Following ref. [21], the 3D CRS traveltime formula can be derived from the well-known two-point paraxial travel time approximations. A multi-covering seismic data set is acquired over a set of homogeneous and isotropic layers, with arbitrary velocities separated by smooth interfaces. In 3D seismic systems, the CRS travel time parameters are generally explained with physical meanings, for example, kinematic wave-field attributes or wave front attributes. The CRS method describes the travel time using eight parameters as follows:

(1) t 2 ( m , h ) = ( t 0 + p z ⋅ m ) 2 + m T Am + h T Bh ,

where t ² is zero-offset traveltime, h is the offset, m is the midpoint displacement, p _z is the 1st order in the midpoint, A is the 2nd order in the midpoint, and B is the 2nd order in the offset.

2.2 Velocity versus azimuth (VVAz) inversion

The method of VVAZ inversion is based on the elliptical NMO equation for Transverse Isotropy (TI) media which has been derived by ref. [22]. Following ref. [10], the velocity used for VVAz inversion is based on the residual variant time of isotropic velocities and azimuthally used to estimate fast and slow NMO velocities and their directions. Thus, the velocity variations with azimuth use an elliptical NMO equation for azimuthal data rather than the conventional NMO to invert fast and slow RMS velocities. The fast direction most likely will indicate the fracture direction, while the fast and slow velocities ratio indicates the HTI anisotropy magnitude. Furthermore, along with fast and slow NMO velocity maps, fracture-induced anisotropy orientation and intensity maps can be created. The advantage of traveltime-based methods is that they measure layer properties rather than the interface or boundary properties [10]. In this case, ref. [23] type interval properties can be calculated to estimate the interval anisotropy for each reservoir interval.

Following ref. [22], in the elliptical NMO equation for TI media, the source-receiver distance or offset does not exceed the depth of the reflector; furthermore, the hyperbolic NMO can be approximated by,

(2) t 2 = t 0 2 + x 2 V NMO 2 ( ϕ ) ,

where,

(3) 1 V NMO ( ϕ ) 2 = 1 V slow 2 cos 2 ( ϕ − β s ) + 1 V fast 2 sin 2 ( ϕ − β s ) ,

where t is the total two-way travel time, t ₀ is the zero-offset two-way travel time, and x is the offset. V _fast and V _slow are the fast and slow NMO velocities, respectively. Furthermore, β _s is the azimuth of the slow NMO velocity, while V _N(ϕ) is the NMO velocity as a function of the source-receiver azimuth.

Equation (3) can be written as follows:

(4) 1 V NMO ( ϕ ) 2 = W 11 cos 2 ( ϕ ) + 2 W 12 cos 2 ( ϕ ) sin 2 ( ϕ ) + 2 W 12 sin 2 ( ϕ ) ,

where W ₁₁, W ₁₂, and W ₂₂ are the ellipse coefficients that are related to the slow and fast NMO velocities.

The azimuthal of the fast NMO velocities and the slow NMO velocity could be written as follows,

(5) 1 V fast 2 = 1 2 [ W 11 + W 22 − ( W 11 − W 22 ) 2 + 4 W 12 2 ] ,

(6) 1 V slow 2 = 1 2 [ W 11 + W 22 + ( W 11 − W 22 ) 2 + 4 W 12 2 ] ,

(7) β s = tan − 1 W 11 − W 22 + ( W 11 − W 22 ) 2 + 4 W 12 2 2 W 12 .

The azimuth of the fast velocity is oriented at 90° to the azimuth of the slow velocities. Furthermore, the total travel equation can be written as follows:

(8) t 2 = t 0 2 + x 2 cos 2 ( ϕ ) W 11 + 2 x cos ( ϕ ) sin ( ϕ ) W 12 + x 2 sin 2 ( ϕ ) W 22 .

Equation (8) can be written as,

(9) d = Gm ,

where d is n-dimensional data vector, m is the 6-dimensional model parameter vector, and G is the n-by-4 data kernel as follows,

(10) t 1 2 t 2 2 ⋮ t n 2 = 1 x 11 2 cos 2 ( ϕ 1 ) x 11 2 cos ( ϕ 1 ) sin ( ϕ 1 ) x 11 2 sin 2 ( ϕ 1 ) 1 x 21 2 cos 2 ( ϕ 2 ) x 21 2 cos ( ϕ 2 ) sin ( ϕ 2 ) x 21 2 sin 2 ( ϕ 2 ) ⋮ ⋮ ⋮ ⋮ 1 x n 1 2 cos 2 ( ϕ n ) x n 1 2 cos ( ϕ n ) sin ( ϕ n ) x n 1 2 sin 2 ( ϕ n ) × t 0 2 W 11 W 12 W 22 .

Equation (10) can be used to solve the model parameters and estimate the anisotropic parameters as follows,

(11) m = ( GTG ) − 1 G T d ,

(12) ε = ( V fast − V slow ) 2 V fast ,

where ε is the fracture intensity.

3 Result and discussion

The seismic data processing schemes of this study are shown in Figure 3, in which the ProMAX software is used. First of all, the geometry and header assignment is carried out using two schemes, namely, all azimuths and four azimuth sectors separately. The all azimuth scheme aims to obtain the RMS velocity and kinematic parameter of the CRS used in the four azimuth sectors. Preconditioning is done to improve the stacking quality of both the synthetic and real data. In the synthetic data, preconditioning is performed by direct wave muting, while in the real data, time-frequency domain noise rejection is performed.

Figure 3

The seismic data processing schemes of this study in which the ProMAX software is used.

In the geometry and header assignment of the azimuth sectors of the synthetic and real data, the relatively narrow azimuth seismic data are divided into four azimuth sectors, with each quadrant at a 45° angle. In dividing the azimuth sector, the regional geological information such as the trend of the main structure (inline direction) or the direction fracture orientation from the well data is used as a reference. In the case of synthetic data, the centers of azimuth values are 0°, 45°, 90°, and 135° because the fracture orientation is N90°E. The centers of azimuth values for the teapot dome data are 359°, 44°, 89°, and 134°. This refers to the main orientation structure of the dome structure, which is N134°E.

Furthermore, the four azimuth sectors’ data are combined to produce seismic data composed of offsets and azimuths called COCA (Common Offset Common Azimuth). The next step is using the NMO of the COCA gather as the input for the VVAz inversion. The Hampson-Russel software is used in the VVAz inversion. The result of the VVAz inversion, the fracture direction (V _fast direction), and the intensity ((V _fast − V _slow)/(2V _fast)) will be validated based on the reference model for the synthetic or well data in case of the real data.

Figure 4 shows the fold distribution of the synthetic data before and after being grouped into four-quadrant azimuth sectors (a). The 64 maximum fold is divided into a four-quadrant data group, and the maximum fold in each quadrant is about 15. The regularization of the CDP gather in each azimuth sector is then carried out using the CRS method. Figure 4b shows the comparison of the before and after CRS of each azimuth center. It shows that the CRS could improve better the nominal fold in each quadrant of the azimuth sector.

Figure 4

The fold distribution of the synthetic data before and after being grouped into four-quadrant azimuth sectors (a). The 64 maximum fold is divided into a four-quadrant data group, and the maximum fold in each quadrant is about 15. The result of the CRS application before and after each azimuth center (b). It shows that the CRS works better to improve the nominal fold in each quadrant of the azimuth sector.

Figure 5 shows the fold distribution of the real data before and after being grouped into four-quadrant azimuth sectors (a). Following the result of the synthetic data, the 65 maximum fold is divided into a four-quadrant data group. The maximum fold in each quadrant is about 15 to 20. Figure 5b shows the comparison of the before and after CRS of each azimuth center. It shows that the CRS works better to improve the nominal fold in each quadrant of the azimuth sector.

Figure 5

The fold distribution of the real data before and after being grouped into four-quadrant azimuth sectors (a). The 65 maximum fold is divided into a four-quadrant data group, and the maximum fold in each quadrant is about 15 to 20. The comparison result of before and after CRS of each azimuth center (b). It shows that the CRS works better to improve the nominal fold in each quadrant of the azimuth sector.

The result of the VVAz inversion of the synthetic data is as follows: the fracture direction is N79.8°W and the fracture intensity is about 0.093. Compared to the model in which the fracture direction is N90°W and inversion is N 79.8°W, a fracture intensity of 0.075 and inversion of about 0.093 still seem reasonable since the dominant frequency used to develop the synthetic data is 5 Hz, and we only used four-quadrant azimuth sectors. In addition to the CRS result for the synthetic data, the COCA after regularization still maintains the AVO class 3 anomaly.

In the real data, the FMI log in the Teapot dome field of the Teensleep sandstone formation is used as validation of the VVAz result. Figure 6a shows that the diameter size of the well is related to the fracture intensity. A well with a large diameter indicates high intensity and vice versa. The five wells have relatively the same fracture direction of about N73°W. The VVAz inversion results after interpolation using the CRS method have a fracture direction and relative intensity distribution that matches pretty well to the five references, as shown in Figure 6a. The high intensity well data are located at seismic data locations with high anomalies, such as the 48-x-28 and 25-1-x-14 wells and vice versa as in the other three wells. The VVAz result of the direction of the fracture shows consistent result referring to the five well data based on the previous study [24].

$Figure 6 The VVAz inversion results after interpolation and regularization using the CRS method and distribution of the well location with the FMI log used as the reference (a). The linear regression between the relative fracture intensity of the well data and the VVAz inversion data (b). The result of the fracture direction and relative intensity distribution matches pretty well compared to the five references. The regression results show that a correlation of 0.9 is obtained, indicating that the trend and the VVAz inversion result’s distribution match the reference well.$

Figure 6

The VVAz inversion results after interpolation and regularization using the CRS method and distribution of the well location with the FMI log used as the reference (a). The linear regression between the relative fracture intensity of the well data and the VVAz inversion data (b). The result of the fracture direction and relative intensity distribution matches pretty well compared to the five references. The regression results show that a correlation of 0.9 is obtained, indicating that the trend and the VVAz inversion result’s distribution match the reference well.

In addition, the linear regression between the relative fracture intensity of the well data and the VVAz inversion data is constructed. The result could be seen in Figure 6b. The regression results show that a correlation of 0.9 is obtained, indicating that the trend and the distribution of the VVAz inversion result in a better match with the reference well.

4 Conclusion

Our work has led us to conclude that CRS can be used as a regularization and interpolation method before VVAz inversion for the relatively narrow azimuth land 3D seismic data. The study found that the application of the CRS method could keep the preservation of the AVO and VVAz anomalies. Based on this study’s result, the narrow azimuth land 3D seismic data with a relatively good offset distribution and a nominal high fold is recommended before applying the VVAz inversion method.

Acknowledgments

The authors wish to thank the RMOTC and the U.S. Department of Energy as the data source. The authors also wish to thank the Global Geophysics Group, Faculty of Mining and Petroleum Engineering, and Bandung Institute of Technology for producing this paper.

Funding information: This research was in part supported by the Riset P3MI ITB 2020 grant funded by the Research and Community Services program (LPPM), Institute of Technology, Bandung (ITB), Indonesia.
Author contributions: Lucky Kriski Muhtar and Wahyu Triyoso developed the method and contributed to the writing of the manuscript. Fatkhan Fatkhan helped and supported the discussion. All authors contributed to the preparation of the manuscript.
Conflict of interest: The authors declare that the research was conducted without any commercial or financial relationships that could be construed as a potential conflict of interest.
Data availability statement: The source data are based on the RMOTC and the U.S. Department of Energy. The data are in the public domain.

References

[1] Roberts G, Wombell R, Gray D, Al-Shamsi A, Suwaina O, Ajlani G, et al. Estimation of fracture orientation offshore Abu Dhabi. Expanded Abstracts of the 63rd EAGE Conference and Technical Exhibition, Expanded Abstracts of the 71st SEG Annual Meeting; 2001.10.3997/2214-4609-pdb.15.O-10Search in Google Scholar

[2] Smith R, McGarrity JP. Cracking the fractures-seismic anisotropy in an offshore reservoir. Lead Edge. 2001;20(1):18–27.10.1190/1.1438870Search in Google Scholar

[3] Gray D, Head KJ. Fracture detection in the Manderson field: a 3D SFD case history. Lead Edge. 2000;19(11):1214–21.10.1190/1.1438508Search in Google Scholar

[4] Ruger A. Variation of P-wave reflectivity with offset and azimuth in anisotropic media. Geophysics. 1998;63:935–47.10.1190/1.1826488Search in Google Scholar

[5] Gray D. Fracture detection in unconventional gas plays using 3D-seismic data. Soc Pet Eng. 2008:114902.10.2118/114902-MSSearch in Google Scholar

[6] Wang J, Zheng Y, Perz M. VVAZ vs AVAZ: practical implementation and comparison of two fracture detection methods. SEG/San Antonio 2007 Annual Meeting; 2007. p. 189–93.10.1190/1.2792408Search in Google Scholar

[7] Li D, Jiao W, Yue F, Xiang X, Pan L. Natural fractures and their effect on reservoir reconstruction in lower cambrian shale, southeast chongqing, china. Energy Explor Exploit. 2015;33:769–84.10.1260/0144-5987.33.6.769Search in Google Scholar

[8] Russell B. Seismic reservoir characterization and pre-stack inversion in resource shale plays. Geoscience Technology Workshop of AAPG Search and Discovery; 2014. Article #41467.Search in Google Scholar

[9] Al Dulaijan K, Margrave GF. Azimuthal velocity analysis of 3D seismic for fracture: Altoment-Bluebell field. CREWES Res Rep. 2015;27:1–15.Search in Google Scholar

[10] Al Dulaijan K. Inversion of azimuthal velocity and amplitude variations for seismic anisotropy (PhD Thesis). University of Calgary, Alberta, Canada; 2017. p. 83–114.Search in Google Scholar

[11] Gulunay N, Chambers RE. Unaliasedf – k domain trace interpolation (UFKI). SEG Technical Program Expanded Abstracts; 1996.10.1190/1.1826390Search in Google Scholar

[12] Gulunay N. Seismic trace interpolation in the fourier transform domain. Geophysics. 2003;68:355–69.10.1190/1.1543221Search in Google Scholar

[13] Gelchinsky B. The common reflecting element (CRE) method (nonuniform asymmetric multifold system). Explor Geophys. 1988;19:71–5.10.1071/EG988071Search in Google Scholar

[14] Hubral P, Hocht G, Jager R. Seismic Illumination. Lead Edge. 1999;18:1268–71.10.1190/1.1438196Search in Google Scholar

[15] Mann J, Jager R, Muller T, Hocht G, Hubral P. Common reflection surface stack a real data example. J Appl Geophys. 1999;42:301–18.10.1016/S0926-9851(99)00042-7Search in Google Scholar

[16] Jager R, Mann J, Hocht G, Hubral P. Common reflection surface stack: image and attributes. Geophysics. 2001;66(1):97–109.10.1190/1.1444927Search in Google Scholar

[17] Eisenberg-Klein G, Pruessmann J, Gierse G, Trappe H. Noise reduction in 2D and 3D seismic imaging by the CRS method. Lead Edge. 2008;27(2):258–65.10.1190/1.2840375Search in Google Scholar

[18] Trappe H, Schuenemann E, Pruessmann J, Suckro S, Linzer H, Schachinger M. CRS azimuth gathers and CRS shot gathers for improved fault imaging in Alpine 3D seismic data. Recorder. 2015;40(4):38–42.Search in Google Scholar

[19] Bergler S, Hubral P, Marchetti P, Cristini A, Cardone G. 3D common-reflection-surface stack and kinematic wavefield attributes. Lead Edge. 2002 Oct;21(10):1010–5.10.1190/1.1518438Search in Google Scholar

[20] Bortfeld R. Geometrical ray theory: rays and travel times in seismic systems second-order approximations of the travel times. Geophysics. 1989;54(3):342–9.10.1190/1.1442659Search in Google Scholar

[21] Hubral P, Schleicher J, Tygel M. Three-dimensional paraxial ray properties part 1: basic relations. J Seismic Explor. 1992;1:265–79.Search in Google Scholar

[22] Grechka V, Tsvankin I. 3-D description of normal moveout in anisotropic inhomogeneous media. Geophysics. 1998;63:1079–92.10.1190/1.1444386Search in Google Scholar

[23] Dix CH. Seismic velocities from surface measurements. Geophysics. 1955;20:68–86.10.1190/1.1438126Search in Google Scholar

[24] Schwartz BC. Fracture pattern characterization of the tensleep formation teapot dome Wyoming. Morgantown, West Virginia: West Virginia University; 2006. p. 40–2.Search in Google Scholar

Received: 2021-04-08

Revised: 2021-05-30

Accepted: 2021-06-01

Published Online: 2021-07-29

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2021-0051

Keywords for this article

fracture direction and intensity; CRS; VVAz inversion

Creative Commons

BY 4.0