FWI of open source marine 2D data (Viking Graben line 12)

[kerim371]

Hi,

I'm trying to obtain FWI results from Viking Graben line 12.
The problem is that the solution doesn't converge.

I would appreciate any help.
My setup: Ubuntu 20.04, 8 threads, 20 Gb RAM, Seismic Unix, OpenSeaSeis, Petrel.

link to processed data

Data description (geometry):
water depth about 330 m
number of samples in raw data: 1500 (in FWI I use only 3600 ms to fit in RAM)
sampling interval: 4 ms
each shot contains 119-120 receivers
about 1000 shots overall (in FWI I use records from SP 300 to SP 700)
SRC spacing = 25 m
GRP spacing = 25 m
distance from source to first receiver = 262 m

The processing:

1. gain: tpow=2
2. band pass filtering: f=1,5,70,80, amps=0,1,1,0
3. Vnmo picking and Radon filtering (remove multiples)
4. NMO & stack

The result (wells converted to time domain using VSP data, without doing Well-Tie: thus wells were transformed to time domain unreliably):

Then using prepared Vnmo, picked horizons, Dix formula and smooth filter I compute starting velocity model (Vinterval):

I cut the model up to about 4 km depth with spacings 12.5 m (equal to CDP step) with 8 Hrz signal or 6.25 m with about 16 Hrz signal to fit in my RAM and I use SEGY records up to 3600 ms.
To compute FWI I use the code:

# FWI on the 2D Overthrust model with L-BFGS from the NLopt library
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using Statistics, Random, LinearAlgebra, Interpolations, DelimitedFiles
using JUDI, NLopt, HDF5, SegyIO, PyPlot, Plots

# prestk_file = "/mnt/HDD_2TB/VikingGraben/line12/s_cdp_tpow_radon_sp_4000ms_sp170_1104.sgy"
prestk_file = "/mnt/HDD_2TB/VikingGraben/line12/s_cdp_tpow_radon_sp_3600ms_sp300_700.sgy"
# model_file = "/mnt/HDD_2TB/VikingGraben/line12/vnmo_model_regr_depth.h5"
model_file = "/mnt/HDD_2TB/VikingGraben/line12/vint_model_depth.h5"
wavelet_file = "/mnt/HDD_2TB/VikingGraben/line12/wavelet/FarField_SourceWavelet.dat"
use_wavelet_file = true
shots_out_name = "shot"
shots_out_dir = "/mnt/HDD_2TB/VikingGraben/line12/shots/"

# Bound constraints
vmin = 1.4f0
vmax = 6.5f0

# JUDI options
buffer_size = 1500f0    # limit model (meters)

# Load data and create data vector
block = segy_read(prestk_file)
srcx = get_header(block, :SourceX)
grpx = get_header(block, :GroupX)
min_x = minimum([minimum(srcx), minimum(grpx)])
max_x = maximum([maximum(srcx), maximum(grpx)])
d_obs = judiVector(block)

# Load starting model
n, d, o, m0 = read(h5open(model_file, "r"), "n", "d", "o", "m")
m0x = o[1]:d[1]:o[1]+(size(m0)[1]-1)*d[1]
# ind = [i for i in eachindex(m0x) if min_x-buffer_size <= m0x[i] <= max_x+buffer_size]
ind = [i for i in eachindex(m0x) if min_x <= m0x[i] <= max_x]
model0 = Model((length(ind), n[2]), (d[1], d[2]), (o[1]+(ind[1]-1)*d[1], o[2]), m0[ind,:])

# Set up wavelet and source vector
src_geometry = Geometry(block; key = "source", segy_depth_key = "SourceDepth")
if use_wavelet_file
    wavelet_raw = readdlm(wavelet_file, skipstart=25)
    itp = linear_interpolation(1:length(wavelet_raw), wavelet_raw[:,1], extrapolation_bc=0f0)
    wavelet = Matrix{Float32}(undef,src_geometry.nt[1],1)
    wavelet[:,1] = itp(0:src_geometry.dt[1]:src_geometry.t[1])
else
    f0 = 0.016     # kHz
    wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], f0)
end
q = judiVector(src_geometry, wavelet)

############################################## FWI #################################################

# JUDI options
global jopt = JUDI.Options(
    save_data_to_disk=true, 
    file_path=shots_out_dir, 
    file_name=shots_out_name,
    limit_m = true,
    buffer_size = buffer_size,
    isic = true)

# optimization parameters
fevals = 10    # allow 10 function evalutations
batchsize = 40
count = 0
fhistory = zeros(Float32, fevals+1)

# NLopt objective function
println("No.  ", "fval         ", "norm(gradient)")
function f!(x, grad)

    # Update model
    model0.m .= x

    global jopt;
    # Select batch and calculate gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, gradient = fwi_objective(model0, q[i], d_obs[i], options=jopt)

    # Reset gradient in water column to zero
    gradient = reshape(gradient, model0.n);
    # gradient[:, 1:21] .= 0f0
    grad[1:end] .= gradient[1:end]

    global count; count += 1
    println("iteration: ", count, "\tfval: ", fval, "\tnorm: ", norm(grad))
    fhistory[count] = fval
    return convert(Float64, fval)
end

# Optimization parameters
opt = Opt(:LD_LBFGS, prod(model0.n))
min_objective!(opt, f!)
mmin = (1f0 ./ vmax)^2   # bound constraints on slowness squared
mmax = (1f0 ./ vmin)^2
lower_bounds!(opt, mmin); upper_bounds!(opt, mmax)
maxeval!(opt, fevals)
(minf, minx, ret) = optimize(opt, copy(model0.m))

# Save results, function values and elapsed time
h5open("result_2D_overthrust_lbfgs.h5", "w") do file
    write(file, "x", sqrt.(1f0 ./ reshape(minx, model0.n)), "fhistory", fhistory)
end

conv_plt = Plots.plot(fhistory/norm(fhistory, Inf), 
    xlabel="Iteration no.",
    ylabel="Normalized residual",
    title="Convergence of FWI w/ L-BFGS")
Plots.savefig(conv_plt,"Convergence.png")

x = o[1]/1000f0:d[1]/1000f0:o[1]/1000f0+(n[1]-1)*d[1]/1000f0
z = o[2]/1000f0:d[2]/1000f0:o[2]/1000f0+(n[2]-1)*d[2]/1000f0
xx = x[ind]
fwi_plt = Plots.heatmap(xx, z, sqrt.(1f0 ./ adjoint(model0.m)).data, c=:rainbow, 
    xlims=(xx[1],xx[end]), 
    ylims=(z[1],z[end]), yflip=true,
    title="Convergence of FWI w/ L-BFGS",
    xlabel="Lateral position [km]",
    ylabel="Depth [km]",
    dpi=600)
Plots.savefig(fwi_plt,"FWI.png")

[mloubout]

Few comments on top of my head.

• The data seems to be preprocessed for imaging rather than FWI as the turning waves are muted (water mute)
• isic is made for imaging so you do not want to use this for FWI
• You probably need a model topmute (mute the water layer in the gradient) because it usually contained such strong amplitudes at the source location it will prevent any convergence

Did you look at the first gradient to see roughly what it's trying to update?

[kerim371]

• The data seems to be preprocessed for imaging rather than FWI as the turning waves are muted (water mute)

Oh, you mean the straight wave that penetrates in the water? If so then yes it was muted by NMO stretch factor (reverse NMO was applied). Completely forgot about that :)

Did you look at the first gradient to see roughly what it's trying to update?

Good idea, haven't done this.

Tomorrow I will try to fix that and launch it again.

How do you think, should I use Ricker wavelet or better Bolt gun signature recorded at farfield (as in the picture below)? And what frequency should I try?

[kerim371]

@mloubout It seems I do something incorrect. If you have ideas...

The gradient at any iteration looks like in picture below:

Calculating single shot at position 5712m with 119 receivers [2500, 5450] gives the following image with the same velocity model and same 16 Hrz Ricker source wavelet:

Also at each iteration I set gradient in water to 0:

    water_depth = Int(350)  # meters
    gradient[:, 1:Int(water_depth/d[2])] .= 0f0

Model grid spacings 6.25

[mloubout]

Looks like a bad reshape that's odd.

[kerim371]

@mloubout thank you for the hint, I forgot to transpose it.
Fist gradient looks like this:

[kerim371]

Also I've noticed that I forgot to multiply depths in my model by 0.5 (two-way travel time staff). The initial model is incorrect. So I do the calculations once again.

[mloubout]

First grad looks as expected that's a good start. Now of you wanna do FWI this is too high frequency so you'll have to work with low-to-high frequency band where for each band you run like 10 iteration with low pass data and source.

Do you still only have the nmo data or do you have the original data with turning waves

[zhangxiaoshuotttt]

@mloubout It seems I do something incorrect. If you have ideas...

The gradient at any iteration looks like in picture below:

Calculating single shot at position 5712m with 119 receivers [2500, 5450] gives the following image with the same velocity model and same 16 Hrz Ricker source wavelet:

Also at each iteration I set gradient in water to 0:

    water_depth = Int(350)  # meters
    gradient[:, 1:Int(water_depth/d[2])] .= 0f0

Model grid spacings 6.25

I also tried to calculate FWI, but my gradient field is all 0, can you give me some hints to help me know what is wrong . Here is my issue: #176