Exercise 7 : Seismic imaging
In this exercise, we conduct a small 2D seismic imaging experiment on a 4-layer model example with JUDI. We suggest you use a Docker image to run the experiment so that software packages are properly installed. To use the docker image, first install docker. Then, in the terminal, do
docker run -p 8888:8888 ddjj1118/judi_eas_project:v4.0
Running this command will produce an output that looks like
Copy/paste this URL into your browser when you connect for the first time,
to login with a token:
http://af637030c092:8888/?token=8f6c664eb945f9c6b7cd72669fef04a6dc70c08194cb87e9
or http://127.0.0.1:8888/?token=8f6c664eb945f9c6b7cd72669fef04a6dc70c08194cb87e9
Copy paste the URL in your browser and replace (af637030c092 or 127.0.0.1) by localhost.
You will then be directed to a jupyter folder that contains the notebooks for the projects.
First, make sure we have the latest version of Devito installed. You can do this by click NEW, then go for TERMINAL in the browser with jupyter folder, then do
pip install devito==4.3
Then, we can create a new julia script and run the experiment in it.
Set up the experiment
## First do using Pkg; Pkg.add(xxx); Pkg.update(xxx) to make sure they are installed and in the latest version
using JUDI, PyPlot, Images, JOLI, IterativeSolvers, LinearAlgebra, Printf, Statistics
Here, we set up a 4-layer model, where the velocities in each layer are 1.5/2.5/3/3.5 km/s. We often refer the first layer as the water layer/column in marine acquisition.
# number of gridpoints
n = (201, 101)
# Grid spacing
d = (10.0, 10.0)
# Origin
o = (0., 0.)
# Velocity [km/s]
v = ones(Float32,n) .+ 0.5f0
v[:,20:50] .= 2.5f0
v[:,51:71] .= 3f0
v[:,71:end] .= 3.5f0
Wave-equation simulations in JUDI are based on squared slowness m. Also, we make up a background model m0 which is a smoothed version of the ground truth one. To make things simple, we assume to know the exact depth of ocean bottom and keep the water column fixed in the background model. The dm, as the difference of m and m0, only contains the sharp reflectors.
# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = convert(Array{Float32,2},imfilter(m, Kernel.gaussian(2)))
m0[:,1:19] = m[:,1:19]
dm = vec(m - m0)
figure();
subplot(1,3,1)
imshow(reshape(m,n)');title("m")
subplot(1,3,2)
imshow(reshape(m0,n)');title("m0")
subplot(1,3,3)
imshow(reshape(dm,n)',cmap="Greys");title("dm")
m and m0 are stored in model structure in JUDI, with physical information.
model0 = Model(n, d, o, m0; nb = 80)
model = Model(n, d, o, m; nb = 80)
We set up the source geometry and receiver geometry as below. Simulations are in 2D so y-direction is always 0. The sources are assumed to fire one-by-one while receivers are fixed.
nsrc = 8 # num of sources
xsrc = convertToCell(range(900f0, stop=1100f0, length=nsrc)) # in [m]
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc)) # in [m]
zsrc = convertToCell(range(6f0, stop=6f0, length=nsrc)) # in [m]
timeS = 1000f0 # recording time [ms]
dtS = 1f0 # time sampling [ms]
srcGeometry = Geometry(xsrc,ysrc,zsrc; dt=dtS, t=timeS)
nrec = 100 # num of receivers
xrec = range(d[1], stop=(n[1]-1)*d[1], length=nrec) # in [m]
yrec = 0f0 # in [m]
zrec = range(10f0, stop=10f0, length=nrec) # in [m]
timeR = 1000f0 # recording time [ms]
dtR = 1f0 # time sampling [ms]
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR, nsrc=nsrc)
Set up the source with Ricker wavelet
f0 = 0.02f0 # 20 Hz wavelet
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)
Set up computational time step
# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)
Set up forward modeling operators in ground truth velocity, background velocity, and migration operator
F = judiModeling(info, model, srcGeometry, recGeometry) # forward modeling w/ true model
F0 = judiModeling(info, model0, srcGeometry, recGeometry) # forward modeling w/ background model
J = judiJacobian(F0,q) # demigration operator (adjoint of J is migration)
Generate data
Generate data in ground truth velocity and background velocity
dobs = F*q
dobs0 = F0*q
Task: plot a single shot record of dobs and dobs0. Choose the same source. What do you see in the shot record? If there are a couple of events, try to match them up with the corresponding reflectors.
# hint: seismic data is as a judiVector
# you can do this by figure();imshow(dobs.data[1],vmin=-0.02*norm(dobs.data[1],Inf),cmap="seismic",vmax=0.02*norm(dobs.data[1],Inf),aspect="auto")
# This plots the shot record in correct velocity generated by the 1st source
# IMPORTANT: always adjust vmin, vmax, aspect so that we can see the events clearly
You can also generate a linearized shot record by
dlin = J*dm
Task: do the same comparison of dlin and dobs-dobs0. What do you see? Why do you observe this theoretically (think about math)?
Reverse time migration (RTM)
From now on, we will only focusing on migrating the linearized data for simplicity. First, we can try reverse time migration by
rtm1 = J'*dlin
To make the RTM image cleaner and with higher quality, Witte et al developed an inverse-scattering imaging condition (ISIC), which has been incorporated in JUDI.
J.options.isic = true
rtm2 = J'*dlin
Task: compare RTM results w/ and w/o ISIC, what do you see? Remember again to clip the vmin/vmax to see all the events in the images.
# hint: rtm1/rtm2 are PhysicalParameters. Do rtm1.data, rtm2.data to access the value. Again take care of vmin, vmax, aspect of plotting.
Least-squares reverse time migration (LS-RTM)
Seismic imaging researchers are not always satisfied with RTM. Seismic imaging basically solves the optimization problem
\min_{\mathbf{\delta m}} \frac{1}{2}\sum_{i=1}^{n_s}\|\mathcal{J}_i\mathbf{\delta m} - \mathbf{\delta d}\|_2^2
where RTM is only taking a full gradient of the objective w.r.t. \delta m. To get better images, we can minimize the objective by LSQR and apply a right preconditioner
# Right Preconditioner
Tm = judiTopmute(model0.n, 19, 2) # Mute water column
S = judiDepthScaling(model0)
Mr = S*Tm
## vanilla LSQR
x1 = 0f0 .* model0.m
lsqr!(x1,J*Mr,dlin;maxiter=2,atol=0f0,verbose=true)
The final solution is given by Mr*x1 (Mr is a right preconditioner so we are actually minimizing \frac{1}{2}\|\mathcal{J}\mathbf{M_r}\mathbf{x}- \mathbf{\delta d}\|_2^2)
Task: compare LS-RTM result with the previous RTM result, what do you see?
Sparsity-promoting least-squares reverse time migration (SPLS-RTM)
The start-of-the-art imaging technique is to do LS-RTM while promoting sparsity of solution in Curvelet domain, see Witte et al. This can be achieved by linearized Bregman iterations, as shown below
# Soft thresholding functions and Curvelet transform
soft_thresholding(x::Array{Float64}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float64, lambda), 0.0)
soft_thresholding(x::Array{Float32}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float32, lambda), 0f0)
n = model0.n
C0 = joCurvelet2D(n[1], 2*n[2]; zero_finest = false, DDT = Float32, RDT = Float64)
function C_fwd(im, C, n)
im = hcat(reshape(im, n), reshape(im, n)[:, end:-1:1])
coeffs = C*vec(im)
return coeffs
end
function C_adj(coeffs, C, n)
im = reshape(C'*coeffs, n[1], 2*n[2])
return vec(im[:, 1:n[2]] .+ im[:, end:-1:n[2]+1])
end
C = joLinearFunctionFwd_T(size(C0, 1), n[1]*n[2],
x -> C_fwd(x, C0, n),
b -> C_adj(b, C0, n),
Float32,Float64, name="Cmirrorext")
src_list = Set(collect(1:nsrc))
batchsize = 2
lambda = 0f0
x2 = 0f0 .* model0.m
z = deepcopy(x2)
niter = 8
# Main loop
for j = 1:niter
# Select batch and set up left-hand preconditioner
length(src_list) < batchsize && (global src_list = Set(collect(1:nsrc)))
i = [pop!(src_list) for b=1:batchsize]
println("LS-RTM Iteration: $(j), imaging sources $(i)")
flush(Base.stdout)
residual = J[i]*Mr*x2-dlin[i]
phi = 0.5 * norm(residual)^2
g = Mr'*J[i]'*residual
# Step size and update variable
t = Float32.(2*phi/norm(g)^2)
# Update variables and save snapshot
global z -= t*g
C_z = C*z
(j==1) && (global lambda = quantile(abs.(C_z), .95)) # estimate thresholding parameter in 1st iteration
global x2 = adjoint(C)*soft_thresholding(C_z, lambda)
@printf("At iteration %d function value is %2.2e and step length is %2.2e \n", j, phi, t)
@printf("Lambda is %2.2e \n", lambda)
end
The solution is given by Mr*x2.