- Email: [email protected]
Stochastic seismic inversion based on an improved local gradual deformation method
Accepted Manuscript Stochastic seismic inversion based on an improved local gradual deformation method Xiuwei Yang, Peimin Zhu PII:
S0098-3004(16)30506-4
DOI:
10.1016/j.cageo.2017.08.010
Reference:
CAGEO 4007
To appear in:
Computers and Geosciences
Received Date: 8 October 2016 Revised Date:
5 August 2017
Accepted Date: 8 August 2017
Please cite this article as: Yang, X., Zhu, P., Stochastic seismic inversion based on an improved local gradual deformation method, Computers and Geosciences (2017), doi: 10.1016/j.cageo.2017.08.010. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
ACCEPTED MANUSCRIPT 1
Stochastic seismic inversion based on an improved local gradual deformation method Xiuwei Yang and Peimin Zhu China University of Geosciences, Institute of Geophysics and Geomatics, Wuhan, China. E-mail:[email protected]; zhupm@ cug.edu.cn
7 8 9
ABSTRACT A new stochastic seismic inversion method based on the local gradual deformation method is proposed, which can incorporate seismic data, well data, geology and their spatial correlations into the inversion process. Geological information, such as sedimentary facies and structures, could provide significant a priori information to constrain an inversion and arrive at reasonable solutions. The local a priori conditional cumulative distributions at each node of model to be inverted are first established by indicator cokriging, which integrates well data as hard data and geological information as soft data. Probability field simulation is used to simulate different realizations consistent with the spatial correlations and local conditional cumulative distributions. The corresponding probability field is generated by the fast Fourier transform moving average method. Then, optimization is performed to match the seismic data via an improved local gradual deformation method. Two improved strategies are proposed to be suitable for seismic inversion. The first strategy is that we select and update local areas of bad fitting between synthetic seismic data and real seismic data. The second one is that we divide each seismic trace into several parts and obtain the optimal parameters for each part individually. The applications to a synthetic example and a real case study demonstrate that our approach can effectively find fine-scale acoustic impedance models and provide uncertainty estimations.
24
SC
14 15 16 17 18 19 20 21 22 23
M AN U
12 13
TE D
10 11
RI PT
2 3 4 5 6
25
Keywords: Stochastic inversion; Heterogeneities; Probability field simulation; Gradual deformation method.
29
1 INTRODUCTION
30 31
Seismic inversion plays an important role in reservoir prediction. The inverted impedance from seismic data can indicate changes in reservoir properties. However, an accurate impedance model is difficult to invert since seismic data is bandwidth limited. Both low frequencies and high frequencies outside of the seismic bandwidth cannot be recovered from seismic data alone. Additional constraints are used to reduce the non-uniqueness of the inversion problem by providing a priori information that is not contained in the seismic data. Geological information, such as sedimentary facies and structures, could provide significant a priori information to constrain the inversion and arrive at reasonable solutions (Gonzalez et al., 2008; Nunes et al., 2016). A conventional seismic inversion, such as generalized linear inversion (Cooke and Schneider, 1983) or sparse spike inversion (Oldenburg et al., 1983), usually gives a best estimate of the impedance that minimizes the difference between synthetic seismic traces and real seismic traces. Such techniques yield a single estimate, which should be the mean of all possible non-unique solutions and relatively smooth. Thus, the deterministic seismic inversion
34 35 36 37 38 39 40 41 42 43
AC C
32 33
EP
26 27 28
2
ACCEPTED MANUSCRIPT
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
RI PT
53 54
SC
51 52
M AN U
49 50
TE D
48
EP
46 47
is only suitable when the reservoir is comprised of relatively thick layers with slow lateral thickness changes and minimal spatial impedance variation (Francis, 2006). It is worth noting that this method loses the small-scale heterogeneities and overestimates the connectivity of the parameters. Geostatistical methods could simulate small-scale variabilities not captured in seismic data. Geostatistical inversion (Doyen, 2007 and Bosch et al., 2010) provides a better strategy for generating a suite of alternative subsurface property realizations conditioned on different kinds of data and incorporates spatial correlation relationships. It is more suitable for use in uncertainty estimations, volumetric calculations or connectivity estimations since it includes the small-scale heterogeneities and spatial variabilities of the geology in the inverted model. Geostatistical inversion using sequential Gaussian simulation was first tested by Bortoli et al. (1993) and Haas and Dubrule (1994). To be constrained to seismic traces, they simulate a whole trace at a time instead of the traditional node-by-node sequential Gaussian simulation, and use a simple Monte Carlo rejection scheme to obtain the optimal trace. Grijalba-Cuenca et al. (2000) offered another geostatistical inversion algorithm, which worked node by node and was optimized using the simulated annealing method. Mukerji et al. (2001) and Eidsvik et al. (2004) introduced statistical rock physics to estimate reservoir parameters from prestack seismic data and evaluated the associated uncertainty. In recent years, some new geostatistical techniques and inversion algorithms have been introduced in seismic inversion. A global iterative geostatistical seismic inversion method (Soares et al., 2007, Nunes et al. 2012, Azevedo et al. 2015), using direct sequential simulation and co-simulations (Soares, 2001), has been proposed. Several models for the entire seismic grid are simulated at once in the first step. Local areas of best fit of the different models are selected and “merged” into a secondary model for the direct co-simulation of the next iteration. Gonzalez et al. (2008) incorporated multiple-point geostatistics and rock physics to characterize previously known geologic information, simulate reservoir models and then select the optimal realization. Grana et al. (2012) estimated the facies and reservoir properties based on the probability perturbation method and tau model (Journel, 2002). Buland and Omre (2003) and Buland et al. (2003) developed a linearized AVA inversion in a Bayesian framework, both in the spatial and Fourier domains. They assume that the elastic parameters are characterized by a log-Gaussian random field and an explicit analytical form of the Gaussian posterior distribution can be efficiently obtained. Grana and Della Rossa (2010) combined this Bayesian linear AVA inversion with statistical rock physics to estimate reservoir properties and the associated uncertainty. The use of Bayesian framework made the link between stochastic inversions and deterministic inversions clearer. The fast Fourier transform moving average (FFT-MA) method (Le Ravalec et al., 2000) is the implementation of the moving average (MA) method (Oliver, 1995) in the frequency domain. The FFT-MA method, combined with the gradual deformation method (GDM) (Le Ravalec et al., 1999; Hu, 2000), has been successfully used in petroleum engineering. This approach uses a perturbation mechanism to gradually change the Gaussian-related reservoir model and guarantees the spatial correlation. Local perturbation (Le Ravalec et al., 1999; Hu, 2000; Marteau et al., 2014) and gradient information (Hu and Le Ravalec-Dupin, 2004) are used to improve the convergence rate. Another common practice to speed up convergence consists in combining several realizations instead of two, which implies the use of several
AC C
44 45
3
ACCEPTED MANUSCRIPT
113
2 METHODOLOGY
114 115 116
The methodology we proposed can integrate the geological a priori information and their spatial correlations in seismic inversion. The workflow of the method is shown in Figure 1. An indicator cokriging (Deutsch and Journel, 1998) is first applied to take into account well data and geological a priori information and to obtain the ccdf of the acoustic impedance (AI) at each grid node of the model to be inverted. Sampled from the locally defined ccdf, the p-field simulation will then generate different realizations consistent with same spatial correlation. A spatially correlated 3D p-field is easily simulated using the FFT-MA technique. Finally, an improved local gradual deformation method (GDM) is applied to perturb the Gaussian white noise locally until there is a match with the seismic data. Two strategies are proposed to make GDM suitable for seismic inversions. Each step of the method used in this research is described in the following subsections.
99 100 101 102 103 104 105 106
117 118 119 120 121 122 123 124 125
SC
97 98
M AN U
95 96
TE D
93 94
EP
92
AC C
90 91
RI PT
107 108 109 110 111 112
deformation parameters (Le Ravalec-Dupin and Noetinger, 2002). GDM has been successfully applied to several synthetic and real datasets (Le Gallo and Le Ravalec-Dupin, 2000; Le Ravalec-Dupin et al. 2004; Pourpak et al., 2009; Gervais-Couplet et al. 2010; Le Ravalec and Mouche, 2012) for history matching, which characterizes a reservoir by matching the static (e.g., porosity, permeability) and dynamic data (e.g., pressure, oil rate) from the wells. Data matching in seismic inversion is different from history matching, since seismic data is defined by millions of points and covers the whole study area. There are problems reducing the mismatch of amplitudes between the numerical seismic responses and the real seismic data exactly using original local GDM. Therefore, new strategies are necessary to develop. In this paper, geological information is included by indicator cokriging to obtain the local conditional cumulative distribution function (ccdf). In order to simulate fine-scale reservoir models consistent with the expected geology, well data and their spatial correlation efficiently, FFT-MA and probability field (p-field) simulation (Srivastava, 1992; Froidevaux, 1993) are used. Then, optimization is performed to match the seismic data via an improved local gradual deformation method. In order to match seismic amplitude at each grid-node, two improved strategies are proposed to make original GDM suitable for seismic inversions. The first strategy is that we select and update local areas of bad fitting between synthetic seismic data and real seismic data. The second one is that we divide each seismic trace into several parts and obtain the optimal parameters for each part individually. We can apply the methodology several times using different random numbers to generate different realizations. Then posterior distribution can be calculated using all these inverted realizations and give uncertainty estimation. We begin by briefly describing the workflow of this approach and then discuss each step in detail, especially the two improved strategies of GDM. The method is first tested on a 3D synthetic case and then on a real data.
88 89
4
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
126 127 128
Figure 1. Workflow of the stochastic inversion.
142
r be a series of cutoff values to discretize I (u ) ; i (uα ; r ) is the indicator transformation of
143
I (u ) ; α stands for the index of AI ( α = 1, 2, ... , n); and y (u; r ) is the a priori probability
144
conditioned on geological information. By solving this kriging system and getting the
145
coefficients λ0 (u , r ) , λα (u; r ) , and v(u; r ) , attributed to F (r ) , i(uα ; r ) and y (u; r ) ,
146
respectively, the a posteriori probability of I (u ) can be obtained using:
137
EP
135 136
AC C
132 133 134
TE D
138 139 140 141
2.1 Local conditional cumulative distribution Generally, the probability function of a continuous variable is often assumed to be Gaussian in modeling or inversion methods. However, the distribution of model parameters in geology is often not close to being Gaussian. Sometimes, the distribution even cannot be expressed using a set of parameters and an explicit analytical form. In this case, the indicator kriging would be selected for a so-called non-parametric estimation of probability distribution of variables. It has the advantage that any distribution can be handled. We first discretize continuous variable AI into K classes by (K-1) cutoff values. Then, a linear interpolation is performed to recover the within-class resolution lost within each separate subclass. K should be a reasonable number that could guarantee each class includes adequate samples from well data. And meanwhile, the inferred ccdf can characterize the main statistical features of AI. In this paper, indicator cokriging using the Markov-Bayes model is used to obtain the local ccdf, conditional on geology and well data. Let I (u ) be the AI value at location u and
129 130 131
n
147
Prob( I (u ) < r | (n + 1)) = λ0 (u , r ) F (r ) + ∑ λα (u; r )i (uα ; r ) + v(u; r ) y (u; r ) α =1
(1)
5
ACCEPTED MANUSCRIPT 148 149
where (n + 1) represents the number of conditioning data including n neighboring AI data and 1 co-located a priori probability. F (r ) is the global a priori probability common to all
150
locations u . These steps only need to be carried out one time before the simulation and inversion; thus, this method is sufficiently simple and fast to handle a large 3D case.
151 152 153 154
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
RI PT
SC
162 163 164 165 166 167
M AN U
160 161
TE D
159
EP
157 158
2.2 FFT-MA and p-field simulation Probability field (p-field) simulation (Srivastava, 1992; Froidevaux, 1993) is a fast conditional simulation method. Contrary to the sequential methods, which need to reconstruct ccdf for each realization at every location, the ccdf of the p-field simulation could remain the same from one realization to another, only conditional on the original data. After we obtained this local ccdfs, the conditional simulations are then obtained by drawing samples from these ccdfs. The probability values used to draw from the local ccdfs constitute a probability field, and are viewed as outcomes of a random function characterized by a uniform distribution and a given covariance function. In this paper, this correlated p-field is quickly constructed by the fast Fourier transform moving average (FFT-MA) method (Le Ravalec et al., 2000). The p-field realizations need not be conditional since the ccdfs are already conditioned on the original data. Since the model generated by FFT-MA follows a Gaussian distribution with a mean of 0, unit variance and autocorrelation ( C ), we have to transform it to a uniform distribution between 0 and 1 in order to ensure the reservoir model satisfies the ccdfs we have obtained. The use of FFT-MA and p-field simulation allows us to generate realizations satisfied by any specified ccdf and spatial correlation. The simulations are fast and stable because construction of the p-field is carried out through FFT, and the ccdf is produced only once, conditional on the original data. Another important feature is that the FFT-MA generator draws the randomness in the spatial domain, which allows us to uncouple the random numbers from the structural parameters. Such features imply that any of these random numbers can be perturbed locally without affecting the others outside of the range of the correlation length. It is potentially of great interest for seismic inversion to modify the AI model locally, based on a comparison between the actual seismic data and synthetic seismic responses of the inverted AI.
AC C
155 156
2.3 Improved local gradual deformation method The model constructed by the p-field simulation has been conditioned on well log data but not conditioned on seismic data. Therefore, we need an inversion to make the result match the seismic data. Generally, a seismic impedance inversion is non-unique. The classical gradient seismic inversion method gives only one smooth result and loses small-scale heterogeneity. Thus, the stochastic inversion is more suitable to obtain different realizations and estimate the uncertainty. The gradual deformation method (GDM) has turned out to be very consistent with FFT-MA and has been successfully used in history matching. Since seismic data is different from production data, two improvements are proposed to make the method suitable for seismic inversion. In geophysics inversion, least-absolute-values criterion (L1 norm) and least-squares
6
ACCEPTED MANUSCRIPT
criterion (L2 norm) are two of the most widely used criterions. The choice of the norm should be driven by the type of data and the purpose of the application. When outliers are suspected in a data set, long-tailed probability density functions should be used to model uncertainties. In this case, L1 norm is suitable since its result is robust to outliers. If uncertainties in the problem can be modeled using Gaussian distributions, L2 norm is a more appropriate choice (Tarantola, 2005). In most case of seismic inversion, the errors of observed data are proved to be a Gaussian distribution or generalized Gaussian distribution. Therefore, we first define the objective function J as the sum of the squared errors of the seismic responses S m and
200
observed data S mobs . J=
201
204 205
1 M ∑ (Sm − Smobs )2 2 m =1
(2)
where M represents the number of sample points of the seismic trace. The lengths of the observed and synthetic seismic data are the same since the inversion is conducted in time domain without the conversion between time and depth. The synthetic seismic trace S is computed by convolving the wavelet w with the reflection coefficients series r :
S = r∗w
SC
202 203
RI PT
192 193 194 195 196 197 198 199
(3) The core idea of GDM is that a new set of independent Gaussian white noise z can be
208
created using two sets of independent Gaussian white noise, denoted z1 and z 2 ,
209
respectively. z1 is the Gaussian white noise for the initial realization and z 2 is another
210
randomly generated one. In addition, their gradual deformation relation is defined as follows:
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232
TE D
216 217
(4)
Then, a search process is implemented to determine the deformation parameters r minimizing Eq. (2). At this point, the fit may not be good enough. A new starting realization chain is constructed from the optimal realization identified at the end of the previous search process. In this way, we only need to create a sequence of independent Gaussian white noise, generate realizations via FFT-MA and p-field simulation, and estimate the optimal r to further minimize the objective function. This process is iterated until a satisfactory matching is achieved. The global GDM does not have a high efficiency because improving the fit in one region may deteriorate the fit in another region. In history matching, a local gradual deformation is developed by partitioning the field into several independent regions and performing the inversion individually for each region. Deterministic partitioning is available in history matching because the observed data (e.g., pressure, oil rate) are scattered in different zones of the field. However, the seismic data are different from the production data because they are defined by millions of points and cover the whole area, aside from some sparse locations. In order to match the seismic amplitude at every grid-node exactly, we developed two improved strategies to make GDM suitable for seismic inversion. The first strategy is to select local areas with bad fits to modify. At the beginning of the inversion, we can modify all of the nodes, and the objective function will decrease quickly. When the convergence slows down, modifying all the nodes at one time is not a good choice. In some regions, if we have obtained a good fit, there is no need to modify the impedance model. Therefore, we only need to update the regions with poor fits. The squared residual at
EP
212 213 214 215
z = z1 cos(π r ) + z 2 sin(π r ), r ∈[-1,1]
AC C
211
M AN U
206 207
7
ACCEPTED MANUSCRIPT
240
z = (1 − H).× z1 + H. × [z1. × cos(π r ) + z 2 . × sin(π r )], r ∈[-1,1]
241 242 243 244 245
where the symbol .× stands for the componentwise product. H is the binary mask. When the convergence slows down, it is difficult to obtain a better model to match the seismic data using this threshold. At this time, we can increase the threshold value T, and the number of nodes that need to be modified is decreased. The modification will focus on the nodes that have the worse fitting and the model will be further refined. The same optimized process is repeated again and again. The modified region decreases to consider the local misfit of the data.
250 251 252 253 254 255 256 257 258 259 260 261 262 263
Figure 2. Binary mask (bottom) defining the modified region obtained from the squared residual image (top) and a threshold.
AC C
249
EP
TE D
246 247 248
(5)
SC
237
M AN U
235 236
RI PT
238 239
each node between the real seismic data and synthetic seismic data stands for the mismatch. Therefore, it can be used to determine the nodes which need to be modified. Thus, at each iteration, we first calculate the squared residual at each node and define a threshold T. When the squared residual of one node is larger than T, this node needs to be modified. We transform the squared residual image into a binary mask using the threshold T, and the modified region is associated with the white pixels (Figure 2). At this time, the GDM equation can be expressed as
233 234
The second strategy is to divide each trace into several parts and to obtain the optimal parameters for each part individually. Even though we have selected local areas with bad fitting to modify, one deformation parameter for all the nodes is still not suitable. A local GDM with multiple parameters is a better choice, but a reasonable partitioning is difficult. We propose a strategy to divide the model and obtain the optimal parameters. In the lateral direction, we assume each trace is independent, since post-stack seismic data has eliminated the trace dependence. As such, each trace is a single subregion with one deformation parameter. In the vertical direction, the amplitude at one node is only influenced by the impedance within the range of the wavelet length. Therefore, we simply divide each trace into several parts. The spatial continuity in the lateral direction and between the different parts in the vertical direction can still be guaranteed with the use of FFT-MA. Therefore, the number
8
ACCEPTED MANUSCRIPT
268 269 270 271
RI PT
266 267
of deformation parameters is n = Nx × Ny × Mz, where Nx represents the number of inline, Ny represents the number of crossline, and Mz is the number of parts in the vertical direction. Then, we solve an optimization problem of the n parameters, t1, t2, …, tn. We choose the golden section search method for all parameters simultaneously and independently based on the local objective function. In this way, the computing time is approximately equivalent to that of one parameter. In fact, when the optimization is done, the optimal parameters are able to divide the whole model into several subregions automatically (Figure 3). The number and shape of the subregions is flexible and different each time.
SC
264 265
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
TE D
281 282 283
Therefore, the optimization of our proposed method is performed within two nested loops. In the outer loop, the whole field is divided into n parts. In each part, the modified nodes are defined based on the squared residuals and threshold T. In the inner loop, we perform an optimization on the n parameters of the gradual deformation method to minimize the objective function. If the objective function of the new model is less than that of the previous model, the new model is accepted. Otherwise, we generate a new Gaussian white noise and repeat the previously described steps. We iterate over this procedure until the objective function is less than a tolerance value.
EP
275 276 277 278 279 280
Figure 3. Optimal parameters of the local GDM (left) are equivalent to the partitions and its parameters in the image to the right. Each rectangle stands for one part, including several samples of one trace in the vertical direction.
2.4 Parallelization In recent years, many parallel computational strategies for sequential simulation have been developed at three distinct levels: the realization level, path level and node level. Multi-core Central Processing Units (CPUs) (Vargas et al., 2007, Nunes and Almeida, 2010) and graphics processing units (GPUs) (Tahmasebi, 2012, Ferreirinha, 2015) are both exploited to improve the calculation efficiency. However, their sequential nature makes them difficult to parallelize very efficiently. The simulation of N nodes cannot be simply parallelized because there are dependences between the sequential and parallel results. The parallelization of our proposed approach is easier than the inversion using a sequential simulation. We exploit the compute unified device architecture (CUDA) to increase the computational efficiency. CUDA (Sanders and Kandrot, 2010) is a parallel computing architecture developed by the technology company Nvidia (2008) that makes it possible to benefit from the ability of GPU in parallel computing with general-purpose computation. In
AC C
273 274
M AN U
272
9
ACCEPTED MANUSCRIPT
302 303 304 305 306 307 308 309
this work, the FFT are the most computationally demanding. The algorithm of the FFT based on CUDA has been successfully designed. In addition, sampling by p-field, generating synthetic seismic data by convolution model, squared residual computation and all the node-based or trace-based calculations could be executed in parallel since the solution for each node or trace is independent of any other. In order to measure the performance of the proposed parallelization approaches, we applied our method in a workstation with Intel Xeon CPU E5-2620 and Nvidia Tesla K20. Several experiments indicate that compared with computation in CPU, the use of GPU reduces the computation time more than 10 times.
RI PT
300 301
3 APPLICATION
311
We first test the algorithm on a synthetic data set and then show an application with real field data. 3.1 Synthetic case The test is presented using a synthetic fluvial channel clastic reservoir model derived from the Stanford V model (Mao and Journel, 1999). The model consists of 100×130×450 grid nodes. The distance between the two traces is 25 m in both the x and y directions, and each trace has 450 samples at intervals of 1 ms. The model has three layers and is characterized by 3 facies (Figure 4): plain facies, sinuous sand-filled channels, and border sands consisting of levee and crevasse splays. In each layer, the channels have different widths, thicknesses and orientations. The spatial correlation of AI is assumed to be stationary and is characterized by an exponential autocorrelation function (equation (6)) with length a =500 m in the x direction, b =250 m in the y direction and c =10 ms in the vertical direction.
322 323
324 325 326 327 328
329 330
M AN U
TE D
319 320 321
1 x2 y 2 z 2 2 R ( x, y, z ) = exp − 2 + 2 + 2 a b c
(6)
A synthetic seismic dataset (Figure 5, left) is generated by the convolutional model with a 25 Hz Ricker wavelet. A Gaussian white noise with a variance α2=0.005 was added to the synthetic data. 15 wells (Figure 5, right) are extracted from real impedance model in the study area.
EP
313 314 315 316 317 318
AC C
312
SC
310
Figure 4. 3D sedimentary facies model and its horizontal cross-sections of each layer.
10
ACCEPTED MANUSCRIPT 331
RI PT
332 333
334
336 337 338 339 340 341
Figure 5. 3D seismic data (left) and well locations (right)
SC
335
356 357
Figure 6. Objective function against the number of iterations.
351 352 353
TE D
EP
348 349 350
AC C
342 343 344 345 346 347
M AN U
354 355
We then performed the inversion procedure for the whole 3D seismic dataset. The wavelet (25 Hz Ricker wavelet), facies model and the autocorrelation function were assumed to be the same as those used to generate the synthetic seismic data. We first used indicator cokriging to obtain the ccdf maps of the AI conditioned to the facies model and 15 known wells. Then, FFT-MA and p-field simulation was applied to generate an initial model of AI. The model is conditioned to the well data, geology and the spatial continuity but does not account for the seismic data. The corresponding synthetic seismic data were calculated, and the objective function was evaluated. An improved local GDM was applied to update the model to minimize the objective function. Two strategies are used together. The nodes which need to be modified are determined by the squared residual at each node between the real seismic data and synthetic seismic data. And we divided each trace into 15 parts. Therefore, the inversion includes 100×130×15 deformation parameters. We could obtain an optimal AI model, whose seismic response has a good match with the real seismic data and requires less than 300 iterations (Figure 6). The procedure is fast when using GPU parallelization, which costs only 5s for one iteration and approximately 1500s for the whole inversion. But when using CPU, it costs about 51s - 54s for one iteration, more than 10 times slower than using GPU. We then applied the methodology 100 times using different random seeds and calculated their posterior mean and variance of model.
11
ACCEPTED MANUSCRIPT
365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382
AC C
EP
383
RI PT
363 364
SC
362
M AN U
360 361
We compared the inverted AI model with the real model and evaluated the uncertainty. Figure 7 displays a 2D section of AI from the 3D inverted model, extracted along the red line in Figure 5 (right). Figure 8 is time slice of the AI at 350 ms. As we can see, the shape of the facies in the inverted model are preserved very well. This means that the geological information plays an important role in the inversion. The mean model is smooth, which recovers the spatially distributed character of AI, except for small-scale variability. It is similar to the result of the deterministic inversion. Each inverted AI realization is well matched with the real model on the seismic scale. Their small-scale heterogeneities were different but had the same statistical characteristics. The variance of model can show the uncertainty of each node. We observed that the traces at the well locations have the lowest variance, which means that our result is conditioned on the well data. At other locations, the uncertainty is related to the facies categories. The AI of the channel facies has the lowest uncertainty, while that of the plain facies is higher. This is because the channel facies have a narrower range of a priori probabilities than the plain facies. The border sands have the highest uncertainties. Since the number of samples of the border sands in the well is so small, we cannot capture the main statistical character of this class to obtain reasonable a priori probability. Therefore, the a priori information is an important factor to influence the uncertainty of the inversion results. Of course, the uncertainty is not only related to a priori information. In each facies category, the variance is not the same which means that the a posteriori probability is influenced by the seismic data as well. We then calculated the percentage decrease of the objective function and correlation coefficient to quantify the match between the synthetic and real seismic data. The value of the objective function decreases by approximately 80% and the correlation coefficient is 0.95, which means that we obtained very good results that are consistent with the seismic data.
TE D
358 359
384
12
ACCEPTED MANUSCRIPT
387 388
Figure 7. 2D section along the red line in Figure 5 of the real AI model and inverted AI model: (top left) real model, (mid-left) mean of 100 realizations, (bottom left) variance of 100 realizations, (right) 3 of 100 different realizations, where the triangles in each section stand for well locations.
389 390 391 392
Figure 8. Time slice (t=350 ms) of the real AI model and inverted AI model: (top left) real model, (top middle) mean of 100 realizations, (top right) variance of 100 realizations, (bottom) 3 of 100 different realizations.
402 403 404 405
406 407 408 409 410
M AN U
TE D
EP
399 400 401
At last, we analyzed the probability distribution of the inverted result. Figure 9 shows a priori probability histogram (top left) and posterior probability histogram (top right) of a grid node (x= 500 m, y= 2500 m, t= 350 ms) from reservoir model. This node belongs to the plain facies. We can observe that the a priori distribution is wide and doesn’t conform to Gaussian distribution. The posterior distribution become narrow after seismic inversion and are still non-Gaussian distribution. Another two histograms (bottom) in Figure 9 shows that the probability distribution of the entire real model (bottom left) is non-Gaussian and our method can reproduce this distribution successfully (bottom right). It is worth noting that there are no assumptions about the a priori or posterior distribution of the variables for each grid node and the distribution of the variables for the entire model being Gaussian distribution before inversion.
AC C
393 394 395 396 397 398
SC
RI PT
385 386
Figure 9. The probability distribution before and after inversion The top two histograms are a priori distribution obtained by indicator cokriging (top left) and posterior distribution after seismic inversion (top right) of one grid node in the reservoir model, respectively. The bottom two histograms are distribution of the real model (bottom left)
13
ACCEPTED MANUSCRIPT and distribution of one inverted realization (bottom right), respectively.
413 414
3.2 Real case The proposed inversion method was then applied to a real dataset in the Gulf of Mexico. The depositional environment is shelf edge delta sediments. The research area contains multiple horizons and a large fault (Figure 10, bottom right). The available data include 3D post-stack seismic data (Figure 10, left) and five wells (Figure 10, top right). The seismic data consist of 45 inlines, with 50 m intervals and 72 crosslines, with 25 m intervals. The vertical extent is from 1800 to 2400 ms with 4 ms intervals. Hydrocarbon bearing sand is out of the Green formation. Four lithofacies can be identified using gamma ray and resistivity: sand, shale, shaley sand, and oil sand. Since the oil sand with low gamma ray values and high resistivity has the lowest AI values (Figure 11), AI inversion can be used to identify the oil sand. For simplicity, we did not consider the influence of the fault and only use the horizons and 5 wells to obtain the ccdf. The wavelet estimation (Figure 10, bottom middle) and well ties were conducted first, using the model-driven wavelet estimation scheme. We consider the spatial correlation of the AI on two different scales. The large-scale correlation lengths are 1000 m in the inline direction, 1500 m in the crossline direction, and 16 ms in the vertical direction. This was considered using indicator cokriging to obtain a priori ccdf. The small-scale correlation lengths are 125 m along the inline direction, 250 m along the crossline direction, and 6 ms in the vertical direction. This was recovered by our proposed seismic inversion. In the inversion process, the squared residual between the real seismic data and synthetic seismic data is used to determine those nodes which need to be modified. The number of deformation parameters used in this application is 45×72×7, which means we divided each trace into 7 parts. Even though the model is smaller than the synthetic model in section 3.1, it still needed approximately 300 iterations (Figure 12) because the real data have lower quality and the parameters used in the inversion are not as perfect as the synthetic case. However, it took very little time, about 0.64 s for one iteration and 240 s to obtain an optimal realization, which is 10 times faster than 6.7 s for one iteration using CPU. We applied the methodology 100 times, using different random seeds, and then evaluated the uncertainty.
422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
AC C
440
SC
420 421
M AN U
418 419
TE D
416 417
EP
415
RI PT
411 412
441 442
Figure 10. Seismic data (right), well locations (top left), estimated wavelet (mid-bottom) and
14
ACCEPTED MANUSCRIPT 443 444
structural model in the geocellular grid (bottom right). The triangle in each section stands for the well location.
Figure 11. Well log data from well 4 (top) and corresponding crossplot (bottom). Oil sand with high resistivity and low gamma ray values can be identified by their lowest AI (points in the red rectangle).
AC C
EP
446 447 448 449
TE D
M AN U
SC
RI PT
445
450 451 452 453 454 455 456
Figure 12. Objective function against the number of iterations. Figure 13 and Figure14 show the inverted AI model in crossline direction and vertical direction, respectively, including the posterior mean, the variance of the inverted 100 realizations and 3 different realizations. They show that the average model is relatively smooth and has lost small-scale heterogeneities. All of the realizations are similar to the
15
ACCEPTED MANUSCRIPT
461 462 463 464 465 466 467 468 469 470
473 474 475 476 477
AC C
EP
TE D
M AN U
471 472
RI PT
459 460
smooth model on the large-scale and contain small-scale heterogeneities. They are statistically similar but are always different in their details. Geological information, such as horizons and faults, constrain the inverted AI, and the difference of the AI in the different horizons is recovered well. The variance of model suggests that different horizons have different uncertainties since the heterogeneities are different. It is obvious that the variance at the well locations is small and becomes bigger gradually as the node is farther. Some areas have high variance which suggests that it may contain thinner beds than seismic resolution. Since we cannot invert thin beds exactly, and these areas have high uncertainty. However, the vertical heterogeneity of each realization is shown to be similar to the well. We can observe that there are two areas with low AI values in the Green formation, which can be identified as oil sand body. Time slices through the Green formation (Figure 14) can help us to determine the hydrocarbon zones by low AI values. We also give a quantitative measurement of the results. The value of the objective function decreases by approximately 78%, and the correlation coefficient of the synthetic and real seismic data is 0.76. Therefore, this can be considered a satisfactory result for a real case.
SC
457 458
Figure 13. Crossline section of the inverted AI model along the red line in Figure 9: (left top) the mean of 100 realizations, (left bottom) variance of 100 realizations, (right) 3 of 100 different realizations, and the triangle in each section stands for the well location.
16
Figure 14. Time slices of the inverted AI model at 36ms below the top of the Green formation. (top left) the mean of 100 realizations, (top right) variance of 100 realizations, (bottom) 3 of 100 different realizations.
M AN U
478 479 480 481
SC
RI PT
ACCEPTED MANUSCRIPT
482
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507
TE D
EP
489 490 491
4 CONCLUSION We presented a new stochastic seismic inversion method based on an improved local GDM that combines seismic data, well log data, geological information and their spatial variability. The use of indicator cokriging allows us to incorporate the geological information and well data. FFT-MA and p-field simulation could generate fine-scaled conditional realizations consistent with spatial correlation information. Two improved strategies make local GDM successful use in seismic inversion. In addition, our proposed method can be used in non-Gaussian reservoir model and can be easily made parallel. The applications in the synthetic and real cases show that our method can generate reasonable AI models and give uncertainty estimations. It also has a high efficiency accelerated by the GPU. In future research, non-stationary natural phenomena should be studied since they are common in geology. Another issue is how to determine the number of parts that we should divide the trace into in the vertical direction. This may be influenced by the wavelet length, the autocorrelation length in the vertical direction and the number of samples in a trace.
AC C
483 484 485 486 487 488
Acknowledgments We gratefully acknowledge the editor and reviewers for the helpful comments and suggestions. REFERENCES Azevedo, L., Nunes, R., Soares, A., Mundin, E. C., Neto, G. S., 2015. Integration of well data into geostatistical seismic amplitude variation with angle inversion for facies estimation. Geophysics 80(6), M113-M128. Bortoli, L. J., Alabert, F. A., Haas, A., and Journel, A. G., 1993. Constraining stochastic images to seismic data. Geostatistics Tróia ’92(1), 325–337.
17
ACCEPTED MANUSCRIPT
519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551
RI PT
517 518
SC
515 516
M AN U
513 514
TE D
512
EP
510 511
Bosch, M., Mukerji, T., Gonzalez, E. F., 2010. Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: a review. Geophysics 75(5), 75A165-75A176. Buland, A., Kolbjørnsen, O., and Omre, H., 2003. Rapid spatially coupled AVO inversion in the Fourier domain. Geophysics 68, 824–836. doi: 10.1190/1.1581035. Buland, A., and Omre, H., 2003. Bayesian linearized AVO inversion. Geophysics 68, 185–198. doi: 10.1190/1.1543206 Cooke, D. A., and Schneider, W. A., 1983. Generalized linear inversion of reflection seismic data. Geophysics 48, 665-676. Deutsch, C. V., Journel, A. G., 1998. Gslib: Geostatistical Software Library and User's Guide: Oxford University Press. Doyen, P., 2007. Seismic Reservoir Characterization: EAGE. Eidsvik, J., Avseth, P., Omre, H., Mukerji, T., Mavko, G., 2004. Stochastic reservoir characterization using prestack seismic data. Geophysics 69, 978–993. doi: 10.1190/1.1778241. Ferreirinha, T., Nunes, R., Azevedo, L., Soares, A., Pratas, F., Tomas, P., Roma, N., 2015. Acceleration of stochastic seismic inversion in openCL-based heterogeneous platforms. Computers & Geosciences 78, 26-36. Francis, A. M., 2006. Understanding stochastic inversion. part 1: First Break 24, 69-77. Froidevaux R., 1993. Probability Field Simulation. Geostatistics Tróia ’92, 73-83. Gervais-Couplet, V., Roggero, F., Feraille, M.,Le Ravalec, M., Seiler, A., 2010. Joint History Matching of Production and 4D-Seismic Related Data for a North Sea Field Case. In: SPE Annual Technical Conference and Exhibition, 19-22 September, Florence, Italy González, E. F., Mukerji, T. and Mavko, G., 2008. Seismic inversion combining rock physics and multiple-point geostatistics. Geophysics 73(1), R11–R21. doi: 10.1190/1.2803748. Grana, D., and Della Rossa, E., 2010. Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion. Geophysics 75(3), O21–O37. doi: 10.1190/1.3386676. Grana, D., Mukerji, T., Dvorkin, J., Mavko, G., 2012. Stochastic inversion of facies from seismic data based on sequential simulations and probability perturbation method. Geophysics 77(4), 53. Grijalba-Cuenca, A., Torres-Verdín, C., and Debeye, H. W. J., 2000. Geostatistical inversion of 3D seismic data to extrapolate wireline petrophysical variables laterally away from the well. In: SPE Annual Technical Conference and Exhibition. Haas, A., and Dubrule, O., 1994. Geostatistical inversion-a sequential method of stochastic reservoir modelling constrained by seismic data. First break 12, 561-569. Hu, L. Y., 2000, Gradual deformation and iterative calibration of Gaussian-related stochastic models. Mathematical Geology 32, 87-108. Hu, L. Y., and Le Ravalec-Dupin, M., 2004. An improved gradual deformation method for reconciling random and gradient searches in stochastic optimizations. Mathematical Geology 36, 703-719. Journel, A. G., 2002. Combining knowledge from diverse sources: an alternative to traditional data independence hypotheses. Mathematical Geosciences 34(5), 573-596. Le Gallo, Y., and Le Ravalec-Dupin, M., 2000. History matching geostatistical reservoir
AC C
508 509
18
ACCEPTED MANUSCRIPT
563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595
RI PT
561 562
SC
559 560
M AN U
557 558
TE D
556
EP
554 555
models with gradual deformation method. In: SPE Annual Technical Conference and Exhibition. Le Ravalec, M., Hu, L. Y., Noetinger, B., 1999. Stochastic reservoir modeling constrained to dynamic data: Local calibration and inference of the structural parameters. In: SPE Annual Technical Conference and Exhibition. Le Ravalec, M., Noetinger, B., and Hu, L. Y., 2000. The FFT moving average (FFT-MA) generator: an efficient numerical method for generating and conditioning Gaussian simulations. Mathematical Geology 32, 701-723. Le Ravalec, M., Noetinger, B., 2002. Optimization with the gradual deformation method. Mathematical Geology 34 (2), 125-142. Le Ravalec-Dupin, M., Coureaud, B., Nicolas, L., Roggero, F., 2004. Conditioning an underground gas storage site to well pressures. Oil & gas science and technology 59 (6), 611-624. Le Ravalec, M., Mouche, E., 2012. Calibrating transmissivities from piezometric heads with the gradual deformation method: An application to the Culebra Dolomite unit at the Waste Isolation Pilot Plant (WIPP), New Mexico, USA. Journal of Hydrology 472, 1-13. Marteau, B., Ding, D. Y., and Dumas, L., 2014. A generalization of the local gradual deformation method using domain parameterization. Computers & Geosciences 72, 233-243. Mukerji, T., Jorstad, A., Avseth, P., Mavko, G., Granli, J. R., 2001. Mapping lithofacies and pore-fluid probabilities in a North Sea reservoir: Seismic inversions and statistical rock physics. Geophysics 66, 988–1001. doi: 10.1190/1.1487078. Nunes, R., Almeida, J. A., 2010. Parallelization of sequential Gaussian, indicator and direct simulation algorithms. Computers & Geosciences 36(8), 1042-1052. Nunes, R., Soares, A., Schwedersky, G., Dillon, L., Guerreiro, L., Caetano, H., Maciel, C., Leon, F., 2012. Geostatistical inversion of prestack seismic data. In: Ninth International Geostatistics Congress, Oslo, Norway, 1-8. Nunes, R., Soares, A., Azevedo, L., Pereira, P., 2016. Geostatistical seismic inversion with direct sequential simulation and co-simulation with multi-local distribution functions. Mathematical Geosciences, 1-19. Oldenburg, D. W., Scheuer, T., Levy, S., 1983, Recovery of the acoustic impedance from reflection seismograms, Geophysics 48, 1318-1337. Oliver, D. S., 1995. Moving averages for Gaussian simulation in two and three dimensions. Mathematical Geology 27, 939-960. Pourpak, H., Bourbiaux, B., Roggero, F., Delay, F, 2009. An integrated method for calibrating a heterogeneous/fractured reservoir model from wellbore flow measurements: Case study. SPE Reservoir Evaluation & Engineer 12, 433-445. Sanders, J., Kandrot, E., 2010. CUDA by Example: An introduction to general-purpose GPU programming. Addison-Wesley Professional. Soares, A., 2001. Direct sequential simulation and cosimulation. Mathematical Geosciences 33(8), 911-926. Soares, A., Diet, J. D., Guerreiro, L., 2007. Stochastic inversion with a global perturbation method. In: EAGE Conference on Petroleum Geostatistics. Srivastava, R. M., 1992. Reservoir characterization with probability field simulation. SPE
AC C
552 553
19
ACCEPTED MANUSCRIPT
598 599 600 601 602 603 604
Annual Technical Conference and Exhibition, 4-7 October, Washington, D.C. Tahmasebi, P., Sahimi, M., Mariethoz, G., Hezarkhani, A., 2012. Accelerating geostatistical simulations using graphics processing units (GPU). Computers & Geosciences 46(3), 51-59. Tarantola, A. 2005. Inverse problem theory and methods for model parameter estimation. Society for Industrial and Applied Mathematics. Vargas, H. S., Caetano, H., Filipe, M., 2007. Parallelization of sequential simulation procedures. In: EAGE Conference on Petroleum Geostatistics.
RI PT
596 597
605 606
AC C
EP
TE D
M AN U
SC
607
ACCEPTED MANUSCRIPT Research highlights A seismic inversion method based on the local gradual deformation method is proposed. The method can incorporate seismic data, well data, geology and spatial correlations.
RI PT
Probability field simulation and FFT moving average is combined to generate models.
AC C
EP
TE D
M AN U
SC
Two improved strategies of gradual deformation method are proposed.
S0098-3004(16)30506-4
DOI:
10.1016/j.cageo.2017.08.010
Reference:
CAGEO 4007
To appear in:
Computers and Geosciences
Received Date: 8 October 2016 Revised Date:
5 August 2017
Accepted Date: 8 August 2017
Please cite this article as: Yang, X., Zhu, P., Stochastic seismic inversion based on an improved local gradual deformation method, Computers and Geosciences (2017), doi: 10.1016/j.cageo.2017.08.010. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
ACCEPTED MANUSCRIPT 1
Stochastic seismic inversion based on an improved local gradual deformation method Xiuwei Yang and Peimin Zhu China University of Geosciences, Institute of Geophysics and Geomatics, Wuhan, China. E-mail:[email protected]; zhupm@ cug.edu.cn
7 8 9
ABSTRACT A new stochastic seismic inversion method based on the local gradual deformation method is proposed, which can incorporate seismic data, well data, geology and their spatial correlations into the inversion process. Geological information, such as sedimentary facies and structures, could provide significant a priori information to constrain an inversion and arrive at reasonable solutions. The local a priori conditional cumulative distributions at each node of model to be inverted are first established by indicator cokriging, which integrates well data as hard data and geological information as soft data. Probability field simulation is used to simulate different realizations consistent with the spatial correlations and local conditional cumulative distributions. The corresponding probability field is generated by the fast Fourier transform moving average method. Then, optimization is performed to match the seismic data via an improved local gradual deformation method. Two improved strategies are proposed to be suitable for seismic inversion. The first strategy is that we select and update local areas of bad fitting between synthetic seismic data and real seismic data. The second one is that we divide each seismic trace into several parts and obtain the optimal parameters for each part individually. The applications to a synthetic example and a real case study demonstrate that our approach can effectively find fine-scale acoustic impedance models and provide uncertainty estimations.
24
SC
14 15 16 17 18 19 20 21 22 23
M AN U
12 13
TE D
10 11
RI PT
2 3 4 5 6
25
Keywords: Stochastic inversion; Heterogeneities; Probability field simulation; Gradual deformation method.
29
1 INTRODUCTION
30 31
Seismic inversion plays an important role in reservoir prediction. The inverted impedance from seismic data can indicate changes in reservoir properties. However, an accurate impedance model is difficult to invert since seismic data is bandwidth limited. Both low frequencies and high frequencies outside of the seismic bandwidth cannot be recovered from seismic data alone. Additional constraints are used to reduce the non-uniqueness of the inversion problem by providing a priori information that is not contained in the seismic data. Geological information, such as sedimentary facies and structures, could provide significant a priori information to constrain the inversion and arrive at reasonable solutions (Gonzalez et al., 2008; Nunes et al., 2016). A conventional seismic inversion, such as generalized linear inversion (Cooke and Schneider, 1983) or sparse spike inversion (Oldenburg et al., 1983), usually gives a best estimate of the impedance that minimizes the difference between synthetic seismic traces and real seismic traces. Such techniques yield a single estimate, which should be the mean of all possible non-unique solutions and relatively smooth. Thus, the deterministic seismic inversion
34 35 36 37 38 39 40 41 42 43
AC C
32 33
EP
26 27 28
2
ACCEPTED MANUSCRIPT
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
RI PT
53 54
SC
51 52
M AN U
49 50
TE D
48
EP
46 47
is only suitable when the reservoir is comprised of relatively thick layers with slow lateral thickness changes and minimal spatial impedance variation (Francis, 2006). It is worth noting that this method loses the small-scale heterogeneities and overestimates the connectivity of the parameters. Geostatistical methods could simulate small-scale variabilities not captured in seismic data. Geostatistical inversion (Doyen, 2007 and Bosch et al., 2010) provides a better strategy for generating a suite of alternative subsurface property realizations conditioned on different kinds of data and incorporates spatial correlation relationships. It is more suitable for use in uncertainty estimations, volumetric calculations or connectivity estimations since it includes the small-scale heterogeneities and spatial variabilities of the geology in the inverted model. Geostatistical inversion using sequential Gaussian simulation was first tested by Bortoli et al. (1993) and Haas and Dubrule (1994). To be constrained to seismic traces, they simulate a whole trace at a time instead of the traditional node-by-node sequential Gaussian simulation, and use a simple Monte Carlo rejection scheme to obtain the optimal trace. Grijalba-Cuenca et al. (2000) offered another geostatistical inversion algorithm, which worked node by node and was optimized using the simulated annealing method. Mukerji et al. (2001) and Eidsvik et al. (2004) introduced statistical rock physics to estimate reservoir parameters from prestack seismic data and evaluated the associated uncertainty. In recent years, some new geostatistical techniques and inversion algorithms have been introduced in seismic inversion. A global iterative geostatistical seismic inversion method (Soares et al., 2007, Nunes et al. 2012, Azevedo et al. 2015), using direct sequential simulation and co-simulations (Soares, 2001), has been proposed. Several models for the entire seismic grid are simulated at once in the first step. Local areas of best fit of the different models are selected and “merged” into a secondary model for the direct co-simulation of the next iteration. Gonzalez et al. (2008) incorporated multiple-point geostatistics and rock physics to characterize previously known geologic information, simulate reservoir models and then select the optimal realization. Grana et al. (2012) estimated the facies and reservoir properties based on the probability perturbation method and tau model (Journel, 2002). Buland and Omre (2003) and Buland et al. (2003) developed a linearized AVA inversion in a Bayesian framework, both in the spatial and Fourier domains. They assume that the elastic parameters are characterized by a log-Gaussian random field and an explicit analytical form of the Gaussian posterior distribution can be efficiently obtained. Grana and Della Rossa (2010) combined this Bayesian linear AVA inversion with statistical rock physics to estimate reservoir properties and the associated uncertainty. The use of Bayesian framework made the link between stochastic inversions and deterministic inversions clearer. The fast Fourier transform moving average (FFT-MA) method (Le Ravalec et al., 2000) is the implementation of the moving average (MA) method (Oliver, 1995) in the frequency domain. The FFT-MA method, combined with the gradual deformation method (GDM) (Le Ravalec et al., 1999; Hu, 2000), has been successfully used in petroleum engineering. This approach uses a perturbation mechanism to gradually change the Gaussian-related reservoir model and guarantees the spatial correlation. Local perturbation (Le Ravalec et al., 1999; Hu, 2000; Marteau et al., 2014) and gradient information (Hu and Le Ravalec-Dupin, 2004) are used to improve the convergence rate. Another common practice to speed up convergence consists in combining several realizations instead of two, which implies the use of several
AC C
44 45
3
ACCEPTED MANUSCRIPT
113
2 METHODOLOGY
114 115 116
The methodology we proposed can integrate the geological a priori information and their spatial correlations in seismic inversion. The workflow of the method is shown in Figure 1. An indicator cokriging (Deutsch and Journel, 1998) is first applied to take into account well data and geological a priori information and to obtain the ccdf of the acoustic impedance (AI) at each grid node of the model to be inverted. Sampled from the locally defined ccdf, the p-field simulation will then generate different realizations consistent with same spatial correlation. A spatially correlated 3D p-field is easily simulated using the FFT-MA technique. Finally, an improved local gradual deformation method (GDM) is applied to perturb the Gaussian white noise locally until there is a match with the seismic data. Two strategies are proposed to make GDM suitable for seismic inversions. Each step of the method used in this research is described in the following subsections.
99 100 101 102 103 104 105 106
117 118 119 120 121 122 123 124 125
SC
97 98
M AN U
95 96
TE D
93 94
EP
92
AC C
90 91
RI PT
107 108 109 110 111 112
deformation parameters (Le Ravalec-Dupin and Noetinger, 2002). GDM has been successfully applied to several synthetic and real datasets (Le Gallo and Le Ravalec-Dupin, 2000; Le Ravalec-Dupin et al. 2004; Pourpak et al., 2009; Gervais-Couplet et al. 2010; Le Ravalec and Mouche, 2012) for history matching, which characterizes a reservoir by matching the static (e.g., porosity, permeability) and dynamic data (e.g., pressure, oil rate) from the wells. Data matching in seismic inversion is different from history matching, since seismic data is defined by millions of points and covers the whole study area. There are problems reducing the mismatch of amplitudes between the numerical seismic responses and the real seismic data exactly using original local GDM. Therefore, new strategies are necessary to develop. In this paper, geological information is included by indicator cokriging to obtain the local conditional cumulative distribution function (ccdf). In order to simulate fine-scale reservoir models consistent with the expected geology, well data and their spatial correlation efficiently, FFT-MA and probability field (p-field) simulation (Srivastava, 1992; Froidevaux, 1993) are used. Then, optimization is performed to match the seismic data via an improved local gradual deformation method. In order to match seismic amplitude at each grid-node, two improved strategies are proposed to make original GDM suitable for seismic inversions. The first strategy is that we select and update local areas of bad fitting between synthetic seismic data and real seismic data. The second one is that we divide each seismic trace into several parts and obtain the optimal parameters for each part individually. We can apply the methodology several times using different random numbers to generate different realizations. Then posterior distribution can be calculated using all these inverted realizations and give uncertainty estimation. We begin by briefly describing the workflow of this approach and then discuss each step in detail, especially the two improved strategies of GDM. The method is first tested on a 3D synthetic case and then on a real data.
88 89
4
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
126 127 128
Figure 1. Workflow of the stochastic inversion.
142
r be a series of cutoff values to discretize I (u ) ; i (uα ; r ) is the indicator transformation of
143
I (u ) ; α stands for the index of AI ( α = 1, 2, ... , n); and y (u; r ) is the a priori probability
144
conditioned on geological information. By solving this kriging system and getting the
145
coefficients λ0 (u , r ) , λα (u; r ) , and v(u; r ) , attributed to F (r ) , i(uα ; r ) and y (u; r ) ,
146
respectively, the a posteriori probability of I (u ) can be obtained using:
137
EP
135 136
AC C
132 133 134
TE D
138 139 140 141
2.1 Local conditional cumulative distribution Generally, the probability function of a continuous variable is often assumed to be Gaussian in modeling or inversion methods. However, the distribution of model parameters in geology is often not close to being Gaussian. Sometimes, the distribution even cannot be expressed using a set of parameters and an explicit analytical form. In this case, the indicator kriging would be selected for a so-called non-parametric estimation of probability distribution of variables. It has the advantage that any distribution can be handled. We first discretize continuous variable AI into K classes by (K-1) cutoff values. Then, a linear interpolation is performed to recover the within-class resolution lost within each separate subclass. K should be a reasonable number that could guarantee each class includes adequate samples from well data. And meanwhile, the inferred ccdf can characterize the main statistical features of AI. In this paper, indicator cokriging using the Markov-Bayes model is used to obtain the local ccdf, conditional on geology and well data. Let I (u ) be the AI value at location u and
129 130 131
n
147
Prob( I (u ) < r | (n + 1)) = λ0 (u , r ) F (r ) + ∑ λα (u; r )i (uα ; r ) + v(u; r ) y (u; r ) α =1
(1)
5
ACCEPTED MANUSCRIPT 148 149
where (n + 1) represents the number of conditioning data including n neighboring AI data and 1 co-located a priori probability. F (r ) is the global a priori probability common to all
150
locations u . These steps only need to be carried out one time before the simulation and inversion; thus, this method is sufficiently simple and fast to handle a large 3D case.
151 152 153 154
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
RI PT
SC
162 163 164 165 166 167
M AN U
160 161
TE D
159
EP
157 158
2.2 FFT-MA and p-field simulation Probability field (p-field) simulation (Srivastava, 1992; Froidevaux, 1993) is a fast conditional simulation method. Contrary to the sequential methods, which need to reconstruct ccdf for each realization at every location, the ccdf of the p-field simulation could remain the same from one realization to another, only conditional on the original data. After we obtained this local ccdfs, the conditional simulations are then obtained by drawing samples from these ccdfs. The probability values used to draw from the local ccdfs constitute a probability field, and are viewed as outcomes of a random function characterized by a uniform distribution and a given covariance function. In this paper, this correlated p-field is quickly constructed by the fast Fourier transform moving average (FFT-MA) method (Le Ravalec et al., 2000). The p-field realizations need not be conditional since the ccdfs are already conditioned on the original data. Since the model generated by FFT-MA follows a Gaussian distribution with a mean of 0, unit variance and autocorrelation ( C ), we have to transform it to a uniform distribution between 0 and 1 in order to ensure the reservoir model satisfies the ccdfs we have obtained. The use of FFT-MA and p-field simulation allows us to generate realizations satisfied by any specified ccdf and spatial correlation. The simulations are fast and stable because construction of the p-field is carried out through FFT, and the ccdf is produced only once, conditional on the original data. Another important feature is that the FFT-MA generator draws the randomness in the spatial domain, which allows us to uncouple the random numbers from the structural parameters. Such features imply that any of these random numbers can be perturbed locally without affecting the others outside of the range of the correlation length. It is potentially of great interest for seismic inversion to modify the AI model locally, based on a comparison between the actual seismic data and synthetic seismic responses of the inverted AI.
AC C
155 156
2.3 Improved local gradual deformation method The model constructed by the p-field simulation has been conditioned on well log data but not conditioned on seismic data. Therefore, we need an inversion to make the result match the seismic data. Generally, a seismic impedance inversion is non-unique. The classical gradient seismic inversion method gives only one smooth result and loses small-scale heterogeneity. Thus, the stochastic inversion is more suitable to obtain different realizations and estimate the uncertainty. The gradual deformation method (GDM) has turned out to be very consistent with FFT-MA and has been successfully used in history matching. Since seismic data is different from production data, two improvements are proposed to make the method suitable for seismic inversion. In geophysics inversion, least-absolute-values criterion (L1 norm) and least-squares
6
ACCEPTED MANUSCRIPT
criterion (L2 norm) are two of the most widely used criterions. The choice of the norm should be driven by the type of data and the purpose of the application. When outliers are suspected in a data set, long-tailed probability density functions should be used to model uncertainties. In this case, L1 norm is suitable since its result is robust to outliers. If uncertainties in the problem can be modeled using Gaussian distributions, L2 norm is a more appropriate choice (Tarantola, 2005). In most case of seismic inversion, the errors of observed data are proved to be a Gaussian distribution or generalized Gaussian distribution. Therefore, we first define the objective function J as the sum of the squared errors of the seismic responses S m and
200
observed data S mobs . J=
201
204 205
1 M ∑ (Sm − Smobs )2 2 m =1
(2)
where M represents the number of sample points of the seismic trace. The lengths of the observed and synthetic seismic data are the same since the inversion is conducted in time domain without the conversion between time and depth. The synthetic seismic trace S is computed by convolving the wavelet w with the reflection coefficients series r :
S = r∗w
SC
202 203
RI PT
192 193 194 195 196 197 198 199
(3) The core idea of GDM is that a new set of independent Gaussian white noise z can be
208
created using two sets of independent Gaussian white noise, denoted z1 and z 2 ,
209
respectively. z1 is the Gaussian white noise for the initial realization and z 2 is another
210
randomly generated one. In addition, their gradual deformation relation is defined as follows:
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232
TE D
216 217
(4)
Then, a search process is implemented to determine the deformation parameters r minimizing Eq. (2). At this point, the fit may not be good enough. A new starting realization chain is constructed from the optimal realization identified at the end of the previous search process. In this way, we only need to create a sequence of independent Gaussian white noise, generate realizations via FFT-MA and p-field simulation, and estimate the optimal r to further minimize the objective function. This process is iterated until a satisfactory matching is achieved. The global GDM does not have a high efficiency because improving the fit in one region may deteriorate the fit in another region. In history matching, a local gradual deformation is developed by partitioning the field into several independent regions and performing the inversion individually for each region. Deterministic partitioning is available in history matching because the observed data (e.g., pressure, oil rate) are scattered in different zones of the field. However, the seismic data are different from the production data because they are defined by millions of points and cover the whole area, aside from some sparse locations. In order to match the seismic amplitude at every grid-node exactly, we developed two improved strategies to make GDM suitable for seismic inversion. The first strategy is to select local areas with bad fits to modify. At the beginning of the inversion, we can modify all of the nodes, and the objective function will decrease quickly. When the convergence slows down, modifying all the nodes at one time is not a good choice. In some regions, if we have obtained a good fit, there is no need to modify the impedance model. Therefore, we only need to update the regions with poor fits. The squared residual at
EP
212 213 214 215
z = z1 cos(π r ) + z 2 sin(π r ), r ∈[-1,1]
AC C
211
M AN U
206 207
7
ACCEPTED MANUSCRIPT
240
z = (1 − H).× z1 + H. × [z1. × cos(π r ) + z 2 . × sin(π r )], r ∈[-1,1]
241 242 243 244 245
where the symbol .× stands for the componentwise product. H is the binary mask. When the convergence slows down, it is difficult to obtain a better model to match the seismic data using this threshold. At this time, we can increase the threshold value T, and the number of nodes that need to be modified is decreased. The modification will focus on the nodes that have the worse fitting and the model will be further refined. The same optimized process is repeated again and again. The modified region decreases to consider the local misfit of the data.
250 251 252 253 254 255 256 257 258 259 260 261 262 263
Figure 2. Binary mask (bottom) defining the modified region obtained from the squared residual image (top) and a threshold.
AC C
249
EP
TE D
246 247 248
(5)
SC
237
M AN U
235 236
RI PT
238 239
each node between the real seismic data and synthetic seismic data stands for the mismatch. Therefore, it can be used to determine the nodes which need to be modified. Thus, at each iteration, we first calculate the squared residual at each node and define a threshold T. When the squared residual of one node is larger than T, this node needs to be modified. We transform the squared residual image into a binary mask using the threshold T, and the modified region is associated with the white pixels (Figure 2). At this time, the GDM equation can be expressed as
233 234
The second strategy is to divide each trace into several parts and to obtain the optimal parameters for each part individually. Even though we have selected local areas with bad fitting to modify, one deformation parameter for all the nodes is still not suitable. A local GDM with multiple parameters is a better choice, but a reasonable partitioning is difficult. We propose a strategy to divide the model and obtain the optimal parameters. In the lateral direction, we assume each trace is independent, since post-stack seismic data has eliminated the trace dependence. As such, each trace is a single subregion with one deformation parameter. In the vertical direction, the amplitude at one node is only influenced by the impedance within the range of the wavelet length. Therefore, we simply divide each trace into several parts. The spatial continuity in the lateral direction and between the different parts in the vertical direction can still be guaranteed with the use of FFT-MA. Therefore, the number
8
ACCEPTED MANUSCRIPT
268 269 270 271
RI PT
266 267
of deformation parameters is n = Nx × Ny × Mz, where Nx represents the number of inline, Ny represents the number of crossline, and Mz is the number of parts in the vertical direction. Then, we solve an optimization problem of the n parameters, t1, t2, …, tn. We choose the golden section search method for all parameters simultaneously and independently based on the local objective function. In this way, the computing time is approximately equivalent to that of one parameter. In fact, when the optimization is done, the optimal parameters are able to divide the whole model into several subregions automatically (Figure 3). The number and shape of the subregions is flexible and different each time.
SC
264 265
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
TE D
281 282 283
Therefore, the optimization of our proposed method is performed within two nested loops. In the outer loop, the whole field is divided into n parts. In each part, the modified nodes are defined based on the squared residuals and threshold T. In the inner loop, we perform an optimization on the n parameters of the gradual deformation method to minimize the objective function. If the objective function of the new model is less than that of the previous model, the new model is accepted. Otherwise, we generate a new Gaussian white noise and repeat the previously described steps. We iterate over this procedure until the objective function is less than a tolerance value.
EP
275 276 277 278 279 280
Figure 3. Optimal parameters of the local GDM (left) are equivalent to the partitions and its parameters in the image to the right. Each rectangle stands for one part, including several samples of one trace in the vertical direction.
2.4 Parallelization In recent years, many parallel computational strategies for sequential simulation have been developed at three distinct levels: the realization level, path level and node level. Multi-core Central Processing Units (CPUs) (Vargas et al., 2007, Nunes and Almeida, 2010) and graphics processing units (GPUs) (Tahmasebi, 2012, Ferreirinha, 2015) are both exploited to improve the calculation efficiency. However, their sequential nature makes them difficult to parallelize very efficiently. The simulation of N nodes cannot be simply parallelized because there are dependences between the sequential and parallel results. The parallelization of our proposed approach is easier than the inversion using a sequential simulation. We exploit the compute unified device architecture (CUDA) to increase the computational efficiency. CUDA (Sanders and Kandrot, 2010) is a parallel computing architecture developed by the technology company Nvidia (2008) that makes it possible to benefit from the ability of GPU in parallel computing with general-purpose computation. In
AC C
273 274
M AN U
272
9
ACCEPTED MANUSCRIPT
302 303 304 305 306 307 308 309
this work, the FFT are the most computationally demanding. The algorithm of the FFT based on CUDA has been successfully designed. In addition, sampling by p-field, generating synthetic seismic data by convolution model, squared residual computation and all the node-based or trace-based calculations could be executed in parallel since the solution for each node or trace is independent of any other. In order to measure the performance of the proposed parallelization approaches, we applied our method in a workstation with Intel Xeon CPU E5-2620 and Nvidia Tesla K20. Several experiments indicate that compared with computation in CPU, the use of GPU reduces the computation time more than 10 times.
RI PT
300 301
3 APPLICATION
311
We first test the algorithm on a synthetic data set and then show an application with real field data. 3.1 Synthetic case The test is presented using a synthetic fluvial channel clastic reservoir model derived from the Stanford V model (Mao and Journel, 1999). The model consists of 100×130×450 grid nodes. The distance between the two traces is 25 m in both the x and y directions, and each trace has 450 samples at intervals of 1 ms. The model has three layers and is characterized by 3 facies (Figure 4): plain facies, sinuous sand-filled channels, and border sands consisting of levee and crevasse splays. In each layer, the channels have different widths, thicknesses and orientations. The spatial correlation of AI is assumed to be stationary and is characterized by an exponential autocorrelation function (equation (6)) with length a =500 m in the x direction, b =250 m in the y direction and c =10 ms in the vertical direction.
322 323
324 325 326 327 328
329 330
M AN U
TE D
319 320 321
1 x2 y 2 z 2 2 R ( x, y, z ) = exp − 2 + 2 + 2 a b c
(6)
A synthetic seismic dataset (Figure 5, left) is generated by the convolutional model with a 25 Hz Ricker wavelet. A Gaussian white noise with a variance α2=0.005 was added to the synthetic data. 15 wells (Figure 5, right) are extracted from real impedance model in the study area.
EP
313 314 315 316 317 318
AC C
312
SC
310
Figure 4. 3D sedimentary facies model and its horizontal cross-sections of each layer.
10
ACCEPTED MANUSCRIPT 331
RI PT
332 333
334
336 337 338 339 340 341
Figure 5. 3D seismic data (left) and well locations (right)
SC
335
356 357
Figure 6. Objective function against the number of iterations.
351 352 353
TE D
EP
348 349 350
AC C
342 343 344 345 346 347
M AN U
354 355
We then performed the inversion procedure for the whole 3D seismic dataset. The wavelet (25 Hz Ricker wavelet), facies model and the autocorrelation function were assumed to be the same as those used to generate the synthetic seismic data. We first used indicator cokriging to obtain the ccdf maps of the AI conditioned to the facies model and 15 known wells. Then, FFT-MA and p-field simulation was applied to generate an initial model of AI. The model is conditioned to the well data, geology and the spatial continuity but does not account for the seismic data. The corresponding synthetic seismic data were calculated, and the objective function was evaluated. An improved local GDM was applied to update the model to minimize the objective function. Two strategies are used together. The nodes which need to be modified are determined by the squared residual at each node between the real seismic data and synthetic seismic data. And we divided each trace into 15 parts. Therefore, the inversion includes 100×130×15 deformation parameters. We could obtain an optimal AI model, whose seismic response has a good match with the real seismic data and requires less than 300 iterations (Figure 6). The procedure is fast when using GPU parallelization, which costs only 5s for one iteration and approximately 1500s for the whole inversion. But when using CPU, it costs about 51s - 54s for one iteration, more than 10 times slower than using GPU. We then applied the methodology 100 times using different random seeds and calculated their posterior mean and variance of model.
11
ACCEPTED MANUSCRIPT
365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382
AC C
EP
383
RI PT
363 364
SC
362
M AN U
360 361
We compared the inverted AI model with the real model and evaluated the uncertainty. Figure 7 displays a 2D section of AI from the 3D inverted model, extracted along the red line in Figure 5 (right). Figure 8 is time slice of the AI at 350 ms. As we can see, the shape of the facies in the inverted model are preserved very well. This means that the geological information plays an important role in the inversion. The mean model is smooth, which recovers the spatially distributed character of AI, except for small-scale variability. It is similar to the result of the deterministic inversion. Each inverted AI realization is well matched with the real model on the seismic scale. Their small-scale heterogeneities were different but had the same statistical characteristics. The variance of model can show the uncertainty of each node. We observed that the traces at the well locations have the lowest variance, which means that our result is conditioned on the well data. At other locations, the uncertainty is related to the facies categories. The AI of the channel facies has the lowest uncertainty, while that of the plain facies is higher. This is because the channel facies have a narrower range of a priori probabilities than the plain facies. The border sands have the highest uncertainties. Since the number of samples of the border sands in the well is so small, we cannot capture the main statistical character of this class to obtain reasonable a priori probability. Therefore, the a priori information is an important factor to influence the uncertainty of the inversion results. Of course, the uncertainty is not only related to a priori information. In each facies category, the variance is not the same which means that the a posteriori probability is influenced by the seismic data as well. We then calculated the percentage decrease of the objective function and correlation coefficient to quantify the match between the synthetic and real seismic data. The value of the objective function decreases by approximately 80% and the correlation coefficient is 0.95, which means that we obtained very good results that are consistent with the seismic data.
TE D
358 359
384
12
ACCEPTED MANUSCRIPT
387 388
Figure 7. 2D section along the red line in Figure 5 of the real AI model and inverted AI model: (top left) real model, (mid-left) mean of 100 realizations, (bottom left) variance of 100 realizations, (right) 3 of 100 different realizations, where the triangles in each section stand for well locations.
389 390 391 392
Figure 8. Time slice (t=350 ms) of the real AI model and inverted AI model: (top left) real model, (top middle) mean of 100 realizations, (top right) variance of 100 realizations, (bottom) 3 of 100 different realizations.
402 403 404 405
406 407 408 409 410
M AN U
TE D
EP
399 400 401
At last, we analyzed the probability distribution of the inverted result. Figure 9 shows a priori probability histogram (top left) and posterior probability histogram (top right) of a grid node (x= 500 m, y= 2500 m, t= 350 ms) from reservoir model. This node belongs to the plain facies. We can observe that the a priori distribution is wide and doesn’t conform to Gaussian distribution. The posterior distribution become narrow after seismic inversion and are still non-Gaussian distribution. Another two histograms (bottom) in Figure 9 shows that the probability distribution of the entire real model (bottom left) is non-Gaussian and our method can reproduce this distribution successfully (bottom right). It is worth noting that there are no assumptions about the a priori or posterior distribution of the variables for each grid node and the distribution of the variables for the entire model being Gaussian distribution before inversion.
AC C
393 394 395 396 397 398
SC
RI PT
385 386
Figure 9. The probability distribution before and after inversion The top two histograms are a priori distribution obtained by indicator cokriging (top left) and posterior distribution after seismic inversion (top right) of one grid node in the reservoir model, respectively. The bottom two histograms are distribution of the real model (bottom left)
13
ACCEPTED MANUSCRIPT and distribution of one inverted realization (bottom right), respectively.
413 414
3.2 Real case The proposed inversion method was then applied to a real dataset in the Gulf of Mexico. The depositional environment is shelf edge delta sediments. The research area contains multiple horizons and a large fault (Figure 10, bottom right). The available data include 3D post-stack seismic data (Figure 10, left) and five wells (Figure 10, top right). The seismic data consist of 45 inlines, with 50 m intervals and 72 crosslines, with 25 m intervals. The vertical extent is from 1800 to 2400 ms with 4 ms intervals. Hydrocarbon bearing sand is out of the Green formation. Four lithofacies can be identified using gamma ray and resistivity: sand, shale, shaley sand, and oil sand. Since the oil sand with low gamma ray values and high resistivity has the lowest AI values (Figure 11), AI inversion can be used to identify the oil sand. For simplicity, we did not consider the influence of the fault and only use the horizons and 5 wells to obtain the ccdf. The wavelet estimation (Figure 10, bottom middle) and well ties were conducted first, using the model-driven wavelet estimation scheme. We consider the spatial correlation of the AI on two different scales. The large-scale correlation lengths are 1000 m in the inline direction, 1500 m in the crossline direction, and 16 ms in the vertical direction. This was considered using indicator cokriging to obtain a priori ccdf. The small-scale correlation lengths are 125 m along the inline direction, 250 m along the crossline direction, and 6 ms in the vertical direction. This was recovered by our proposed seismic inversion. In the inversion process, the squared residual between the real seismic data and synthetic seismic data is used to determine those nodes which need to be modified. The number of deformation parameters used in this application is 45×72×7, which means we divided each trace into 7 parts. Even though the model is smaller than the synthetic model in section 3.1, it still needed approximately 300 iterations (Figure 12) because the real data have lower quality and the parameters used in the inversion are not as perfect as the synthetic case. However, it took very little time, about 0.64 s for one iteration and 240 s to obtain an optimal realization, which is 10 times faster than 6.7 s for one iteration using CPU. We applied the methodology 100 times, using different random seeds, and then evaluated the uncertainty.
422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
AC C
440
SC
420 421
M AN U
418 419
TE D
416 417
EP
415
RI PT
411 412
441 442
Figure 10. Seismic data (right), well locations (top left), estimated wavelet (mid-bottom) and
14
ACCEPTED MANUSCRIPT 443 444
structural model in the geocellular grid (bottom right). The triangle in each section stands for the well location.
Figure 11. Well log data from well 4 (top) and corresponding crossplot (bottom). Oil sand with high resistivity and low gamma ray values can be identified by their lowest AI (points in the red rectangle).
AC C
EP
446 447 448 449
TE D
M AN U
SC
RI PT
445
450 451 452 453 454 455 456
Figure 12. Objective function against the number of iterations. Figure 13 and Figure14 show the inverted AI model in crossline direction and vertical direction, respectively, including the posterior mean, the variance of the inverted 100 realizations and 3 different realizations. They show that the average model is relatively smooth and has lost small-scale heterogeneities. All of the realizations are similar to the
15
ACCEPTED MANUSCRIPT
461 462 463 464 465 466 467 468 469 470
473 474 475 476 477
AC C
EP
TE D
M AN U
471 472
RI PT
459 460
smooth model on the large-scale and contain small-scale heterogeneities. They are statistically similar but are always different in their details. Geological information, such as horizons and faults, constrain the inverted AI, and the difference of the AI in the different horizons is recovered well. The variance of model suggests that different horizons have different uncertainties since the heterogeneities are different. It is obvious that the variance at the well locations is small and becomes bigger gradually as the node is farther. Some areas have high variance which suggests that it may contain thinner beds than seismic resolution. Since we cannot invert thin beds exactly, and these areas have high uncertainty. However, the vertical heterogeneity of each realization is shown to be similar to the well. We can observe that there are two areas with low AI values in the Green formation, which can be identified as oil sand body. Time slices through the Green formation (Figure 14) can help us to determine the hydrocarbon zones by low AI values. We also give a quantitative measurement of the results. The value of the objective function decreases by approximately 78%, and the correlation coefficient of the synthetic and real seismic data is 0.76. Therefore, this can be considered a satisfactory result for a real case.
SC
457 458
Figure 13. Crossline section of the inverted AI model along the red line in Figure 9: (left top) the mean of 100 realizations, (left bottom) variance of 100 realizations, (right) 3 of 100 different realizations, and the triangle in each section stands for the well location.
16
Figure 14. Time slices of the inverted AI model at 36ms below the top of the Green formation. (top left) the mean of 100 realizations, (top right) variance of 100 realizations, (bottom) 3 of 100 different realizations.
M AN U
478 479 480 481
SC
RI PT
ACCEPTED MANUSCRIPT
482
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507
TE D
EP
489 490 491
4 CONCLUSION We presented a new stochastic seismic inversion method based on an improved local GDM that combines seismic data, well log data, geological information and their spatial variability. The use of indicator cokriging allows us to incorporate the geological information and well data. FFT-MA and p-field simulation could generate fine-scaled conditional realizations consistent with spatial correlation information. Two improved strategies make local GDM successful use in seismic inversion. In addition, our proposed method can be used in non-Gaussian reservoir model and can be easily made parallel. The applications in the synthetic and real cases show that our method can generate reasonable AI models and give uncertainty estimations. It also has a high efficiency accelerated by the GPU. In future research, non-stationary natural phenomena should be studied since they are common in geology. Another issue is how to determine the number of parts that we should divide the trace into in the vertical direction. This may be influenced by the wavelet length, the autocorrelation length in the vertical direction and the number of samples in a trace.
AC C
483 484 485 486 487 488
Acknowledgments We gratefully acknowledge the editor and reviewers for the helpful comments and suggestions. REFERENCES Azevedo, L., Nunes, R., Soares, A., Mundin, E. C., Neto, G. S., 2015. Integration of well data into geostatistical seismic amplitude variation with angle inversion for facies estimation. Geophysics 80(6), M113-M128. Bortoli, L. J., Alabert, F. A., Haas, A., and Journel, A. G., 1993. Constraining stochastic images to seismic data. Geostatistics Tróia ’92(1), 325–337.
17
ACCEPTED MANUSCRIPT
519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551
RI PT
517 518
SC
515 516
M AN U
513 514
TE D
512
EP
510 511
Bosch, M., Mukerji, T., Gonzalez, E. F., 2010. Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: a review. Geophysics 75(5), 75A165-75A176. Buland, A., Kolbjørnsen, O., and Omre, H., 2003. Rapid spatially coupled AVO inversion in the Fourier domain. Geophysics 68, 824–836. doi: 10.1190/1.1581035. Buland, A., and Omre, H., 2003. Bayesian linearized AVO inversion. Geophysics 68, 185–198. doi: 10.1190/1.1543206 Cooke, D. A., and Schneider, W. A., 1983. Generalized linear inversion of reflection seismic data. Geophysics 48, 665-676. Deutsch, C. V., Journel, A. G., 1998. Gslib: Geostatistical Software Library and User's Guide: Oxford University Press. Doyen, P., 2007. Seismic Reservoir Characterization: EAGE. Eidsvik, J., Avseth, P., Omre, H., Mukerji, T., Mavko, G., 2004. Stochastic reservoir characterization using prestack seismic data. Geophysics 69, 978–993. doi: 10.1190/1.1778241. Ferreirinha, T., Nunes, R., Azevedo, L., Soares, A., Pratas, F., Tomas, P., Roma, N., 2015. Acceleration of stochastic seismic inversion in openCL-based heterogeneous platforms. Computers & Geosciences 78, 26-36. Francis, A. M., 2006. Understanding stochastic inversion. part 1: First Break 24, 69-77. Froidevaux R., 1993. Probability Field Simulation. Geostatistics Tróia ’92, 73-83. Gervais-Couplet, V., Roggero, F., Feraille, M.,Le Ravalec, M., Seiler, A., 2010. Joint History Matching of Production and 4D-Seismic Related Data for a North Sea Field Case. In: SPE Annual Technical Conference and Exhibition, 19-22 September, Florence, Italy González, E. F., Mukerji, T. and Mavko, G., 2008. Seismic inversion combining rock physics and multiple-point geostatistics. Geophysics 73(1), R11–R21. doi: 10.1190/1.2803748. Grana, D., and Della Rossa, E., 2010. Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion. Geophysics 75(3), O21–O37. doi: 10.1190/1.3386676. Grana, D., Mukerji, T., Dvorkin, J., Mavko, G., 2012. Stochastic inversion of facies from seismic data based on sequential simulations and probability perturbation method. Geophysics 77(4), 53. Grijalba-Cuenca, A., Torres-Verdín, C., and Debeye, H. W. J., 2000. Geostatistical inversion of 3D seismic data to extrapolate wireline petrophysical variables laterally away from the well. In: SPE Annual Technical Conference and Exhibition. Haas, A., and Dubrule, O., 1994. Geostatistical inversion-a sequential method of stochastic reservoir modelling constrained by seismic data. First break 12, 561-569. Hu, L. Y., 2000, Gradual deformation and iterative calibration of Gaussian-related stochastic models. Mathematical Geology 32, 87-108. Hu, L. Y., and Le Ravalec-Dupin, M., 2004. An improved gradual deformation method for reconciling random and gradient searches in stochastic optimizations. Mathematical Geology 36, 703-719. Journel, A. G., 2002. Combining knowledge from diverse sources: an alternative to traditional data independence hypotheses. Mathematical Geosciences 34(5), 573-596. Le Gallo, Y., and Le Ravalec-Dupin, M., 2000. History matching geostatistical reservoir
AC C
508 509
18
ACCEPTED MANUSCRIPT
563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595
RI PT
561 562
SC
559 560
M AN U
557 558
TE D
556
EP
554 555
models with gradual deformation method. In: SPE Annual Technical Conference and Exhibition. Le Ravalec, M., Hu, L. Y., Noetinger, B., 1999. Stochastic reservoir modeling constrained to dynamic data: Local calibration and inference of the structural parameters. In: SPE Annual Technical Conference and Exhibition. Le Ravalec, M., Noetinger, B., and Hu, L. Y., 2000. The FFT moving average (FFT-MA) generator: an efficient numerical method for generating and conditioning Gaussian simulations. Mathematical Geology 32, 701-723. Le Ravalec, M., Noetinger, B., 2002. Optimization with the gradual deformation method. Mathematical Geology 34 (2), 125-142. Le Ravalec-Dupin, M., Coureaud, B., Nicolas, L., Roggero, F., 2004. Conditioning an underground gas storage site to well pressures. Oil & gas science and technology 59 (6), 611-624. Le Ravalec, M., Mouche, E., 2012. Calibrating transmissivities from piezometric heads with the gradual deformation method: An application to the Culebra Dolomite unit at the Waste Isolation Pilot Plant (WIPP), New Mexico, USA. Journal of Hydrology 472, 1-13. Marteau, B., Ding, D. Y., and Dumas, L., 2014. A generalization of the local gradual deformation method using domain parameterization. Computers & Geosciences 72, 233-243. Mukerji, T., Jorstad, A., Avseth, P., Mavko, G., Granli, J. R., 2001. Mapping lithofacies and pore-fluid probabilities in a North Sea reservoir: Seismic inversions and statistical rock physics. Geophysics 66, 988–1001. doi: 10.1190/1.1487078. Nunes, R., Almeida, J. A., 2010. Parallelization of sequential Gaussian, indicator and direct simulation algorithms. Computers & Geosciences 36(8), 1042-1052. Nunes, R., Soares, A., Schwedersky, G., Dillon, L., Guerreiro, L., Caetano, H., Maciel, C., Leon, F., 2012. Geostatistical inversion of prestack seismic data. In: Ninth International Geostatistics Congress, Oslo, Norway, 1-8. Nunes, R., Soares, A., Azevedo, L., Pereira, P., 2016. Geostatistical seismic inversion with direct sequential simulation and co-simulation with multi-local distribution functions. Mathematical Geosciences, 1-19. Oldenburg, D. W., Scheuer, T., Levy, S., 1983, Recovery of the acoustic impedance from reflection seismograms, Geophysics 48, 1318-1337. Oliver, D. S., 1995. Moving averages for Gaussian simulation in two and three dimensions. Mathematical Geology 27, 939-960. Pourpak, H., Bourbiaux, B., Roggero, F., Delay, F, 2009. An integrated method for calibrating a heterogeneous/fractured reservoir model from wellbore flow measurements: Case study. SPE Reservoir Evaluation & Engineer 12, 433-445. Sanders, J., Kandrot, E., 2010. CUDA by Example: An introduction to general-purpose GPU programming. Addison-Wesley Professional. Soares, A., 2001. Direct sequential simulation and cosimulation. Mathematical Geosciences 33(8), 911-926. Soares, A., Diet, J. D., Guerreiro, L., 2007. Stochastic inversion with a global perturbation method. In: EAGE Conference on Petroleum Geostatistics. Srivastava, R. M., 1992. Reservoir characterization with probability field simulation. SPE
AC C
552 553
19
ACCEPTED MANUSCRIPT
598 599 600 601 602 603 604
Annual Technical Conference and Exhibition, 4-7 October, Washington, D.C. Tahmasebi, P., Sahimi, M., Mariethoz, G., Hezarkhani, A., 2012. Accelerating geostatistical simulations using graphics processing units (GPU). Computers & Geosciences 46(3), 51-59. Tarantola, A. 2005. Inverse problem theory and methods for model parameter estimation. Society for Industrial and Applied Mathematics. Vargas, H. S., Caetano, H., Filipe, M., 2007. Parallelization of sequential simulation procedures. In: EAGE Conference on Petroleum Geostatistics.
RI PT
596 597
605 606
AC C
EP
TE D
M AN U
SC
607
ACCEPTED MANUSCRIPT Research highlights A seismic inversion method based on the local gradual deformation method is proposed. The method can incorporate seismic data, well data, geology and spatial correlations.
RI PT
Probability field simulation and FFT moving average is combined to generate models.
AC C
EP
TE D
M AN U
SC
Two improved strategies of gradual deformation method are proposed.