Support for Joint Surface and Borehole Acquisition in LS-RTM (2D Marmousi Example)

[WYJLCYWHZ]

Hi, thank you for the great work on JUDI! I'm currently experimenting with the 2D Marmousi LS-RTM example (lsrtm_marmousi_sgd.jl) and I would like to extend it to support joint surface and borehole (VSP) acquisition.

My idea is to concatenate both surface and borehole datasets into a single judiVector, and create a unified modeling operator (judiModeling) and Jacobian (judiJacobian) for both source-receiver geometries. Then, I would randomly select sources from the joint dataset during each SGD iteration.

Before I proceed, I have a few questions:


1. Does JUDI natively support this kind of joint surface + borehole modeling and inversion setup?
   
2. Are there any existing examples or best practices for this use case using JUDI?

[mloubout]

JUDI does not make any "surface" assumptions you can put the source and receiver location anywhere you want.

[WYJLCYWHZ]

JUDI does not make any "surface" assumptions you can put the source and receiver location anywhere you want.

Thank you. Are there any existing examples for joint surface and borehole (VSP) RTM example using JUDI?

[mloubout]

There is no specific example no. Is there any issue you ran into trying it?

[WYJLCYWHZ]

There is no specific example no. Is there any issue you ran into trying it?

Hi! My initial approach was to concatenate both surface and borehole datasets into a single judiVector, then create a unified modeling operator (judiModeling) and Jacobian (judiJacobian) for both source-receiver geometries. During each SGD iteration, randomly selected sources from the joint dataset. However, this approach didn't perform as well as I expected. It performed far worse than using only the surface data.

Additionally, the current implementation does not effectively incorporate the correlation between the borehole and surface data in the inversion process. I believe that there should be a significant correlation between the wavefields of the borehole and the surface. I would like to ask for your advice on how to better address this issue. Thank you.

The code：

using Statistics, Random, LinearAlgebra
using JUDI, SegyIO, HDF5, PyPlot

#################################################################################################

n, d, o, m0 = read(h5open("marmousi_model.h5", "r"), "n", "d", "o", "m0")
model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)

#################################################################################################
block_surf = segy_read("datacube_in.segy")
geom_surf = Geometry(block_surf; key = "source", segy_depth_key = "SourceDepth")
wavelet_surf = ricker_wavelet(geom_surf.t[1], geom_surf.dt[1], 0.03)
d_surf = judiVector(block_surf)
q_surf = judiVector(geom_surf, wavelet_surf)

block_vsp = segy_read("marmousi_2D1_vsp.segy")
geom_vsp = Geometry(block_vsp; key = "source", segy_depth_key = "SourceDepth")
wavelet_vsp = ricker_wavelet(geom_vsp.t[1], geom_vsp.dt[1], 0.03)
d_vsp = judiVector(block_vsp)
q_vsp = judiVector(geom_vsp, wavelet_vsp)

#################################################################################################

nsrc = q_surf.nsrc

d_joint = judiVector()
q_joint = judiVector()

for i in 1:nsrc
   
    geom_joint = Geometry([geom_surf[i]; geom_vsp[i]])
    wavelet_joint = wavelet_surf  # or wavelet_vsp, assumed same
    

    data_joint = vcat(d_surf[i], d_vsp[i])
    

    push!(q_joint, judiVector(geom_joint, wavelet_joint))
    push!(d_joint, data_joint)
end

#################################################################################################

opt = Options(optimal_checkpointing = true)
M = judiModeling(model0, q_joint.geometry, d_joint.geometry; options = opt)
J = judiJacobian(M, q_joint)

Mr = judiTopmute(model0; taperwidth = 10)
Ml = judiDataMute(q_joint.geometry, d_joint.geometry)

#################################################################################################

x = zeros(Float32, prod(model0.n))
niter = 20
batchsize = 1
fval = zeros(Float32, niter)

#################################################################################################

for j = 1:niter
    println("Iteration: ", j)
    

    i = randperm(nsrc)[1:batchsize]


    r = Ml[i]*J[i]*Mr*x - Ml[i]*d_joint[i]
    println("Residual size: ", size(r))

    g = adjoint(Mr) * adjoint(J[i]) * adjoint(Ml[i]) * r########chengxiang
    g = vec(g)
    println("Gradient size: ", size(g))

    fval[j] = 0.5f0 * norm(r)^2
    t = norm(r)^2 / norm(g)^2
    global x -= t * g#diedaiyouhua
    println("Model updated, size: ", size(x))
end

[mloubout]

then create a unified modeling operator (judiModeling)

Yes that sounds like a good approach.

However, this approach didn't perform as well as I expected. It performed far worse than using only the surface data.

Would you have pictures/... to support that?

Ml = judiDataMute(q_joint.geometry, d_joint.geometry)

I honestly have no idea how that will behave with VSP data, I would assume very poorly. I would plot a muted shot to check what it does.

The issue with VSP is that you will have a strong imprint of the vsp in the gradient, it might mess up the inversion. I would look at the first gradient to have an idea