#### Abstract

Based on the physical quantity of log data, the accurate identification of oil- and gas-bearing properties may be caused by the prestack inversion of fluid prediction, which will affect the success rate of exploration and development. Prestack data contain more information of amplitude and frequency. Using the frequency-dependent viscoelastic impedance equation and Bayesian inversion framework, the objective function of frequency-dependent elastic impedance inversion can be established to realize the frequency-dependent impedance inversion at different angles. According to the elastic impedance equation of the frequency-varying viscoelastic fluid factor, the relationship between elastic impedance and the frequency-dependent viscoelastic fluid factor is established, and the prestack seismic inversion method of the frequency-dependent viscoelastic fluid factor is studied. However, one of the important factors easily neglected is that we have been using logging data to establish fluid-sensitive parameters and the lithophysical version for fluid identification, so there are differences between logging and seismic frequency bands for fluid identification. The indicator factors with higher sensitivity to fluid can be selected by laboratory measurements. This article applies this method on Luojia oilfield data and verifies this method with log interpretation results, based on the sample of rock physics obtained in a low-frequency rock physics experiment; the technique of dispersion and fluid-sensitive parameters is studied, and the fluid prediction technology of a multifrequency band rock physics template is adopted, which can build the relationship between rock physical elastic parameters and fluid properties by the multifrequency broadband impedance method.

#### 1. Introduction

Seismic waves are affected by the pore fluid properties (fluid type, saturation, etc.) and the characteristics of rock skeleton (porosity, bedrock modulus, etc.) when they propagate in the pore medium containing fluid, which will produce seismic wave attenuation and dispersion. Scholars have made a great deal of research on fluid property theory about dispersion and seismic wave propagation in the pore medium. Russell et al. [1] defined the parameter as the fluid factor; then, Russell et al. [2] presented the fluid modulus with the rock physics theory of the porous elastic medium. The construction of the fluid factor works as the key problem in seismic fluid discrimination [3, 4]. The attenuation and dispersion of seismic waves in a porous medium are mainly caused by the equivalent flow of the pore fluid. According to the scale of the fluid, it can be divided into a macroscopic scale [5], microscopic scale, and mesoscopic scale fluid flow. The law of mass conservation of fluid in the pore medium, generalized Newton’s theorem, and elastic constitutive relation are deduced by Biot. The Lagrangian equation was deduced in the Biot theory without considering inertial coupling between the fluid and solid. The Biot theory contains biphasic conditions and a multiphase media elastic wave theory. The fluid attenuation mechanism is that when seismic waves propagate in a fluid-containing medium; wave crests and troughs are balanced, which result in friction, that is, fluid flow on a macroscopic scale. The fluid attenuation scale is equivalent to that of the seismic wavelength. Under the assumption that the pore medium contains a small amount of gas, for the first time, the theory of Biot elastic fluctuation at the macroscopic scale and the attenuation mechanism of “jet flow” at the microscopic scale are combined, which established the theory of BISQ elastic fluctuation. The relation between predicted permeability, seismic wave attenuation, and dispersion velocity is consistent with the experimental measurement [6]. However, the theory of attenuation and dispersion of seismic waves only take effect at macroscopic and microscopic scales, which cannot explain the attenuation and dispersion of seismic waves in the seismic frequency band in detail. White and others first put forward the concept of the mesoscopic scale petrophysical theory [7], and two White models based on the Biot theory framework were established. One is the model of alternating patchy saturation of horizontally layered media containing two different fluids; the other is the saturation model of spherical bubbles containing periodic arrangement in saturated media. The model successfully predicted the attenuation and dispersion of the seismic frequency band when a seismic wave propagates in the pore medium. The White model is mainly used to study the nonuniform infiltration of heterogeneous fluids into porous rocks, which is accepted as the patchy saturation model.

Then, Dutta and Ode use a more rigorous elastic mechanics theory to solve the seismic wave dispersion and attenuation of the spherical periodic patch saturation model. The influence of water saturation and spherical radius on attenuation and dispersion of seismic waves is analyzed. Dutta and Seriff further improve it; the degree of agreement between the White model and the Gassmann theory at the zero frequency limit is improved [8]. Gurevich and Lopatnikov derived a plane wave modulus with random variation of characteristic thickness in the White model using a statistical smoothing theory, finding that the mesoscopic attenuation mechanism in random media has a larger frequency band range. The reason is that the change of medium thickness characteristic will differentiate the seismic wave attenuation. By studying the permeability and fluid pressure, Shapiro and Müller considered that the mesoscopic fluid flow is the main mechanism of seismic frequency band energy attenuation and velocity dispersion [9]. Johnson improved the bubble White model by introducing two new geometric parameters, which have irregular shape. Considering the effect of jet flow between adjacent pores and individual microscale fractures when seismic waves propagate in microfractures and pore media, the effect of large-scale directional (centimeter-meter) cracks is considered. The theory of elastic fluctuation of anisotropic fracture-pore microstructure is established; fluid type, saturation, porosity, and fracture density can also be analyzed by the theory [10]. The fracture development scale has great influence on seismic wave attenuation and dispersion. Gurevich established a coin-shaped fissure connecting the pore unit, skeleton bulk modulus, and shear modulus with frequency variation. According to Gurevich’s research, the low frequency result is consistent with the Gassmann theory and the high frequency result is consistent with the Mavko-Jizba theory. Rubino used a simple numerical simulation process to link the wave equation with the mesoscopic scale attenuation mechanism; at the same time, he showed the attenuation and dispersion results.

In addition to the inhomogeneous petrophysical model caused by infiltration with a multiphase fluid and heterogeneous saturation, the heterogeneity of the rock skeleton structure, pore shape, and distribution can also lead to inhomogeneity, among which the most representative is the double pore medium model. A number of scholars have studied the seismic response of the double pore medium model. Berryman first studied the definition and calculation method for each parameter of the equation. The elastic wave propagation and attenuation in a double pore double permeability medium were further studied later. Based on the volume average approximation, Pride derived the fluid seepage control equation in an isotropic double pore medium. In his research, there are two fluid seepage modes. Pride studied the attenuation and velocity dispersion of seismic waves in detail based on three longitudinal wave attenuation models simulated by a double pore medium model. It is pointed out that the mesoscale inhomogeneous porous solid skeleton model and the heterogeneous fluid inhomogeneous saturation model can produce longitudinal wave dispersion in the seismic frequency band range. The wave equations of Biot phase-to-phase displacement and Biot absolute displacement are derived. The numerical simulation of a mesofluid is realized based on a pseudospectral method.

Many seismic inversion methods in the transform domain and time domain inversion methods have developed rapidly because of their excellent antinoise ability and stability. In recent years, time domain seismic inversion has developed rapidly in two aspects: in view of the inherent band-limited characteristics of seismic signals, prior information of the regularization constraint can reduce the space of model parameters and retrieve elastic parameter information. Considering that the reflection coefficient sequence of underground media is sparse, the constraints such as L0/L1 norm, minimum entropy criterion, Huber norm, and Cauchy criterion are usually used as regular terms in the inversion process. Zhang built a wedge-shaped reflection coefficient dictionary. Nguyen studied the matching tracking reflection coefficient inversion and orthogonal matching tracking prestack seismic inversion methods, respectively. In recent years, scholars have also done a lot of research on Bayesian sparse seismic prestack AVO inversion methods. First of all, for low frequency model constraints, Yin and Zong studied point constraints, petrophysical model constraints, and smoothing model constraints for prestack seismic inversion, which improve the stability of the inversion. Gelderblom proposed a model inversion algorithm, which can impose a geological prior constraint in the model. Hamid developed a multichannel seismic impedance inversion method based on stratigraphic lateral constraints. The model test and data processing showed that this method has good effect in suppressing random noise and fidelity of formation boundary. During the recursive inversion of the electromagnetic model, Kim deduced the iterative weighted least square algorithm under the constraint of a model parameter. The reliability of the inversion results is improved by adding the low bound and high bound constraints of the model parameters in the inversion process. Zhu and Engel proposed a recursive inversion algorithm based on linear inequality constraints for the regularization least square inverse problem. In actual data processing, using logging data, drilling data, and geological stratification information, it is often possible to obtain background information on the model parameters of underground media. Therefore, how to use prior constraint information to improve the reliability of inversion results has become a hot and difficult point.

Considering the stability of multiparameter inversion, the nonlinear seismic inversion algorithm proposed a fast random inversion based on FFT-MA spectrum simulation. Compared with the conventional linear inversion algorithm, the resolution is improved obviously. However, considering the computation time consumption, stability, and reliability, the nonlinear inversion algorithm has not been widely used in practical processing. Although the Zoeppritz exact equation is a forward operator, which will greatly increase the degree of nonlinearity of the inverse problem, the process will increase the computation and time consumption of inversion; also, the stability and reliability of the inversion results are reduced; as a result, this direction has not been widely used in practical data processing.

Compared with the conventional time domain seismic inversion, the inversion of the reflection coefficient in the frequency domain has higher formation resolution [11] and is applied to distinguish a thin layer thickness interpretation within the tuning thickness. The spectral inversion algorithm is based on the parity decomposition of the reflection coefficient. Then, the time-window seismic data matching and the spectrum of the stationary seismic wavelet are obtained by using the short-time Fourier transform. Secondly, the inversion target functional is constructed based on the least square strategy in the frequency domain. The inversion method of the high resolution reflection coefficient is realized. Rubino developed a prestack AVA spectral inversion method based on the fast simulation algorithm, which also shows that the frequency domain inversion is helpful to distinguish thin layers in tuning thickness.

All the above-mentioned petrophysical models attempt to explain attenuation and dispersion phenomena. At the same time, it lays a key theoretical foundation [12] for the application of attenuation and frequency variation attributes in fluid identification. Taner et al. found that the frequency below the reservoir moves to a low frequency direction when an instantaneous frequency attribute is extracted by complex seismic channel analysis, which can be called as a low frequency shadow [13]. Castagana used an instantaneous spectral decomposition technique to detect low frequency shadow characteristics related to a reservoir hydrocarbon, which shows a good application prospect of the low frequency shadow in seismic data analysis and interpretation. However, the formation mechanism of the low frequency shadow is not clearly given. Geophysicists have been focusing on the causes of “low frequency shadows” and numerical simulations.

In this paper, the viscoelastic fluid factor inversion equation and target functional are derived. The sensitivity of the viscoelastic fluid factor is analyzed by the experimental model, and the method is applied in the Luojia oilfield. The results are compared with the conventional fluid factor to verify that the method can effectively improve the accuracy of fluid detection.

#### 2. Materials and Methods

The time and frequency statistics of seismic signals vary with time during the propagation of underground complex media; that is, the seismic signal itself is a nonstationary or time-varying signal. In order to reveal more time and frequency information contained in seismic data, scholars have developed many time-frequency analysis methods for nonstationary seismic signals, such as the short-time Fourier transform (STFT), continuous wavelet transform (CWT), S transform (ST), Wigner-Ville distribution (WVD), and time-frequency characterization (MP). Stratigraphic data interpretation, thin layer thickness estimation, and reservoir fluid detection have been widely used recently.

However, at present, among many spectral decomposition algorithms, the time-frequency resolution of matching tracking spectral decomposition algorithm is the highest; the classical matching tracking adaptive decomposition algorithm based on Gabor time-frequency atoms is proposed by Mallat and Zhang. Aiming at the problem of huge computation of the matching tracking algorithm, scholars put forward a large number of fast algorithms; the transition from the greedy algorithm to the fast dynamic algorithm proposed a complex domain double parameter fast matching tracking algorithm based on the Morlet wavelet dictionary. Then, the number of dynamic matching search parameters is reduced; as a result, calculation efficiency is improved obviously.

Velocity dispersion and energy decay usually occur during the propagation of seismic signals in underground complex media. Research shows that this property becomes particularly prominent in hydrocarbon-bearing reservoirs. In order to make full use of seismic wave attenuation and dispersion properties to describe reservoir characteristics, with the geophysical exploration progress, it is necessary to develop an effective method to predict reservoir distribution and evaluate the reservoir hydrocarbon. In the face of a more complex geological background, a viscoelastic frequency-dependent fluid factor inversion method based on broadband impedance is developed in combination with amplitude and frequency information of seismic data, which provides reliable basis for complex fluid identification and well location selection. When the seismic wave propagates in the oil and gas reservoir, a dispersion phenomenon will occur [14–17]. How to use the seismic wave attenuation and the dispersion attribute to describe reservoir fluid characteristics becomes one of the key research points [18]. Wu combined the RSPWVD algorithm to achieve the frequency AVO inversion method; Zhang et al. studied the calculation method of the longitudinal wave velocity dispersion attribute based on poststack inversion [19]. Based on the Shuey approximate formula, the concept of the velocity dispersion gradient of the longitudinal wave is proposed; Cheng et al. applied the frequency conversion AVO fluid identification technology to actual seismic data processing [20].

##### 2.1. Methodology

Therefore, considering the attenuation of the medium, the amplitude and frequency anomaly information contained in the prestack seismic data should be fully utilized to construct the frequency variable viscoelastic fluid factor, which is helpful to predict the distribution state of the fluid and further improve the reliability of the reservoir description.

The constitutive matrix of the inhomogeneous medium can be used as the superposition of the homogeneous background, isotropic disturbance, and viscosity disturbance:

is the background medium, is the isotropic medium, and is the fluid factor.

The viscoelastic fluid factor can be defined as the sum of the elastic part and the fluid factor:

The viscoelastic fluid factor can be expressed as

Substituting equation (3) into equation (1) and simplifying the equation:

where

Using Taylor expansion at a certain frequency to solve [5]

After using the sparse pulse inversion method developed by this project to obtain the data body of the frequency variable elastic impedance, the corresponding frequency variable elastic parameters should be extracted based on the relationship between the frequency variable elastic impedance and the frequency variable viscoelastic fluid factor. According to the equation of the viscoelastic frequency-varying elastic impedance derived in the previous chapter, there is a nonlinear exponential relationship between the viscoelastic frequency-varying fluid factor and impedance. In order to simplify the solution method, the formula is logarithmic transformed into a linear form. VEI is the broadband impedance, and and are broadband impedance-related coefficients. According to equation (6), two or more than two angle elastic impedances are needed. Under the influence of the velocity ratio of longitudinal and shear waves and different frequency selection, the equations are solved directly in The values of coefficients and applicable to the study area can be obtained by solving the above matrix coefficients corresponding to different angles and frequencies, respectively, which can be obtained by calculating for different incident angles and different frequency choices. By substituting coefficients into the above equations, the following equations can be obtained:

##### 2.2. Model Tests

In Figure 1, the indicator factors with higher sensitivity to fluid can be tested by laboratory model measurements. The fluid sensitivity of frequency variation attribute depends on the difference of dispersion degree caused by different fluid types in the seismic frequency band. Compared with the elastic medium fluid factor, the viscoelastic fluid factor can better simulate the propagation process of the seismic wave in the complex medium. The dispersion degree of eight sensitive parameters including the viscoelastic fluid factor was tested by laboratory model data. Figure 1(a) represents the dispersion test results of normalized sensitive parameters with frequency change. Then, the change rate parameters are obtained by the dispersion results of eight sensitive parameters. It can be seen that the sensitivity order is the highest . The viscoelastic fluid factor is tested according to the model. The validity of the fluid factor construction process is proved by the model test.

**(a) Dispersion of eight sensitive parameters**

**(b) Sensitivity of eight parameters**

Viscoelastic fluid factor inversion can be summed up into three parts as in Figure 2:
(1)*Multiscale decomposition*: selecting the reference frequency and a series of dominant frequencies; based on the cwt scale decomposition method, a series of equal frequency band seismic data is obtained at the dominant frequency, and the variation characteristics of each equal frequency profile with frequency are analyzed.(2)*Broadband impedance inversion*: selecting the reference frequency and a series of dominant frequencies; based on the cwt scale decomposition method, a series of equal frequency band seismic data is obtained at the dominant frequency, and the variation characteristics of each equal frequency profile with frequency are analyzed.(3)*Frequency-dependent viscoelastic fluid factor extraction*: according to the elastic impedance curve of the well side channel and the viscoelastic frequency variable fluid factor curve, the matrix equation of the dispersion property is constructed, and the weighted coefficient of the frequency angle is predicted.

#### 3. Results and Discussion

According to the characteristic equation of the viscoelastic prestack frequency variable reflection, combined with the results of the frequency variable viscoelastic fluid factor construction, a frequency variable viscoelastic fluid factor estimation strategy based on prestack seismic inversion is proposed. The frequency variable viscoelastic fluid factor is an extension of the conventional frequency variable seismic inversion. The mapping equation and target function are constructed by combining Bayesian estimation and a prior model regularization term to improve the reliability and noise resistance of the frequency variable viscoelastic property inversion. The application of the viscoelastic frequency variable fluid factor in reservoir fluid identification in a practical work area is helpful to improve the prediction accuracy of the reservoir pore fluid spatial distribution range.

In order to further verify the effectiveness of this method, this article uses the Luojia oilfield to verify frequency-dependent inversion based on broadband impedance. To verify the practicability of the frequency variable viscoelastic fluid factor method, based on the theoretical model to verify the feasibility of the method, the authors deal with the actual seismic data of the exploration area in the eastern part of China. The actual seismic profile is superimposed on the partial angle of the l m well as shown in Figure 3, where the frequency is 12 Hz, 18 Hz, 24 Hz, and 40 Hz. The sampling interval is 2 ms; the black curve in the section shows the resistivity of the logging curve (the shallow bilateral resistivity logging results are abnormally high in the reservoir development position), as indicated by the black arrow in the figure. In the 15 Hz profile, there are obvious “low frequency shadow” anomalies under the reservoir (white elliptical region). With the increasing frequency component, the low frequency strong energy mass under the reservoir gradually disappears. The phenomenon of the “low frequency shadow” provides evidence for the inversion of the variable viscoelastic fluid factor, but there are many “false bright spots” in the single frequency profile, so it will interfere with the pore fluid identification only by the “low frequency shadow.” The frequency-varying viscoelastic fluid factor inversion profile highlights reservoir characteristics and reduces the false appearance of single frequency attribute fluid identification.

Figure 4 shows the fluid factor profile compared with the viscoelastic fluid factor profile of wells A and B in the Luojia oilfield. According to well A within the black dotted circle, the oil layer is predicted by the fluid factor profile; however, the viscoelastic fluid factor profile has no oil response, which can verify the reliability of the viscoelastic fluid factor. As for well B at the red arrow, when the fluid layer is drilled, the conventional fluid prediction section is the oil layer, but the reservoir is not developed in the viscoelastic fluid factor prediction section, and the prediction result is accurate. The results of these two wells effectively test the effectiveness of the method. Compared with the results of the conventional fluid factor inversion, the frequency variable viscoelastic fluid can effectively identify the development position of oil reservoirs in seismic inversion. The dispersion degree of the viscoelastic fluid factor will be greatly improved by the sensitivity of reservoir indication, which indicates that the properties of viscoelastic fluid properties can better characterize the hydrocarbon properties in the reservoir than those of conventional dispersion fluid factors.

#### 4. Conclusions

The prestack inversion of the frequency variable viscoelastic fluid factor is an extension of the conventional frequency viscoelastic fluid factor inversion, which takes into account the influence of the underground medium viscosity on the seismic wave viscoelastic fluid factor reflection characteristics. The viscoelastic frequency-dependent viscoelastic fluid factor is approximately the key theoretical basis of the frequency-dependent viscoelastic fluid factor inversion. The viscoelastic approximation equation derived in this paper is compared with the exact Zoeppritz equation by the theoretical model, which verifies that the frequency-dependent viscoelastic fluid factor approximation formula has high accuracy. In terms of the inversion algorithm, the inversion mapping matrix and target functional are constructed by combining the Bayesian formula and the prior model regularization term. The model test shows the feasibility and noise resistance of the proposed method in extracting viscoelastic dispersion attributes. Finally, the frequency-dependent viscoelastic fluid factor inversion method is applied to the actual seismic data processing and reservoir fluid identification test. By comparing with the application effect of single frequency seismic profile and conventional frequency variable fluid factor, it is found that the inversion result of the frequency variable viscoelastic fluid factor has less “false bright spot” interference and will have higher prediction accuracy in reservoir pore hydrocarbon detection. The viscoelasticity of the medium can better simulate the seismic wave propagation in the complex medium according to the research in this article. Viscoelastic property has close relation with the fluid reservoir, which is a key attribute from qualitative to quantitative. The frequency changing reservoir characteristic is the foundation to establish the viscoelastic property model. The viscoelastic medium can better simulate the wave propagation process in a complex medium, in which the viscoelastic fluid factor can be extracted based on the relation between amplitude and frequency. Also the study of the frequency-dependent viscoelastic fluid factor prestack seismic inversion method is helpful to improve the reliability of fluid identification.

#### Data Availability

All data included in this study are available upon request through contact with the corresponding author.

#### Conflicts of Interest

The authors declare no conflict of interest.

#### Acknowledgments

We gratefully acknowledge the Wave Phenomena and Intelligent Inversion Imaging (WPI) group for the financial support. The authors thank the sponsors of the WPI group for their financial support and help. WPI’s research works are also financially supported by the National Key R&D Program of China (2018YFA0702503, 2019YFC0312004), National Natural Science Foundation of China (41774126), and National Science and Technology Major Project of China (2016ZX05024-001, 2016ZX05006-002).