EXERCISE — Synthetic seismic: wedge

Q. Try changing the rock properties.

Q. Try plotting unresolved and resolved traces.

Q. try using a different wavelet, e.g. spike, square, sinc, Ormsby, Butterworth, etc.

Q. Try making a different geometry.

Make an earth model



In [1]:

    
import matplotlib.pyplot as plt
import numpy as np
% matplotlib inline



In [116]:

    
import matplotlib as mpl

img = mpl.image.imread('../data/Hubbard_etal_2014.png')



In [118]:

    
meanimg = np.mean(img, axis=-1)

model = np.digitize(meanimg, bins=[0.5, 0.6, 0.7, 0.8])



In [119]:

    
np.unique(model)









    Out[119]:





array([0, 1, 2, 3, 4])



In [100]:

    
length = 80  # x range
depth = 200  # z range



In [91]:

    
model = 1 + np.tri(depth, length, -depth//3, dtype=int)

plt.imshow(model)
plt.show()



In [92]:

    
model[100]









    Out[92]:





array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

Now set the upper part of the model — above the wedge — to zero (blue).



In [93]:

    
model[:depth//3,:] = 0

plt.imshow(model)
plt.colorbar()
plt.show()

Now we can make some Vp-rho pairs (rock 0 and rock 1).



In [120]:

    
rocks = np.array([[2600, 2700],
                  [2700, 2750],
                  [2400, 2450],
                  [2550, 2650],
                  [2800, 3000]])

I was using np.take here, but 'fancy indexing' is shorter and more intuitive. We are just going to index rocks using the integers in model. That is, if model has a 1, we take the second element, [2400, 2450], from rocks. We'll end up with an array containing the rocks corresponding to each element of earth.



In [121]:

    
earth = rocks[model]

Now apply np.product to those Vp-rho pairs to get impedance at every sample.



In [122]:

    
imp = np.apply_along_axis(np.product, -1, earth)



In [123]:

    
plt.imshow(imp)
plt.colorbar()









    Out[123]:





<matplotlib.colorbar.Colorbar at 0x7f9699e365f8>

Model seismic reflections

Now we have an earth model — giving us acoustic impedance everywhere in this 2D grid — we define a function to compute reflection coefficients for every trace.



In [124]:

    
def make_rc(imp):
    upper = imp[ :-1, :]
    lower = imp[1:  , :]
    
    return (lower - upper) / (lower + upper)

rc = make_rc(imp)

plt.imshow(rc, aspect=1)
plt.show()

We'll use scipy.signal.ricker to make a wavelet:



In [125]:

    
from scipy.signal import ricker

f = 25
wavelet = ricker(100, a=5)

plt.plot(wavelet)
plt.show()



In [132]:

    
import bruges
w = bruges.filters.ricker(duration=0.1, dt=0.001, f=60)



In [133]:

    
plt.plot(w)









    Out[133]:





[<matplotlib.lines.Line2D at 0x7f9699295630>]

Now apply 1D convolution to every trace:



In [134]:

    
def convolve(trace):
    return np.convolve(trace, w, mode='same')

synth = np.apply_along_axis(convolve,
                            axis=0,
                            arr=rc)

plt.figure(figsize=(12,6))
plt.imshow(synth, cmap="RdBu", aspect=1)
plt.colorbar()
plt.show()

Measure amplitude at the peak



In [53]:

    
apparent = np.amin(synth, axis=0)
t_apparent = np.argmin(synth, axis=0)

actual = synth[depth//3, :]
t_actual = np.ones_like(t_apparent) * depth//3



In [54]:

    
plt.plot(apparent, label="apparent")
plt.plot(actual, label="actual")
plt.text(6, 0, "actual", color='g')
plt.text(1, -0.18, "apparent", color='b')
# plt.legend()
plt.show()



In [55]:

    
# t_apparent = t_apparent.astype(np.float)
# t_apparent[0] = np.nan



In [56]:

    
plt.figure(figsize=(12,6))

# Plot seismic
plt.imshow(synth, cmap="Greys", aspect=0.2)
plt.colorbar()

# Plot horizons
plt.plot(t_apparent, 'r')
plt.plot(t_actual, 'c')

plt.xlim([0, 79])
plt.ylim([199, 0])
plt.title("Wedge model, {} Hz wavelet".format(f))
plt.show()

Find the tuning thickness

We can find the thickness at which the amplitude peak starts... this is called the tuning thickness.



In [58]:

    
tuning = np.argmin(actual)
vp = 2400
wavelength = 2 * 1000 * (vp / f) / vp  # in milliseconds TWT

print("Tuning thickness is {:.1f} ms, or 1/{:.1f} wavelength".format(tuning, wavelength/tuning))









    



Tuning thickness is 18.0 ms, or 1/4.4 wavelength



In [ ]:



In [ ]: