Seismic AVO statistical inversion incorporating poroelasticity

Seismic amplitude variation with offset (AVO) inversion is an important approach for quantitative prediction of rock elasticity, lithology and fluid properties. With Biot–Gassmann’s poroelasticity, an improved statistical AVO inversion approach is proposed. To distinguish the influence of rock porosity and pore fluid modulus on AVO reflection coefficients, the AVO equation of reflection coefficients parameterized by porosity, rock-matrix moduli, density and fluid modulus is initially derived from Gassmann equation and critical porosity model. From the analysis of the influences of model parameters on the proposed AVO equation, rock porosity has the greatest influences, followed by rock-matrix moduli and density, and fluid modulus has the least influences among these model parameters. Furthermore, a statistical AVO stepwise inversion method is implemented to the simultaneous estimation of rock porosity, rock-matrix modulus, density and fluid modulus. Besides, the Laplace probability model and differential evolution, Markov chain Monte Carlo algorithm is utilized for the stochastic simulation within Bayesian framework. Models and field data examples demonstrate that the simultaneous optimizations of multiple Markov chains can achieve the efficient simulation of the posterior probability density distribution of model parameters, which is helpful for the uncertainty analysis of the inversion and sets a theoretical fundament for reservoir characterization and fluid discrimination.


Introduction
Subsurface rock is composed of solid mineral matrix, dry pores/fractures and various fluid mixtures and set foundation for geophysical interpretation (Biot 1956;Han and Batzle 2004;Russell et al. 2011;Yin et al. 2015;Ding et al. 2019a;Li et al. 2020). Biot-Gassmann's theory describes the propagation of elastic wave in saturated porous media and establishes the quantitative relationship between rock modulus and seismic wave velocity, which is the physical basis of seismic hydrocarbon discrimination (Gassmann 1951;Biot 1956;Ding et al. 2019b;Liu et al. 2020;Li et al. 2020). Prestack seismic AVO (amplitude variation with offset) inversion is an important approach for quantitative characterization of rock elasticity, physical properties, lithology and fluid properties (Smith and Gidlow 1987;Rutherford and Williams 1989;Goodway et al. 1997;Gray et al. 1999;Russell et al. 2011;Zhang et al. 2012;Yin et al. 2015;Zong et al. 2015, Zong andGrana 2016). The sensitivity of model parameters to fluid types is a key factor affecting the accuracy of fluid discrimination in hydrocarbon-bearing reservoirs. With the guidance of rock physics, scholars define different characteristic parameters related to pore fluid properties as fluid indicators and establish various quantitative relationships between pore fluid types and elastic parameters (Ostrander 1984;Rutherford and Williams 1989;Goodway et al. 1997;Castagna et al. 1998;Han and Batzle 2004;Yin and Zhang 2014;Yin et al. 2015). Russell et al. (2003Russell et al. ( , 2011 derived the Gassmann fluid term f that mainly reflected the fluid type and proposed a linear AVO approximate equation characterized by Gassmann Edited by Jie Hao fluid term. Yin and Zhang (2014) proposed the fluid bulk modulus K f as the fluid factor based on Gassmann equation and critical porosity model (Nur et al. 1998) and developed the solid-fluid decoupling prestack seismic fluid identification technology. Lang and Grana (2018) derived the partial derivatives of rock bulk modulus, shear modulus and density with rock-matrix, fluid modulus and porosity, respectively, and developed a Bayesian AVO linear inversion based on the explicit solution of posterior probability density function (PDF) and sequential simulation.
For seismic AVO inversion algorithms in different domains (Buland 2003;Yuan et al. 2015;Yin et al. 2013Yin et al. , 2016Li et al. 2017a, b;Zong et al. 2017), the timecomplex frequency-domain AVO inversion is helpful to consider the attenuation effects during the propagation of seismic waves and improve inversion resolution. Li et al. (2017b) deduced the seismic AVO forward operator in the complex frequency domain and estimated the low-frequency information of the model parameters with AVO stepwise inversion in complex frequency domain. Zong et al. (2017) studied the influence of attenuation coefficients and frequency components in complex frequency-domain AVO inversion on mining the low-frequency information of model parameters. Since the high-frequency information contained in well-logging data is utilized in statistical inversion incorporating statistical learning theory, the resolution and identification of thin layers can be improved effectively (Eidsvik et al. 2004;Braak 2006;Bosch et al. 2007; Grana and Della Rossa 2010;Figueiredo et al. 2018). The statistical inversion based on the Markov chains Monte Carlo model (MCMC) is usually suggested to solve the seismic AVO inverse problem with strong nonlinearity and evaluate the uncertainty of the estimated model parameters simultaneously (Hansen et al. 2006;Grana and Della Rossa 2010;Alemie and Sacchi 2011;Yuan et al. 2015;Yin and Zhang 2014;Yin et al. 2016;Li et al. 2017aLi et al. , b, 2019. However, for MCMC-based seismic AVO inversion, it is difficult to give the proposed distribution and the optimization direction of candidate states. Furthermore, the computational efficiency of the independent simulation driven by single-Markov chain AVO inversion is lower than that of multi-chains AVO inversion (Li et al. 2019(Li et al. , 2020. When seismic waves propagate in hydrocarbon-bearing porous media, there is a strong nonlinear relationship between P-wave AVO reflectivity and the rock-skeleton moduli, rock-matrix moduli, porosity and fluid modulus (Li et al. 2020). AVO inversion and fluid identification methods are usually developed with linearized P-wave AVO approximate reflectivity including first-order partial derivatives or statistical rock physical models (Grana and Della Rossa 2010;Russell et al. 2011;Yin and Zhang 2014;Zong et al. 2015;Figueiredo et al. 2018;Lang and Grana 2018). From the perspective of Biot-Gassmann's poroelasticity, an improved AVO statistical inversion method driven by differential evolution-Markov chain Monte Carlo model (DE-MCMC) is introduced. A novel AVO reflectivity equation in terms of porosity, rock-matrix moduli, rock density and pore fluid bulk modulus is initially derived. Furthermore, the differential evolution-Markov chain Monte Carlo model is introduced into the stochastic simulation of the seismic AVO inverse problem (Li et al. 2020). And then, the statistical estimation of rock porosity, rock-matrix modulus, density and fluid modulus is achieved through time-complex frequency-domain Bayesian AVO stepwise inversion. The simultaneous optimizations of multiple Markov chains in DE-MCMC model can realize the multiple simulations of model parameters, which is more efficient than the conventional MCMC-based AVO inversion. Finally, models and field data example demonstrate that the statistical inversion of porosity, rock-matrix moduli and pore fluid modulus is helpful in reservoir characterization and fluid discrimination. Biot (1941Biot ( , 1956) established a static mechanical model of fluid-bearing porous media and described the theory of elastic wave propagation in saturated porous media. The elastic characteristics of homogeneous isotropic saturated rocks are described as Gassmann theory (1951). Biot and Gassmann theories are widely applied in seismic reservoir characterization and fluid identification. P-wave and S-wave velocity (V p , V s ) of seismic wave are expressed in terms of dry-rock skeleton (K d , μ d ), rock-matrix moduli (K m , μ m ), rock porosity (ϕ) and fluid modulus K f as,

AVO reflectivity equation with poroelasticity
where K sat is the bulk modulus of saturated rocks. μ sat is the shear modulus of saturated rocks. ρ sat is the density. K d is the dry-rock bulk modulus. α 2 is Biot coefficient. M is a constant related to effective porosity ϕ. K m is rock-matrix bulk modulus. α 2 is Biot coefficient and K f is fluid bulk modulus. μ d is the dry-rock shear modulus. And, With critical porosity model of clastic rocks (Nur et al. 1998), the bulk modulus and shear modulus of dry-rock skeleton in terms of rock-matrix moduli, effective porosity and critical porosity ϕ 0 are where K m and μ m are bulk modulus and shear modulus of the rock-matrix, respectively. ϕ 0 is the critical porosity of clastic rocks, depending on the sorting and rounding of solid particles during the rock-deposition process. ϕ c is the relative porosity and the ratio of rock porosity to critical porosity ( c = 0 ). To implement the simultaneous estimation of rock-matrix moduli (K m and μ m ), dry pores (ϕ) and pore fluids (K f ) with seismic AVO inversion, the model parameters to be inverted are set as = K f , K m , , m , sat T . Substituting Eqs.
A three-layer theoretical model as shown in Table 1 is set up to verify the accuracy of AVO Eq. (8). The top, middle and bottom layers are mudstone, gas-bearing sandstone and mudstone, respectively. The preset rock physical parameters and elastic parameters of this model are shown in Table 1. P-wave AVO reflection coefficients of exact Zoeppritz equation, Gray equation and the proposed Eq. (8) are displayed in Figs. 1 and 2, respectively. And they compare the AVO reflection coefficients and the approximate errors of different AVO equations at the top and bottom reflection interfaces, respectively, which demonstrate that derived Eq. (8) has higher consistency with the exact Zoeppritz equation and Gray linear AVO approximation. With the increase in incidence angles θ, the approximate error also increases gradually. However, when the P-wave incident angle less than 35°, the approximate error of AVO reflection coefficients is still well controlled, and this meets the requirements of seismic AVO inversion.

Sensitivity analysis of porosity, rock-matrix moduli and fluid modulus
The inversion feasibility of the model parameters is discussed with Eq. (8) with models. In Figs. 3, 4, 5, 6 and 7, the sensitivities of model parameters (including porosity ϕ, rock-matrix bulk modulus K m , rock-matrix shear modulus μ m , rock density ρ sat and fluid bulk modulus K f ) to AVO reflection coefficients at the top reflection interface are displayed, respectively. It can be seen that the porosity ϕ has the greatest influence on AVO reflection coefficients, followed by the rock-matrix moduli (K m and μ m ) and density, and fluid bulk modulus K f has the least influences on AVO reflection coefficients. Besides, the porosity, rock-matrix moduli and density of different rocks (mudstone and gas-bearing sandstone) have the same magnitude of influence on AVO reflection coefficients. However, the influence of fluid bulk modulus K f of the top mudstone on the AVO reflection coefficients is different from that of the gas-bearing sandstone in Fig. 7a, b, respectively. Because the porosity of gas-bearing sandstone (ϕ = 0.20) is higher than that of mudstone (ϕ = 0.10). Therefore, the change of the pore fluid modulus illustrates a greater impact on the elastic modulus (K sat ) and density (ρ sat ) of gas-bearing sandstone with high porosity, and a smaller impact on the elastic modulus and density of the low-porosity mudstone leads to a greater impact of fluid bulk modulus of the gasbearing sandstone on the AVO reflection coefficients. From those, it can be concluded that the feasibility of the inversion decreases in turn for porosity, rock-matrix moduli, density and fluid bulk modulus.

Method of statistical AVO inversion
Incorporating seismic wavelets with Eq. (8) as forward solver, an improved AVO statistical inversion by utilizing where ̃ PP is AVO reflection coefficients at different incident angles, respectively.
M is the number of incident angles. R m is the reflectivity of (9) PP =̃ ( ) ⋅ the model parameters to be inverted. ̃ ( ) is the mapping relationship between the AVO reflection coefficients (Li et al. 2017b). Since seismic inversion in complex frequency domain has the advantage of low-frequency recovery, the AVO forward model in the complex frequency domain is introduced as (Zong et al. 2017;Li et al. 2017a, b), where + i is complex frequency. + i , i is the complex spectrum of seismic AVO data.
is the angle-dependent wavelet. Substituting Eq. (9) into Eq. (10) yields the mapping equation between the complex frequency-domain seismic AVO data and the reflectivity of model parameters as, Reflectivity = 1.9 g/cm 3 = 2.0 g/cm 3 = 2.1 g/cm 3 = 2.2 g/cm 3 = 2.3 g/cm 3 = 2.4 g/cm 3 = 2.5 g/cm 3 ρ ρ ρ ρ ρ ρ ρ sat = 1.9 g/cm 3 sat = 2.0 g/cm 3 sat = 2.1 g/cm 3 sat = 2.2 g/cm 3 sat = 2.3 g/cm 3 sat = 2.4 g/cm 3 sat = 2.5 g/cm 3 ρ ρ ρ ρ ρ ρ ρ sat sat sat sat sat sat sat Fig. 6 Influences of density on P-wave AVO reflection coefficients at the top reflection interface. a The influences of density of the top mudstone on AVO reflection coefficients and b the influences of density of the gas-bearing sandstone on AVO reflection coefficients where l is number of the input frequency components ω. k is the number of attenuation coefficients σ. N is the number of sampling points.
i is the complex-frequency spectrum at incident angle i is the wavelet matrix in complex frequency domain with i n c i d e n t a n g l e θ i .
where F is the AVO forward operator in the complex frequency domain. When the input attenuation coefficients σ = 0, Eq. (12) is simplified as the frequency-domain AVO forward model. The time-domain seismic AVO data are also the conditional data in AVO statistical inversion. The timedomain AVO forward model is Markov chain as an example, and defining the current iteration number as i and the current state as i (k) , the candidate state (i+1) * (k) of the kth Markov chain at (i + 1)th iterations is (Li et al. 2019(Li et al. , 2020, where i (x) and i (z) are two different Markov chains to be selected (i.e., the xth and zth Markov chains). N 0, r 2 is a Gaussian PDF with mean of 0 and variance of r 2 that is a small symmetric disturbance introduced in the differential mutation operation. F i+1 k =F min ⋅ 2 exp 1− Q Q+1−i is the mutation scale of the (i + 1)th iteration (Li et al. 2019(Li et al. , 2020. F min is the minimum mutation scale F Q k of the kth Markov chain at the end of iteration. The mutation scale with exponential decay not only ensures the diversity of candidate states, but also avoids the destruction of the optimal state of Markov chain by large mutation scale. The posterior PDF contains exponential and continuous multiplication operation. Therefore, the calculation stability and accuracy of p | , are low. With the natural logarithm on both sides of Eq. (14) and the hard constraints of model parameters, Eq. (14) can be simplified as follows, where O | , is the equivalent objective function of posterior PDF p | , . is the hard constraint data of relative model parameters.
= ⋅ . U is the hard constraint forward operator, where P is the integration matrix of model parameters. With Eq. (17), the acceptance probability of the candidate state (i+1) * (k) is given as, Then, a random variable u is generated based on a uniform distribution Uniform(0, 1) . If ln (u) ≤ , the candidate state is updated i+1 Through the simultaneous optimizations of multiple Markov chains, multiple optimal solutions are obtained simultaneously, and the statistical characteristics including the mean, variance and confidence interval are obtained. Thus, the optimal solution and uncertainty of the model parameters can be evaluated effectively.
In view of the posterior PDF or objective function shown in Eq. (16), AVO inversion in different domains is discussed

Model examples
To demonstrate the feasibility of the proposed statistical AVO inversion method, the resampling logging data (ϕ, K m , μ m , ρ sat and K f ) are utilized to implement the model tests, as shown by the blue lines in Fig. 9a. Figure 8a synthetic AVO data with no noise and S/N = 2 that are synthesized by a zero-phase 35 Hz Ricker wavelet, respectively. In Fig. 8, black lines are the seismic AVO data calculated using real logging data and red lines are the seismic AVO data calculated using the estimated model parameters.
Due to the lack of low-frequency component in seismic data, it is necessary to compensate the low-frequency background of the model parameters in AVO inversion. Furthermore, a reliable low-frequency background is required to estimate the absolute model parameters of subsurface media. Firstly, seismic AVO inversion in the complex frequency domain driven by the DE-MCMC model is developed to estimate the low-frequency background of model parameters. The input frequency components are f ∈ [1.95, 11.72 ]Hz and ∈ 0, 3.6 , as shown in Figs. 9a and 10a, respectively. From the figures, we can see that the estimated lowfrequency model parameters (gray lines, 50-simulation results) under different signal-to-noise ratios (Figs. 9a, no noise, and 10a, S/N = 2) maintain good consistency with the large-scale information of the real logging data. The mean solutions (red dashed lines) of 50 simulations can reflect the large-scale background of model parameters effectively, which provides reliable prior information for the estimation of absolute model parameters through time-frequency joint domain seismic AVO inversion.
Based on the estimated low-frequency model parameters in the complex frequency-domain AVO inversion, the absolute model parameters are obtained by time-frequency joint domain AVO inversion, as shown in Figs. 9b, c and 10b, c, respectively. In the Figs. 9b and 10b, blue lines are the resampling logging data, gray lines are the 50 simulation results, red lines are the mean solution of 50 simulations, and black lines are the 95% confidence intervals of estimated model parameters. The figures show that the final estimated results of model parameters (the mean of 50 simulations) are in good agreement with the well-logging data. Furthermore, the overall inversion errors shown in Figs. 9c and 10c are well controlled. The relative errors only exist in the areas where the model parameters change sharply, which might be caused by the approximate errors of AVO reflection coefficients. The black arrows in Figs. 9b and 10b indicate two gas-bearing sandstones with high porosity. The anomalies of high porosity ϕ, low density ρ sat and low-fluid bulk modulus K f are also illustrated in the estimated results, which verify the effectiveness of the proposed method in quantitative characterization of hydrocarbon reservoirs.
The inversion accuracy of porosity ϕ, rock-matrix moduli (K m and μ m ), rock density (ρ sat ) and pore fluid moduli (K f ) is further compared and analyzed. From Figs. 9b and 10b, the variances of porosity (ϕ) and rock-matrix moduli (K m and μ m ) are smaller than those of rock density (ρ sat ) and pore fluid modulus (K f ). The dynamic ranges of the estimated porosity (ϕ) and rock-matrix moduli (K m and μ m ) are relatively concentrated, and the confidence intervals are relatively narrow, and the uncertainties of inversion results are small. In addition, the inversion accuracy of porosity and rock-matrix moduli is almost the same with each other. Similarly, the accuracy of density is the same as that of pore fluid modulus. The uncertainty of pore fluid modulus K f is the strongest among these model parameters. Figure 11 shows the stochastic simulation processes and convergence curves of the model parameters at 1700 ms in gas-bearing sandstone with high porosity. In this test, 50 Markov chains are set up and the stable simulation of the posterior PDFs could be realized through 6000 iterations of model parameters. In terms of the convergence rates, porosity is the fastest, followed by rock-matrix moduli, and density and fluid bulk modulus are the slowest among these model parameters. From the simulation results of the posterior PDFs (Fig. 11), the same conclusion can be drawn that the posterior means of porosity and rock-matrix moduli are more accurate than those of density and fluid bulk modulus. From those, the theoretical analysis results of the seismic AVO reflection coefficients are verified effectively, which is consistent with the previous analyses of the sensitivities of model reflectivity to seismic AVO reflection coefficients.

Field data example
A filed data example is demonstrated to verify the application effect of the proposed inversion approach in reservoir characterization and fluid identification. In Fig. 12, the incident angles of partial-angle stacking seismic profiles after AVO amplitude-preserved process are 3°-10° (5°), 6°-18° (12°), 16°-24° (20°), 20°-32° (26°) and 26°-36° (31°), respectively. When the P-wave incident angle is less than 35°, the approximate errors of the proposed AVO equation are small. And the real seismic data meets the requirements of incident angles range for seismic AVO inversion.
The estimation of the low-frequency models is firstly implemented. The input frequency components are Fig. 9 Inversion results and errors of model parameters without noise interfere (including porosity ϕ, rock-matrix bulk modulus K m , rockmatrix shear modulus μ m , density ρ sat and fluid bulk modulus K f ). a The estimated results of low-frequency model parameters based on complex Laplace domain AVO inversion, blue lines are the resampling logging data, green lines are the initial low-frequency models, gray lines are the 50 simulation results, red dashed lines are the mean solution of 50 simulations, and black lines are the 95% confidence intervals. b The estimated absolute model parameters based on timefrequency joint domain AVO inversion, blue lines are the resampling logging data, gray lines are the 50 simulation results, red dashed lines are the mean solution of 50 simulations, and black lines are the 95% confidence intervals. c The inversion errors of absolute model parameters gray lines are the errors of 50 simulations, and red lines are the errors of mean solutions  . Figure 13a-e shows the low-frequency model results of porosity ϕ, rock-matrix moduli (K m , μ m ), density ρ sat and pore fluid modulus K f , respectively. Three developing wells in the inversion profiles are well L, well M and well R from left to right, and the logging data is 15 Hz low-frequency filtering result. It can be seen from Fig. 13 that the inversion low-frequency models maintain good consistencies with the large-scale information of subsurface media. Especially, the inversion results of boreholeside traces are in good agreement with the low-frequency filtering data (f c < 15 Hz) of the real logging data, which can reflect the overall longitudinal trends of model parameters. These low-frequency models can provide a reliable priori for the next step of time-frequency joint domain AVO statistical inversion. The low-frequency prior models in Fig. 13 are utilized as the initial models of time-frequency joint domain AVO inversion in the second step and the final inversion results of model parameters are shown in Fig. 14a-e, respectively. In the inversion profiles, the displayed logging data are 160 Hz low-frequency filtering data. The black arrows indicate the positions of gas-bearing sandstones. With the constraints of the prior low-frequency models in Fig. 13, it can be seen from Fig. 14 that the inversion results of model parameters restore the small-scale information of subsurface media effectively and maintain good lateral continuity and inversion resolution. The inversion results of the model parameters of borehole-side traces are in good agreements with the 160 Hz low-frequency filtering results of developing wells. In the gas-bearing sandstone reservoirs, the final inversion results show obvious anomalies of high porosity ϕ, high rock-matrix moduli (K m and μ m ), low density (ρ sat ) and low pore fluid modulus (K f ), which are consistent with the logging interpretation results and rock physical analysis.
To elaborate the accuracy of statistical AVO inversion, the posterior probability density distributions and inversion errors of model parameters at well L, M and R are established by 20 stochastic simulations, as shown in Figs. 15 and 16. The black lines are the elastic parameters calculated with the real logging data, the red lines indicate the 95% confidence intervals of model parameters, and gray lines are the inversion errors of 20 simulations. It can be seen that the posterior PDFs have the highest probabilities at the location of real logging curves, and the posterior mean of the posterior PDFs matches well with the logging data. In the gas-bearing sandstone reservoirs indicated by a series of white arrows, the porosity ϕ shows abnormally high value, which illustrates that the physical properties of reservoirs are good. Furthermore, the fluid bulk modulus K f shows low value and indicates the fluid-bearing properties well of the high-porosity sandstones. In Fig. 15, the real logging data of fluid bulk modulus K f at gas-bearing reservoirs can fall within the 95% confidence interval. In Fig. 16, the inversion uncertainties of the model parameters could be quantified intuitively. Besides, the inversion accuracy of fluid bulk modulus at gas-bearing reservoirs is higher than that of fluid modulus at shale stones. In summary, the proposed AVO inversion can be utilized for the uncertainty evaluations and reservoir characterization.

Conclusions
Based on Biot-Gassmann's poroelasticity, the DE-MCMC model-based statistical AVO inversion method is proposed and it can realize the simultaneous estimation of rock porosity, rock-matrix moduli, density and fluid modulus, which can be utilized in reservoir characterization and pore fluid discrimination. The main conclusions are as follows: 1. With Biot-Gassmann's theory, the seismic AVO equation characterized by porosity, rock-matrix moduli, density and fluid bulk modulus is derived.   3. Incorporating Biot-Gassmann's theory, critical porosity model and DE-MCMC model, the proposed AVO inversion approach realizes the simultaneous estimation of porosity, rock-matrix modulus, density and fluid modulus. The inversion of porosity and rock-matrix moduli is helpful to reservoir characterization, and the direct statistical inversion of fluid bulk modulus is helpful to guide seismic fluid identification. 4. The selection of initial population size and the number of iterations of Markov chains need to be set in advance to ensure that the model parameter with the slowest convergence speed has high convergence accuracy. Some theoretical tests are essential to determine the number of the selected Markov chains when outputting the final inversion results.
Based on the Zoeppritz equation and Aki-Richards linear AVO approximation, Gray et al. (1999) proposed an approximate formula for AVO reflection coefficients characterized by bulk moduli, shear moduli and density of saturated rock, as following, where a( ) = 1 4 − 1 3 2 sec 2 , b( ) = 1 2 sec 2 3 − 2 sin 2 , c( ) = 2−sec 2 4 , θ is the incident angle, γ is the ratio of P-wave and S-wave velocity. Substituting Eqs. (26) and (27)  are the fluid modulus reflectivity, rock-matrix bulk modulus reflectivity, porosity reflectivity, rock-matrix shear modulus reflectivity and density reflectivity at the reflection interface, respectively. The angle-dependent coefficients are written as, Fig. 12 Partial-angle stacking seismic profiles acquired in an exploration area. a Seismic profile with incident angles of 3°-10°. b Seismic profile with incident angles of 6°-18°. c Seismic profile with incident angles of 16°-24°. d Seismic profile with incident angles of 20°-32°. e Seismic profile with incident angles of 26°-36° the joint probability density distribution of observation data (S and d). The prior Laplace probability model p ; R mi , 2 R mi is given as, where R mi is the model parameter at the ith sampling point,

R mi
and 2 2 R mi are the prior mean and variance of the Laplace PDF, respectively, R mi determines the center position of Laplace PDF, and R mi controls the steepness and confidence