Bayesian Markov chain Monte Carlo inversion for anisotropy of PP- and PS-wave in weakly anisotropic and heterogeneous media

A single set of vertically aligned cracks embedded in a purely isotropic background may be considered as a long-wavelength effective transversely isotropy (HTI) medium with a horizontal symmetry axis. The crack-induced HTI anisotropy can be characterized by the weakly anisotropic parameters introduced by Thomsen. The seismic scattering theory can be utilized for the inversion for the anisotropic parameters in weakly anisotropic and heterogeneous HTI media. Based on the seismic scattering theory, we first derived the linearized PP- and PS-wave reflection coefficients in terms of P- and S-wave impedances, density as well as three anisotropic parameters in HTI media. Then, we proposed a novel Bayesian Markov chain Monte Carlo inversion method of PP- and PS-wave for six elastic and anisotropic parameters directly. Tests on synthetic azimuthal seismic data contaminated by random errors demonstrated that this method appears more accurate, anti-noise and stable owing to the usage of the constrained PS-wave compared with the standards inversion scheme taking only the PP-wave into account.


Introduction
Aligned vertical cracks, which are embedded in a purely isotropic background and treated as a vital storage space and migration pathway of reservoirs, may be considered as a long-wavelength effective transversely isotropy (HTI) medium with a horizontal symmetry axis (Mallick et al. 1998;Thomsen 2002;Tsvankin and Grechka 2011;Liu and Martinez 2012;Crampin and Gao 2014). The crack-induced HTI anisotropy can be characterized by the weakly anisotropic parameters introduced by Thomsen (Thomson 1986;Tsvankin 1996;Rüger 1997Rüger , 1998. Linearized PP-and PS-wave reflection coefficients in weakly anisotropic and heterogeneous media are difficult to be derived due to the complexity of seismic wave propagation in such media. Based on the assumption of weak elastic anisotropy and anisotropic notations (Thomson 1986;Tsvankin 1996), linearized PP-wave reflection coefficients were derived to characterize the weakly anisotropic HTI media (Rüger 1997(Rüger , 1998Pšenčik and Vavryčuk 1998;Pšenčik and Martins 2001). However, the inversion results obtained just by taking the PP reflectivity into account were less accurate than the ones obtained by taking the PS reflectivity into account. Thus, the conventional PP-wave pre-stack inversion methods are poorly constrained, resulting in high uncertainties to estimate the crack parameters in HTI media. As a consequence, individual parameters of the cracked reservoir cannot be recovered uniquely, or the recovered quantities may not be sufficiently accurate. Taking the fact that S-waves do not propagate through fluid into account, the study of combining PP-and PS-wave seismic data may increase our ability to derive lithology and fluid information from seismic data, reducing the uncertainty in reservoir characterization (Veire and Landro 2006). Using the firstorder perturbation theory, Jilek (2002) derived the general closed-form approximations of azimuthally dependent weak-contrast, weak-anisotropy converted PS-wave reflection coefficients at a horizontal weak-contrast interface separating two weakly anisotropic half-spaces. Tang et al. (2007) studied the PP-and PS-wave seismic reflection response in cracked tight gas reservoir. Garotta and Granger (2010) proposed PP-and PS-wave reflection approximations and applied the multi-component seismic data to complete the reservoir characterization. Based on the closed-form PS-wave approximations of weakly HTI media (Jilek 2002), Cui et al. (2010) derived a more simple form of PS-wave reflection coefficient approximation. These derived equations mentioned above, however, may be restricted to the horizontally layered medium assumption, because the seismic waves may scatter when coming across the cracks in the heterogeneous media during the propagation (Shaw and Sen 2006;Burns et al. 2007;Yin et al. 2013). Meanwhile, first-order PP-and PS-wave reflections usually lose their accuracy on account of the neglected high-order terms (Shaw and Sen 2004). In this paper, linearized PP-and PS-wave reflection coefficients in weakly anisotropic and heterogeneous HTI media are derived using seismic scattering theory.
Seismic rock-physics models are the basis of reservoir quantitative prediction, which build the bridge between seismic response and fracture properties such as fracture density, pore aspect ratio and pore-filled fluids; thus, fractured anisotropic rock-physics models contribute to the estimate of fracture properties (Avseth et al. 2005). There are three commonly used cracked rock-physics models, including the linear slip model (Schoenberg 1980(Schoenberg , 1983, the aligned penny-shaped crack model (Hudson 1980(Hudson , 1981 and the fractured porous medium model (Thomsen 1995;Hudson et al. 1996). Studies show that reasonable anisotropic rock-physics models help to evaluate the elastic and anisotropic parameters (Zhang et al. 2013;Chen et al. 2015).
Pre-stack seismic inversion can be used for the estimate of cracked parameters. Mallick et al. (1998) was the first scholar to study the pre-stack seismic inversion based on the azimuthal seismic data. Bachrach et al. (2009) reconstructed the cracked reservoir anisotropic elastic parameters and crack characteristics and inverted the anisotropic parameters based on the PP-wave properties of amplitude variation with the angle and azimuth (AVAZ). Downton and Roure (2010) studied the inversion of the crack weaknesses based on the linear slip model. Chen et al. (2014) estimated the crack weaknesses based on the simplified small-angle approximate AVAZ inversion and the rock-physics model. In this paper, we proposed a novel Bayesian Markov chain Monte Carlo (MCMC) inversion method (Zhang et al. 2011) of PP-and PS-wave for anisotropy in weakly anisotropic and heterogeneous HTI media to characterize the crack properties. Tests on synthetic azimuthal seismic data contaminated by random errors demonstrated that this method appears more accurate, anti-noise and stable owing to the usage of the constrained PS-wave compared with the standards inversion scheme taking only the PP-wave into account.

Theory
2.1 Effective stiffness tensor in weakly anisotropic and heterogeneous HTI media Perturbation theory mainly used for wave propagation problems in weakly anisotropic elastic media (Thomson 1986;Vavryčuk 1997Vavryčuk , 2003Pšencik and Vavrycuk 1998) can be also generalized to the parameterization in heterogeneous media, which is introduced as a small perturbation of a purely isotropic elastic background. Thus, the effective elastic stiffness tensor of heterogeneous media can be then considered to be the sum of the stiffness tensor of the homogeneous isotropic background and the perturbations (Pšencik and Vavrycuk 1998; Shaw and Sen 2004;Yin et al. 2013;Zong et al. 2015): where C, C 0 and DC are the effective elastic stiffness tensor of the heterogeneous media, the homogeneous isotropic background and the weak perturbation, respectively. Similarly, the effective elastic stiffness tensor of weakly anisotropic and heterogeneous media could be thus written as: where DC iso and DC ani represent the stiffness tensor of the isotropic perturbation and the anisotropic perturbation, respectively. Before proceeding to obtain the effective stiffness tensor of weakly anisotropic and heterogeneous HTI media, the Pand S-wave velocities of the background media are first addressed. According to Pšenčik and Gajewski (1998), the best choice for the background P-and S-wave velocities of reflection seismic is the vertical P-wave velocity and the vertical S-wave velocity with polarization in the plane of isotropy, respectively. Moreover, we found that such a preference for the background medium in our approach only has an effect on the required algebra and the structure of the isotropic and anisotropic part of the reflection coefficients. In this paper, we choose q 0 a 2 0 ¼ c 0 33 and q 0 b 2 0 ¼ c 0 55 , where a 0 , b 0 and q 0 represent the background P-and S-wave velocity, and density.
In the case of isotropic background and isotropic perturbation, the effective elastic stiffness tensor is given by: where C 33 ¼ C 0 33 þ DC 33 and C 55 ¼ C 0 55 þ DC 55 . The two parts in Eq. (3) represent isotropic background and isotropic perturbations, respectively.
According to the perturbation theory, the effective elastic stiffness tensor of weakly anisotropic and heterogeneous HTI media can be also parameterized as a perturbation over a homogeneous isotropic background given as:  þ C 11 À C 33 C 13 À C 33 þ 2C 55 C 13 À C 33 þ 2C 55 0 0 0 The three parts in Eq. (4) represent isotropic background, weakly isotropic perturbations and weakly anisotropic perturbations introduced by the aligned vertical cracks, respectively. Following Pšenčik and Vavryčuk (1998), three weakanisotropy (WA) parameters (Thomson 1986) are introduced as: and Tsvankin (1996) and Rüger (1997Rüger ( , 1998 used the WA parameters introduced by Thomson (1986), which can be written as: and According to Sayers (1994), the differences between the above WA parameters are of second orders. Thomsen's WA parameters can be expressed in terms of the anisotropic parameters in Eq. 5a-c in the following ways: and We rewrite the perturbations of stiffness tensor in Eq. (4) in terms of introduced anisotropic parameters e, d and c for the case of a perturbation over a homogeneous isotropic background: 2.2 Linearized PP-and PS-wave reflection coefficients based on seismic scattering theory With the Born integral and stationary phase method, the linearized PP-and PS-wave reflection coefficients for arbitrary anisotropic media embedded in an isotropic background media are given by (Shaw and Sen 2004) and where h and u are the angles of incidence and reflection, respectively; r 0 is the point on the horizontal interface; S r 0 ð Þ is scattering function (Eaton and Stewart 1994;Burridge et al. 1998;Č ervený 2001), and its expression is given by: where and p and t represent the slowness vectors and the polarization vectors, respectively. The prime represents the scattered wave.
The linearized PP-and PS-wave reflection coefficients in terms of crack weaknesses for the case of an interface separating two HTI medium are then derived as (see Appendix): and where I P0 = q 0 a 0 and I S0 = q 0 b 0 represent the background P-and S-wave impedances; DI P and DI S indicate the changes in P-and S-wave impedances across the interface, respectively; h and u are the angles of incidence and reflection, respectively; / = / obs -/ sym is the azimuthal difference angle between the observation azimuth angle / obs and the crack symmetry axis azimuth angle / sym , and in order to obtain the linear expression of Eqs. 16a and 16b, the crack orientation / sym can be first extracted using the amplitude inversion method of small incident angle data easily and accurately (Jenner 2002;Mahmoudian and Margrave 2012). But in this paper we did not discuss the azimuth angle.

PP-and PS-wave nonlinear scattering inversion for anisotropic parameters based on MCMC method
The inverse problem in our study can be expressed as: where d PPÀPS is the PP-and PS-wave azimuthal seismic data; e represents noise, and it is simulated to be the additional part and independent of the model parameters; f(Á) represents forward operators, that is, PP-and PS-wave reflection coefficients. Based on the Bayesian framework, the posterior probability density function (pdf) of six elastic and anisotropic parameters can be calculated by the a priori pdf and the likelihood function: where m ¼ I P ; I S ; q; d; e; c ½ T represent the six elastic and anisotropic parameters inverted.
The first term at the right side of Eq. 14 is the a priori pdf, which is obtained from other sources of information such as the well log data or core data to assure that the inverted parameters can contain the low-frequency information. Under the assumption of I P , I S , q, d, e and c are independent of each other and the a priori pdf is consistent with Gaussian distribution, whose means and standard a 2 0 sin 2 / 2 cos 2h þ sin 2h cot u À sin 2h tan u ð Þ Dc deviations are l I P ; l I S ; l q ; l d ; l e ; l c and r I P ; r I S ; r q ; r d ; r e ; r c , respectively. So the a priori pdf of six elastic and anisotropic parameters is given by: where N is the number of the inverted parameters; r 2 I P ; r 2 I S ; r 2 q ; r 2 d ; r 2 e ; r 2 c are the corresponding variances of six elastic and anisotropic parameters, which can be gotten by the statistical analysis of well information.
The second term of Eq. 14 is referred to the likelihood function, which describes the relationship between the observed seismic data and six parameters. In the case of observable PP-and PS-wave seismic data containing a Gaussian distribution random noise e with zero mean and r PP , r PS standard deviation, thus the likelihood function p d PPÀPS jm ð Þis given by: To obtain the posterior pdf p mjd PPÀPS ð Þ , we use the MCMC algorithm to generate the Markov chains (Zhang et al. 2011, which avoid the calculation of the complicated marginal integral in the Bayesian rules owing to the algorithm structure of MCMC method.
The most common algorithms to construct transmission kernel are Metropolis-Hastings algorithm and Gibbs sampler. We choose Metropolis-Hastings algorithm in our study. The way to generate the Markov Chain is shown in Fig. 1 Þand a uniform random number U between 0 and 1 to decide whether to accept the transmission. If U a m i ; m Ã ð Þ, we accept the candidate value m Ã , otherwise remain the precious value m i not to change. The acceptance probability that can be proved that the transmission kernel satisfies the detailed balance condition is as follows: Then, the Markov chain has the stationary distribution p m ð Þ. Metropolis algorithm whose proposal distribution is a symmetrical random walk is a special case of Metropolis-Hastings algorithm, that is where d represents the step size. q m Ã ; m i ð Þ¼q m i ; m Ã ð Þ; thus, the acceptance probability a m i ; m Ã ð Þ in our study is given by: Combing Eqs. 15 and 16, we can obtain the posterior pdf given as:  ð21Þ and w m ð Þ ¼ À In order to get the Markov chain converging to the posterior pdf p mjd PPÀPS ð Þ , the stationary distribution p m ð Þ must be the posterior pdf. Thus, the stationary distribution p m ð Þ is the product of a prior probability density and likelihood function.
So the final acceptance probability can be given as: Finally, we can get the Markov chains which converge to the posterior pdf, and we can gain the final inverted results of six elastic and anisotropic parameters by the statistical analysis of Markov chains. But what I want to emphasize is that the direct inversion of six elastic and anisotropic parameters can eliminate the errors, resulting from the trace integration or the recurrence of the reflection coefficients of elastic and anisotropic parameters. So the cumulative error may be too small to ignore.

Well log data from cracked carbonate reservoir
There is one carbonate reservoir in Western China, and the crack orientation is nearly NE-SW strike. Studies indicate that the reservoir shows AVAZ anisotropic characteristics. In this paper, the method of Bayesian Markov chain Monte Carlo inversion for anisotropy of PP-and PS-wave in weakly anisotropic and heterogeneous media is first tested on synthetic data of real well A from the carbonate reservoir.
According to the process of building anisotropic cracked rock-physics models (Zhang et al. 2013), the six elastic and anisotropic parameters of well A are calculated. Figure 2 shows the flow of anisotropic rock-physics models. Figure 3 shows the estimated elastic and anisotropic parameters of well A based on the built rock-physics models.
Step 1 Step 2 Step 3 Step  Step 1 Estimate the elastic modulus of mineral mixture; Step 2 Append intergranular pores into the rock matrix and calculate the elastic modulus of noncracked dry rock skeleton; Step 3 Append cracks into the non-cracked dry rock skeleton based on the linear slip model (Schoenberg 1980(Schoenberg , 1983 and then calculate the elastic modulus of cracked dry rock; Step 4 Implement the anisotropic fluid substitution in the cracked saturated rock (Brown and Korringa 1975) and finally calculate the elastic and anisotropic parameters from the stiffness matrix of the cracked saturated rock The synthetic seismic data were generated using the convolution model based on the zero-phase Ricker wavelet, the derived HTI reflection coefficient Eq. 12a, b and the estimated elastic and anisotropic parameters shown in Fig. 3 driven by the built anisotropic rock-physics model. Synthetic seismic data with signal-to-noise ratio (SNR) being 8 are shown in Fig. 4, which is convolved with 45 Hz Ricker wavelet. Only one angle trace of 0°azimuth is displayed, and other four azimuthal angle traces are 45°, 90°, 135°and 180°. The comparison between the estimated results and the well log data is shown in Fig. 5, where the green lines indicate the initial model and blue dashed ellipses show the cracked reservoir based on the known logging interpretation. It shows that the six elastic and anisotropic parameters estimated using the method of PP-and PS-wave nonlinear scattering inversion based on MCMC can meet inversion requirements preferably, so it verifies the feasibility of the method proposed in this paper.  Now, random noise is added into the synthetic seismic data, which is displayed in Fig. 6, and only one angle trace of 0°azimuth is displayed and other four azimuthal angle traces are also 45°, 90°, 135°and 180°. Figure 7 shows that the inverted elastic and anisotropic parameters with SNR being 2 are also consistent with the well log data, where the cracked reservoir can be identified obviously by using the inverted parameters, so it further validates the availability of cracked reservoir identification and detection with the method of PP-and PS-wave nonlinear scattering inversion for six elastic and anisotropic parameters based on MCMC.

Complex 2D overthrust model
We select a complex 2D overthrust model constructed by Aminzadeh et al. (1994) and Mulder et al. (2006) to further test our method, which can be considered as a HTI medium. We use the convolution formula to produce the PP-and PS-wave synthetic seismic gathers with 0°azimuth, 45°azimuth, 90°azimuth, 135°azimuth and 180°a zimuth, and random noise is then added to the PP-and PS-wave azimuthal seismic gathers. The original models are shown in Fig. 8. Figures 9 and 10 show the inverted elastic and anisotropic parameters of 2D overthrust model and seismic trace (CDP100) with SNR being 8, and Figs. 11 and 12 show the same results but with SNR being 2. The noisy inverted results show that the 2D overthrust model values estimated using PP-and PS-wave nonlinear scattering inversion based on MCMC are consistent with the true values roughly, and the estimated values could more accurately describe the stratigraphic features of 2D over-thrust model, where its lateral continuity is also well, so it verifies the reliability and accuracy of the method proposed in this paper. The results of the single seismic trace are also consistent with the true values, so it further validates the feasibility of the method (Figs. 11, 12).

Conclusions
The crack-induced anisotropy in weakly anisotropic and heterogeneous HTI media can be characterized by the seismic scattering theory. The linearized PP-and PS-wave reflection coefficients in terms of P-and S-wave impedances, density as well as three anisotropic parameters in weakly anisotropic and heterogeneous HTI media were first derived based on the seismic scattering theory, and then, we proposed a method of Bayesian MCMC inversion for anisotropy of PPand PS-wave in weakly anisotropic and heterogeneous media. Tests on synthetic azimuthal seismic data show that the method proposed in this paper may provide an alternative way to identify the crack-induced anisotropy. Next, we may take the symmetry axis azimuth into account together. and the Fundamental Research Funds for the Central Universities (15CX08002A) for their funding in this research. We also thank the anonymous reviewers for their constructive suggestions.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.