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

[zhangxiaoshuotttt]

Discussed in #175

^{Originally posted by zhangxiaoshuotttt March 17, 2023}
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:

[kerim371]

@zhangxiaoshuotttt hi,

First of all are you sure that your gradient is all zero?
I would try to check that by clipping displayable range of values coded by color (vmin/vmax):

imshow(reshape(g, model0.n)',vmin=-maximum(g)/10f0,vmax=maximum(g)/10f0,  extent=(0,10,4,0), cmap="jet")

[zhangxiaoshuotttt]

@zhangxiaoshuotttt hi,

First of all are you sure that your gradient is all zero? I would try to check that by clipping displayable range of values coded by color (vmin/vmax):

imshow(reshape(g, model0.n)',vmin=-maximum(g)/10f0,vmax=maximum(g)/10f0,  extent=(0,10,4,0), cmap="jet")

Yes,here is the valure of g that I outputed and the image

[zhangxiaoshuotttt]

Is it possible that there is a problem with my coordinate setting? In the seismic data,I set the y-coordinate of the source point, the receiver point and the CMP point to 0.

The X-coordinate of my source,receiver and CMP points do not coincide. And I set model geometry by:

min_source_x = minimum(get_header(block,:SourceX))

shape = (942,361) # Number of gridpoints nx, nz
spacing = (12.5,12.5) # #n meters here
origin = (min_source_x, 0) # In meters as well

d_model = judiVector(VpModel)/1000;
d_model_T = d_model.data[1]'
m0 = 1f0./d_model_T.^2f0
model0 = Model(shape,spacing,origin,m0)

[kerim371]

From the first look I don't see any mistakes.
To find the cause of that I recommend to do the following:

1. define somewhere above a function save_results(model, gradient, count). In that function save plotting result to a .png image (the same way you call imshow(reshape(g, model0.n)',vmin=-maximum(g)/10f0,vmax=maximum(g)/10f0, extent=(0,10,4,0), cmap="jet") but save it to the .png).
2. within function objf!(x, grad) call the function save_results(model0.m, grad, count). Thus you will see what grad you have after each iteration.

Also I would pay attention if grad gets updated within this block:

    if length(grad) > 0
        grad[1:end] = vec(gradient)
    end

And also what returns fval, gradient = fwi_objective(model0, q[i], d_obs[i],options=jopt)

Saving these results to file will get you understanding what is wrong.

And for testing purposes you don't need to run all the shots, probably 1 shot at each iteration would be enough for debugging purposes.

[zhangxiaoshuotttt]

From the first look I don't see any mistakes. To find the cause of that I recommend to do the following:

1. define somewhere above a function save_results(model, gradient, count). In that function save plotting result to a .png image (the same way you call imshow(reshape(g, model0.n)',vmin=-maximum(g)/10f0,vmax=maximum(g)/10f0, extent=(0,10,4,0), cmap="jet") but save it to the .png).
2. within function objf!(x, grad) call the function save_results(model0.m, grad, count). Thus you will see what grad you have after each iteration.

Also I would pay attention if grad gets updated within this block:

    if length(grad) > 0
        grad[1:end] = vec(gradient)
    end

And also what returns fval, gradient = fwi_objective(model0, q[i], d_obs[i],options=jopt)

Saving these results to file will get you understanding what is wrong.

And for testing purposes you don't need to run all the shots, probably 1 shot at each iteration would be enough for debugging purposes.

Thank you for your advice Sir. I will check it as you said

[kerim371]

@zhangxiaoshuotttt yes the problem may be related to the coordinates.
Why you have Source/Rec X position about millions and your model is 50 km?
Also your gradient picture is only 10 km. If your testing shots out of these 10 KM then you won't see gradient changes. Probably leave extent without modification.

You have to check coordinate scaler trace header as well. To be sure that you coordinates are fine please do the forward modeling for a single shot with the same model. If your coordinates are completely wrong then forward modeling will return empty. I see that RecSourceScalar=-10 that means that your coordinates will be mulltiplied by -10 or by 0.1 I don't remember exactly how JUDI handle this.

More over JUDI options for FWI should have IC = "fwi".

[mloubout]

I think there is a scaling issue accessing th le min x through the headers as it should be rescaled by 1e4.

[mloubout]

You can double check by making sure the origin coincides with Geometry(q[1].geometry).xloc[1]

[zhangxiaoshuotttt]

You can double check by making sure the origin coincides with Geometry(q[1].geometry).xloc[1]
Ok. I'll check. Thanks a lot.

[zhangxiaoshuotttt]

@zhangxiaoshuotttt yes the problem may be related to the coordinates. Why you have Source/Rec X position about millions and your model is 50 km? Also your gradient picture is only 10 km. If your testing shots out of these 10 KM then you won't see gradient changes. Probably leave extent without modification.

You have to check coordinate scaler trace header as well. To be sure that you coordinates are fine please do the forward modeling for a single shot with the same model. If your coordinates are completely wrong then forward modeling will return empty. I see that RecSourceScalar=-10 that means that your coordinates will be mulltiplied by -10 or by 0.1 I don't remember exactly how JUDI handle this.

More over JUDI options for FWI should have IC = "fwi".

I rezeroed the origin by uniformly subtracting the x origin coordinate. Next, I will make a forward model to verify whether there are still problems with my coordinates. Thank you for your advice.

[zhangxiaoshuotttt]

You can double check by making sure the origin coincides with Geometry(q[1].geometry).xloc[1]

My "Function val" is not zero , it is "3.47902e+07".

[zhangxiaoshuotttt]

@kerim371 @mloubout
Thanks for your help, I solved the problem, and here's my gradient field.

But I still have a question, why does my gradient field looks so shallow? The depth of my OBN nodes is floating around a depth of 90m.

By the way, I used statistical wavelet.

[mloubout]

I think your wavelet has a time shift so there is very little illumination. Check what JUDI default tucker wavelet looks like you'll see it need to start around time=0

[zhangxiaoshuotttt]

Thank you very much, this advice is really helpful!