I tried FWI with JUDI on the actual marine data, but my gradient was always 0

[zhangxiaoshuotttt]

Hi,

I tried to obtain FWI results from marine data, but I can't get gradient, it is always 0.
I would appreciate some help or tips.

My setup: Ubuntu22.04;12 threads; 16GB RAM; Seismic Unix

Data Description
Water depth about 100 m
number of samples in raw data: 1950(I used only 3900ms; about 4.5km)
sample intervals: 2 ms
each shot contains 230 receivers
In FWI I use 300 shots.
src spacing: 37.5m

The processing:

1. I removed the surface waves and abnormal amplitudes.
2. band pass filtering: f=1,5,70,80

To compute FWI I used the code:

using SegyIO, HDF5, PyPlot, JUDI, NLopt, Random, LinearAlgebra, Printf,Distributed

##Setting  up the initial model geometry
shape = (942, 361) # Number of gridpoints nx, nz
spacing = (12.5, 12.5) # #n meters here
origin = (0.0, 5001) # In meters as well

##Load the model
VpModel = segy_read("Z_Velocity_Field_Depth.sgy");
d_model = judiVector(VpModel)/1000;  #its minvalue:1.5 ; maxvalue: 3.5
d_model_T = d_model.data[1]'
m0 = 1f0./d_model_T.^2f0
model0 = Model(shape,spacing,origin,m0)

##Plot
extent = [0, 50, 4.5, 0]
figure(figsize=(7, 7))
title("Starting OBN model")
imshow(sqrt.(1f0./m0)', cmap="jet", extent=extent, aspect=4)
xlabel("Lateral position");
ylabel("Depth [km]");

##load marine data
block = segy_read("Z.sgy");
d_obs = judiVector(block);

##Plot
extent = [0, 10, 3.9, 0]
figure(figsize=(14, 14))
subplot(221)
imshow(d_obs.data[1], vmin=-0.1, vmax=0.1, cmap="PuOr", extent=extent, aspect=4, interpolation="hamming")
xlabel("Receiver position(km)")
ylabel("Time(s)")
subplot(222)
imshow(d_obs.data[75], vmin=-0.1, vmax=0.1, cmap="PuOr", extent=extent, aspect=4, interpolation="hamming")
xlabel("Receiver position(km)")
ylabel("Time(s)")
subplot(223)
imshow(d_obs.data[150], vmin=-0.1, vmax=0.1, cmap="PuOr", extent=extent, aspect=4, interpolation="hamming")
xlabel("Receiver position(km)")
ylabel("Time(s)")
subplot(224)
imshow(d_obs.data[225], vmin=-0.1, vmax=0.1, cmap="PuOr", extent=extent, aspect=4, interpolation="hamming")
xlabel("Receiver position(km)")
ylabel("Time(s)")
tight_layout()

##Set up wavelet
src_geometry = Geometry(block; key="source",segy_depth_key="SourceDepth");
src_data = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.008f0);
q = judiVector(src_geometry, src_data);

batchsize = 30;
count = 0;

# NLopt objective function
function objf!(x, grad)
    if count == 0
        @printf("%10s %15s %15s\n","Iteration","Function Val","norm(g)")
    end
    # Update model
    model0.m .= Float32.(reshape(x, model0.n))
    #JUDI Options
    jopt = JUDI.Options(limit_model_to_receiver_area=true,IC='as')
    # Seclect 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:8] .= 0f0
    if length(grad) > 0
        grad[1:end] = vec(gradient)
    end
    global count += 1
    @printf("%10d %15.5e %15.5e\n",count, fval, norm(g))
    return convert(Float64, fval)
end

g = zeros(prod(model0.n))
f0 = objf!(vec(model0.m), g)
# Reset count
global count = 0;

imshow(reshape(g, model0.n)',vmin=-1e3,vmax=1e3,  extent=(0,10,4,0), cmap="jet")
title("FWI first gradient")
xlabel("Lateral position [km]");
ylabel("Depth [km]");

And here is my gradient:

[mloubout]

Could you check your source positions? The origin might not be zero putting all your sources outside the domain.

Is fval zero too?

[zhangxiaoshuotttt]

Yes Sir, I also noticed the problem. I reset the geometry here, but it still doesn't work #176