The seismic inversion method combined with multipoint geostatistics theory has begun to receive attention, but the acquisition accuracy and calculation efficiency of 3D training image still need more optimization. This paper presents a novel method of 3D multipoint geostatistical inversion based on 2D training images directly. The 2D training image was scanned by the data template to acquire the multipoint statistical probability in 2D direction. The probability fusion method is used to fuse the 2D multipoint probability into 3D multipoint probability. The rock facies types and patterns of the simulated points are obtained by random sampling. On this basis, the elastic parameters are extracted from the statistical rock physics model, and the seismic records are convoluted. Then, the synthetic records and the actual records were compared under a given threshold. If the error exceeds the given threshold, the iterative adaptive spatial sampling method will be used to repeat the process above-mentioned, so as to ensure that the error is below the threshold. Because the 2D training image is easy to obtain and evaluate, the demand problem of 3D training image is solved. The 2D training image scanning, probability storage and access are more convenient, and the adaptive spatial sampling method is more efficient than the reject sampling, so as to ensure the operation efficiency. The model from the Stanford Center for Reservoir Forecasting is selected to test the effectiveness of this newly designed method.

Seismic inversion is an important petrophysical way to predict the subsurface reservoir during the exploration and development phases of almost every oil and gas field [1, 2]. Early inversion often gives unique result, which terms deterministic inversion; later, the stochastic inversion is introduced to give multiple results for risk assessment. Due to the limitation of seismic wave frequency, the simplicity of forward modeling, and strong reservoir heterogeneity, even if there is no noise in seismic records, the inversion of seismic records into reservoir attributes is uncertain, and its reasonable interpretation is a challenge [3, 4]. The establishment of statistical rock physics model can characterize the multisolution between petrophysical parameters and reservoir attributes and the uncertainty of inversion to a certain extent. Constrained by the limited seismic bandwidth, the conventional inversion can only reflect the exact resolution of seismic data. It is usually difficult to recover full frequency information. In order to improve the resolution of inversion results, it has become a consensus to add multiple constraints of logging data and geological data in the inversion process. Thus, an inversion method integrating statistical rock physics and geostatistics theory has been gradually formed [3, 5, 6]. Bosch et al. divided these inversion methods into two categories [3]: stepwise inversion method and unified inversion method based on Bayesian formula. The stepwise inversion method does not output elastic parameters and reservoir attributes at the same time, but carries out inversion first, and then applies inversion parameters with clustering or costimulating method to reservoir attribute prediction. Because the elastic parameters and reservoir attribute prediction are studied separately in this method, how to combine the uncertainties in different steps and reflect them in the final result would be a problem. The unified inversion method based on Bayesian formula integrates seismic rock physics and reservoir attributes through a single formula, which better expresses the uncertainty and has become a common seismic inversion method gradually. Obtain the optimal reservoir model by seeking the maximum likelihood function or a posteriori probability or sampling multiple models from a posteriori probability to meet the uncertainty evaluation [7, 8]. The adopted sampling method is usually the reject sampling [9]. However, reject sampling needs evaluating a lot of forward models to obtain a posteriori probability in the statistical sense, and its efficiency is relatively low. An iterative inversion method realizes part of the data through sampling the previous simulation as the next inversion conditional data, which can achieve convergence much faster, but it is easy to fall into local optimization. As the sampling is random, the points with unexpected errors will be accepted and these errors have a great impact on the accuracy of subsequent inversion. At the same time, it also increases the number of iterations which reduces the calculation efficiency. Jeong et al. [8] proposed an adaptive spatial sampling method to solve this problem. By comparing the error between the inversion results and the actual observation data, the sampling points are controlled in the area with less errors, which accelerates the posterior probability to reach the steady state, while the inversion model is closer to the reality.

The prior geological model in the Bayesian inversion framework is mainly obtained by geostatistics. The traditional two-point statistical method describes the characteristics of parameter space by variogram function and obtains the conditional probability by solving Kriging equation. The algorithm is efficient, but it is impossible to describe the complex curvature reservoir shape and the spatial superimposition mode due to the single point (track) sampling [6, 10]. In the field of geological modeling, a variogram-based modeling method is gradually replaced by multipoint geostatistics method. Multipoint geostatistics reveals the spatial variability and correlation of reservoir through three dimensional (3D) training images. This prior information can be obtained through multipoint data template scanning, which can better reproduce the reservoir shape and superposition style. Replacing the variogram-based geostatistical inversion by multipoint geostatistics one has become a research hotspot. The first developed multipoint geostatistical inversion algorithm has achieved good results in theoretical and practical application [4, 7, 10]. However, the multipoint geostatistics method faced the unavoidable problems, that is, the accurate acquisition of 3D training images and computational efficiency [3, 10]. Although scholars have carried out lots of research on the acquisition, evaluation, and computational efficiency of 3D training images [11, 12], the complexity of subsurface reservoir structure makes the training image matching more difficult, and the 3D image scanning greatly increases the time consumed on computation.

In geological research, the distribution and superposition style of geological bodies are often reflected by two-dimentional (2D) plane facies map and section facies map. The compilation of these maps not only considers well data, seismic data, dynamic response, and other data but also makes full use of the constraints provided by modern sedimentation, ancient outcrops, similar sedimentary areas, UAV detection, and other achievements. Finally, based on the experience of geologists, a conceptual model of the distribution of reservoirs is formed, which is the deep combination of objective data and subjective cognition. So the 2D training image can represent the subsurface structure. Reconstruction of geological model based on 2D training images using multipoint geostatistics has become a hot spot in geological modeling and has been tried and applied in reservoir modeling. The reconstruction of 3D geological model using 2D images was first applied for the construction of pore throat model by Okabe and Blunt [13, 14]. Comunian et al. [15] investigated the 2D probability fusion mechanism in detail and designed the s2Dcd method for subsurface reservoir prediction. Chen et al. [16] proposed the reconstruction method of the local partition and search strategy on the basis of s2Dcd to meet the local stationarity during probability fusion. Wang et al. [17, 18], aiming at solving the problem of poor continuity of 2D reconstructed geological body, introduced the adaptive spatial sampling method to retain the data with small observation error as the conditional data for the next modeling, so as to improve the connectivity of geological body and the accuracy of model. But using the 2D constructing method for 3D seismic inversion has not been proposed.

In this paper, a 3D multipoint geostatistical inversion method based on 2D sedimentary facies profile is proposed. For each simulated grid, the given data template is used to scan the 2D sedimentary facies section to obtain the 2D multipoint probability. The probability fusion strategy is used to generate the 3D multipoint geostatistical model. Combined with the statistical rock physics model, the sedimentary facies and its corresponding elastic parameters are sampled. The seismic records are synthesized by convolution theory formula, and a comparison with the actual seismic records is conducted. The sedimentary facies and elastic parameters are sampled, and the best elastic parameters are retained through adaptive spatial sampling method, so as to realize the synchronous output of lithology and elastic parameters and achieve the purpose of seismic inversion. The method is compared and tested using the data set from the Stanford Center for Reservoir Forecasting.

Seismic inversion is to sample the elastic parameters and to calculate synthetic seismic records by convolution, which matches actual seismic records under the given threshold. The most known inversion method can be expressed by the Bayesian formula [19].
(1)σMm=cγMmγDfm,
where c is a constant for normalizing the posterior probability, γMm is the prior probability of reservoir property, and γDfm is the likelihood function to measure the bias between the synthetic results and the observation records. M is the working area, and m is the initial property, such as the elastic parameters, the facies categories. fm is the forward operator, which may be a simple function with an analytic form or computational algorithm with no simple analytic expression (3). σMm is the posterior probability which is updated from the prior probability.

The solution approaches of seismic inversion are divided into two types, that is, to search the maximum posterior probability density or to sampling the best properties from the posterior probability density [3]. As the elastic parameters are generally sampled from the rock physics function and the rock type is forecasted by the prior probability γMm, how to predict γMm is a key point in seismic inversion [20]. Traditional inversion used a variogram-based method such as sequential Gaussian simulation to simulate γMm, and then, the elastic parameters are sampled and the best matching are determined. This is the form of the most initial variogram-based geostatistical inversion [19]. But the variogram function based on two-point statistics is difficult to characterize the correlation characteristics between complex curved geological bodies and attributes. Therefore, using multipoint geostatistics for probability inference has become a natural choice. However, multipoint geostatistics needs to scan 3D training images to obtain multipoint probability. The acquisition of 3D training images is very difficult, and repeated scanning leads to low computational efficiency, which has become the main problem restricting the practical application of multipoint geostatistical inversion. The newly developed method of reconstructing 3D multipoint probability using 2D profiles has played a role in geological modeling. It is a useful attempt to extend it to multipoint geostatistical inversion.

Based on the development of the reconstruction algorithm of 3D model using 2D profiles, a new multipoint geostatistics inversion method based on 2D profile is designed. The multipoint probabilities of different facies combinations are obtained by 2D profile scanning, and the possible rock facies types are obtained by sampling. Based on the facies model, the elastic parameters are sampled using the established rock physics model. The synthetic records are obtained by a forward operator and compared with the input seismic records, and the most matched sampling elastic parameters are retained as the inversion results. This sampling optimization theory ensures that the prior lithofacies model comes from the sampling of equal probability model and meets the requirement that the error between the inversion synthetic record and the actual observation data is within a certain range under a given statistical rock physics model [10]. The method includes four parts: 2D probability acquisition, 2D pattern probability fusion, statistical rock physics relationship, and spatial sampling.

2.1. The 2D Multipoint Probability

The conventional multipoint geostatistical inversion first adopts the multipoint simulation method to scan the 3D training image for the acquisition of probability using data template. However, the 3D training image is difficult to construct, and the calculation efficiency of each scanning is very low. The 2D training image solves this problem well with a probability fusion method. Firstly, the geological analysis section is easy to obtain, and it is the possible situation of subsurface geological conditions by geologists on the basis of data analysis, which is representative and avoid the consideration of accuracy evaluation of training images. Secondly, the scale of 2D training image is smaller than that of 3D training image, and the demand for memory is low. The probability can be stored by search tree to meet the demand of fast access. Finally, the 3D probability can be calculated by probability fusion method. In this paper, the Snesim method [21] is used to calculate the conditional probability on 2D section. For the given direction, 2D data template is used to obtain the conditional data points to form the special 2D data event. Then, the 2D training images are scanned with the special 2D data event in the parallel direction to find the matching one. The numbers of the repetition data event, namely, Ri,jZk, at the evaluated point belonging to the property Zk are counted. Finally, the frequency of each data event is used to represent the 2D conditional probability of occurrence of facies Zk, that is,
(2)Pi,jZk=Ri,jZkk=1vRi,jZk.

Considering the variability of sedimentary patterns at different locations, pattern probability extraction needs to be carried out within confined regions to guarantee the stationary assumption. Chen et al. [16] proposed a local search strategy to solve this problem (Figure 1). The space is divided into 27 subzones by six cross-sections in 3 orthogonal directions. Each direction exists 2 parallel cross-sections. As can be seen in Figure 1, the subzone is surrounded by N subsections; N is between 3 and 6 because in some subzones the area may not be closed. As described by local strategy, the statistical information of the simulated grid is only calculated from the adjacent N subsections, which guarantee the stationarity assumption, and so the statistical information represented the local characteristics than the whole region scanning strategy.

2.2. The Probability Fusion of 2D Section

The probability fusion of 2D section is to fuse the probabilities of sections in different directions. Many scholars have carried out research on the method of aggregate 2D vertical section probability to form 3D probability [14, 15, 22]. The probability aggregation method was divided into addition-based aggregation and multiplication-based aggregation function [22]. The Linear Pooling addition formula was used to aggregate the probability from two independent parallel sections in the same direction due to the similarity of the patterns, and the probability aggregation in different directions uses the Log-Linear Pooling multiplication formula to reflect anisotropy. Wang et al. [17] and Chen et al. [16] depicted the fusion method in detail. The brief idea is as follows.

For the simulated grid, there are many slices in m directions, m is parallel to X or Y or Z direction and is no more than 3. In each direction (X,Y,Z), n adjacent sections are existed (1n2). Then, the probability of point x can be calculated by scanning the sections (see Figure 2). First the conditional data events around point x of the given direction is obtained, then the 2D training images adjacent to that direction are scanned to obtain matching data events, thirdly the conditional probability density function PZx of the point x is calculated using the occurrence frequency approximation, and finally, the aggregate probability PtZx is obtained using the addition formula, as shown by [17]
(3)Pt,iZX=j=1nωi,jPi,jZXi=1,2,,m,
where ωi,j=li,j/j=1nli,j as the probability in the parallel direction is fused using equation (3). The probabilities in orthogonal directions then are fused using a multiplication formula, and the final 3D conditional probability density is forecasted, as
(4)PfZXP0ZX1i=1mwii=1mPt,iZXwi.
P0ZX is a priori probability, and wi is the weight values in different directions, which can be equal or different. In real case, the inverse distance function is usually adopted to assign the weight value, and the formula is as follows:
(5)wi=1Lii=1mLi,
where Li=minli,1,li,2,,li,n.

2.3. Statistical Rock Physics Model

A rock physics model plays an important role to connect elastic parameters with reservoir properties. Theoretical and empirical rock physics model believes that rock modulus is affected by mineral type, porosity, pore throat type, cement, mineral type, sorting, and pressure. This paper does not discuss the types of rock physics models, nor does it analyzes the uncertainty caused by the difference between the reservoir attributes interpreted from well data and the scale of elastic parameters obtained from seismic data. The quantitative relationship between rock physics parameters and reservoir attributes is given only from the perspective of statistics. Using this quantitative relationship, the sampling and inversion of elastic parameters under the constraint of reservoir attributes can be realized. The general statistical rock physics method is mainly to extract the elastic parameters of the well bypass, carry out statistical fitting analysis with interpreted sedimentary facies from the well, and establish the statistical function of elastic parameters and sedimentary facies. Figure 3 shows the distribution characteristics of elastic parameters of different sedimentary facies obtained from the test data of Stanford University. In inversion prediction, once the sedimentary facies are obtained through 2D profile prediction, the probability relationship between sedimentary facies and elastic parameters can be used to sample elastic parameters.

2.4. Spatial Sampling

In the process of seismic inversion, iterative spatial resampling is often used to meet the stability requirements of final sampling. This study still adopts the idea of iteration. After establishing the 3D geological model by the 2D sections, sampling elastic parameters is carried out according to the statistical relationship between sedimentary facies and elastic parameters. Then, the forward convolution is performed to get the synthetic seismic records. A comparison between the synthetic one and the actual seismic records is executed and the error is calculated. If the error is larger than the given threshold, iteration is executed and elastic parameters are sampled again. At this time, the idea of adaptive spatial sampling is introduced; that is, in iteration, a certain number of simulation results in the area with small error are retained as conditional data, and they are no longer involved in modeling and inversion; iterative modeling and inversion are only carried out for the regions with a large number of errors, so that the inversion error can quickly reach the given range. The implementation steps are illustrated as follows:

  • (Step 1)

    The cross-sections of the 2D training images are acquired, and the simulation parameters are input, such as the data template, the maximum search radius, the maximum number of conditional data, the error threshold to determine the termination, the sampling area, the retain maximum sampled data points, sedimentary facies, and elastic parameters relation functions

  • (Step 2)

    Select the points to be estimated on the random path, search the conditional points with the template corresponding to each direction, obtain the conditional data events, search the matching pattern in the adjacent training images in the same direction, and establish the probability density function of the corresponding direction

  • (Step 3)

    In the process of probability fusion, the same direction is processed first. If there are two adjacent training images, the probability of each data events in this direction is fused by formula (3). The pattern probability in different directions is fused by formula (4). After obtaining the pattern probability, the sedimentary facies value of the point to be evaluated is determined by random sampling. Repeat steps 2 and 3 until all grid nodes are simulated to obtain the sedimentary facies model

  • (Step 4)

    The elastic parameters are obtained by sampling according to the relationship between statistical rock physics, that is, the relationship between elastic parameters and sedimentary facies, and the forward convolution generates synthetic seismic records, which are compared with the actual seismic records. If the error is less than the given threshold, the sedimentary facies and elastic parameter model are output to complete the inversion. Otherwise, the data points are randomly sampled from the area where the error is less than the threshold as the conditional points for the next simulation, and the simulation returns back to step 2. The workflow is shown in Figure 4. The code is compiled by Intel Fortran, and the result is shown in the commercial software such as Petrel

3.1. Basic Data

The Stanford VI reservoir model was often used to test newly designed geostatistical or inversion methods [10]. The model consists of three different sedimentary system. The sediment source is from the south and propagate to the northern basin. The top layer is deltaic deposits, the middle is meandering channel deposits, and the bottom is sinuous channel deposits. In this study, the top channel deposits were selected as the test object, in which channels and point bars were the main sand body deposits. The model region was gridded to 150×200×80 where the grid size is 25 m in the horizontal direction and 1 m in the vertical direction [18]. 68 pseudowells and 14 2D training image positions were designed (where Nx=3,34,75,125;Ny=21,62,99,142,182;andNZ=3,20,37,54,71); the lithofacies data in the wells and 14 2D cross sections were acquired directly from the open model and were used both as the input data and as 2D training images (Figure 5). 13 test wells were selected for cross-validation to evaluate the precision of the inversion results.

3.2. Comparison with the Existing Model

The method designed in this paper is used for the inversion and prediction of elastic parameter field. The 68 pseudowell data and 14 training images are used as input data to constrain the simulation. The search radius is set to 9×9×5. The maximum conditional points were set to 60. The subzone ratio to be scanned was 0.8, and the termination of iteration is to realize that the error of the proportion of the channel between simulation and public model is within 5%. During resampling, the area with the error between the synthetic seismogram and the original seismogram less than 1‰ was taken as the sampling point area, and the proportion of sampled points occupied the model was less than 25%. The simulation process went through two iterations. The adaptive sampling area is defined as the point set with the distance between the selected data event, and the conditional one was smaller than 3% and the aggregated 3D probability excessed 95%.

From the inversion results, the inversion method based on 2D profiles reproduces the distribution of elastic parameter field well. On the whole, the distribution of density, velocity, and impedance data obtained by inversion is close to the theoretical model (Figure 6). From the cross-section, the distribution of impedance value in the theoretical model is controlled by sedimentary facies. The impedance of point bar is generally low, while the impedance of channel is relatively high and the impedance of mudstone is the highest. There is an obvious high value of impedance data on the contact interface between point bar and channel. In the inversion model, the low impedance of point bar, the relatively high impedance of channel, and the high impedance of mudstone are well reflected, and the connection between point bar and channel is also clear. It reflects that the inversion results match the reality.

From the comparison of statistical distribution, the elastic parameter distribution obtained by inversion is in good agreement with the theoretical model (Figure 7). It shows that this new method is faithful to statistical characteristics. The inverse impedance is subtracted from the actual impedance to analyze its absolute error (Figure 8). The absolute error is largely distributed near the value of 0, but in local areas, there are areas with large errors. On the whole, the impedance error between the inversion model and the actual model is analyzed and counted. The proportion of absolute error less than 0.1 is 87.8%, and the proportion of points with relative error less than 10% is 94%. It shows that inversion can reveal the distribution of actual impedance and has high accuracy.

3.3. Cross-Validation by Wells

It is common to check the accuracy of reservoir prediction through blind wells. The accuracy of 2D profile inversion method is checked by 13 randomly selected wells. From the cross-profile (Figure 9), the impedance of the inversion is very close to the impedance data in actual model. The absolute errors are mainly distributed near 0, but there are a small number of high value areas in local positions on the well. Statistics show that the absolute error less than 0.1 is accounting for 92.9%. The relative error less than 10% is accounting for 96.2%. The statistical results of blind well prediction show that the 2D profile inversion method can predict the distribution of elastic parameters between wells with high accuracy. Therefore, this method can be applied to practical seismic reservoir prediction.

In geological modeling, the prediction of reservoir attributes mainly adopts the prediction method of facies controlled petrophysical simulation; that is, the physical properties of porosity and permeability are controlled by sedimentary facies. Before modeling, it is necessary to analyze the statistics of physical parameters of different sedimentary facies according to well data and establish the distribution function of physical parameters using the control of sedimentary facies. In the petrophysical modeling, the physical property parameters are obtained by sampling from the established physical property parameter distribution function of different facies. Similarly, elastic parameters are affected by physical properties, pore throat structure, clay minerals, and fluid properties, which are mainly controlled by sedimentary facies. Therefore, establishing the relationship between elastic parameters and sedimentary facies for elastic parameter prediction can be adopted in the seismic inversion, which can realize the optimal matching quickly between the forward synthetic records and actual one, and can improve the operation efficiency and the accuracy of inversion prediction.

The accuracy of sedimentary facies prediction has a huge impact on the inversion results. Although 3D multipoint modeling using 2D profile reconstruction avoids the difficulty of obtaining 3D training image, the statistical information in 3D training image is much more abundant than 2D profile. Therefore, 2D profile is at the expense of losing information, which may affect the accuracy of sedimentary facies model. Wang et al. [18] discussed the influence of section numbers on the geological model. As the section number increases, the channel shape and continuity become better. When the real 3D model is used as the training image, the modeling result is the best. However, for the delta reservoir, the outcome of geological modeling directly using the 3D delta training image is not ideal (Figure 10), because the delta reservoir does not have statistical stability; that is, the stationary assumption is not applicable. However, the section based geostatistical modeling just splits the research area, so as to ensure the local statistical stability. Therefore, for nonstationary reservoirs, the effect of section based modeling should be significantly improved. For geostatistical inversion, the traditional method does not carry out zoning simulation, resulting in low inversion quality when the reservoir is complex and nonstationary. The inversion based on 2D profile will potentially solve the problem of nonstationary and can predict the complex reservoir.

As the resampling is conducted by the error between the real seismic record and the synthetic one, the quality of seismic data is very important for inversion. If the quality of seismic record is low, the inversion process may be very slow and even not workable. So in practicality, the quality of the seismic data must be insured and a larger threshold is set for lower resolution data, which guarantees the convergence. However, the uncertainty will be increased, and the model must be handled carefully.

This method is a sampling optimization method. Although the posterior probability distribution is not given, sampling is carried out to produce inversion results. However, previous studies have shown that there is little difference between the sampling optimization method and MCMC inversion results [8], so it does not affect the final inversion.

A novel method of multipoint geostatistical inversion using 2D training images is proposed in this paper. The 3D multipoint probability is generated by 2D training images scanning and probability fusion method, and the sedimentary facies model is constructed by sampling. According to the facies model, the elastic parameter sampled and forward modeling synthesis mechanism are carried out. Then, a comparison between the synthetic record and the actual one is performed to select the best matching elastic parameters. Adaptive spatial sampling is used to accelerate the iterative inversion process. The test is completed by using the theoretical data of Stanford University. The inversion results have high similarity with the theoretical model, and the absolute error is low. The simulation points with the relative error of impedance less than 10% account for 94%. The results of 13 blind wells show that the relative error less than 10% is accounting for 96.2%. It shows that the seismic inversion result of the new design method is reliable.

Based on the theoretical model test, the influence of complex factors such as noise, the relationship between statistical rock physics parameters and reservoir attributes, time depth conversion, and the scale integration of well and seismic data are excluded. In practical subsurface reservoir seismic prediction, these factors have an important impact on the inversion quality. Although this novel method has achieved success in theoretical model, how to carry out the actual reservoir inversion remains to be further tested and studied.

All the data can be accessible from the corresponding author.

We confirm that we have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

We are thankful for the financial support from the National Natural Science Foundation of China (grant nos. 42130813 and 41872138) and the Open Foundation of Top Disciplines in Yangtze University (grant no. 2019KFJJ0818029).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).