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 scipy.io import loadmat

from pysit import *
from pysit.gallery import horizontal_reflector

from cp9_sol import load_challenge_data

#	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()

shots = load_challenge_data(mesh)					 
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})
g = objective.compute_gradient(shots, m0)

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})
g = objective.compute_gradient(shots, m0)

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})
g = objective.compute_gradient(shots, m0)

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})
g = objective.compute_gradient(shots, m0)

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})
g = objective.compute_gradient(shots, m0)

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,150], label="Trace 150")
plt.plot(shot.receivers.data[:2500,155], label="Trace 155") # ~middle index
plt.legend();
In [8]:
x = (155-150)*mesh.x.delta
max_150 = np.where(shot.receivers.data[:,150] == shot.receivers.data[:,150].max())[0][0]
max_155 = np.where(shot.receivers.data[:,155] == shot.receivers.data[:,155].max())[0][0]
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})
g = objective.compute_gradient(shots, m0)

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||
                        'gradient_frequency'		: 1,  # g_k
                        'gradient_length_frequency' : 1,  # ||g_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
  gradnorm 5.07702344148
  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
  gradnorm 4.33496675138
  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
  gradnorm 0.967884193673
  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
  gradnorm 0.559995336607
  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
  gradnorm 0.472170679645
  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
  gradnorm 0.377573437476
  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
  gradnorm 0.371795877604
  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
  gradnorm 0.606441236078
  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
  gradnorm 0.442643193555
  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
  gradnorm 0.434964752642
  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
  gradnorm 1.1968640548
  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
  gradnorm 0.355419930403
  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
  gradnorm 0.251523410369
  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
  gradnorm 0.21781550684
  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
  gradnorm 0.178373908195
  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
  gradnorm 0.187583315682
  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
  gradnorm 0.298839497386
  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
  gradnorm 0.16455734996
  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
  gradnorm 0.141174581672
  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
  gradnorm 0.158852835356
  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
  gradnorm 0.148615411008
  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
  gradnorm 0.276291964574
  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
  gradnorm 0.0942254819576
  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
  gradnorm 0.068088634427
  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
  gradnorm 0.0725120679971
  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()
vis.plot(invalg.gradient_history[0].data, mesh, clim=(-0.01,0.01))
plt.colorbar();

plt.figure()
vis.plot(invalg.gradient_history[12].data, mesh, clim=(-0.01,0.01))
plt.colorbar();

plt.figure()
vis.plot(invalg.gradient_history[24].data, mesh, clim=(-0.001,0.001))
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:

Properties of the data:

Known properties of the medium:

Helpful PySIT routines: