Tag Archives: visualization

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.

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).

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.

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

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”.

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.

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.

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 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

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

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.

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

Do you know any cool apps?

Posted on by
8

I’d like to pick up my Apps page, which I sort of abandoned a while back.

If you have any great app to recommend, I’d love to hear about it so please add them in the comment section to this post. I am looking for Apps for Android and iPhone/iPad in the following categories – ideally free or very low-cost, possibly open-source:

Geology

Geophysics

Cartography and mapping

Planetary Science

Image Processing

Visualization

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

Google Earth and a 5 minutes book review: Geology Illustrated

Posted on by
2

A few years ago I bought on e-bay Geology Illustrated – by John S. Shelton, for just 10 US dollars. Every time I look at, and inside the book I can’t but think those were the best 10 dollars I ever invested in books.

There are already reviews and plenty of praise for this book out there – no need to repeat any of that if not briefly. My take is that the geology is clear and well explained. A bit simple, but simple is not always bad. And Shelton himself in the preface recommends this book as a “point of departure rather than something to lean on…” but that is perfect if you are a teacher looking for material, a first year college student, or a non-geologist looking for a high quality introduction.

But the photographs are priceless, and Shelton, who was also a pilot, took them all himself. Again, the author reminds us that nothing can replace field experience, and having  been trained as a field geologist (an average one, but that’s another story) I cannot but agree. However, lacking access or time to go to the field, or both, I find looking at a book like this can be an extraordinary substitute. That is especially true if you combine the reading with using Google Earth (particularly if you are a visual-spatial learner) and that is exactly what I did.

I already praised Google Earth for visualisation in this post. This program is a fantastic tool for learning geology, and today, to reinforce the point, I want to show you a couple of examples of Google Earth views replicating almost exactly figures from Chapter 14 of Geology Illustrated: The works of streams and rivers.

The first view is a replica of Figure 130 in the book, showing a fantastic example of a stream (the Colorado River) deepening its valley at the Marble Canyon.

The second view is a replica of Figure 135, showing many excellent examples of stream capture by headward erosion. Notice that in the 60s, when the photo was taken by Shelton, the highway (US Highway 101 north of San Juan Capistrano, California) was the only visible evidence of human activity.

The last view is a replica of Figure 137 in the book, showing the meander belt of the Animas River a few miles from Durango, Colorado. Looking at this was by far my favourite as it gave me the opportunity to create my own time lapse: a repeat snapshots of the same landscape nearly 50 years apart. Tis is priceless: 50 years are nothing in geological time scale, and yet there are already some significant differences in the two images. For example, it looks like the meander cutoff  in the lower left portion of the image had ‘just’ happened in the 60s, whereas at the time the imagery used by Google Earth was acquired (I imagine in the last few years), the remnant oxbow lake is more clearly defined. Another oxbow lake in the center has nearly disappeared.

I found that this process of looking for and replicating the photos in the book, zooming in and out, then in again changing view was a fantastic way to see the geological features as part of the larger geological context, visualize them, see the interconnection with other elements of the landscape, observe how erosion and deposition, and human processes have modeled the landscape in just a few decades (as in the second and third examples).  As a geophysicist, sitting in the office away from the outcrops, this is  invaluable, and a great aid in finding analogs in support of seismic interpretations. And really you don’t need a book in your lap to start the process.

In a future post I will show my results at creating similar views using HD lidar data, which can be downloaded from the National Center for Airborne Laser Mapping, as done in this blog post on Quest.

Resources

John Shelton’s obituary, August 2008

Geomorphology from space

Visualization tips for geoscientists: Matlab, part III

Posted on by
5

 Introduction

Last weekend I had a few hours to play with but needed a short break from writing about color palettes, so I decided to go back and finish up (for now) this series on geoscience visualization in Matlab. In the first post of the series I expanded on work by Steve Eddins at Mathworks on overlaying images using influence maps and demonstrated how it could be used to enhance the display of a single geophysical dataset.

Using transparency to display multiple data sets an example

At the end of the second post I promised I would go back and show an example of using transparency and influence maps for other tasks, like overlaying of different attributes. Here’s my favorite example in Figure 1. The image is a map in pastel colors of the Bouguer Gravity anomaly for the Southern Tuscany region of Italy, with three other layers superimposed using the techniques above mentioned.

It is beyond the objectives of this post to discuss at length about gravity exploration methods or to attempt a full interpretation of the map. I will go back to it at a later time as I am planning a full series on gravity exploration using this data set, but if you are burning to read more about gravity interpretation please check these excellent notes by Martin Unsworth, Professor of Physics at the Earth and Atmospheric Sciences department, University of Alberta, and note 4 at the end of this post. Otherwise, and for now, suffice it to say that warm colors (green to yellow to red) in the Bouguer gravity map indicate, relatively speaking, excess mass in the subsurface and blue and purple indicate deficit of mass in the subsurface.

The black and grey lines are lineaments extracted from derivatives of the Bouguer gravity data using two different methods [1]. The semitransparent, white-filled polygons show the location of some of the  basement outcrops (the densest rocks in this area).

Lineaments extracted from gravity data can correspond to contacts between geological bodies of different density, so a correlation can be expected between basement outcrops and some of the lineaments, as they are often placed in lateral contact with much lesser dense rocks. This is often exploited in mineral exploration in areas such as this where mineralization occurs at or in the vicinity of this contacts. As an example, I show in Figure 2 the occurrences (AGIP – RIMIN, unpublished industry report, 1989) of silicization (circles) and antimony deposits (triangles), superimposed on the distance from one of the set of lineaments (warm colors indicate higher distance) from Figure 1.

The fact that different methods give systematically shifted results is a known fact, often due the trade-off between resolution and stability, whereby the more stable methods are less affected by noise, but often produce smoother edges over deeper contacts, and their maxima may not correspond. This is in addition to the inherent ambiguity of gravity data, which cannot, by themselves, be interpreted uniquely. To establish which method might be more correct in this case (none is a silver bullet) I tried to calibrate the results using basement outcrops (i.e. does either method more closely match the outcrop edges?). Having done that, I would have more confidence in making inferences on possible other contacts in the subsurface suggested by lineament. I would say the black lines do a better overall job in the East, the gray perhaps in the West. So perhaps I’m stuck? I will get back to this during my gravity series.

Figure 1

Figure 2

Matlab code

As usual I am happy to share the code I used to make the combined map of Figure 1. Since the data I use is in part from my unpublished thesis in Geology and in part from Michele di Filippo at the University of Rome, I am not able to share it, and you will have to use your own data, but the Matlab code is simply adapted. The code snippet below assume you have a geophysical surface already imported in the workspace and stored in a variable called “dataI”, as well as the outcrops in a variable called “basement”, and the lineaments in “lnmnt1” and “lnmnt2”. It also uses my cube1 color palette.

 
% part 1 - map gravity data
figure; imagesc(XI,YI,dataI); colormap(cube1); hold on;
%
% part 2 - dealing with basement overlay
white=cat(3, ones(size(basement)), ones(size(basement)),...
 ones(size(basement)));
ttt=imagesc(Xb,Yb,white); % plots white layer for basement
%
% part 3 - dealing with lineaments overlays
black=cat(3, zeros(size(lnmnt1)), zeros(size(lnmnt1)),...
 zeros(size(lnmnt1)));
grey=black+0.4;
basement_msk=basement.*0.6;
kkk=imagesc(XI,YI,black); % plots black layer for lineament 1
sss=imagesc(XI,YI,gray); % plots gray layer for lineament 2
hold off
%
% part 4 - set influence maps
set(ttt, 'AlphaData', basement_msk); % influence map for basement
set(kkk, 'AlphaData', lnmnt1); % influence map for linement 1
set(sss, 'AlphaData', lnmnt2); % influence map for linement 2
%
% making it pretty
axis equal
axis tight
axis off
set(gca,'YDir','normal');
set(gcf,'Position',[180 150 950 708]);
set(gcf,'OuterPosition',[176 146 958 790]);

Matlab code, explained

OK, let’s break it down starting from scratch. I want first to create a figure and display the gravity data, then hold it so I can overlay the other layers on top of it. I do this with these two commands:

figure;imagesc(XI,YI,dataI);

hold on;

The layer I want to overlay first is the one showing the basement outcrops. I make a white basement layer covering the full extent of the map, which is shown in Figure 3, below.

Figure 3

I create it and plot it with the commands:

white=cat(3, ones(size(basement)), ones(size(basement)), ones(size(basement)));

ttt=imagesc(Xb,Yb,white);

The handle  ttt is to be used in combination with the basement influence map to produce the partly transparent basement overlay: remember that I wanted to display the outcrops in white color, but only partially opaque so the colored gravity map can still be (slightly) seen underneath. I make the influence map, shown in Figure 4, with the command:

basement_msk=basement.*0.6;

Since the original binary variable “basement” had values of 1 for the outcrops and 0 elsewhere, whit the command above I assign an opacity of 0.6 to the outcrops, which will be applied when the next command, below, is run, achieving the desired result.

set(ttt, ‘AlphaData’, basement_msk); % uses basement influence map

Figure 4

For the lineaments I do things in a similar way, except that I want those plotted with full opacity since they are only 1 pixel wide.

As an example I am showing in Figure 5 the black layer lineament 1 and in Figure 6 the influence map, which has values of 1 (full opacity) for the lineament and 0 (full transparency) for everywhere else.

Figure 5

Figure 6

Now a few extra lines to make things pretty, and this is what I get, shown below in Figure 7: not what I expected!

Figure 7

The problem is in these two commands:

white=cat(3, ones(size(basement)), ones(size(basement)), ones(size(basement)));

ttt=imagesc(Xb,Yb,white);

I am calling the layer white but really all I am telling Matlab is to create a layer with maximum intensity (1). But the preceding colormap(cube1) command assigned a salmon-red color to the maximum intensity in the figure, and so that is what you get for the basement overlay.

Again, to get the result I wanted, I had to come up with a trick like in the second post examples. This is the trick:

I create a new color palette with this command:

cube1edit=cube1; cube1edit(256,:)=1;  

The new color palette has last RGB triplet actually defined as white, not salmon-red.

Then I replace this line:

figure; imagesc(XI,YI,dataI); colormap(cube1); hold on;

with the new line:

figure; imagesc(XI,YI,dataI, [15 45]); colormap (cube1edit); hold on;

The highest value in dataI is around 43. By spreading the color range from [15 43] to [15 45], therefore exceeding max(dataI) I ensure that white is used for the basement overlay but not in any part of the map where gravity is highest but there is no basement outcrop. In other words, white is assigned in the palette but reserved to the overlay.

Please let me know if that was clear. If it isn’t I will try to describe it better.

Notes

[1] One method is the total horizontal derivative. The other method is the hyperbolic tilt angle – using Matlab code by Cooper and Cowan (reference). This is how I produced the two overlays:  first I calculated the total horizontal derivative and the tilt angle, then I found the maxima to use as the overlay layers. This is similar to Figure 3e in Cooper and Cowan, but I refined my maxima result by reducing them to 1-pixel-wide lines (using a thinning algorithm).

Reference

Cooper, G.R.J., and Cowan, D.R. (2006) – Enhancing potential field data using filters based on the local phase  Computers & Geosciences 32 (2006) 1585–1591

Related posts (MyCarta)

Visualization tips for geoscientists: Surfer

Visualization tips for geoscientists: Matlab

Visualization tips for geoscientists: Matlab, part II

Image Processing Tips for Geoscientists – part 1

What is a colour space?

Posted on by
7

Since I am deconstructing the rainbow/spectrum in RGB, HSV color spaces, I will show it in dichromat “color space”, and finally will then make a new one in CIE L*a*b* color space I thought I should include a review of the basics. So what is a color space? This good post answers that question.

Color Use Guidelines for Mapping and Visualization

Posted on by
7

I find the Color Use Guidelines by Cynthia Brewer (of Color Brewer fame) very well done and extremely useful. Below is a screen captures of the main page. To see an explanation of and example for each color scheme visit the interactive guidelines.

The rainbow is dead…long live the rainbow! – Perceptual palettes, part 1

Posted on by
33

Introduction

This is the first  post in a series on the rainbow and similar color palettes. My goal is to demonstrate it is not a good idea to use these palettes to display scientific data, and then answer these two questions: (1) is there anything we can do to “fix” the rainbow, and (2) if not, can we design a new one from scratch.

The rainbow is dead…some examples

In a previous post I showed a pseudo-3D rendering of my left hand x-ray using intensity (which is a measure of bone thickness) as the elevation. I mapped the rendering to both grayscale and rainbow color palettes, and here I reproduced the two images side by side:


I used this example to argue (briefly) that the rainbow obscures some details and confuses images by introducing artifacts. Notice that in this case it clearly reduces the effectiveness of the pseudo-3D rendering in general. It also introduces inversions in the perception of elevation. The thick part in the head of the radius bone, indicated by the arrow, looks like a depression, whereas it is clearly (and correctly) a high in the grayscale version.

Continue reading →