A new prediction method of reservoir porosity based on improved Kriging interpolation

In the process of oil and gas field development, the pore characteristics of reservoir are the main factors that affect the reservoir capacity and production efficiency. The physical property parameters and pressure distribution of layers have changed greatly, and the pressure system of multi-pressure layer is formed underground of multi-layer reservoir and heterogeneous sandstone oil reservoir after water injection and chemical injection. Accurate prediction of reservoir physical properties and formation pressure distribution after changes is significance for efficient production and safe operation of reservoirs. In this paper, based on geostatistics theory and percolation mechanics theory, combined with reservoir geology, drilling, development and testing data of Bohai Oilfield, the prediction method of reservoir porosity based on improved Kriging interpolation and effective thickness has been established. This prediction method can provide certain guidance for the formulation of oilfield injection and production distribution plans at the oilfield and the shut-in and adjustment plans for injection-production wells during drilling and completion operations.


Introduction
In the process of oil and gas field development, the main factor that affects the storage capacity and production efficiency of reservoir fluids is the pore structure of the reservoir rock. After the development of multi-oil layer and heterogeneous sandstone oil fields with water injection and chemical injection [1,2], the physical parameters and pressure distribution of the formation have changed greatly. Taking Bohai Oilfield as an example, the development of fault blocks leads to formation stress concentration, the risk of borehole collapse during drilling operations is relatively high, and there are many complicated downhole conditions. Moreover, the strata in each block are very different and affected by ground stress. The multi-layer pressure system brings the difficulties of easy to leak, easy to spray, and easy to jam the broach. Therefore, the prediction of reservoir physical properties and formation pressure is directly related to the success or failure of drilling operations.
The physical properties of the adjustment well reservoir mainly refer to the thickness, porosity, permeability, and oil and gas content of the adjustment well formation. From the perspective of reservoir engineering, porosity is of great significance to the exploration and development of oil and gas fields, and it is also one of the main contents of reservoir research. After years of exploration and development, there have been many methods to obtain porosity from seismic data. The main methods can be roughly  [3][4][5][6][7][8][9][10]: Wyllie formula method, well-constrained inversion and non-wellconstrained inversion, function approximation method and Kriging and its improvement method. The Wyllie formula method is simple in principle and easy to use. It is one of the methods widely used at present, but its calculation accuracy is not high. The accuracy of well-constrained inversion method is higher than that of non-well-constrained inversion method, which is suitable for blocks with sparse well points. The prerequisites of these two methods are difficult to meet, so the application effect is not good. The function approximation method establishes the functional relationship between porosity and seismic characteristics based on the fitting of seismic data and logging data beside the well, and uses the function approximation function to predict porosity. This method is only suitable for areas with uniform porosity distribution, and not suitable for areas without wells or wells but extremely uneven porosity distribution.
Kriging method is a method of spatial interpolation of irregularly distributed log porosity by using correlation function. It takes into account the spatial position relationship between porosity values. It is currently one of the most accurate and widely used methods. Kriging interpolation is the core of geostatistics [11,12]. It mainly uses covariance function and variogram to study regional variables, including reservoir physical parameters such as thickness, porosity, and permeability obtained from well profiles. The randomness and regionality of this kind of space have correlation. In this paper, we first take the value of the well points around the point to be sought, obtain the local variogram, then calculate the covariance function, and finally perform the improved Kriging interpolation method of ordinary Kriging interpolation estimation. This method can well characterize the geological conditions within a small range of the points to be sought, and the calculation accuracy is high.

Ordinary Kriging interpolation principle
Kriging interpolation method is a spatial prediction method proposed by South African geomineralist P G. Krige. It is based on the theoretical analysis of variogram, and is an unbiased and optimal method for estimating the value of regional variables in a limited area. Different from the traditional interpolation method, this method not only considers the spatial position between the adjacent observation data point and the point to be interpolated when estimating the sample value, but also considers the positional relationship between the neighboring points near each estimated sample value [13]. Therefore, the estimated value is more accurate and more realistic than the traditional method.
Ordinary Kriging interpolation is the most commonly used unbiased optimal estimation method in Kriging interpolation. It considers the spatial relationship between the point to be measured and the known point or between known points. In view of the geological and well pattern conditions of the Bohai Oilfield, in terms of inter-well valuation, the Kriging interpolation method has relatively high valuation accuracy and can reflect the objective laws of geology. It is suitable for the case where the mathematical expectation of the estimated regionalization variable is constant. The calculation formula is similar to that of the inverse distance weighted interpolation method: Among them, ( ) (i = 0,1,2, . . . , n) is the characteristic value of the data point , 0 is the unknown point, the other is the known point, and is the weight coefficient. Different from the inverse distance weighted interpolation, the weight coefficient here is not only determined by the distance, but also judged according to the calculation result of the variogram. The conditions for calculating the variogram are under the unbiased condition and the minimum variance. This algorithm embodies the advantages of ordinary Kriging interpolation.
Assume that * ( 0 ) still adopts the form in equation (1). Different from the simple Kriging interpolation method, in order to ensure that * ( 0 ) is an unbiased estimate of * ( 0 ), it is assumed that 0 = 0, which can be derived 1 − ∑ = 0 −1 , so that the unbiased linear estimation of ( 0 ) can be written as follows: If 0 = 1, = − (i = 1,2, . . . , n), then the same as the derivation process of simple Kriging interpolation equations, the estimated variance has the following expression: In order to make the estimated variance of (i = 1,2, . . . , n) obtained under the constraints of ∑ = 1 by minimizing variance of estimation errors, the Lagrangian multiplier method is used to open the constructor ) . If a set of values (i = 1,2, . . . , n) meets the above requirements, then this set of values can satisfy the following formula: The above is the ordinary Kriging interpolation equations. According to formula (4), the ordinary Kriging interpolation equations in matrix form can be obtained as: Among them: In the calculation of ordinary Kriging interpolation estimation, the estimated variance obtained is the minimum estimated variance. This minimum estimated variance is called ordinary Kriging interpolation variance and can be expressed by the following formula: In the ordinary Kriging interpolation method, the obtained weighting coefficient is called the Kriging interpolation weighting coefficient. It can be seen from the form of the Kriging interpolation equations that these Kriging interpolation weighting coefficients not only depend on the functional relationship between , and the covariance function ( , ) of the random function ( ), but also it depends on the relative positional relationship between the known points , and the relative positional relationship between the estimated point 0 and each known point .

Variogram analysis
Variogram is one of the basic tools of geostatistics, which can describe the changes of regionalized variables in the spatial structure, and can describe the random changes of regional changes. Here, the variation function is used to find the covariance function in the Kriging interpolation method [14], so as to solve the Kriging interpolation equations.
Suppose there is a stationary random function whose mean value m and variance 2 are known. The average value and variance have nothing to do with their location, that is to say, m( ) = and 2 ( ) = 2 in the entire study area. The average value m usually includes regional and vertical trends, which are processed by the deterministic model of the average value and work together with the residuals of the local variable averages. The variogram is defined as: The semivariogram γ(h) is half of the variogram 2γ(h). The variogram is a measure of variation, and the value of the variogram increases when the sample point gap becomes larger. Covariance is a numerical quantity used to measure correlation: By definition, the covariance C(0) at ℎ = 0 is the variance 2 . When the scattered ℎ values are nonlinearly related, the value of C(0) is zero.
Expanding the square in formula (6), the following relationship between semivariogram and covariance is obtained: The prerequisite for the establishment of this relationship is that the mean and variance of the model are both constant and independent of the position. That is: (1) the base value of the variogram is the variance, and the variogram at this time is zero-correlated; (2) Where 0 is the nugget effect constant, referred to as the nugget constant, c is the arch height, 0 + is the abutment value, and is the range. Because when h = 0, γ(h) = 0; h > a, γ(h) = 0 + , so we only need to consider the case of 0 < h < a. So there is: Set 0 = 0 , 1 = 3 2 , 2 = − 1 2 3 , = (h), 1 = ℎ, 2 = ℎ 3 . Get the linear function from (10): Suppose the value of the experimental data ℎ is ℎ 1 , ℎ 2 , ⋯ , ℎ , and the experimental variation function is calculated accordingly as: * (h 1 ), * (h 2 ), * (h 3

), ⋯ , * (h )
And substituting them into equation (14), the linear equations is: Where, 1 = ℎ , 2 = ℎ 3 , = * (h ), = 1,2, ⋯ , , is the coefficient matrix:  Analyzing the relationship between the geological variation and the observed behavior of the variogram [15,16] is the basis for accurate deduction and simulation of the variogram. Figure 1 shows the relationship between the vertical and lateral variograms under different geological conditions. Each graph is a corresponding variogram in the vertical and horizontal directions. In reality, the detailed graph of the variable is unknown, but the behavior of the variogram must be directly deduced from the graph of the variogram. The initial variogram behaviors include: irregular or lack of spatial connections, spatial relationships that decrease with increasing distance, geological structure direction, plane orientation, stratigraphic stratification, geological periodicity, etc. The actual variogram is usually a comprehensive reflection of the behavior of these different variograms.
There is uncertainty in the modeling of the variogram. One way to capture this uncertainty is to include geological information in the variogram model. The geological information to be added includes, but is not limited to, the geological environment, sequence stratigraphy, pore properties, isoporosity maps, and isopermeability maps. The interpretation of this information can derive the main continuous direction, horizontal expansion and anisotropy index. Geological information and reservoir description Due to the sufficient vertical data, the estimation of the vertical variogram is very easy. On the contrary, the estimation of the lateral variation function is a difficult problem [17,18]. Therefore, most of the calculation process is just pure nugget variation function. The pure nugget variogram is a variogram that has no correlation on all lags. Simulating this variation function will only get random distribution. Due to the spatial relationship of geological data, the pure nugget model is not suitable for geological models [5,19]. Due to the scarcity of data, pure nugget variograms in the horizontal direction are often unavoidable.
This paper uses the error averaging technique proposed by Bahar and Kelkar to solve this problem. The main reason for using this technique is that this technique is based on a continuous geological inference of contribution, such as considering only a low frequency information and filtering out the high frequency part of the data. In this technology, the sill value of the variogram is taken from the vertical variogram, and the range of the variogram is calculated after calculating the vertical mean of the data to generate 2D data, and calculate the variogram from these 2D data. Figure 2 is a comparison chart of pure nugget variogram and error average technical results.  As a measure of spatial variability, the variogram is used in geostatistical reservoir modeling. In a statistical sense, the 3D model is constructed based on the variogram. Therefore, whether the variogram can represent the heterogeneity of the entire study area is very important. It represents the modeler's quantitative understanding of the variability of the entire desired area and other relevant geological information.

Improved Kriging interpolation method to predict reservoir physical properties
In terms of well spacing and geological conditions in the Bohai Oilfield, the ordinary Kriging interpolation method can better reflect the objective laws of geology, and its estimation accuracy is also higher. In terms of selecting the variation function model, the spherical model is suitable for threedimensional conditions and is a commonly used method in geological estimation. Therefore, the  (1) Establish a spherical model to find the variation parameters 0 and Spherical model: Due to the complex geological conditions in the target area, it is impossible to fit a variogram model that is applicable to the entire block. Therefore, when calculating the variogram parameters, you can first use the nearby known point data to find γ(h) = ( − ( +ℎ) ) 2 /2 , where Z is the value of the known point. Substitute the γ(h) and h of each point taken into (ℎ) = 0 + ( 3 2 ) ℎ + (− 2 3 ) ℎ 3 , take the diameter of the circle in the selected point range for , find c and 0 greater than zero, and determine the variation function calculation formula.
According to the data of known well points A, B, C, and D, the c and 0 values (c=11.15, 0 =0.327) obtained by inverse calculation, a is set to 782.58. c=11.15, 0 =0.327, and =782.58 can be used to predict the physical properties of the D zone. Then the variation function is: (2) Ordinary Kriging interpolation method to solve the weight function is the covariance function, and (h) = (0) − (ℎ), ( , ) is the covariance between two different points. Here ( , ) is abbreviated as When = j , it is the same point itself, without variation function (ℎ) = 0 , and because of (h) = (0) − (ℎ), get ( , ) = (0) = 2 = 0 + = 11.48.

When
≠ j , = (| − |) = 0 + − (| − |) , where (ℎ) = 0 + ( 3 2 ) ℎ + (− 2 3 ) ℎ 3 , | − | is the distance between the two points to be calculated, thus Out all the values of ( , ). And bring it into the matrix [K][λ] = [ 2 ] to calculate the parameters and the values of D h , D k , D ϕ . By comparing with the actual value measured on the oil field, the calculation error is obtained. The calculation result is shown in Table 1

Conclusion
According to the characteristics of the Bohai Oilfield, a prediction method for reservoir physical properties predicting based on the ordinary Kriging interpolation method with the improved variogram is given. By using this method, the parameter of the unknown point can be calculated by the local variation function, which can first obtained by the reservoir properties of the well point around the unknown point. The thickness deviation rate of single well is less than 10%, porosity deviation rate is less than 5%, and permeability deviation rate is less than 15% of reservoir physical property prediction by using the improved Kriging interpolation method. This method can well represent the geological situation of the unknown point with a small deviation, and the calculation accuracy is higher. This prediction method can provide certain guidance for the formulation of oilfield injection and production distribution plans at the oilfield and the shut-in and adjustment plans for injection-production wells during drilling and completion operations.