import matplotlib.pyplot as plt import numpy as np % matplotlib inline
We'll start off with an earth model --- an array of 'cells', each of which has some rock properties.
Line 1 sets up some basic variables, then in line 2 I've used a little matrix-forming trick,
np.tri(m, n, k), which creates an m × n matrix with ones below the kth diagonal, and zeros above it. The
dtype specification just makes sure we end up with integers, which we need later for the indexing trick.
Then line 3 just sets every row above
// is integer division, because NumPy prefers integers for indexing arrays), to 0.
length, depth = 40, 100 model = 1 + np.tri(depth, length, -depth//3, dtype=int) model[:depth//3,:] = 0
We'll have a quick look with some very basic plotting commands.
plt.imshow(model, cmap='viridis', aspect=0.2) plt.show()
Now we can make some Vp-rho pairs (rock 0, rock 1, rock 2) and select from those with
np.take. This works like
vlookup in Excel --- it says "read this array,
model in this case, in which the values i are like 0, 1, ... n, and give me the ith element from this other array,
rocks in this case.
rocks = np.array([[2700, 2750], # Vp, rho [2400, 2450], [2800, 3000]])
Edit: 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 = rocks[model]
np.product to those Vp-rho pairs to get impedance at every sample.
This might look a bit magical, but we're just telling Python to apply the function
product() to every set of numbers it encounters on the last axis (index
-1). The array
earth has shape (100, 40, 2), so you can think of it as a 100 row x 40 column 'section' in which each 'sample' is occupied by a Vp-rho pair. That pair is in the last axis. So product, which just takes a bunch of numbers and multiplies them, will return the impedance (the product of Vp and rho) at each sample location. We'll end up with a new 100 x 40 'section' with impedance at every sample.
imp = np.apply_along_axis(np.product, -1, earth)
We could have saved a step by taking from
np.product(rocks, axis=1) but I like the elegance of having an earth model with a set of rock properties at each sample location. That's how I think about the earth --- and it's similar to the concept of a geocellular model.
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.
I love this indexing trick though I admit it looks weird the first time you see it. It's easier to appreciate for a 1D array. Let's look at the differences:
>>> a = np.array([1,1,1,2,2,2,3,3,3]) >>> a[1:] - a[:-1] array([0, 0, 1, 0, 0, 1, 0, 0])
This is equivalent to:
>>> np.diff(a, axis=0)
But I prefer to spell it out so it's analogous to the sum on the denominator.
rc = (imp[1:,:] - imp[:-1,:]) / (imp[1:,:] + imp[:-1,:]) plt.imshow(rc, cmap='Greys', aspect=0.2) plt.show()
We'll use a wavelet function from
bruges. This is not cheating! Well, I don't think it is... we could use
scipy.signal.ricker but I can't figure out how to convert frequency into the 'width' parameter that function wants. Using the Ricker from
bruges keeps things a bit simpler.
import bruges w = bruges.filters.ricker(duration=0.100, dt=0.001, f=40)
Let's make sure it looks OK:
Now one more application of
apply_along_axis. We could use a loop to step over the traces, but the rule of thumb in Python is "if you are using a loop, you're doing it wrong.". So, we'll use
It looks a bit more complicated this time, because we can't just pass a function like we did with
product before. We want to pass in some more things, not just the trace that
apply_along_axis is going to send it. So we use Python's 'unnamed function creator',
lambda (in keeping with all things called
lambda, it's a bad name that no-one can quite explain).
synth = np.apply_along_axis(lambda t: np.convolve(t, w, mode='same'), axis=0, arr=rc) plt.imshow(synth, cmap="Greys", aspect=0.2) plt.show()
That's it! And it only needed 9 lines of Python! Not incldung boring old imports and plotting stuff.
Here they are so you can count them:
length, depth = 40, 100 model = 1 + np.tri(depth, length, -depth//3) model[:depth//3,:] = 0 rocks = np.array([[2700, 2750], [2400, 2450], [2800, 3000]]) earth = np.take(rocks, model.astype(int), axis=0) imp = np.apply_along_axis(np.product, -1, earth) rc = (imp[1:,:] - imp[:-1,:]) / (imp[1:,:] + imp[:-1,:]) w = bruges.filters.ricker(duration=0.100, dt=0.001, f=40) synth = np.apply_along_axis(lambda t: np.convolve(t, w, mode='same'), axis=0, arr=rc)
© Agile Geoscience 2016