In this short post I show how to adapt Agile Scientific‘s Python tutorial x lines of code, Wedge model and adapt it to make 100 synthetic models in one shot: X impedance models times X wavelets times X random noise fields (with I vertical fault).
You can download the notebook with the full code from GitHub. N.B. the code is optimized for Python 2.7.
I begin by making a 6-layer model, shown in Figure 1:
model = numpy.zeros((50,49), dtype=np.int) model[8:16,:] = 1 model[16:24,:] = 2 model[24:32,:] = 3 model[32:40,:] = 4 model[40:,:] = 5
next I make some Vp-rho pairs (rock 0, rock 1, … , rock5):
rocks = numpy.array([[2700, 2750], # Vp, rho [2400, 2450], [2600, 2650], [2400, 2450], [2800, 3000], [3100, 3200],])
and then create 10 slightly different variations of the Vp-rho pairs one of which are is shown in Figure 2:
rnd = numpy.random.rand(10,6,2)*0.2 manyrocks = np.array([rocks + rocks*rn for rn in rnd], dtype=np.int) earth = manyrocks[model]
at which point I can combine Vp-rho pairs to make 10 impedance models, then insert a vertical fault with:
impedances = [np.apply_along_axis(np.product, -1, e).astype(float) for e in earth]# Python 2 faulted = copy.deepcopy(impedances) for r, i in zip(faulted, np.arange(len(faulted))): temp = np.array(r) rolled = np.roll(np.array(r[:,:24]), 4, axis = 0) temp[:,:24]=rolled faulted[i]=temp
next I calculate reflection coefficients (Figure 4)and convolve them with a list of 10 Ricker wavelets (generated using Agile’s Bruges) to make synthetic seismic models, shown in Figure 5.
rc = [(flt[1:,:] - flt[:-1,:]) / (flt[1:,:] + flt[:-1,:]) for flt in faulted] ws = [bruges.filters.ricker(duration=0.098, dt=0.002, f=fr) for fr in [35, 40, 45, 50, 55, 60, 65, 70, 75, 80]] synth = np.array([np.apply_along_axis(lambda t: np.convolve(t, w, mode='same'), axis=0, arr=r) for r in rc for w in ws ])
The last bit is the addition of noise, with the result is shown in Figure 6:
blurred = sp.ndimage.gaussian_filter(synth, sigma=1.1) noisy = blurred + 0.5 * blurred.std() * np.random.random(blurred.shape)
The notebook with the full code is on GitHub, let me know if you find this useful or if you come up with more modeling ideas.
UPDATE, April 22, 2019.
I received an email from a reader, Will, asking some clarification about this article and the making of many impedance models. I’m re-posting here what I wrote back.
I think the key to understand this is how we multiply velocity by density in each of the possible earth model.
print (np.shape(earth)) >>> (10, 50, 49, 2)
impedances = [np.apply_along_axis(np.product, -1, e).astype(float) for e in earth]
numpy.apply_along_axis to multiply Vp by Rho along the last dimension,
-1 , using
numpy.product, and the list comprehension
[... for e in earth] to do it for all models in the array