How to generate reverse time migration

[1106-nie]

Sorry to interrupt, but I was wondering if we could get a RTM if we changed“n_epochs” to 1

[1106-nie]

[ar4]

Hello and thank you for your question.

Yes, if you perform only one update step then it is equivalent to RTM. It depends on how you wish to define RTM and what you wish the physical meaning of your image amplitudes to be, but if the gradient of the FWI loss with respect to the velocity model is acceptable, then you can get slightly better computational performance by using the regular scalar propagator instead of the Born propagator, and extracting the result from the ".grad" attribute of the velocity model.

For example:

v.requires_grad_()
out = deepwave.scalar(v, ...)
loss = loss_fn(out[-1], observed_data)  # You may wish to mask direct arrivals here
loss.backward()
image = v.grad.detach()

[1106-nie]

Dear professer:
I'm sorry to bother you again.The first image is the RTM I made and the second image is theLSRTM. But I feel that this LSRTM is not very clear and I want to know if this image result is correct or not.

[ar4]

Hello again,

I'm glad to see that you got both RTM and LSRTM to work successfully.

LSRTM does seem to have made some improvements to the RTM image: the deep part and sides of the image have better amplitude, and some artifacts have been reduced. The amplitude of the shallow layers has been reduced relative to some of the deeper layers, but that is probably correct for this model as I think the velocity difference between the shallow layers is small. There are some artifacts remaining, though, which I suspect you are referring to when you say that the LSRTM image is not clear.

Many of the remaining artifacts are probably caused by multiples. RTM/LSRTM assume that the dataset only contains primary reflections, so multiples will cause artifacts. Trying to remove multiples, especially internal multiples, is unfortunately difficult. If you would like to cheat a little bit, you could try slightly smoothing the true model before you use it to generate the observed data, which might help.

May I ask what you used as your migration velocity model? Also, how did you remove the direct arrival? Did you run forward modeling with the migration velocity model (with the max_vel parameter set appropriately to ensure the same time step size was used as when generating the observed data) and subtract the result from the observed data to produce the target data used in RTM/LSRTM? (If you do not understand some of the things that I said, please let me know.) Which optimiser did you use, and how many iterations did you run it for?

[1106-nie]

Dear professer:
I am glad that you have given your valuable comments. Regarding the migration model, I used a Gaussian function to smooth the velocity model to produce the migration model. I am not sure if my perception of the migration model is correct. In order to remove the direct waves, I modeled the original and smoothed velocities orthogonally with scalars in the same observing system, and then subtracted the original velocity seismic data from the smoothed velocity seismic data. For the number of iterations of LSRTM, I chose 10 times, and in practice, it does not work well when the number is 3 times. Below is my code, I hope you can guide me. Thank you very much!
lsrtm.txt

[ar4]

Thank you for your code.

You seem to have an outer loop over shots. That is unusual. If this is intended to reduce memory requirements, then it is more typical to have the loop over shots within the loop over epochs, as in the reducing memory usage example.

I have made that change, and the change I proposed of slightly smoothing the true velocity model to try to reduce multiples, in the code below. Please let me know if it produces any improvement in the result.

import torch
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
import deepwave
from deepwave import scalar, scalar_born
import numpy as np

device = torch.device('cuda:0' if torch.cuda.is_available()
                      else 'cpu')
ny = 2301
nx = 801
dx = 4.0
nt = 750
dt = 0.004
n_shots = 115
freq = 26

v_true = torch.ones(ny, nx) * 1500
v_true[:, -751:] = torch.from_file('marmousi_vp.bin',
                         size=ny * 751).reshape(ny, 751)

# Smooth v_true slightly to reduce multiples
v_true = torch.tensor(1 / gaussian_filter(1 / v_true.numpy(), 5))

smooth2 = (torch.tensor(1 / gaussian_filter(1 / v_true.numpy(), 5))
           .to(device))

v_true = v_true.to(device)

n_sources_per_shot = 1
d_source = 20  # 20 * 4m = 80m
first_source = 10 # 10 * 4m = 40m
source_depth = 2  # 2 * 4m = 8m

n_receivers_per_shot = 384
d_receiver = 6  # 6 * 4m = 24m
first_receiver = 0  # 0 * 4m = 0m
receiver_depth = 2  # 2 * 4m = 8m

peak_time = 1.5 / freq

source_locations = torch.zeros(n_shots, n_sources_per_shot, 2,
                               dtype=torch.long, device=device)
source_locations[..., 1] = source_depth
source_locations[:, 0, 0] = (torch.arange(n_shots) * d_source +
                             first_source)

receiver_locations = torch.zeros(n_shots, n_receivers_per_shot, 2,
                                 dtype=torch.long, device=device)
receiver_locations[..., 1] = receiver_depth
receiver_locations[:, :, 0] = (
    (torch.arange(n_receivers_per_shot) * d_receiver +
     first_receiver)
    .repeat(n_shots, 1)
)

source_amplitudes = (
    (deepwave.wavelets.ricker(freq, nt, dt, peak_time))
    .repeat(n_shots, n_sources_per_shot, 1).to(device)
)

# Use the same max_vel value for each propagation to ensure consistent
# internal time steps (for better data match)
observed_true_data = scalar(v_true, dx, dt, source_amplitudes=source_amplitudes,
                            source_locations=source_locations,
                            receiver_locations=receiver_locations,
                            pml_freq=freq,
                            max_vel=v_true.max().item())[-1]
observed_smooth2_data = scalar(smooth2, dx, dt, source_amplitudes=source_amplitudes,
                            source_locations=source_locations,
                            receiver_locations=receiver_locations,
                            pml_freq=freq,
                            max_vel=v_true.max().item())[-1]
observed_scattered_data = observed_true_data - observed_smooth2_data

scatter = torch.zeros_like(smooth2)
scatter.requires_grad_()

optimiser = torch.optim.LBFGS([scatter])
loss_fn = torch.nn.MSELoss()

n_epochs = 10

for epoch in range(n_epochs):
   def closure():
       optimiser.zero_grad()
       total_loss = 0
       for shotidx in range(n_shots):
           out = scalar_born(
               smooth2, scatter, dx, dt,
               source_amplitudes=source_amplitudes[shotidx:shotidx+1],
               source_locations=source_locations[shotidx:shotidx+1],
               receiver_locations=receiver_locations[shotidx:shotidx+1],
               pml_freq=freq,
               max_vel=v_true.max().item()
           )
           loss = 1e20 * loss_fn(out[-1][0], observed_scattered_data[shotidx])
           loss.backward()
           total_loss += loss.item()
       print(epoch, total_loss)
       return total_loss

   optimiser.step(closure)