Category Archives: Tutorial

Machine learning in geoscience with scikit-learn. Part 1: checking, tidying, and analyzing the dataset

Posted on by
2

The idea behind this series of articles is to show how to predict P-wave velocity, as measured by a geophysical well log (the sonic), from a suite of other logs: density, gamma ray, and neutron, and also depth, using Machine Learning.

The log suite is from the same well that Alessandro Amato del Monte used in the Seismic Petrophysics Notebook accompanying his Geophysical tutorial article on The Leading Edge.

I will explore different Machine Learning methods from the scikit-learn Python library and compare their performances.

To wet your appetites, here’s an example of P-wave velocity, Vp, predicted using a cross-validated linear model, which will be the benchmark for the performance of other models, such as SVM and Random Forest:

multilinear

In the first notebook, which is already available on GitHub here, I show how to use the Pandas and Seaborn Python libraries to import the data, check it, clean it up, and visualize to explore relationships between the variables. For example, shown below is a heatmap with the pairwise Spearman correlation coefficient between the variables (logs):

heatmap

Stay tuned for the next post / notebook!

PS: I am very excited by the kick-off of the Geophysical Tutorial (The Leading Edge) Machine Learning Contest 2016. Check it out here!

Machine learning in geoscience and planetary science with scikit-learn: series outline

Posted on by
Reply
  • Machine learning in geoscience with scikit-learn. Part 3: the SEG ML contest
  • Machine Learning in Geoscience with Scikit-learn. Part 4: TBE

sketch2model – sketch image enhancements

Posted on by
7

This is the second post of in a series of collaborative articles about sketch2model, a project from the 2015 Calgary Geoscience Hackathon organized by Agile Geoscience.

The first post was written by Elwyn Galloway and published on both his Scibbatical blog and here on MyCarta. In that article Elwyn mentioned the need for an adaptive image conditioning workflow for binarization of photos with geological sketches in images. Binarization is the process of converting a natural image to a binary image (check this simple but awesome interactive demonstration of binarization), which in our case is necessary to separate the sketch from the background.

The following is a demonstration of the preliminary image processing operations applied to the input photo when sketch2model is run. The full code listing of this demonstration is available as a Jupyter notebook on GitHub. Also in GitHub you can find a Jupyter Notebook with the fully documented version of sketch2model.

First we import one of the photos with sketches and convert it to a grayscale image.

im = io.imread('paper_breaks.png')
im = color.rgb2gray(im[0:-1:2,0:-1:2])

paper_breaks_post

Next we enhance the grayscale image with a couple of cascaded processes. But before we do that, let’s graph the intensity values across the image to understand the degree of contrast between sketch edges and background, which ultimately will determine our success in separating them from it. We show this in the figure below, on the left, for one column of pixels (y direction). The black line across the input image on the right shows the location of the column selected. It is fairly obvious from the plot on the left that the intensity of the background is not uniform, due to variable light conditions when the photo was taken, and towards the right (e.g. bottom of the photo) it gets closer to that of the edges. In some images it might even become less than the intensity of the edge. This highlights the need for (preemptively) applying the enhancements illustrated in the remainder of the post.

section_raw

The first enhancement is called compressor, or limiter. I read many years ago that it is used in electronics to find hard edges in data: the idea is to square each element in the data (image, or other type of data), smooth the result (enough to remove high frequency variations but not so much as to eliminate variability), take the square root, and finally divide each element in the input by the square root result.

I experimented with this method (at the time using Matlab and its Image Processing Toolbox) using the same gravity dataset from my 2015 geophysical tutorial on The Leading Edge (see the post Mapping and validating geophysical lineaments with Python). An example of one such experiments is shown in the figure below where: the top left map is the Bouguer data; the centre top map is the squared data; the top right is the result of a Gaussian blur; the bottom left the result of square root, and centre right is the final output, where the hardest edges in the original data have been enhanced.

Sketch2model_compressor_campi

The most important parameter in this process is the choice of the smoothing or blur; using a Gaussian kernel of different size more subtle edges are enhanced, as seen in the bottom right map (these are perhaps acquisition-related gridding artifacts).

In our sketch2model implementation the size of the Gaussian kernel is hardcoded; it was chosen following trial and error on multiple photos of sketches and yielded optimal results in the greatest majority of them. We were planning to have the kernel size depend on the size of the input image, but left the implementation to our ‘future work’ list.’

Here’s the compressor code from sketch2model:

# compressor or limiter (electronics): find hard edges in data with long 
# wavelength variations in amplitude
# step 1: square each element in the image to obtain the power function
sqr = im**2
# step 2: gaussian of squared image
flt2 = sp.ndimage.filters.gaussian_filter(sqr,21)
# step 3: divide the intensity of each original pixel by the square root 
# of the smoothed square
cmprs= im/(np.sqrt(flt2))

and a plot of the result (same column of pixels as in the previous one):

section_compressed

From the plot above we see that now the background intensity is uniform and the contrast has been improved. We can maximize it with contrast stretching, as below:

# contrast stretching
p2, p98 = np.percentile(cmprs, (2, 98))
rescale = exposure.rescale_intensity(cmprs, in_range=(p2, p98))

section_stretched

We now have ideal contrast between edges and background, and can get a binary image with the desired sketch edges using a scalar threshold:

# binarize image with scalar threshold
binary = ~(color.rgb2gray(rescale) > 0.5)

section_binary

Bingo!

Simulating seismic surveys using King Tut’s CAT scan

Posted on by
2

The remote sensing used to study the human body is very similar to the remote sensing used to study the subsurface. Apart from a scaling factor (due to the different frequencies of the signals used) the only major difference between the two methods of investigation is in that radiologists and doctors looking at an x-ray, ultrasound, or CAT scan image know what to look for in those images, as bones, tissues, and anomalies, have known characteristics, whereas the subsurface is always to a large extent unknown.

In this short visual post I am going to use a CAT scan of King Tut’s skull to explore the effect on the image quality of progressive decimation of the data followed by upsampling it back to the initial size. I will also look at the effect of these manipulations on the results of edge detection.

With this I want to simulate the progressive reduction in imaging quality that happens when going from high density 3D seismic acquisition to medium density 3D seismic to high quality, but sparse 2D seismic lines.

Here’s the input image in Figure 1.

tut20bone20frag

Figure 1. CAT scan of King Tut’s skull – Supreme Council of Antiquities. guardians.net/hawass/press_release_tutankhamun_ct_scan_results.htm

 

In Figure 2 I am showing the image after import into a Jupiter Notebook and conversion to grayscale, and the result of edge detection using the Sobel filter. Notice the excellent quality of the edge detection result.

ground_Tut

Figure 2. Original image, or ground truth for the experiment,  and edge detection result.

To simulate a high-resolution 3D seismic acquisition I decimated the original image by a factor of 4 in both directions. The resulting image (no interpolation) is, shown in Figure 3, is of good quality, and so is the edge detection result.

highres_3D

Figure 3. Simulated high-resolution 3D survey and edge detection result.

The image in Figure 4 results from a further decimation by a factor of 2 of the image in Figure 3, then interpolation to upsample to the same size as the image in Figure 4. The image and the edge detection are still of fair quality overall, but some of the smaller features have either disappeared, merged, or faded.

Figure 4. Simulated medium resolution 3D survey and edge detection result.

Figure 4. Simulated medium resolution 3D survey and edge detection result.

Now look at Figure 5: this is the equivalent of a high quality (in one direction) 2D dataset. Although we can still guess at what this represents, I would argue this is a result of our a priori knowledge of what it is supposed to represent – a human skull; and yet I don’t think anybody would want their doctor to make a diagnosis  based on this image.

Figure 5. Simulated set of very high-resolution 2D lines.

Figure 5. Simulated set of very high-resolution 2D lines.

The image in Figure 6 results from 2D interpolation (my intention is to simulate the result we would get by gridding 2D data to get a continuous image. We can now definitely interpret this as a skull, but the edge detection result is very unsatisfactory.

Figure 6. Simulated interpolation of 2D lines.

Figure 6. Simulated interpolation of 2D lines.

In  future post we will explore the effects of adding periodic noise (similar to seismic acquisition footprint) on these images and on the edge detection results. I will also show you how to remove it using 2D FFT filters, as promised (now more than a year ago) in my post Moiré Patterns.

If you would like to play with the code, get the Jupiter Notebook here.

A study of Ricker wavelets in MS Excel

Posted on by
Reply

There, I said it, and aloud: ‘EXCEL’!

I got this idea of making a modern (no, I am NOT kidding) educational tool to interactively construct and study Ricker wavelets after reading William Ashcroft’s A Petroleum Geologist’s Guide to Seismic Reflection. The book is a gem, in my opinion (I discovered it thanks to a blog post by Matt Hall a couple of years ago) and comes with some interesting tools on a DVD, which are unfortunately a bit outdated. So I went back to my old Geophysics Bible (Yilmaz’ book Seismic Data Analysis), and mashed a few ideas up; you can Download the Excel file here. Please feel free to use it, peruse it, and abuse it, in your classes, labs, nightmares, etcetera…

The tool is fairly simple. In Sheet 1 the user enters the dominant frequency of the desired Ricker wavelet, as shown in the middle of Figure 1. From that informatin the wavelet is constructed using the equation A = g^2 * 1/exp g^2 where g is the ration between frequency f  (in increments of 5 Hz up to an arbitrary 125 Hz – but this could be easily changed!) and the dominant frequency f1 just entered.  frequencies. The frequency spectrum of the wavelet is shown as a graph.

sheet 1 input and amplitude spectrum 30 hz

Figure 1

In Figure 2 I show the same Sheet, but for a wavelet of dominant frequency equal to 50 Hz.

sheet 1 input and amplitude spectrum 50 hz

Figure 2

A number of ancillary graphs, shown in Figure 3, display individual building blocks of the formula A = g^2 * 1/exp g^2.

Figure 3

Figure 3

In Sheet 2 the user can view a plot of the wavelet in the time domain,  add Constant phase shifts and Linear phase shifts, and experiment with different combinations of them, as shown in Figures 4-7, below.

sheet 2 zero phase

Figure 4

sheet 2 linear phase shift 180

Figure 5

sheet 2 constant phase shift 90

Figure 6

sheet 2 constant 90 and linear 180 phase shifts

Figure 7

Finally, Sheet 3 displays a plot of the wavelet, and of the individual frequency components that make it up. Have fun!!!

sheet 3 component sinusoids

Figure 8

References and further playing

A Petroleum Geologist’s Guide to Seismic Reflection, William Ashcroft, 2011, Wiley.

Seismic Data Analysis, Öz Yilmaz, 2001, Society of Exploration Geophysicists.

To plot a wavelet [in Python]. blog post by Evan Bianco.

Seismic terrain displays

Posted on by
8

Introduction

A couple of years ago I stumbled in a great 2001 paper by Beyer [1] on The Leading Edge. Being interested in visualization techniques I was drawn by the display in Figure 1 (which is a low resolution copy from Figure 2 in the paper). But what really amazed me was the suggestion that a display like this could be created in a few minutes, without doing any interpretation, by just manipulating instantaneous phase slices. With the only condition of having data of fair quality, this promised to be an awesome reconnaissance tool.

Beyer_Fig2_lowres

Figure 1 – Copyright SEG

The theory

The idea that instantaneous phase is a great attribute for interpretation has been around for a long time. There is, for example, a 1989 Exploration Geophysics paper by Duff and Mason [2]. These two authors argue that amplitude time slices are a suboptimal choice, and that instantaneous phase slices should be preferred. They give three reasons:

1) on amplitude time slices only relatively strong events remain above the bias level after gain and scaling. Weak events are submerged below the bias and remain unmappable.  Although it is the topic of a future post, it is worth mentioning I think this effect is exacerbated by the common but unfortunate choice of a divergent color palette with white in the middle. White is so bright (I call it white hole) that even more low amplitude events become indiscernible.

2) discrete boundaries corresponding to unique positions on the wavelet are displayed on instantaneous phase slices – this intra wavelet detail is lost on amplitude slices.

After-Duff-Mason 2

Figure 2 – after Duff and Mason, Figure 2

3) instantaneous time slices give DIRECTLY the sense of time dip for dipping events. In Figure 2 I show 2 parallel dipping reflectors, represented by 5 (non consecutive) traces, and 2 (non-consecutive) instantaneous phase slices (at arbitrary t1 an t2). I marked 5 discrete phase events for the top dipping reflector. The sense of time dip is given (with appropriate color palette) by the sense of color transition. Conversely, this intra wavelet detail would be lost on the amplitude time slices, with amplitudes between the black center trace and the red traces, and amplitudes between the red traces and the green traces lost within a single broad zone. The difference is probably not as dramatic nowadays with the increase in dynamic ranges available, but using instantaneous phase slices still remains advantageous for detailed mapping.

Beyer’s seismic terrain is just a natural extension of the instantaneous time slice as ( quoted from [1]):”… it then follows that the instantaneous phase (-180 deg to +180 deg) can simply be rescaled to the wavelength in ms of pseudoseismic two-way time… Seismic terrain can be thought of as a type of instantaneous wavelength generated from instantaneous phase along a time slice”. With reference to the top dipping reflector in Figure 3, the method allows generating converting the brown phase segment to the dipping blue segment (and similarly the yellow phase segment to the dipping green segment for bottom dipping event).

After-Beyer

Figure 3

The practice – Petrel

Let’s see how we can create a terrain display similar to that in Figure 1 using Petrel.

Raw data

The process starts with migrated seismic data, from which we need generate both the phase and frequency component to get us the instantaneous wavelength.

For this tutorial I use a public seismic dataset (BPA9901) available on the Norwegian Public Data Portal. In Figure 4 below I am showing an amplitude time slice (above the Chalk) from the migrated seismic volume.

time_slice

Figure 4

Step 1 – generate phase component

The first step is to generate an instantaneous phase attribute volume. This is found in the volume attributes. In Figure 5 below I am showing the instantaneous phase slice corresponding to the amplitude time slice of Figure 4.

*** N.B. *** If significant regional dips are observed in the seismic data, care should be taken in some cases  it may be beneficial (please see comment section) to remove them through flattening prior to the terrain generation.

Figure 5

Figure 5

Step 2 – generate frequency component

This is the trickiest part. In theory to get the instantaneous wavelength we would have to calculate the instantaneous frequency and divide the  instantaneous frequency attribute can be very noisy and can have spurious values in areas of low amplitude in the input data. A good practical alternative is to measure a single value of the dominant wavelet period T in an area of relatively flat reflections near the zone of interest as I am showing in Figure 6.

For the more avid readers, this is all explained quite nicely in Beyer (quoted from [1]): “Complex trace relationships dictate that the wavelength is the phase component divided by the frequency component. Thus one may be compelled to derive the seismic terrain by dividing the extracted instantaneous phase by the  extracted instantaneous frequency (carefully applying unit conversion of 1000 ms/ 360 deg or 2.78). However extracted instantaneous frequency tends to include spurious values in low amplitudes (approaching infinity according to literature and practice which correspond to poor data quality zones. Instantaneous frequency or even averaged instantaneous frequency renders the seismic terrain noisy, unrealistic, and misleading. Years of extensive use have shown that  single value of visually estimated dominant wavelet period (i.e. cycles per second) produces a very high-quality seismic terrain that closely fits the seismic events over wide areas”.

period

Figure 6

Step 3 – generate the instantaneous wavelength (seismic terrain)

Having estimated the dominant wavelet period T (in ms) , I can now use it to generate the instantaneous wavelength. We are essentially converting the data from the range [-180 180] deg to the range [-T/2 T/2] ms.

In Petrel I do it in the calculator with a formula of the type:

terrain=(((instantaneous phase +180)*T/2)/180)

The actual formula used is shown on the top row of Figure 7 below. You will notice that it isn’t exactly the same as the above formula. I added 1 to 180 to avoid division of zero values, e.g.:

(([-180 180]+180)*28/180) = ([0 180]*28/180) = ([0 5040]/180) % not good!

whereas:

(([-180 180]+181)*28/180) = ([1 181]*28/180) = ([28 5068]/180) % good!

Once the division is performed I subtract T/2 again.

Notice from Figure 7 that because we added 181 but divided by 180 there is a small adjustment to be made by hand. I get this small adjustment by double clicking on the output volume to get the statistics. In this case it is +-0.16 ms, so I run a second time the formula (bottom row, Figure 7), this time subtracting T/2 +0.16 instead of just T/2.

formula

Figure 7

Step 4 – display seismic terrain.

There are two options in Petrel to display the resulting seismic terrain volume:

Option 1 – display bump mapped terrain slices in 2D or 3D window

This is my preferred option for scanning up and down through the terrain slices. The bump mapping effect is done by double clicking on the terrain survey with a 2D or 3D window open and selected, in the Style tab>Intersection  tab. In Figure 8 I am showing the bump mapped terrain slice corresponding to the instantaneous phase time slice of Figure 5.

terrain_slice

Figure 8

Option 2 – display selected horizons of interest in 3D window

One may want to create a display such as the one in Figure 1, which for me is intended for a later stage, when integrating perhaps with extracted amplitude or attribute anomalies highlighting hydrocarbon presence.

As often in Petrel there are different ways of achieving the same result. This is how I do it. First, I create a flat surface, with a TWT value corresponding to the slice I am interested in (as in Figure 9, left panel). Then I extract and append to this surface the terrain values from that slice of interest (as in Figure 9, right panel).

Figure 9

Figure 9

Figure 10

Figure 10

Finally, in the calculations tab, I add a constant time shift corresponding to the time associated with the slice  of interest: notice the difference in Z value between the left panel in Figure 10 (before the calculation), and the right panel (after the calculation). It is also necessary to use the extracted value as visual vertical position as illustrated in Figure 11.

Figure 11

Figure 11

I am showing the result in Figure 12. This is the same terrain slice as in Figure 8.

terrain_slice extracted

Figure 12

Discussion

As a quick QC I am displaying in Figure 13 a vertical section (corresponding to the thick black line in Figure 12) from the input seismic data, with the extracted surface drawn as a thin black line.

The terrain deteriorates to the far left as we approach the edge of the survey, with fold decreasing and noise increasing, and there is a cycle skip towards the far right. But all in all I think  this is a very good result: it captures the faults well, and the whole process took less than 20 minutes with no picking.

arbitrary

Figure 13

Limitations

This method has one limitation: the maximum fault throws or stratigraphic relief (in milliseconds) that can be mapped is equal to the period T.

Acknowledgements

I wish to thank DONG Energy for agreeing to the publication of the seismic images, which were generated using company licensed Petrel.

** UPDATE **

A few readers asked clarifications on what the benefits and potential uses are of using this attribute. The short answer is that this is in fact a pseudo-horizon that tracks dipping  seismic events accurately within the range of the period of the dominant frequency period (as seen in Figure 13), which makes it an excellent reconnaissance tool.  A good quality first pass map can be made in minutes in areas where detailed mapping can take days.  An excellent example is Figure 16 in the original paper [1].  More details can be found in the last two paragraphs of the paper:  Reservoir-scale structural data from seismic terrain and Fast 3-D screening. 

References

[1] Beyer, L. (2001). – Rapid 3-D screening with seismic terrain: deepwater Gulf of Mexico examples. The Leading Edge, 20 (4), 386–395.

[2] Duff, B.A., and Mason, D.J. (1989) – The instantaneous-phase time slice: A crucial display for enhancing 3-D interpretation. Exploration Geophysics 20 (2) 213 – 217

Visualization tips for geoscientists – series outline

Posted on by
2

Visualization tips for geoscientists: Surfer

Visualization tips for geoscientists: Matlab

Visualization tips for geoscientists: Matlab, part II

Visualization tips for geoscientists: Matlab, part III

The rainbow is dead…long live the rainbow! – The rainbow is dead…long live the rainbow! – Perceptual palettes, part 4 – CIE Lab heated body

Posted on by
9

  In my last post I discussed the two main issues with the rainbow color palette from the point of view of human color vision, and concluded one of these issues is insurmountable.

But before I move to presenting alternative color palettes, let me give you one last example of how bad the rainbow is. It was sent to me by Antony Price, a member of the LinkedIn group Worldwide Geophysicists. Antony created a grayscale and a rainbow-colored version – using the same data range and number of intervals – of the satellite altimeter derived free-air gravity map of the world [1].  I am showing the two maps below.

Continue reading →