Phase and the Hilbert transform
|This tutorial originally appeared as a featured article in The Leading Edge in October 2014 — see issue.|
The concept of phase permeates seismic data processing and signal processing in general, but it can be awkward to understand, and manipulating it directly can lead to surprising results. It doesn't help that the word phase is used to mean a variety of things, depending on whether we refer to the propagating wavelet, the observed wavelet, poststack seismic attributes, or an entire seismic data set. Several publications have discussed the concepts and ambiguities   .
In this tutorial, we will focus on aspects of phase relevant to the interpreter. We will look at how to manipulate the phase of a seismic trace by manipulating the phase of its Fourier transform and will use that idea to generate the well-known instantaneous-phase poststack attribute. We will also check how close to zero phase our test data set is.
The plots included in this tutorial were created using GNU Octave. The code to produce them, some of which is included here, is available in the SEG tutorials GitHub repository. The data we are using are the Penobscot 3D survey from the Open Seismic Repository, which is openly licensed CC-BY-SA by the Canada Nova Scotia Offshore Petroleum Board, and dGB Earth Sciences.
The Hilbert transform
The Fourier transform is complex. Taking the transform of any real signal will result in a set of complex coefficients. Complex numbers are essentially 2D vectors, meaning they have two components: magnitude and phase angle. Most of the time when dealing with Fourier transforms, we concentrate on magnitude, which tells us about the distribution of signal energy through frequency. However, every signal also has a phase spectrum, and the phase encodes the signal's structure - the distribution of the signal energy through time. We do not often examine phase spectrum because it is difficult to interpret, but we can manipulate Fourier phase to change the structure of our signal without affecting its amplitude spectrum.
The Hilbert transform is a linear operator that produces a 90° phase shift in a signal, and it is a good first step in our exploration of phase. It is also commonly used in poststack seismic analysis to generate the analytic signal from which we can compute the standard complex trace attributes such as envelope, instantaneous phase, and instantaneous frequency.
The definition of the Hilbert transform is rather cryptic; it is much easier to consider in terms of its Fourier transform definition. The Hilbert transform H of a signal u is related to the Fourier transform F like this:
So we apply a Hilbert transform by multiplying all negative frequencies by i and all positive frequencies by –i, leaving any DC component untouched. This is perhaps not intuitive, but we can gain some additional insight by thinking about it geometrically.
The multiplication of two complex numbers a (such as a Fourier coefficient) and b (such as above) can be thought of as a rotation in the complex plane. Euler's formula helps us to see that to rotate any complex number by a specific angle , we must multiply it by the complex number . We see that multiplication by i alone, as in the equation above, is equivalent to a rotation by 90°. When we take the inverse Fourier transform, the result is a phase-shifted version of our signal.
Indeed, we can generalize the definition of the Hilbert above to produce a phase shift to any angle, α:
Phase shifting in GNU Octave
Let us use this to phase-shift a seismic trace. We will take the fast Fourier transform (FFT) of our signal, manipulate the coefficients, and apply the inverse FFT. The following code implements this for an array x (see fftshifter.m in the repo):
% Create an array to apply to the coefficients % positive and negative frequencies appropriately R0 = exp(i∗phase_shift_in_radians); N = length(x); R = ones([1 N]); M = ceil((N+1)/2); R(1:M) = R0; R(M+1;N) = conj(R0); % Apply the phase shift in the frequency domain Xshifted = R.∗fft(x); % Recover the shifted time domain signal y = real(ifft(Xshifted));
Complex trace attributes
Now that we have this ability to phase-shift, what can we do with it? We can compute some seismic attributes. A whole set of commonly used attributes depends on the Hilbert transform:
% Load a seismic trace; this trace was % exported from OpendTect as a Simple File filename = ‘data/trace_il1190_xl1155.trace’ [s, t] = load_simple_trace(filename); % Compute the analytic (complex) trace z = s + i∗fftshifter(s,-pi/2); envelope = abs(z); phase = angle(z);
The code above generates the Hilbert transform, envelope, and phase traces for a trace from the Penobscot data set. A section of this trace is shown in Figure 1 (see plot_complex_attributes_on_a_trace.m for full plotting code).
Once you have the complex trace representation, z, there are many additional attributes that you can compute, such as instantaneous frequency, phase acceleration, and envelope-weighted frequency, merely by playing around with the complex trace.
Checking the phase of your seismic
The complex attributes generated from a seismic trace and its Hilbert transform are useful in seismic interpretation. However, the quality of our seismic interpretation depends on another aspect of phase, that is, the phase of our seismic data set and whether the observed wavelet has been corrected to a zero-phase response during processing.
Roden and Sepúlveda (1999) give a practical way to assess the phase of seismic data as a quality check for seismic interpretation. They point out that where we have a strong reflection from a sharp geologic interface that we can confidently expect to produce a zero-phase response, we can measure the actual response in the seismic data set on that horizon to determine the amount of phase error present.
The principle that Roden and Sepúlveda (1999) apply is that at strong reflections in zero-phase data, we would expect peaks (or troughs) in the trace to align with peaks in the envelope. Therefore, any horizon that we pick on zero-phase data could also be picked on the envelope of that data.
Figure 2 shows a segment of our trace from Figure 1 with envelope peaks marked as spikes (see find_peaks.m). In Figure 2b, I have isolated those events and plotted the actual signal phase against the phase we would expect for a zero-phase data set. The difference between the two could be systematic residual phase. If we think a particular reflector has minimal interference and no AVO or other effects, then the rest is residual phase. The strong trough at the Abenaki reflector in Figure 2 suggests that there could be 0.64 radians, or about 37°, of residual phase in the data. That is something that we could now correct for by phase-shifting all our traces.
The Hilbert transform opens up a world of seismic attributes, some of which have everyday application for the interpreter. To see how we can extend them to 3D data and extract volumes of residual phase, check the expanded notes and full code at the SEG tutorials GitHub page.
- ↑ 1.0 1.1 1.2 1.3 Roden, R., and H. Sepúlveda, 1999, The significance of phase to the interpreter: Practical guidelines for phase analysis: The Leading Edge, 18, no. 7, 774–777, http://dx.doi.org/10.1190/1.1438375.
- ↑ Liner, C. L., 2002, Phase, phase, phase: The Leading Edge, 21, no. 5, 456–457, http://dx.doi.org/10.1190/1.1885500.
- ↑ Simm, R., and R. E. White, 2002, Tutorial: Phase, polarity and the interpreter’s wavelet: First Break, 20, no. 5, 277–281.
- Madagascar workflow - reproducible using Madagascar open-source software
- Corresponding author: Steve Purves,