Seismic pre-stack AVA inversion scheme based on lithology constraints

Seismic pre-stack AVA inversion using the Zoeppritz equation and its approximations as a forward engine yields P- and S-wave velocities and density. Due to the presence of seismic noise and other factors, the solution to seismic inversion is generally ill-posed and it is necessary to add constraints to regularize the algorithm. Moreover, since pre-stack inversion is a nonlinear problem, linearized optimization algorithms may fall into false local minima. The simulated annealing (SA) algorithm, on the other hand, is capable of finding the global optimal solution regardless of the initial model. However, when applied to multi-parameter pre-stack inversion, standard SA suffers from instability. Thus, a nonlinear pre-stack inversion method is proposed based on lithology constraints. Specifically, correlations among the elastic parameters are introduced to establish constraints based on a Bayesian framework, with special intention of mitigating the ill-posedness of the inversion problem as well as addressing the lithological characteristics of the formations. In particular, to improve the stability, a multivariate Gaussian distribution of elastic parameters is incorporated into the model updating the SA algorithm. We apply the algorithm to synthetic and field seismic data, indicating that the proposed method has a good resolution and stability performance.


Introduction
Seismic inversion can be classified into post-stack impedance and pre-stack multi-parameter inversions. Since the post-stack class only makes use of stacked data and extracts acoustic impedance information, it is relatively useful in obtaining information about the presence of hydrocarbons (Ghosh 2000;Liu et al. 2012). On the other hand, the pre-stack inversion takes advantage of angle gathers that contain amplitude variations with incident angle or offset (AVA/AVO) theory (Aki & Richards 1980), and is able to estimate the parameter group, including P-and S-wave velocities as well as density, constituting one of the most powerful techniques for reservoir prediction and fluid identification (Chiappa & Mazzotti 2009;Zong & Yin 2017;Luo et al. 2019;Ba et al. 2017Ba et al. , 2019Pang et al. 2019;Zhang et al. 2019;Zhao et al. 2014Zhao et al. , 2015Zhao et al. , 2017. Various optimization algorithms have been introduced to solve the pre-stack inversion. Although linearized optimization algorithms, such as the damped least squares (Marquardt 1963) and conjugate gradient (Bae et al. 2012) methods, are robust and highly efficient, they may fall into local minima. Since pre-stack inversion is highly nonlinear, other methods such as genetic algorithms (Mallick 1992), particle swarm optimization (Shaw & Srivastava 2007;Yuan et al. 2009) and the multi-mutation differential evolution algorithm (Gao et al. 2016) have been employed. In particular, simulated annealing (SA), a commonly used nonlinear inversion algorithm, was initially developed by Kirkpatrick et al. (1983). Rothman (1986) first used the SA algorithm to the problem of residual static corrections. Then, Yao (1995) improved the computational efficiency of SA by using fuzzy prior information to determine the minimum temperature and modify the cost function. Endowed with fast cooling speed and high efficiency, fast simulated annealing (FSA) (Szu & Hartley 1987), which is a variant of SA, has been widely used in seismic inversion (Misra & Sacchi 2008;Guo et al. 2018). In FSA, the updating of the model parameters is a key factor of the algorithm. However, pre-stack seismic inversion is a multi-parameter inversion problem, leading to instability of multiple inversion results. The instability can be caused by the update in the conventional FSA, and a subspace method may overcome this issue (Wang & Houseman 1994, 1995Wang 2003;Liu & Wang 2020).
Noises and inadequacy of observed data cause the ill-posed solution to pre-stack inversion (Tarantola 2012). By introducing a regularization technique (Tikhonov & Arsenin 1977), i.e. adding constraints into objective function, the ill-posedness can be mitigated and stability can be achieved (Li & Peng 2017;Li & Zhang 2017). Zhang et al. (2013) combined generalized linear AVO inversion with Bayesian theory, and introduced a three-variable Cauchy distribution to constrain the inversion. Specifically, the approaches that add prior constraints into the objective function based on the Bayesian framework seem to work quite well. Among them, the joint probability distribution is effective (Downton & Ursenbach 2006;Alemie & Sacchi 2011), by which multiple parameters are correlated and stabilized based on statistics, i.e. expectation and covariance. However, these attributes are usually the same for all the layers in standard approaches (Downton & Ursenbach 2006), without considering the discrepancies due to lithology.
We address the two problems mentioned above, by implementing a nonlinear three-parameter pre-stack seismic inversion method based on lithology constraints. Specifically, different prior terms (expectations and covariances) are used for different lithologies in the objective function. Moreover, the multivariate Gaussian probability distribution is incorporated into the model update of the FSA.
The paper is arranged as follows. First, a model is introduced and the objective function is based on lithology constraints. Then, the FSA algorithm, aided by a multivariate probability distribution is proposed, including a discussion of the model update. Finally, the method is tested with synthetic and log data, and applied to field data.

AVA forward modeling
To estimate elastic parameters more accurately, the reflection coefficients of P-waves can be calculated by the exact Zoeppritz equation as (Aki & Richards 1980) , and where V P and V S represent the longitudinal and shear wave velocities, denotes bulk density, 1 denotes the angle of P-wave incidences, 1 is the reflected angle of converted waves and 2 and 2 denote the refracted angles of the P-and S-transmissions, respectively. The subscripts 1 and 2 represent the upper and lower layers, respectively. Snell's law expresses the horizontal slowness p as (2)

The objective function
Under the assumption that subsurface rock consists of a series of lateral layers attributed by discretized elastic (model) parameters, z, the seismic responses can be simulated based on the convolutional model, in which the observed data are obtained by convolving the reflectivity series r with the source wavelet w, i.e. (Alemie & Sacchi 2011;Aleardi & Salusti 2019), where z denote the elastic parameters that represent the properties of subsurface rock; for the three-parameter pre-stack inversion problem, z consists of the discretized V P , V S and ; denotes the incidence angle of the P-wave; r(z, ) denotes the reflection coefficient of the P-wave, which is a function of z and according to equation (1); e denotes the noise resulting from measurement errors and '*' denotes the convolution operator. Compared with post-stack inversion, the elastic-parameter groups, e.g. P-wave velocity, S-wave velocity, density and Poisson's ratio, P-wave impedance, density, etc., can be simultaneously estimated through inversion technique from pre-stack angle gathers. However, inversion by using pre-stack datasets is an ill-posed problem due to instabilities and multimodal solutions. By introducing prior information based on a Bayesian framework (Buland & Omre 2003;Kjønsberg et al. 2010), we can mitigate the ill-posed problem. The posterior probability of the model parameters Prob(m | d obs ) is where Prob(d obs ) denotes the marginal distribution of the observed data and it is a constant value given the known observed gathers d obs ; Prob(m) denotes the priori information on m, which can be assigned as prior information for the unknown model parameters to be inverted and Prob(d obs | m), the likelihood function, corresponds to the forward relationship connecting observations and model parameters. Actually, Prob(d obs | m) represents the similarity between the responses of simulated and observed data during inversion process. Assuming that the noise e follows a zero-mean Gaussian distribution and is independent of m, the likelihood function is expressed as where e 2 denotes the variance of noise, and Sca d and Sca r are scale factors, which reflect the magnitudes of the measured amplitude and the synthesized amplitude, respectively. Given the large amount of observed data, computational efficiency is important for a nonlinear inversion algorithm. Both the forward (modeled) response and observed data represent the relative amplitude of the seismic response, so it is necessary to normalize the modeled response and observed data in each iteration. However, standard normalization approaches usually have low efficiency. Here, the normalization is realized by dividing the (modeled and observed) amplitudes by the corresponding scale factors (i.e. the mean value of single-channel sampling points) in the likelihood function, which is also intended to improve the forward efficiency.
In Bayesian theory, prior information is often used to constrain inversion process, making the obtained posterior probability closer to the real result. Different prior distributions can produce different constraint effects. In the actual inversion, the distribution of model parameters is considered to be Gaussian distribution in accordance with the general law of signal distribution. Since V P , V S and are statistically correlated by log data (Gardner et al. 1974;Todoeschuck et al. 1990;Greenberg & Castagna 1992), we add a prior term to stabilize the multi-parameter inversion. Based on the bivariate Gaussian probability function, the prior term takes the form of (Downton & Ursenbach 2006) where m and C m denote the prior mean value (expectation) and the covariance matrix of m, respectively, and N is the dimension of m. By combining equations (5) and (6), the posterior probability is As for the problem of minimizing the corresponding negative logarithmic posteriori of maximum a posteriori (MAP) estimation of m, the objective function has the following expression:

Lithology-dependent prior constraints
For a three-parameter pre-stack inversion, the covariance matrix C m in the prior constraints is a [3 × 3] matrix (Downton & Ursenbach 2006), which takes the form where the diagonal elements are the variances of V P , V S and , and the nondiagonal elements are the covariances, which describe the statistical correlation of V P and V S , V P and , as well as V S and . Therefore, the prior term, based on multivariate Gaussian distribution is Actually, the statistical characteristic of the elastic parameters is for different lithologies. If the same C m is used for all the layers, the inversion results are erroneous. Therefore, we propose lithology-dependent prior constraints, i.e. different C m values are extracted from logging data and used to formulate the prior term of the specific formation.
The multilayer-model properties are shown in Table 1, from which the first three layers (L1-L3) illustrate the different prior probability distributions for the lithologies, with different mean values and covariances of the elastic parameters. The sampling rate of the model is 1 ms. Figure 1a shows the prior probability of the multivariate Gaussian distribution of V P and V S (equation (6)) of all the layers (without considering the lithological differences) and figure 1 parts b-d show the prior probability of V P and V S for each lithology. Figure 1a indicates that the distributions of all the layers lies between 2500 and 4500 m s −1 and 1200 and 2700 m s −1 , respectively, and the Gaussian ellipse is flat and nearly strip-shaped, indicating the P-wave velocity has a very strong correlation with the S-wave velocity. Compared with figure 1b-d, the statistics are quite different and the different shape of the Gaussian ellipse also indicates that the statistical relationships between V P and V S are quite different. Thus, prior lithological constraints should be used to better represent the statistical properties of the subsurface properties.

Conventional FSA
The goal of the inversion problem is to minimize the objective function or maximize the posterior probability distribution. Therefore, to obtain the estimation of the model parameters, an optimization algorithm is usually employed to achieve the minimization. Derived from the principle of solid annealing that simulates the state of solid from melting to crystallization, the SA algorithm is a commonly used nonlinear optimization algorithm (Sen & Stoffa 1991). As a variant of SA, FAS has the advantage of fast convergence and improved model updating, making it widely applied in practical inversions (Mosegaard & Vestergaard 1991). Different from conventional SA, the acceptance probability of a new solution in FSA obeys a Cauchy-like distribution (Alemie & Sacchi 2011), in which the tail of the distribution is flat and wide, indicating that the probability of accepting large-scale disturbance is large and it is likely to avoid a local minimum. The new solution is (Szu & Hartley 1987): where m (k) denotes the current value and m (k+1) denotes the updated value of m; denotes a random number in the range of 0 and 1; m denotes the search range of m, and sign(.) denotes the sign function. In particular, the generation of new solution is controlled by the temperature T, which is gradually annealing based on the cooling schedule.
In each iteration of the optimization algorithm, the worse solution still has a chance to be accepted when the cost function value increases, so as to prevent the solution from falling into local minima, and a model update so that the cost function decreases is also accepted (Alemie & Sacchi 2011). The acceptance probability P of the updated value is where E denotes the energy variation between the objective function values after the update, and t is a constant value.
The cooling schedule has a great influence on the convergence and efficiency of the conventional fast simulated annealing. The exponential annealing method rapidly cools down in the early iterations, and slowly cools down in the later iterations, finally ensuring global convergence. The cooling equation is expressed by (Geman & Geman 1993;Triki et al. 2005): where T (k) denotes the current temperature that begins with the initial temperature T (0) , and c is the damping factor, which is a constant value.

FSA based on multivariate Gaussian distribution
The three-parameter pre-stack inversion attempts to invert V P , V S and simultaneously. Therefore, the parameter perturbation (update) in FSA involves multiple parameters. In the conventional FSA algorithm, the three parameters are usually perturbed at the same time according to equation (10), which causes instabilities. Liang et al. (2017) proposed to improve the using relationships among the three parameters. Specifically, the P-wave velocity is perturbed first, then the S-wave velocity is evaluated based on a constant P-to S-wave velocity ratio, and density is obtained by the generalized Gardner equation (Gardner et al. 1974) in the parameter perturbation. To achieve a better agreement with the statistical characteristics of the actual formation elastic parameters, a multi-parameter perturbation approach that combines multiple probability density functions is proposed here. We start with the 2D probability density function (multivariate Gaussian distribution function) of V P and V S : The correlation coefficient q is introduced to expand the covariance matrix and equation (13) can be rearranged as where V P and V S denote the mean values of V P and V S , respectively, and q = VpVs / Vp Vs .
V P is updated according to equation (10) as In standard FSA, V P , V S and density are usually perturbed simultaneously, according to equation (10), leading to instability when generating new solutions. To illustrate the proposed model perturbation, we take figure 2 as an example. When V P is specified, a 1D Gaussian (normal) distribution of V S can be obtained from the 2D Gaussian distribution of V P -V S . Substituting the new value of V P into the 2D probability density function P (V P , V S ), and assuming V P -V P = V, equation (13) can be rewritten as Omitting the constant term, the 1D Gaussian distribution of V S (given V P ) can be obtained as Therefore, when V P is specified, equation (16) turns out to be the 1D Gaussian distribution of V S with the expectation (mu s ) and variance (va s ), which can be obtained as 7 Figure 3. The multilayer model listed in Table 1 and the corresponding angle gather.
By combining mu s and va s with equation (10), a new mode disturbance of V S combined with the multivariate probability distribution can be obtained as follows: Similarly, the new value V P is introduced into the 2D probability density function P (V P , ): When V P is specified, equation (19) is also a 1D Gauss distribution of . The expectation of P ( | V P ) is mu and the variance is va . By combining these two values with equation (10), the mode updating of , combined with a multivariate probability distribution, can be obtained as The proposed update method of the three-parameter (referred as to multivariate probability FSA hereafter) is based on a multivariate probability distribution, in which the statistical correlation of multi-parameters is effectively incorporated into the FSA algorithm. The synthetic data test verifies the method in the next section.

Synthetic data
The multilayer model listed in Table 1 is employed again to verify the proposed inversion method. Table 1 shows the mean values of elastic parameters for each lithologic layer, and ±5% of the mean value is taken as the variance of elastic parameters of each layer. The pre-stack angle gather profile generated from this multilayer model is shown in figure 3.
Three sets of inversion test are performed, whose results are shown in figures 4-6. Figure 4 shows the three-parameter results by conventional FSA without any prior terms in the objective function. The update of the V P , V S and follows from equation (10). It can be seen that the accuracy of is poor and the results are unstable, since the discrepancies with the true model are large. Figure 5 shows results by multivariate probability FSA. The prior term is introduced in the objective function with uniform covariance (extracted from all the layers). The update the of S-wave velocity and density follows from equations (18) and (20). It can be seen that the accuracy is better than before, since the discrepancies are smaller. Figure 6 shows by the proposed method, which is based on lithology-dependent prior terms and multivariate probability FSA. The update of the S-wave velocity and density follows from equations (18) and (20). The accuracy and resolution are higher than those shown in figures 4 and 5, especially the density.
To test the robustness of the proposed method, white Gaussian noises are added to the synthesized data. Figure 7 shows the seismic data with noises, of which signal-to-noise (S/N) ratio is 20 dB. Figure 8 shows the three-parameter inversion results by the method. The result accuracy is high and stable, indicating the acceptable stability of the method.    Through the synthetic data tests, the multivariate probability FSA based on lithology-dependent prior constraints shows higher stability and accuracy. Next, real log data is used to gain further confidence in the proposed method.

Log data
In order to verify the algorithm, real well-log data (Well #1) from an area in the northern part of South China Sea is selected. The elastic-parameter model of the target layer from Well #1 is shown in figure 9. This layer is divided into three formations. From top to bottom, there are mudstone, limestone and sandstone formations. Based on a lithological profile, the prior probability of V P -Vs and V P -(den) 2D Gaussian distributions can be obtained. Figure 10 shows the V P -V S and V P -density multivariate Gaussian distribution prior probability for the whole model. The Gaussian ellipse shows that the statistical distribution ranges for the V P , V S and density are 2500-5500, 1000-3000 and 2200-2600 kg m −3 , with the mean values of 3700 m s −1 , 2300 m s −1 and 2450 kg m −3 , respectively. Figures 12-14 show the prior probability of the multivariate Gaussian distribution for different lithologies. It can be seen that the statistical distribution ranges and mean values of the elastic parameters for different lithologies are quite different from those of the whole model, and the shapes of the Gaussian ellipse are different, which indicates that lithology-dependent prior constraints can better represent the statistical properties of the different formations. Figures 11 and 15 show the inversion results of the angle gather along Well #1, which is based on uniform prior constraints and lithology-dependent prior constraints, respectively. In figures 11 and 15, the updated equations of V P , V S and density    are equations (15), (18) and (20), respectively. It can be seen that using uniform prior constraints, the result is unstable and different from the log curves (figure 11). In particular, the actual oil-gas-bearing layer is the middle limestone layer, whose resolution is not high. In contrast, the results based on lithology-dependent prior constraints (figure 15) are more stable and more consistent with the logging profiles, especially for the density. With the intention to further validate the effectiveness of the proposed method, a field data is used in the next section.

Field-data test
The field data was acquired at the northern part of the South China Sea, from which a seismic line was selected for the test (the stacked section is shown in figure 16). There is one well (Well #2) in the study area which is located at the 121st CDP. The line has 212 CDPs, each of which has 11 traces with angle range of 3-36°, interval of 3°and sampling rate of 0.002 s. The target zone is identified and its top and bottom interfaces are depicted by dotted lines in figure 16, involving three groups of formations; i.e. mudstone (top layer), limestone (the target layer) and sandstone (bottom layer).
Three sets of inversion tests are designed here, and we use the initial model with a low frequency trend (established according to the background velocities). Correlation coefficient between the pre-stack seismic angle gathers synthesized with the initial model and the actual seismic data is 0.4. In particular, the model update in the proposed multivariate probability FSA is compared with the conventional approach in Liang et al. (2017). Figure 17 shows the results without any prior constraints, where the update of the P-wave velocity is based on equation (10). According to well-log data, the correlation of the three   parameters is analyzed, and relationships between the velocities and/or density are derived. The S-wave velocity and density update is in agreement with Liang et al. (2017), where V p k+1 is the updated value. Figures 18 and 19 show the results by the multivariate probability FSA with uniform prior constraints and lithologydependent prior constraints, respectively. The lithology-dependent constraints were estimated from the data of Well #2 and the model update of the three parameters is based on equations (15), (18)   bottom of the middle limestone are discontinuous. Moreover, the density does not follow the formation structure very well, and shows a poor agreement with the log curve. The results obtained with the multivariate probability FSA (figure 18) with uniform prior constraints are better, but those of the lithology-dependent prior constraints ( figure 19) are the best. However, in certain intervals, such as between 1775 and 1790 ms in figure 22, the results show some anomalies, which may come from a cumulative error caused by the stochastic algorithm. 16 Figure 21. Comparison of the results in CDP121 shown in figure 18 with the well-log profiles of V P , V S and in Well #2.  figure 19 with the well-log profiles of V P , V S and in Well #2.

Conclusions
We propose a three-parameter pre-stack seismic inversion method using lithology-dependent prior constraints of the objective function based on a Bayesian approach. It is shown that these constraints can better represent the statistical properties of the different lithologies, and the results have a high accuracy and stability. The method uses a multivariate Gaussian distribution of model parameters, and the problem is solved with fast simulated annealing. Compared with the conventional approach, the new method effectively improves the stability of the multi-parameter inversion. Tests with synthetic, logging data and 2D field data verify the resolution of the algorithm.