Amplitude variation with incident angle inversion for fluid factor in the depth domain

The development of Pre-stack depth migration makes the imaging of the subsurface structure in the depth possible, which set a foundation for the study of amplitude variation with incident angle (AVA) inversion. This leads to the increasing demanding of the seismic inversion methods in the depth domain for guiding reservoir characterization. However, the conventional seismic inversion methods in the time domain are not suitable in the depth domain due to the seismic wavelet in the depth domain is depth-variant and depending on medium velocity. To address this issue, we proposed a pragmatic seismic inversion method for fluid factor in the depth domain with amplitude variation with incident angle gathers by using a true-depth wavelet on the process of seismic inversion. This wavelet is estimated by converting the time wavelet to the depth wavelet with seismic velocity. To guide the fluid discrimination, the proposed method directly estimates the fluid factor from the pre-stack seismic data and all the process of the method is implemented in the depth domain. To deal with the weak nonlinearity induced by the velocity, the Bayesian inference, the prior information and model constraint are introduced in this seismic inversion method. Tests on synthetic data show that the fluid factor can be well estimated reasonably even with moderate noise. The field data example illustrates the feasibility and efficiency of the proposed method in application and the estimated fluid factor and shear modulus are in good agreement with the drilling results.


INTRODUCTION
The development of pre-stack depth migration has propelled to the forefront in investigations of the imaging of subsurface structure in the depth domain. This renders the possibility of inversion for reservoir properties in depth domain. The development of reservoir models and subsequent interpretation requires the transfer of the reservoir properties from time domain to the depth domain. And well data are usually recorded in depth domain. Therefore, seismic inversion has dramatically shaped queries on depth domain in recent years to avoid the depth-to-time conversion and time-to-depth conversion. Zhang et al. [2006] presented a rock property parameter inversion method by using well logs and seismic data in the depth domain with a neural networks algorithm, which was a mapping function between seismic data and rock property parameters. Fletcher et al. [2012] suggested an approach to implement the seismic amplitude inversion directly in the depth domain and imaging by using the practical methods of optimization. Singh [2012] developed a deterministic inversion of seismic data in the depth domain. He utilized the pseudo-depth transformation with constant velocities to handle the wavelet estimation. Letki [2015] also documented a formulation of seismic inversion in the depth domain by using the gird of practical methods of optimization to capture the dip dependent spatial wavelet. Although this method could estimate more balanced acoustic impedance compared to the results in the time domain, the practical methods of optimization are not valid anymore in presence of strong boundaries and the results will not be reliable. Zhang et al. [2016] developed a seismic inversion method in depth domain by using the compressed sensing technique with a constant wavelet by using the seismic data in the depth domain with convolution model, and this may lead to a changing frequency and phase of the seismic data according to the variation of depth. Compared to conventional inversion in time domain, these methods obtained some advantages in the resolution and continuity. However, their methods avoid the estimation of wavelet or indirectly extract the wavelet in the depth domain. In this condition, the prevailing methods such as AVO inversion method cannot be applied to estimate the reservoir properties from pre-stack seismic data. In order to implement convolution and inversion in depth domain, a stable algorithm for estimating the wavelet is important.
The convolution model known as linear time-invariant system is usually preferred as the convolution of reflection coefficients and seismic wavelet and is widely accepted as forward solver in seismic inversion. This framework solely pertains to seismic modeling in time domain. However, seismic wavelet is depth-variant due to different velocities in depth domain [Hu et al., 2007;Singh, 2012]. To satisfy the linear depth-invariant, Hu et al. [2007] presented a formulation of interval velocity and thickness transformation to pseudo-depth by using constant velocity. However, pseudo-depth model was resampled in a constant depth interval, which might omit some reflection coefficients. To address this issue, Wei et al. [2017] proposed an analytic seismic wavelet to implement a convolution algorithm for seismogram synthesis in depth domain in terms of weighted superposition principle. On the basis of previous studies, we directly give the expression of depth wavelet and utilize it to implement the inversion in depth domain.
Considering the importance of the fluid discrimination with various fluid factors in seismic exploration and reservoir characterization, we try to utilize the wavelet in the depth domain to implement the amplitude variation with incident angle inversion for fluid factors from pre-stack seismic data. The prevailing methods for estimating the fluid factors in time domain from pre-stack seismic data renders the rapid development of technology for reservoir forecasting and fluid discrimination [Smith and Gidlow, 2000;Quakenbush et al., 2006]. Smith and Gidlow [1987] initially proposed to estimate the fluid factor and a pseudo Poisson ratio to guide lithology and fluid discrimination, based on a weighted stack method that used prestack seismic data. Biot [1941Biot [ , 1956 and Gassmann [1951] suggested the poroelasticity theories, which set the foundation of modeling physical rock-fluid interaction in situ and fluid discrimination in exploration geophysics. Russell et al. [2003] further estimated the fluid component by applying these theories and it was the combination of P-wave impedance, S-wave impedance, and a constant. They added the dry and saturated properties of in-situ reservoir rock to the expression for the P-wave velocity. That was a general extension of the amplitude variation with offset (AVO) approximate equation as in the Goodway et al. [1997]. The introduction of poroelasticity theory made the fluid discrimination more possible, at least in theory, when compared to the LMR technology. In this study, we focus on the estimation of the fluid term from pre-stack seismic data in depth domain. There are two mainly steps, depth wavelet estimation and AVA inversion in depth domain. The method can directly extract the fluid factor from seismic data and steadily estimate the depth wavelet in depth domain.
Tests on synthetic data show that the fluid factor is still estimated reasonable with moderate noise. A test on a real data set demonstrates that the estimated results are in good agreement with the results of drilling results.

SEISMIC WAVELET ESTIMATION IN THE DEPTH DOMAIN
The first issue of seismic inversion is the estimation of wavelet. The time wavelet and depth wavelet are equal in value, while sampling interval satisfy this relationship, where t ∆ is the time sampling interval, h ∆ is the depth sampling interval and v is the interval velocity. Given the frequency of time wavelet, the relationship between time and depth wavelet can be written as, We suppose that seismic velocity models in the time domain are Where, M is the number of sampling point and a is the reflection coefficient in depth domain.

AVA INVERSION IN THE DEPTH DOMAIN FOR FLUID FACTOR
Poroelasticity theories [Biot, 1956;Gassmann, 1951] have been increasingly utilized for fluid discrimination and reservoir prediction because the abundant amplitude variation with offset information of pre-stack seismic data. The relationship between velocity data and fluid properties is useful in reasonably explaining the underground fluid migration from seismic data [Mavko and Mukerji 1995]. From Biot-Gassmann theory, the P-wave velocity p V and S-wave velocity s V are expressed as: where sat K , µ , and sat ρ are the bulk modulus, shear modulus, and density of fluid-saturated rock, respectively. And where d K , s K , and f K are the equivalent bulk modulus of dry rock, rock matrix, and the saturated fluid, respectively. φ is the effective porosity of the rock. d µ is the shear modulus of dry rock. Assuming the saturated fluid consists of oil, gas, and water, where 0 S , g S , and w S are the saturation of oil, gas, and water, respectively. 0 K , g K , and w K are the bulk modulus of oil, gas, and water, respectively. Russell et al. (2011) defined the fluid/pore term f as, Combining Eq. (5), (6) and (8), the P-wave velocities may be rewritten as, With Eq. (9), it is easy to get, Simplifying the Eq. (10) where ( ) Convoluting the Eq. (12) with the depth wavelet as the forward solver, the fluid term can be estimated under Bayesian inversion scheme. In this paper, Cauchy distribution is utilized for a priori probability distribution of the model parameters [Sacchi and Ulrych 1995;Alemie and Sacchi 2011], to enhance the inversion resolution, and the likelihood function obey the Gaussian distribution [Buland and Omre 2003]. To weaken the strong degree of correlation among the three parameters in Eq. (12), an operator is utilized to decorrelate the parameters via the diagonalization of the covariance matrix. Furthermore, the low-frequency models from the well logs are incorporated into the objective function to enhance the stability and accuracy of the inversion.
Simplifying Eq. (12) as, where ( ) The posteriori probability distribution of the estimated model parameters And is the constraint coefficient for the f , µ , and density, respectively.
where U is the decorrelation matrix, and it is defined as follows. The covariance matrix r C of parameters to estimate is defined as,

Synthetic Example
We test the proposed amplitude variation with incident angle (AVA) inversion method in the depth domain with the synthetic seismic traces from real well logs which shows in the figure 2 and the sampling interval is 4 m. The seismic forward modeling is implemented using a convolution of the exact Zoeppritz equation with a depth Ricker wavelet, whose main frequency is 40 Hz. Figure 3 shows the synthetic common midpoint profile from real well logs (red) and the generated synthetic data in the depth domain (black). Figure 4 shows the original (in blue), initial (in green) and inverted (in red) fluid factor, shear modulus and density of well A with different S/N.
From Figure 4a, we see that all the parameters can be inverted well when there is no noise in the synthetic data even with smooth initial fluid factor, shear modulus and density. To further test the stability of the proposed inversion method, we add random Gaussian noise to the synthetic traces when the ratio of the signal to noise (S/N) is 5:1 and 2:1, respectively (as displayed in figure 4b and 4c). It is easy to demonstrate that the fluid factor can also be estimated well even with moderate noise, and the shear modulus can also be estimated reasonably. However, the density shows more bias from the true value due to the small contribution of density to P-wave reflection coefficients [Zong et al. 2015b]. Fluid factor and shear modulus are the main contribution to the reflectivity curve at small angles. Therefore, Fluid factor and shear modulus can be estimated more accurately than density.

Real data example
Real data is utilized to validate the feasibility of the proposed inversion method in the depth domain for application. True-amplitude processing has been implemented before the inversion of the real seismic data. The seismic data is processed by a contractor and the processing sequence is defined to ensure that the final pre-stack amplitudes should image the reflection strength of the subsurface interfaces as correctly as possible [Zong et al. 2015b]. Three partial angle-stack seismic profiles in time domain and depth domain through well A are shown in figure 6a-6c, the gamma ray of well A is shown in profile. The maximum incident angle is around 24 degrees for target areas around 2500m. A well was drilled to verify that two reservoirs develop.
One is the gas reservoir with water as labeled in red dashed circle at 2430m, and other is the gas reservoir as labeled in red dashed circle at 2560m in figure 6 and 7. Figure 7 displays the velocity field of the profiles and the lateral variation of velocity in work area is not obvious. Therefore, the differences between profiles of seismic data plotted in time-and depth-domains are small. But the seismic data will be different between time and depth domain when the area occur severe lateral velocity variation. In this condition, pre-stack depth migration image is an ideal seismic data for the interpretation of complex structure. We can directly implementing the inversion in the depth domain, which avoids the conversion from depth to time domain. Figure 5 displays the depth seismic wavelet of near (in black), middle (in blue), and far (in red) angle range in one depth point. These depth wavelet local at around 2463m. The initial models of the fluid factor and shear modulus are shown in figure 8, which were obtained by longitudinally smoothing the interpolation model which were generated by well log interpolation. The inverted results of fluid factor and shear modulus with smooth initial models by utilizing the proposed method are displayed in figure 9-10, respectively. The gas reservoir is labeled in red with water labeled in blue at 2430m, and the other is labeled in red at 2560m. From figure 9 and 10, we can see that the fluid factor shows anomalously low value in the gas reservoirs in well A and the shear modulus also shows low abnormality because the development of a large amount of cracks and pores leads to a weak anti-shear capacity of the framework in the gas reservoirs. Those results are consistent with the drilling and rock physics analysis results. Furthermore, curves of filtered well curves in the depth domain and the inverted results at known location are displayed in figure 11. The sampling interval of well curves is 4m, which is the same as seismic data. From those inversion results, at the zones where the reservoirs develop, the inverted f reveals anomaly with low value, which coincides with the results of rock physics analysis well.

CONCLUSIONS
We