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)

Done!

**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]

you use `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 `earth`

.

Pingback: Adding normal and reverse faults to an impedance model with Python | MyCarta·

hello, i want to ask about the fault on the figure.

1. is there a code to make to fault not vertical?

2. is there a code to count the slip of the fault?

thank you

Hello Phillip

I show how to make a fault at ~45 degrees in here: https://mycarta.wordpress.com/2017/10/03/adding-normal-and-reverse-faults-to-an-impedance-model-with-python/

Knowing the shift in samples, you could convert it to meters.

For a more general solution to non-vertical faults, I would look at Pynoddy: https://pynoddy.readthedocs.io/en/latest/notebooks/4-Create-model.html

Hi Matteo,

very useful notebook, I was looking for something similar, but I’m still struggling to save the models to bin files, I tried to do this:

earth = np.asarray(earth, order=’C’)

for w in range(1,11):

filename = “figures/model{}.bin”.format(w)

with open(filename, ‘wb’) as outfile:

outfile.write(earth)

It seems to save the earth models, but when I use Seismic Unix to open the bin files, that should have as first dimension =50, it give me only a black figure with zeros. Is it wrong what I’m doing?

Hi Flavio

Unfortunately I’ve never worked with that format, or Seismic Unix for what matters.for that matter. But I can almost guarantee if you ask in the Software Underground Slack group someone will come up with something…

https://softwareunderground.org/slack

Thanks Matteo,

I’ll ask there.