# Nonlinear FAVO Dispersion Quantification Based on the Analytical Solution of the Viscoelastic Wave Equation

Authors: Yuanqiang Li; Jingye Li; Xiaohong Chen; Jian Zhang; Chen Zhou; Wei Tang
Journal: Geofluids
 (2020)

---

## 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.


---

## Body

## 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 [1–3]. 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 [10–12]. 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 [19–21]. 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 [29–31], 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 [33–36], 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 ofN horizontal layers, the total PP-reflectivity RPPω,p can be derived by the transfer function
(1)νn=Δ−RSPΔ−RSSΔRPPΔRPSΔRΔT,where Δ is the determinant associated with the system matrix, the value of which does not affect the calculation result. RPP is PP-reflectivity, RPS is PS-reflectivity, and others are similar. R=RPPRSSRPSRSP does not have physical meaning.To compute the total transfer functionv0 of the horizontal attenuated layers, we redefine the wave propagator matrix Qn for more convenient calculations:
(2)νn=Qnνn+1,where Qn=Tn+EnTn−; matrix En is the phase shift of reflection:
(3)En=diage−iωΔhqp+qs1e−iωΔhqp−qseiωΔhqp−qs1eiωΔhqp+qs.AndTn+, Tn−represents downward and upward interface energy distribution matrices, respectively. The matrix Tn+ contains 16 independent elements:
(4)Tn+=−p2+qpqs/μ−2pqp/μ−p2−qpqs/μp2−qpqs/μ−2pqs/μp2+qpqs/μiqs/v~sω20−iqs/v~sω2−iqs/v~sω20−iqs/v~sω2−ipΓ+2qpqs−4ip2qp−ipΓ−2qpqsipΓ−2qpqs−2iΓqs−ipΓ+2qpqs−ipΓ+2qpqs−2iΓqp−ipΓ−2qpqsipΓ−2qpqs−4ip2qs−ipΓ+2qpqs−iqp/v~sω20−iqp/v~sω2−iqp/v~sω20iqp/v~sω2−μΓ2+4p2qpqs−4μΓpqp−μΓ2−4p2qpqsμΓ2−4p2qpqs−4μΓpqs−μΓ2+4p2qpqs.The matrixTn− can be expressed by elements of matrix Tn+:
(5)Tn−=t61t51t31t31t21t11−t650−t45−t350−t15−t63−t51−t33−t33t12−t13t63−t51t33t33t21t13−t620−t42−t320−t12t61−t51t31t31−t21t11,where Δh is the thickness, qp=1/v~pω2−p2 and qs=1/v~sω2−p2 are the vertical slowness of the P-wave and S-wave, respectively. Γ=2p2−1/v~sω2, μ=ρv~sω2. The horizontal slowness is given as p=sinθ/v~pω. It should be noted that v~pω and v~sω, 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
(6)v~pω=νp1+1πQprlnωωr+i2Qpr,v~sω=νs1+1πQsrlnωωr+i2Qsr.νp and νs are reference velocities. Qpr and Qsr are quality factors. The parameter ωr 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), 1/πQrlnω/ωr is the dispersion term, which controls the phase of the wavelet, and i/2Qr 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 Qsr=Qpr.There is an elastic half-space below the base of the medium.νN can be expressed as
(7)νN=100000T.Snell’s lawp=sinθ/v~pω 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 v0:
(8)RPPω,θ=ν04ν01=RPPΔΔ.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 En.
(9)En=diage−iωΔh1/v~pω+1/v~sω1e−iωΔh1/v~pω−1/v~sωeiωΔh1/v~pω−1/v~sω1eiωΔh1/v~pω+1/v~sω.We then obtain the angle gather seismogramgt,θ by Inverse Fourier transform:
(10)gt,θ=12π∫−∞∞WωRPPω,θeiwtdω,where Wω 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 pointf0 by its expansion in a Taylor series. In this study, the nonlinear reflectivity function can be expressed as:
(11)Rf=Gmf,where G is the nonlinear operator; mf represents model parameters at different frequencies. The first-order approximation of the function in f0 is:
(12)Rf≈Gmf0+∂Gmf∂mfmf=mf0∂mf∂ff=f0f−f0,where Gmf0=Rf0. The derivative of the reflection coefficient with respect to the model is written as: Fr=∂Gmf/∂mfmf=mf0. Dispersion property, the derivative of the model parameter with respect to the frequency, can be written as: I=∂mf/∂ff=f0=IpIsIrhoT. Usually, we can use the P-wave dispersion property Ip=∂Vpf/∂ff=f0 as the fluid indicator. For prestack seismic data, equation (12) can be rewritten as:
(13)Rθ,t,f−Rθ,t,f0=Fr∗I∗f−f0.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 tracesSθ,t,f can be calculated using spectral decomposition techniques. Then using the frequency domain wavelet Wf, we get
(14)Sθ,t,f=Rθ,t,f∗Wf.According to the equations (13) and (14), the equation (15) is used to eliminate the overprint effect of wavelet.
(15)Sθ,t,f∗Wf0−Sθ,t,f0∗Wf=Wf0∗Wf∗f−f0∗Fr∗I.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,
(16)Rθ,t,f=Sθ,t,fweθ,f,where weθ,f=∑τS2θ,t,f0/∑τS2θ,t,f; τ represents the time window.When the operatorG 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 Fr is an essential step. Under the viscoelastic medium, applying the chain rule, we can obtain:
(17)∂R∂mn∗=∂R∂m~n∗∗∂m~n∗∂mn∗,where m~n∗=m~nv~pv~sρ denotes complex P-wave velocity, complex S-wave velocity, and density of layer n, respectively, mn∗=mnvpvsρ 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):
(18)∂R∂m~n∗=∂ν04/∂m~n∗ν01−ν04∂ν01/∂m~n∗ν012.∂ν0∂m~n∗=Q0Q1⋯∂Qn∂m~n∗⋯QN−1νN,∂Qn∂m~n∗=∂Tn+∂m~n∗EnTn−+Tn+∂En∂m~n∗Tn−+Tn+En∂Τn−∂m~n∗.∂m~nvp∂mnvp=∂m~nvs∂mnvs=1+1πQrlnωωr+i2Qr,∂m~nρ∂mnρ=1,so the derivatives of the viscoelastic medium can be expressed as:
(19)∂R∂mnvp=1+1πQrlnωωr+i2Qr∂R∂m~nv~p,∂R∂mnvs=1+1πQrlnωωr+i2Qr∂R∂m~nv~s,∂R∂mnρ=∂R∂m~nρ.The model parametermf0 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.Figure1 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.Figure 1
The nonlinear FAVO dispersion estimation workflow.


## 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 ofN horizontal layers, the total PP-reflectivity RPPω,p can be derived by the transfer function
(1)νn=Δ−RSPΔ−RSSΔRPPΔRPSΔRΔT,where Δ is the determinant associated with the system matrix, the value of which does not affect the calculation result. RPP is PP-reflectivity, RPS is PS-reflectivity, and others are similar. R=RPPRSSRPSRSP does not have physical meaning.To compute the total transfer functionv0 of the horizontal attenuated layers, we redefine the wave propagator matrix Qn for more convenient calculations:
(2)νn=Qnνn+1,where Qn=Tn+EnTn−; matrix En is the phase shift of reflection:
(3)En=diage−iωΔhqp+qs1e−iωΔhqp−qseiωΔhqp−qs1eiωΔhqp+qs.AndTn+, Tn−represents downward and upward interface energy distribution matrices, respectively. The matrix Tn+ contains 16 independent elements:
(4)Tn+=−p2+qpqs/μ−2pqp/μ−p2−qpqs/μp2−qpqs/μ−2pqs/μp2+qpqs/μiqs/v~sω20−iqs/v~sω2−iqs/v~sω20−iqs/v~sω2−ipΓ+2qpqs−4ip2qp−ipΓ−2qpqsipΓ−2qpqs−2iΓqs−ipΓ+2qpqs−ipΓ+2qpqs−2iΓqp−ipΓ−2qpqsipΓ−2qpqs−4ip2qs−ipΓ+2qpqs−iqp/v~sω20−iqp/v~sω2−iqp/v~sω20iqp/v~sω2−μΓ2+4p2qpqs−4μΓpqp−μΓ2−4p2qpqsμΓ2−4p2qpqs−4μΓpqs−μΓ2+4p2qpqs.The matrixTn− can be expressed by elements of matrix Tn+:
(5)Tn−=t61t51t31t31t21t11−t650−t45−t350−t15−t63−t51−t33−t33t12−t13t63−t51t33t33t21t13−t620−t42−t320−t12t61−t51t31t31−t21t11,where Δh is the thickness, qp=1/v~pω2−p2 and qs=1/v~sω2−p2 are the vertical slowness of the P-wave and S-wave, respectively. Γ=2p2−1/v~sω2, μ=ρv~sω2. The horizontal slowness is given as p=sinθ/v~pω. It should be noted that v~pω and v~sω, 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
(6)v~pω=νp1+1πQprlnωωr+i2Qpr,v~sω=νs1+1πQsrlnωωr+i2Qsr.νp and νs are reference velocities. Qpr and Qsr are quality factors. The parameter ωr 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), 1/πQrlnω/ωr is the dispersion term, which controls the phase of the wavelet, and i/2Qr 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 Qsr=Qpr.There is an elastic half-space below the base of the medium.νN can be expressed as
(7)νN=100000T.Snell’s lawp=sinθ/v~pω 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 v0:
(8)RPPω,θ=ν04ν01=RPPΔΔ.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 En.
(9)En=diage−iωΔh1/v~pω+1/v~sω1e−iωΔh1/v~pω−1/v~sωeiωΔh1/v~pω−1/v~sω1eiωΔh1/v~pω+1/v~sω.We then obtain the angle gather seismogramgt,θ by Inverse Fourier transform:
(10)gt,θ=12π∫−∞∞WωRPPω,θeiwtdω,where Wω 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 pointf0 by its expansion in a Taylor series. In this study, the nonlinear reflectivity function can be expressed as:
(11)Rf=Gmf,where G is the nonlinear operator; mf represents model parameters at different frequencies. The first-order approximation of the function in f0 is:
(12)Rf≈Gmf0+∂Gmf∂mfmf=mf0∂mf∂ff=f0f−f0,where Gmf0=Rf0. The derivative of the reflection coefficient with respect to the model is written as: Fr=∂Gmf/∂mfmf=mf0. Dispersion property, the derivative of the model parameter with respect to the frequency, can be written as: I=∂mf/∂ff=f0=IpIsIrhoT. Usually, we can use the P-wave dispersion property Ip=∂Vpf/∂ff=f0 as the fluid indicator. For prestack seismic data, equation (12) can be rewritten as:
(13)Rθ,t,f−Rθ,t,f0=Fr∗I∗f−f0.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 tracesSθ,t,f can be calculated using spectral decomposition techniques. Then using the frequency domain wavelet Wf, we get
(14)Sθ,t,f=Rθ,t,f∗Wf.According to the equations (13) and (14), the equation (15) is used to eliminate the overprint effect of wavelet.
(15)Sθ,t,f∗Wf0−Sθ,t,f0∗Wf=Wf0∗Wf∗f−f0∗Fr∗I.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,
(16)Rθ,t,f=Sθ,t,fweθ,f,where weθ,f=∑τS2θ,t,f0/∑τS2θ,t,f; τ represents the time window.When the operatorG 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 Fr is an essential step. Under the viscoelastic medium, applying the chain rule, we can obtain:
(17)∂R∂mn∗=∂R∂m~n∗∗∂m~n∗∂mn∗,where m~n∗=m~nv~pv~sρ denotes complex P-wave velocity, complex S-wave velocity, and density of layer n, respectively, mn∗=mnvpvsρ 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):
(18)∂R∂m~n∗=∂ν04/∂m~n∗ν01−ν04∂ν01/∂m~n∗ν012.∂ν0∂m~n∗=Q0Q1⋯∂Qn∂m~n∗⋯QN−1νN,∂Qn∂m~n∗=∂Tn+∂m~n∗EnTn−+Tn+∂En∂m~n∗Tn−+Tn+En∂Τn−∂m~n∗.∂m~nvp∂mnvp=∂m~nvs∂mnvs=1+1πQrlnωωr+i2Qr,∂m~nρ∂mnρ=1,so the derivatives of the viscoelastic medium can be expressed as:
(19)∂R∂mnvp=1+1πQrlnωωr+i2Qr∂R∂m~nv~p,∂R∂mnvs=1+1πQrlnωωr+i2Qr∂R∂m~nv~s,∂R∂mnρ=∂R∂m~nρ.The model parametermf0 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.Figure1 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.Figure 1
The nonlinear FAVO dispersion estimation workflow.


## 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 Figure2, 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.Figure 2
Unequal thickness interlayer models (P-wave velocity, S-wave velocity, and density).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. Figure3 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.Figure 3
Angle gathers comparison based on the VAE and the ASM. Part1 and Part2 respectively, correspond to angle gathers by VAE, ASM. Part3 is the residual between VAE and ASM; red box area highlights attenuation and transmission loss. Blue box area highlights internal multireflections.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 Figure3) 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.Figure 4
(a) VAE reflection coefficient variation curves with angle at different frequencies. (b) ASM reflection coefficient variation curves with angle at different frequencies.
(a)(b)


### 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 ofvp, vs, ρ, and Q are, respectively, shown in Figure 5. The vp, vs, and ρ data were taken from field well logging after Backus-averaging, and the Q value was calculated by empirical formula [46, 47] Q=avpb (in this actual study area, a=5.7 and b=2.25). 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°).Figure 5
P-wave velocity, S-wave velocity, density, and Q value curves.Figure 6
The synthetic angle gathers were generated by ASM using well logs.Following the workflow shown in Figure1, 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: Hf=exp−f∫π/Qvdl on seismic amplitude [38]. So we can use a low-frequency trend of the Q value instead of the real Q value curve to perform compensation or inversion [48]. Using the low-frequency trend of the Q 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 Q 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.Figure 7
(a) Inversion results by PVWI using low-frequency trend of Q value curve. (b) Inversion results by single-interface AVO inversion. Blacks are true parameters, reds are inverted parameters, and blues are initial parameters and low-frequency trend of Q.
(a)(b)After spectral decomposition and spectral balance, the P-wave dispersion was estimated by the nonlinear frame. As shown in Figure8(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.Figure 8
(a) The comparison between real dispersion and estimated P-wave dispersion by the nonlinear FAVO. The black line indicates the real dispersion. The red line represents estimated P-wave dispersion. (b) The comparison of the estimated P-wave dispersion between the nonlinear FAVO and the single-interface FAVO inversion (Wilson-based). The black line indicates the estimated P-wave dispersion by single-interface FAVO inversion. The red line represents estimated P-wave dispersion by the nonlinear FAVO.
(a)(b)


## 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 Figure2, 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.Figure 2
Unequal thickness interlayer models (P-wave velocity, S-wave velocity, and density).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. Figure3 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.Figure 3
Angle gathers comparison based on the VAE and the ASM. Part1 and Part2 respectively, correspond to angle gathers by VAE, ASM. Part3 is the residual between VAE and ASM; red box area highlights attenuation and transmission loss. Blue box area highlights internal multireflections.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 Figure3) 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.Figure 4
(a) VAE reflection coefficient variation curves with angle at different frequencies. (b) ASM reflection coefficient variation curves with angle at different frequencies.
(a)(b)


## 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 ofvp, vs, ρ, and Q are, respectively, shown in Figure 5. The vp, vs, and ρ data were taken from field well logging after Backus-averaging, and the Q value was calculated by empirical formula [46, 47] Q=avpb (in this actual study area, a=5.7 and b=2.25). 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°).Figure 5
P-wave velocity, S-wave velocity, density, and Q value curves.Figure 6
The synthetic angle gathers were generated by ASM using well logs.Following the workflow shown in Figure1, 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: Hf=exp−f∫π/Qvdl on seismic amplitude [38]. So we can use a low-frequency trend of the Q value instead of the real Q value curve to perform compensation or inversion [48]. Using the low-frequency trend of the Q 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 Q 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.Figure 7
(a) Inversion results by PVWI using low-frequency trend of Q value curve. (b) Inversion results by single-interface AVO inversion. Blacks are true parameters, reds are inverted parameters, and blues are initial parameters and low-frequency trend of Q.
(a)(b)After spectral decomposition and spectral balance, the P-wave dispersion was estimated by the nonlinear frame. As shown in Figure8(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.Figure 8
(a) The comparison between real dispersion and estimated P-wave dispersion by the nonlinear FAVO. The black line indicates the real dispersion. The red line represents estimated P-wave dispersion. (b) The comparison of the estimated P-wave dispersion between the nonlinear FAVO and the single-interface FAVO inversion (Wilson-based). The black line indicates the estimated P-wave dispersion by single-interface FAVO inversion. The red line represents estimated P-wave dispersion by the nonlinear FAVO.
(a)(b)


## 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. Figure9 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°).Figure 9
(a) The wavelet estimated based on statistical. (b) The wavelet spectrum for selecting the key frequency points for calculation.
(a)(b)First, we extract the wavelet (Figure10(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 f0=22 Hz. Then, we selected the stagnation and inflection points of the amplitude spectrum function around f0 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 f1,f2,f3,f4=13,18,28,34, respectively.Figure 10
Seismic poststack section of the interest.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 formulaQ=avpb, where vp is from the initial model. Through the regional attenuation background survey, the coefficient was corrected using logging data (a=5.9, b=2.18). 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.Figure 11
The initial model of P- and S-wave velocities, density, and the low-frequency trend of the Q: (a) P-wave velocity; (b) S-wave velocity; (c) density; (d) quality factor.
(a)(b)(c)(d)Figure 12
Inverted parameters by PVWI including results of P-wave velocity (a), S-wave velocity (b), and density (c).
(a)(b)(c)Figure 13
Inverted parameters by single-interface AVO inversion including results of P-wave velocity (a), S-wave velocity (b), and density (c).
(a)(b)(c)After spectral decomposition and spectral balancing, we estimated P-wave dispersion by the nonlinear FAVO (Figure14(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.Figure 14
The comparison of estimated P-wave dispersion between the (a) single-interface FAVO inversion and (b) the nonlinear FAVO. The red ellipse represents the target reservoir.
(a)(b)Figure 15
The spectral decomposition frequency sections of the poststack data: (a) 16 Hz; (b) 22 Hz; (c) 28 Hz; (d) 34 Hz.
(a)(b)(c)(d)


## 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.



---

*Source: 7616045-2020-06-05.xml*