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();
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 ####.
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)
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:
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:
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...
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();
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();
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.
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();
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
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,
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.
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.
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
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.
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)
First, we can check the convergence:
obj_vals = np.array([v for k,v in invalg.objective_history.items()])
plt.figure()
plt.semilogy(obj_vals)
Next lets us look at our proposed solution:
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')
We can also look at the gradient at the first, a middle, and last step:
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.
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, :])
Properties of the physical and computational domain:
mesh.x.delta and mesh.z.delta)Properties of the data:
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:
Helpful PySIT routines:
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])
The migration operator can be applied to the residual with the code (where C0 is the model about which you wish to do migration m0 = solver.ModelParameters(mesh,{'C': C0}) g = objective.compute_gradient(shots, m0)
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)
To view a single depth slice of the result, you can reshape the resulting velocity vector and plot it C_reshape = reshape_pysit_vector_to_grid(result.C, mesh) plt.plot(C_reshape[105, :]) # Reshaped vector is indexed (x_idx,z_idx)
cp9.py has an example for plotting the value of the objective function