# Seismic Imaging Blind Challenge¶

## Resources and Code Setup¶

The file cp9.py is provided to you at the computing resources website. This file implements a significant portion of the code you will need for today.

You will also need to download the seismic data for the challenge problem. cp9.py has the code required to load it. Be sure to put it in the same directory as your cp9.py script.

In cp9.py there is a section with comments describing how to set up an initial velocity model. Tunable parameters are marked with #### TUNABLE PARAMETER #### comments. Helpful plotting code is marked with #### PLOT ####.

In [1]:
import time
import math

import numpy as np
import matplotlib.pyplot as plt

from pysit import *
from pysit.gallery import horizontal_reflector

#	Define Domain
pmlx = PML(0.1, 100)
pmlz = PML(0.1, 100)

x_fac = 6
z_fac = 3
x_config = (0.0, 0.35*x_fac, pmlx, pmlx)
z_config = (0.0, 0.50*z_fac, pmlz, pmlz)

domain = RectangularDomain(x_config, z_config)

nx = 35*x_fac
nz = 50*z_fac
mesh = CartesianMesh(domain, nx, nz)
XX, ZZ = mesh.mesh_coords()

shot = shots[0]

# Define and configure the wave solver
trange = shot.trange

solver = ConstantDensityAcousticWave(mesh,
formulation='scalar',
model_parameters={'C': np.ones(mesh.shape())},
spatial_accuracy_order=2,
trange=trange,
use_cpp_acceleration=True,)

objective = TemporalLeastSquares(solver)

(210L, 150L) (31500, 1)


## Questions¶

In this problem, you will attempt to recover an unknown medium from the provided seismic data.

Use the theoretical and computational tools (PySIT) you have developed and used over the last two weeks to attempt to answer the following questions about the unknown domain.

Be sure to look at the details section for important information.

Question 1: What is the background velocity of the medium? How did you come to this conclusion? What indications do you have that you are correct?

Answer 1: The background velocity is $c_0(x) = 1 - 0.05\pi \approx 0.8429203673205103$, which would be very difficult to recover on its own. There are a few ways that an approximation could be recovered.

The first method is examining the gradient for a few different background velocities. We know from the given information that we have a horizontally layered medium. Given that the source is in the center, for the correct constant background we expect the gradient show reflectors as flat oscillitory things where there is sufficient aperture. When the velocity is too slow, these reflectors will appear to "smile" and when the velocity is too fast, they will "frown".

For example, for $c_0(x) = 0.65$, we have:

In [2]:
C0 = 0.65*np.ones_like(ZZ)
m0 = solver.ModelParameters(mesh,{'C': C0})

plt.figure()
vis.plot(g.data, mesh, clim=(-0.005,0.005))
plt.colorbar();


For $c_0(x) = 1.1$, we have:

In [3]:
C0 = 1.1*np.ones_like(ZZ)
m0 = solver.ModelParameters(mesh,{'C': C0})

plt.figure()
vis.plot(g.data, mesh, clim=(-0.01,0.01))
plt.colorbar();


Continuing with a binary search...

In [4]:
C0 = 0.875*np.ones_like(ZZ)
m0 = solver.ModelParameters(mesh,{'C': C0})

plt.figure()
vis.plot(g.data, mesh, clim=(-0.01,0.01))
plt.colorbar();

In [5]:
C0 = 0.7625*np.ones_like(ZZ)
m0 = solver.ModelParameters(mesh,{'C': C0})

plt.figure()
vis.plot(g.data, mesh, clim=(-0.01,0.01))
plt.colorbar();

In [6]:
C0 = 0.81875*np.ones_like(ZZ)
m0 = solver.ModelParameters(mesh,{'C': C0})

plt.figure()
vis.plot(g.data, mesh, clim=(-0.01,0.01))
plt.colorbar();


I'd be fairly satisfied with this background: $c_0(x) = 0.81875$.

Another approach would be to use travel time tomography. This would normally be more difficult, but because you know that the medium is constant, you can look directly at the data and identify features in each one. For example, comparing the times at which the maximum point in two traces is obtained to the distance between the receivers should be a good approximation. Due to numerical dispersion, we must pick nearby receivers.

Lets take receiver 150 and 155 and (to make the visualization more useful) look at the first 2500 data points.

In [7]:
plt.figure()
plt.plot(shot.receivers.data[:2500,155], label="Trace 155") # ~middle index
plt.legend();

In [8]:
x = (155-150)*mesh.x.delta
print max_150, max_155
t = (max_155-max_150)*shot.dt
print x, t
print "Estimated speed: ", x/t

930 1011
0.05 0.0593550652581
Estimated speed:  0.84238808908


The estimated velocity is very close to the true velocity $c_0(x) = 0.8429203673205103$. This is good enough, though you could average over a number of small windows to come up with something more accurate. This is exhibited in the function fit_travel_times in the solution file, which produces $c_o(x) = 0.84293408738$. To compare, the gradient looks like,

In [9]:
C0 = 0.84293408738*np.ones_like(ZZ)
m0 = solver.ModelParameters(mesh,{'C': C0})

plt.figure()
vis.plot(g.data, mesh, clim=(-0.01,0.01))
plt.colorbar();


which is only barely better than the crude search approach.

Note that I deliberately oversampled the data so cannot use the data's $\Delta t$ to back out the maximum velocity from the model, but you could use this for a bound.

Question 2: How many reflectors are there and what is their location? How did you come to this conclusion?

Answer 2: When the background velocity is correct, the position and number of reflectors is evident. At first glance, it appears that there are three reflectors, but when we re-examine the scale, we see at more appear.

In [10]:
plt.figure()
vis.plot(g.data, mesh, clim=(-0.01,0.01))
plt.colorbar();


This can also be seen by plotting the center depth slice.

In [11]:
from cp9_sol import reshape_pysit_vector_to_grid

plt.figure()
g_reshaped = reshape_pysit_vector_to_grid(g.data, mesh)
plt.plot(np.arange(mesh.z.n)*mesh.z.delta,g_reshaped[mesh.x.n/2, :])
plt.ylim(-0.01, 0.01); # rescale to be able to see things

Out[11]:
(-0.01, 0.01)

It turns out that most of these wiggles in the image are due to 'multiples' which result in 'ghost' reflectors. There are in fact three reflectors, at $z = 0.6, 0.9,$ and $1.05$, though it is not 100% clear at this point.

Question 3: What is the amplitude of the reflectors? How did you come to this conclusion?

Answer 3: After running a few iterations of gradient descent, you can get a reasonable image.

In [12]:
invalg = LBFGS(objective)

print('Running FWI...')
tt = time.time()

nsteps = 25
status_configuration = {'value_frequency'		    : 1,  # m_k
'residual_frequency'		: 1,  # r_k
'residual_length_frequency' : 1,  # ||r_k||
'objective_frequency'	    : 1,  # J[m_k]
'step_frequency'			: 1,  # s_k
'step_length_frequency'	    : 1,  # ||s_k||
'run_time_frequency'		: 1,  # obvious
'alpha_frequency'		    : 1   # line search length
}

result = invalg(shots, m0, nsteps,
line_search='backtrack',
status_configuration=status_configuration, verbose=True)

print '...run time:  {0}s'.format(time.time()-tt)

Running FWI...
Iteration 0
objective 0.000700760976375
residual 0.0374369062925
Starting:  2718.63918202 0.000700760976375
Pass 1: a:0.663730269049; 0.000842213785596 ?<= 0.000699625438603
Pass 2: a:0.530984215239; 0.000565722833393 ?<= 0.000700034232201
alpha 0.530984215239
run time 26.118999958s
Iteration 1
objective 0.000565722833393
residual 0.0336369687514
Starting:  1.0 0.000565722833393
Pass 1: a:1.0; 0.000304862288474 ?<= 0.000565270867657
alpha 1.0
run time 22.3930001259s
Iteration 2
objective 0.000304862288474
residual 0.0246926016642
Starting:  1.0 0.000304862288474
Pass 1: a:1.0; 0.000282349553079 ?<= 0.000304832558139
alpha 1.0
run time 21.1879999638s
Iteration 3
objective 0.000282349553079
residual 0.0237633984556
Starting:  1.0 0.000282349553079
Pass 1: a:1.0; 0.000264244388033 ?<= 0.000282327570993
alpha 1.0
run time 18.3580000401s
Iteration 4
objective 0.000264244388033
residual 0.0229888837499
Starting:  1.0 0.000264244388033
Pass 1: a:1.0; 0.000249278564295 ?<= 0.000264225632682
alpha 1.0
run time 18.7809998989s
Iteration 5
objective 0.000249278564295
residual 0.0223283928797
Starting:  1.0 0.000249278564295
Pass 1: a:1.0; 0.000233098990097 ?<= 0.000249258137776
alpha 1.0
run time 19.0780000687s
Iteration 6
objective 0.000233098990097
residual 0.0215916182857
Starting:  1.0 0.000233098990097
Pass 1: a:1.0; 0.000215339082454 ?<= 0.000233052934141
alpha 1.0
run time 21.2669999599s
Iteration 7
objective 0.000215339082454
residual 0.020752786919
Starting:  1.0 0.000215339082454
Pass 1: a:1.0; 0.000195766918232 ?<= 0.000215306358199
alpha 1.0
run time 23.0469999313s
Iteration 8
objective 0.000195766918232
residual 0.0197872139642
Starting:  1.0 0.000195766918232
Pass 1: a:1.0; 0.000179244685589 ?<= 0.000195742977229
alpha 1.0
run time 23.6029999256s
Iteration 9
objective 0.000179244685589
residual 0.0189338155473
Starting:  1.0 0.000179244685589
Pass 1: a:1.0; 0.000176409376427 ?<= 0.000179204579423
alpha 1.0
run time 27.5870001316s
Iteration 10
objective 0.000176409376427
residual 0.0187834702026
Starting:  1.0 0.000176409376427
Pass 1: a:1.0; 0.000156729894724 ?<= 0.000176374280248
alpha 1.0
run time 29.4530000687s
Iteration 11
objective 0.000156729894724
residual 0.0177047956624
Starting:  1.0 0.000156729894724
Pass 1: a:1.0; 0.000150571340883 ?<= 0.000156722617084
alpha 1.0
run time 29.7219998837s
Iteration 12
objective 0.000150571340883
residual 0.0173534631059
Starting:  1.0 0.000150571340883
Pass 1: a:1.0; 0.000142571110692 ?<= 0.0001505594308
alpha 1.0
run time 30.1700000763s
Iteration 13
objective 0.000142571110692
residual 0.0168861547246
Starting:  1.0 0.000142571110692
Pass 1: a:1.0; 0.00013600002291 ?<= 0.00014256219471
alpha 1.0
run time 31.8150000572s
Iteration 14
objective 0.00013600002291
residual 0.0164924238916
Starting:  1.0 0.00013600002291
Pass 1: a:1.0; 0.000129851114458 ?<= 0.000135989028627
alpha 1.0
run time 31.7809998989s
Iteration 15
objective 0.000129851114458
residual 0.016115279362
Starting:  1.0 0.000129851114458
Pass 1: a:1.0; 0.000126937895954 ?<= 0.000129843124314
alpha 1.0
run time 32.0469999313s
Iteration 16
objective 0.000126937895954
residual 0.0159334802196
Starting:  1.0 0.000126937895954
Pass 1: a:1.0; 0.000123896444495 ?<= 0.000126933117117
alpha 1.0
run time 28.8510000706s
Iteration 17
objective 0.000123896444495
residual 0.0157414385934
Starting:  1.0 0.000123896444495
Pass 1: a:1.0; 0.000121534905499 ?<= 0.000123893104464
alpha 1.0
run time 28.8229999542s
Iteration 18
objective 0.000121534905499
residual 0.0155906962961
Starting:  1.0 0.000121534905499
Pass 1: a:1.0; 0.000118594260593 ?<= 0.000121530677709
alpha 1.0
run time 28.5299999714s
Iteration 19
objective 0.000118594260593
residual 0.0154009259847
Starting:  1.0 0.000118594260593
Pass 1: a:1.0; 0.000114705432394 ?<= 0.00011858912611
alpha 1.0
run time 28.3970000744s
Iteration 20
objective 0.000114705432394
residual 0.0151463152215
Starting:  1.0 0.000114705432394
Pass 1: a:1.0; 0.000112635657885 ?<= 0.000114697967432
alpha 1.0
run time 29.486000061s
Iteration 21
objective 0.000112635657885
residual 0.0150090411343
Starting:  1.0 0.000112635657885
Pass 1: a:1.0; 0.000109656901616 ?<= 0.000112630428994
alpha 1.0
run time 36.3550000191s
Iteration 22
objective 0.000109656901616
residual 0.0148092472203
Starting:  1.0 0.000109656901616
Pass 1: a:1.0; 0.00010888262822 ?<= 0.000109655741033
alpha 1.0
run time 34.0869998932s
Iteration 23
objective 0.00010888262822
residual 0.014756871499
Starting:  1.0 0.00010888262822
Pass 1: a:1.0; 0.000107837195428 ?<= 0.000108881311633
alpha 1.0
run time 33.507999897s
Iteration 24
objective 0.000107837195428
residual 0.0146858568308
Starting:  1.0 0.000107837195428
Pass 1: a:1.0; 0.000106370083533 ?<= 0.000107834785978
alpha 1.0
run time 33.7380001545s
...run time:  688.184999943s

C:\Anaconda\lib\site-packages\pysit\solvers\model_parameter.py:74: RuntimeWarning: invalid value encountered in sqrt
def unlinearize(cls, data): return np.sqrt(1./data)


First, we can check the convergence:

In [13]:
obj_vals = np.array([v for k,v in invalg.objective_history.items()])
plt.figure()
plt.semilogy(obj_vals)

Out[13]:
[<matplotlib.lines.Line2D at 0x6b1c208>]

Next lets us look at our proposed solution:

In [15]:
C = result.C
clim = C.min(),C.max()

plt.figure()
vis.plot(C0, mesh, clim=clim)
plt.gca().set_aspect('equal')
plt.title('Initial Model')

plt.figure()
vis.plot(result.C, mesh, clim=clim)
plt.gca().set_aspect('equal')
plt.title('Reconstruction')

Out[15]:
<matplotlib.text.Text at 0x4d31bc88>

We can also look at the gradient at the first, a middle, and last step:

In [21]:
plt.figure()
plt.colorbar();

plt.figure()
plt.colorbar();

plt.figure()
plt.colorbar();


Unfortunately, at this resolution, it is difficult to recover the the reflectors exactly (well, it is always difficult...). But we can at least infer that the last reflector and the first reflector are similar in amplitude (in fact the last one is more intense, but this is perhaps muffled by the weaker waves that hit it). The second is definitely the weakest of the three, and in fact it is difficult to distinguish from the ghosts surrounding it.

## True Model¶

In [34]:
C_orig, C0_orig = horizontal_reflector(mesh, reflector_depth=[0.4, 0.6, 0.7], reflector_scaling = [0.6, 0.3, 0.8])
# Depths given as percentage of the domain
speed_factor = (1-0.05*np.pi)
C0_true = speed_factor*C0_orig
dC_refl = C_orig - C0_orig

C = C0_true + dC_refl # + dC_grad
C0 = C0_true

print [d*1.5 for d in [0.4, 0.6, 0.7] ]

plt.figure()
vis.plot(C,mesh)

plt.figure()
C_reshaped = reshape_pysit_vector_to_grid(C, mesh)
plt.plot(np.arange(mesh.z.n)*mesh.z.delta,C_reshaped[mesh.x.n/2, :])

(210L, 150L) (31500, 1)
[0.6000000000000001, 0.8999999999999999, 1.0499999999999998]

Out[34]:
[<matplotlib.lines.Line2D at 0x77e9da0>]

## Details¶

Properties of the physical and computational domain:

• Domain is defined on $[0.0, 2.1) \times [0.0, 1.5)$
• Computational mesh is $210 \times 150$
• $\Delta x = \Delta z = 0.01$ (accessible in the code via the variables mesh.x.delta and mesh.z.delta)

Properties of the data:

• There is one shot of data at $(x,z) = (1.05, 0.05)$
• There are 210 recievers, one at each pixel at depth $z=0.05$

• Data is recorded from $T=0$ until T=4\$

• There are 6824 data points per trace

• $\Delta t \approx 0.00073277$ (available exactly in the variable shot.dt)

• The PySIT solver will automatically interpolate this data as needed

• The reciever data array shot.receivers.data has shape (nt, nr)

Known properties of the medium:

• The medium has a constant background velocity between $0.65$ and $1.1$
• The medium is horizontally layered

• The reciever data can be plotted with pysit.vis.display_seismogram(shot)
• A single trace can be plotted with plt.plot(shot.receivers.data[:,x_idx])
You can plot the resulting migrated image with pysit.vis.plot(g.data, mesh) Note that g is not a numpy array, but g.data is. Additionally, you may want to pass an optional color limit argument clim=(lower, upper) to pysit.vis.plot to help see small details. Good values for upper and lower can be found by glancing at the colorbar obtained with plt.colorbar().
• After running an inversion, the resulting velocity is found in result.C, which can be visualized by plotting with pysit.vis.plot(result.C, mesh)
• cp9.py has an example for plotting the value of the objective function