A seismic diffraction extraction method for the study of discontinuous geologies using a regularisation algorithm

Key words: diffraction extraction, local slope estimation, plane-wave equation, regularisation method.

In exploration of oil and gas, accurately locating small-scale geological features, such as tiny fault blocks, pinch outs, wedge outs, reef edges or any abrupt change in facies (Kanasewich and Phadke, 1988; Moser and Howard, 2008), can be a challenging issue. Identifying these small-scale subsurface geologies is necessary for prediction of reservoir condition. In the view of seismic wave theory, diffractions including diffracted/scattered waves carry information about these small-scale geologies, and can be used for high resolution imaging.

Using seismic diffractions for studying faults and scatterers has been emphasised in several studies (Krey, 1952; Hagedoorn, 1954; Kunz, 1960; Bansal and Imhof, 2005). Utilising the different characteristics of reflections and diffractions, various methods to separate reflections from diffractions have been developed. Landa et al. (1987) successfully attenuated reflections and enhanced diffractions in the common offset section by employing the double-square-root traveltime moveout equation. As reflections and diffractions have different kinematic characteristics with respect to a point source, plane-wave destruction (PWD) methods have also been investigated for separation of diffraction (Taner et al., 2006; Fomel, 2002; Fomel et al., 2007). Based on different focusing geometries of diffraction and reflection, the focusing-defocusing diffraction imaging approach has been proposed by Khaidukov et al. (2004). Several publications also discuss modifying Kirchhoff migration by a weighting factor to enhance diffractions (Kozlov et al., 2004; Zhang, 2004; Moser and Howard, 2008). In order to increase the signal-to-noise (S/N) ratio of diffractions in separation, multi-focusing (Berkovitch et al., 2009) and common reflection surface (CRS) techniques (Dell and Gajewski, 2011; Asgedom et al., 2011) were studied using super-gather stacking. Furthermore, a windowed or steered version of the multiple signal clarification method for overcoming limited apertures and finite bandwidth of seismic data was developed by Gelius et al. (2013).

In fact, during the process of seismic imaging, diffracted/scattered waves will constructively stack with each other to become a strong event. In a mathematical sense, seismic reflections with a smooth property can be described by a regularisation theory. Therefore, in this paper, we develop an L₂-norm diffraction extraction method by subtracting the predictable reflections. To enhance stability, two regularising constraints for adjusting the smooth property of reflection events are imposed on the minimisation model.

The paper is organised as follows. First, we briefly introduce the method of diffraction extraction. Then a new technique based on regularisation for estimating of reflection dip angle is proposed. Applications of the proposed method are given in both synthetic and field examples.

Identifying and removing reflections using the plane-wave destruction technique has been considered by Claerbout (1992) and Fomel (2002). Generally, the plane-wave destruction filter is defined by the local plane-wave differential equation

where P(t, x) is seismic wave field and σ denotes the local slope varied with time and space. The general solution of the differential Equation 1 is in the form of a plane-wave

with f (x) a waveform function. If the slope σ is irrelevant to time, then the Equation 1 can be transformed into the frequency domain

The solution of the plane-wave differential Equation 3 in the frequency domain is

where is the Fourier transform of P(t, x), the complex part of the exponential in Equation 4 represents t-trace shifting. With the Z transform of Equation 4 and the Taylor expansion of the complex exponential e^iωσ, estimation of the slope σ reduces to solving a nonlinear equation (Fomel, 2002)

where d is a vector of observed seismic data, the operator C(σ) is the convolution of a 2D filter with seismic data and r is the destruction residual. The form of the operator C(σ) and the seismic data d can be written as

where P_i,j(σ_i) means prediction of trace j from trace i according to the plane-wave equation, I is the identity operator, s_i, i =  1, 2, ··· N represent multi-traces. Assume the slope σ is known, Equation 5 becomes a forward problem. Our aim is to invert the slope σ, which is an inverse problem.

Since the matrix operator C(σ) contains the local slope information, the slope information can be obtained by solving a least-squares problem: ‖‖r(σ)‖‖ → min. It is a nonlinear optimisation problem. The nonlinearity comes from the operator C(σ), where the filter comes from a series of Z transforms. Generally, the linear iterative method (e.g. Gauss-Newton algorithm) is the most practical method and the solution can be solved through:

where Cʹ(σ_k) is the differentiation of the filter operator C(σ_k) with respect to the local slope. Therefore, the updated value Δσ_k can be obtained by solving Equation 6, and the slope information can be estimated using iterative methods.

It is evident that the residual r should approach zero if the estimation of the slope is accurate enough. Therefore, instead of Equation 6, we have the approximate equation

The estimation of the slope can be updated using the formula: σ_k+1 = σ_k + Δσ_k.

Solving for the slope σ is an inverse problem. It can be proved that C(σ) is a compact operator. In reality, seismic wavefields may be contaminated with noises and the plane-wave assumption may be invalid where faults and conflicting boundaries exist. Therefore, direct solution of the Equation 5 may be unstable, i.e. ill-posedness holds. To tackle the ill-posed nature, we consider using the Tikhonov regularisation technique.

Note that our aim is to minimise the norm of the residual ‖‖r(σ)‖‖. To ensure stability, we establish the following constrained minimisation model

In Equation 8, α > 0 is the regularisation parameter, Ω(σ) is the stabiliser. We consider the form of the stabiliser Ω(·) as the Sobolev penalty functional: . Using integration by parts, we have (Wang, 2007)

Here, denotes the outward unit normal to the boundary ∂Ω, and L₁ denotes the negative Laplacian:

With homogeneous Neumann boundary conditions, Equation 9 reduces to the following form

where (·,·) denotes the inner product over the L²(Ω ) space, and L₁ is the negative Laplacian.

Considering the complex structure of layers beneath the earth, we further enforce another constraint on the solution, i.e. the difference operator which is in the form

Then, our minimisation problem becomes

With proper choice of the regularisation parameters α and β, the two constraints can ensure the computational values remain in the feasible region of the true solutions, and hence the problem represented by Equation 11 will lead to a stable approximate solution of the desired solution (Wang et al., 2008). The values of α and β are used to balance the fitness of the predicted traces to the observed traces. Best values of parameters α and β should be calculated using the discrepancy principle, however, much more computational effort will be added. Therefore, in this paper, we only consider a priori values of α and β. The values of α and β are set in an a priori manner, i.e. we choose α ≤ β < 1.

The gradient of J^α,β(σ) in Equation 11 can be calculated as

where

We have the following proposition: σ* is a critical point for Equation 12 if and only if g_J^α,β(σ) = 0. Therefore, if σ* is a minimiser of Equation 11, then

for any ξ > 0, where P_Γ(·) denotes a projection onto the feasible set of the slope.

Equation 14 gives us a fixed point iterative method. The formula reads as

Note that in computation, σ_k is a vector to be updated which represents an updated slope field. Therefore, we propose the following iterative scheme:

• Step 1: Input the initial guess value of the slope σ₀; set the iteration index k := 0;

• Step 2: Iterate for k until convergence:

• Step 3: Compute the negative gradient s_k := −g_J^α,β (σ_k )

• Step 4: Line search:

• Step 5: Update the solution: σ_k+1 = P_Γ (σ_k + ξ_ks_k)

• Step 6: Set k := k + 1, GOTO Step 2.

It can be proved that the iteration sequence {σ_k}_k converges to the minimiser of Equation 11 (Wang, 2007). The parameter ξ_k is the step size, which can be obtained by the Wolfe line search criterion, i.e. the parameter ξ_k should satisfy the following two conditions (Yuan, 2008)

where , ω and are two constant numbers requiring . In our numerical experiments, we fix values of ω as 0.1 and as 0.4.

As mentioned above, in order to extract seismic diffractions, we develop a stable reflection slope estimation method by considering the regularisation technique and nonlinear iteration algorithm. For clarity, the whole workflow is illustrated in Figure 1. This diffraction extraction method can be carried out in both the time and depth image, and also applied to prestack or poststack data. To apply our proposed method, the only requirement is that diffractions in the data should not be severely attenuated.

Fig. 1.  The workflow diagram of seismic diffraction extraction. ‘PWD method’ refers to the plane-wave destruction method.

In this section, we generate synthetic seismograms and test our diffraction extracting algorithm on two typical geologic models: a tiny fault-scatterer model and a basin-cavity model.

The first model includes a simple geology with four vertical faults and three diffractors. These faults are located, respectively, at the horizontal positions of 0.6 km, 3.6 km, 6.6 km and 9.6 km, with the scale as 1 wavelength, 1/2 wavelength and 1/4 wavelength and 1/8 wavelength. The diffractors are placed at the horizontal locations of 2.1 km, 5.1 km and 8.1 km, with scale as 1/2 wavelength and 1/4 wavelength and 1/8 wavelength, respectively. The velocities used between the two layers are 2 km/s and 2.6 km/s. Synthetic data were generated by a Kirchhoff forward modelling algorithm with a 20 Hz Ricker wavelet. The synthetic seismic record (Figure 2) shows that there exist diffractions at the edges of vertical faults and scatterers with different scales. The imaging results of the seismic synthetic records are shown in Figure 3, which includes both large-scale and small-scale geological features. However, small-scale features seem invisible because of their weak energies. In order to separate the discontinuous diffraction information, we apply our method to calculate the slopes, and the reflection dip angle field is shown in Figure 4. From the dip field section, we can easily see that faults and diffractors have different features. With this reflected dip angle information, the accurately estimated reflections are displayed in Figure 5. Then following the procedure in the Workflow for extracting seismic diffractions section above, we obtain the diffraction section (Figure 6) by directly subtracting the prediction reflection data (Figure 5) from the seismic image data (Figure 3).

Fig. 2.  Synthetic seismic record of tiny fault-scatterer model.

Fig. 3.  Image result of tiny fault-scatterer model.

Fig. 4.  Reflection dip angle of tiny fault-scatterer model.

Fig. 5.  Reflections estimated by the proposed method.

Fig. 6.  Extracted diffractions of tiny fault-scatterer model.

In real applications, geologists often need to solve different scale problems, like macro-scale layer interpretation, micro-scale fracture locating, fault tracking or reservoir analysis, which exist in all stages of oil-gas exploration and development. Reflections with strong energy have the advantage of revealing macro-scale structures. Diffractions have weak energy but high resolution in enhancing small-scale discontinuous objects. Their resolution can even be smaller than a quarter of the seismic wavelength in horizontal direction and one-sixteenth of the seismic wavelength in vertical direction (Phadke and Kanasewich, 1990). This synthetic experiment illustrates the effectiveness of the proposed method in predicting reflections and extracting diffractions.

The second model is a cavity-like geological structure relevant to carbonate reservoirs. We perform diffraction extraction on prestack depth migration image in this example. Six strings of cavities are successively placed at the shallow parts from CMP 3.7 km to CMP 5.4 km and another two strings are respectively placed at CMP 5.2 km and CMP 5.6 km. This model includes cavities with scales from 5 m to 20 m in diameter. At the bottom, we design a basin-like structure, which is usually the geological exploration object. There are also several tiny faults in the central and boundary of the basin. For this experiment, the main task is to reveal cavities, tiny faults and diffractors, which control the flow property of oil and gas. The synthetic data was generated by a Kirchhoff algorithm and a shot gather is shown in Figure 7.

Fig. 7.  Synthetic seismic record of basin-cavity model.

Using prestack Kirchhoff depth migration to shot gathers, we get the seismic image shown in Figure 8. The imaging result demonstrates large-scale structures, such as the horizontal and dipping layers. However, cavity structures and small-scale scatterers cannot be clearly identified. In seismic exploration, reflections are generally used to interpret large-scale layers and the stratigraphic features of the subsurface. For reservoir research and oil/gas development, the resolution of seismic reflections limited by Rayleigh criterion cannot satisfy the geological requirements. Diffractions with higher resolution but weak energy are generally masked in an invisible way by reflections. Therefore, removing reflections is necessary for diffraction analysis. Using the proposed regularisation method, we can accurately estimate the reflection local slopes in Figure 9. Correspondingly, the reflections can be predicted in Figure 10. Although some cavities still have their shadows in the reflection prediction image, their amplitudes differ from those in prestack depth migration image. Following the procedure in Figure 1, we obtain the diffraction extraction result illustrated in Figure 11. To be clear, we mark the cavities with red circles and faults with black arrows. The diffraction extraction profile (Figure 11) clearly shows edges, cavities and inner details of a complex geological model.

Fig. 8.  Conventional image profile of basin-cavity model.

Fig. 9.  Reflection dip angle of basin-cavity model.

Fig. 10.  Reflections estimated by the proposed method.

Fig. 11.  Extracted seismic diffraction profile of basin-cavity model. Black arrows and the red circle indicate faults and cavities, respectively.

To further test the robustness of our regularisation method, we add random noise into the simulated data. The noisy data and the corresponding conventional imaging result are shown, respectively, in Figures 12 and 13. With our regularisation method, the results of reflection dip angle and extracted diffraction are displayed in Figures 14 and 15 (the marks having the same meanings as those in Figure 11). The simulation results indicate that our method can extract correct diffraction information even with noisy data. This verifies the robustness of our method in extracting seismic diffractions.

Fig. 12.  Synthetic seismic record of basin-cavity model with random noise added.

Fig. 13.  Conventional image section of basin-cavity model (noisy case).

Fig. 14.  Reflection dip angle of basin-cavity model (noisy case).

Fig. 15.  Extracted diffraction profile of basin-cavity model (noisy case). Black arrows and the red circle indicate faults and cavities, respectively.

Small-scale faults and cavities are critical in oil and gas exploration. In the following, we present a field data example to demonstrate the effectiveness of applying the proposed diffraction extraction method to identify and locate distributions of the faults and cavities.

The 3D land seismic data, for studying carbonate reservoirs, is from Western China. The Ordovician layer, displayed between 2.3 s and 2.7 s in vertical direction of Figure 16, proves to be the main oil-gas generating and accumulating object. There are many faults and cavities developed in this layer because of multi-stage tectonic movements and underground rivers in the long history of geological evolution. These discontinuous geologies not only influence the connectivity property but also can provide a storage space. After denoising and migration velocity analysis, we get an inline prestack time migration image by Kirchhoff migration in Figure 16. We emphasise that no attenuation of seismic diffractions is performed during the noise-removal process. The prestack time image can illustrate macro-scale geological layers. However, tracking faults and identifying small-scale cavities are still difficult because of the limited resolution of reflections by Rayleigh criterion. Therefore, in order to improve seismic resolution, diffractions should be extracted to reveal discontinuous geologies. We apply the proposed regularisation method to obtain a reflection dip angle section in Figure 17, where the smooth and discontinuous events behave differently. The smooth events generally imply macro-scale layers with stable amplitudes and the discontinuous events are more likely to indicate faults and cavities. Using the extracted local slopes, we get diffraction profile in Figure 18 definitely displaying small-scale discontinuous objects.

Fig. 16.  Prestack time image.

Fig. 17.  Reflection dip angle estimated by the regularisation method.

Fig. 18.  Seismic diffraction profile obtained by the regularisation method. The red circle indicates karsts, the red square indicates small-scale cavities and black arrows indicate tiny faults.

At the top of the diffraction section (e.g. marked by the red circle), many karsts that have developed in sedimentary layers emerge, with the possibility of becoming storing spaces for petroleum. In the deep layers, many small-scale cavities show themselves in diffraction section (e.g. marked by the red square). Tiny faults can be also traced in the diffraction section shown by black arrows. Also, we can observed a phenomena that cavity often develops with faults. This interesting story seems to be untold in the conventional image (Figure 16) but quite clear in the diffraction image (Figure 18).

From the field data example, we conclude that removing reflections can really reveal higher resolution of diffractions. This field application demonstrates the efficiency of the proposed regularised diffraction extracting method in revealing small-scale cavities and tiny faults.

This paper addresses the challenging issues of how to effectively extract diffractions. From a practical view, we require the seismic image without suppressing diffractions. That is to say, a seismic image with both reflections and diffractions migrated properly, which means that diffractions are not severely attenuated in the denoising process and the imaging kernels are not designed against diffractions. We remark that extraction of the diffraction energy can be performed on unmigrated seismic data or the migrated seismic data. However, for complex subsurface structures, seismic diffractions are generally coupled with each other and distributed in a large region. These features make separation of diffraction hard to perform in unmigrated data.

It should be pointed out that straightforward use of diffraction cannot be fully made without the removal of reflections. To resolve the fine details beyond reflection imaging, we have presented a regularisation method for extracting diffractions. The method incorporates two constraints, the negative Laplacian operator and the difference operator to the solution, and can ensure a stable and unique solution. In the slope computation, we solve a nonlinear minimisation problem via fixed point iteration. Applications of the method in both synthetic and field examples demonstrate the ability of the proposed method in identifying small-scale discontinuous geological features.