Advance Wave Modeling and Diffractions for High-Resolution Subsurface Seismic Imaging

Seismic modeling and Imaging for the small-scale feature in a complex subsurface geology such as salt deposit, fracture reservoir, and Carbonate is not casual because of propagated wave affected by many objects once it hits the geologic structure in the subsurface. The principal goal of newly developed seismic modeling & imaging is to get a subsurface image of structural features with greatest sharpness or resolution. Using model dataset the Sigsbee and Marmousi, we illustrate the accuracy of conventional and advance wave modeling techniques. However, in conventional a Finite difference (FD) algorithm is used to generate the data and in advanced wave modeling, the low-rank (LR) approximation is used to acquire zero-offset configuration data. A field dataset from Malaysian basin is re-processed and imaged using diffraction imaging which shows an enhancement in structural interpretation. Furthermore, the results gained from the proposed modeling and imaging approach significantly enhance the bandwidth of the imaged data. Finally, a frequency spectrum shows a recovery of low-frequency from 0 to 60 Hz which is an optimal resolution of seismic imaging.


Malaysian basin
Malaysian basins are structurally complicated because their assessment through different phases of continental accretion, rifting and mountain building. The geology of the Malay Basin is very old and different that the other part of Sarawak and Sabah basin [1]. Geophysical and geological challenges include, fine sand imagery, often beyond seismic resolution; imaging under gas chimney and under carbonates; diffraction imaging; imaging of the internal architecture of the basement; understand the propagation of waves in efficient media and the associated anisotropy; speed analysis and anisotropy; Gas cloud imaging using complete waveform inversion; and multiple elimination [2]. The focus of research in this part of the world is to enhance the image quality of the faults, fracture and karst using diffraction imaging. Figure 1a shows the geographical location of the area which is near to the Borneo main land and a part of Sarawak Basin. Figure 1b is the cross section taken from Malay Basin subsurface geometry which is an extensional regime and because of the tectonic activity the normal faulting can be seen.

Sarawak basin
Diffraction imaging analyzes were carried out in the Sarawak Basin, located north-west of Borneo, forming the southern boundary of the Oligocene-Recent Basin of the South China Sea; its tectonic evolution has been almost matched by seabed rifting and spread in the marginal basin of the South China Sea [4]. The majority of the Sarawak Basin, including this part of Malaysia, is composed of carbonates, making seismic imagery difficult. The pre-carbonated Oligocene deposit in the lower Miocene and later the terrestrial deposit filled the shelter zone in Cycle I and Cycle II. The cycles are described as follows: Cycles I and II (Upper Eocene to Lower Miocene) have been interpreted as channel sands, lowered clays, and deposited coal. Cycle III (lower-middle Miocene) contains fine limestone shale and sandstone, cycle IV (middle Miocene) is composed of limestone with few mixed clastics, cycle V (middle to upper Miocene) is recognized as limestone, cycle  [3] and (b) a cross section of Malay Basin subsurface structure. Fractured basement located at a depth of 2.0-5.0 km, varying with lateral extension [2]. VI to cycle VIII (upper Miocene to Pleistocene) is composed of open marine and coastal clays and sand, respectively. The progradante sediments from cycle V1 to cycle VIII have gradually stifled the large accumulations of carbonates to the present day. In the center of Luconia, carbonate deposits began at the beginning of the Miocene (cycle III) and increased significantly in the middle and late Miocene cycles IV and V [3]. The structure is dominated by a simple extension phase NS trend faults. These faults affected the deposits during cycles I and II and served as foci for the development of carbonate reefs of cycle IV and V and their accumulation. The dominance of intensive pre-carbonation structuring has therefore resulted in poor seismic image quality.

Diffraction
Diffraction hyperbolic patterns occur frequently in recorded seismic data, particularly in carbonate reservoirs due to abrupt lateral changes in impedance contrast and discontinuity of subterranean layers. However, a very serious doubt to the application of classical theory has claimed that stacked seismic data is not true zero separation data because it is not clear that the results of the stacking of data recorded on a wide range of source-receiver separations is practically close to the results of real zero separation recording, since diffraction amplitudes are concerned [5]. Berryhill also explained the concept and compared the theory of zero separation with source-geophone distance without zero separation, concluding that the diffraction amplitudes at a source-receiver separation different from zero and well known with respect to the wave and geometric propagation paths, a point diffractor gives rise to a hyperbolic diagram on the stacked section. Hyperbola is explained as a symmetrical open curve formed by the intersection of a circular cone with a plane at a smaller angle with its axis than the sides of the cone. Diffraction only be considered as hyperbolic if the recovery is homogeneous, which not a very natural supposition to make is. The curvature of the diffraction hyperbola depends on the speed of the medium, the apex being an indicator of the location of the defect. Diffraction imaging is a dare in seismic processing and adopted by a workflow by means of the common reflection surface by [6][7][8].
In reflection seismology applied to exploration geophysics, such wave propagation phenomena are employed to estimate the properties of the Earth's subsurface reflector. Furthermore, the diffraction phenomenon is also concerned with reflection because of the properties of the subsurface as defined above. The acoustic (seismic) impedance, Z=, Where V is the seismic wave velocity and ρ is density. Although seismic migration is now the one of the primary imaging tools employed in the field, the earliest analogue seismic records took the form of simple single-fold illustration [4,[9][10][11][12]. These recordings were characterized by diffracted energy and random noise, but they nonetheless provided a useful interpretation of the subsoil of the Earth. Later, mechanical migration eliminated the structural deformation of early seismic data, with CMP stacks condensing the amount of random noise as the diffracted energy is preserved. Seismic reflection and diffraction waves are essentially different physical phenomena. Most seismic processors tune and image by improving the seismic reflection data, and do not deliberate the diffracted waves present in the processed data, which carry most of the information regarding minor but important subsurface events [13]. Such small-scale underground events (such as faults, fractures, channels, karsts, and salty body edges) occur as diffracted waves in the seismic data, which can be captured [14,15].
The plane wave destruction filter (PWD) was originally introduced by Claerbout [28] for the characterization of seismic images using the superposition of local plane waves. This PWD filter was based on the plane wave differential equation, after the original plane wave destruction filter with the same approximation proved poor when applied to spatial folding data [16]. On the other hand, the dip frequency filtering approach is applied in the f-k domain. Here we use the Fourier transform to convert time data to the frequency domain, with a filter designed to eliminate reflections based on wave cycles per kilometer. In this paper, we develop a workflow that captures these events on a small scale through separate diffractions based on the regularity and continuity of the slope of the local event that corresponds to the reflection event. We compare the two techniques of dipfrequency filtering and wave-destruction filtering, before integrating the two approaches and performing comparative analyzes on the optimal preservation of diffractions.

Diffraction imaging and resolution
Seismic resolution and Diffraction imaging have a direct relationship to each other, while the frequency of the seismic wave also affects seismic image, a high frequency providing a high resolution image. The relationship between seismic wave frequency and resolution, wavelength, penetration, and hyperbolic diffraction curve is shown in Table 1 [17,18].
The resolution of seismic data is typically around λ/4, which means that high frequency waves illuminate a small object; however, as the increase in frequency will also affect penetration depth, it is usually defined by the purpose of the research. Seismic surveys typically use lower frequencies for subsurface imaging. As a result, seismic images have a slightly lower resolution (λ=4 to λ=8 dependent on data quality) but a deeper penetration (10 km) than those produced by biomedical imaging, requiring a lower depth of penetration (in feet) but higher resolution on the acquired frequency of the data.

Diffraction imaging in depth and time
Although many studies have been conducted on diffractions based on both depth and time domains, the latter does not provide an accurate indication of the discontinuity of the subsurface feature examined in this paper. In the field of depth, the existence of diffraction is a frequent indicator of high complexity and a strongly inhomogeneous trend [19]. As a result, pre-stacking depth imaging methods require more work and computation than those using time imaging [20]. In addition, the success of diffraction identification and isolation in the depth domain depends on the accuracy of the velocity model used to provide the full wave depth image.

Multi-focusing diffraction imaging (MFDI)
MFDI is a novel temporal displacement correction method that effectively describes diffraction events, with the ideal sum of diffracted events and the attenuation of specular reflections allowing the creation of an image mainly comprising diffraction energy. The time correction, which is based on the multiple focus method, depends on two parameters, the emergence angle and the radius of curvature of the diffracted wavefront. The above parameters are calculated from the seismic traces of the pre-stack [19]. The result of the IFM is therefore a highresolution full-azimuth seismic image that comprises optimally stacked diffraction events. The diffraction section contains important data on local heterogeneities and discontinuities in the geology of the subsoil, which can be used to improve horizontal drilling and determine the optimal location of exploration wells.
The multi-focus diffraction imaging makes it possible to: • Identify naturally fractured reservoirs • Avoid unwanted liquid pathways • Cartographic sources of anisotropy of speed • Contour defects and salt body

Finite difference modeling
Finite difference methods (FDM) are extensively used in seismic modeling and migration. In this chapter, a conventional FDM model was used for modeling, with model input being velocity and density values, and the model produces seismic data. FDM are numerical methods for solving differential equations by approximating them with differential equations, in which finite differences approach derivatives. In seismic wave modeling, FD methods are used to propagate the wave in the subsoil. This method has no immersion limitation and produces all events associated with the wave equation such as multiple reflection, head waves and elastic wave equation, anisotropic effects, and wave conversion data [21]. Therefore, the modeling of the F-D wave equations is the ideal way to produce the synthetic seismic data. However, the ultimate goal of migration is to get the image of the real earth using seismic data, which is difficult to test the accuracy of the migration methods with the desired results. In case of seismic inversion, the entry includes the traces and the output of the structural image. For this purpose, the quantities below must be calculated [22]: ∅ wave propagation angle τ two way travel time β angle of incident from source or receiver ∂β/∂x geometrical spreading parameter For each source or receiver, the above quantities must satisfy the equations [22].

Plane-wave destructors (PWD)
Plane wave destruction filter is derived from the local plane wave model for characterizing the seismic data. This filter operates in the time domain (T-X), such as the time distance, and can be extended to the frequency domain. PWD is constructed using an imbedded finite difference scheme for the local plane wave equation as described [24].
For the characterization of several plane waves, it is possible to flow several filters similar to that of the following equation: where Z1, Z2, …, ZN are the zeros of the polynomial. The Taylor series method (assimilating the coefficients of the expansion of the Taylor series to the zero frequency) gives the expression For a three-point centered filter B3 (Zt) For a five-point centered filter B5 (Zt), the derivation of Eqs. (7) and (8) can be found in Fomel [24].
The filter used in the present work is a modified version of filter A (Zt, Zx): This filter avoids the requirement of a polynomial partition. In the case of the three-point filter Eq. (7) and the two-dimensional filter Eq. (7), there are six coefficients consisting of two columns; in each column there are three coefficients and the second column is an inverse copy of the first one. However, the decomposition algorithm is significantly more expensive than the FD. This algorithm is extended to Song's anisotropic wave propagation in 2013 by involving the Eigen function rather than rows and columns of the original extrapolation matrix [25] ( Figure 2).

Slope estimation
Slope estimation is a necessary step in applying the FD plane-wave filters to real data [27], although estimating dissimilar dual slopes, σ 1 and σ 2 , in the available data is complicated than estimating a single slope [24].
The regularization condition should thus be applied to both Δσ 1 and Δσ 2 as follows: The solutions of the above equations depend on the primary values of slop. 1 and slop. 2, which should not be equal, but can be prolonged to the numerical equation with respect to the grid number of the dataset. However, this equation is used here to calculate slopes for the given data set. In the current study, we used a reformed and better version of the plane wave destruction method for seismic diffraction parting based on Claerbout [28].

Diffraction separation methods from seismic full-wave data
One of the best practices and methods of diffraction preservation is the plane wave destruction (PWD) filter initially introduced by Claerbout [28] for the description of seismic images using local plane wave superposition. This PWD filter is based on the plane-wave differential equation, after the original wavedestroying filter plane with the same approximation showed poor performance when applied to spatially aliased data in comparison with other filters. Frequencydistance prediction error [16]. Planar Wave Destruction Filter, which can be considered as a time-distance (TX) analog of the frequency-distance (FX) prediction error filter, derived from a local plane wave model is used to characterize the seismic data [24]. Unfortunately, the first experiments in applying flat-wave destruction for spatially-aliased data interpolation [28] have shown poor results compared to standard frequency distance prediction error (FX) filters [16]. A workflow that uses plane wave destruction for diffraction imaging is shown below. Flat-wave destruction of common shift data may have difficulty in extracting diffractions in regions with complex geological and velocity variations [11] (Figure 3).

2D examples of a complex faulted model: Marmousi
The proposed methodology is tested on a model called Marmousi, created in 1988 by the French Petroleum Institute (IFP) [30]. This model contains 158 horizons with horizontal layers and series of normal faults, which makes it complex especially at the center of the model. This model length is 9.2 km and the depth is 3 km. Figure 4(a) shows the Marmousi model, in which we assume that our task is to imagine the structure of the anticline below and consider it as a reservoir. We use a Ricker wavelet at a point source with a dominant frequency of 40 Hz. The size of the horizontal grid x is 4 m and z is 4 m. Figure 4(b) shows the result of the results of the advanced wave modeling using a low rank approximation. This modeling technique is capable of carrying out different designed studies, but for this we have selected the zero shift data record.
This imaging workflow is a new development for high-resolution imaging, as shown in Figure 3. We tested the workflow by separating the diffraction and the residue and representing it separately. Figure 5(a) shows the diffraction separated using plane wave destruction filtering and Figure 5(b) is the residue after diffraction separation. Figure 6(a) shows the zero shift data migrated with full wave imaging and (b) is the imaged section combining both reflection migration and diffraction migration. This shows an improvement in the resolution especially in the amplitude of faults and small discontinuities that are not resolved by the conventional imaging method.  For quantitative interpretation of the results, Figure 6(c) shows a frequency spectrum for both conventional (red) and diffractive (green) imaging. Higher diffraction response. In addition, higher frequency data of 50-60 Hz is enhanced for high resolution imaging.

Application of method in field data: Sigsbee (Gulf of Mexico)
Most information obtained from high resolution or super-resolution images is the result of careful processing. Consideration of diffraction is necessary in advanced wave modeling and processing because in conventional processing, diffractions are suppressed intentionally or implicitly [32]. In addition, the separation of these diffractions using plane wave destruction filtering could preserve the  diffracted amplitude. To prove the proposed modeling and diffraction migration results, we examined the Sigsbee 2A model to study wave extrapolation in a complex velocity field model containing a sedimentary sequence fragmented by a number of normal faults and overlapping (Figure 7a). In addition, there is a complex salt structure in the model that causes illumination problems as part of the current processing and imaging approach. These imaging problems allowed us to develop appropriate algorithms for better imaging. The Sigsbee 2A model has an absorbent free surface condition and a lower than normal water body reflection. These properties do not generate the effect of free surface multiples and lower than normal internal multiples.
For seismic imaging purposes, a zero shift design was chosen to acquire the zero offset seismic data for simplicity, as shown in Figure 7b. Advanced modeling of low-rank waves is used to obtain non-dispersive seismic. The comparison of two waves modeling the finite difference and the lower row is performed for high resolution imaging purposes. Figure 5 shows an evaluation of the results of conventional Figure 8(a) and advanced Figure 8(b) wave modeling and imaging results that indicate that saline body edges are unresolved. Modeling and imaging.
The plane wave destruction method is used to separate the diffractions from the full-wave reflection data. It is important to rely on the estimated dip of the data and to recognize an accurate determination of the sinking wave as it is an essential parameter for plane wave destruction filtering to separate diffraction and reflection. In this research, we separated the diffractions from the complete wave and tried to imagine these diffractions with a correct velocity model (Figure 9b). Figure 9a is the result of a depth-migrated image using conventional wave modeling results, which are poor compared to the results of advanced wave modeling and diffraction  imaging. Figure 9c shows the frequency spectra of the migrated data. The original offset using classical modeling and full wave migration with Fourier migration in several steps is indicated in green, and advanced waveform modeling with diffraction migration is displayed in purple. Improvements in diffraction imaging, especially for small scale events, are highlighted in the spectrum, allowing us to achieve high resolution images with low frequency data recovery after migration. Low frequency signals from long-lived seismic data are critical values for many areas of exploration seismology and hydrocarbon prediction.

Application of method in field data: Malaysian basin
The angular stack is calculated to measure the reflectivity of a particular angle, the angle stacked data are used to observe the amplitude over offset (AVO) for the Direct Hydrocarbon Indicators (DHI) and the inversion in the oil and gas industry [33,34]. The stacking of angles also applies to the overall combination of interceptions and gradients. These angle stacks typically have a near, mid, and far angle, but an angle stack may have more than three angle stacks with a minimum of 1°. Here, for the diffraction studies, a stack of 3 angles has been realized, as shown in Figure 10. For the comparison on the near and far, the stacking is performed with two angles of 4.5 and 31.5 degrees, as shown in Figure 11. A near angle provides the best diffraction amplitude rather than a far angle stack in this study area.
As the stacking angle rises, there is a loss of amplitude in the data, as can be seen in the lower section ( Figure 10). A point diffraction that can be observed in the data (2150 ms) shows the behavior of the diffraction change; the diffraction flanks are not wide in the stack of nearly staggered angles, while in the offsets, the flanks of the curves are wider. This diffraction would require a greater migration aperture to stack the energy of the top diffraction flanks (Figure 11).  Figure 12 illustrates the diffraction with different range of offset. Theoretically, the sides of the diffraction curve are influenced by the velocity and time/depth of the point diffractor [35]. However, Figure 13 shows that the diffraction also depends on the data shift because the near offset data has less diffraction response than the far-lag data. The reason for this is that the far-lag data covers more offsets and therefore recorded the wider angle data set.

Methodology for diffraction separation for high-resolution imaging: Real field data
Non-migrated offset collection data was on condition that for this project. The initial processing, including the sorting of the common midpoint offset (CMP) collection, was performed to obtain the seismic section of the stack. The following processes were implemented in the sorting process:  1. Window selection around the desired structure with the maximum diffraction response.
2. Select the inline or crossline 3D data for 2D analysis.    6. Stack the data for diffraction analysis in the full stack data set.
7. Estimate the dip components from the data using Eqs. (10) and (11) given above.
8. Remove reflections and preserve diffractions via PWD filtering (Plane-wave Destruction) [24]. Figure 14 shows a 2D line extracted which is an unmigrated seismic from a carbonated field in the Sarawak Basin. It has been carefully treated using preimaging procedures. The diffraction separation method was then extended to preserve the diffractions in the actual data. Figure 15a shows the estimated dip components of the data, which help identify tilt faults and pinches, while Figure 15b shows the corresponding texture [36] obtained by convolving a random number field with the inverse of The destruction by plane wave filters The latter was built  using helical filtering techniques [37,38]. The advantage of displaying the texture is to visualize the characteristics of the local plane in the data with the dip. Figure 16 shows the separate diffraction which is the input for diffraction imaging. This diffraction data is migrated separately and merged with the residual data migration shown in Figure 17. The final result of the migration including the diffraction and reflection data which improved the resolution of the data. Inside the red circle on the left side, a major fault is imaged and can be interpreted, further on the right side in the red circle small scale faults are illuminated after imaging.

Conclusion
We present the significance of advance wave modeling and significance of careful pre-image processing and diffraction imaging in the multifaceted earth like fractured zones, fault edges, small-scale faults and pinch-outs. The effect of the Angles and offsets also be incorporated in the study by concluding that a far stack data has a higher diffraction response than the near stack seismic data. The analysis of different offset gather data shows the observation of the diffraction with the angle of an offset. The results of the far offset data provide the higher diffraction response because of the coverage of the long offset in the data.
Plane-wave destruction filters with an improved finite-difference design have added value in processing for maintaining diffraction. These diffraction data contains the information of the diffracted events which are not recorded in the reflected wave, and further migrate it separately and merge with the reflected data. Finally, considering diffraction in the processing and imaging is very useful for the high-resolution diffraction imaging.
A.P Wan Ismail, Khairul Ariffin, Abdul Halim, Siti Nur Fathiyah, Luluan Lubis. My special appreciation goes to Prof. Sergey Fomel and Luke Decker from The University of Texas Austin & Texas Consortium for Computational Seismology for their valuable advice and guidance over the course of this study.
I would also like to thank my colleagues in Center for Seismic Imaging (CSI) who assisted and supported me in the completion of this project including Dr. Yaser, Dr. Sajid, Dr. Iftikhar, Dr. Maman, Amir, Annur, Liu, Teresa, Hammad, and Irfan and to all CSI members. Many thanks for discussion and sharing of knowledge during the research work.

Author details
Yasir Bashir* and Deva Prasad Ghosh Centre of Excellence in Subsurface Seismic Imaging and Hydrocarbon Prediction (CSI), Geosciences Department, Universiti Teknologi Petronas (UTP), Seri Iskandar, Perak, Malaysia *Address all correspondence to: dryasir.bashir@live.com © 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.