Petrophysical seismic inversion based on lithofacies classification to enhance reservoir properties estimation: a machine learning approach

For estimation of petrophysical properties in industry, we are looking for a methodology which results in more accurate outcome and also can be validated by means of some quality control steps. To achieve that, an application of petrophysical seismic inversion for reservoir properties estimation is proposed. The main objective of this approach is to reduce uncertainty in reservoir characterization by incorporating well log and seismic data in an optimal manner. We use nonlinear optimization algorithms in the inversion workflow to estimate reservoir properties away from the wells. The method is applied at well location by fitting nonlinear experimental relations on the petroelastic cross-plot, e.g., porosity versus acoustic impedance for each lithofacies class separately. Once a significant match between the measured and the predicted reservoir property is attained in the inversion workflow, the petrophysical seismic inversion based on lithofacies classification is applied to the inverted elastic property, i.e., acoustic impedance or Vp/Vs ratio derived from seismic elastic inversion to predict the reservoir properties between the wells. Comparison with the neural network method demonstrated this application of petrophysical seismic inversion to be competitive and reliable.


Introduction
Different approaches for petrophysical properties prediction from elastic properties are routinely employed in oil and gas reservoirs (Doyen 2007;Bosch et al. 2010;Grana and Della Rossa 2010;Figueiredo et al. 2018). Empirical equations, geostatistical methods, multi-attribute regression and neural network, petrophysical seismic inversion and co-simulation after stochastic inversion are the predominant procedures (Doyen 1988;Bortoli et al. 1993;Russell et al. 2011;Lang and Grana 2018). Each methodology for reservoir properties prediction possesses deficiency in terms of demonstrating proper relations between elastic and petrophysical properties per each lithofacies class. For instance, fitting an experimental polynomial equation to acoustic impedance (AI) versus porosity (Phi) cross-plot, will not always result in a proper match between the measured and the predicted porosity due to low correlation in some lithofacies classes, e.g., shale and coal. To overcome this issue in the industry methods such as geostatistical methods and neural networks are employed (Doyen 2007;Grana 2020).
Machine learning (ML) is an application of artificial intelligence using data which enables system to learn without explicit programming (Alpaydin 2010;Zimek and Schubert 2017). In the proposed workflow, seismic inversion and petroelastic modeling are incorporated to estimate reservoir properties in a machine learning approach.
The seismic forward modeling is the process of transferring subsurface geological properties into the seismic response (Lang and Grana 2018). Seismic data inversion is a technique through which an elastic property such as acoustic impedance is estimated from seismic response (Tarantola 2005;Sen 2006;Sayers and Chopra 2009;Sen and Stoffa 2013). This approach can be extended to petrophysical properties prediction derived from elastic properties which is called petrophysical seismic inversion (Doyen 2007;Bosch et al. 2010;Grana and Della Rossa 2010). This method is known by other terms such as petroelastic inversion or inversion of inversion (Mukerji et al. 2001;Coléou et al. 2005;Bornard et al. 2005). The process starts with an initial fine-scale model and updated iteratively. Petrophysical seismic inversion has been developed by applying analytical Bayesian techniques and using different rock physics models (Bosch et al. 2010;Grana et al. 2012Grana et al. , 2016Grana et al. , 2018Figueiredo et al. 2018;Lang and Grana 2018;Grana 2020).
Petroelastic modeling (PEM) is used to convert petrophysical properties to the elastic ones (Avseth et al. 2005;Mavko et al. 2009;Dvorkin et al. 2014). Therefore, at the first step, we need to find a reliable petroelastic model at a well location in order to establish a nonlinear experimental relation for solving an inverse problem algorithm. In this study, PEM is optimized via calibrating the predicted elastic logs with the measured logs by selecting the most suitable pore space stiffness value. This technique aids to reduce uncertainty and validate the results (Babasafari et al. 2020). Next, the seismic elastic inversion is performed. In this paper, we present a practical approach in petrophysical seismic inversion using a stable nonlinear optimization algorithm in which a lithofacies volume is interactively involved.

Petrophysical seismic inversion
Since getting petrophysical properties from seismic data directly is impractical, elastic properties as a key link between seismic and petrophysics should be predicted first. There is an approach  for petrophysical seismic inversion in which the elastic model is generated from the petrophysical model through petroelastic modeling and then seismic forward modeling used to create synthetic seismic data from the elastic model. Recorded and synthetic seismic data are compared via solving an inverse problem algorithm, and the petrophysical model is updated iteratively. Coléou et al. (2005) proposed a detailed technique which is started from an initial fine-scale geomodel defined from a 3D stratigraphic grid. Afterward, a petroelastic model (PEM) is employed to calculate elastic properties at each cell of the geomodel from porosity, volume of clay and water saturation volumes. Using the Zoeppritz equation at each trace location, the angle-dependent reflectivity series is calculated from the elastic properties. 3D angle-dependent synthetic reflectivity series are produced after wavelet convolution. The degree of match between the synthetic and the recorded angle stacks is optimized through introducing perturbations of the properties of the geomodel using a simulated annealing algorithm. After convergence of errors, the final geomodel honors the observed seismic amplitudes and also is consistent with the user-defined PEM . It is worth noting that even minor changes in the initial model such as the small-scale distribution of rock type or change in PEM result in different solutions. The final models due to the nonuniqueness problem in seismic inversion results then represent alternative solutions that are in compliance with the seismic data.
This approach utilizes a robust petroelastic relation as a constraint in an inversion workflow that is finalized on the basis of minimizing the residual error between recorded and synthetic angle-dependent reflectivity series. It is therefore known to be an accurate and reliable technique in reservoir properties prediction.
In the proposed methodology in this research, a workflow is amended as follows; first, the elastic model is derived from an initial petrophysical model via petroelastic modeling which is known as modeled elastic property. As a parallel process, seismic elastic inversion is implemented to convert seismic data into elastic data which is known as inverted elastic property. Elastic properties in this study are acoustic impedance and V p V s ratio, and desired petrophysical models are porosity and clay volume. For instance, in porosity estimation, the modeled acoustic impedance originated from petroelastic model and inverted acoustic impedance derived from elastic inversion are then compared for matching via solving an inverse problem algorithm and the porosity model is updated iteratively. Once the residual error between the inverted and the modeled acoustic impedance reaches the defined criteria for minimizing misfit function, the porosity model is finalized. In this research two different nonlinear optimization algorithms are employed and compared: the Newton-Raphson and the Secant methods. Figure 1 illustrates the workflow of petrophysical seismic inversion for porosity estimation.

Reservoir properties estimation
The petrophysical seismic inversion workflow for porosity estimation is carried out in two main stages: First nonlinear petroelastic experimental relations, i.e., porosity versus acoustic impedance, are fitted at well location per each identified lithofacies class separately. Petrophysical seismic inversion is performed using two different optimization algorithms (Newton-Raphson and Secant). Afterward, the results of predicted porosity are compared with those from measured porosity through cross-plot and statistical analysis.
Once a good match between measured and predicted porosity is attained, the petrophysical seismic inversion is applied to the inverted acoustic impedance produced from seismic elastic inversion to predict porosity between the wells. For this purpose, a lithofacies volume is interactively employed (Ghosh et al. 2018). In the end, results are validated through the use of a blind well and compared with those from the probabilistic neural network algorithm. A similar process is applied to The advantages of the proposed application of petrophysical seismic inversion compared to the Coleou's approach are as follows: (a) Instead of comparing recorded and synthetic seismic data via solving an inverse problem algorithm in the matching process, the modeled and inverted elastic properties are utilized and compared. This helps to mitigate the uncertainty of using seismic forward modeling in Coleou's method in which noise, multiples and migration artifacts should be taken into consideration in order to relate the residual errors more to the petrophysical properties estimation and not to the processing issues. The simulated annealing algorithm was employed to solve the inverse problem in Coleou's approach that is replaced with the Newton-Raphson and Secant algorithms in the present case. (b) Since the petroelastic relations vary from facies to facies, lithofacies classification is interactively utilized to assign an identified lithofacies class to each cell of the geomodel. To this end, first elastic seismic inversion is conducted. The outputs are incorporated in the workflow: acoustic impedance, In this inversion workflow, each petrophysical property is predicted using the more appropriate elastic property in terms of getting a higher correlation coefficient in petroelastic cross-plot, e.g., porosity versus acoustic impedance. Although each petrophysical property is independently estimated, this research tries to demonstrate that the basic reservoir parameters such as porosity and volume of clay can be attributed to the corresponding elastic parameters, i.e., acoustic impedance and V p V s ratio, respectively.
The application of petrophysical seismic inversion in this study is not a stochastic approach. It is a deterministic technique. However, it proposes a reliable tool for crossvalidating of results and reducing uncertainties in reservoir characterization which is the residual error between inverted and modeled elastic property.

Optimization algorithms
Two different nonlinear optimization algorithms are selected and used to minimize the misfit function of the petroelastic inverse problem: Newton-Raphson and Secant.
Newton-Raphson is an approximation method which is based on the Taylor series and requires an initial value, and also a derivative of the nonlinear equation should be known. This iterative procedure can be generalized by writing following equation Eq. (1), where i represents the iteration number:

Results and discussion
Data from a field located in the Malay basin are used for this study. This research focuses on reservoir layer E that is mainly composed of clastic sediments. 3D pre-stack seismic data and four wells including elastic and petrophysical logs are used. Lithofacies classes are created utilizing petrophysical logs and available lithological data (Grana and Della Rossa 2010;Babasafari et al. 2020).
Lithofacies information at the well location is obtained from petrophysical logs (porosity, water saturation and mineral grain volumes). Four lithofacies classes are classified including shale, wet sand, gas sand and coal. Elastic seismic inversion was conducted, obtaining 3D inverted ratio and lithofacies class. First, the results of reservoir properties prediction are illustrated at well location and then the final results are displayed on 2D cross sections. The results consist of porosity and clay volume prediction based on lithofacies classification by utilizing nonlinear optimization algorithm.
The results of predicted porosity are compared with the measured porosity through cross-plot and statistical analysis. Figures 3 and 4 display quality controlling steps. Figure 3a shows overlapped zones between different lithofacies classes, and Fig. 3b-e is used as quality control (QC) for the results. Figure 3b illustrates the overlay comparison cross-plot of measured and predicted porosity to assess the prediction accuracy at well location. Figure 3c-e displays statistical analysis of measured and predicted porosity. Figure 4a illustrates the elastic logs, raw logs and petrophysical interpretation logs. In Fig. 4b the QC is conducted in the depth domain by superimposing measured porosity on the predicted porosity derived from the acoustic impedance log using petrophysical seismic inversion within selected interval. However, in Fig. 6c the QC is performed in two-way time (TWT) domain by superimposing measured porosity on the predicted porosity derived from inverted acoustic impedance that might be influenced by the well to seismic tie in the inversion process. Figure 5 shows the seismic elastic properties, i.e., acoustic impedance and V p V s ratio in 2D cross sections as the input of petrophysical seismic inversion. Figure 6 represents lithofacies classification and predicted porosity. In the cross section one constrained well is located in common depth point (CDP) no. 215 and blind well in CDP no. 31. According to Fig. 6c the correlation coefficient between predicted and measured porosity logs is 71%. In Fig. 7 the residual error between the modeled and the inverted acoustic impedance is displayed. Apart from the yellow color and the dark blue zones in Fig. 7b that are attributed to shale lithology, the rest including the reservoir interval demonstrate a residual error value close to zero. High values of residual error in the mentioned locations (yellow and dark blue color zones) are due to the divergence of the residual error miscalculation (overfitting) where inverted acoustic impedance value does not fit the selected nonlinear equation (porosity vs. acoustic impedance) for shale lithology. This can represent uncertain zones for the mentioned method. Meanwhile, the accumulation of errors is quantified in Fig. 7c where the Newton-Raphson and Secant methods are compared. It is not a stochastic approach; therefore, the uncertainty is not specifically quantified based on an objective criterion of a couple of realizations and only the average error is evaluated. Since the Newton-Raphson method reaches the Fig. 2 Approximation for the Newton-Raphson method (Urroz 2004) minimum criteria of misfit function via a lower number of iterations than the Secant method, it was selected to be applied for the final petrophysical seismic inversion (Fig. 7c). The result of the blind well in this method was compared with the probabilistic neural network method and validated with higher correlation coefficient and lower root-mean-square error (RMSE) values in Table 1. The correlation coefficient between measured porosity log and predicted one at blind well using two different methods, petrophysical seismic inversion and probabilistic neural Fig. 3 a Cross-plot of acoustic impedance (m/s * gr/cc) versus porosity (%) at well location color coded by lithofacies classes. b Overlay comparison cross-plot of acoustic impedance versus porosity, meas-ured porosity in black and predicted porosity in red color. c Predicted porosity versus measured porosity. d and e Distribution histogram of measured and predicted porosity, respectively 1 3 network that are 48.27% and 33.78%, respectively, is illustrated in Table 1. The matching quality is not acceptable for both methods at blind well. However, this is the only available well for cross-validation and the seismic noise effect might not be neglected. In the next step of petrophysical seismic inversion process, a nonlinear petroelastic experimental relation, i.e., clay volume vs. V p V s ratio, is fitted at well location per each identified lithofacies class separately (Fig. 8). Finally petrophysical seismic inversion is applied to the inverted V p V s ratio extracted from seismic elastic inversion to predict volume of clay between the wells (Fig. 9). In the same process, the results of predicted clay volume are compared with the measured clay volume through cross-plot and statistical analysis. Figure 8 displays the quality control steps at well location. Figure 9 demonstrates predicted clay volume in 2D cross section.
This method of petrophysical seismic inversion needs some requirements as follows: • The experimental relationship between elastic and petrophysical properties derived from well logs should be defined at well location. • Lithofacies volume needs to be involved.
Compared to the other approaches, the proposed application of petrophysical seismic inversion has its own advantages as indicated below: • Statistical analysis is used to attain a similar probability distribution function (PDF) trend of measured objective (reservoir property) for the predicted reservoir property. Fig. 4 a Elastic logs, raw logs and petrophysical interpretation logs. b Measured porosity superimposed on predicted porosity derived from petrophysical seismic inversion color coded by identified lithofacies classes within selected interval of well • Unlike neural network method, the residual error between inverted and modeled elastic property away from the wells validates the result as a remarkable tool for quality control of optimization algorithm.
• This application of petrophysical seismic inversion benefits from seismic elastic inverted results; however, unlike the co-simulation of petrophysical properties in geostatistical inversion, there is no variography effect in outcome. The main disadvantage of the proposed workflow is that the application of this approach is not a multivariate joint inversion of petrophysical properties in which a set of nonlinear equations are incorporated and reservoir properties are simultaneously estimated.

Conclusion
This research demonstrated how the proposed application of petrophysical seismic inversion can be employed to achieve useful petrophysical properties which are then utilized in Fig. 6 a Cross section of lithofacies classification. b Cross section of predicted porosity using petrophysical seismic inversion method and c well log porosity superimposed on predicted porosity using petro-physical seismic inversion method extracted from model representing by white arrow in cross section (a) and (b) reservoir modeling. The method reduces uncertainty in reservoir characterization. The residual error between inverted and modeled elastic property is a reliable tool for cross-validating of results and reducing uncertainties. Unlike the co-simulation of petrophysical properties in geostatistical inversion, the variography effect is not observed in the outcome. Estimation of the reservoir properties (porosity and clay volume), using this approach, shows above 70% correlation with measured logs. The calculated RMSE values at a blind well are 0.017 Fig. 7 Residual error between modeled and inverted acoustic impedance a cross-section of first residual error b cross section of final residual error and c average residual error for different iteration numbers using two methods of Newton-Raphson and Secant at well location and 0.027 for petrophysical seismic inversion and probabilistic neural network techniques, respectively. The outcome of the proposed technique compared to the results obtained from available advance methods in the industry, e.g., probabilistic neural network, is proven to be competitive and reliable. However, this method is not a multivariate joint inversion of petrophysical properties. The petrophysical seismic inversion can be extended to multivariate joint inversion workflow and simultaneous prediction of multiple reservoir properties (porosity and clay volume) through using required equations in petroelastic modeling.