Abstract

Sparsity constraint inverse spectral decomposition (SCISD) is a time-frequency analysis method based on the convolution model, in which minimizing the l1 norm of the time-frequency spectrum of the seismic signal is adopted as a sparsity constraint term. The SCISD method has higher time-frequency resolution and more concentrated time-frequency distribution than the conventional spectral decomposition methods, such as short-time Fourier transformation (STFT), continuous-wavelet transform (CWT) and S-transform. Due to these good features, the SCISD method has gradually been used in low-frequency anomaly detection, horizon identification and random noise reduction for sandstone and shale reservoirs. However, it has not yet been used in carbonate reservoir prediction. The carbonate fractured-vuggy reservoir is the major hydrocarbon reservoir in the Halahatang area of the Tarim Basin, north-west China. If reasonable predictions for the type of multi-cave combinations are not made, it may lead to an incorrect explanation for seismic responses of the multi-cave combinations. Furthermore, it will result in large errors in reserves estimation of the carbonate reservoir. In this paper, the energy and phase spectra of the SCISD are applied to identify the multi-cave combinations in carbonate reservoirs. The examples of physical model data and real seismic data illustrate that the SCISD method can detect the combination types and the number of caves of multi-cave combinations and can provide a favourable basis for the subsequent reservoir prediction and quantitative estimation of the cave-type carbonate reservoir volume.

1. Introduction

Spectral decomposition is a widely used technique in seismic data processing and interpretation. It transforms seismic data from the time domain to the time-frequency domain and utilizes the sensitivities of different frequency data to different geological conditions to describe various geological anomalies. The technique can be used to evaluate stratigraphic thickness variations quantitatively and to depict the discontinuity of a geological anomalous body (e.g. Liu et al2015, Chopra and Marfurt 2016). Moreover, it can also overcome the resolution limitations of seismic data to some extent (Partyka 1999, Han 2013).

The application of spectral decomposition techniques is mainly focused on thin bed prediction, special stratigraphic characteristics analysis and direct hydrocarbon detection, etc Partyka et al (1999) adopted the short-time Fourier transformation (STFT) to produce a series of discrete frequency energy cubes for stratigraphic and structural visualization and reservoir characterization. To investigate internal changes within thin-bed reservoirs, Marfurt and Kirlin (2001) developed a suite of seismic attributes, such as the amplitude and frequency of the spectral peak, from spectral decomposition volumes. Castagna et al (2003) applied the instantaneous spectral analysis to detect the low-frequency shadows associated with hydrocarbons. Portniaguine and Castagna (2004) proposed a method that decomposed a seismic trace by solving an inverse problem with a sparse spike constraint. Sinha et al (2005) proposed a new methodology to detect frequency shadows caused by hydrocarbons and to identify subtle stratigraphic features for reservoir characterization, which transformed the time-scale map into a time-frequency map by taking the Fourier transform of the inverse continuous-wavelet transform (CWT). Wang (2007) used the matching pursuit method to analyse lithology and detect gas reservoirs. Bonar and Sacchi (2010) introduced an inverse spectrum decomposition method adopting complex Ricker wavelets as a wavelet library. Rodriguez et al (2012) presented a noncoherent noise attenuation technique based on a constrained time-frequency transform, and applied it to microseismic data denoising. Puryear et al (2012) proposed a new method called constrained least-squares spectral analysis (CLSSA), which is an inversion-based algorithm for calculating the time-frequency spectrum of seismic data and through which good results are obtained using both synthetic and real seismic data.

Since the pore-cave fracture reservoir is microscopic relative to the whole carbonate reservoir, the time-frequency information would be very helpful for describing the reservoir heterogeneity. Therefore, the spectral decomposition technique is widely used in vuggy carbonate reservoir prediction. Li and Zheng (2008) applied the smoothed Wigner–Ville distribution-based spectral decomposition for carbonate reservoir characterization. Li et al (2011) used the Wigner–Ville distribution-based spectral decomposition to find high-frequency anomalies in a carbonate reservoir, and noted that high-frequency anomaly may be used as a more reliable indicator for a good carbonate reservoir than low-frequency shadows in real applications. Lu and Li (2013) proposed a modified deconvolutive STFT spectrogram to describe the characterization of a cavern carbonate reservoir.

In this paper, we review the theory of sparsity constraint inverse spectral decomposition (SCISD) and apply it to carbonate reservoirs for the identification of multi-cave combinations. First, through different synthetic seismic trace examples, we compare SCISD with the traditional spectral decomposition methods, such as STFT, CWT and S-transform. We show that the time-frequency energy and phase spectra of SCISD have higher time and frequency resolutions and better time-frequency focusing property, which can more clearly reflect the local characteristics of seismic signals. Then, we apply SCISD to physical model data and real seismic data, and complete the identification of multi-cave combinations.

2. Sparsity constraint inverse spectral decomposition

The seismic convolution model (Robinson and Treitel 1980) can be described as follows. A seismic trace s(t) is composed of the convolution of a seismic wavelet w(t) and the reflectivity sequence r(t), with additive random noise n(t):
1
where represents the convolution operator. As shown above, equation (1) does not show signal frequency changes over time. Therefore, equation (1) could be expanded to a more complex expression, in which the seismic record is obtained by the convolution of a series of seismic wavelets with different frequency bands and the corresponding frequency-dependent pseudo reflection coefficient sequence, as well as with additive random noise. It is called the non-stationary seismic convolution model and reads:
2
where wk(t) is the seismic wavelet, the dominant frequency of which corresponds to the subscript k, rk(t) is the corresponding frequency-dependent reflection coefficient, k can be considered to be linearly related to the dominant frequency of the wavelet, and K represents the total number involved in the calculation of wavelets or reflection coefficient vectors (Rodriguez et al2012).
According to linear algebra theory, equation (2) can be rewritten as
3
or
4
where Wk and rk(k[1K]) refer to the convolutional matrix for complex Ricker wavelet wk(t) with central frequency k and their associated reflectivity sequences, respectively. The matrix D denotes the convolution matrix library, and m presents the vector containing all pseudo-reflectivity sequences. Consequently, the seismic trace can be decomposed into different frequency components at a specific time by using a dictionary of complex Ricker wavelets that are uniquely defined by a series of central frequencies. Finally, the frequency-dependent reflection coefficient matrix m is provided, which can be considered to represent the time-frequency spectrum of the seismic record (e.g. Bonar and Sacchi 2010, Bonar et al2010).
Because the number of elements in unknown m is often far more than that in seismic records, equation (4) is an underdetermined problem. To solve equation (4), we turned the ill-posed equation (4) into a well-posed problem by imposing an additional constraint on the solution of the inversion problem (e.g. Yuan and Wang 2013, Yuan et al2015). In this paper, we used a sparseness constraint implemented by an l1 norm to produce the sparseness solution (Bonar and Sacchi 2010, Bonar et al2010). The resulting objective function is given by
5

The first term inside the square brackets in equation (5) is the misfit, which matches the solutions with the seismic data. The second term inside the square brackets in equation (5) is the constraint, which shapes the resulting distribution of reflection coefficients. λ is a regularization parameter that controls the sparsity desired in the solution and the stability of the solution.

There are many algorithms for solving an l1 norm problem. In this paper, we use the fast iterative shrinkage-threshold algorithm (FISTA) to solve the objective function equation (5) (Beck and Teboulle 2009). The analytical solution of equation (5) is obtained by
6
where α is a constant that should be no less than the maximum eigenvalue of DHD and yi is an update of the model estimate mi reducing computational time (Bonar and Sacchi 2010, Bonar et al2010).

3. Synthetic data examples

We generated a synthetic trace to analyse the differences among the STFT, S-transform, CWT and SCISD. Figure 1(a) shows a synthetic seismic trace which contains several seismic wavelets with different phases and different dominant frequencies. The dominant frequency of the first Ricker wavelet is 10 Hz with 0° phase. The dominant frequency of the second is 30 Hz with  -90° phase. The dominant frequency of the third is 45 Hz with 180° phase. The dominant frequencies of the last two wavelets are both 60 Hz with a 20 ms two-way travel time. Figures 1(b)(e) show the results obtained by using STFT, S-transform, CWT and SCISD, respectively. It can be observed that the STFT magnitude spectrum has a high time resolution while it loses frequency resolution due to the window effect. The resulting decomposition of CWT suffers from resolution limitations; it has high frequency resolution at low frequencies and high time resolution at high frequencies. The S-transform has a character similar to that of CWT. Compared with these three conventional spectral decomposition methods, the resulting spectrogram obtained from SCISD has better focus-ability and higher time and frequency resolutions. Figures 2(a)(d) are the phase spectra from STFT, S-transform, CWT and SCISD, respectively. It is obvious that the phase spectra of STFT, S-transform and CWT are complex, and it is therefore difficult to extract the desired seismic information directly. Nevertheless, the phase spectrum of SCISD clearly shows the distribution of seismic time-frequency energy and the phase variation of the seismic signal.

A comparison of results for STFT, S-transform, CWT and SCISD applied to a synthetic trace. (a) The synthetic seismic trace, (b) the STFT spectrogram, (c) the S-transform spectrogram, (d) the CWT spectrogram, and (e) the SCISD spectrogram.
Figure 1.

A comparison of results for STFT, S-transform, CWT and SCISD applied to a synthetic trace. (a) The synthetic seismic trace, (b) the STFT spectrogram, (c) the S-transform spectrogram, (d) the CWT spectrogram, and (e) the SCISD spectrogram.

A comparison of phase spectra for STFT (a), S-transform (b), CWT (c) and SCISD (d) applied to a synthetic trace.
Figure 2.

A comparison of phase spectra for STFT (a), S-transform (b), CWT (c) and SCISD (d) applied to a synthetic trace.

We constructed a layered blocky impedance synthetic model that is obtained by the convolution of a 30 Hz Ricker wavelet with a series of positive reflectivity dipoles, and it is adopted to demonstrate the ability of the SCISD method to identify the thin layers. The time thickness of the layers varies from 2 ms in the shallow part to 34 ms in the deep part (shown in figure 3(b)) with a time distance between layers of 100 ms. We applied STFT, S-transform, CWT and SCISD to the synthetic trace (figure 3(a)). The results (figures 3(c)(f)) show that it is hard to identify the event of the thin layer in the frequency domain for STFT (figure 3(c)). In the S-transform spectrogram (figure 3(d)), low-frequency streaking phenomena can be observed due to the interferences. The interferences exist not only between the boundary surfaces of individual thin layers but also between nearby layers. In this example, the CWT spectrum (figure 3(e)) displays low time resolution at low-frequency components because it suffers from the use of wavelets to analyse the low frequency components. In contrast, the result of the SCISD (figure 3(f)) has high time resolution, which resolves layers as thin as approximately 14 ms (the red arrow), and shows a relatively good ability to distinguish thin layers. Figure 4 displays the phase spectra of STFT, S-transform, CWT and SCISD, which correspond to figures 3(c)(f), respectively. Some phase changes are observed in the phase spectra of STFT, S-transform and CWT, but it is difficult to determine which specific phase changes correspond to which seismic signals; it is also hard to identify thin layers from the phase spectra. The phase variation with time can be seen in the phase spectrum of SCISD, and it can be observed that the phase changes from negative to positive and from negative to positive again at 400 ms (in which time thickness is 14 ms), which correspond to two acoustic impedance interfaces.

A comparison of the abilities of STFT, S-transform, CWT and SCISD to identify thin layers. (a) The synthetic seismic trace, (b) the reflection coefficient sequence with increasing thickness as a function of time, (c) STFT spectrogram, (d) S-transform spectrogram, (e) CWT spectrogram, and (f) SCISD spectrogram.
Figure 3.

A comparison of the abilities of STFT, S-transform, CWT and SCISD to identify thin layers. (a) The synthetic seismic trace, (b) the reflection coefficient sequence with increasing thickness as a function of time, (c) STFT spectrogram, (d) S-transform spectrogram, (e) CWT spectrogram, and (f) SCISD spectrogram.

The phase spectra corresponding to figures 3(c)–(f) (a) STFT phase spectrogram, (b) S-transform phase spectrogram, (c) CWT phase spectrogram, and (d) SCISD phase spectrogram.
Figure 4.

The phase spectra corresponding to figures 3(c)(f) (a) STFT phase spectrogram, (b) S-transform phase spectrogram, (c) CWT phase spectrogram, and (d) SCISD phase spectrogram.

4. The identification of multi-cave combinations of carbonate reservoirs

The carbonate fractured-vuggy reservoir is the main hydrocarbon reservoir in the Halahatang area of the Tarim Basin, north-west China. In this area, the carbonate reservoir is deeply buried and strongly heterogeneous. Seismic wave scattering is very strong, and the inter-layer reflection of the deep formation is weak. In addition, due to complex surface conditions and subsurface structure of this region, deep seismic data have low resolution and a low signal-to-noise (S/N) ratio. Because of the acoustic impedance difference between the reservoir and the surrounding rock, the carbonate caves on the seismic profile (figure 5) display anomalous seismic-amplitude bright spots called the ‘string of beads response’ (SBR) (Xu et al2014), which is a significant seismic phenomenon that distinguishes carbonate pore-cavity reservoirs. These issues present serious challenges for the identification and quantitative prediction of carbonate reservoirs.

A real seismic section containing carbonate caves.
Figure 5.

A real seismic section containing carbonate caves.

4.1. Seismic physical model data examples

According to the carbonate reservoir geology of the Halahatang area (figure 5) and the similarity theorem, we designed a complex 3D physical model (figure 6) that is proportional to the realistic geological model. The ratio of similitude is
where vm is the velocity of model material, v is actual formation velocity, lm is the size of model, l is the spatial scale of the actual formation, fm is the dominant frequency of the seismic source in the physical model experiments, and f is the dominant frequency of real seismic data. For instance, 1 mm in the physical model corresponds to 20 m at field scale. The velocity of the model material is 1480 m s-1, which corresponds to 2960 m s-1 at field scale. The dominant frequency of the ultrasonic transducer is approximately 250 kHz, which corresponds to 25 Hz at field scale.
Vertical profile sketch of the physical model (the left part is north, and the right part is south).
Figure 6.

Vertical profile sketch of the physical model (the left part is north, and the right part is south).

The major features of the model are as follows (figure 6): (1) The top formations above the target layer are simplified to horizontal formations. (2) Three nearby overlying strata above the target layer gradually thin out from the south to the north. (3) The fracture-cave reservoir is distributed in the formation between 6300 m and 7000 m. (4) A total of 14 zones are set on a plane, where different types of fracture cave reservoirs are placed.

The 3D seismic data for this physical model are acquired using a seismic physical model experimental simulator. The ultrasonic transducers are used as ultrasonic transmitters and receivers in the seismic data acquisition system. The parameters of the 3D seismic geometry are given in table 1. Conventional processing including observation system loading, trace editing, true amplitude recovery, filtering, multiple suppression, deconvolution, acquisition footprints suppressing and prestack time migration, etc is performed for the collected data. In this paper, we choose two 2D lines which contain different cave-combination types to analyse the identification ability of SCISD for the multi-cave combinations.

Table 1.

Parameters of seismic observation system.

Physical modelGeological model
Observation system12L-30S-240R12L-30S-240R
Longitudinal observation system317.5-20-2.5-20-317.5 mm6350-400-50-400-6350 m
Line spacing25 mm500 m
Perpendicular offset25 mm500 m
Trace interval2.5 mm50 m
Shot interval2.5 mm50 m
Minimum off-inline distance1.25 mm25 m
Maximum off-inline distance173.75 mm3475 m
Bin dimension1.25 mm  ×  1.25 mm25 m  ×  25 m
Folds12 (longitudinal)  ×  6 (horizontal)12 (longitudinal)  ×  6 (horizontal)
Azimuth angle0.550.55
Physical modelGeological model
Observation system12L-30S-240R12L-30S-240R
Longitudinal observation system317.5-20-2.5-20-317.5 mm6350-400-50-400-6350 m
Line spacing25 mm500 m
Perpendicular offset25 mm500 m
Trace interval2.5 mm50 m
Shot interval2.5 mm50 m
Minimum off-inline distance1.25 mm25 m
Maximum off-inline distance173.75 mm3475 m
Bin dimension1.25 mm  ×  1.25 mm25 m  ×  25 m
Folds12 (longitudinal)  ×  6 (horizontal)12 (longitudinal)  ×  6 (horizontal)
Azimuth angle0.550.55
Table 1.

Parameters of seismic observation system.

Physical modelGeological model
Observation system12L-30S-240R12L-30S-240R
Longitudinal observation system317.5-20-2.5-20-317.5 mm6350-400-50-400-6350 m
Line spacing25 mm500 m
Perpendicular offset25 mm500 m
Trace interval2.5 mm50 m
Shot interval2.5 mm50 m
Minimum off-inline distance1.25 mm25 m
Maximum off-inline distance173.75 mm3475 m
Bin dimension1.25 mm  ×  1.25 mm25 m  ×  25 m
Folds12 (longitudinal)  ×  6 (horizontal)12 (longitudinal)  ×  6 (horizontal)
Azimuth angle0.550.55
Physical modelGeological model
Observation system12L-30S-240R12L-30S-240R
Longitudinal observation system317.5-20-2.5-20-317.5 mm6350-400-50-400-6350 m
Line spacing25 mm500 m
Perpendicular offset25 mm500 m
Trace interval2.5 mm50 m
Shot interval2.5 mm50 m
Minimum off-inline distance1.25 mm25 m
Maximum off-inline distance173.75 mm3475 m
Bin dimension1.25 mm  ×  1.25 mm25 m  ×  25 m
Folds12 (longitudinal)  ×  6 (horizontal)12 (longitudinal)  ×  6 (horizontal)
Azimuth angle0.550.55

Different multi-cave combinations can result in some special SBRs. Figure 7 shows the physical model modules photo (upper part) and the corresponding seismic response (lower part) of various multi-cave combinations. The first spherical cave leads to long SBR, the combination of two vertical ellipsoidal caves is characterized by SBR-like chops, the combination of four spherical caves causes the waved SBR in the seismic profile, and a single ellipsoidal cave results in a short-wide SBR. Although the types of SBRs are different, it is difficult to determine seismic reflection from a single cave or multiple-cave combinations depending on SBRs. In addition, the information from the SBRs is not sufficient to describe the details of the multi-cave combinations, such as the number and combination types of the caves. In the actual exploration and production, these will give rise to difficulty in describing the inner structure and the distribution of the carbonate reservoir and will influence the prediction accuracy of recoverable reserves.

The physical model photo (upper part) and the corresponding seismic response (lower part) of multi-cave combinations.
Figure 7.

The physical model photo (upper part) and the corresponding seismic response (lower part) of multi-cave combinations.

From the synthetic seismic data example, it can be recognized that SCISD is a method for providing high time and frequency resolution and effective phase information. In this section, SCISD is adopted to identify multi-cave combinations in a carbonate reservoir. Figure 8 shows the amplitude spectrum of the target formation in the physical model data. It can be observed that the dominant frequency of the target formation is approximately 20 Hz. The spectral components of two 2D lines from 10 to 30 Hz are calculated. In the next section, to illustrate the effect of the recognition of multi-cave combinations, only the spectral component with 20 Hz is displayed.

The amplitude spectrum of the physical model data from the target area.
Figure 8.

The amplitude spectrum of the physical model data from the target area.

4.1.1. Regular cave combinations.

Figure 9(a) displays photos of multi-cave combination modules placed in the physical model. The combination forms of the caves in figure 9(a) are a regular triangle, an asymmetric diamond, a left triangle, a rhombus, two near-vertical caves, a square, two vertical caves and two horizontal caves from left to right. The diameter of each single spherical cave at field scale is 60 m, and the P-wave velocity is 2512 m s-1. Figure 9(b) displays seismic responses of multi-cave combinations corresponding to figure 9(a). The triangle and diamond combinations show the tilted wavy SBRs in the seismic section, and the tilt of the SBRs is related to the position of the caves. Seismic responses of vertical double-cave combinations are long SBRs in the seismic section, and long SBRs will gradually be separated into two short SBRs with an increase in the two-cave spacing. The fifth SBR is inclined, which is created by two nearly vertical caves (i.e. two caves have a little tilt). The seismic response of two horizontal caves is a chop-shaped SBR. The seismic response of a single cave is a short SBR, and two caves are not sufficiently separated from the point of view of seismic exploration due to the small distance between these caves. The result is two SBRs stacked together. From figure 9(b) it can be recognized that different combinations can form different seismic reflection characteristics, and the spatial occurrence of the pore-cave reservoir can be approximately determined by the observed seismic reflection characteristics. However, it is difficult to make a correct and reasonable judgement for distribution details of the underground caves only by observing seismic reflection characteristics. Figures 9(c)(f) show STFT with a 40 ms Hanning window, S-transform, CWT and SCISD isofrequency sections at 20 Hz, respectively. CWT shows the outline of multi-cave combinations. It is obviously not enough to describe the details of the combinations (figure 9(e)). Like CWT, S-transform just reflects the whole structure of the combinations and has poor time resolution in the low-frequency area (figure 9(d)). Compared with the S-transform and CWT, STFT yields better results for identifying multi-cave combinations and higher time resolution, and the energy is also relatively concentrated, all of which are achieved at the cost of frequency resolution (figure 9(c)). In general, STFT, CWT and S-transform can only show the contour of multi-cave combinations but cannot be used to identify the combination types and the number of caves. In addition, the seismic response of the vuggy reservoir is obviously stretched in the vertical direction in STFT, S-transform and CWT isofrequency sections, which makes the apparent height4

4

The apparent height is defined as the distance between the center of the upper cave and the center of the bottom cave.

of multi-cave combinations appear to be higher than the true height. If these attributes are adopted to estimate the vuggy reservoir volume, it will give rise to a considerable error. Compared with these three conventional spectral decomposition methods, SCISD (figure 9(f)) can directly determine the combination form and the number of the caves, and the energy is also more concentrated. The apparent heights of cave combinations are closer to the real height (see table 2, which corresponds to seven combinations in figure 9(f) from left to right). Therefore, it can be concluded that SCISD is a very useful tool for vuggy reservoir identification and height estimation.

The physical model photo (a), the corresponding seismic response (b) and 20 Hz isofrequency sections are obtained by using STFT (c), S-transform (d), CWT (e), and SCISD (f).
Figure 9.

The physical model photo (a), the corresponding seismic response (b) and 20 Hz isofrequency sections are obtained by using STFT (c), S-transform (d), CWT (e), and SCISD (f).

Table 2.

Comparison between apparent height and true height of multi-cave combinations.

1234567
Apparent height (m)134.9173.9137.9155.8119.9143.9203.8
True height (m)120170140150100130200
Relative error0.12420.02290.01500.03870.19900.10690.0190
1234567
Apparent height (m)134.9173.9137.9155.8119.9143.9203.8
True height (m)120170140150100130200
Relative error0.12420.02290.01500.03870.19900.10690.0190
Table 2.

Comparison between apparent height and true height of multi-cave combinations.

1234567
Apparent height (m)134.9173.9137.9155.8119.9143.9203.8
True height (m)120170140150100130200
Relative error0.12420.02290.01500.03870.19900.10690.0190
1234567
Apparent height (m)134.9173.9137.9155.8119.9143.9203.8
True height (m)120170140150100130200
Relative error0.12420.02290.01500.03870.19900.10690.0190

Figures 10(a)(d) display the phase spectra of STFT, S-transform, CWT and SCISD, which correspond to figures 9(c)(f), respectively. It is difficult to obtain effective information about the cave reservoirs from STFT, S-transform and CWT phase spectra. Nevertheless, the phase spectrum of SCISD clearly shows the time-frequency energy distribution of cave combinations. It also contains the information of local phase and the phase change of cave combinations in the horizontal and vertical direction obviously, which is very beneficial for precise interpretation of multi-cave combinations. In the time-frequency energy spectrogram of SCISD (figure 9(f)), the energy of the upper and lower caves of the second combination and the lower cave of the fourth combination is weak, which will cause some trouble for multi-cave combination explanation, but in the phase spectrogram (figure 10(d)), their energy is strong and the distribution details of caves are clear; additionally, the phase change in the longitudinal direction is highlighted. The energy responses of the third combination are mixed together in the time-frequency energy spectrum (figure 9(f)); by contrast, in the time-frequency phase spectrum, the cave in the upper part and two caves in the lower part are clearly separated, and the phase near the boundary changes obviously from negative to positive.

The phase spectra of STFT (a), S-transform (b), CWT (c) and SCISD (d), which correspond to the figures 9(c)–(f), respectively.
Figure 10.

The phase spectra of STFT (a), S-transform (b), CWT (c) and SCISD (d), which correspond to the figures 9(c)(f), respectively.

In the time-frequency phase spectrum of SCISD, it can be seen that the phase change of each cave in the combinations is obvious. These phase changes can more clearly show the external outline and longitudinal distribution characteristics of the combinations, which is conducive to determining the number of caves in the longitudinal direction and also provides a further basis for the estimation of the apparent heights of cave combinations. When the energy of some caves in combination is relatively weak in the time-frequency energy spectrum, the high-resolution phase information displayed on the phase spectrum of SCISD can be used as an important supplement to the time-frequency energy spectrum. Both complement each other, which allows us to infer the distribution of multi-cave combinations and estimate the apparent height.

4.1.2. The columnar caves and layered vertical combinations.

The sketch map of the columnar caves (the lateral five caves on the left) and layered vertical combinations (the lateral three combinations on the right) is illustrated in figure 11(a). The diameter of all columnar caves is equal to 80 m, and the height changes from 25 m to 400 m. The layered vertical combinations consist of multiple small caves in the vertical direction, which are evenly spaced and stacked over a depth range of 200 m. The height of each small columnar cave in the layered vertical combination is 10 m. The first layered vertical combination contains two caves, and the distance between the two adjacent caves is 180 m. The second one contains three caves, and the distance between any two adjacent caves is 85 m. The third one contains five caves, and the distance between any two adjacent caves is 35 m. Figure 11(b) shows seismic responses of the columnar caves and layered vertical combinations. For the layered vertical combinations, the number of SBRs in seismic section increases with the number of caves in the combination. When the spacing between two adjacent caves is larger, the SBR of each cave is separated and forms two short SBRs, whose characteristics are all ‘one peak and one trough’ (The black represents peak and the red represents trough). When the distance between two adjacent caves is 85 m, the seismic responses of each cave mix together and form one long SBR, the energy of which is strong. When the distance between two adjacent caves is 35 m, the energy of the long SBR is significantly weaker because of the destructive interference of the diffracted waves. For columnar caves, when the height of the caves is less than 100 m, SBRs are similar and display ‘one peak, one trough and one peak’; additionally, the top and bottom reflections of the caves are mixed together and are therefore hard to distinguish due to the resolution limitations of seismic data. When the height of the caves is more than 100 m, the top and bottom reflections of the caves are separated, and the top and bottom interfaces can be easily identified. In the same height range (200 m), whether it is a single cave or a multi-cave vertical combination, seismic responses displays all long SBRs.

The sketch map, seismic section, 20 Hz isofrequency section and phase spectrum of the columnar caves and layered vertical combinations. (a) The sketch map, and (b) seismic section of the columnar caves and layered vertical combinations, (c) the 20 Hz SCISD isofrequency section, and (d) the phase spectrum corresponding to (c).
Figure 11.

The sketch map, seismic section, 20 Hz isofrequency section and phase spectrum of the columnar caves and layered vertical combinations. (a) The sketch map, and (b) seismic section of the columnar caves and layered vertical combinations, (c) the 20 Hz SCISD isofrequency section, and (d) the phase spectrum corresponding to (c).

Figures 11(c) and (d) are the 20 Hz SCISD isofrequency profile and phase spectrum profile, respectively, which correspond to figure 11(b). As you can see in figures 11(c) and (d), when the height of a single columnar cave is higher than or equal to 100 m, the top and bottom of caves are obviously separated in both the time-frequency energy spectrum and the time-frequency phase spectrum. The phase on top of the caves change from positive to negative; the phase on bottom of the caves change from negative to positive. However, the effective spectrum information in the middle part of the caves is lost due to large cave height. The 25 m or 50 m columnar cave gives a strong energy group in the time-frequency energy spectrum, and the phase changes from positive to negative to positive again in the time-frequency phase spectrum. By comparing these with their phase changes on top and bottom of the large scale columnar caves, we have made a joint analysis and inferred that the local phase mutation indicates the top and bottom interface of the small scale columnar cave. For the layered vertical combinations, when the spacing between two adjacent caves is 180 m or 85 m, the single cave can be clearly identified in the time-frequency energy spectrum and phase spectrum. When the distance between two adjacent caves is 35 m, the single cave of combination has an obvious phase change, but it is difficult to be recognized in the time-frequency energy spectrum. In the same height range (200 m), layered vertical combinations of two caves in the time-frequency energy spectrum show two energy groups, and so does the single columnar cave, which is because the energy in the middle part of the single cave is lost. To distinguish whether the two energy groups represent a single cave or a two-cave combination, we also need to analyse the phase spectrum. When it is a single cave, the phase changes of two energy groups are from positive to negative and from negative to positive. When it is a two-cave vertical combination, the phase changes of both energy groups are from positive to negative. Combining the energy spectrum and the phase spectrum therefore provides a more accurate portrayal of the SBRs of the vuggy carbonate reservoir and also provides a fairly exact description of the carbonate reservoir distribution and a quantitative estimation of the reservoir volume.

4.2. Field data example

The field data are from the Halahatang area of Tarim Basin, north-west China, in which the Upper Ordovician marine carbonate reservoir is the main hydrocarbon reservoir. In this area, the target formation is buried at a depth of more than 5500 m. Because of the complex surface conditions and subsurface structure in this region, the deep seismic data have low resolution and a low S/N ratio, yielding poor seismic imaging results.

The effect of STFT is better than that of S-transform and CWT, as shown in the physical model data example, so we chose STFT to compare with SCISD. The field data have a dominant frequency of approximately 20 Hz. Similarly, the spectral components of the field data from 10 to 30 Hz are calculated, and only the 20 Hz isofrequency section is displayed in order to illustrate the effect of the recognition of multi-cave combinations. Figure 12(a) is a typical seismic section in the target area, in which we can see some cave-related SBRs, such as long SBR (number 1), waved SBR (number 2), and short SBR (number 3) (Numbers have been marked). Figures 12(b) and (c) present the 20 Hz time-frequency energy spectrum sections, obtained by STFT and SCISD, respectively. Figures 12(d) and (e) show the phase spectra, corresponding to figures 12(b) and (c), respectively. In general, higher resolution and better S/N ratio are provided by SCISD. Compared with STFT, SCISD delineates the locations and extents of carbonate caves more clearly. Although there are some changes in the phase spectrum of STFT, it is still difficult to extract the effective phase information in order to identify the carbonate caves. From the details, long SBR (number 1) in the 20 Hz time-frequency energy spectrum section of SCISD shows two energy groups, but the phase changes from positive to negative three times. We find that these characteristics are similar to those of the three-cave vertical combination in figure 11. Combing the time-frequency spectra of SCISD with the characteristics of the multi-cave combinations in the physical model data, we infer that long SBR (number 1) is a joint response of three caves that are vertically distributed. The waved SBR (number 2) is similar to the seismic response of the third multi-cave combination in figure 9(a). It shows two energy groups in the energy and phase spectra of SCISD, the whole of which is inclined and divided into an upper group and a lower group. We also find the same time-frequency characteristics in the physical model example (the third multi-cave combination in figures 9(f) and 10(d)). Combining the time-frequency energy spectrum and phase spectrum of SCISD with the physical model data example (figures 9(a), (b) and (f) and 10(d)), we make a comprehensive analysis and conclude that it is a combination of three caves similar to the left triangle. The short SBR (number 3) is considered to be a single cave according to the results of SCISD and the characteristics of the single cave in the physical model data, and its vertical dimension is much closer to the actual situation than it is with STFT. The phase changes from positive to negative to positive.

20 Hz isofrequency sections and phase spectra of a real seismic section. (a) The real seismic section, (b) the 20 Hz STFT isofrequency section, (c) the 20 Hz SCISD isofrequency section, (d) the phase spectrum corresponding to (b), (e) the phase spectrum corresponding to (c).
Figure 12.

20 Hz isofrequency sections and phase spectra of a real seismic section. (a) The real seismic section, (b) the 20 Hz STFT isofrequency section, (c) the 20 Hz SCISD isofrequency section, (d) the phase spectrum corresponding to (b), (e) the phase spectrum corresponding to (c).

5. Discussion

Our modelling and field data examples show noteworthy improvement in the identification of multi-cave combinations in the carbonate reservoir that results from the application of SCISD. By adding a sparsity constraint to the algorithm, the results of SCISD are sparse and display good time-frequency focusing, and the energy is more concentrated. In the synthetic data cases, the spectra of SCISD observed have narrower bandwidths than STFT, S-transform and CWT spectra due to the absence of the window effect in SCISD method, which also results in most of the energy of the spectrum being included in a narrow frequency band. In other words, the energy is concentrated within the effective frequency band (In this article, the effective bandwidth of the physical model data and the field data are all 10–30 Hz). Thus, to illustrate the effect of SCISD on the identification of multi-cave combinations, the 20 Hz spectral component was chosen to characterize the internal structures of the carbonate cave complexes.

For the interpretation of multi-cave combinations, the results of applying SCISD to physical model data validate that the SCISD method is more effective at describing the external contour and the internal structure of multi-cave combinations than conventional spectral decomposition methods. Not only can the energy spectrum of the SCISD method reflect more detail of the internal structure of the multi-cave combinations, but the phase spectrum also has the same high-resolution level with the time-frequency energy spectrum. This method also has relatively good results for the field data. Although there is no well data for the selected field data to support the inference on the combination forms of the multiple caves, we can find similar time-frequency features in the physical model data examples. Thus, the SCISD method and the interpretation achievement of physical model data are combined as the basis for our inference. In such situations, our speculation on multi-cave combinations in the field data example may be available.

6. Conclusion

In this paper, we conducted the sparsity constrained inversion spectral decomposition (SCISD) for carbonate cave reservoirs to identify the multi-cave combinations. This approach is a time-frequency analysis method based on the convolution model, in which sparsity is imposed on the time-frequency spectrogram of the seismic signal by adding a regularization condition. Synthetic seismic trace examples illustrate that the energy spectrum of this method has high time and frequency resolutions and good time-frequency focusing property. The phase spectrum of the SCISD clearly shows the distribution of the seismic time-frequency energy and the phase variations of the seismic signal and also contains the information on the local phase.

In the physical model data examples, for the regular cave combinations, the combination types and the number of caves can be directly identified by the 20 Hz SCISD isofrequency section and the energy of the caves is more concentrated. Additionally, the apparent heights of the cave combinations are closer to their real heights, which indicates that SCISD is useful for the height estimation of a vuggy reservoir. The 20 Hz phase spectrum not only is able to distinctly show the energy distribution of the multi-cave combinations but also contains the corresponding phase information. Combined usage of the energy and phase spectra is beneficial for more accurate interpretation of the combination form and the number of the caves, and can provide a good basis for the quantitative estimation of the volume of the caves in the subsequent interpretation. For the columnar caves and layered vertical combinations, it is helpful for determining the combination forms as well as the boundary reflection of the caves by combining the energy spectrum with the phase spectrum; it can also provide an effective basis for subsequent accurate depiction of the cave reservoir distribution.

In the field data example, compared with STFT, the higher resolution and better noise immunity brought by SCISD highlights the locations and extent of the carbonate caves more clearly. The phase changes of the single cave and multi-cave combinations are clearly shown in the phase spectrum of SCISD. The combination of energy spectrum and phase spectrum can help us to deduce the real structure of the cave combinations.

Acknowledgments

We are very grateful to Dr Han Li for offering help to us. In addition, we are grateful for the support from Petro China Tarim Oilfield Company Institute for the experiment of the physical model. We also thank the Northwest Branch of Petro China Exploration and Development Research Institute for the data processing. The work was supported by the Key State Science and Technology Project of China (2011ZX05007-006). Special thanks go to the editors and reviewers for their constructive suggestions, which significantly improved the paper.

References

Beck
A
Teboulle
M
,
2009
,
A fast iterative shrinkage-thresholding algorithm for linear inverse problems
Siam J. Imaging Sci.
, vol.
2
(pg.
183
-
202
)
Bonar
D C
Sacchi
M D
,
2010
,
Complex spectral decomposition via inversion strategies
SEG Technical Program Expanded Abstracts
, vol.
vol 29
(pg.
1408
-
1412
)
Bonar
D C
et al.
,
2010
,
Time-frequency analysis via deconvolution with sparsity constraints
Geo-Canada Expanded Abstracts
(pg.
1
-
5
)
Castagna
J
Sun
S J
Siegfried
R W
,
2003
,
Instantaneous spectral analysis: detection of low-frequency shadows associated with hydrocarbons
Leading Edge
, vol.
22
(pg.
120
-
127
)
Chopra
S
Marfurt
K J
,
2016
,
Spectral decomposition and spectral balancing of seismic data
Leading Edge
, vol.
35
(pg.
176
-
179
)
Elboth
T
Presterud
I V
Hermansen
D
,
2010
,
Time-frequency seismic data denoising.
Geophys. Prospect.
, vol.
58
(pg.
441
-
453
)
Gardner
T J
Magnasco
M O
,
2006
,
Sparse time-frequency representations
Proc. Natl Acad. Sci. USA
, vol.
103
(pg.
6094
-
6099
)
Han
L
Sacchi
M D
Han
L G
,
2013
,
Spectral decomposition and denoising via time-frequency and space-wavenumber reassignment
Geophys. Prospect.
, vol.
62
(pg.
244
-
257
)
Liu
Y
et al.
,
2015
,
Inversion spectral decomposition with channel complex identification
85th Annual Int. Meeting, SEG, Expanded Abstracts
(pg.
1628
-
1632
)
Li
Y D
Zheng
X D
,
2008
,
Spectral decomposition using Wigner–Ville distribution with applications to carbonate reservoir characterization
Leading Edge
, vol.
27
(pg.
1050
-
1057
)
Li
Y D
Zheng
X D
Zhang
Y
,
2011
,
High-frequency anomalies in carbonate reservoir characterization using spectral decomposition
Geophysics
, vol.
76
(pg.
V47
-
V57
)
Lu
W K
Li
F Y
,
2013
,
Seismic spectral decomposition using deconvolutive short-time Fourier transform spectrogram
Geophysics
, vol.
78
(pg.
V43
-
V51
)
Marfurt
K J
Kirlin
R L
,
2001
,
Narrow-band spectral analysis and thin-bed tuning
Geophysics
, vol.
66
(pg.
1274
-
1283
)
Partyka
G
Gridley
J
Lopez
J
,
1999
,
Interpretational applications of spectral decomposition in reservoir characterization
Leading Edge
, vol.
18
(pg.
353
-
360
)
Portniaguine
O
Castagna
J
,
2004
,
Inverse spectral decomposition
SEG Technical Program Expanded Abstracts
, vol.
vol 23
(pg.
1786
-
1789
)
Puryear
C I
Portniaguine
O N
Cobos
C M
Castagna
J P
,
2012
,
Constrained least-squares spectral analysis: application to seismic data
Geophysics
, vol.
77
(pg.
V143
-
V167
)
Rodriguez
I V
Bonar
D
Sacchi
M
,
2012
,
Microseismic data denoising using a 3C group sparsity constrained time-frequency transform
Geophysics
, vol.
77
(pg.
V21
-
V29
)
Robinson
E A
Treitel
S
,
1980
Geophysical Signal Analysis
Englewood Cliffs, NJ
Prentice-Hall
Sejdic
E
Djurovic
I
Jiang
J
,
2009
,
Time-frequency feature representation using energy concentration: an overview of recent advance
Digit. Signal Process.
, vol.
19
(pg.
153
-
183
)
Shafi
I
Ahmad
J
Shah
S I
Kashif
F M
,
2009
,
Techniques to obtain good resolution and concentrated time-frequency distributions: a review
EURASIP J. Adv. Signal Process.
, vol.
2009
(pg.
1
-
43
)
Sinha
S
Routh
P S
Anno
P D
Castagna
J P
,
2005
,
Spectral decomposition of seismic data with continuous-wavelet transform
Geophysics
, vol.
70
(pg.
P19
-
P25
)
Wang
Y H
,
2007
,
Seismic time-frequency spectral decomposition by matching pursuit
Geophysics
, vol.
72
(pg.
V13
-
V20
)
Wang
H L
Siu
K
Ju
K
Chon
K H
,
2006
,
A high resolution approach to estimating time-frequency spectra and their amplitudes
Ann. Biomed. Eng.
, vol.
34
(pg.
326
-
338
)
Williams
G
Chadwick
A
,
2012
,
Quantitative seismic analysis of a thin layer of CO2 in the Sleipner injection plume
Geophysics
, vol.
77
(pg.
R245
-
R256
)
Xu
C
Di
B R
Wei
J X
,
2014
,
A seismic physical modeling study of cavernous carbonate reservoirs at seismic scale
SEG Technical Program Expanded Abstracts
(pg.
3579
-
3583
)
Yuan
S Y
Wang
S X
,
2013
,
Edge-preserving noise reduction based on Bayesian inversion with directional difference constraints
J. Geophys. Eng.
, vol.
10
025001
Yuan
S Y
Wang
S X
Luo
C M
He
Y X
,
2015
,
Simultaneous multitrace impedance inversion with transform-domain sparsity promotion
Geophysics
, vol.
80
(pg.
R71
-
R80
)