Application of Reversible Jump Markov Chain Monte Carlo Algorithms to Elastic and Petrophysical Amplitude-Versus-Angle Inversions

Abstract

We infer the elastic and petrophysical properties from pre-stack seismic data through a transdimensional Bayesian inversion. In this approach the number of model parameters (i.e. the number of layers) is treated as an unknown, and a reversible jump Markov Chain Monte Carlo (rjMCMC) algorithm is used to sample the variable-dimension model space. This inversion scheme provides a parsimonious solution, and reliably quantifies the uncertainties affecting the estimated model parameters. Parallel tempering, which employs a sequence of interacting Markov chains in which the likelihood function is successively relaxed, is used to improve the efficiency of the probabilistic sampling. In addition, the delayed rejection updating scheme is employed to speed up the convergence of the rjMCMC algorithm to the stationary regime. Both elastic and petrophysical inversions invert the amplitude versus angle responses and employ a convolutional forward modelling based on the exact Zoeppritz equations. First, synthetic tests are used to assess the reliability of the implemented rjMCMC algorithms, then their applicability is demonstrated by inverting field seismic data acquired onshore. In this case the inversion was aimed at inferring the elastic and petrophysical properties around a gas-saturated reservoir hosted in a shale-sand sequence. In this case, the final outcomes provided by the rjMCMC algorithms are also compared with the predictions of linear Bayesian elastic and petrophysical inversions. The synthetic and field data examples demonstrate that the implemented algorithms can successfully estimate model uncertainty, model dimensionality and subsurface parameters with an affordable computational cost.

Appendices

Appendix: Acceptance Probabilities for Elastic AVA Inversion

Here we shortly derive the acceptance probabilities for each move in the implemented rjMCMC algorithm. Since we adopt the same MCMC recipe for both the elastic and petrophysical inversions, in the following we limit the attention to the elastic inversion. Similar results would be derived for the petrophysical inversion: the only difference is that the elastic property vector ${\mathbf{e}}$ would be replaced by the petrophysical property vector ${\mathbf{r}}$. For a more detailed mathematical description we refer the interested reader to Bodin and Sambridge (2009) and Dosso et al. (2014).

The Acceptance Probabilities for Fixed Dimensional Moves

In the elastic property move or interface move the number of layers remain fixed. In the former case the elastic properties are determined by applying the principal-component perturbation described in the text. For the latter we consider a Gaussian proposal with a given standard deviation $\sigma_{z}$ centred around the current interface location. Then, the proposed interface location vector will be derived as:

$${\mathbf{z}}^{\prime } = {\mathbf{z}} + g\sigma_{z} {\mathbf{u}}_{{\mathbf{i}}} ,$$

(24)

where $g$ is a random number drawn from the normal distribution $N\left( {0,1} \right)$ with zero mean and standard deviation 1, and ${\mathbf{u}}_{{\mathbf{i}}}$ is the unit vector along the i-th model space dimension that identifies the perturbed interface. For these moves, the prior ratio and the proposal ratio are equal to 1, and the acceptance probability simply becomes:

$$\begin{aligned} p\left( {{\mathbf{m}}^{\prime } ,n^{\prime } {|}{\mathbf{m}},n} \right) & = {\min}\user2{ }\left[ {1,\frac{{p\left( {{\mathbf{d}}{|}{\mathbf{m}}^{\prime } ,n^{\prime } } \right)}}{{p\left( {{\mathbf{d}}{|}{\mathbf{m}},n} \right)}}} \right] \\ & = {\text{min }}\left[ {1,{\text{ exp}}\left( { - \frac{{{{\varphi }}\left( {{\mathbf{m}}^{\prime } ,n^{\prime } } \right) - {{\varphi }}\left( {{\mathbf{m}},n} \right)}}{2T}} \right)} \right]. \\ \end{aligned}$$

(25)

The Acceptance Probabilities for Transdimensional Moves

3.1 Prior Ratio

The prior ratio can be written as follows:

$$\frac{{p\left( {{\mathbf{m}}^{\prime } ,n^{\prime } } \right)}}{{p\left( {{\mathbf{m}},n} \right)}} = \frac{{p({\mathbf{e}}^{\prime } |n^{\prime } )}}{{p({\mathbf{e}}|n)}}\frac{{p({\mathbf{z}}^{\prime } |n^{\prime } )}}{{p({\mathbf{z}}|n)}}\frac{{p\left( {n^{\prime } } \right)}}{p\left( n \right)},$$

(26)

where $n^{\prime} = n + 1$ for a birth move and $n^{\prime} = n - 1$ for a death move. By extending Eq. 14.1 for a birth move we get:

$$\begin{aligned} \frac{{p\left( {{\mathbf{m}}^{\prime } ,n^{\prime } } \right)}}{{p\left( {{\mathbf{m}},n} \right)}} & = \frac{{\left( {n + 1} \right)!\left( {N - \left( {n + 1} \right)} \right)!}}{{N! \Delta n \mathop \prod \nolimits_{i = 1}^{3} \left( {\Delta {\mathbf{e}}_{i}^{n + 1} } \right)}} \times \frac{{N! \Delta n \mathop \prod \nolimits_{i = 1}^{3} \left( {\Delta {\mathbf{e}}_{i}^{n} } \right)}}{{n!\left( {N - n} \right)!}} \\ & = \frac{n + 1}{{\left( {N - n} \right)\mathop \prod \nolimits_{i = 1}^{3} \left( {\Delta {\mathbf{e}}_{i} } \right)}} \\ \end{aligned}$$

(27)

where N is the number of discrete time samples forming the considered vertical interval.

Similarly, for a death move we obtain:

$$\begin{aligned} \frac{{p\left( {{\mathbf{m}}^{\prime } ,n^{\prime } } \right)}}{{p\left( {{\mathbf{m}},n} \right)}} & = \frac{{\left( {n - 1} \right)!\left( {N - \left( {n - 1} \right)} \right)!}}{{N! \Delta n \mathop \prod \nolimits_{i = 1}^{3} \left( {\Delta {\mathbf{e}}_{i}^{n - 1} } \right)}} \times \frac{{N! \Delta n \mathop \prod \nolimits_{i = 1}^{3} \left( {\Delta {\mathbf{e}}_{i}^{n} } \right)}}{{n!\left( {N - n} \right)!}} \\ & = \frac{{\mathop \prod \nolimits_{i = 1}^{3} \left( {\Delta {\mathbf{e}}_{i} } \right) \left( {N - n + 1} \right) }}{n} \\ \end{aligned}$$

(28)

3.2 Proposal Ratio

The proposal ratio can be decomposed as:

$$\frac{{q({\mathbf{m}},n|{\mathbf{m}}^{\prime } ,n^{\prime } )}}{{q({\mathbf{m}}^{\prime } ,n^{\prime } |{\mathbf{m}},n)}} = \frac{{q({\mathbf{z}}|{\mathbf{z}}^{\prime } ,n^{\prime } )}}{{q({\mathbf{z}}^{\prime } |{\mathbf{z}},n)}}\frac{{q({\mathbf{e}}|{\mathbf{e}}^{\prime } ,n^{\prime } )}}{{q({\mathbf{e}}^{\prime } |{\mathbf{e}},n)}} .$$

(29)

For a birth move the probability of generating new property values at the (n + 1)-th layer is:

$$q\left( {{\mathbf{e^{\prime}}}{|}{\mathbf{e}},n} \right) = \frac{1}{{\mathop \prod \nolimits_{i = 1}^{3} \sqrt {2\pi } \sigma_{{{\mathbf{e}}_{i} }}^{{}} }}{\exp}\left( { - \mathop \sum \limits_{i = 1}^{3} \frac{{\left( {{\mathbf{e}}_{i}^{\prime } - {\mathbf{e}}_{i} } \right)^{2} }}{{2\sigma_{{{\mathbf{e}}_{i} }}^{2} }}} \right).$$

(30)

where $\sigma_{{{\mathbf{e}}_{i} }}^{{}}$ is the standard deviation for the proposal distribution for the i-th elastic property, and ${\mathbf{e}}_{i}$ is the i-th elastic property in the current model within the time interval associated to the new (n + 1)-th layer. The probability of the reverse step (removing the elastic properties when a layer is deleted) is:

$$q({\mathbf{e}}|{\mathbf{e^{\prime}}},n^{\prime } ) = 1.$$

(31)

The probability to generate the (n + 1)-th layer results:

$$p\left( {{\mathbf{z^{\prime}}}{|}{\mathbf{z}},n} \right) = \frac{1}{N - n}.$$

(32)

The probability of the reverse step in a birth move (deleting the layer positions) is:

$$p\left( {{\mathbf{z^{\prime}}}{|}n^{\prime}} \right) = \frac{1}{n + 1}.$$

(33)

By combining equations Eqs. 29–32 for a birth step we found:

$$\frac{{q({\mathbf{m}},n|{\mathbf{m^{\prime}}},n^{\prime})}}{{q({\mathbf{m^{\prime}}},n^{\prime}|{\mathbf{m}},n)}} = \frac{{\left( {N - n} \right)\mathop \prod \nolimits_{i = 1}^{3} \sqrt {2\pi } \sigma_{{{\mathbf{e}}_{i} }}^{{}} }}{n + 1}{\exp}\left( {\mathop \sum \limits_{i = 1}^{3} \frac{{\left( {{\mathbf{e}}_{i}^{\prime } - {\mathbf{e}}_{i} } \right)^{2} }}{{2\sigma_{{{\mathbf{e}}_{i} }}^{2} }}} \right)$$

(34)

With similar mathematical derivation the proposal for the death move is:

$$\frac{{q({\mathbf{m}},n|{\mathbf{m^{\prime}}},n^{\prime})}}{{q({\mathbf{m^{\prime}}},n^{\prime}|{\mathbf{m}},n)}} = \frac{n}{{\left( {N - n + 1} \right)\mathop \prod \nolimits_{i = 1}^{3} \sqrt {2\pi } \sigma_{{{\mathbf{e}}_{i} }}^{{}} }}{\exp}\left( { - \mathop \sum \limits_{i = 1}^{3} \frac{{\left( {{\mathbf{e}}_{i}^{\prime } - {\mathbf{e}}_{i} } \right)^{2} }}{{2\sigma_{{{\mathbf{e}}_{i} }}^{2} }}} \right).$$

(35)

3.3 Acceptance Probability

By combining the prior ratio (Eq. 27), the proposal ratio (Eq. 35) and with the likelihood ratio (Eq. 25) for a birth move we get:

$$\left[ {p\left( {{\mathbf{m^{\prime}}},n^{\prime}{|}{\mathbf{m}},n} \right)} \right]_{birth} = \min \left[ {1,\frac{{\mathop \prod \nolimits_{i = 1}^{3} \sqrt {2\pi } \sigma_{{{\mathbf{e}}_{i} }}^{{}} }}{{\mathop \prod \nolimits_{i = 1}^{3} \Delta {\mathbf{e}}_{i} }}{\exp}\left( {\left( {\mathop \sum \limits_{i = 1}^{3} \frac{{\left( {{\mathbf{e}}_{i}^{\prime } - {\mathbf{e}}_{i} } \right)^{2} }}{{2\sigma_{{{\mathbf{e}}_{i} }}^{2} }}} \right) - \frac{{{{\varphi }}\left( {{\mathbf{m^{\prime}}},n} \right) - {{\varphi }}\left( {{\mathbf{m}},n} \right)}}{{2\sigma_{d}^{2} T}}} \right)} \right].$$

(36)

For a death move we obtain:

$$\left[ {p\left( {{\mathbf{m^{\prime}}},n^{\prime}{|}{\mathbf{m}},n} \right)} \right]_{death} = {\min}\left[ {1,\frac{{\mathop \prod \nolimits_{i = 1}^{3} \Delta {\mathbf{e}}_{i} }}{{\mathop \prod \nolimits_{i = 1}^{3} \sqrt {2\pi } \sigma_{{e_{i} }}^{{}} }}{\exp}\left( {\left( {\mathop \sum \limits_{i = 1}^{3} - \frac{{\left( {{\mathbf{e}}_{i}^{\prime } - {\mathbf{e}}_{i} } \right)^{2} }}{{2\sigma_{{{\mathbf{e}}_{i} }}^{2} }}} \right) - \frac{{{{\varphi }}\left( {{\mathbf{m^{\prime}}},n} \right) - {{\varphi }}\left( {{\mathbf{m}},n} \right)}}{{2\sigma_{d}^{2} T}}} \right)} \right].$$

(37)

After determining the acceptance probability, we draw a random number $\alpha$ from the uniform distribution $U\left( {0, 1} \right)$. If $p\left( {{\mathbf{m^{\prime}}},n^{\prime}{|}{\mathbf{m}},n} \right) \ge \alpha$ then $\left( {{\mathbf{m}},n} \right) = ({\mathbf{m}}^{^{\prime}} ,n)\user2{ }$ otherwise the current model is repeated in the chain and another perturbation is proposed.