 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

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
seismic = line

trace = seismic[100, :]

vmin, vmax = np.min(line), np.max(line)
vrng = max(abs(vmin), vmax)

print(f"inline range {f.ilines}, {f.ilines[-1]}")
print(f"xline range {f.xlines}, {f.xlines[-1]}")
print(f"sample range {f.samples}, {f.samples[-1]}")
print(f"data range {vmin}, {vmax}")
import matplotlib.pyplot as plt
%matplotlib inline

plt.figure(figsize=(18,11))
plt.imshow(seismic.T, vmin=-vrng/2, vmax=vrng/2, cmap='seismic', aspect='auto', interpolation='bilinear')
_=plt.title(f"Penobscot - partial inline 1020")
plt.figure(figsize=(16,3))
plt.plot(trace)
plt.grid('on')
plt.title('Trace from F3')
plt.ylim(-7000,7000)

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

# Create an array to apply the coefficient to
# positive and negative frequencies appropriately
N = len(x);
R = np.ones(x.shape, dtype=np.complex);
R=0.;
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.figure(figsize=(16,3))
plt.plot(trace)
plt.plot(yp,':')
plt.grid('on')
plt.title(f"Trace and {rotation_angle_degrees} degree phase shifted version")
plt.ylim(-7000,7000)

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) + 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):
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))

il_idx=100
t = np.arange(400,1104, 4.0)
axs.plot(t,seismic[il_idx, :])
axs.plot(t, hilbert[il_idx, :], ':')
axs.plot(t, env[il_idx, :], 'k', linewidth=1)
axs.plot(t, -env[il_idx, :], 'k', linewidth=1)
axs.set_title('Trace Amplitudes')
axs.grid('on')
axs.set_xlim(400,1100)
axs.legend(labels=['seismic','hilbert','envelope'])

axs.plot(t,phase[il_idx, :], 'k', linewidth=1)
axs.set_title('Instantaneous Phase')
axs.grid('on')
axs.set_xlim(400,1100)

### #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')
ax.set_title(title)
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
cbar = fig.colorbar(im,
ax=ax,
ticks=ticks,
orientation='horizontal',
cbar.ax.set_xticklabels(ticklabels)

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

plt.tight_layout()

## #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
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():
warnings.simplefilter("ignore")
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.plot(t,seismic[il_idx, :])
axs.plot(t,env[il_idx, :], linewidth=1)
axs.plot(t,env[il_idx, :],'k', linewidth=1)
axs.plot(t,-env[il_idx, :],'k', linewidth=1)
axs.stem(t,penv, 'C1--','C1o',':')
axs.stem(t,-penv, 'C1--','C1o',':')
axs.grid('on')
axs.set_xlim(400,1100)
axs.legend(labels=['seismic','peaks','envelope'])
axs.set_title('Envelope Peaks')

axs.stem(t, idealised_phase,'C1:','C1o',':')
axs.stem(t, pphase,'C0:','C0x',':')
axs.grid('on')
axs.set_xlim(400,1100)
axs.legend(labels=['idealised','actual'])
axs.set_title('Ideal & Actual Phase at Envelope Peaks')

plt.subplots_adjust(hspace=0.25)

## #Plot phase error on a section

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

for x in tqdm(range(line.shape)):
_, pphase = find_peaks(env[x, :], phase[x, :])
idealised_phase = get_idealised_phase(pphase)
phase_error[x,:] = np.abs(idealised_phase-pphase)
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')
References
1. Purves, S. (2014). Phase and the Hilbert transform. The Leading Edge, 33(10), 1164–1166. 10.1190/tle33101164.1