ML IN SEISMIC INTERPRETATION

Facies analysis

Spectral decomposition

Stepan Goryachev
Data Analysis Center
9 min readJul 19, 2021

--

Preparing your data to train a neural network is no less important than the training process itself. The better features you have — the easier it is for a model to make decent predictions. That is why exploring data and generating features are important steps in composing your ML experiment.

In seismic analysis, features are often called attributes.

In the first part of the Facies analysis article, we have demonstrated the basic attributes that highlight anomalous horizon areas, each in its own way. Now when you are familiar with seismic data representation concepts, let’s dig in a bit deeper and talk about spectral decomposition. It’s a procedure commonly used by seismic experts during facies interpretation.

In this article, we will briefly explain the core idea of spectral decomposition and introduce the procedure of applying Fourier and Wavelet transform to seismic horizons.

Signal spectrum — what is it?

Seismic traces might look messy and chaotic at a first glance. But in fact, like many other signals, they are composed of multiple simple components and a bit of noise (sometimes not a bit, but a lot). Therefore one can approach seismic signals analysis by separating them into those simple components. The most common approaches to do this are Fourier Transform and Wavelet Transform.

Both procedures produce results that describe original signal properties as a function of its components. For Fourier transform, these components are sine and cosine waves corresponding to various frequencies; for Wavelet transform — wavelets themselves corresponding to the waveforms of various scales.

The spectrum obtained by decomposition into mentioned components is a sequence of numbers that describes the contribution of every component into the original signal.

This might sound a bit abstract, so I’ll throw in a few graphs. Here goes the seismic trace.

To capture the nature of the signal in both the big picture and details, one needs to apply decomposition that considers the impact of all its components. From a spectral analysis point of view (Fourier/Wavelet), these components are harmonics/wavelets of various frequencies/scales. And that is exactly how spectral attributes differ from Hilbert transform-related ones: while the latter accounts for only local information, spectral decomposition attributes account for both.

Decomposition of seismic trace into primitives that correspond to different frequencies/scales makes sense due to the nature of the signal itself. Properties of seismic reflection in a specific underground point depend on rock properties at that point. More specifically, the velocity of a seismic wave is quasi-constant in a homogenous stratum (yeah, “assume the cow is a sphere…”). Therefore rock layers of various properties produce reflections of various intensities, that are imprinted into resulting seismic trace.

Fourier power spectrum

So, how do you calculate the impact of various harmonics in your signal? That piece of code below might show you the way.

Results in this example will differ from those on the plot below since the signal is generated randomly.

If you are unfamiliar with numpy methods in use, better go read the docs!

And that is how the power spectrum of a seismic trace obtained via such transform looks like.

The spectrum shows how components of different frequencies contribute to the trace on the whole. It can be seen that the component’s contribution outside of the 10–60 Hz interval is negligible. As for components inside this interval, they are more perceptible.

And here is another view on the same spectrum data, where intensities are encoded with color. The brighter — the bigger the frequency contribution. As one might notice, the Fourier transform does not show us where exactly in signal the contribution of a specific frequency is greater and where it is smaller. The spectrum just describes the given signal in general, i.e., what harmonics is it made of.

Sure, there is a way to achieve time resolution with harmonics decomposition — short-time Fourier transform (SFTF). Apply the same procedure as above to consequent signal parts of the same length, and you’ll get data more detailed in time. But due to the very nature of Fourier transform you’ll lose in resolution in the frequency domain.

I encourage you to try SFTF out yourself (either modify an example provided above or use scipy.signal.stft) to observe described time-frequency resolution tradeoff with your own eyes!

Wavelet decomposition

There is a way to keep both time resolution and decompose signal in detail. Actually, that is the very aim of the Wavelet transform.

To apply the transform, one should choose the wavelet function and its scales. In seismic interpretation, a Ricker wavelet (also known as Mexican hat) function is used because of its similarity to the seismic waveform. As for the choice of wavelets’ widths, there is no simple rule of thumb for that. It depends greatly on signal properties, so the best advice for beginners would be “try different ranges”.

This piece of code uses scipy.signal.cwt method to calculate wavelet transform.

Results in this example will differ from those on the plot below since the signal is generated randomly.

Let’s look at the visualization of the transformed result. Contrary to Fourier transform, there is a time resolution here. The color of every cell in a row encodes the contribution of a specific wavelet at a specific time moment.

That was the most general overview of results produced by Fourier and Wavelet transforms correspondingly, if applied to a single signal. If you are completely new to these two procedures then for sure you must have more questions rather than answers. That's totally fine. If you wish to learn more about spectral analysis, I suggest the following nice introductions to Fourier transform (intuitive-friendly with great visualizations!) and wavelets (also great, just with fewer animations).

And if you are a glutton for punishment, try these articles [1], [2] to know more about the math basics behind the scenes.

Fourier transform attribute

To obtain an attribute map for a whole reflective surface, one needs to transform a fragment of a fixed length of every trace that goes through that surface (just like in the first part of the article). The image below shows the basic “amplitudes” attribute — a slice of seismic cube values cut along the horizon.

To get a Fourier transform result with decent frequency resolution one needs to obtain a slice at least a couple dozen counts wide along the time axis. Applied to that data transform gives a 3d array. Along the first two dimensions, it has a shape of an input data slice along inlines and crosslines. As for the third dimension, it has the number of channels equal to the number of components extracted by transform, which is half the signal’s length.

Fourier transform result always covers the frequency range from 0 to half of the sampling frequency. However, the shorter signal you use as a transform input, the less frequency-detailed picture you get.

The first way to visualize obtained attribute is to display transform components separately. They depict various structures on the horizon depending on the frequency they correspond.

The second approach to visualizing a multiple-channel array is to reduce the number of channels somehow. Seismic experts traditionally approach this problem by taking three specific channels (e.g., 15, 25, 35 Hz) and treating them as RGB channels.

Such an approach requires manual channels choice since various objects of interest might be present in different parts of the spectrum. Instead of following this procedure of supervised choice, we apply dimensionality reduction (via PCA) to the result of the Fourier transform. Reducing it to a single component gives the following image.

The interpretation of color can be really messy after PCA is applied to data. What is important on the attributes below, is how the color changes from one point to another.

As one might see, there are two noticeable differences between the “amplitudes” attribute and “Fourier” ones. The first is that “Fourier” is more contrast. The second is it blurs some details and sharpens others.

The mask labeled by the expert does not match the distinctive bright on the Fourier attribute part perfectly. That can mean either that this attribute depicts facies edges worse than raw amplitudes on that horizon (i.e. the feature is not very useful) or that the dimensionality reduction is not an ideal way to visualize such data (probably) or that the expert just did not take the spectral attribute into account (i.e. the provided labeling is subjective). The last case happens more often than one might suspect 🤷.

Wavelet transform attribute

As for wavelet transform, being applied to a slice of amplitudes along the horizon, it produces a 4d result — the first 3 dimensions correspond to the shape of the amplitude slice itself. The 4th dimension corresponds to the wavelet scale parameter.

Method for the wavelet transform from scipy is slow, when applied trace-wise. Therefore I suggest a better approach, that implements the same transform via 3d convolution applied to a whole amplitudes array. It works fast but can be improved even more with numba.cuda.

It’s assumed that the time dimension is along the last array axis.

Visualizing 4d data is hard, so firstly, one needs to reduce the number of dimensions. Slicing the wavelet transform result along the depth axis gives a 3d array that corresponds to a specific time. After that available approaches are the same as with Fourier-decomposed data.

On separate plots of wavelet-decomposed seismic data, various regions are highlighted depending on the wavelet scale parameter chosen for decomposition.

As for attribute, obtained by reducing several wavelet components into a single via PCA, it’s somehow similar to the “Fourier” attribute — image is also more contrast compared to “amplitudes”, also blurs some horizon regions and sharpens the other. However, it looks less extreme, i.e. contrast and details are changed not so intensely as on the “Fourier” attribute.

Facies gallery

In the illustration below introduced spectral attributes are compared to base cube values for various horizons from different fields. As one might notice, depending on a specific facies object different attributes highlight alluvial fans contours better or worse. That might suggest that there is no universal attribute that allows to easily segment facies on any seismic cube. Instead, various attributes might be useful depending on the situation.

Conclusion

Spectral attributes for sure highlight some areas of interest on reflective surfaces. However, it will be too bold if from here we just say that they are definitely will be helpful during training of the model that segments facies.

For sure, their ability to represent information about various objects along the horizon in different channels across frequency/scale dimension is interesting. Combined with attention mechanism in neural networks, usage of such features might allow making more detailed and accurate facies predictions depending on the frequency range the anomaly is present in.

However, such attributes might be also excessive for the model. The only valid way to know it is to make an experiment that checks various scenarios — training of multiple models with different features combinations and comparison of those models' predictions in terms of quality metrics.

To be continued

In the following articles of the series, we will share our approaches to the facies segmentation task itself. Firstly, we build a baseline that loads the data correctly and makes decent-quality predictions. After that goes the research of everything — features, data load defaults, model parameters, etc. To get familiar with a tool, that allows making such experiments read this article. Based on the results of that research we build an even better model.

Stay tuned!👋

--

--