Brian Russell , Dan Hampson, and Todor Todorov Hampson-Russell Software Calgary, Alberta.

1. Brian Russell, Dan Hampson, and Todor TodorovHampson-Russell SoftwareCalgary, Alberta. Combining geostatistics and multiattribute transforms – A channel sand case study

2. Introduction • In this talk, I will look at a new approach to integrating well log and seismic data, which involves post-stack inversion, geostatistics, and multiattribute transforms. • This method will be applied to data slices extracted from multiple 3D volumes. • I will illustrate this approachusing the Blackfoot dataset from Alberta, Canada.

3. Full Blackfoot survey Portion of data used in the present talk. N 09-08 This map is an amplitude slice of the complete Blackfoot survey, averaged below the channel sand top, showing the area used in this talk, and well 09-08 for reference.

4. Base map from portion of Blackfoot survey N Xline 18 This map shows the location of the 12 wells in the area within the seismic base map, and Xline 18. Note that we have rotated the data slice by 900 for better visual presentation.

5. Seismic line from volume This figure shows Xline 18 from the seismic volume, with sonic logs from two intersecting wells spliced in.

6. Line from inverted seismic volume This figure shows Xline 18 from the inverted volume, with sonic logs from two intersecting wells spliced in. The color key indicates impedance.

7. Correlating the logs with the seismic data This figure shows the correlation of well 14-09 with the seismic data, where the synthetic trace is in blue the and the seismic trace is in red. The sonic and porosity logs are on the right. The top and base of sand are also shown.

8. Introduction to geostatistics • Map-based geostatistics involves producing three types of maps: • Optimal maps (Best Linear Unbiased Estimates) from sparse well data (kriging). • Maps that incorporate both sparse well data and a secondary seismic attribute (cokriging and kriging with external drift (KED)). • Simulations of a range of possible scenarios, which will not be discussed today. • All these maps are based on the variogram, which will be discussed shortly. • A good place to start is with a crossplot between the primary (well) and secondary (seismic) datasets.

9. Wells showing average porosity This map shows the distribution of the wells with the average porosity over the zone of interest in color.

10. Acoustic impedance slice This map shows the arithmetic average of the acoustic impedance over a 10 ms window below the channel top event.

11. Initial crossplot This is the crossplot between the well porosities and acoustic impedance values. The red line is the regression fit, and the correlation is –0.65.

12. Regression applied to inversion slice This map shows the application of a the regression line from the previous slide to the inversion slice. Unfortunately, the wells do not tie!

13. Kriging and the spatial covariance • Kriging involves creating map values (f0) that are a weighted sum of the input well values (f1through fN): • The weights are determined by solving a set of linear equations involving the spatial covariance, C(hij), where hij is the distance between the ith and jth porosity values:

14. Well to well variogram The spatial covariance is derived from the spatial variogram, g(h), by the equation: C(h) = C(0) -g(h). The figure above shows the variogram from the 12 wells on the map.

15. Seismic to seismic variogram This figure shows the seismic to seismic variogram calculated from the inversion slice shown earlier. Using the Markov-Bayes linear assumption, we can use this variogram for all of our maps.

16. Kriged result This map shows the result of applying kriging to the 12 wells within the 3D seismic survey, using the scaled seismic to seismic variogram.

17. Estimated kriging error To see the error associated with kriging, we often display the error variance. But this is simply the “theoretical” error, and will go to zero as the variance between the input values goes to zero.

18. Kriging cross-validation error A more meaningful error is the cross-validation error, found by leaving each well out in turn, and predicting its value using kriging. The validation error for the kriged map is shown above, displaying the absolute error at each well in % porosity. The standard deviation is 3.25%.

19. Cokriging and KED • Cokriging is an extension of kriging, in which the original equation is extended to include a weighted sum of secondary (seismic) values. • In collocated cokriging, only the seismic value at the output location is used (i.e., we have 3D seismic data). • In kriging with external drift (or KED for short) the drift, or trend, from the seismic data is added to the kriged result. • As with kriging, an error variance plot can be made for both the cokriged and KED maps. • The results of cokriging and KED applied to our dataset are shown in the next four slides.

20. KED result This map shows the initial KED result, using the inversion slice as the secondary variable. The standard deviation of the validation error in this map is 3.10 %, which is less than for kriging. Note the “imprint” of the kriged map.

21. Cokriged result This map shows the initial collocated cokriging result, using the inversion slice as the secondary variable. Again, note the imprint of the kriged map, which is more dominant than for cokriging. However, standard deviation of the validation error is 2.91%, which is better than for kriging or KED.

22. The multiattribute transform • The multiattribute transform involves combining multiple attributes to predict a reservoir parameter such as porosity, by training at the well locations. • Initially, this method was applied to seismic volumes, using each well logs sample over a window. Here, we will use maps and averaged well log values. • Two transforms are available: multilinear regression and the probabilistic neural network (PNN). • For the map case, we felt that there were not enough points to train adequately PNN, so this method will not be discussed here. • In the first stage we train the multiattribute transform. • In the second stage we apply the results of the training to the multiple map attributes.

23. The multiattribute map transform Attribute map 1 X Y Attribute map 2 Attribute map M This figure shows the multiattribute map transform approach in schematic form, where we need to compute the weights wi which, when multiplied by the attribute values, will produce the log value.

24. Multilinear regression • In the multiattribute transform using multilinear regression, we compute M+1 weights such that the log value L(X,Y) at a particular map value is a weighted sum of M attributes: • If we have N input values for each attribute, the forward problem can be written L=AW, where A is an N x M matrix of observations, and the solution is W = (ATA)-1ATL, which involves the covariance matrix shown below:

25. Attribute slices • We have already seen one of the slices that will be used in our multiattribute transform: the impedance slice. • The next three slides will show the other slices used. Each attribute (except trace length) was derived from the seismic volume by taking an RMS average over a 10 msec window below the picked top of sand. • The following attribute slices were extracted: • Seismic amplitude • Amplitude envelope • Instantaneous phase • Instantaneous frequency • Integrated seismic trace • Trace length – the total length of the trace over the window

26. Seismic amplitude slices (a) (b) The map in (a) shows the RMS average of the seismic amplitude over a 10 ms window below the channel top event, whereas the map in (b) shows the RMS average of the amplitude envelope over the same window. Notice that the two slices are very close in appearance.

27. Instantaneous phase and frequency slices (a) (b) The map in (a) shows the RMS average of the instantaneous phase over a 10 ms window below the channel top event, whereas the map in (b) shows the RMS average of the instantaneous frequency over the same window.

28. Integrated trace and trace length slices (a) (b) The map in (a) shows the RMS average of the integrated trace over a 10 ms window below the channel top event, whereas the map in (b) shows the total trace length over the same window.

29. Computing the multilinear regression result • We are now ready to compute the multilinear regression map, which will be some linear combination of the previous maps. • The next slide shows the single correlation coefficients for each of the attribute slices. • We then compute the best combination of attributes using a technique called step-wise regression. • Finally, we decide which attributes are significant using a validation technique in which the target well is left out in jackknife fashion. • This is shown in the subsequent five slides.

30. Correlation coefficients for all the slices This table shows the correlation coefficients between the well porosity values and all of the attribute slices, sorted by decreasing correlation coefficient.

31. Validation error plot This is the validation error for the 5 attributes used in the multiattribute process. The red line leaves out the target well and shows that only the first 3 attributes should be used.

32. Validation error and weights Only use the first three attributes based on validation error. Weights The table at the top shows the numerical values from the previous slide and also shows that the best non-linear fit is between the square root of porosity and the inverse attributes. The bottom table shows the weights.

33. Multilinear regression result This map is the result of applying the multilinear regression weights, shown in the previous slide, to the attributes. Note that the result is in pseudo-porosity.

34. New crossplot This is the new crossplot between the well porosity and the pseudo-porosity from the multiattribute transform. Note that the correlation coefficient has gone up to 0.81.

35. Combining multilinear regression with geostatistics • We are now at the final stage of the process, and will combine the multilinear regression result with the well values using geostatistics. • That is, the multiattribute transform will replace the inversion slice as the secondary variable. • The first step is to re-compute the seismic to seismic variogram. • We will then compute new KED and cokriging results. • Finally, we will do a statistical analysis of the results.

36. New seismic variogram This is the seismic variogram used in the final KED and cokriging processes. (Note for the geostatistician: this is an exponential fit, rather than the spherical fit used in the earlier variograms.)

37. KED with the multiattribute transform This map shows the result of applying the KED process using the multiattribute transform as the secondary variable. The standard deviation of the validation error is 2.45%, which is the better than the previous maps.

38. Cokriging with the multiattribute transform This map shows the result of applying the collocated cokriging process using the multiattribute transform as the secondary variable. The standard deviation of the validation error is now 2.3%, which is the best of all the maps.

39. KED and cokriged maps (a) (c) Std. Dev. = 3.1 % Std. Dev. = 2.45 % (b) (d) Std. Dev. = 2.91 % Std. Dev. = 2.33 % This map shows a summary of the four key maps, where (a) is KED with inversion, (b) is cokriging with inversion, (c) is KED with multiattributes, and (d) is cokriging with multiattributes.

40. Kriging through cokriging with multiattributes (a) Std. Dev. = 3.25 % (b) Std. Dev. = 2.91 % (c) Note the increase in geological information as we move from (a) the kriged map with wells alone to (b) cokriging with inversion, and finally to (c) cokriging with the multiattribute transform. Note that all three maps match the wells. Std. Dev. = 2.33 %

41. Conclusions • In this paper, we have combined geostatistics with multiattribute map analysis. • Traditional geostatistics uses a single secondary variable for cokriging and KED. We found that, in this case, using inversion as the secondary attribute produced results that showed a strong “imprint” from the wells. • By using a multiattribute transform based on multilinear regression, we were able to get a much better fit between the wells and the final pseudo-porosity map. • Following up with a second pass of geostatistics, the new KED and cokriged maps were much more feasible from and exploration point of view. • Statistically, the standard deviation of the errors at the wells showed a constant decrease, and were best using cokriging with the multiattribute transform.

42. Acknowledgements • We wish to thank the University of Calgary CREWES consortium and their sponsors, especially PanCanadian Petroleum, for permission to use the Blackfoot dataset. • We also want to thank Jim Schuelke and John Quirein for originally suggesting the multiattribute approach.