A High Resolution Method for Fluid Prediction Based on Geostatistical Inversion

In order to predict the fluid in thin layer precisely, this paper proposed a high-resolution method for fluid prediction. The method used geostatistical inversion with lithology masks to calculate water saturation. We applied this method to theoretical model and real data.Theresult was compared with that of prestack AVA simultaneous inversion for fluid prediction. It showed that this method had high resolution both in vertical and lateral directions for fluid prediction and could also predict the fluid in thin layer efficiently.


Introduction
Nowadays the main target of geological exploration has changed from a simple constructed oil and gas reservoir to complex structure reservoir and lithology reservoir.The need of precise exploration is increasing.The precise prediction for thin-layer fluid becomes an urgent problem for oil and gas exploration.Conventional fluid prediction is mainly based on AVO technology, including AVO analysis [1,2], prestack AVA simultaneous inversion [3,4], elastic impedance inversion [5,6], and other methods.However, these methods do analysis mainly based on seismic data, and their vertical resolution is limited.So it cannot precisely predict the fluid in thin layer and cannot meet the needs of precise exploration.In this paper, in order to solve this issue, we presented a high-resolution method for fluid prediction.The method used geostatistical inversion with lithology masks to calculate water saturation.The result was compared with that of prestack AVA simultaneous inversion.

Geostatistical Inversion with Lithology Masks.
Geostatistical inversion is based on geological information (including seismic, drilling, and logging).It applies random function theory and geostatistics (histogram analysis, variogram analysis, etc.) and combines with traditional seismic inversion technique.It can generate multiple optional inversion results with the same probability [7][8][9].It includes stochastic simulation and stochastic inversion.Stochastic simulation includes sequential Gaussian simulation, sequential indicator simulation, and other simulation methods.The geostatistical inversion with lithology masks mentioned in this paper contained two parts: sequential Gaussian simulation with lithology masks and stochastic inversion.

Sequential Gaussian Simulation with Lithology Masks.
The basic principle of sequential Gaussian simulation is the application of Kriging algorithm for locally estimating the geological variables, which is using a set of values at grid points where it is known (it is generally the value of well logging) to estimate the value at any grid point where it is not known (Figure 1).The key point of Kriging algorithm is to find the weight of known points to unknown points.For example, the formula used by the simple Kriging algorithm to estimate any unknown point using known points is where  * (  ) is the value of any unknown point, (  ) is the value of known points around,   is the weight of known points to the unknown point,  is the number of known points Figure 1: Sketch of Kriging algorithm (the black dots are known points and the white dot is the unknown point needed to be estimated).
around, and  is the mean of all the known data.The formula of calculating the weight is where  is covariance function, it can be obtained by variogram.
Variogram is the most commonly used method to measure reservoir space during geostatistics methods.The mathematical expression of variogram is where ℎ is lag distance, which is the distance between two points;  is measured value; (ℎ) is the variogram value with lag distance ℎ; (ℎ) is the number of data pair with distance ℎ.The relation between variogram and covariance is where (ℎ) is the covariance value with lag distance ℎ.
The basic steps of sequential Gaussian simulation are [10][11][12] as follows: (1) sequential Gaussian simulation requires the data with normal distribution, but most of known data distribute nonnormally in the application.So we need to do histogram analysis for all known data, then we do normal transformation for the data and calculate the mean  and variance  2 ; (2) according to some random order, do the following simulation to all the unknown grid points needed to be simulated; (3) at each unknown grid point, firstly we calculate the covariance function  according to formula (4) after analyzing the variogram; secondly, we use the Kriging algorithm to calculate the weight  of the known points around to the unknown point according to formula (2); thirdly, we obtain the mean   and the variance  2  of the unknown grid point according to formula (5); finally, a local probability density function (pdf) is built, and it is converted to cumulative distribution function (cdf): (4) draw the simulated value randomly from the cdf and add the simulated value to the set of the known points; (5) back to Step 3 until all unknown grid points are simulated; (6) do normal inverse transformation, transform the data to the original space, and finish one simulation; (7) redo another realization according to another random order.
Sequential Gaussian simulation can generate many data volumes with the same probability, which can meet the known conditions.We call them realizations.
Lithology mask is a method that according to lithology data volume, the lithology which is nonreservoir is masked, and sequential Gaussian simulation is done only in reservoir not in nonreservoir.The lithology data volume can be obtained from sequential indicator simulation using lithology data from well logs.
Sequential indicator simulation is another stochastic simulation method.This method has similar theories and procedures to those of sequential Gaussian simulation.Their differences are that sequential Gaussian simulation uses the simple Kriging method for continuous data, and sequential indicator simulation uses the indicator Kriging to estimate the local probability distribution for discrete data, such as lithology.The following are the steps for sequential indicator simulation [13,14]: (1) do indicator transformation to data.It assumes that (  ) is a group of discrete measured data at the points   ( = 1, 2, . . ., ), and there are  types of lithologies in research area.K types of indicator variables can be defined as where (  ) is the indicator variable at the point    is the lithology, and  = 1, 2, . . ., .The expectation value of indicator variable for lithology  at the point where {  (  ;  | (  ),  = 1, . . ., )} is the probability distribution of lithology  at the point   ; (2) according to some random order, do the following simulation to all the unknown grid points needed to be simulated; (3) at each unknown grid point, do variogram analysis and indicator Kriging calculation to every lithology, and calculate the weights of them to get their expectation values.Then, we can obtain the probability distributions of them at the grid point according to formula (7), and cumulative distribution function can be built; (4) the following steps are the same as those of sequential Gaussian simulation.
We cannot produce a continuous probability density function (pdf) for lithology; so we can only estimate the probability of each lithology and use them to build a cumulative probability distribution function (cdf).We need to notice that the sum of the probability values of all lithologies should be 100%.
Lithology mask is done after lithology data volume is obtained by doing sequential indicator simulation to lithology data from well logs.Then, sequential Gaussian simulation is done for water saturation data from well logs only in reservoir according to the steps of sequential Gaussian simulation.This method studies the fluid properties in reservoir directly and can avoid the case that non-reservoir with low water saturation is mistaken as containing oil and gas.It can make the simulation result more accurate.

Stochastic Inversion.
Reflection coefficients from the stochastic simulation results and wavelet were used to generate synthetic seismograms then, compare the synthetic seismograms with the original seismic data.If the result cannot meet the requirements for accuracy, we should resimulate the grid points of the channel until getting a good match between synthetic seismograms and original seismic data.Choose the value of the grid point having the best match with synthetic seismograms as the result of inversion.After all unknown points in a realization are completed, the inversion of another realization will go on [15][16][17].The inversion result is many data volumes with the same probability.The number of them is the same as that of stochastic simulation.During optimization, the most intuitive and effective way is to fully understand the geological background of the study area.Through comparing every data volume, we chose one result matching the geological knowledge or used the mean or weighted average of multiple realizations as the final inversion result.
Finally, according to a variety of water saturation threshold, we can distinguish between oil, gas, and water and predict the fluids.

Prestack AVA Simultaneous Inversion.
Prestack AVA simultaneous inversion is based on Aki-Richards formula and the idea of constrained sparse spike inversion.It carries out the simultaneous inversion constrained by well logs for the seismic data volume of multiple angle gathers [18,19].The low-frequency information of P impedance, S impedance, and density from well logs is a constrained condition.Calculate a minimization problem to the target function [20]: where

Theoretical Modeling
To illustrate the resolution and the effect of above methods, we built a multilayer model interbedded with a shale wall (Table 1).We did geostatistical inversion with lithology masks for the theoretical model and compared the result with that of prestack AVA simultaneous inversion.
The top depth, the total thickness, and the total width of the model were 1600 m, 144 m and 2200 m, respectively.It included shale, gas bearing sand, and water bearing sand with thickness of 30 m, 10 m, and 3 m.There was a shale wall with thickness of 144 m, and width of 640 m in the middle of the model.The shale wall was located between 780 m and 1420 m in lateral direction and penetrated the whole model in vertical direction.The shale wall was homogeneous medium.The parameters of each layer on both sides of the shale wall were the same.The distributions and the detailed parameters of the shale layers, the sand layers, and the shale wall were shown in Table 1.It was assumed that the water saturation of shale, water bearing sand, and gas bearing sand was 1, 0.8, and 0.2 respectively.The threshold value of gas was that the   /  being less than 1.712 or the water saturation being less than 0.55.Two wells were built on each side of the shale wall, well01 and well02, respectively.
Poststack seismic data and shot gathers were obtained by wave equation forward modeling for the theoretical model.The poststack seismic data was used for geostatistical inversion.The shot gathers were processed to get three angle gathers (near, middle, and far angle gathers), and they were used for prestack AVA simultaneous inversion.
In the result sections of theoretical modeling, the well logs of both wells were logging interpretation results.They were totally set according to the theoretical model in Table 1.The red was gas layer.The blue was water layer.The brown was shale.The white represented that there was no interpretation result in this layer (Figure 2).The parameters in both sides of the shale wall were the same, so the logging interpretation results of two wells were also the same.
Figure 3 shows the   /  profile of prestack AVA simultaneous inversion.In this figure, the value of   /  in golden yellow regions and red regions is less than 1.712, corresponding to gas layer.It can be seen from the figure that the gas layer with thickness of 30 m can be distinguished very clearly in vertical direction.It corresponded with well logs very well.The gas layer with thickness of 10 m can also be distinguished and corresponded with well logs, but the thickness is different from well logs.When the thickness of gas layer is 3 m, it is completely impossible to be distinguished.The thickness of the three gas layers in the inversion profile are larger than 10 m, which cannot correspond with well logs.The thin gas layer cannot be distinguished.In lateral direction, each thin layer was cut off by the shale wall except the 30 m layer.The shale wall is very clear and the boundary is also consistent with the model.Figure 4 shows the water saturation profile of geostatistical inversion with lithology masks.The water saturation of the yellow zone is less than 0.55, which is gas layer.The blue area is water layer.It can be seen from the figure that the shale on both sides of the shale wall has been masked in vertical direction.The gas layers with thickness of 30 m, 10 m, and 3 m can be clearly resolved and have a good agreement with well logs.The thin gas layer can be distinguished.In lateral direction, the gas layers with thickness of 10 m and 3 m and water layer are completely disconnected, and only the gas layer with 30 m thickness and water layer in the shale wall are not completely disconnected.The boundary of the shale wall is very clear and the prediction result for gas layer is close to that of prestack AVA simultaneous inversion.Geostatistical inversion with lithology masks can get a better prediction result for gas layer both in vertical direction and in lateral direction.

Analysis
From the results of theoretical modeling, it can be clearly seen that in lateral direction prestack AVA simultaneous inversion can reflect the change of formation, but it cannot distinguish thin gas layer in vertical direction.This is because prestack AVA simultaneous inversion uses the thought of constrained sparse spike inversion, which is based on direct conversion on seismic data.Its result retains most original features of seismic data.In lateral direction, seismic data has high density, so it has high lateral resolution.But it was constrained by the International Journal of Geophysics   seismic bandwidths and its basic assumption is that the strong reflection coefficient is sparse, only considering the impedance of strong reflections; so the details cannot be reflected and its vertical resolution is limited [21].However, the geostatistical inversion with lithology masks not only has high vertical resolution, but also has high lateral resolution.In vertical direction, thin gas layer can be well resolved.In the lateral direction, it can accurately reflect the lateral variation of the formation and gets the similar result with that of prestack AVA simultaneous inversion.This is because geostatistical inversion with lithology masks needs to go through sequential indicator simulation, sequential Gaussian simulation, and stochastic inversion.Sequential indicator simulation and sequential Gaussian simulation both take advantage of well data to locally estimate the geological variables based on Kriging algorithm.Their results also retain the basic characteristics of well data: high vertical resolution and low lateral resolution [22].The vertical resolution can theoretically reach the accuracy of well logs.Stochastic inversion generates synthetic seismograms with the result of stochastic simulation.Through the iteration contrast with original seismic data, it can get the final inversion result, which is mostly consistent with the original   seismic data.This process will make the results of stochastic simulation have the characteristics of seismic data.Although the vertical resolution will drop slightly, the lateral resolution has been greatly improved.Geostatistical inversion combined stochastic simulation with stochastic inversion, which can bring their advantage into play sufficiently.

Application
Here, we carried out geostatistical inversion with lithology masks on a two-dimensional line of some gas field.The result was also compared with that of prestack AVA simultaneous inversion.The threshold value of gas was that the   /  being less than 1.65 or the water saturation being less than 0.5.In the sections of the inversion results for real data, the well log of Well01 was still logging interpretation result.The red was gas layer, and the white represented that there was no interpretation result in this layer.Other types of layers were not interpreted in this well (Figure 5). Figure 6 shows the   /  profile of prestack AVA simultaneous inversion.In this figure, the value of   /  in golden yellow regions is less than 1.65, corresponding to gas layer.It can be seen from the figure that the vertical distribution of gas layer is consistent with the well logs basically, and only the thickness is different from well logs.In the inversion profile, the thickness of the continuous gas layer on the top of Well01 is 15 m, and the thickness of the gas layer below is 10 m.Their thickness are significantly higher than those of the gas layers shown in the well logs.Thin gas layer has not been accurately predicted.Figure 7 shows the water saturation profile of geostatistical inversion with lithology masks.The shale in nonreservoir has been masked.The water saturation of the yellow zone is less than 0.5, corresponding to that of gas layer.The blue area is water layer.It can be seen from the figure that the vertical distribution of gas layer is in good agreement with the well logs.The thickness is only a little larger than that of the interpretation result.In the profile, the thickness of the gas layer on the top of Well01 is 12 m, and the thickness of the thin gas layer below is only 4 m.Compared to the result of prestack AVA simultaneous inversion, this result is very close to that of the well logs.Thin gas layer has been accurately predicted.In the lateral direction, both two sections seemed more natural.The gas layer is neither too much nor too continuous.There are obvious changes in the lateral direction.By comparing these two pictures, it can be clearly seen that calculating water saturation by geostatistical inversion with lithology is a fluid prediction method which has high vertical and high lateral resolution.

Conclusion
Geostatistical inversion utilized well data as simulation control points and was constrained by seismic data between the wells.The inversion result can reflect the character of well data and seismic data at the same time.It gave full play to the advantage of high vertical resolution of well data and high lateral resolution of seismic data.Using geostatistical inversion with lithology masks for water saturation can directly research the properties of fluids in reservoir.It can avoid the case that non-reservoir with low water saturation is mistaken as containing oil and gas.It can reduce the error of the inversion result.It can get the fluid prediction result with high vertical and lateral resolution and precisely predict thin layer fluid.From the theoretical modeling and real data application, it can be found that this high-resolution method of fluid prediction gets a good prediction result for thin layer fluid.

Figure 2 :
Figure 2: The logging interpretation results of well01 and well02 in theoretical model.

Figure 3 :
Figure 3: The   /  profile of prestack AVA simultaneous inversion for the theoretical model (30 m, 10 m, and 3 m are the thickness of the gas layers).

Figure 4 :
Figure 4: The water saturation profile of geostatistical inversion with lithology masks for the theoretical model (30 m, 10 m, and 3 m are the thickness of the gas layers).

Figure 5 :
Figure 5: The logging interpretation result of Well01 in real data.

Figure 7 :
Figure 7: The water saturation profile of geostatistical inversion with lithology masks for the real data (12 m and 4 m are the thickness of the gas layers).
Through iterations, the stable data volumes of P impedance   , S impedance   , and density  can be inversed simultaneously.The parameters which can reflect fluid properties, such as the ratio of P velocity to S velocity   /  , Poisson's ratio , bulk modulus , shear modulus , and lame coefficient , can be calculated by the relationships among all kinds of the parameters.

Table 1 :
The parameters of multi-layered model with a shale wall.
Figure 6: The   /  profile of prestack AVA simultaneous inversion for the real data (15 m and 10 m are the thickness of the gas layers).