-
Notifications
You must be signed in to change notification settings - Fork 29
/
fwi_example_2D.jl
69 lines (54 loc) · 2.08 KB
/
fwi_example_2D.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
# 2D FWI on Overthrust model using minConf library
# Author: Philipp Witte, [email protected]
# Date: December 2017
#
using Statistics, Random, LinearAlgebra
using JUDI, SlimOptim, HDF5, SegyIO, PyPlot
# Load starting model
n,d,o,m0 = read(h5open("$(JUDI.JUDI_DATA)/overthrust_model.h5","r"), "n", "d", "o", "m0")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)
# Bound constraints
v0 = sqrt.(1f0 ./ m0)
vmin = ones(Float32, model0.n) .* 1.3f0
vmax = ones(Float32, model0.n) .* 6.5f0
vmin[:,1:21] .= v0[:,1:21] # keep water column fixed
vmax[:,1:21] .= v0[:,1:21]
# Slowness squared [s^2/km^2]
mmin = vec((1f0 ./ vmax).^2)
mmax = vec((1f0 ./ vmin).^2)
# Load data
block = segy_read("$(JUDI.JUDI_DATA)/overthrust_shot_records.segy")
d_obs = judiVector(block)
# Set up wavelet
src_geometry = Geometry(block; key="source")
wavelet = ricker_wavelet(src_geometry.t[1],src_geometry.dt[1],0.008f0) # 8 Hz wavelet
q = judiVector(src_geometry,wavelet)
############################### FWI ###########################################
F0 = judiModeling(deepcopy(model0), src_geometry, d_obs.geometry)
# Optimization parameters
niterations = parse(Int, get(ENV, "NITER", "10"))
batchsize = 16
fhistory_SGD = zeros(Float32, niterations)
# Projection operator for bound constraints
proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)]; dims=2),model0.n)
ls = BackTracking(order=3, iterations=10, )
# Main loop
for j=1:niterations
# get fwi objective function value and gradient
i = randperm(d_obs.nsrc)[1:batchsize]
fval, gradient = fwi_objective(model0, q[i], d_obs[i])
p = -gradient/norm(gradient, Inf)
println("FWI iteration no: ",j,"; function value: ",fval)
fhistory_SGD[j] = fval
# linesearch
function ϕ(α)
F0.model.m .= proj(model0.m .+ α * p)
misfit = .5*norm(F0[i]*q[i] - d_obs[i])^2
@show α, misfit
return misfit
end
step, fval = ls(ϕ, 1f-1, fval, dot(gradient, p))
# Update model and bound projection
model0.m .= proj(model0.m .+ step .* p)
end
figure(); imshow(sqrt.(1f0./adjoint(model0.m))); title("FWI with SGD")