Visualizing normal moveout (NMO) in Python

Recently I was reading in Yilmaz - Seismic Data Analysis (1987) about normal moveout (NMO) correction, a conventional method for velocity analysis based on the assumption of hyperbolic moveout. A key figure explaining the method is reproduced in Figure 1.

Figure 1  NMO correction involves mapping nonzero-offset traveltime t onto zero-offset traveltime  t₀. (a) Before and (b) after NMO correction. (From Yilmaz (1987), Fig. 3.1-3).

I thought it would be instructive to reproduce this figure using Python. Following Yilmaz' example, there is a single shallow event at t0 = 0.8 s and a medium velocity of v = 2264 m/s. The offset range is 0 to 3150 m, with a 50-m trace separation. The wavelet is assumed to be a 20-Hz Ricker wavelet. Here is the code:

import numpy as np
import matplotlib.pyplot as plt
from math import pi
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show

def Ricker(f0,t): # Ricker wavelet
return (1-2*(pi*f0*t)**2)*np.exp(-(pi*f0*t)**2)

f0=20. # Dominant frequency of Ricker wavelet (in Hz)
t0=0.8 # Arrival time at zero offset (in s)
v=2264. # Velocity of the subsurface (in m/s)

Xmax=3150. # Maximum offset (in m)
dX=5. # Trace spacing (in m)
Tmax=2. # Maximum time (in s)
dT=0.001 # Sampling period (in s)

t=np.arange(0,Tmax,dT) # Coordinate vector for time
x=np.arange(0,Xmax,dX) # Coordinate vector for offset

(X,T)=np.meshgrid(x,t) # Coordinate matrices for time and offset

def tr(x): # Traveltime from source to receiver (in s)
return np.sqrt(t0**2+(x/v)**2) # Normal moveout with constant velocity

def w(x,t): # Wavefield
return Ricker(f0,t-tr(x)) # Ricker wavelet from a single reflector at t=tr

def w_NMO(x,t0):
tx=np.sqrt(t0**2+x**2/v**2)
return w(x,tx)

W=w(X,T) # Wavefield
w_NMO=w_NMO(X,T) # NMO-corrected wavefield

plt.close("all")

# Travel time as function of offset
plt.figure()
plt.plot(x,tr(x))
plt.ylim([0.,Tmax])
plt.xlim([0.,Xmax])
plt.xlabel('Offset (m)')
plt.ylabel('Two-way time (s)')
plt.gca().invert_yaxis()

# Wavefield
plt.figure()
im = imshow(W,cmap=cm.RdBu,aspect='auto',extent=[0,Xmax,Tmax,0])
plt.xlabel('Offset (m)')
plt.ylabel('Two-way time (s)')

# NMO-corrected wavefield
plt.figure()
im = imshow(w_NMO,cmap=cm.RdBu,aspect='auto',extent=[0,Xmax,Tmax,0])
plt.xlabel('Offset (m)')
plt.ylabel('Two-way time (s)')

plt.show()

The figures generated by this script are shown below. Figure 2 shows the traveltime as a function of offset, and Figure 3 shows the seismic data before and after NMO correction.

Figure 2  Traveltime as a function of offset for Yilmaz' example, showing hyperbolic moveout.

Figure 3  Wavefield before (above) and after (below) NMO correction. Blue and red colors represent positive and negative amplitudes, respectively.

As seen from Figure 3, in the NMO-corrected common midpoint (CMP) gather, the reflector appears to 'stretch' with increasing offset. This phenomenon is (aptly) called "NMO stretch".

In short, this example shows how to implement NMO correction for analytically defined seismic data and display the results in a color plot using Python.

Acknowledgement

I followed an example by The Glowing Python.