In [1]:
%pylab inline
pylab.rcParams['figure.figsize'] = (10.0, 8.0)
import toolbox
import numpy as np
from IPython.display import HTML
print dir(toolbox)
import su
Seismic exploration uses seismic waves to explore the subsurface. The resulting raw data does not have much in common with the actual subsurface, and also has unwanted signals (noise). Seismic processing attempts to remove noise, and create a seismic image which can be related to the geology of the subsurface.
The kind of model we are going to build is based upon the "Convolutional Model"
This is only one type of model. Another type of model might be a finite difference model, which uses the partial differential wave equation to model seismic waves propogating through a medium.
The convolutional model of the seismic trace states that the trace we record is the result of the earth's reflectivity (what we want) convolved with the source wavelet, the recording system and some noise.
ie.
$$E * S * R + N$$where
E = earth reflection series
S = source wavelet
R = receiver function
N = noise
* = convolution
In [2]:
HTML('<iframe src=http://en.wikipedia.org/wiki/Convolution#Definition width=1000 height=240</iframe>')
Out[2]:
To convolve is to roll together, so convolution can be though of as "rolling" two signals together, and is usually denoted by (*).
I'm building the source wavelet by band-limiting a spike in the frequency domain. The result is known as a "klauder wavelet", a wavelet which is produced when you use vibroseis as a source.
There is some good geophysical signal processing theory involved in this, but I dont expect you to write this code, the function
wavelet = build_wavelet(low,high)
has been supplied in the file toolbox.py
In [3]:
from toolbox import build_wavelet
low = 5. #hz
high = 60. #hz
wavelet = np.convolve(np.ones(3)/3., build_wavelet(low, high)[75:125], mode='same')
xaxis = np.arange(50)/1000.
pylab.plot(xaxis, wavelet)
pylab.ylabel('amplitude')
pylab.xlabel('time')
pylab.title('Line plot of klauder wavelet')
pylab.figure()
pylab.plot(xaxis, wavelet, 'bo')
pylab.ylabel('amplitude')
pylab.xlabel('time')
pylab.title('Discrete plot of klauder wavelet')
Out[3]:
Seismic waves travel from the source, bounce off the reflection point and travel back up to the geophone.
So the two things we need to build our synthetic record is the time from the source to the refelction point to the geophone (two way travel time) and the amplitude of the wave at the geophone.
In [4]:
def traveltime(sx, gx, v, z):
h = gx - sx
v = float(v)
distance = np.sqrt(h*h + z*z)
time = 2.0*distance/v
return time
In [5]:
def diverge(distance):
return np.abs(1.0/(distance*distance))
The ratio of amplitude of the reflected wave to the incident wave, or how much energy is reflected.
If the wave has normal incidence, then its reflection coefficient can be expressed as:
$$ R = \frac{(\rho_2V_2 - \rho_1V_1)}{(\rho_2V_2 + \rho_1V_1)} $$$\rho V$ is often referred to as the accoustic impedance "Z".
One regular use is to generate seismograms
At non-normal incidence, the reflection coefficient defined as a ratio of amplitudes depends on other parameters, such as the shear velocities, and is described as a function of incident angle by the Zoeppritz equations.
In [6]:
HTML('<iframe src=http://subsurfwiki.org/wiki/Zoeppritz_equation width=1000 height=800</iframe>')
Out[6]:
We can examine the effect of angle of incidence using the CREWES Zoeppritz Explorer.
The Zoeppritz equations are fairly complex to implement and not very intuitive. There are a number of approximations available. You can see the Aki-Richards approximation in action in the CREWES Zoeppritz Explorer. An even simpler approximation is the Shuey equation. You are going to use the Shuey equation to calculate the amplitudes in the synthetic model.
Create a simple 2-layer model. Define Vp, Vs and density for each layer. Implement the Shuey Equation for a single shot-geophone pair. Plug in the same numbers into Zoeppritz explorer to confirm correctness.
In [7]:
def reflection_coefficient(z0, z1):
z = ((z1 - z0)**2)/(z1+z0)**2
if z1 < z0: z*= -1
return z
def transmission_coefficient(z0, z1):
return (4*z0*z1)/((z1+z0)**2)
Assume we have a single shot, and a line of geophones. The angle to each geophone will be different, resulting in a different solution to the Shuey Equations. Take the code from Exersize 1 and loop it over a range of geophones (still in the simple 2 layer case).
In [8]:
#define survey geometry, ie shot and reciever points
sx_coords = [120] #np.arange(99)+0.5
rx_coords = np.arange(500)
#define an output array
nx = 500 #number of samples in x direction (assume 1m cells)
ns = 1000 #number of samples in z direction (assume 1ms cells)
output = np.zeros((nx, ns), 'f')
In [9]:
import toolbox
model = toolbox.build_model()
pylab.imshow(model['model']['vp'].T, aspect='auto')
Out[9]:
In [10]:
direct = np.zeros_like(output)
#direct. assume speed of sound (330m/s)
v = 330.0
x_locs = np.abs(np.arange(500) - sx_coords[0])
t_locs = (x_locs/v)*1000
t = np.floor(t_locs[t_locs < 1000]).astype(np.int)
x = np.arange(500)[t_locs < 1000].astype(np.int)
amps = np.ones_like(x, 'f')*diverge(x_locs[t_locs < 1000])
amps[np.argmax(amps)] = amps[np.argmax(amps)+1]
amps *= 0.05
direct[x,t]= amps
#pylab.imshow(output.T, aspect='auto')
sutype = su.typeSU(1000)
data = np.zeros(500, dtype=sutype)
data['gx'] = range(1, 501)
data['tracl'] = range(1,501)
data['sx'] = 120
data['offset'] = data['gx'] - data['sx']
data['cdp'] = (data['sx']+data['gx'])/2
data['trace'] = direct.astype(np.float32)
data['ns'] = 1000
data['dt'] = 1000
su.writeSU(data, "direct.su")
In [11]:
#refraction
refraction = np.zeros_like(output)
def refract(x, v0, v1, z0):
ic = np.arcsin(v0/v1)
t0 = 2.0*z0*np.cos(ic)/v0
t = t0 + x/v1
return t
z0 = model['dz'][0]
v0 = model['vp'][0]
v1 = model['vp'][1]
x_locs = np.abs(np.arange(500) - sx_coords[0])
tt = np.floor((refract(x_locs, v0, v1, z0)*1000.0)).astype(np.int)
xx = np.arange(500).astype(np.int)
ampr = np.ones_like(xx, 'f')*diverge(x_locs[tt < 1000])
ampr[np.argmax(ampr)] = 0
ampr *= 0.1
refraction[xx,tt] += ampr
tmin = np.ceil((z0*2.2/v0)*1000).astype(np.int)
refraction[:,:tmin] = 0
#pylab.imshow(refraction.T, aspect='auto')
sutype = su.typeSU(1000)
data = np.zeros(500, dtype=sutype)
data['gx'] = range(1, 501)
data['tracl'] = range(1,501)
data['sx'] = 120
data['offset'] = data['gx'] - data['sx']
data['cdp'] = (data['sx']+data['gx'])/2
data['trace'] = refraction.astype(np.float32)
data['ns'] = 1000
data['dt'] = 1000
su.writeSU(data, "refraction.su")
In [12]:
#calculate traveltimes and amplitudes for each reflection point
num = 100
vp = model['model']['vp']
rho = model['model']['rho']
R = model['model']['R']
reflection = np.zeros_like(output)
gxs = np.arange(0,500)
sx = 120
sz = 0
gz = 0
print R.shape
for gx in gxs:
cmpx = np.floor((gx + sx)/2.).astype(np.int)
h = cmpx - sx
for cmpz in (np.nonzero(R[cmpx,:])[0]):
ds = np.sqrt(cmpz**2 + (h)**2)/float(num) # line step distance
#predefine outputs
amp = 1.0
time = 0.0
#traveltime from source to cdp
vp_down = toolbox.find_points(sx, sz, cmpx, cmpz, num, vp)
time += np.sum(ds/vp_down)
#traveltime from cdp to geophone
vp_up = toolbox.find_points(cmpx, cmpz, gx, gz, num, vp)
time += np.sum(ds/vp_up)
#loss due to spherical divergence
amp *= diverge(ds*num*2)
#transmission losses from source to cdp
rho_down = toolbox.find_points(sx, sz, cmpx, cmpz, num, rho)
z0 = rho_down * vp_down
z1 = np.roll(z0, shift=1)
z1[0] = z0[0]
z1[-1] = z0[-1]
correction = np.cumprod(transmission_coefficient(z0, z1))[-1]
amp *= correction
#amplitude loss at reflection point
correction = R[cmpx,cmpz]
amp *= correction
#transmission loss from cdp to source
rho_up = toolbox.find_points(cmpx, cmpz, gx, gz, num, rho)
z0 = rho_up * vp_up
z1 = np.roll(z0, shift=-1)
z1[0] = z0[0]
z1[-1] = z0[-1]
correction = np.cumprod(transmission_coefficient(z0, z1))[-1]
amp *= correction
sample = np.floor(time*1000).astype(np.int)
gx = np.floor(gx).astype(np.int)
reflection[gx, sample] += amp
#pylab.figure()
#pylab.imshow(record[:,:].T, aspect='auto')
#colorbar()
sutype = su.typeSU(1000)
data = np.zeros(500, dtype=sutype)
data['gx'] = range(1, 501)
data['tracl'] = range(1,501)
data['sx'] = 120
data['offset'] = data['gx'] - data['sx']
data['cdp'] = (data['sx']+data['gx'])/2
data['trace'] = reflection.astype('f4')
data['ns'] = 1000
data['dt'] = 1000
su.writeSU(data, "reflection.su")
In [13]:
#check traveltimes d/v
#print 1000*20./800
#print record[np.nonzero(record)]
In [14]:
noise = np.random.normal(0.0, 1e-9, size=(output.shape))
output = direct + refraction + reflection + noise
output /= np.amax(output)
record = np.apply_along_axis(lambda m: np.convolve(m, wavelet, mode='same'), axis=1, arr=output)
In [15]:
sutype = su.typeSU(1000)
data = np.zeros(500, dtype=sutype)
data['gx'] = range(1, 501)
data['tracl'] = range(1,501)
data['sx'] = 120
data['offset'] = data['gx'] - data['sx']
data['cdp'] = (data['sx']+data['gx'])/2
data['trace'] = record.astype('f4')
data['ns'] = 1000
data['dt'] = 1000
data = data[::2]
su.writeSU(data, "record.su")
In [18]:
#top of coals will be 2*(time(a) + time(b))
sx = 120
gx = np.arange(122,300)
z0, z1, z2, z3 = 20., 40., 20., 100.
v0, v1, v2, v3 = 800., 2200., 1800., 2400.
h = (gx - sx)/2.
i1 = np.arctan(h/(z0+z1))
l_a = z0/np.cos(i1)
l_b = z1/np.cos(i1)
t_ab = (l_a/v0) + (l_b/v1)
tt_1 = t_ab*2
i2 = np.arctan(h/(z0+z1+z2))
l_c = z0/np.cos(i2)
l_d = z1/np.cos(i2)
l_e = z2/np.cos(i2)
t_cde = (l_c/v0) + (l_d/v1) + (l_e/v2)
tt_2 = t_cde*2
plot(tt_1)
plot(tt_2)
ylim([0.4,0])
print tt_1[-1], tt_2[-1]