AVO Detuning Effect Analysis Based on Sparse Inversion

Abstract

The wave field characteristics of thin reservoirs are extremely complex due to the tuning and interference between the top and bottom interfaces of the reservoirs, which leads to large uncertainty in thin layer AVO (Amplitude Versus Offset) analysis. In order to reduce the uncertainty of thin layer AVO analysis, we study the uncertainty dominant factors of the effect of thin layer on the AVO response characteristics from the aspects of theoretical derivation and forward simulation. Based on the research results, we use the AVO fitting forward method with offset and tuning utility as the joint inversion operator to establish an AVO detuning effect method, based on the sparse fitting inversion strategy, and study the objective function of the fitting inversion method. We optimize the sparsity constraints and the sparsity method to reduce the non-independence of multiparameter variables and seismic data, and the noise of inversion. Through the verification analysis of the model using actual data, the AVO detuning effect method studied in this paper has a correct and reasonable technical theory and obvious application effect.

1. Introduction

With the gradual deepening of exploration and development, the chance of encountering thin or thin interlayer reservoirs is increasing. The tuning effect of thin interlayers affects the AVO analysis, leading to the difficulty of AVO application in thin interlayer reservoirs. Many scholars have carried out a lot of studies on the seismic response characteristics, frequency characteristics and thickness relationship of thin reciprocal reservoirs. Widess [1] studied the relationship between thin layer thickness and tuning amplitude using the zero-phase wavelets. When the thickness of thin layer is one-quarter wavelength, the reflected waveform is the first-order derivative of the incident waveform, which breaks the limit of finding the thickness of thin layer, meanwhile defining strata with thickness less than one-quarter wavelength as thin layer. Since then, the relationship between the thickness of the thin interlayer and the seismic response has revolved around the fact that when thin interlayer exists, the tuning effect, multiple waves between layers, and attenuation can all lead to significant differences between reflection characteristics and those of single interface, especially for the AVO response [2,3]. Simmons and Backus [4] pointed out that the conversion wave has some influence on the seismic response of the thin layer and the real AVO response of the thin layer cannot be accurately obtained with the Zoeppritz equation alone, and also analyzed the influence of conversion waves on the AVO response of thin layers. Liu and Schmitt [5] derived an exact solution for the AVO response of thin reservoir amplitude with offset, analyzed the influence of reservoir thickness and Poisson’s ratio on the AVO response. Pan et al. [6] pointed out that Zoeppritz equation is not suitable for analyzing the reflection and transmission response of thin layers, and analyzed the AVO characteristics of thin layers in elastic media using the Brekhovski equation, while there are differences in the response of thin layer in acoustic and elastic wave media, and that the AVO curves of strata with severe attenuation are close to the results of the Zoeppritz equation. Wang [7] and Lu [8] simplified the Brekhovski equation and used the simplified equation to give the characteristics of AVO response in different thin layer models, and proved through simulated annealing inversion that the elastic parameters obtained from the simplified equation are more reasonable than those obtained from the Zoeppritz equation, and they analyzed the AVO gradient changes of seismic reflection amplitudes of the target layer in different frequency bandwidths, which solved the AVO identification problem of thin interbedded reservoirs. Guo [9] studied the relationship between the physical parameters of thin layers and their seismic reflection characteristics using the Gassmann equation and the AVO technique based on thin layer reflection dynamics and viscoelastic anisotropy theory. Liu et al. [10] analyzed the variation of seismic reflection characteristics with the cumulative thickness of thin interbedded sandstone, the thickness of sandstone and mudstone interlayers by orthorectified simulation of a wedge model through the seismic convolution method. However, there are fewer studies involving AVO tuning correction.

We study the characteristics of AVO and its correction method from the aspects of theoretical derivation and forward simulation. The expression of thin layer coupled AVO reflection is derived and combined with the influence of thin layer tuning and interference. Conceptual models for different thicknesses of thin layers are constructed, and the wavefield orthorectification is carried out using the reflectivity method to investigate the effects of thin layer structures on AVO response characteristics. Based on the simplified model, the AVO detuning effect method based on the sparse fitting inversion strategy is investigated. The AVO fitting orthorectification method with offset and tuning effect as joint orthorectification operators is discussed, and sparse constraints methods are introduced to reduce the multi-solution of the inversion in order to address the non-independence of multi-parameter variables and the noise problem of seismic data in the inversion problem. The validity of the research method is demonstrated by the validation analysis of the model and the actual data.

2. Geological Setting

The Qiongdongnan Basin (QDNB) is a Cenozoic basin located on the continental shelf in the northern part of the South China Sea (Figure 1). The basin is separated from the Yinggehai Basin by the Deep Great Rift in the west, from the Zhujiangkou Basin by the Shenhu Uplift in the east, the Hainan Uplift in the north and the Xisha Uplift in the south. The QDNB comprises five secondary tectonic units: the Hainan Uplift, the Northern Depression, the Central Uplift, the Central Depression and the Southern Uplift. The QDNB has undergone two phases of tectonic evolution, the Paleoproterozoic Faulting Stage and the Neoproterozoic Depressional Stage, forming a typical ‘double-layered structure’ of the faulted basin. The basin is filled from bottom to top by the Eocene, YC and Lingshui formations in the Paleocene rifting stage, and the Sanya, Meishan, Huangliu and Yinggehai formations in the Neoproterozoic depression stage. The QDNB is an important target for deep-water exploration in the South China Sea, and AVO technology has been used to achieve good results in this area [11,12]. Fluid-sensitive information extracted from pre-stack information can effectively reduce the multi-solution of strong amplitude highlights, and several medium and large oil and gas fields in LS17-2, LS18-1 and LS25-1 have been discovered, demonstrating good exploration prospects [13,14]. However, with the progressive exploration, the large reservoir targets were exhausted and then shifted to the turbidite waterway sandstone with more complex stratigraphic assemblage, which is limited by the high cost of deep-water exploration, few drilling wells and the uncertainty of AVO in the thin layer, putting forward higher requirements on AVO analysis for deep-water exploration. In this study, the LS31 channel in the QDNB in the western part of the South China Sea is taken as the target area.

3. Methods

The theory of thin layer reflection coefficient, which discusses the reflection and transmission propagation of a plane wave in the layered medium from the perspective of elastic wave dynamics, is based on the Brekhovski equation. Suppose two semi-infinite elastic media are sandwiched by the horizontal media (Figure 2). If the thickness of the stratum is less than one wavelength, the fluctuations at the interface between the top and bottom of this stratum will interfere and superimpose significantly, so that the fluctuations formed at the top surface of the thin layer will have different characteristics from those of the thick layer in all aspects of intensity variation and waveform characteristics. The result of the superimposed reflection depends mainly on the thickness of the thin layer, the propagation velocity of the thin layer and the angle of incidence of the plane wave.

3.1. Theory of Thin Layer Reflectance Coefficient

Let the plane wave be incident from the top layer into the below layer, and the angle of incidence of the P-wave into each layer is $\theta$, then the displacement function expressions for each wave in the top layer [15] can be expressed as

$${\varphi }_e=A_D^2\exp \left[i\left(\omega t-k_{x}x-k_{dz}z\right)\right]$$

(1)

$${\phi }_e=B_D^2\exp \left[i\left(\omega t-k_{x}x-k_{sz}z\right)\right]$$

(2)

$${\varphi }_r=A_U^2\exp \left[i\left(\omega t-k_{x}x+k_{dz}z\right)\right]$$

(3)

$${\phi }_r=B_U^2\exp \left[i\left(\omega t-k_{x}x+k_{sz}z\right)\right]$$

(4)

where ${\varphi }_e$, ${\varphi }_r$ is respectively the displacement functions of the incident and the reflected P-wave; ${\phi }_e$, ${\phi }_r$ is respectively the displacement functions of the incident and reflected S-wave; $A$, $B$ is respectively the amplitudes of P-wave and S-wave; $D$, $U$ is respectively upward and downward wave; $\omega$ is angular frequency; $t$ is time; $k_x$, $k_{dz}$, $k_{sz}$ is respectively the number of horizontal waves in the longitudinal direction, the number of vertical waves in the longitudinal direction and the number of vertical waves in the transverse direction.

The displacement-stress relationship can be expressed as

$$\left[\begin{array}{c}{u}^{\left(1\right)}\\ w^{\left(1\right)}\\ {\sigma }_{zz}{}^{\left(1\right)}\\ {\sigma }_{zx}{}^{\left(1\right)}\end{array}\right]=A_{ij}^{\left(2\right)}\left[\begin{array}{c}{u}^{\left(3\right)}\\ w^{\left(3\right)}\\ {\sigma }_{zz}{}^{\left(3\right)}\\ {\sigma }_{zx}{}^{\left(3\right)}\end{array}\right]$$

(5)

where superscripts 1, 2 and 3 is respectively the first, second and third layer of the medium; $u$, $w$ is respectively the displacement function;${\sigma }_{zz}$ is the vertical stress; ${\sigma }_{zx}$ is the horizontal stress; $A_{ij}$ is the propagation matrix; superscript 2 is the propagation matrix $A_{ij}$ related to medium 2.

The formulae for calculating the reflection and transmission coefficients [16] can be expressed as

$$R_{pp}=\frac{A_2^1}{A_1^1}, R_{ps}=\frac{B_2^1}{A_1^1}\frac{{\alpha }_1}{{\beta }_1}, T_{pp}=\frac{A_1^3}{A_1^1}\frac{{\alpha }_1}{{\alpha }_3}, T_{ps}=\frac{B_1^3}{A_1^1}\frac{{\alpha }_1}{{\beta }_3}$$

(6)

where $R_{pp}$, $R_{ps}$, $T_{pp}$, $T_{ps}$ is respectively longitudinal reflection coefficient, transverse reflection coefficient, longitudinal transmission coefficient and transverse transmission coefficient; $\alpha$, $\beta$ is P-wave and S-wave velocity.

Reflected waves from the thin interlayer in multilayered media are composite waves formed by the superposition of a single wave in the thin interlayer after different numbers of reflections, so recursive algorithm can be used to calculate the displacement and stress relationship between the top and bottom of the thin interlayer, so that the reflection transmission coefficient of each layer can be found.

$$\left[\begin{array}{c}{u}^{\left(1\right)}\\ w^{\left(1\right)}\\ {\sigma }_{zz}{}^{\left(1\right)}\\ {\sigma }_{zx}{}^{\left(1\right)}\end{array}\right]=A_{ij}^{\left(2\right)}\times A_{ij}^{\left(3\right)}\times \cdots \times A_{ij}^{\left(n\right)}\left[\begin{array}{c}{u}^{\left(n+1\right)}\\ w^{\left(n+1\right)}\\ {\sigma }_{zz}{}^{\left(n+1\right)}\\ {\sigma }_{zx}{}^{\left(n+1\right)}\end{array}\right]$$

(7)

where $A_{ij}^{\left(2\right)}=B_{ij}^{\left(2\right)}{C}_{ij}^{\left(2\right)}$,$B_{ij}^{\left(2\right)}$, $C_{ij}^{\left(2\right)}$ can be expressed as

$$C_{ij}^{\left(2\right)}=\left[\begin{array}{cccc}\frac{i2k_x}{k_s^2}& 0& -\frac{1}{{\mu }_2{k}_s^2}& 0\\ 0& -\frac{i\left[k_{sz}^2-k_x^2\right]}{k_s^2{k}_{dz}}& 0& \frac{2k_x}{k_s^2{k}_{dz}}\\ \frac{i\left[\frac{\lambda }{{\mu }_2}{k}_d^2+2k_{dz}^2\right]}{k_s^2{k}_{sz}}& 0& \frac{k_x}{{\mu }_2{k}_s^2{k}_{sz}}& 0\\ 0& \frac{i2k_x}{k_s^2}& 0& \frac{2}{k_s^2}\end{array}\right]$$

(8)

$$B_{ij}^{\left(2\right)}=\left[\begin{array}{cccc}i{k}_x\cos p& -k_x\sin p& -i{k}_{sz}\cos q& -k_{sz}\sin q\\ k_{dz}\sin p& i{k}_{dz}\cos p& -k_x\sin q& -i{k}_x\cos q\\ -M\cos p& M\sin p& 2{\mu }_2{k}_{xz}\cos q& -i2{\mu }_2{k}_x{k}_{xz}\sin q\\ -i{k}_x{k}_{dz}\sin p& k_x{k}_{dz}\cos p& -iN\sin q& N\cos q\end{array}\right]$$

(9)

where $p=k_{dz}H$, $q=k_{sz}H$, $M={\lambda }_2{k}_x^2+\left({\lambda }_2+{\mu }_2\right)k_{dz}^2$, $N=1/2\left(k_{sz}^2-k_x^2\right)$.$k_x$ is slowness in the horizontal direction; $k_s$, $k_d$ is respectively P-wave and SV-wave slowness; $k_{sz}$, $k_{dz}$ is respectively P-wave and SV-wave vertical slowness; H is the thickness of the stratum.

3.2. Thin Layer AVO Curves Versus Stratigraphic Thickness

When the thickness of the thin layer is equal to one quarter of the seismic wavelength, the seismic wave field appears to interfere in the same phase, making the amplitude appear extremely large, and this phenomenon is called the tuning effect of the thin layer. In order to analyze the relationship between the thickness of the thin layer, the main seismic frequency and its pre-stack AVO characteristics, a single sand wedge model is designed (Figure 3), and the elastic parameters of the wedge model are shown in

Table 1. The elastic parameters of the wedge model.

AVO	Mudstone Enclosures			Sandstone			Zoeppritz Calculation
	P-Wave (m·s−1)	S-Wave (m·s−1)	Density (g cm−3)	P-Wave (m·s−1)	S-Wave (m·s−1)	Density (g cm−3)	Intercept	Gradient
III	2810	1480	2.53	2750	1750	2.12	−0.023	−0.011

. The results of the pre-stack AVO simulation are shown in Figure 4. When the stratigraphic thickness is greater than the tuning thickness, the simulated results show that the top surface of the sandstone exhibits III AVO characteristics, which is basically the same as the actual AVO characteristics of the stratigraphic layer, and the variation of amplitude with the angle is more-or-less stable without much different, indicating that the simulation results of the thin layer reflection coefficient formula and the Zoeppritz equation are the same when the stratigraphic thickness is thicker than the tuning thickness.

We performed rendezvous analysis of the intercept and gradient of the simulation results for different dominant frequencies (Figure 5). As a whole, there is a spiral distribution in the intercept and gradient rendezvous, with different dominant frequencies having different abilities to identify the thin layer. For the dominant frequency of 10 Hz (where half-wavelength = 137.5 m), when the thickness of the layer is less than a half-wavelength, the intercept and gradient of the top surface of sandstone is less than a half-wavelength deviation from their true values (red dots are the values of the Zoeppritz equation calculation when the thickness of the layer is half-wavelength) due to the interference between the top and bottom reflection interfaces. When the dominant frequency is 20 Hz (half-wavelength = 68.75 m), the intercept gradient converges to the same position for sandstone thicknesses greater than a half-wavelength, and basically does not change with the thickness of the formation. However, when the sandstone thickness is less than a half-wavelength, the intercept and gradient are principally influenced by the thickness of the formation: when the sandstone thickness is less than the tuning thickness, the intercept increases with the sandstone thickness; when the sandstone thickness is greater than the tuning thickness, the intercept decreases with the sandstone thickness and the gradient follows a similar pattern of change characteristics, but the position of the inflection point is delayed. As the dominant frequency gradually increases, the wavelength progressively decreases (45.8 m at 30 Hz half-wavelength, 34.4 m at 40 Hz half-wavelength, and 27.5 m at 50 Hz half-wavelength, using the layer P-wave velocity of 2750 m/s), and the intercept and gradient converge to the true value of the layer when the layer thickness exceeds a half-wavelength, and the higher the main frequency the better the recognition of the thin layer.

3.3. AVO Detuning Effect Analysis Based on Sparse Inversion

In the tuning thickness range, how to correct the amplitude and phase data, remove the tuning effect and the stretching effect, and restore the influence of the offset factors on the waveform, is the key to the successful application of AVO technology [17,18,19]. The previous studies show that the top-bottom interference in the tuning thickness range has a large impact on the thin layer’s true AVO relationship, and here a sparse inversion strategy is introduced to the detune effect method.

Theoretically, zero-offset seismic traces have only wavelet and stratigraphic tuning effects, while non-zero offset traces have offset tuning effects in addition to wavelet and stratigraphic tuning effects. The fitted AVO detuning inversion uses the AVO attribute as the model’s spatial parameter, and uses the Shuey approximation formula to obtain the attribute parameter model without the tuning effect by fitting the approximation to the actual pre-stack seismic trace set data with the tuning effect.

The binomial approximation of AVO reflection coefficient [20] can be expressed as

$$R\left(\theta \right)=P+G{\sin }^2\left(\theta \right)$$

(10)

where $R$ is the reflection coefficient, $P$ is the intercept, $G$ is the gradient, $\theta$ is the angle. The above equation can also be expressed in matrix form

$$R\left(\theta \right)=CM$$

(11)

where $M={\left[P,G\right]}^T$, $C=\left[1,{\sin }^2(\theta )\right]$.

The detuning effect inversion problem can be expressed [21] as

$$I_{min}=\arg u\left[\min \left(C{M}^{\prime }-d\right)\right]$$

where $d=CM$ is actual pre-stack seismic data. $M^{\prime }={\left[P^{\prime },G^{\prime }\right]}^T$ is the intercept and gradient matrix after detuning effects.

Whether it is the acoustic impedance tuning effect or the elastic tuning effect, the parameters are not independent, and the seismic data itself is noisy (especially the far offset seismic data, which is noisier) which further increases the degree of ambiguity of the parameter variables. The inversion problem is severely unsuitable due to multiple parameters, non-independence between parameters, and noise. In order to adapt the inversion problem to reduce the multi-solution nature of the inversion, constrained inversion is generally performed [22,23]. In this paper, we propose the use of the constraint sparse inversion method of model parameters and seismic data. The inverse problem can be expressed as

$$I_{min}=\arg u{\left[\min \left(C{M}^{\prime }-d\right)\right]}^{,s.t.}{\Vert d\Vert }_0$$

(12)

where ${\Vert d\Vert }_0$ is sparse optimization results for reflection coefficients or reflection amplitudes.

The seismic pre-stack gather is the convolution of the reflection coefficient and wavelet, and the seismic record is projected onto the wavelet to achieve sparse transformation. The degree of sparsification is related to the band width of the wavelet: the wider the band width of the wavelet, the narrower the wavelet pulses in the time domain, and the greater the degree of sparsification of the seismic record after the sparse transformation by wavelet projection [24,25]. The main parts of this sparsification process include the creation of a wavelet dictionary library and hierarchical sparsification.

Firstly, for wavelet dictionary library building, wavelet vector can be expressed as

$$\mathbf{W}={\left[w_0,w_1,w_2,\dots ,w_{N-1}\right]}^T$$

(13)

Secondly, hierarchical sparsification. A basic idea of sparse reconstruction is to arrange and intercept data from the largest to the smallest, and the sparsity is the number of selected data from the largest to the smallest. If a signal is continuous and has a certain distribution, the complete signal can be reconstructed with fewer samples. The reflection coefficients are not continuously predictable in the longitudinal direction, and the corresponding seismic records are also not continuously predictable in the longitudinal direction [26,27,28]. To address these problems, we propose a hierarchical sparsification reconstruction method. The method is based on the following idea: for a section of a seismic record, the maximum and minimum values (absolute values) are obtained, and the data are divided into multiple sequences of different levels from large to small at a certain segmentation spacing, and then each sequence is reconstructed by sparse reconstruction, and finally the reconstruction results of each level are then weighted and summed.

In order to verify the reliability of this technical method, this paper randomly designed a geological model with different thin layer combinations for pre-stack forward simulation (Figure 6), where the stratigraphic real shows III AVO response characteristics, and only the stratigraphic combination of the thin layer is changed. The analysis of the forward simulation shows that the actual stratigraphy has both III and IV AVO characteristics due to the tuned influence of the thin stratigraphic assemblage, indicating that the tuning effect has changed the real AVO of the actual stratigraphy, resulting in some uncertainty in the AVO analysis. The sparse inversion strategy proposed in this paper is used to detune this effect, and the correlation analysis with the theoretical orthodromic shows that (i) the de-coupling effect and the theoretical orthodromic AVO types basically match each other, showing III AVO response characteristics, and (ii) the correlation coefficients with the theoretical orthodromic AVO can reach above 90%, which can better reflect the real AVO response characteristic of thin gas-bearing sandstones and effectively reduce the multi-solvability of AVO analysis.

4. Results

LS31 channel, which is located in the QDNB of the South China Sea, is geological conditioned for developing large and medium-sized oil and gas fields [29,30,31]. A statistical analysis of the AVO types of sandstone reservoirs in this area at different depths shows that the reservoir in the central canyon channel in the deep-water area has good physical properties: the gas-bearing sandstones are mostly IV AVO, the gas-bearing sandstones away from the canyon channel are mostly III AVO, most of the water-bearing sandstones show IV AVO characteristics [32,33], and the gas-water difference is relatively obvious. Consider the post-stack seismic profile and the pre-stack gather across the LS31 channel as an example (Figure 7). In the post-stack seismic profile, the channel features are obvious, with the top surface of the channel showing negative-phase trough reflections and the bottom surface showing positive-phase reflections and obvious undercutting features, with the middle thickening and the sides gradually thinning. However, in the pre-stack gather, the AVO characteristics of the top surface of the LS31 channel vary greatly due to the influence of the change in the thickness from the edge to the main body of the channel, and III AVO and IV AVO existing between each other, especially at the position of the tuned thickness; therefore, the AVO curve shows an obvious V-shape variation, with III AVO characteristics at small angles in the near gathers and IV AVO at large angles in the far gathers. This makes the analysis of AVO in this area very difficult.

We use the detuning method proposed in this paper to compare and analyze the PG attribute profiles of the high and low parts of the LS31 channel, and the AVO change characteristics before and after the detuning (Figure 8). Before the correction of the PG attribute (Figure 8a), due to the interference of sandstone thickness, the high and low parts of the target area are mostly dominated by III AVO and IV AVO coexistence, and the AVO anomaly of the sandstone in the main body of the LS31 channel is not clear. We also show the detuned profile using the conventional method (waveform correction) due to the influence of stretching factors, and also the near and far channel frequency change due to the gun inspection distance, and thus the tuning effect of AVO is not effectively improved (Figure 8b). After the correction of the detuning effect via sparse inversion (Figure 8c), the AVO anomalies in the high and low parts of the target area are relatively consistent, all showing mainly III AVO, and the AVO variation characteristics before and after detuning are basically consistent with the actual stratigraphic AVO, and the V-shape feature on the single point AVO gather is practically eliminated (Figure 9).

From the AVO attribute extracted before and after detuning (Figure 10), it can be seen that the pre-detuning AVO shows the coexistence of III and IV AVO, and after correction by sparse inversion detuning, the influence of the detuned V-shape on AVO is basically eliminated, which is basically consistent with the actual drilled wells. The L-1 well encountered 18 m of water-bearing sandstone in the Yinggehai Formation, which exhibited IV AVO characteristics, which correlates with what was observed after detuning; however, before detuning the AVO characteristics were not obvious, and both III and IV AVO characteristics were observed. Combining the detuned PG attributes, drilling well L-2 was deployed in the main part of the Yinggehai Formation, producing a high-quality gas reservoir of 41.8 m with a porosity of 20.8%, adding approximately 9.6 billion cubic meters of proven reserves. Subsequently well L-3 is to be drilled to progress and expand the scale of reserves. This technology has opened up a new dimension of exploration in the shallow Yinggehai Formation in the Lingshui Depression, and has effectively promoted the exploration and development process of the deep-water gas fields in the western part of the South China Sea.

5. Conclusions

The forward simulation analysis of tuned AVO in thin reservoirs reveals that the difference in the tuning effect of near gathers and far gathers in the tuning thickness range is a key factor leading to the multiresolution of thin AVO. In order to solve this problem, the AVO detuning effect method based on the sparse fitting inversion strategy is studied, which can better solve the problem of multi-solutions of multivariate coupled inversion. At the same time, the fitted AVO inversion method can effectively remove the AVO tuning effect and reduce the uncertainty of AVO analysis in thin reservoirs, which is more beneficial to the AVO inscription analysis of thin reservoirs. The DF gas field and LD gas field in Yinggehai Basin are also faced with the uncertainty of AVO analysis of thin reservoirs. The technical method proposed in this paper has been effectively verified and applied in the small well area of Yinggehai Basin, and could be further employed to solve the uncertainties experienced at other wells, both in this Basin and beyond.