Seismic terrain displays

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 examplesThe 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 interpretationExploration Geophysics 20 (2) 213 – 217

NASA Worldview satellite image browser adopts MyCarta perceptual rainbow

I was thrilled this week to learn from Ryan Boller that his team at NASA’s ESDIS Project included MyCarta’s perceptual rainbow (the CubicYF) as one of the palettes for the Worldview satellite imagery browser.

If you’d like to try it, once on the viewer you can load an overlay and then you can choose from among several color palettes. The perceptual rainbow palette is listed here as “Rainbow 2”.

I am including below an example using the Land surface temperature for April 13 2013 from MODIS Aqua mission:

Land_surf_temp_130413

This is really exciting news as NASA’s adoption will increase the palette’s exposure and its chances of becoming more mainstream. This is also as close as I will ever get to realizing my childhood dream of becoming an astronaut. Thanks ESDIS, and thanks Ryan, on both accounts.

Colour: A quick guide to its use in informative graphics

A terrific tutorial on color / colour use in infographics from Darrell Wilkinson Geo-computing blog.

An introduction to color palettes for seismic amplitude – teaser

Introduction

In a future posts I will take a look at some of the color palettes used for seismic amplitude display, and discuss ways we can design more perceptual and more efficient ones.

For now, I would like to ask readers to look at two sets of seismic images and answer the survey questions in each section. Far from being exhaustive sets, these are meant as a teaser to get a conversation started and exchange opinions and preferences.

Stratigraphic interpretation

The seismic line below is inline 424 from the F3 dataset, offshore Netherlands from the Open Seismic Repository (licensed CC-BY-SA).

I generated an animation, played at 0.5 frames/second, where 8 different color palette are alternated in sequence.  Please click on the image to see a full resolution animation. I also generated a 0.25 frame/second version and a 1 frame/second version.

05

Fault interpretation

The images used to create the panel below are portions of seismic displays kindly provided by Steve Lynch of 3rd Science Solutions, generated using data released by PeruPetro. I am grateful to both.

faults_sm

Thanks to Matt Hall and Evan Bianco of Agile Geoscience for their suggestions.

Color palettes for seismic structure maps and attributes

I created three color palettes for structure maps (seismic horizons, elevation maps, etcetera) and seismic attributes. To read about the palettes please check these previous blog posts:

The rainbow is dead…long live the rainbow! – Part 5 – CIE Lab linear L* rainbow
The rainbow is dead series – Part 7 – Perceptual rainbow palette – the method
The rainbow is dead series – Part 7 – Perceptual rainbow palette – the goodies

The palettes are available as plain ASCII files and also formatted for a number of platforms and software products:

Geosoft
Hampson-Russell
Kingdom
Madagascar
Matlab
OpendTect
Petrel
Seisware
Surfer
Voxelgeo

Please download them from my Color Palettes page and follow instructions therein.

Enjoy!

linearlfb

Image courtesy of Sergey Fomel of the Madagascar Development blog

Perceptual rainbow palette – the goodies

Perceptual rainbow palette – Matlab function and ASCII files

In my last post I introduced cubeYF, my custom-made perceptual lightness rainbow palette. As promised there, I am sharing the palette  with today’s post. For the Matlab users, cube YF, along with the other palettes I introduced in the series, is part of the Matlab File Exchange submission Perceptually improved colormaps.

For the non-Matlab users, please download the cubeYF here (RGB, 256 samples). You may also be interested in cube1, which has a slightly superior visual hue contrast, due to the addition of a red-like color at the high lightness end but at the cost of a modest deviation from 100% perceptual. I used cube 1 in my Visualization tips for geoscientists series.

Perceptual rainbow palette – preformatted in various software formats

The palettes are also formatted for a number of platforms and software products: Geosoft, Hampson-Russell, SMT Kingdom, Landmark Decision Space Geoscience, Madagascar, OpendTect, Python/Matplotlib, Schlumberger Petrel, Seisware, Golden Software Surfer, Paradigm Voxelgeo. Please download them from my Color Palettes page and follow instructions therein.

Another example

In Comparing color palettes I used a map of South America [1] to compare a linear lightness palette to some common rainbow palettes using  grayscale as a perceptual benchmark. Below, I am doing the same for the cubeYF colormap.

South_America_maps_CubeYF_rainbow

Comparison of South America maps using, from left to right: ROYGBIV (from this post) , classic rainbow, cubeYF, and grayscale

Again, there is little doubt in my mind that cubeYF does a superior job compared to the other two rainbow palettes as it is free of artefacts [2] and more similar to grayscale  (with the additional benefit of color).

The ROYGBIV and cubeYF map have been included in Marek Kultys’ excellent tutorial Visual Alpha-Beta-Gamma: Rudiments of Visual Design for Data Explorers, recently published  on Parsons Journal for information mapping, Volume V, Issue 1.

An online palette testing tool

Both cubeYF and cube1 feature in the colormap evaluation tool by the Data Analysis and Assessment Center at the Engineer Research and Development Center. If you want to quickly evaluate a number of palettes, this is the right tool. The tool has a collection of many palettes, organized by categories, which can be used on 5 different test image, and examined in terms of RGB components and human perception. Below here is an example using cube YF.

hpc_terrain

An idea for a palette’s mood test

A few weeks ago, thanks to Matt Hall (@kwinkunks on twitter),  I discovered Colour monitor, a great online tool by Richard Weeler (@Zephyris on twitter). You supply an image; Colour monitor analyses its colors in terms of hue, saturation and luminance and produces a graphical representation of the image’s mood [3]. I thought, what a wonderful idea!

Then I wondered: what if I used this to tell me something about a color palette’s mood? The circular histogram of colors reminded me of the Harmonic templates [4] on the hue wheel from this paper And so I created fat colorbars using the three  palettes I used in the last post, saved them as images, and run the monitor with them. Here below are the results for Matlab jet, Industry Spectrum, and cubeYF. Looking at these palettes in terms of harmony I would say that jet is not very harmonic (too large a portion of the hue circle; the T template, which is the largest, spans 180 degrees), and that the spectrum is terrible.

CubeYF is also exceeding a bit 180 degrees, but looks very close to a T template rotated by 180 degrees (rotations are allowed). So perhaps I could trim it a bit? But to me it looks a lot nicer and gives me a vibe of really good mood, and reminds me of one of those beautiful central american headdresses, like Moctezuma’s crown.

jet-clrmp-mood

Jet mood

Spectrum-industry-clrmp-mood

Spectrum mood

cubeYF-clrmp-mood

cubeYF mood

Notes

[1] Created with data from the Global Land One-km Base Elevation Project at the National Geophysical Data Center.

[2] Looking at the intensity of the colorbars may help in the assessment: the third and fourth colorbars are very similar and both look perceptually linear, whereas the first and second do not.

[3] Quoted from Richard’s blog post: “… in the middle is a circular histogram of the colours (spectral shades) in the image, and gives an idea of how much of each colour there is. Up the left is a histogram of image brightness (lightness of colour), and up the right is a histogram of colour saturation (vibrancy)”.

[4] Quoted from the paper’s abstract: “Harmonic colors are sets of colors that are aesthetically pleasing in terms of human visual perception. If you are interested in this idea there is a set of slides and a video on the author’s website

Related posts

Perceptual rainbow palette – the method

With this post I would like to introduce my new, perceptually balanced rainbow color palette. I used the palette for the first time in How to assess a colourmap, an essay I wrote for 52 Things You Should Know About Geophysics, edited by Matt Hall and Evan Bianco of Agile Geoscience.

In my essay I started with the analysis of the spectrum color palette, the default  in some seismic interpretation softwares, using my Lightness L* profile plot and Great Pyramid of Giza test surface (see this post for background on the tests and to download the Matlab code). The profile and the pyramid are shown in the top left image and top right image in Figure 1, from the essay.

spectrum vs cubeYF

Figure 1

In the plot the value of L* varies with the color of each sample in the spectrum, and the line is colored accordingly. This erratic profile highlights several issues with spectrum: firstly, the change in lightness is not monotonic. For example it increases from black (L*=0) to magenta [M] then drops from magenta to blue [B], then increases again and so on. This is troublesome if spectrum is used to map elevation because it will interfere with the correct perception of relief, particularly if shading is added. Additionally, the curve gradient changes many times, indicating a nonuniform perceptual distance between samples. There are also plateaus of nearly flat L*, creating bands of constant color (a small one at the blue, and a large one at the green [G]).

The Great Pyramid has monotonically increasing elevation (in feet – easier to code) so there should be no discontinuities in the surface if the color palette is perceptual. However, clearly using the spectrum we have introduced many artificial discontinuities that are not present in the data. For the bottom row in FIgure 1 I used my new color palette, which has a nice, monotonic, compressive Lightness profile (bottom left). Using this palette the pyramid surface (bottom right) is smoothly colored, without any perceptual artifact.

This is how I created the palette: I started with RGB triplets for magenta, blue, cyan, green, and yellow (no red), which I converted to L*a*b* triplets using Colorspace transformations, a Matlab function available on the Matlab File Exchange. I modified the new L* values by fitting them to an approximately cube law L* function (this is consistent with Stevens’ power law of perception), and adjusted a* and b* values using Lab charts like the one in Figure 2 (from CIELab Color Space by Gernot Hoffmann, Department of Mechanical Engineering, University of Emden)  to get 5 colors moving up the L* axis along an imaginary spiral (I actually used tracing paper). Then I interpolated to 256 samples using the same ~cube law, and finally reconverted to RGB [1].

L*50_RGBval

Figure 2

There was quite a bit of trial and error involved, but I am very happy with the results. In the animations below I compare the spectrum and the new palette, which I call cubeYF, as seen in CIELab color space. I generated these animations with the method described in this post, using the 3D color inspector plugin in ImageJ:

I also added Matlab’s default Jet rainbow – a reminder that defaults may be a necessity, but in many instances not the ideal choice:

OK, the new palette looks promising, insofar as modelling is concerned. But how would it fare using some real data? To answer this question I used a residual gravity map from my unpublished thesis in Geology at the University of Rome. I introduced this map and discussed the geological context and objectives of the geophysical study in a previous post, so please refer to that if you are curious about it. In this post I will go straight to the comparison of the color palettes; if you are unfamiliar with gravity data, try to imagine negative residuals as elevation below sea level, and positive residuals as elevation above seal level – you won’t miss out on anything.

In Figures 3 to 6 I colored the data using the above three color palettes, and grayscale as benchmark. I generated these figures using Matlab code I shared in my post Visualization tips for geoscientists: Matlab, and I presented three of them (grayscale, Spectrum, and cubeYF) at the 2012 convention of the Canadian Society of Exploration Geophysicists in Calgary (the extended abstract, which I co-authored with Steve Lynch of 3rd Science, is available here).

In Figure 3, the benchmark for the following figures, I use grayscale to represent the data, assigning increasing intensity from most negative gravity residuals in black to most positive residuals in white (as labeled next to the colorbar). Then, I used terrain slope to create shading: the higher the slope, the darker the shading that is assigned, which results in a pseudo-3D display that is very effective (please refer to Visualization tips for geoscientists: Surfer, for an explanation of the method, and Visualization tips for geoscientists: Matlab for code).

Figure 3 - Grayscale benchmark

Figure 3 – Grayscale benchmark

In Figure 4 I color the pseudo-3D surface with the cubeYF rainbow. Using this color palette instead of grayscale allows viewers to appreciate smaller changes, more quickly assess differences, or conversely identify areas of similar anomaly, while at the same time preserving the peudo-3D effect. Now compare Figure 4 with Figure 5, where we use the spectrum to color the surface: this palette introduces several artefacts (sharp edges and bands of constant hue) which confuse the display and interfere with the perception of pseudo-relief, all but eliminating the effect.  For Figure 6 I used Matlab’s default Jet color palette, which is better that the spectrum, and yet the relief effect is somewhat lost (due mainly to a sharp yellow edge and cyan band).

campi cube YF

Figure 4 – cube YF rainbow

campi spectrum

Figure 5 – Industry spectrum

campi jet

Figure 6 – Matlab Jet

It looks like both spectrum and jet are poor choices when used for color representation of a surface, with the new color palette a far superior alternative. In the CSEG convention paper mentioned above (available here) Steve and I went further by showing that the spectrum not only has these perceptual artifacts and edges, but it is also very confusing for viewers with deficient color vision, a condition that occurs in about 8% of Caucasian males. We did that using computer software [2] to simulate how viewers with two types of deficient color vision, Deuteranopia and Tritanopia, would see the two colored surfaces, and we compare the results. In other words, we are now able to see the images as they would see them. Please refer to the paper for a full discussion on these simulation.

In here, I show in Figures  7 to 9 the Deuteranope simulations for cubeYF, spectrum, and jet, respectively. In all three simulations the hue discrimination has decreased, but while the spectrum and jet are now even more confusing, the cubeYF has preserved the relief effect.

Deuteranope Simulation of campi cube YF

Deuteranope Simulation of cube YF

Deuteranope Simulation of campi spectrum

Deuteranope Simulation of Industry spectrum

Deuteranope Simulation of campi jet

Deuteranope Simulation of Matlab Jet

That’s it for today. In my next post, to be published very shortly, you will get the palette, and a lot more.

References

A more perceptual color palette for structure maps, CSEG/CSPG 2012 convention, Calgary

How to assess a colourmap, in 52 Things You Should Know About Geophysics

Notes

[1] An alternative to the method I used would be to start directly in CIELab color space, and use a some kind of spiral *L lightness profile programmatically.  For example:

– Using 3D helical curves from: http://www.mathworks.com/matlabcentral/fileexchange/25177-3d-curves

– Using Archimedes spiral

– Expanding on code by Steve Eddins at Mathworks (A path through L*a*b* color space) in this article , one could create a spiral cube lightness with something like:

%% this creates best-fit pure power law function 
%  Inspired by wikipedia - http://en.wikipedia.org/wiki/Lightness
l2=linspace(1,power(100,0.42),256); 
L2=(power(l2,1/0.42))'; 

%% this makes cielab real cube function spiral 
radius = 50; 
theta = linspace(0.6*pi, 2*pi, 256).'; 
a = radius * sin(theta); b = radius * cos(theta); 
Lab1 = [L2, a, b]; RGB_realcube=colorspace('RGB<-Lab',(Lab1));

[2] The simulations are created using ImageJ, an open source image manipulation program, and the Vischeck plug-in. I later discovered Dichromacy, anther ImageJ plug-in for these simulations, which has the advantage of being an open source plugin. They can also be performed on the fly (no upload needed) using the online tool Color Oracle.

Related posts