Synthetic seismograms can be created by doing basic calculus on travel-time functions. Integrating slowness (the reciprocal of velocity) can a yield a time-to-depth relationship for making a well tie. Differentiating an acoustic impedance (velocity times density) log yields a reflectivity function along the borehole, which can be convolved with a seismic wavelet. In effect, the integral tells us where a given interval of rock is in the time-domain, the derivative tells us what it should look like; its seismic character.
My goal is to demonstrate the robustness of simple mathematical operations of travel time,
The code and synthetic data required to reproduce the results and figures can be found here in an accompanying IPython Notebook. For this tutorial, we will make use of the open-source Python libraries NumPy and Matplotlib, and there is a handy script for parsing LAS files here: http://wiki.scipy.org/Cookbook/LASReader.
As an example, I will use the sonic and density log from the Penobscot L-30 well. This well sits within the same volume of 3D seismic data that Matt Hall used in his horizon smoothing tutorial, it's a vertical well, so we don't have to worry about deviated trajectories through time and space.
Make sure you are in the correct directory, or cd to the correct one. Then let's check the files.
In [1]:
ls
Import NumPy, and the LASReader module from las.py. The las.py can be downloaded from the SciPy Cookbook (URL above), but we also give it here in the same directory. We create a well object called L30 by passing the 'L-30.las' file into the LASReader module.
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
In [3]:
import numpy as np
from las import LASReader
L30 = LASReader('L-30.las', null_subs=np.nan)
Let's take a look at which curves are in the file,
In [4]:
print L30.curves.names
Let's take a look at the measured depth units for this log.
In [5]:
print L30.curves.DEPTH.units
The units for this log are in feet (yuck). It turns out that we will have to convert to metres a number of times so let's write a function to do this,
In [6]:
def f2m(item_in_feet):
"converts feet to meters"
try:
return item_in_feet / 3.28084
except TypeError:
return float(item_in_feet) / 3.28084
return converted
For the curves that we will be accessing in the code, let's given them some shorter names. Also, some of the curves need to be converted into SI units,
In [7]:
z = f2m(L30.data['DEPTH']) # convert feet to metres
GR = L30.data['GRS']
IL8 = L30.data['LL8']
ILM = L30.data['ILM']
ILD = L30.data['ILD']
NPHISS = L30.data['NPHISS']
NPHILS = L30.data['NPHILS']
ILD = L30.data['ILD']
DT = L30.data['DT']*3.28084 # convert usec/ft to usec/m
RHOB = L30.data['RHOB']*1000 # convert to SI units
To deal with the shallow section, we first need to adjust the depths relative to sea level by subtracting the KB elevation (30.2 m) from the measured depths,
If we now integrate the sonic log, we see that time t = 0 corresponds to a depth of 347 m TVDss. To position the top of the log at the correct travel time on the seismic section, we need to use the thickness, and replacement velocities for both the water column, and the section above the log.
The replacement velocity for water = 1480 m/s
The replacement velocity for section above log = 1600 m/s
Update - missing variables now defined (Jan 5, 2014)
In [8]:
KB_elev = f2m(L30.well.KB.data) # Kelly Bushing elevation(ft)
water_depth = f2m(L30.well.GL.data) # Depth to Sea Floor below sea level (ft)
top_of_log = f2m(L30.start) # top of log (ft) relative to KB (actually 1150 ft)
print "KB elevation [m]: ", f2m(L30.well.KB.data) # Kelly Bushing (ft)
print "Seafloor elevation [m]: ", f2m(L30.well.GL.data) # Depth to sea floor below sea level (ft)
repl_int = f2m(L30.start) - f2m(L30.well.KB.data) + f2m(L30.well.GL.data)
water_vel = 1480 # velocity of sea water [m/s]
EGL_time = 2.0 * np.abs(f2m(L30.well.KB.data)) / water_vel
print "Ground Level time above SRD :", EGL_time
water_twt = 2.0*abs(water_depth + EGL_time) / water_vel
print "water_time: ", water_twt
print "Top of sonic log [m]: ", f2m(L30.start) # top of log (ft) (actually 1150 ft)
print "replacement interval [m]: ", repl_int
repl_vel = 1600 # m/s
repl_time = 2.0 * repl_int / repl_vel
print "two-way-replacement time: ", repl_time
log_start_time = water_twt + repl_time
print 'log_start_time:', log_start_time
def tvdss(md):
"assumes a vertical well"
return md - f2m(L30.well.KB.data)
top_log_TVDss = f2m(L30.well.KB.data) - f2m(L30.well.GL.data)
The starting time for our log is :
In [9]:
print log_start_time, ' seconds'
We would like to clip spikes. But np.clip() clips absolute values. This won't work for us because there's a trend in the data. Instead, we'd like to compare to the neighbourhood, and clip locally spurious things.
This is going to be a bit like, but not exactly like, the conservative filter.
In [10]:
def rolling_window(a, window):
shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
strides = a.strides + (a.strides[-1],)
rolled = np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
return rolled
Now we can smooth with a rolling median filter (preserves edges, only attacks spikes). Note we are not actually going to use the smoothed version, but we are going to use the diff to decide where the spikes are.
In [11]:
window = 13 # the length of filter is 13 samples or ~ 2 metres
rho_sm = np.median(rolling_window(RHOB,window), -1)
rho_sm = np.pad(rho_sm, window/2, mode='edge')
In [12]:
plt.figure(figsize=(18,4))
plt.plot(z,rho_sm,'b', z, RHOB, 'k', alpha=0.5)
plt.plot(z, rho_sm)
plt.show()
Now we will find spikes by comparing the original data with the median-smoothed data. If they differ by more than 100 (arbitrary), then it's a spike. Alternately, we could use the DRHO curve as a threshold for this:
In [13]:
def despike(curve, curve_sm, max_clip):
spikes = np.where(curve - curve_sm > max_clip)[0]
spukes = np.where(curve_sm - curve > max_clip)[0]
out = np.copy(curve)
out[spikes] = curve_sm[spikes] + max_clip # Clip at the max allowed diff
out[spukes] = curve_sm[spukes] - max_clip # Clip at the min allowed diff
return out
Now we know the indices of the spiky data, we can replace with a clipped version of the same spike.
In [14]:
rho = despike(RHOB, rho_sm, max_clip = 100)
Let's examine a segment of the log to see what effect it had,
In [15]:
start = 13000
end = 14500
plt.figure(figsize=(18,4))
plt.plot(z[start:end], RHOB[start:end],'k')
plt.plot(z[start:end], rho_sm[start:end],'b')
plt.plot(z[start:end], rho[start:end],'r')
plt.title('de-spiked density')
Out[15]:
Do the same thing for the sonic log,
In [16]:
dt_sm = np.median(rolling_window(DT,window), -1)
dt_sm = np.pad(dt_sm, window/2, mode='edge')
dt = despike(DT, dt_sm, max_clip = 10)
In [17]:
start = 13000
end = 14500
plt.figure(figsize=(18,4))
plt.plot(z[start:end], DT[start:end],'k')
plt.plot(z[start:end], dt_sm[start:end],'b')
plt.plot(z[start:end], dt[start:end],'r')
plt.title('de-spiked sonic')
Out[17]:
The time-to-depth relationship is obtained by scaling the sonic log by the sample interval (6" or 0.1524 m) and by calling np.cumsum() on it.
In [18]:
# two-way-time to depth relationship
scaled_dt = 0.1524 * np.nan_to_num(dt) / 1e6
tcum = 2 * np.cumsum(scaled_dt)
tdr = tcum + log_start_time
In [19]:
Z = (1e6/dt) * rho
In [20]:
RC = (Z[1:] - Z[:-1]) / (Z[1:] + Z[:-1])
Let's make a dictionary from the file 'tops.txt'
In [21]:
# Get well tops for plot annotation
tops = {}
with open('tops.txt') as f:
for line in f.readlines():
if not line.startswith('#'):
temp = line.strip().split('\t')
tops[temp[-1].replace('_',' ')] = float(temp[1])
tops
Out[21]:
In [22]:
tops.items() , tops.values()
Out[22]:
Find top relative to nearest sample in time domain
In [23]:
def find_nearest(array, value):
idx = (np.abs(array - value)).argmin()
return idx
For plotting tops in the time-domain, make a second entry in the dictionary for the travel-time.
In [24]:
tops_twt = {}
for key, val in tops.iteritems():
tops_twt[key] = tdr[find_nearest(z, val)]
tops_twt
Out[24]:
Plot the logs in the depth domain
In [25]:
f1 = plt.figure(figsize = (10,8))
ax1 = f1.add_axes([0.1, 0.1, 0.18, 0.8])
ax1.plot( DT, z,'b', alpha=0.5)
ax1.set_title('a) function', style = 'italic')
ax1.set_ylabel('measured depth ' + '$[m]$', fontsize = '12' )
ax1.set_xlabel('P-wave slowness '+ r'$[\mu s/m]$', fontsize = '12')
ax1.set_ylim(0, 4500)
ax1.set_xticklabels('')
ax1.invert_yaxis()
ax1.grid()
ax2 = f1.add_axes([0.3 , 0.1, 0.18, 0.8])
ax2.plot(tcum, z, 'b', alpha = 0.5)
ax2.set_title('integral', style = 'italic')
ax2.set_xlabel('two-way time ' + '$[s]$', fontsize = '12')
ax2.invert_yaxis()
ax2.set_yticklabels('')
ax2.set_xticklabels('')
ax2.set_ylim(4500, 0)
ax2.grid()
ax3 = f1.add_axes([0.5, 0.1, 0.18, 0.8])
ax3.plot( Z, z, 'k', alpha=0.5)
ax3.set_title('impedance', style = 'italic')
ax3.set_xlabel(r'$kg/m^2s^2$', fontsize = '12')
ax3.invert_yaxis()
ax3.set_yticklabels('')
ax3.set_xticklabels('')
ax3.set_ylim(4500, 0)
ax3.grid()
ax4 = f1.add_axes([0.7, 0.1, 0.18, 0.8])
ax4.plot( RC, z[:-1], 'k', alpha=0.5)
ax4.set_title('derivative', style = 'italic')
ax4.set_xlabel('reflectivity', fontsize = '12')
ax4.invert_yaxis()
ax4.set_yticklabels('')
ax4.set_xticklabels('')
ax4.set_ylim(4500, 0)
ax4.grid()
for i in range (4):
for top in tops.values() :
f1.axes[i].axhline( y = float(top), color = 'b', lw = 1,
ls = ':',
alpha = 0.5, xmin = 0.05, xmax = 0.95 )
for top, depth in tops.iteritems():
ax4.text( x = max(ax4.xaxis.get_data_interval())*1.0,
y = float(depth), s = top,
alpha=0.75, color='k',
fontsize = '10',
horizontalalignment = 'left',
verticalalignment = 'center',
bbox=dict(facecolor='white', alpha=1.0, lw = 0.25),
weight = 'light')
ax5 = f1.add_axes([0.8, 0.1, 0.18, 0.8])
ax5.invert_yaxis()
ax5.set_ylim(4500, 0)
ax5.set_axis_off()
But we need the functions represented by travel time
In [26]:
# RESAMPLING FUNCTION
dt = 0.004
maxt = 3.0
t = np.arange(0, maxt, dt)
Z_t = np.interp(x = t, xp = tdr, fp = Z)
RC_t = (Z_t[1:] - Z_t[:-1]) / (Z_t[1:] + Z_t[:-1])
In [27]:
# Define a wavelet.
def ricker(f, length, dt):
t = np.linspace(-length / 2, (length-dt) / 2, length / dt)
y = (1. - 2.*(np.pi**2)*(f**2)*(t**2))*np.exp(-(np.pi**2)*(f**2)*(t**2))
return t, y
In [28]:
RC_t = np.nan_to_num(RC_t)
tw, w = ricker (f=25, length = 0.512, dt = 0.004)
synth = np.convolve(w, RC_t, mode='same')
In [29]:
f2 = plt.figure(figsize=[10,12])
ax1 = f2.add_axes([0.05, 0.1, 0.2, 0.9])
ax1.plot(Z, z,'k', alpha=0.75)
ax1.set_title('impedance')
ax1.set_ylabel('measured depth ' + '$[m]$', fontsize = '12' )
ax1.set_xlabel(r'$kg/m^2s^2$ ', fontsize = '16')
ax1.set_ylim(0, 4500)
ax1.set_xticks( [0.0e7, 0.5e7, 1.0e7, 1.5e7, 2.0e7 ] )
ax1.invert_yaxis()
ax1.grid()
ax2 = f2.add_axes([0.325, 0.1, 0.2, 0.9])
ax2.plot(Z_t, t,'k', alpha=0.75)
ax2.set_title('impedance')
ax2.set_ylabel('two-way time ' + '$[s]$', fontsize = '12' )
ax2.set_xlabel(r'$kg/m^2s^2$ ', fontsize = '16')
ax2.set_ylim(0, 3)
ax2.set_xticks( [0.0e7, 0.5e7, 1.0e7, 1.5e7, 2.0e7 ] )
ax2.invert_yaxis()
ax2.grid()
ax3 = f2.add_axes([0.675, 0.1, 0.1, 0.9])
ax3.hlines(t[:-1], 0, RC_t, color='k', lw = 1) # Stems
ax3.plot([0, 0], [t.min(), t.max()], '-', c='k', alpha = 0.5) # Middle bar
ax3.set_title('reflectivity')
ax3.set_xlabel('', fontsize = '10')
ax3.set_ylim(0, 3)
ax3.set_xlim(-0.5, 0.5)
ax3.invert_yaxis()
ax3.set_yticklabels('')
ax3.set_xticks([-0.3, 0, 0.3] )
ax3.grid()
ax4 = f2.add_axes([0.8, 0.1, 0.2, 0.9])
ax4.plot(synth, t[:-1],'k')
ax4.fill_betweenx(t[:-1], synth, 0, synth > 0.0, color='k', alpha = 1.0)
ax4.set_title('synthetic')
ax4.set_xlabel('', fontsize = '10')
ax4.set_ylim(0, 3)
ax4.set_xlim(-0.05, 0.05)
ax4.invert_yaxis()
ax4.set_yticklabels('')
ax4.set_xticks([-0.6, -0.3, 0, 0.3, 0.6 ] )
ax4.grid()
for i in range(1):
for top, depth in tops.iteritems():
f2.axes[i].axhline( y = float(depth), color = 'b', lw = 2,
alpha = 0.5, xmin = 0.05, xmax = 0.95 )
f2.axes[i].text( x = 1e7, y = float(depth)-0.015, s = top,
alpha=0.75, color='k',
fontsize = '12',
horizontalalignment = 'center',
verticalalignment = 'center',
bbox=dict(facecolor='white', alpha=0.5, lw = 0.5),
weight = 'light')
for i in range(1,4):
for twt in tops_twt.values():
f2.axes[i].axhline( y = float(twt), color = 'b', lw = 2,
alpha = 0.5, xmin = 0.05, xmax = 0.95)
for i in range(1,2):
for top, twt in tops_twt.iteritems():
f2.axes[i].text( x = 2.75e7, y = float(twt), s = top,
alpha=0.75, color='k',
fontsize = '12',
horizontalalignment = 'center',
verticalalignment = 'center',
bbox=dict(facecolor='white', alpha=0.5, lw = 0.5),
weight = 'light'
)
real data is located at trace 77 of 500
In [30]:
traces = np.loadtxt('PenobXL_1155.txt')
# rearrange traces for plotting
traces = np.fliplr(np.transpose(traces))
In [31]:
import matplotlib.gridspec as gridspec
f3 = plt.figure(figsize=[16,12])
gs = gridspec.GridSpec(1, 6, width_ratios = [1.25, 1.25, 0.5, 1.25, 0.75, 4])
# depth domain
axa = plt.subplot(gs[0])
axa.plot( DT, z,'k', alpha=0.75, lw=0.25)
axa.set_title('travel time')
axa.set_ylabel('measured depth ' + '$[m]$', fontsize = '12' )
axa.set_xlabel('slowness 'r'$[\mu s/m]$', fontsize = '12')
axa.set_ylim(0, 4500)
axa.set_xticklabels('')
axa.invert_yaxis()
axa.grid()
gs.update(wspace=0.1)
axb = plt.subplot(gs[1])
axb.plot(tcum, z, 'k', alpha = 0.75)
axb.set_title('integrated travel time')
axb.set_xlabel('two-way time ' + '$[s]$', fontsize = '12')
axb.invert_yaxis()
axb.set_yticklabels('')
axb.set_xticklabels('')
axb.set_ylim(4500, 0)
axb.grid()
for top, depth in tops.iteritems():
axa.axhline( y = float(depth), color = 'b', lw = 2,
alpha = 0.5, xmin = 0.05, xmax = 0.95 )
axb.axhline( y = float(depth), color = 'b', lw = 2,
alpha = 0.5, xmin = 0.05, xmax = 0.95 )
axb.text( x = 1.25, y = float(depth)-0.015, s = top,
alpha=0.75, color='k',
fontsize = '10',
horizontalalignment = 'center',
verticalalignment = 'center',
bbox=dict(facecolor='white', alpha=1.0, lw = 0.5),
weight = 'light')
# time domain
#white space between depth and time plots
axoff = plt.subplot(gs[2])
axoff.set_axis_off()
axc = plt.subplot(gs[3])
axc.plot( Z_t, t, 'k', alpha=0.75)
axc.set_title('impedance')
axc.set_ylabel('two-way time '+ '$[s]$', fontsize = '12' )
axc.set_xlabel(r'$kg / m^2s^2$ ', fontsize = '12')
axc.set_ylim(0, 3)
axc.set_xlim(0, 2e7)
axc.invert_yaxis()
axc.grid()
for top, twt in tops_twt.iteritems():
axc.axhline( y = float(twt), color = 'b', lw = 2,
alpha = 0.5, xmin = 0.05, xmax = 0.95)
axc.text(x = 1e7, y = float(twt), s = top,
alpha=0.75, color='k',
fontsize = '10',
horizontalalignment = 'center',
verticalalignment = 'center',
bbox=dict(facecolor='white', alpha=0.5, lw = 0.5),
weight = 'light')
axd = plt.subplot(gs[4])
axd.hlines(t[:-1], 0, RC_t, color='k', lw = 1) # Stems
axd.plot([0, 0], [t.min(), t.max()], '-', c='k', alpha = 0.5) # Middle bar
axd.set_title('reflectivity')
axd.set_xlabel('', fontsize = '12')
axd.set_ylim(0, 3)
axd.set_xlim(-0.5, 0.5)
axd.invert_yaxis()
axd.set_yticklabels('')
axd.set_axis_off()
axd.grid()
for depth in tops_twt.values():
axd.axhline(y = float(depth), color = 'b', lw = 2,
alpha = 0.5, xmin = 0.25, xmax = 0.75)
axe = plt.subplot(gs[5])
axe.imshow( traces[:750,100:], cmap = 'Greys',
vmin = -7000, vmax = 7000, alpha = 0.8,
aspect = 'auto'
)
axe.set_axis_off()
axe.set_title('seismic')
# put synthetic on seismic
bottom = axe.get_position().get_points()[0][1]
top = axe.get_position().get_points()[1][1]
axf = axe.figure.add_axes([0.79, bottom, 0.05, top-bottom])
gain_synth = 1
axf.plot(gain_synth * synth, t[:-1], 'k', alpha = 0.9)
axf.fill_betweenx(t[:-1], gain_synth * synth, 0,
gain_synth * synth > 0.0,
color = 'k', alpha = 0.5)
axf.set_title('synthetic')
axf.set_xlabel('', fontsize = '12')
axf.set_ylim(0, 3)
axf.set_xlim(-0.5 ,0.5)
axf.invert_yaxis()
axf.set_yticklabels('')
axf.set_axis_off()
axf.grid()
for top, twt in tops_twt.iteritems():
axf.axhline( y = float(twt), color = 'b', lw = 2,
alpha = 0.25, xmin = 0.25, xmax = 0.75)
axf.text(x = -0.4, y = float(twt), s = top,
alpha=0.75, color='k',
fontsize = '10',
horizontalalignment = 'right',
verticalalignment = 'center',
bbox=dict(facecolor='white', alpha=0.25, lw = 0.5),
weight = 'light')
#f3.savefig('Figure_1.png', facecolor = 'white', dpi = 300)