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

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, Swave 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.


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][2][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 depen-dent) 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][11][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 . 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][20][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][30][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][34][35][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.

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 N horizontal layers, the total PP-reflectivity R PP ðω, pÞ 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. R PP is PP-reflectivity, R PS is PS-reflectivity, and others are similar. jRj = R PP R SS R PS R SP does not have physical meaning.
To compute the total transfer function v 0 of the horizontal attenuated layers, we redefine the wave propagator matrix Q n for more convenient calculations: where Q n = T + n E n T − n ; matrix E n is the phase shift of reflection: And T + n , T − n represents downward and upward interface energy distribution matrices, respectively. The matrix T + n contains 16 independent elements: The matrix T − n can be expressed by elements of matrix T + n : where Δh is the thickness, q p = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1/ṽ p ðωÞ 2 − p 2 q and q s = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1/ṽ s ðωÞ 2 − p 2 q are the vertical slowness of the P-wave and S-wave, respectively. Γ = 2p 2 − 1/ṽ s ðωÞ 2 , μ = ρṽ s ðωÞ 2 . The horizontal slowness is given as p = sin θ/ṽ p ðωÞ. It should be noted thatṽ p ðωÞ andṽ 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 ν p and ν s are reference velocities. Q pr and Q sr 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/πQ r Þ ln jω/ω r j is the dispersion term, which controls the phase of the wavelet, and i/2Q r 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 Q sr = Q pr .
There is an elastic half-space below the base of the medium. ν N can be expressed as Snell's law p = sin θ/ṽ 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 v 0 : 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 E n .
We then obtain the angle gather seismogram gðt, θÞ by Inverse Fourier transform: where WðωÞ is the source wavelet in the frequency domain.

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 f 0 by its expansion in a Taylor series. In this study, the nonlinear reflectivity function can be expressed as: where G is the nonlinear operator; mðf Þ represents model parameters at different frequencies. The first-order approximation of the function in f 0 is: where Gðmðf 0 ÞÞ = Rð f 0 Þ. The derivative of the reflection coefficient with respect to the model is written as: Fr = ∂Gðmð f ÞÞ/∂mðf Þj mð f Þ=mð f 0 Þ . Dispersion property, the derivative of the model parameter with respect to the frequency, can be written as: I = ∂mð f Þ/∂f j f =f 0 = ½I p I s I rho T . Usually, we can use the P-wave dispersion property I p = ð∂V p ðf Þ/∂f Þ j f =f 0 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 Sðθ, t, f Þ can be calculated using spectral decomposition techniques. Then using the frequency domain wavelet Wð f Þ, 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 When the operator G 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: wherem * n =m fṽ pṽs ρg n denotes complex P-wave velocity, complex S-wave velocity, and density of layer n, respectively, m * n = m fv p v s ρg n 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 mð f 0 Þ 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  4 Geofluids frequency points for calculation, the velocity dispersion properties can be estimated.

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 atten-uation 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.

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 v p , v s , ρ, and Q are, respectively, shown in Figure 5. The v p , v s , and ρ data were taken from field well logging after Backus-   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: QvÞdlÞ on seismic amplitude [38]. So we can use a lowfrequency 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   Figure 6: The synthetic angle gathers were generated by ASM using well logs. 7 Geofluids assumed to be the real seismic data in which the propagation effect is not removed, and then, we performed singleinterface 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.   Blacks are true parameters, reds are inverted parameters, and blues are initial parameters and low-frequency trend of Q.

Geofluids
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 highfrequency 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.

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°a nd 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 Amp Figure 9: (a) The wavelet estimated based on statistical. (b) The wavelet spectrum for selecting the key frequency points for calculation. 9 Geofluids Figure 10(b). From the wavelet spectrum, we can estimate the wavelet dominant frequency and use it as the reference frequency f 0 = 22 Hz. Then, we selected the stagnation and inflection points of the amplitude spectrum function around f 0 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 ½ f 1 , f 2 , f 3 , f 4 = ½13, 18, 28, 34, respectively.
After horizon picking, the initial model of P-and Swave velocities and density were generated by interpolation with four wells. The quality factor was determined by the empirical formula Q = av p b , where v p 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 lowfrequency 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 singleinterface 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 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 dril-ling verification confirms that the red ellipse area is a favorable hydrocarbon reservoir.

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 singleinterface 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 lowfrequency 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 frequencydependent 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 zeroangle 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.