Abstract

Wave-induced fluid flow is the main cause of seismic attenuation and dispersion. So the estimated velocity dispersion information can be used to identify reservoir fluid and effectively reduce the risk of reservoir drilling. Using equivalence of dispersion and attenuation between poroelastic and viscoelastic media, we developed the method of FAVO (frequency-dependent amplitude variation with offset) dispersion quantitative estimation based on the analytical solution of 1D viscoelastic wave equation. Compared with the current single-interface velocity dispersion estimation method, the new nonlinear approach uses the analytical solution of 1D viscoelastic wave equation as the forward modeling engine. This method can conveniently handle the attenuation and generate the full-wave field response of a layered medium. First, the compound matrix method (CMM) was applied to rapidly obtain the analytical solution by vectorization. Further, we analyzed the seismic response characteristics through the model data to clarify the effectiveness of the forward modeling method. Then, the more reliable P-wave velocity, S-wave velocity, and density were recovered based on prestack viscoelastic waveform inversion (PVWI). Combining with the inversion results, the derivative matrix was calculated to perform nonlinear velocity dispersion estimation. Finally, the new estimation method was tested with the model and actual data. The experiments show that the developed method is clearly superior to the single-interface dispersion estimation method in accuracy and resolution. This approach can be used as a new choice reservoir fluid identification.

1. Introduction

AVO technology is a common method for fluid identification and is widely used in hydrocarbon exploration industries. This approach uses the information of amplitude variation with offset/angle to quantitatively obtain the subsurface elastic parameters including P-wave velocity, S-wave velocity, and density [13]. Then, the inversion results and their derived parameters (P-wave impedance, VP/VS ratio, and Poisson’s ratio) can be used for fluid identification [4]. With the deepening of seismic exploration, the conventional offset dependent only technology can no longer meet the demand of seismic fluid identification and reservoir prediction [5]. And using frequency-dependent seismic information for fluid identification has become a research hot topic [6, 7]. As an organic combination of AVO inversion (offset dependent) and frequency-dependent fluid identification technology, FAVO can make full use of the information of reflectivity variation with offset and frequency to perform more effective fluid identification [8, 9].

In the past few decades, the theory of FAVO has been gradually improved. On the one hand, the study on rock physics mechanisms of the attenuation effect laid a theoretical foundation for identifying fluid by using velocity dispersion characteristics [1012]. White [12] proposed the theory of patchy-saturated equivalent medium at the mesoscopic scale to explain the attenuation and dispersion in the seismic exploration frequency band (1–100 Hz). After theoretical analysis [13, 14] and laboratory experiments [15, 16], it is widely accepted that the wave-induced fluid flow at mesoscopic scale is the main cause of seismic wave intrinsic attenuation in porous media. On the other hand, the FAVO forward method transits from the early single-interface assumption [17, 18] to the layered medium assumption [1921]. This makes the seismic reflection response more reasonably related to the attenuation and dispersion [22]. The velocity dispersion can be quantitatively estimated from the information of seismic reflection [8], and velocity dispersion is related to fluid flow. So it can be used as a fluid indicator.

Based on the Smith-Gidlow approximation, Wilson et al. [23] firstly proposed the FAVO inversion for quantitative estimation dispersion properties from seismic reflection data and used model data to verify the effectiveness of the method. The implementation process of the FAVO inversion method can be briefly summarized as the following key steps: ① calculating the time-frequency spectrum of seismic traces using spectral decomposition techniques such as short-time Fourier transform [24], wavelet transform [25], S transform [26]; ② performing the spectral balancing technique to eliminate the overprint effect of wavelets; thus, we can obtain the time-frequency spectrum of the reflection coefficient [27]; ③ by selecting the key frequency points for calculation, the velocity dispersion properties can be estimated by the Smith-Gidlow approximation. Since then, based on different approximations and spectral decomposition methods [5, 9, 28], the ability of FAVO inversion to identify fluids has been verified on actual data.

From a theoretical point of view, Wilson-based dispersion estimation methods are not satisfactory and remain controversial in many important aspects. First, these methods are based on the single-interface assumption. Seismic data consist of primary reflections only and is not “contaminated” by transmission loss effect, intermultiples, and interconversion waves. This will lead to incorrect fluid identification. Second, based on small-angle approximation [2931], the calculated reflection coefficient is unreliable at large angles (>30°) [32]. Unfortunately, the large-angle seismic data is even more critical in fluid identification. Third, these methods only consider the attenuation and dispersion at the interface and ignore the attenuation and dispersion of propagation in the medium. In fact, the propagation attenuation effect is far greater than the interface attenuation effect. Therefore, it is necessary to completely eliminate the propagation attenuation effect before the dispersion estimation. Obviously, this is difficult to achieve. Based on the equivalence of dispersion and attenuation between poroelastic and viscoelastic media [3336], we can use the analytical solution of the 1D viscoelastic wave equation to address above issues.

In this study, we developed the nonlinear FAVO dispersion quantitative estimation method. Based on the Kolsky–Futterman (KF) attenuation model [37, 38], the compound matrix method (CMM) was used to rapidly obtain the analytical solution of the viscoelastic wave equation. After prestack viscoelastic waveform inversion (PVWI), we obtained the more reliable P-wave velocity, S-wave velocity, and density. Combining with the inversion results, the derivatives matrix was calculated to perform the nonlinear velocity dispersion estimation. Finally, the validity and practicability of the developed method was verified by the model and actual data.

2. Theory

2.1. Forward Modeling

In general, we can use the propagator matrix algorithm [39, 40] to obtain the analytical solution. However, because of numerical difficulties [41], the algorithm is always unstable. In order to rapidly generate the synthetic seismograms without the numerical difficulties, CMM (Phinney, et al., 2010) was used to solve the 1D wave equation. The detailed description of the vectorization method is as follows.

For a stack of horizontal layers, the total -reflectivity can be derived by the transfer function where is the determinant associated with the system matrix, the value of which does not affect the calculation result. is PP-reflectivity, is PS-reflectivity, and others are similar. does not have physical meaning.

To compute the total transfer function of the horizontal attenuated layers, we redefine the wave propagator matrix for more convenient calculations: where ; matrix is the phase shift of reflection:

And , represents downward and upward interface energy distribution matrices, respectively. The matrix contains 16 independent elements:

The matrix can be expressed by elements of matrix : where is the thickness, and are the vertical slowness of the P-wave and S-wave, respectively. , . The horizontal slowness is given as . It should be noted that and , respectively, represent the complex P- and S-wave velocity; the real part of which represents the corresponding phase velocity.

Based on the KF model, the P- wave and S-wave velocity complex expression of frequency [29] is

and are reference velocities. and are quality factors. The parameter is reference frequency. Same as inverse Q filtering [42] and for calculation convenience, we choose the main frequency of the wavelet as the reference frequency. In equation (6), is the dispersion term, which controls the phase of the wavelet, and is the absorption term only decreasing seismic amplitudes. It is noteworthy that S-wave attenuation has less impact on PP-seismic records [12]. For the convenience of calculation, we can let .

There is an elastic half-space below the base of the medium. can be expressed as

Snell’s law is also applicable in viscoelastic media [43]. Substituting equation (6) into recursive equation (2), the total PP-reflectivity with attenuation in domain can be given by :

Oliveira et al. [44] used a simple mathematical transformation to obtain the angle gather seismogram, but the stretch occurs at large angles. To avoid the stretching, we rewrite the phase shift matrix .

We then obtain the angle gather seismogram by Inverse Fourier transform: where is the source wavelet in the frequency domain.

2.2. Nonlinear Inversion Method

According to Wilson FAVO inversion theory, the information of reflectivity variation with frequency is directly related to the dispersion of velocity. The reflectivity function needs to be first-order approximated around a frequency point by its expansion in a Taylor series. In this study, the nonlinear reflectivity function can be expressed as: where is the nonlinear operator; represents model parameters at different frequencies. The first-order approximation of the function in is: where . The derivative of the reflection coefficient with respect to the model is written as: . Dispersion property, the derivative of the model parameter with respect to the frequency, can be written as: . Usually, we can use the P-wave dispersion property as the fluid indicator. For prestack seismic data, equation (12) can be rewritten as:

The left term of the equation is related to the time-spectrum of the reflection coefficient. This means we need to obtain the time-frequency spectrum of the reflection coefficient from seismic traces. First, the time-frequency spectrum of seismic traces can be calculated using spectral decomposition techniques. Then using the frequency domain wavelet , we get

According to the equations (13) and (14), the equation (15) is used to eliminate the overprint effect of wavelet.

This step is called the spectral balance. When the wavelet is inaccurate or unknown, Wu et al. [9] removed the effect of the source wavelet by designing a suitable weight function, where ; represents the time window.

When the operator is constructed by Smith-Gidlow approximation, the operator is linear, and the equation (13) is the famous Wilson FAVO equation. For nonlinear dispersion estimation method based on the analytical solution of the wave equation, obtaining the derivatives is an essential step. Under the viscoelastic medium, applying the chain rule, we can obtain: where denotes complex P-wave velocity, complex S-wave velocity, and density of layer , respectively, denotes P-wave reference velocity, S-wave reference velocity, and density. The elements for the partial derivatives of reflectivity can be derived by equation (8): so the derivatives of the viscoelastic medium can be expressed as:

The model parameter can be more accurately estimated by PVWI. This method is the compromise solution between full-waveform inversion and AVO inversion, and the resolution can live up to expectations of reservoir prediction. The two key steps about PVWI: forward modeling method and Fréchet derivatives have been presented in the above theory. Like other waveform inversion methods [45], PVWI can be easily implemented based on these two steps.

Figure 1 shows the nonlinear FAVO inversion workflow. Different from FAVO inversion based on single interface, we need to perform PVWI inversion to obtain accurate and reliable elastic parameters. Then, the inversion result is used to calculate the derivatives matrix. Finally, selecting the key frequency points for calculation, the velocity dispersion properties can be estimated.

3. Synthetic Example

3.1. Model Analysis

An unequal thickness interlayer model was used to demonstrate the effectiveness of the analytical solution method (ASM) in forward modeling. As shown in Figure 2, thin interbeds of nearly 400 ms act as protruding interlayer reflections (interlayer multiples and interlayer converted waves), and thick layers for highlighting transmission loss and attenuation effects.

With a range of incidence angle is 0–45° and the angle interval is 5°, the synthetic data were generated by a 40-Hz Ricker wavelet. The quality factor Q is a constant of 60. Figure 3 shows the comparison of normalized angle gathers that are, respectively, generated by the single-interface viscoelastic AVA equation (VAE) [17] and the ASM. Part1 and Part2 show angle gathers by the VAE and ASM, respectively, and Part3 is the residual between VAE and ASM. Since only the attenuation of the interface is considered, angle gathers synthesized by VAE have no significant resolution change and waveform distortion. This is not consistent with the attenuation phenomenon. The blue frame in Figure 3 reveals the apparent interlayer reflection, which often acts as effective weak reflections. With the increase of incident angle, the reflection coefficient value becomes larger. As shown in the red dotted box, the transmission loss and effect of attenuation have an increasing influence on the primary reflections with depth. It is worth noting that those propagation effects are difficult to completely eliminate by the current geophysical processing approach.

Using spectral decomposition technique (short-time Fourier transform), the variation curves of seismic amplitude with angle and frequency were calculated. After eliminating the overprint effect of wavelet, the reflection coefficient variation with angle and frequency at the reflection (indicated by the red arrow in Figure 3) was shown in Figure 4. Figure 4(a) shows the VAE reflection coefficient variation curves with angle at different frequencies (10–70 Hz). Since the attenuation of the interface is not obvious, the reflection coefficient changes slightly with frequency. The ASM reflection coefficient variation curves are shown in Figure 4(b). Different from VAE, the attenuation effects of interface and propagation are all considered. With the increase of frequency, the reflection coefficient decreases, but the AVO trend appears consistent. In summary, the propagation attenuation effect is far greater than the interface attenuation effect. Therefore, using ASM instead of VAE as the forward engine is a more reasonable choice.

3.2. Model Testing of the Nonlinear Dispersion Estimation Method

The nonlinear dispersion estimation method is then tested by using a group of well data. The curves of , , , and Q are, respectively, shown in Figure 5. The , , and data were taken from field well logging after Backus-averaging, and the Q value was calculated by empirical formula [46, 47] (in this actual study area, and ). With the 40-Hz Ricker wavelet, the synthetic angle gathers (Figure 6) were generated by ASM (angles range of 0–50° and intervals of 5°).

Following the workflow shown in Figure 1, the angle gathers were taken as the input gathers, and we performed the PVWI. From the mathematical point of view, the attenuation is expressed as an integral effect: on seismic amplitude [38]. So we can use a low-frequency trend of the value instead of the real value curve to perform compensation or inversion [48]. Using the low-frequency trend of the value, the recovered elastic parameters are shown in Figure 7(a). We can found that the inversion results are quite consistent with the true elastic parameters in numerical accuracy. So it is effective to apply the low-frequency for inversion, which will facilitate actual data inversion. For comparison, the angle gathers were assumed to be the real seismic data in which the propagation effect is not removed, and then, we performed single-interface AVO inversion. Because of the propagation effects, especially attenuation, the inverted elastic parameters (Figure 7(b)) are significantly deficient in accuracy and resolution. As the depth increases, the inversion bias also increases.

After spectral decomposition and spectral balance, the P-wave dispersion was estimated by the nonlinear frame. As shown in Figure 8(a), there is a certain bias between the real dispersion obtained by the KF model and the estimated dispersion. After reviewing the workflow, we thought there were several reasons for the bias. First, according to Heisenberg’s uncertainty criterion, we cannot obtain spectral decomposition results with high time resolution and high-frequency resolution (point spectrum). Second, the wavelet overprinting cannot be completely eliminated (just like deconvolution). Third, the method fits the dispersion properties at the reference frequency through reflection information of other frequencies, so the estimation result is influenced by frequency selection. Please note that whatever the linear or nonlinear estimation methods all remain these limitations. Nevertheless, the estimated dispersion can reflect the mid-low frequency trend of the real dispersion. Compared with the current single-interface FAVO inversion (Wilson-based) results (Figure 8(b)), the nonlinear method is clearly superior in accuracy and resolution. In a word, even if these restrictions remain, this method can be used as a higher-precision reservoir fluid identification technology in real data.

4. Real Data Example

The practicability of the developed method was verified by the actual data. The seismic data was obtained from work area X in northern China. Figure 9 shows the interest seismic poststack section with highly developed sag (time window 2.3–2.8 s and interval 2 ms). The seismic data for the FAVO inversion consists of 200 angle gathers (angle range 3–43° and angle interval 8°).

First, we extract the wavelet (Figure 10(a)) based on a statistical method. For selecting the key frequency points for calculation, the wavelet spectrum was calculated as shown in Figure 10(b). From the wavelet spectrum, we can estimate the wavelet dominant frequency and use it as the reference frequency Hz. Then, we selected the stagnation and inflection points of the amplitude spectrum function around as the frequency calculation points. In theory, in addition to the reference frequency, only one calculation frequency is needed to estimate dispersion. In order to improve the accuracy and stability of the dispersion estimation method, we generally select 4-6 key frequency points. As shown in Figure 10(b), the four key frequencies are , respectively.

After horizon picking, the initial model of P- and S-wave velocities and density were generated by interpolation with four wells. The quality factor was determined by the empirical formula , where is from the initial model. Through the regional attenuation background survey, the coefficient was corrected using logging data (, ). As shown in Figure 11, Figures 11(a)11(d), respectively, represent the initial model of P- and S-wave velocities, density, and the low-frequency trend of the Q. Then, we performed the PVWI. The inverted parameters are shown in Figures 12(a)12(c). For a comparison, the parameters estimated by single-interface AVO inversion are shown in Figures 13(a)13(c). As a result of the combination of attenuation compensation with inversion, the inversion resolution of PVWI is significantly higher than AVO inversion. We can conclude that PVWI can provide more accurate and higher-resolution inversion results.

After spectral decomposition and spectral balancing, we estimated P-wave dispersion by the nonlinear FAVO (Figure 14(a)). The dispersion of Wilson-based FAVO inversion is shown in Figure 14(b). Due to the difference in the forward method, there is a large difference in the resolution of the two dispersion properties. The spectral decomposition frequency sections of the poststack data at 16 Hz, 22 Hz, 28 Hz, and 34 Hz are shown in Figures 15(a)15(d). The low shadow (indicated by the red ellipse in Figure 15) is consistent with the fluid factor indicated by the high value of the dispersion property (the red ellipse in Figure 14(a)). This again verified the validity of the nonlinear dispersion estimation. Unlike the qualitative low-frequency shading technique, FAVO inversion can quantitatively indicate the fluid. Later drilling verification confirms that the red ellipse area is a favorable hydrocarbon reservoir.

5. Conclusions and Discussion

We developed the method of FAVO dispersion quantitative estimation based on the analytical solution of the 1D viscoelastic wave equation. Compared with the current single-interface forward method, this method can conveniently handle the attenuation and generate the full-wave field response of a layered medium. The model and actual data experiments show that the developed method is clearly superior to the single-interface dispersion estimation method in accuracy and resolution. Unlike the qualitative low-frequency shading technique, FAVO inversion can quantitatively indicate fluid. This approach can be used as a new choice reservoir for fluid identification.

As described in the model inversion section, FAVO inversion remains some limitations. And current geophysical technologies cannot handle these limitations. In addition, the several issues remain to be discussed. First, the attenuation effect is the main cause of the frequency-dependent phenomenon, but it is not unique. The frequency-dependent phenomenon also occurs due to the stretching of normal moveout, partial stacking, and multiple-wave interference [7]. On the one hand, we should be careful to avoid the frequency-dependent caused by numerical reasons during the seismic processing. On the other hand, the estimated dispersion results can be combined with conventional fluid factors to reduce the risk of false fluid indications.

Second, the physical meaning of the dispersion attribute is the derivative of the velocity dispersion function at a reference frequency. In theory, choosing different reference frequencies can obtain the derivatives at different frequencies. Using the Petrophysical model [10], this frequency-dependent dispersion property is directly related to the reservoir’s physical parameters. This provides theoretical support for using seismic data to obtain reservoir physical parameters.

Third, the main seismic attenuation mechanism in the seismic frequency band is mainly due to P-wave impact [12, 13]. For PP systems (P-wave source, P-wave reception), the P-wave velocity dispersion can be estimated from the poststack profile. Compared with FAVO inversion, the AVF (amplitude varies with frequency) inversion calculation efficiency will be greatly improved. At this time, the zero-angle analytical solution of the viscous acoustic wave can be used as the forward method.

Data Availability

The model data included in this study are available upon request by contact with the corresponding author. Due to confidentiality, actual data is only available upon request from the data provider (maintain secrecy).

Conflicts of Interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Acknowledgments

This work is financially supported by the China National Key Research Project (NO: 2019YFC0312003), National Natural Science Foundation of China (41774129, 41774131), and Science Foundation of China University of Petroleum, Beijing (2462018QZDX01).