For the risk exploration of carbonate rocks in the west of China, clarifying the small-scale cavities, fractures and other geological structures plays an important role in reservoir analysis. Thus, performing research on detecting these discontinuous and inhomogeneous geological structures has an important value in evaluating the conditions of oil and gas. In the petroleum industry, there have been many available methods for fault detection and among them the seismic coherence technique provides a solid foundation for automatic identification. The first generation of coherent volume algorithm (Bahoric et al., 1998) is based on the classical normalized crosscorrelation theory but with the requirements of a higher S/N of seismic data. In order to deal with the low S/N data, Marfurt et al., (1998) proposed the second generation of coherent volume algorithm for improving noise immunity that introduces a multi-channel similarity of covariance matrix. However, the averaging effect will reduce the lateral resolution of seismic interpretation. The third generation coherent volume algorithm was developed (Gersztenkorn et al., 1999) by calculating the eigenvalues of covariance matrix. This algorithm can improve seismic interpretation by eliminating noise and increasing lateral resolution. Some other seismic discontinuity detection methods inspirited by the coherent volume technology, also grow up rapidly, such as the seismic attribute method (Li et al., 2005 ; Zhang et al., 2010 ; Li et al., 2011 ; Fang et al., 2013 ; Huo et al., 2014 ; Wang et al., 2014 ; Chen et al., 2014), the spectral decomposition method (Liao et al., 1990 ; Wei, 2009 ; Zhu et al., 2009 ; Chen et al., 2010 ; Chen et al., 2011 ; Huang et al., 2012 ; Yu et al., 2013 ; Cai et al., 2014 ; Zhang et al., 2014 ; Wang et al., 2015), the ant tracking technique (Zhang, 2010 ; Wang et al., 2013 ; Huo et al., 2014). In order to eliminate the scattering and diffraction effect generated from fault zones, some scholars (Wang, 2009) even introduced the chaotic theory for measuring the regularity and chaos properties of seismic amplitudes. Liao et al.(1990) proposed a fault identification method based on the spectral decomposition. Zhang et al.(2006) applied the spectral decomposition technique to study the carbonate reservoir.

In fact, in the stage of analyzing carbonates, in addition to caring about tiny faults and other discontinuities, the inhomogeneity related to the inner conditions of reservoirs, such as cavities should also be paid attention to. Therefore, in this paper, we propose a method based on the L₁ non-smooth norm that can provide a mathematical mode for dealing with both the discontinuous and inhomogeneous information.

The seismic strong reflections from large-scale structures generally will mask the existence of small-scale discontinuous and inhomogeneous geological structures, because the energy difference between them can be one or two orders of magnitude. Thus, the hidden small-scale information can be extracted if we remove the strong reflections. The plane wave destruction filter that uses a local plane wave model to characterize seismic structure has a good potential in predicting linear events and will be adopted for estimating reflections. After subtracting the predicted reflections from seismic imaging data, the residual data will include nonlinear structural features and noises. Then, a L₁ sparse inversion method will be developed for inversing the inhomogeneous information.

The reflections with stable amplitudes are corresponding to large-scale smooth structural characteristics of seismic imaging data. The plane wave destruction method (PWD)(Fomel, 2002) can take a job for estimating reflections and is given by

(1)

where d is a vector of seismic residual data, s = [s₁, s₂, …, s_N, ]^T is the given seismic imaging data, R (σ) is a nonlinear plane wave prediction operator, σ is the local slope of seismic imaging data. The formula (1) can be expressed as

(2)

where I is the identity matrix, σ is the local slope vector, the filter P_{i, j} (σ_i) means a prediction of trace j from trace i according to the plane wave equation and local slopes.

In solving seismic local slopes, the initial local slopes will be given and then the conjugate gradient method is used to update the local slopes. For prediction reflections, the low-order term of the filter P_{i, j} (σ_i) is used because only the smooth information of seismic imaging data is considered, and its form is given as

(3)

where Z_t, Z_x are the Z transform of t, x respectively.

As seismic local slopes corresponding to the locations of faults, fractures, cavities and other geological anomalies cannot be estimated accurately by the PWD method, the discontinuous and inhomogeneous geological information will be left into the seismic residual data d.

In order to extract the useful discontinuous and inhomogeneous information from the seismic residual data, the nonlinear structure-enhancing filtering method (Liu et al., 2014 ; Liu et al., 2010) is required. Here, the operator D (σ) can be a Gaussian similar mean filter (Liu et al., 2009) or lower-upper-middle filter (Hardie et al., 1993). In this paper we apply the lower-upper-middle filter (Liu et al., 2010):

(4)

where, r* is the result of the lower-upper-middle filter, k is the smoothness parameter, l is the sharpening parameter, N is the number of sample point, s* is the intermediate value of sample point.

The magnitude of seismic energy from small-scale discontinuities and inhomogeneities is usually comparable with noises, therefore the inversion method will be considered for accomplishing a higher S/N image. Since the discontinuous and inhomogeneous information has a characteristic of sparsity, a L₁ norm model will be constructed and the minimization problem can be written as

(5)

where m is the inversion model. The above equation is equivalent to solving an unconstrained problem using the Lagrange operator:

(6)

where α is an adjustable regularization parameter, and Ω(m) is defined as (Wang et al., 2011a)

(7)

where ${{h}_{{\tilde{\varepsilon }}}}$ is defined as a piecewise function:

(8)

where $\tilde{\varepsilon }$ is a small constant parameter. Obviously, as $\tilde{\varepsilon }$ → 0, Ω(m) approximates the L₁ norm very well.

Supposing H (_m) is the Hessian matrix of the function Jα(m), quasi-Newton methods can be used to solve this minimization problems. The benefit of the quasi-Newton method lies in only the requirement of the first-order partial derivatives to ensure a fast convergence, which is more appropriate for solving the large-scale computation problem.

The objective function Jα(m) is a nonlinear function and its gradient can be calculated:

(9)

The above equation can be rewritten as

(10)

where the function K (m) is defined as

(11)

where and ${{h}_{{\tilde{\varepsilon }}}}$(m_i) is defined as

(12)

Supposing B ≈ H^-1(m), the gradient of the function Jα(m) at the k-th iteration is defined as g_k = g (m_k), and the direction for updating inversion parameters is defined as s_k =(m) _k+1 -(m) _k. Therefore, the next inversion parameters can be updated by the following equations (Wang et al., 2011b):

(13)

where the matrix B_k is calculated:

(14)

The parameter ω_k is the step size for a line search, which can be solved by the Armijo-Goldstein criterion (Yuan and Sun, 1997).

The iteration steps of the quasi-Newton algorithm can be described as follows:

Step 1: Given initial inversion parameters m₀, a constant parameter r and a symmetric positive definite matrix B₀ = I, and set k = 0.

Step 2: If k 6 ≠ 0, compute m_k and B_k.

Step 3: Compute g (m_k). If , where ξ < 0.1 STOP; otherwise, GOTO

Step 4.

Step 5: Set k = k + 1, and update m_k and B_k; GOTO Step 3.

Remark: Stopping criteria is a key issue for an inversion algorithm. For general inverse problems, the accuracy of a solution should be dictated by the error of input data (Wang, 2007). If we consider the conventional gradient energy value as a stopping criterion, it will fall into many iterations and lead to the accumulation of computation errors. Therefore, we choose is satisfied, then terminate the iteration' as a stopping criterion, where ξ is a user-controlled parameter and ξ < 0.1 is usually required.

In order to test the performance of the proposed seismic sparse detection method for detecting discontinuities and inhomogeneities, a fracture-cavity model is designed that includes breakpoints, cavities and rough reflectors. Seismic diffractions will be generated from this kind of small-scale geological structures or elements and will be imaged after migration. However, compared with the reflections, diffractions are generally much weaker. So it is difficult to observe the weak diffraction information in conventional migration images. The conventional migration result of the fracture-cavity model is shown in Fig. 1, where the large-scale structures are clearly imaged but the small-scale fractures and cavities seem invisible because of their weak energies. Therefore, strong reflections corresponding to large-scale structures must be removed for releasing the high-resolution diffraction details. The estimated reflections by the plane-wave prediction method are shown in Fig. 2. After subtracting the estimated reflections, we get the seismic residual data that concludes the non-smooth information. The sparse inversion result of the seismic residual data is shown in Fig. 3. The four series of cavities in the shallow parts, three series of cavities in the deep part, and discontinuous information of the breakpoints and faults can be clearly seen in Fig. 3. In this sparse inversion process, the regularization parameter α in the objective function plays a key role: if α is chosen larger, the effect of sparsity constraint will become stronger; if α is chosen smaller, it will become a smoother constraint model. The purpose of using L₁ sparsity norm is to protect the discontinuous and inhomogeneous information. Compared with the conventional migration result, the proposed method has advantages in imaging the breakpoints, fractures and cavities. The migration noises also get suppressed in this sparse inversion process. Although there are some linear events left in Fig. 3, the amplitudes along them change dramatically, which in some way demonstrate the sensitivity of the proposed method to amplitude variations. After removing the smooth large-scale geological information, the sedimentary discontinuities, faults, cavities, rough reflection interfaces and other information are separated and enhanced. This kind of small-scale geological information can be used for studying the inner conditions of oil and gas reservoirs.

Fig. 1

Fig. 1 Depth migration of a fracture-cavity model

Fig. 2

Fig. 2 Reflection information predicted by plane-wave method

Fig. 3

Fig. 3 Discontinuous and inhomogeneous geological information extracted by the seismic sparse inversion method

The testing result of the fracture-cavity model shows that the proposed seismic sparse inversion method can clearly reveal small-scale geological structures or elements that are related to the migration channel or accumulation of oil and gas.

The high-capacity carbonate reservoirs in abroad are generally related to homogeneous reefs, complete structures of dolomite weathering crusts and are less involved with heterogeneity properties. However, in China, a typical feature of the majority of effective carbonate reservoirs is strong discontinuity or inhomogeneity, such as cavities or fractures, as they suffered from a long history of tectonic movements and weathering. In this case, detecting karst fractures or cavity groups has a significant value in discovering a large-scale oil and gas reservoir, especially for the risk exploration in the west of China. In this given example, we will try to clarify the distributions of such small-scale discontinuous and inhomogeneous geological information along the Ordovician layer in the Tarim basin.

One inline image of the 3D prestack time migration result is shown in Fig. 4, where the vertical location of weathering crusts is at 3.5 seconds and the vertical position corresponding to bead-shaped cavities is at 4 seconds. Although Fig. 4 can show the changes of amplitudes, the shapes of weathering crusts and even the shadows of cavities in a certain extent, it is not clear enough to describe the small-scale karst cavities in the deep part. We provide the linear prediction result of Fig. 4 by the plane-wave destruction method in Fig. 5 that mainly describes the macro-scale smooth events but not for the small-scale fractures, karst cavities and other geological units. The discontinuous and inhomogeneous geological information extracted by the proposed seismic sparse inversion method is shown in Fig. 6. After removing the smooth and strong reflections, the tiny faults and small-scale cavities in the deep part can be clearly showed and even for the breaking points along the weathering crusts. More important is its ability to depict the development of the karst zones along the Ordovician layers.

Fig. 4

Fig. 4 Three-dimensional pre-stack time migration profile in inline direction

Fig. 5

Fig. 5 Smooth reflection profile obtained by the plane wave destruction method

Fig. 6

Fig. 6 Discontinuous and inhomogeneous geological information extracted by the seismic sparse inversion method

In order to further study the planar distribution of fractures and karst cavities, some seismic attribute analysis is performed along the Ordovician layer. One of the important seismic attributes is the coherent volume that can highlight the non-correlation information, corresponding to faults, fractures and so on, by transforming seismic imaging data into similar coefficients and the result is shown in Fig. 7a. The shape of an ancient river can be clearly shown in Fig. 7a, but the small-scale geological bodies such as karst cavities and tiny faults are difficult to be found. After calculating the root-mean-square (RMS) amplitudes of the seismic sparse inversion data, we present Fig. 7b that can clearly show the planner distribution of the small-scale karst cavities and even the development characteristics of fractures.

Fig. 7

Fig. 7 Seismic attribute maps along the Ordovician layer (a) Coherence attribute;(b) RMS amplitude attribute of seismic sparse inversion data.

The 3D seismic field application in the Tarim basin confirms the ability of the proposed seismic sparse inversion method for detecting the discontinuous structures, such as tiny faults and fractures, and small-scale inhomogeneous elements, such as karst cavities.

By constructing a L₁ sparsity-constraint inversion model, this paper presents a method for detecting the small-scale discontinuous and inhomogeneous geological information. The proposed method demonstrates a good performance in extracting the small-scale geological structures or elements, such as tiny faults and cavities. For practical application of the method, some valuable notes should be considered:

(1) In this paper, the plane-wave destruction filter is used to predict the linear and smooth reflection information. Some other methods, such as F-K filter or median filter, can achieve a similar technical effect.

(2) In order to improve the S/N of the discontinuous and inhomogeneous geological information, this paper constructs a sparsity-constraint model and then uses the piecewise function approximation method and the quasi Newton algorithm to achieve a high-preformation computation. Thus, the proposed method has a capability for processing 3D high-density seismic imaging data.

(3) One application condition of the proposed method is the requirement of seismic imaging data to have the discontinuous and inhomogeneous information imaged. Otherwise, some additional seismic data processing should be performed.

(4) In the future, an issue that how to individually discriminate the small-scale discontinuous or inhomogeneous information, and then only use the extracted small-scale inhomogeneous information to study "beaded" reservoir, can be further explored.

This work is supported by National Key Research and Development Program (2016YFC0501102), National Science and Technology Major Project (2016ZX05066-001), Coal United Project of National Natural Science Foundation (U1261203), Shanxi Natural Science Funds Project (2013012001). We are grateful to the editor and reviewers for kind comments and suggestions which greatly improved the quality of our paper.