Juyter Logo


Steve Purves

#Phase and the Hilbert Transform

The first verision of this tutorial was published in 2014 in where GNU Octave was used for coding and plotting.

A python version is long overdue and thank's to the notebook we can now bring all the example code inline easily.

#setup (conda)

conda env create -f environment.yml
conda activate seg-tutorial-phase

#Load Data

We have 1 SEGY data file available:

  1. data/f3_seismic.sgy a subset of the Netherlands F3 Seismic Survey
import segyio
import numpy as np

with segyio.open('./data/f3_seismic.sgy') as f:
    line = f.xline[1020]
    seismic = line
    trace = seismic[100, :]
    vmin, vmax = np.min(line), np.max(line)
    vrng = max(abs(vmin), vmax)
    print(f"inline range {f.ilines[0]}, {f.ilines[-1]}")
    print(f"xline range {f.xlines[0]}, {f.xlines[-1]}")    
    print(f"sample range {f.samples[0]}, {f.samples[-1]}")
    print(f"data range {vmin}, {vmax}")
inline range 200, 650
xline range 700, 1200
sample range 400.0, 1100.0
data range -25080.0, 21983.0
import matplotlib.pyplot as plt
%matplotlib inline

plt.imshow(seismic.T, vmin=-vrng/2, vmax=vrng/2, cmap='seismic', aspect='auto', interpolation='bilinear')
_=plt.title(f"Penobscot - partial inline 1020")
<Figure size 1296x792 with 1 Axes>
plt.title('Trace from F3')
(-7000.0, 7000.0)
<Figure size 1152x216 with 1 Axes>

#Phase Shifting

Our phase shifting approach works by manipulating the FFT coefficents of the trace. The code below has been ported from the original Octave code in the TLE version and can be used to produce an arbitrary phase shift.

from math import ceil,exp
from scipy.fftpack import fft, ifft

def fftshifter(x, phase_shift_in_radians):
    # Create an array to apply the coefficient to 
    # positive and negative frequencies appropriately
    N = len(x);
    R0 = np.exp(-1j*phase_shift_in_radians);
    R = np.ones(x.shape, dtype=np.complex);
    R[1:N//2] = R0;
    R[N//2:] = np.conj(R0);

    # Apply the phase shift in the frequency domain
    Xshifted = R*fft(x);
    # Recover the shifted time domain signal
    y = np.real(ifft(Xshifted));
    return y
rotation_angle_degrees = 90
yp = fftshifter(trace, np.pi*rotation_angle_degrees/180)

plt.title(f"Trace and {rotation_angle_degrees} degree phase shifted version")
(-7000.0, 7000.0)
<Figure size 1152x216 with 1 Axes>

Which can then of course be used to perform a Hilbert Transform y shifting by 90 degrees.

But! of course we are now in python with scipy available and can benefit from a built in hilbert transform functions. Note there are 2!

that have different calling semantics but a hilbert transform is a hilbert transform and below we can see they are equivalent to the fftshifter at 90 degrees.

TODO: Currently there are differences, as we are not applying a window to our trace prior to FFT computation in the function above. Needs to be addressed, check the scipy implementations to see what they are using.

#Complex Trace Attributes

Computing the Hilbert Transform allows us to use the Analytic Trace representation

z=s(t)+iH(s(t))z = s(t) + i \cdot H(s(t))

to compute a variety of well known seimsic trace attributes. For example, envelope and instantaneous phase.

from scipy.signal import hilbert as hilbert_transform
hilbert = np.zeros_like(seismic)

for x in range(seismic.shape[0]):
    hilbert[x,:] = np.imag(hilbert_transform(seismic[x,:]))

z = seismic + 1j*hilbert
env = np.abs(z)
phase = np.angle(z)

#On a Trace

fig, axs = plt.subplots(2,1, figsize=(16,6))

t = np.arange(400,1104, 4.0)
axs[0].plot(t,seismic[il_idx, :])
axs[0].plot(t, hilbert[il_idx, :], ':')
axs[0].plot(t, env[il_idx, :], 'k', linewidth=1)
axs[0].plot(t, -env[il_idx, :], 'k', linewidth=1)
axs[0].set_title('Trace Amplitudes')

axs[1].plot(t,phase[il_idx, :], 'k', linewidth=1)
axs[1].set_title('Instantaneous Phase')
(400.0, 1100.0)
<Figure size 1152x432 with 2 Axes>

#On a Section

fig, axs = plt.subplots(2,2, figsize=(16,10))

def plot_it(ax, data, vmin, vmax, title, cmap, ticks, ticklabels, alpha=1.0, interpolation='bilinear'):
    im = ax.imshow(data.T, vmin=vmin, vmax=vmax, cmap=cmap, aspect='auto', alpha=alpha, interpolation=interpolation)
    ax.axvline(il_idx, 0, 250, color='black')
    cbar = fig.colorbar(im, 
                        fraction=0.03, pad=0.03)

plot_it(axs[0][0], seismic, -vrng/2, vrng/2, 'Seismic', 'seismic', [-vrng/2, 0, vrng/2], ['-ve', '0', '+ve'])
plot_it(axs[0][1], hilbert, -vrng/2, vrng/2, 'Hilbert', 'seismic', [-vrng/2, 0, vrng/2], ['-ve', '0', '+ve'])
plot_it(axs[1][0], env, 0, vrng/2, 'Envelope', 'viridis', [0, vrng/2], ['0', 'max'])
plot_it(axs[1][1], phase, -np.pi, np.pi, 'Instantaneous Phase', 'twilight', [-np.pi, 0, np.pi], ['tau/2', '0', 'tau/2'])

<Figure size 1152x720 with 8 Axes>

#Checking the phase of the seismic

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.

We can use this principle to check for residual phase at given reflectors

from math import ceil
from scipy.ndimage.filters import maximum_filter1d
import warnings

def find_peaks(x, a):
    y = np.full(x.shape, np.nan)
    ya = np.full(x.shape, np.nan)
    xmax = maximum_filter1d(x, 3)
    max_mask = x >= xmax
    y[max_mask] = x[max_mask]
    ya[max_mask] = a[max_mask]
    return y, ya

def get_idealised_phase(phase_at_peaks):
    # flip to nearest 'zero' phase point
    ideal_phase = np.full(phase_at_peaks.shape, np.nan);
    with warnings.catch_warnings():
        ideal_phase[np.abs(phase_at_peaks) < np.pi/2] = 0.;    
        ideal_phase[np.abs(phase_at_peaks - np.pi) < np.pi/2] = np.pi;
        ideal_phase[np.abs(phase_at_peaks + np.pi) < np.pi/2] = -np.pi;
    return ideal_phase
penv, pphase = find_peaks(env[il_idx, :], phase[il_idx, :])
idealised_phase = get_idealised_phase(pphase)
fig, axs = plt.subplots(2,1, figsize=(16,7))

t = np.arange(400,1104, 4.0)

axs[0].plot(t,seismic[il_idx, :])
axs[0].plot(t,env[il_idx, :], linewidth=1)
axs[0].plot(t,env[il_idx, :],'k', linewidth=1)
axs[0].plot(t,-env[il_idx, :],'k', linewidth=1)
axs[0].stem(t,penv, 'C1--','C1o',':')
axs[0].stem(t,-penv, 'C1--','C1o',':')
axs[0].set_title('Envelope Peaks')

axs[1].stem(t, idealised_phase,'C1:','C1o',':')
axs[1].stem(t, pphase,'C0:','C0x',':')
axs[1].set_title('Ideal & Actual Phase at Envelope Peaks')

<Figure size 1152x504 with 2 Axes>

#Plot phase error on a section

from tqdm.notebook import tqdm
phase_error = np.zeros_like(seismic)

for x in tqdm(range(line.shape[0])):
    _, pphase = find_peaks(env[x, :], phase[x, :])
    idealised_phase = get_idealised_phase(pphase)
    phase_error[x,:] = np.abs(idealised_phase-pphase)
HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=226.0), HTML(value='')))

fig, ax = plt.subplots(1,1, figsize=(16,8))

plot_it(ax, phase_error, 0, np.pi/2, 'Phase Error', 'plasma', [0, np.pi/2], ['0', '45deg'], interpolation='none')
<Figure size 1152x576 with 2 Axes>
  1. Purves, S. (2014). Phase and the Hilbert transform. The Leading Edge, 33(10), 1164–1166. 10.1190/tle33101164.1