Questions on Implementing Conventional RTM and Wavefield Decomposition in DeepWave

[LiZHENGYui]

Dear Professor,

Thank you very much for developing and sharing the open-source DeepWave toolkit. I am a PhD student working on medical ultrasound imaging, and I have found DeepWave’s automatic differentiation capability to be extremely valuable for my research.

Currently, I am exploring velocity inversion by constraining RTM images. For this purpose, I would like to implement a conventional RTM based on wavefield cross-correlation. From my understanding, this might be approximated by performing a single iteration of LSRTM and extracting the gradient. However, I am uncertain how to apply illumination compensation to the resulting image in this case. Could you advise whether there is a recommended way to achieve this within DeepWave?

In addition, I am also interested in improving RTM image quality through wavefield decomposition. Would this require modifying the underlying wavefield propagation code, or is there an existing functionality in DeepWave that could support this?

I greatly appreciate your time and any guidance you can provide.

Best regards,

[ar4]

Dear @LiZHENGYui,

Thank you for your kind words. I am glad to hear that you find Deepwave useful.

Thank you also for your interesting questions.

You are correct that for an RTM image you can perform a single iteration of LSRTM. You can also just use the gradient from one iteration of FWI, which is more computationally efficient, and that is what I would normally recommend, but I don't know whether your desire to perform velocity inversion by constraining the RTM images affects that. On the one hand, you might need the separate background and scattered waves of the scalar born propagator to implement your idea (perhaps something like the joint migration-inversion example), but on the other hand maybe the scalar propagator's ability to perform double back propagation (so you can calculate gradients with respect to gradients) is what you need (double back propagation is not implemented for the other propagators). I think your current question, however, is about applying illumination compensation. The most recent release of Deepwave adds the ability to pass callback functions that are called during propagation. These are an efficient way to provide more flexibility to users, and might be useful for your needs. As part of the recent release, I rewrote the custom imaging condition example so that it now demonstrates how to use these callbacks to implement source-side illumination compensation. Does that enable you to do what you want?

For your second question, it sounds like you wish to modify the imaging condition further to use wavefields that you have decomposed in some way. I am hopeful that the new callback feature will also be useful for this application. You will probably need to store snapshots of the forward wavefield using a forward callback, and then use a backward callback to load the appropriate snapshot and then combine it with the current backward wavefield to implement your custom imaging condition. You will accumulate the gradient yourself, and then in the final backward time step (step 0) you can overwrite Deepwave's standard gradient with your own gradient. Does that sound like it will work, or will you need additional flexibility that the current implementation of the callback feature does not provide?

[LiZHENGYui]

Dear Professor, @ar4

Thank you so much for your timely and generous help—it really means a lot to me. Following your guidance, I tried out the custom_imaging_condition script you provided, and it almost perfectly resolved the issue I had raised! I’m truly excited about this. Of course, whether it can fully accomplish the work I need still requires further experimental validation, and I’ll make sure to share updates here as progress is made.

As you mentioned, my work is indeed similar to what you referred to as Joint Migration Inversion, or WEMVA in the seismic domain. This represents a rather cutting-edge direction in medical ultrasound. Conventional ultrasound imaging typically assumes a homogeneous background velocity model, whereas we are now exploring ways to improve the velocity distribution by correcting image discrepancies, thereby enhancing image quality. For students in ultrasound, deriving the adjoint-state equations and writing a near-perfect forward modeling code can be extremely challenging. That’s why I believe your DeepWave project is such a remarkable contribution—it will be an invaluable resource for researchers in medical ultrasound imaging.

Once again, thank you for your prompt and effective response. I truly hope DeepWave continues to grow and thrive—it is undoubtedly a great achievement. Wishing you good health and continued success in your work!

Best regards.

[ar4]

I am very happy to hear that the new callback feature is helpful to you. Please let me know if you encounter any more problems. I wish you luck in your research.

[LiZHENGYui]

Dear Professor, @ar4

Sorry to bother you again, but I’ve run into some issues with scalar wave equation forward modeling and would really appreciate your help.

Here’s my setup: I placed 128 sensors in a line along the top of the domain. Each forward simulation uses one element as the transmitter while all 128 elements act as receivers. I defined two velocity models, v0 (a homogeneous medium with velocity 1540 m/s) and v (the same background plus a single point target). My expectation was that by running forward modeling in both v and v0, and then subtracting the results, I could cancel out the direct arrivals (inter-element interference). However, the results are not as expected: the direct wave in the early part of the traces is not fully canceled.

From my understanding, the deepwave uses finite-difference simulation for forward modeling—so why doesn’t this subtraction completely remove the direct wave? Or is there something wrong with my code?

Another issue is that I observe diagonal artifacts in the channel data. My initial suspicion is that this comes from imperfect PML absorption. After extracting the signals, I found that the amplitude of these artifacts is still quite noticeable compared to the echo from the point target. Since my target application is medical ultrasound imaging (where velocity contrasts are relatively small), I’m concerned this interference could have a significant impact. Should I try increasing the PML thickness further, or are there other approaches you would recommend?

I’ve attached my code below.

import math
import numpy as np
import torch
import matplotlib.pyplot as plt
import scipy.io as sio
from scipy.signal import cheby1, lfilter
import deepwave
from deepwave import scalar
from skimage.transform import resize
from scipy.ndimage import gaussian_filter
from scipy.io import savemat

print(">>> Loaded from:", deepwave.__file__)

def tone_burst(f0, cycles, nt, dt, start_time=None, window="hann", peak_norm=True):
    t = np.arange(nt, dtype=np.float64) * dt
    T = cycles / f0
    if start_time is None:
        start_time = 1.0 / f0

    burst = np.zeros_like(t, dtype=np.float64)
    on = (t >= start_time) & (t < start_time + T)
    if np.any(on):
        tau = (t[on] - start_time) / T  # 0..1
        if window.lower() == "hann":
            win = 0.5 * (1.0 - np.cos(2.0 * np.pi * tau))
        else:
            win = np.ones_like(tau)
        burst[on] = np.sin(2.0 * np.pi * f0 * (t[on] - start_time)) * win

    if peak_norm and np.max(np.abs(burst)) > 0:
        burst = burst / np.max(np.abs(burst))
    return burst.astype(np.float32)

device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

n_points = 128
freq = 2e6                     # Hz
ImageZ = 0.08                  # m
ImageX = 0.04                  # m
v_init = 1540.0                # m/s
pitch = 0.0003                 # m

l_lambda = v_init / freq
dx = l_lambda / 10.0
dz = dx
nz = math.floor(ImageZ / dz)
nx = math.floor(ImageX / dx)

#  Velocity Model
# data = sio.loadmat("sound_speed_map.mat")
# sound_speed = data["sound_speed"].astype(np.float32)
# v_numpy = resize(sound_speed, (nz, nx), order=1, mode="reflect", anti_aliasing=True)
# v = torch.tensor(v_numpy, dtype=torch.float32, device=device)

# v_numpy = v.detach().cpu().numpy()
# v_smooth_numpy = gaussian_filter(v_numpy, sigma=10)
# v0 = torch.tensor(v_smooth_numpy, dtype=torch.float32, device=device)

#  Test Velocity Model: Point
z0, x0 = int(0.35*nz), int(0.5*nx)
Z = torch.arange(nz, device=device)[:,None]
X = torch.arange(nx, device=device)[None,:]
sigma = 1.2  # 1??1.5?
G = torch.exp(-0.5*((Z-z0)**2 + (X-x0)**2)/sigma**2)

v0 = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

v = v0.clone()
contrast = 0.03  # +3%
v = v * (1 + contrast*G) 

#  Test Velocity Model: Two Layer 
# v = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)
# v[nz//2 : , :] = 1600

# v0 = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

#  Plot Velocity Models
fig0, ax0 = plt.subplots(1, 2, figsize=(10, 5), sharey=True)

vmin = min(v.min(), v0.min())
vmax = max(v.max(), v.max())

im0 = ax0[0].imshow(v.detach().cpu().numpy(), cmap='jet', aspect='auto',
                    extent=[0, ImageX*1000, ImageZ*1000, 0],
                    vmin=vmin, vmax=vmax)
ax0[0].set_title("True Velocity Model (v)")
ax0[0].set_xlabel("X (mm)")
ax0[0].set_ylabel("Z (mm)")

im1 = ax0[1].imshow(v0.detach().cpu().numpy(), cmap='jet', aspect='auto',
                    extent=[0, ImageX*1000, ImageZ*1000, 0],
                    vmin=vmin, vmax=vmax)
ax0[1].set_title("Initial Velocity Model (v0)")
ax0[1].set_xlabel("X (mm)")

cbar = fig0.colorbar(im1, ax=ax0.ravel().tolist(), location="right", fraction=0.046, pad=0.04)
cbar.set_label("Velocity (m/s)")

# plt.tight_layout()
plt.show()

#  Array & Time 
n_shots = n_points
n_sources_per_shot = 1
d_source = math.floor(pitch / dx)

n_receivers_per_shot = n_points
d_receiver = math.floor(pitch / dx)

first_source   = (nx - 1 - (n_points - 1) * d_source)  // 2
first_receiver = (nx - 1 - (n_points - 1) * d_receiver) // 2

source_depth   = 10
receiver_depth = 10

dt = 0.5 * min(dx, dz) / (v.max().item() * np.sqrt(2.0))
nt = int(np.ceil(2.2 * ImageZ / v.min().item() / dt))

#  Geometry & Source 
# [n_shots, n_src_per_shot, 2] (x,z)
source_locations = torch.zeros(n_shots, n_sources_per_shot, 2, dtype=torch.long, device=device)
source_locations[..., 0] = source_depth
source_locations[:, 0, 1] = (torch.arange(n_shots, device=device) * d_source + first_source)

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

# Tone-burst
# cycles = 2
# burst_np = tone_burst(freq, cycles, nt, dt, start_time=1.0/freq, window="hann", peak_norm=True)
# source_amplitudes = torch.from_numpy(burst_np).to(device).repeat(n_shots, n_sources_per_shot, 1)

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


#  Forward Modeling
out  = scalar(v,  dx, dt, source_amplitudes=source_amplitudes,
              source_locations=source_locations, receiver_locations=receiver_locations,
              accuracy=8, pml_freq=freq, pml_width=80)
out0 = scalar(v0, dx, dt, source_amplitudes=source_amplitudes,
              source_locations=source_locations, receiver_locations=receiver_locations,
              accuracy=8, pml_freq=freq, pml_width=80)

# [n_shots, n_rx, nt]
rf_data = out[-1] - out0[-1]
# rf_data = out[-1]
# rf_data = rf_data/(dx*dx)

#  Save 
rf_data_np = rf_data.detach().cpu().numpy()
src_loc_np = source_locations.detach().cpu().numpy().astype(np.int32)
rcv_loc_np = receiver_locations.detach().cpu().numpy().astype(np.int32)
v_np  = v.detach().cpu().numpy()
v0_np = v0.detach().cpu().numpy()
src_amp_np = source_amplitudes.detach().cpu().numpy().astype(np.float32)

save_dict = {
    "rf_data": rf_data_np,              # [n_shots, n_receivers, nt]
    "source_locations": src_loc_np,     # [n_shots, n_src_per_shot, 2]
    "receiver_locations": rcv_loc_np,   # [n_shots, n_rx_per_shot, 2]
    "freq": np.array([freq], dtype=np.float32),
    "dt":   np.array([dt], dtype=np.float32),
    "nt":   np.array([nt], dtype=np.int32),
    "v":    v_np,                       # ????
    "v0":   v0_np,                       # ????
    "nz": np.array([nz], dtype=np.int32),
    "nx": np.array([nx], dtype=np.int32),
    "dx": np.array([dx], dtype=np.float32),
    "source_amplitudes": src_amp_np
}

savemat("Kidney_simulation_data.mat", save_dict)

#  Visuability: A-scan 
shot_idx = 64  
rx_idx   = 64  

trace = rf_data[shot_idx, rx_idx, :].cpu().numpy()
t = np.arange(nt) * dt * 1e6  # ??? us

plt.figure(figsize=(6,4))
plt.plot(t, trace)
plt.xlabel("Time (us)")
plt.ylabel("Amplitude")
plt.title(f"Shot {shot_idx}, Receiver {rx_idx}")
plt.show()
#  Visuability: Seismogram
shot_idx = 64
seismogram = rf_data[shot_idx].cpu().numpy()  # shape = [n_receivers, nt]

plt.figure(figsize=(10,6))
plt.imshow(seismogram.T, aspect='auto', cmap='jet',
           extent=[0, nt*dt*1e6, n_points, 0],
           vmin = -5e-6, vmax = 5e-6)  # time ??,?????
plt.xlabel("Time (us)")
plt.ylabel("Receiver index")
plt.title(f"Shot gather {shot_idx}")
plt.colorbar(label="Amplitude")
plt.show()

Looking forward to your reply, and thank you again for your time!

Best regards.

[ar4]

Dear @LiZHENGYui

Thank you for your questions. They are not a bother.

The reason for the incomplete subtraction is because if you do not specify the propagator's max_vel parameter, then Deepwave will use the actual maximum velocity in the velocity model to determine what time step interval it needs to use internally to satisfy the CFL condition for stability. Since your two velocity models have different maximum velocities, it is likely using different time step intervals for the two propagations, which leads to the small differences you see. If you specify the same value for the max_vel parameter for both of the propagations (I suggest v.max().item() in this case), then it will use the same time step interval for both and you should get better cancellation.

Your second question is more difficult to resolve. A small gaussian anomaly will cause very little reflection. The artifacts you see are indeed caused by incomplete absorption by the PML absorbing boundaries, but in reality the artifacts are quite small and only look large because the reflection is so small. When you use the max_vel parameter to get better subtraction of the direct arrivals, I think it should also subtract a lot of the artifacts in this case, so that might help. If you are able to make your anomaly more reflective, or to use a transmission setup (where the receiver is on the other side of the anomaly from the source), then that might also help. For simplicity I currently do not expose all of the PML parameters, but I will consider changing this in a future release, so that users who need to get optimal absorption can have greater control.

[LiZHENGYui]

Dear Professor, @ar4

Thank you very much for your timely and detailed reply. Your explanation resolves the issue I was facing, and I truly appreciate the clear guidance. The results now look much clearer.

I also look forward to seeing Deepwave continue to grow with even more powerful and versatile features in the future. Your efforts are making the tool increasingly valuable for research and applications.

Best regards！

[LiZHENGYui]

Dear Professor, @ar4

Sorry to bother you again, and thank you—your previous suggestions resolved my earlier issue. I’m now trying to produce more detailed RTM images, but I’m unsure whether the limitations I’m seeing come from the underlying imaging principle or from how I’m using DeepWave. If it’s not too much trouble, I would really appreciate your advice.

In fact, my goal is to use RTM to image a simple multi-layer velocity model. I tried two strategies: (1) using the FWI gradient directly as the RTM image, and (2) using the gradient with respect to m = 1/v² as the imaging condition. In both cases I also applied illumination compensation. Additionally, during forward modeling I generated data from a smooth/background velocity model to suppress direct arrivals.

However, the resulting images consistently exhibit a prominent arc-shaped artifact ahead of the first interface, with additional arc-like shadows at subsequent interfaces.

In medical ultrasound we usually view results in dB scale, but the image quality here is clearly suboptimal. By comparison, an RTM workflow I recently built on top of k-Wave forward modeling looks noticeably better (although, due to the nature of pseudospectral methods, removing direct arrivals is more challenging).

This issue has puzzled me for quite some time and I haven’t been able to resolve it. May I ask whether this is an inherent limitation of my RTM setup or whether I might be misusing DeepWave? Below are my forward-modeling and RTM scripts; any guidance you can share would be greatly appreciated. Thank you for your time again!

Forward Modeling Code

import math
import numpy as np
import torch
import matplotlib.pyplot as plt
import scipy.io as sio
from scipy.signal import cheby1, lfilter
import deepwave
from deepwave import scalar
from skimage.transform import resize
from scipy.ndimage import gaussian_filter
from scipy.io import savemat
import matplotlib.animation as animation


print(">>> Loaded from:", deepwave.__file__)

# ========================= Tone-burst (Hann) =========================
def tone_burst(f0, cycles, nt, dt, start_time=None, window="hann", peak_norm=True):
    t = np.arange(nt, dtype=np.float64) * dt
    T = cycles / f0
    if start_time is None:
        start_time = 1.0 / f0

    burst = np.zeros_like(t, dtype=np.float64)
    on = (t >= start_time) & (t < start_time + T)
    if np.any(on):
        tau = (t[on] - start_time) / T  # 0..1
        if window.lower() == "hann":
            win = 0.5 * (1.0 - np.cos(2.0 * np.pi * tau))
        else:
            win = np.ones_like(tau)
        burst[on] = np.sin(2.0 * np.pi * f0 * (t[on] - start_time)) * win

    if peak_norm and np.max(np.abs(burst)) > 0:
        burst = burst / np.max(np.abs(burst))
    return burst.astype(np.float32)

# ========================= Basic Config ==============================
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

n_points = 128
freq = 2e6                     # Hz
ImageZ = 0.04                  # m
ImageX = 0.04                  # m
v_init = 1540.0                # m/s
pitch = 0.0003                 # m

l_lambda = v_init / freq
dx = l_lambda / 10.0
dz = dx
nz = math.floor(ImageZ / dz)
nx = math.floor(ImageX / dx)

# ========================= Test Velocity Model: Multi Layers ============================
v = torch.ones(nz, nx, device=device) * v_init
for interface_idx in range(1, 5):
    v[interface_idx * int(nz/5) :] += 20.0

v0 = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

# ========================= Plot Velocity Models ======================
fig0, ax0 = plt.subplots(1, 2, figsize=(10, 5), sharey=True)

vmin = min(v.min(), v0.min())
vmax = max(v.max(), v.max())

im0 = ax0[0].imshow(v.detach().cpu().numpy(), cmap='jet', aspect='auto',
                    extent=[0, ImageX*1000, ImageZ*1000, 0],
                    vmin=vmin, vmax=vmax)
ax0[0].set_title("True Velocity Model (v)")
ax0[0].set_xlabel("X (mm)")
ax0[0].set_ylabel("Z (mm)")

im1 = ax0[1].imshow(v0.detach().cpu().numpy(), cmap='jet', aspect='auto',
                    extent=[0, ImageX*1000, ImageZ*1000, 0],
                    vmin=vmin, vmax=vmax)
ax0[1].set_title("Initial Velocity Model (v0)")
ax0[1].set_xlabel("X (mm)")

cbar = fig0.colorbar(im1, ax=ax0.ravel().tolist(), location="right", fraction=0.046, pad=0.04)
cbar.set_label("Velocity (m/s)")

# plt.tight_layout()
plt.show()

# ========================= Array & Time ==============================
n_shots = n_points
n_sources_per_shot = 1
d_source = math.floor(pitch / dx)

n_receivers_per_shot = n_points
d_receiver = math.floor(pitch / dx)

first_source   = (nx - 1 - (n_points - 1) * d_source)  // 2
first_receiver = (nx - 1 - (n_points - 1) * d_receiver) // 2

source_depth   = 0
receiver_depth = 0

dt = 0.5 * min(dx, dz) / (v.max().item() * np.sqrt(2.0))
# nt = int(np.ceil(2.2 * ImageZ / v.min().item() / dt))
nt = 3400

# ========================= Geometry & Source =========================
# [n_shots, n_src_per_shot, 2] (x,z)
source_locations = torch.zeros(n_shots, n_sources_per_shot, 2, dtype=torch.long, device=device)
source_locations[..., 0] = source_depth
source_locations[:, 0, 1] = (torch.arange(n_shots, device=device) * d_source + first_source)

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

# Tone-burst
cycles = 2
burst_np = tone_burst(freq, cycles, nt, dt, start_time=1.0/freq, window="hann", peak_norm=True)
source_amplitudes = torch.from_numpy(burst_np).to(device).repeat(n_shots, n_sources_per_shot, 1)

callback_frequency = 10
n_frames = (nt + callback_frequency - 1) // callback_frequency  # ceil(nt / freq)

forward_snapshots = torch.zeros(n_frames, nz, nx, dtype=torch.float32)
# forward_snapshots = torch.zeros(nt // callback_frequency, nz, nx)
class ForwardCallback:

    def __init__(self):
        """Initialises the ForwardCallback."""
        self.step = 0  # Forward propagation starts at time step 0

    def __call__(self, state):
        """A function called during the forward pass."""
        # We only need the wavefield for the first shot in the batch
        forward_snapshots[self.step] = (
            state.get_wavefield("wavefield_0")[0].cpu().clone()
        )
        self.step += 1

# ========================= Forward Modeling ==========================
out  = scalar(v,  dx, dt, source_amplitudes=source_amplitudes, max_vel = v.max(),
              source_locations=source_locations, receiver_locations=receiver_locations,
              accuracy=8, pml_freq=freq, pml_width=80,freq_taper_frac=0.2,
              time_pad_frac=0.2,
              forward_callback=ForwardCallback(), callback_frequency=callback_frequency)
out0 = scalar(v0, dx, dt, source_amplitudes=source_amplitudes, max_vel = v.max(),
              source_locations=source_locations, receiver_locations=receiver_locations,
              accuracy=8, pml_freq=freq, pml_width=80,freq_taper_frac=0.2,
              time_pad_frac=0.2)

# [n_shots, n_rx, nt]
rf_data_true = out[-1]
rf_data_back = out0[-1]

# ========================= Save ==========================
rf_data_np_true = rf_data_true.detach().cpu().numpy()
rf_data_np_back = rf_data_back.detach().cpu().numpy()
src_loc_np = source_locations.detach().cpu().numpy().astype(np.int32)
rcv_loc_np = receiver_locations.detach().cpu().numpy().astype(np.int32)
v_np  = v.detach().cpu().numpy()
v0_np = v0.detach().cpu().numpy()
src_amp_np = source_amplitudes.detach().cpu().numpy().astype(np.float32)

save_dict = {
    "rf_data_true": rf_data_np_true,              # [n_shots, n_receivers, nt]
    "rf_data_back": rf_data_np_back,              # [n_shots, n_receivers, nt]
    "source_locations": src_loc_np,     # [n_shots, n_src_per_shot, 2]
    "receiver_locations": rcv_loc_np,   # [n_shots, n_rx_per_shot, 2]
    "freq": np.array([freq], dtype=np.float32),
    "dt":   np.array([dt], dtype=np.float32),
    "nt":   np.array([nt], dtype=np.int32),
    "v":    v_np,                       # ????
    "v0":   v0_np,                       # ????
    "nz": np.array([nz], dtype=np.int32),
    "nx": np.array([nx], dtype=np.int32),
    "dx": np.array([dx], dtype=np.float32),
    "source_amplitudes": src_amp_np
}

savemat("Sim_MultiLayer_Data.mat", save_dict)

# ========================= Visuability: A-scan ==========================
shot_idx = 64  
rx_idx   = 64  

rf_data = rf_data_np_true - rf_data_np_back

trace = rf_data[shot_idx, rx_idx, :]
t = np.arange(nt) * dt * 1e6  # ??? us

plt.figure(figsize=(6,4))
plt.plot(t, trace)
plt.xlabel("Time (us)")
plt.ylabel("Amplitude")
plt.title(f"Shot {shot_idx}, Receiver {rx_idx}")
plt.show()
# ========================= Visuability: Seismogram ==========================
shot_idx = 64
seismogram = rf_data[shot_idx]  # shape = [n_receivers, nt]

plt.figure(figsize=(10,6))
plt.imshow(seismogram.T, aspect='auto', cmap='gray',
           extent=[0, n_points, nt*dt*1e6, 0])
           # vmin = -5e-6, vmax = 5e-6)  # time ??,?????
plt.xlabel("Receiver index")
plt.ylabel("Time (us)")
plt.title(f"Shot gather {shot_idx}")
plt.colorbar(label="Amplitude")
plt.show()

# ========================= Visuability: Animation ==========================
fig, ax = plt.subplots(figsize=(12, 4))  # ??:ax ????

# ??????(????max??)
v_abs = forward_snapshots.float().abs()
v_sample = v_abs[::2, ::4, ::4]  # ????,?????
vmax_forward = torch.quantile(v_sample.reshape(-1), 0.99).item()
# ???????:vmax_forward = v_abs.max().item()

im1 = ax.imshow(
    forward_snapshots[-1].cpu().numpy(),  # ? numpy ??
    vmin=-vmax_forward,
    vmax=vmax_forward,
    cmap="seismic",
    animated=True,
)
ax.set_xticks([]); ax.set_yticks([])

def init():
    im1.set_data(forward_snapshots[-1].cpu().numpy())
    return (im1,)  # ????????,??? im1 ??

def update(frame):
    idx = len(forward_snapshots) - 1 - frame
    im1.set_data(forward_snapshots[idx].cpu().numpy())
    ax.set_title(f"Forward ({idx * callback_frequency * dt:.3f} s)")
    return (im1,)   # ??????????

ani = animation.FuncAnimation(
    fig,
    update,
    frames=len(forward_snapshots),
    init_func=init,
    interval=10,
    blit=True,
)

ani.save("Forward_animation_new.gif", writer="pillow")

RTM Code

import math
import matplotlib.pyplot as plt
import torch
from scipy.ndimage import gaussian_filter
import deepwave
import numpy as np
import scipy.io as sio
from scipy.signal import cheby1, lfilter
import deepwave
from deepwave import scalar, scalar_born
from skimage.transform import resize
from scipy.ndimage import gaussian_filter
from scipy.io import savemat
from scipy.signal import hilbert


def _pad_edge(a, k=1):
    """replicate ????,?????? k ????"""
    if k == 0: 
        return a
    return np.pad(a, [(k, k)] * a.ndim, mode='edge')

def _iter_smooth_axis(x, axis, nit):
    """??????? nit ???: 0.25*? + 0.5*? + 0.25*??"""
    nit = int(nit)
    if nit <= 0:
        return x
    for _ in range(nit):
        x = 0.25 * (np.roll(x, -1, axis=axis) + np.roll(x,  1, axis=axis)) + 0.5 * x
    return x

def highpass_numpy(d_in, nit1, nit2=None, nit3=None, crop_mode="fixed"):
    """
    ????(1D/2D/3D)?????????(???? nit*),????? d_in - smooth?
    - crop_mode="fixed": ???????(??/???)
    - crop_mode="matlab": ?????????:3D ?????,3rd ???
    """
    a0 = np.asarray(d_in)
    nd = a0.ndim
    if nd not in (1, 2, 3):
        raise ValueError(f"Only 1D/2D/3D supported, got ndim={nd}")

    if nd == 1:
        n1 = int(nit1)
        if n1 == 0:
            return a0.copy()
        a = _pad_edge(a0, 1)
        a = _iter_smooth_axis(a, 0, n1)
        smooth = a[1:-1]  # 1D ????
        return a0 - smooth

    if nd == 2:
        n1 = int(nit1)
        n2 = n1 if nit2 is None else int(nit2)
        if n1 == 0 and n2 == 0:
            return a0.copy()
        a = _pad_edge(a0, 1)
        a = _iter_smooth_axis(a, 0, n1)
        a = _iter_smooth_axis(a, 1, n2)
        smooth = a[1:-1, 1:-1]
        return a0 - smooth

    # nd == 3
    n1 = int(nit1)
    n2 = n1 if nit2 is None else int(nit2)
    n3 = n1 if nit3 is None else int(nit3)
    if n1 == 0 and n2 == 0 and n3 == 0:
        return a0.copy()
    a = _pad_edge(a0, 1)
    a = _iter_smooth_axis(a, 0, n1)
    a = _iter_smooth_axis(a, 1, n2)
    a = _iter_smooth_axis(a, 2, n3)

    if crop_mode == "matlab":
        # ?????????? 3D ?????:?????,???????
        smooth = a[1:-1, 1:-1, :]
    else:
        # ???/???????:????
        smooth = a[1:-1, 1:-1, 1:-1]
    return a0 - smooth


# ========================= Basic Config ==============================
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

# data = sio.loadmat("Kidney_simulation_data.mat")
# data = sio.loadmat("Sim_Point_Data.mat")
data = sio.loadmat("Sim_MultiLayer_Data.mat")
print(data.keys())

n_shots = 128
freq = float(data["freq"].squeeze())
nz   = int(data["nz"].squeeze())
nx   = int(data["nx"].squeeze())
nt   = int(data["nt"].squeeze())
dx   = float(data["dx"].squeeze())
dt   = float(data["dt"].squeeze())
v_init = 1540

v_true = torch.tensor(data["v"],  dtype=torch.float32, device=device)
v_mig  = torch.tensor(data["v0"], dtype=torch.float32, device=device)
v0     = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

observed_data1 = torch.tensor(data["rf_data_true"], dtype=torch.float32, device=device)
background_data = torch.tensor(data["rf_data_back"], dtype=torch.float32, device=device)
observed_data = observed_data1 - background_data

source_locations   = torch.tensor(data["source_locations"],   dtype=torch.long, device=device)
receiver_locations = torch.tensor(data["receiver_locations"], dtype=torch.long, device=device)

source_amplitudes = torch.tensor(data["source_amplitudes"], dtype=torch.float32, device=device)

# ========================= Illumination Function ==============================
# def forward_cb(state):
#     nonlocal illum_batch
#     w = state.get_wavefield("wavefield_0", view="full").detach()  # (bs, nz, nx)
#     e = (w * w).sum(dim=0)  # -> (nz, nx)
#     illum_batch = e if illum_batch is None else (illum_batch + e)

# def backward_cb(state):
#     if state.step == 0:
#         eps = 1e-3 * illum_batch.max()
#         g = state.get_gradient("v", view="full")  # (nz, nx)
#         g /= (illum_batch + eps)

illumination = None
def forward_callback(state):
    global illumination
    if illumination is None:
        illumination = (
            state.get_wavefield("wavefield_0", view="full").detach() ** 2
        ).sum(dim=0)
    else:
        illumination += (
            state.get_wavefield("wavefield_0", view="full").detach() ** 2
        ).sum(dim=0)


def backward_callback(state):
    if state.step == 0:
        gradient = state.get_gradient("v", view="full")
        gradient /= illumination + 1e-3 * illumination.max()
# ========================= RTM ============================== 
v = v_mig.clone().detach().requires_grad_(True)
# # v_true.requires_grad_(True)
n_batch = 128  # The number of batches to use
n_shots_per_batch = (n_shots + n_batch - 1) // n_batch

# # (1) Without illumination compensation
# for batch in range(n_batch):
#     batch_start = batch * n_shots_per_batch
#     batch_end = min(batch_start + n_shots_per_batch, n_shots)
#     if batch_end <= batch_start:
#         continue
#     s = slice(batch_start, batch_end)
#     out = scalar(v, dx, dt,
#                  source_amplitudes=source_amplitudes[s],
#                  source_locations=source_locations[s], 
#                  receiver_locations=receiver_locations[s],
#                  accuracy=8, pml_freq=freq, pml_width=80, max_vel = v_true.max(),freq_taper_frac=0.2,
#                  time_pad_frac=0.2)
#     loss = -(out[-1] * observed_data[s]).sum() # rf_data = out[-1] - out0[-1]
#     # loss = torch.nn.MSELoss()(out[-1], observed_data[s]) # rf_data = out[-1]
#     loss.backward()

# fwi_grad = v.grad.detach().cpu().clone()
# v.grad.zero_()

# # (2) With illumination compensation
# loss_fn = torch.nn.MSELoss(reduction="mean")
# shots_per_batch = max(1, math.ceil(n_shots / n_batch))

# for b in range(n_batch):
#     start = b * shots_per_batch
#     end   = min(start + shots_per_batch, n_shots)
#     if end <= start:
#         continue
#     s  = slice(start, end)

#     holder = {"illum": None}  # ????,?? nonlocal

#     def forward_cb(state, holder=holder):
#         w = state.get_wavefield("wavefield_0", view="full").detach()  # (bs, nz, nx)
#         e = (w * w).sum(dim=0)                                        # (nz, nx)
#         holder["illum"] = e if holder["illum"] is None else holder["illum"] + e

#     def backward_cb(state, holder=holder):
#         if state.step == 0:
#             illum = holder["illum"]
#             # eps   = 1e-3 * illum.max()
#             g = state.get_gradient("v", view="full")                  # (nz, nx)
#             eps = (1e-3 * illum.max()).clamp_min(1e-12)   # ?? illum ?????
#             g /= (illum + eps)
#             # g /= illum

#     out = scalar(
#         v, dx, dt,
#         source_amplitudes=source_amplitudes[s],
#         source_locations=source_locations[s],
#         receiver_locations=receiver_locations[s],
#         accuracy=8, pml_freq=freq, pml_width=80, max_vel=v_true.max(),
#         forward_callback=forward_cb,
#         backward_callback=backward_cb,
#         callback_frequency=1,
#     )
    
#     # ????? backward ??????? == ???? shots ???????
#     # loss = loss_fn(out[-1], observed_data[s])
#     # loss = -(out[-1] * observed_data[s]).sum() # rf_data = out[-1] - out0[-1]
#     # loss = -((out[-1]-background_data[s]) * observed_data[s]).sum() 
#     loss = -((out[-1]-background_data[s]) * observed_data1[s]).sum() 
#     loss.backward()

# fwi_grad_comp = v.grad.detach().cpu().clone()
# v.grad.zero_()

# (3) Real RTM without illumination compensation
m = (1.0 / (v_mig**2)).clone().detach().requires_grad_(True)
# v_from_m = 1.0 / torch.sqrt(m)
# loss_fn = torch.nn.MSELoss(reduction="mean")
for batch in range(n_batch):
    batch_start = batch * n_shots_per_batch
    batch_end = min(batch_start + n_shots_per_batch, n_shots)
    if batch_end <= batch_start:
        continue
    s = slice(batch_start, batch_end)
    v_from_m = 1.0 / torch.sqrt(m)
    out = deepwave.scalar(
        v_from_m, dx, dt,
        source_amplitudes=source_amplitudes[s],
        source_locations=source_locations[s],
        receiver_locations=receiver_locations[s],
        accuracy=8, pml_freq=freq, pml_width=80, max_vel = v_true.max(),freq_taper_frac=0.2,
        time_pad_frac=0.2)
    loss = -(out[-1] * observed_data[s]).sum()
    # loss = loss_fn(out[-1], observed_data[s])
    loss.backward()
    
rtm_img_m = m.grad.detach().cpu().clone()
m.grad.zero_()

#  # (4) Born scatter without illumination compensation
# scatter = torch.zeros_like(v_mig)
# scatter.requires_grad_()
# loss_fn = torch.nn.MSELoss()
# n_epochs = 1
# optimiser = torch.optim.SGD([scatter], lr=1e9)
# for _epoch in range(n_epochs):
#     epoch_loss = 0
#     optimiser.zero_grad()
#     for batch in range(n_batch):
#         batch_start = batch * n_shots_per_batch
#         batch_end = min(batch_start + n_shots_per_batch, n_shots)
#         if batch_end <= batch_start:
#             continue
#         s = slice(batch_start, batch_end)
#         out = scalar_born(
#             v_mig, scatter, dx, dt,
#             source_amplitudes=source_amplitudes[s],
#             source_locations=source_locations[s],
#             receiver_locations=receiver_locations[s],
#             accuracy=8, pml_freq=freq, pml_width=80, max_vel = v_true.max(),freq_taper_frac=0.2,
#             time_pad_frac=0.2)
#         # out = out0[-1] - out1[-1]
#         # loss = loss_fn(out[-1], observed_data[s])
#         loss = -(out[-1] * observed_data[s]).sum()
#         epoch_loss += loss.item()
#         loss.backward()
#     optimiser.step()

# scatter_img = scatter.detach().cpu().clone()

# ========================= Visualization ============================== 
fig, ax = plt.subplots(figsize=(6, 5))
# im = ax.imshow(fwi_grad, cmap="gray", origin="upper", aspect="auto")
# im = ax.imshow(fwi_grad_comp, cmap="gray", origin="upper", aspect="auto")
im = ax.imshow(rtm_img_m, cmap="gray", origin="upper", aspect="auto")
# im = ax.imshow(scatter_img, cmap="gray", origin="upper", aspect="auto")
ax.set_title("RTM (wrt v)")
ax.set_xlabel("x"); ax.set_ylabel("z")
plt.colorbar(im, ax=ax, label="?L/?v (normalized range)")
plt.tight_layout()
plt.show()

# ========================= Visualization: B-mode ============================== 
from scipy.signal import hilbert

# img_temp = fwi_grad.cpu().numpy()  # ? fwi_grad_comp;?? numpy
# img = img_temp[:,68:449]

# img_temp = fwi_grad_comp.cpu().numpy()  # ? fwi_grad_comp;?? numpy
# img = img_temp[:,68:449]

img_temp = rtm_img_m.cpu().numpy()  # ? fwi_grad_comp;?? numpy
img = img_temp[:,68:449]

# img_temp = scatter_img.cpu().numpy()  # ? fwi_grad_comp;?? numpy
# img = img_temp[:,68:449]


# env = np.abs(img)
# plt.figure(figsize=(6,5))
# plt.imshow(env, cmap="gray", origin="upper", aspect="auto")
# plt.xlabel("x"); plt.ylabel("z")
# plt.colorbar(label="normalized intensity")
# plt.tight_layout(); plt.show()

img = highpass_numpy(img, 10,2)
env = np.abs(hilbert(img, axis=0))  # ? z ???
plt.figure(figsize=(6,5))
plt.imshow(env, cmap="gray", origin="upper", aspect="auto")
plt.xlabel("x"); plt.ylabel("z")
plt.colorbar(label="normalized intensity")
plt.tight_layout(); plt.show()

# eps = 1e-12
# env = env / (env.max() + eps)
# bmode_db = 20.0 * np.log10(env + eps)

# eps = 1e-12

env = env / env.max() 
bmode_db = 20.0 * np.log10(np.clip(env, 1e-12, None))

plt.figure(figsize=(6,5))
plt.imshow(bmode_db, cmap="gray", origin="upper", aspect="auto")
plt.xlabel("x"); plt.ylabel("z")
plt.colorbar(label="normalized intensity")
plt.tight_layout(); plt.show()


DR = 60.0                              # [-60, 0] dB ??
plt.figure(figsize=(6,5))
im = plt.imshow(bmode_db, cmap="gray", origin="upper", aspect="auto",
                vmin=-DR, vmax=0.0)   # ??:vmin/vmax ?? dB ??
plt.title(f"B-mode (DR={DR} dB)")
plt.xlabel("x"); plt.ylabel("z")
cb = plt.colorbar(im, label="Intensity (dB)")  # ??? dB
plt.tight_layout(); plt.show()

[ar4]

Hello again. I am sorry to hear that you have been having problems. Thank you for your thorough description of the issue, and for your code which looks nice.

I am unfortunately not an expert on getting good results with inversion, and your question also unfortunately comes at a time when I am especially busy, but I will try to look at your code in more detail soon.

My first comment is that the artifacts you see are probably a result of incomplete cancellation of the direct arrivals. The direct arrivals will not only be affected by the velocity at the surface, but will also be sensitive to the velocity near the surface. The layers near the surface will therefore change the direct arrivals and so the prediction will not cancel them completely. Using a smoothed version of the true velocity model, rather than a constant velocity model, to predict the direct arrivals might help to further attenuate the artifacts. Attempting to remove non-reflected arrivals in this way is quite a simplistic and crude approach, however. If you wish to obtain better results, there are more complicated approaches. One that is slightly more complicated, so probably still not good enough for your application, is described in this paper. It looks like you might have already tried something like it, with your high pass approach, but you should perhaps try the paper's method to see if it is better.

My second comment was going to be to suggest using a few iterations of LSRTM, but it looks like you might have already tried that as I see what looks like an LSRTM optimisation in the commented-out code. Did it improve the results at all?

Getting good results is difficult. It often requires careful regularisation and other advanced techniques that are unfortunately not topics I have much experience with.

[ar4]

I apologise again for the delay. I think I know what the problems might be.

1. You should subtract the predicted direct arrivals from the forward data in the loss function. You should change your line

loss = -(out[-1] * observed_data[s]).sum()

to

loss = -((out[-1] - background_data[s]) * observed_data[s]).sum()

2. I also think I was wrong about the cause of the artifacts you observed. I now think that they are probably finite aperture artifacts. You can attenuate these by applying a taper across the receiver data to taper the edges of the receiver data to zero. This can be achieved by creating a taper:

taper = torch.hann_window(
            observed_data.shape[1],
            periodic=False,
            device=device,
        ).reshape(1, -1, 1)

and then applying it to the data:

loss = -((out[-1] - background_data[s]) * observed_data[s] * taper).sum()

or, if you use another loss function, such as MSELoss:

loss = loss_fn((out[-1]-background_data[s])*taper, observed_data[s]*taper)

Please let me know if this helps.

[LiZHENGYui]

Dear Professor, @ar4

Thank you so much for the prompt and thoughtful guidance—fantastic! I really admire your clear reasoning and the practical way you approached the problem.

Following your suggestions, I (1) removed the predicted direct arrivals from the forward data and (2) applied an apodization taper across the receiver aperture (I used a window along the receiver dimension). This immediately improved the results: the persistent arc-shaped artifacts have essentially disappeared. It’s the first time I’ve seen, so clearly, how receiver-side weighting of a rectangular aperture benefits RTM—many thanks for pointing me in the right direction!

I’ll paste the updated code and the corresponding figures next. In the dB images, the arcs are gone, though there remain some faint horizontal line artifacts in the shallow region. My current guess is that they’re related to multiples or illumination effects; I’m still experimenting with mitigation strategies. Any brief pointers would be very much appreciated, but please don’t feel obliged.

import math
import matplotlib.pyplot as plt
import torch
from scipy.ndimage import gaussian_filter
import deepwave
import numpy as np
import scipy.io as sio
from scipy.signal import cheby1, lfilter
import deepwave
from deepwave import scalar, scalar_born
from skimage.transform import resize
from scipy.ndimage import gaussian_filter
from scipy.io import savemat
from scipy.signal import hilbert


def _pad_edge(a, k=1):
    """replicate ????,?????? k ????"""
    if k == 0: 
        return a
    return np.pad(a, [(k, k)] * a.ndim, mode='edge')

def _iter_smooth_axis(x, axis, nit):
    """??????? nit ???: 0.25*? + 0.5*? + 0.25*??"""
    nit = int(nit)
    if nit <= 0:
        return x
    for _ in range(nit):
        x = 0.25 * (np.roll(x, -1, axis=axis) + np.roll(x,  1, axis=axis)) + 0.5 * x
    return x

def highpass_numpy(d_in, nit1, nit2=None, nit3=None, crop_mode="fixed"):
    """
    ????(1D/2D/3D)?????????(???? nit*),????? d_in - smooth?
    - crop_mode="fixed": ???????(??/???)
    - crop_mode="matlab": ?????????:3D ?????,3rd ???
    """
    a0 = np.asarray(d_in)
    nd = a0.ndim
    if nd not in (1, 2, 3):
        raise ValueError(f"Only 1D/2D/3D supported, got ndim={nd}")

    if nd == 1:
        n1 = int(nit1)
        if n1 == 0:
            return a0.copy()
        a = _pad_edge(a0, 1)
        a = _iter_smooth_axis(a, 0, n1)
        smooth = a[1:-1]  # 1D ????
        return a0 - smooth

    if nd == 2:
        n1 = int(nit1)
        n2 = n1 if nit2 is None else int(nit2)
        if n1 == 0 and n2 == 0:
            return a0.copy()
        a = _pad_edge(a0, 1)
        a = _iter_smooth_axis(a, 0, n1)
        a = _iter_smooth_axis(a, 1, n2)
        smooth = a[1:-1, 1:-1]
        return a0 - smooth

    # nd == 3
    n1 = int(nit1)
    n2 = n1 if nit2 is None else int(nit2)
    n3 = n1 if nit3 is None else int(nit3)
    if n1 == 0 and n2 == 0 and n3 == 0:
        return a0.copy()
    a = _pad_edge(a0, 1)
    a = _iter_smooth_axis(a, 0, n1)
    a = _iter_smooth_axis(a, 1, n2)
    a = _iter_smooth_axis(a, 2, n3)

    if crop_mode == "matlab":
        # ?????????? 3D ?????:?????,???????
        smooth = a[1:-1, 1:-1, :]
    else:
        # ???/???????:????
        smooth = a[1:-1, 1:-1, 1:-1]
    return a0 - smooth


# ========================= Basic Config ==============================
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

# data = sio.loadmat("Kidney_simulation_data.mat")
# data = sio.loadmat("Sim_Point_Data.mat")
data = sio.loadmat("Sim_MultiLayer_Data.mat")
print(data.keys())

n_shots = 128
freq = float(data["freq"].squeeze())
nz   = int(data["nz"].squeeze())
nx   = int(data["nx"].squeeze())
nt   = int(data["nt"].squeeze())
dx   = float(data["dx"].squeeze())
dt   = float(data["dt"].squeeze())
v_init = 1540

v_true = torch.tensor(data["v"],  dtype=torch.float32, device=device)
v_mig  = torch.tensor(data["v0"], dtype=torch.float32, device=device)
v0     = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

observed_data1 = torch.tensor(data["rf_data_true"], dtype=torch.float32, device=device)
background_data = torch.tensor(data["rf_data_back"], dtype=torch.float32, device=device)
observed_data = observed_data1 - background_data

source_locations   = torch.tensor(data["source_locations"],   dtype=torch.long, device=device)
receiver_locations = torch.tensor(data["receiver_locations"], dtype=torch.long, device=device)

source_amplitudes = torch.tensor(data["source_amplitudes"], dtype=torch.float32, device=device)

# ========================= RTM ============================== 
v = v_true.clone().detach().requires_grad_(True)

grad_total = torch.zeros_like(v)
illum_total = torch.zeros_like(v)

n_batch = 128  # The number of batches to use
n_shots_per_batch = (n_shots + n_batch - 1) // n_batch

loss_fn = torch.nn.MSELoss(reduction="mean")
shots_per_batch = max(1, math.ceil(n_shots / n_batch))

taper = torch.hann_window(
            observed_data.shape[1],
            periodic=False,
            device=device,
        ).reshape(1, -1, 1)
# taper = torch.blackman_window(observed_data.shape[1], periodic=False, device=device).reshape(1, -1, 1)

for b in range(n_batch):
    start = b * shots_per_batch
    end   = min(start + shots_per_batch, n_shots)
    if end <= start:
        continue
    s  = slice(start, end)

    holder = {"illum": None}
    
    def forward_cb(state, holder=holder):
        w = state.get_wavefield("wavefield_0", view="inner").detach()
        e = (w * w).sum(dim=0)
        holder["illum"] = e if holder["illum"] is None else holder["illum"] + e
        
    def backward_cb(state, holder=holder):
        pass

    out = scalar(
        v, dx, dt,
        source_amplitudes=source_amplitudes[s],
        source_locations=source_locations[s],
        receiver_locations=receiver_locations[s],
        accuracy=8, pml_freq=freq, pml_width=80, max_vel=v_true.max(),
        forward_callback=forward_cb,
        backward_callback=backward_cb,
        callback_frequency=1,
    )
    
    loss = -((out[-1] - background_data[s]) * observed_data[s] * taper).sum()
    # loss = -((out[-1]) * observed_data[s] * taper).sum()
    # loss = loss_fn((out[-1])*taper, observed_data1[s]*taper)
    loss.backward()
    
    with torch.no_grad():
        grad_total += v.grad
        v.grad.zero_()   

        if holder["illum"] is not None:
            illum_total += holder["illum"]

with torch.no_grad():
    eps = 1e-2 * illum_total.max()
    grad_total = grad_total / (illum_total + eps)

plt.imshow(grad_total.cpu(), cmap="gray", origin="upper", aspect="auto")
plt.show()

# fwi_grad_comp = v.grad.detach().cpu().clone()
# v.grad.zero_()

# # ========================= Visualization ============================== 
# fig, ax = plt.subplots(figsize=(6, 5))
# # im = ax.imshow(fwi_grad, cmap="gray", origin="upper", aspect="auto")
# im = ax.imshow(fwi_grad_comp, cmap="gray", origin="upper", aspect="auto")
# # im = ax.imshow(rtm_img_m, cmap="gray", origin="upper", aspect="auto")
# # im = ax.imshow(scatter_img, cmap="gray", origin="upper", aspect="auto")
# ax.set_title("RTM (wrt v)")
# ax.set_xlabel("x"); ax.set_ylabel("z")
# plt.colorbar(im, ax=ax, label="?L/?v (normalized range)")
# plt.tight_layout()
# plt.show()

# # ========================= Visualization: B-mode ============================== 
# from scipy.signal import hilbert

img_temp = grad_total.cpu().numpy()  # ? fwi_grad_comp;?? numpy
img = img_temp[:,82:717]


# # img_temp = fwi_grad.cpu().numpy()  # ? fwi_grad_comp;?? numpy
# # img = img_temp[:,82:717]

# img_temp = fwi_grad_comp.cpu().numpy()  # ? fwi_grad_comp;?? numpy
# img = img_temp[:,82:717]

# # img_temp = rtm_img_m.cpu().numpy()  # ? fwi_grad_comp;?? numpy
# # img = img_temp[:,68:449]

# # img_temp = scatter_img.cpu().numpy()  # ? fwi_grad_comp;?? numpy
# # img = img_temp[:,68:449]


# # env = np.abs(img)
# # plt.figure(figsize=(6,5))
# # plt.imshow(env, cmap="gray", origin="upper", aspect="auto")
# # plt.xlabel("x"); plt.ylabel("z")
# # plt.colorbar(label="normalized intensity")
# # plt.tight_layout(); plt.show()

img = highpass_numpy(img, 10,2)
env = np.abs(hilbert(img, axis=0))  # ? z ???
# env = np.abs(img)  # ? z ???
plt.figure(figsize=(6,5))
plt.imshow(env, cmap="gray", origin="upper", aspect="auto")
plt.xlabel("x"); plt.ylabel("z")
plt.colorbar(label="normalized intensity")
plt.tight_layout(); plt.show()

env = env / env.max() 
bmode_db = 20.0 * np.log10(np.clip(env, 1e-12, None))

plt.figure(figsize=(6,5))
plt.imshow(bmode_db, cmap="gray", origin="upper", aspect="auto")
plt.xlabel("x"); plt.ylabel("z")
plt.colorbar(label="normalized intensity")
plt.tight_layout(); plt.show()


DR = 60.0                              # [-60, 0] dB ??
plt.figure(figsize=(6,5))
im = plt.imshow(bmode_db, cmap="gray", origin="upper", aspect="auto",
                vmin=-DR, vmax=0.0)   # ??:vmin/vmax ?? dB ??
plt.title(f"B-mode (DR={DR} dB)")
plt.xlabel("x"); plt.ylabel("z")
cb = plt.colorbar(im, label="Intensity (dB)")  # ??? dB
plt.tight_layout(); plt.show()

Thanks again for your detailed, patient reply and for sharing such practical tips. This has been tremendously helpful to me.

[ar4]

I'm glad to hear that there was an improvement. I have several more suggestions.

I think the biggest problem at the moment is that you seem to be using the same velocity model for migration and to generate the true data. The forward modelled data will therefore be almost identical to the observed data, and so there will be no residual for RTM to image.

I suggest that you create a migration model by smoothing the true model. You should use this migration model both to generate the background direct arrivals that you subtract from the true data, and as the velocity model for migration.

v0 = torch.tensor(1 / gaussian_filter(1 / v.numpy(), 20))

I suggest that you use the standard MSE loss function unless there is a good reason to use something else. If you do this, you should change the loss function calculation to:

loss = loss_fn((out[-1] - background_data[s])*taper, observed_data[s]*taper)

I hope this provides further improvements. Please let me know how it goes.

[LiZHENGYui]

Dear Professor, @ar4

Thank you again for your prompt and thoughtful reply.

Following your suggestions, I tried the following:

1. Background (direct-wave) prediction.
   I generated the background model using your slowness smoothing recipe

v0 = torch.tensor(1 / gaussian_filter(1 / v.numpy(), 20))

and used it to synthesize background data for direct-arrival removal. (The corresponding velocity model is shown in the figure.)

2. RTM loss with apodization.
   In the RTM step I used

loss = loss_fn((out[-1] - background_data[s]) * taper, observed_data[s] * taper)

Unfortunately, this correction did not materially change the 60 dB display—the main structure remained essentially the same.

However, when I expand the dynamic range to 160 dB, new behavior becomes visible: there are faint horizontal line artifacts between otherwise well-defined interfaces. Interestingly, these lines do not appear in k-Wave results (k-Wave images, without apodization on the observations, still show the familiar arc artifacts instead).

My current interpretation is that the two implementations produce systematically different low-energy tails: the regular 60 dB view hides them, but a 160 dB view reveals a clear divergence. I’m curious whether this stems from the pseudospectral vs. finite-difference discretization (e.g., dispersion, implicit periodicity, boundary treatment/PML) or from my use of DeepWave (e.g., how I apply apodization and illumination compensation).

I plan to investigate this more carefully and will share any progress here. If you have a brief comment on whether such shallow horizontal lines are a known consequence of the pseudospectral scheme or might indicate a configuration issue on my side, that would be extremely helpful—though please don’t feel obliged.

Thank you again for your timely and patient guidance. It has been tremendously helpful.

These are my ForwardModeling Code:

import math
import numpy as np
import torch
import matplotlib.pyplot as plt
import scipy.io as sio
from scipy.signal import cheby1, lfilter
import deepwave
from deepwave import scalar
from skimage.transform import resize
from scipy.ndimage import gaussian_filter
from scipy.io import savemat
import matplotlib.animation as animation


print(">>> Loaded from:", deepwave.__file__)

# ========================= Tone-burst (Hann) =========================
def tone_burst(f0, cycles, nt, dt, start_time=None, window="hann", peak_norm=True):
    t = np.arange(nt, dtype=np.float64) * dt
    T = cycles / f0
    if start_time is None:
        start_time = 1.0 / f0

    burst = np.zeros_like(t, dtype=np.float64)
    on = (t >= start_time) & (t < start_time + T)
    if np.any(on):
        tau = (t[on] - start_time) / T  # 0..1
        if window.lower() == "hann":
            win = 0.5 * (1.0 - np.cos(2.0 * np.pi * tau))
        else:
            win = np.ones_like(tau)
        burst[on] = np.sin(2.0 * np.pi * f0 * (t[on] - start_time)) * win

    if peak_norm and np.max(np.abs(burst)) > 0:
        burst = burst / np.max(np.abs(burst))
    return burst.astype(np.float32)

# ========================= Basic Config ==============================
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

n_points = 128
freq = 2e6                     # Hz
ImageZ = 0.04                  # m
ImageX = 0.04                  # m
v_init = 1540.0                # m/s
pitch = 0.0003                 # m

l_lambda = v_init / freq
# dx = l_lambda / 10.0
dx = 5e-05
dz = dx
nz = math.floor(ImageZ / dz)
nx = math.floor(ImageX / dx)

# ========================= Velocity Model ============================
# data = sio.loadmat("sound_speed_map.mat")
# sound_speed = data["sound_speed"].astype(np.float32)
# v_numpy = resize(sound_speed, (nz, nx), order=1, mode="reflect", anti_aliasing=True)
# v = torch.tensor(v_numpy, dtype=torch.float32, device=device)

# v_numpy = v.detach().cpu().numpy()
# v_smooth_numpy = gaussian_filter(v_numpy, sigma=10)
# v0 = torch.tensor(v_smooth_numpy, dtype=torch.float32, device=device)

# v = np.full((nz, nx), v_init, dtype=np.float32)
# anomaly_center = (nz // 2, nx // 2)
# anomaly_radius = nx // 8
# high_speed = 1560.0
# for iz in range(nz):
#     for ix in range(nx):
#         if (iz - anomaly_center[0])**2 + (ix - anomaly_center[1])**2 < anomaly_radius**2:
#             v[iz, ix] = high_speed
# v = torch.from_numpy(v).to(device)

# v0 = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

# ========================= Test Velocity Model: Point ============================
# z0, x0 = int(0.35*nz), int(0.5*nx)
# Z = torch.arange(nz, device=device)[:,None]
# X = torch.arange(nx, device=device)[None,:]
# sigma = 1.2  # 1??1.5?
# G = torch.exp(-0.5*((Z-z0)**2 + (X-x0)**2)/sigma**2)

# v0 = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

# v = v0.clone()
# contrast = 0.03  # +3%
# v = v * (1 + contrast*G) 

# ========================= Test Velocity Model: Multi Layers ============================
v = torch.ones(nz, nx, device=device) * v_init
for interface_idx in range(1, 5):
    v[interface_idx * int(nz/5) :] += 20.0

v0 = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

with torch.no_grad():
    inv = (1.0 / v).detach().cpu().numpy()                  # CPU numpy
    inv_smooth = gaussian_filter(inv, sigma=20, mode="nearest")  # SciPy ? CPU ????
    v0 = torch.as_tensor(1.0 / inv_smooth,
                         device=v.device,
                         dtype=v.dtype)                      # ????? device/dtype

# ========================= Test Velocity Model: Point ============================
# v = torch.ones(nz, nx, device=device) * v_init
# v[int(2*nz/3)-2:int(2*nz/3)+2,int(nx/2)-2:int(nx/2)+2] = 1600

# v0 = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

# ========================= Plot Velocity Models ======================
fig0, ax0 = plt.subplots(1, 2, figsize=(10, 5), sharey=True)

vmin = min(v.min(), v0.min())
vmax = max(v.max(), v.max())

im0 = ax0[0].imshow(v.detach().cpu().numpy(), cmap='jet', aspect='auto',
                    extent=[0, ImageX*1000, ImageZ*1000, 0],
                    vmin=vmin, vmax=vmax)
ax0[0].set_title("True Velocity Model (v)")
ax0[0].set_xlabel("X (mm)")
ax0[0].set_ylabel("Z (mm)")

im1 = ax0[1].imshow(v0.detach().cpu().numpy(), cmap='jet', aspect='auto',
                    extent=[0, ImageX*1000, ImageZ*1000, 0],
                    vmin=vmin, vmax=vmax)
ax0[1].set_title("Initial Velocity Model (v0)")
ax0[1].set_xlabel("X (mm)")

cbar = fig0.colorbar(im1, ax=ax0.ravel().tolist(), location="right", fraction=0.046, pad=0.04)
cbar.set_label("Velocity (m/s)")

# plt.tight_layout()
plt.show()

# ========================= Array & Time ==============================
n_shots = n_points
n_sources_per_shot = 1
d_source = math.floor(pitch / dx)

n_receivers_per_shot = n_points
d_receiver = math.floor(pitch / dx)

first_source   = (nx - 1 - (n_points - 1) * d_source)  // 2
first_receiver = (nx - 1 - (n_points - 1) * d_receiver) // 2

source_depth   = 0
receiver_depth = 0

# dt = 0.5 * min(dx, dz) / (v.max().item() * np.sqrt(2.0))
# # nt = int(np.ceil(2.2 * ImageZ / v.min().item() / dt))
# nt = 3400

dt = 9.270833333333323e-09
nt = 5000

# ========================= Geometry & Source =========================
# [n_shots, n_src_per_shot, 2] (x,z)
source_locations = torch.zeros(n_shots, n_sources_per_shot, 2, dtype=torch.long, device=device)
source_locations[..., 0] = source_depth
source_locations[:, 0, 1] = (torch.arange(n_shots, device=device) * d_source + first_source)

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

# Tone-burst
cycles = 2
burst_np = tone_burst(freq, cycles, nt, dt, start_time=1.0/freq, window="hann", peak_norm=True)
source_amplitudes = torch.from_numpy(burst_np).to(device).repeat(n_shots, n_sources_per_shot, 1)

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

callback_frequency = 10
n_frames = (nt + callback_frequency - 1) // callback_frequency  # ceil(nt / freq)

forward_snapshots = torch.zeros(n_frames, nz, nx, dtype=torch.float32)
# forward_snapshots = torch.zeros(nt // callback_frequency, nz, nx)
class ForwardCallback:

    def __init__(self):
        """Initialises the ForwardCallback."""
        self.step = 0  # Forward propagation starts at time step 0

    def __call__(self, state):
        """A function called during the forward pass."""
        # We only need the wavefield for the first shot in the batch
        forward_snapshots[self.step] = (
            state.get_wavefield("wavefield_0")[0].cpu().clone()
        )
        self.step += 1

# ========================= Forward Modeling ==========================
out  = scalar(v,  dx, dt, source_amplitudes=source_amplitudes, max_vel = v.max(),
              source_locations=source_locations, receiver_locations=receiver_locations,
              accuracy=8, pml_freq=freq, pml_width=80,freq_taper_frac=0.2,
              time_pad_frac=0.2,
              forward_callback=ForwardCallback(), callback_frequency=callback_frequency)
out0 = scalar(v0, dx, dt, source_amplitudes=source_amplitudes, max_vel = v.max(),
              source_locations=source_locations, receiver_locations=receiver_locations,
              accuracy=8, pml_freq=freq, pml_width=80,freq_taper_frac=0.2,
              time_pad_frac=0.2)

# [n_shots, n_rx, nt]
rf_data_true = out[-1]
rf_data_back = out0[-1]

# rf_data = out[-1]
# rf_data = rf_data/(dx*dx)

# ========================= Save ==========================
rf_data_np_true = rf_data_true.detach().cpu().numpy()
rf_data_np_back = rf_data_back.detach().cpu().numpy()
src_loc_np = source_locations.detach().cpu().numpy().astype(np.int32)
rcv_loc_np = receiver_locations.detach().cpu().numpy().astype(np.int32)
v_np  = v.detach().cpu().numpy()
v0_np = v0.detach().cpu().numpy()
src_amp_np = source_amplitudes.detach().cpu().numpy().astype(np.float32)

save_dict = {
    "rf_data_true": rf_data_np_true,              # [n_shots, n_receivers, nt]
    "rf_data_back": rf_data_np_back,              # [n_shots, n_receivers, nt]
    "source_locations": src_loc_np,     # [n_shots, n_src_per_shot, 2]
    "receiver_locations": rcv_loc_np,   # [n_shots, n_rx_per_shot, 2]
    "freq": np.array([freq], dtype=np.float32),
    "dt":   np.array([dt], dtype=np.float32),
    "nt":   np.array([nt], dtype=np.int32),
    "v":    v_np,                       # ????
    "v0":   v0_np,                       # ????
    "nz": np.array([nz], dtype=np.int32),
    "nx": np.array([nx], dtype=np.int32),
    "dx": np.array([dx], dtype=np.float32),
    "source_amplitudes": src_amp_np
}

savemat("Sim_MultiLayer_Data.mat", save_dict)

# ========================= Visuability: A-scan ==========================
shot_idx = 64  
rx_idx   = 64  

rf_data = rf_data_np_true - rf_data_np_back

trace = rf_data[shot_idx, rx_idx, :]
t = np.arange(nt) * dt * 1e6  # ??? us

plt.figure(figsize=(6,4))
plt.plot(t, trace)
plt.xlabel("Time (us)")
plt.ylabel("Amplitude")
plt.title(f"Shot {shot_idx}, Receiver {rx_idx}")
plt.show()
# ========================= Visuability: Seismogram ==========================
shot_idx = 64
seismogram = rf_data[shot_idx]  # shape = [n_receivers, nt]

plt.figure(figsize=(10,6))
plt.imshow(seismogram.T, aspect='auto', cmap='gray',
           extent=[0, n_points, nt*dt*1e6, 0])
           # vmin = -5e-6, vmax = 5e-6)  # time ??,?????
plt.xlabel("Receiver index")
plt.ylabel("Time (us)")
plt.title(f"Shot gather {shot_idx}")
plt.colorbar(label="Amplitude")
plt.show()

# ========================= Visuability: Animation ==========================
fig, ax = plt.subplots(figsize=(12, 4))  # ??:ax ????

# ??????(????max??)
v_abs = forward_snapshots.float().abs()
v_sample = v_abs[::2, ::4, ::4]  # ????,?????
vmax_forward = torch.quantile(v_sample.reshape(-1), 0.99).item()
# ???????:vmax_forward = v_abs.max().item()

im1 = ax.imshow(
    forward_snapshots[-1].cpu().numpy(),  # ? numpy ??
    vmin=-vmax_forward,
    vmax=vmax_forward,
    cmap="seismic",
    animated=True,
)
ax.set_xticks([]); ax.set_yticks([])

def init():
    im1.set_data(forward_snapshots[-1].cpu().numpy())
    return (im1,)  # ????????,??? im1 ??

def update(frame):
    idx = len(forward_snapshots) - 1 - frame
    im1.set_data(forward_snapshots[idx].cpu().numpy())
    ax.set_title(f"Forward ({idx * callback_frequency * dt:.3f} s)")
    return (im1,)   # ??????????

ani = animation.FuncAnimation(
    fig,
    update,
    frames=len(forward_snapshots),
    init_func=init,
    interval=10,
    blit=True,
)

ani.save("Forward_animation_new.gif", writer="pillow")

And this is my rtm code

import math
import matplotlib.pyplot as plt
import torch
from scipy.ndimage import gaussian_filter
import deepwave
import numpy as np
import scipy.io as sio
from scipy.signal import cheby1, lfilter
import deepwave
from deepwave import scalar, scalar_born
from skimage.transform import resize
from scipy.ndimage import gaussian_filter
from scipy.io import savemat
from scipy.signal import hilbert


def _pad_edge(a, k=1):
    """replicate ????,?????? k ????"""
    if k == 0: 
        return a
    return np.pad(a, [(k, k)] * a.ndim, mode='edge')

def _iter_smooth_axis(x, axis, nit):
    """??????? nit ???: 0.25*? + 0.5*? + 0.25*??"""
    nit = int(nit)
    if nit <= 0:
        return x
    for _ in range(nit):
        x = 0.25 * (np.roll(x, -1, axis=axis) + np.roll(x,  1, axis=axis)) + 0.5 * x
    return x

def highpass_numpy(d_in, nit1, nit2=None, nit3=None, crop_mode="fixed"):
    """
    ????(1D/2D/3D)?????????(???? nit*),????? d_in - smooth?
    - crop_mode="fixed": ???????(??/???)
    - crop_mode="matlab": ?????????:3D ?????,3rd ???
    """
    a0 = np.asarray(d_in)
    nd = a0.ndim
    if nd not in (1, 2, 3):
        raise ValueError(f"Only 1D/2D/3D supported, got ndim={nd}")

    if nd == 1:
        n1 = int(nit1)
        if n1 == 0:
            return a0.copy()
        a = _pad_edge(a0, 1)
        a = _iter_smooth_axis(a, 0, n1)
        smooth = a[1:-1]  # 1D ????
        return a0 - smooth

    if nd == 2:
        n1 = int(nit1)
        n2 = n1 if nit2 is None else int(nit2)
        if n1 == 0 and n2 == 0:
            return a0.copy()
        a = _pad_edge(a0, 1)
        a = _iter_smooth_axis(a, 0, n1)
        a = _iter_smooth_axis(a, 1, n2)
        smooth = a[1:-1, 1:-1]
        return a0 - smooth

    # nd == 3
    n1 = int(nit1)
    n2 = n1 if nit2 is None else int(nit2)
    n3 = n1 if nit3 is None else int(nit3)
    if n1 == 0 and n2 == 0 and n3 == 0:
        return a0.copy()
    a = _pad_edge(a0, 1)
    a = _iter_smooth_axis(a, 0, n1)
    a = _iter_smooth_axis(a, 1, n2)
    a = _iter_smooth_axis(a, 2, n3)

    if crop_mode == "matlab":
        # ?????????? 3D ?????:?????,???????
        smooth = a[1:-1, 1:-1, :]
    else:
        # ???/???????:????
        smooth = a[1:-1, 1:-1, 1:-1]
    return a0 - smooth


# ========================= Basic Config ==============================
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

data = sio.loadmat("Sim_MultiLayer_Data.mat")
print(data.keys())

n_shots = 128
freq = float(data["freq"].squeeze())
nz   = int(data["nz"].squeeze())
nx   = int(data["nx"].squeeze())
nt   = int(data["nt"].squeeze())
dx   = float(data["dx"].squeeze())
dt   = float(data["dt"].squeeze())
v_init = 1540

v_true = torch.tensor(data["v"],  dtype=torch.float32, device=device)
v_mig  = torch.tensor(data["v0"], dtype=torch.float32, device=device)
v0     = torch.full((nz, nx), v_init, dtype=torch.float32, device=device)

observed_data1 = torch.tensor(data["rf_data_true"], dtype=torch.float32, device=device)
background_data = torch.tensor(data["rf_data_back"], dtype=torch.float32, device=device)
observed_data = observed_data1 - background_data

source_locations   = torch.tensor(data["source_locations"],   dtype=torch.long, device=device)
receiver_locations = torch.tensor(data["receiver_locations"], dtype=torch.long, device=device)

source_amplitudes = torch.tensor(data["source_amplitudes"], dtype=torch.float32, device=device)

# ========================= RTM ============================== 
with torch.no_grad():
    inv = (1.0 / v_true).detach().cpu().numpy()                  # CPU numpy
    inv_smooth = gaussian_filter(inv, sigma=20, mode="nearest")  # SciPy ? CPU ????
    v0 = torch.as_tensor(1.0 / inv_smooth,
                         device=v_true.device,
                         dtype=v_true.dtype)                      # ????? device/dtype
# v0 = torch.tensor(1 / gaussian_filter(1 / v_true.numpy(), 20))
# v = v_mig.clone().detach().requires_grad_(True)
v = v0.clone().detach().requires_grad_(True)
# # v_true.requires_grad_(True)
n_batch = 128  # The number of batches to use
n_shots_per_batch = (n_shots + n_batch - 1) // n_batch

fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharex=True, sharey=True,
                         constrained_layout=True)
extent = [0.0, nx * dx, nz * dx, 0.0]
im0 = axes[0].imshow(v_true.detach().cpu().numpy(),
                     cmap="gray", origin="upper", aspect="auto",
                     vmin=v0.min(), vmax=v0.max(), extent=extent)
axes[0].set_title("v_true (m/s)")

im1 = axes[1].imshow(v0.detach().cpu().numpy(),
                     cmap="gray", origin="upper", aspect="auto",
                     vmin=v0.min(), vmax=v0.max(), extent=extent)
axes[1].set_title("v0 (m/s)")

# ??:? fig.colorbar,?? ax=axes ??(?? colorbar)
cbar = fig.colorbar(im1, ax=axes, location="right", shrink=0.9, pad=0.02)
cbar.set_label("Velocity (m/s)")

plt.show()

# (2) With illumination compensation
loss_fn = torch.nn.MSELoss(reduction="mean")
shots_per_batch = max(1, math.ceil(n_shots / n_batch))

taper = torch.hann_window(
            observed_data.shape[1],
            periodic=False,
            device=device,
        ).reshape(1, -1, 1)
# taper = torch.blackman_window(observed_data.shape[1], periodic=False, device=device).reshape(1, -1, 1)

for b in range(n_batch):
    start = b * shots_per_batch
    end   = min(start + shots_per_batch, n_shots)
    if end <= start:
        continue
    s  = slice(start, end)

    holder = {"illum": None}  # ????,?? nonlocal

    def forward_cb(state, holder=holder):
        w = state.get_wavefield("wavefield_0", view="full").detach()  # (bs, nz, nx)
        e = (w * w).sum(dim=0)                                        # (nz, nx)
        holder["illum"] = e if holder["illum"] is None else holder["illum"] + e

    def backward_cb(state, holder=holder):
        if state.step == 0:
            illum = holder["illum"]
            # eps   = 1e-3 * illum.max()
            g = state.get_gradient("v", view="full")                  # (nz, nx)
            eps = (1e-3 * illum.min()).clamp_min(1e-12)   # ?? illum ?????
            g /= (illum + eps)
            # g /= illum

    out = scalar(
        v, dx, dt,
        source_amplitudes=source_amplitudes[s],
        source_locations=source_locations[s],
        receiver_locations=receiver_locations[s],
        accuracy=8, pml_freq=freq, pml_width=80, max_vel=v_true.max(),
        forward_callback=forward_cb,
        backward_callback=backward_cb,
        callback_frequency=1,
    )
    
    loss = loss_fn((out[-1] - background_data[s])*taper, observed_data[s]*taper)
    loss.backward()

fwi_grad_comp = v.grad.detach().cpu().clone()
v.grad.zero_()

# ========================= Visualization ============================== 
fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(fwi_grad_comp, cmap="gray", origin="upper", aspect="auto")
ax.set_title("RTM (wrt v)")
ax.set_xlabel("x"); ax.set_ylabel("z")
plt.colorbar(im, ax=ax, label="?L/?v (normalized range)")
plt.tight_layout()
plt.show()

# ========================= Visualization: B-mode ============================== 
from scipy.signal import hilbert

img_temp = fwi_grad_comp.cpu().numpy()  # ? fwi_grad_comp;?? numpy
img = img_temp[:,82:717]

img = highpass_numpy(img, 10,2)
env = np.abs(hilbert(img, axis=0))  # ? z ???
plt.figure(figsize=(6,5))
plt.imshow(env, cmap="gray", origin="upper", aspect="auto")
plt.xlabel("x"); plt.ylabel("z")
plt.colorbar(label="normalized intensity")
plt.tight_layout(); plt.show()

env = env / env.max() 
bmode_db = 20.0 * np.log10(np.clip(env, 1e-12, None))

plt.figure(figsize=(6,5))
plt.imshow(bmode_db, cmap="gray", origin="upper", aspect="auto")
plt.xlabel("x"); plt.ylabel("z")
plt.colorbar(label="normalized intensity")
plt.tight_layout(); plt.show()


DR = 60.0                              # [-60, 0] dB ??
plt.figure(figsize=(6,5))
im = plt.imshow(bmode_db, cmap="gray", origin="upper", aspect="auto",
                vmin=-DR, vmax=0.0)   # ??:vmin/vmax ?? dB ??
plt.title(f"B-mode (DR={DR} dB)")
plt.xlabel("x"); plt.ylabel("z")
cb = plt.colorbar(im, label="Intensity (dB)")  # ??? dB
plt.tight_layout(); plt.show()

[ar4]

That is disappointing. It is also surprising, as when I run (approximately) your code I do not observe these artifacts:

I made some small modifications to the code to reduce the memory usage and runtime. I changed dt to dt = 1/(4*freq)and uncommented your code to calculate the number of time steps using nt = int(np.ceil(2.2 * ImageZ / v.min().item() / dt)). I also removed unnecessary options to the propagator, using code such as

out  = scalar(v,  dx, dt, source_amplitudes=source_amplitudes, max_vel = v.max(),
              source_locations=source_locations, receiver_locations=receiver_locations,
              pml_freq=freq,)

I tried it with and without illumination compensation, and did not see the artifacts in either case.

I hope that making these changes fixes the problem for you. If it does, then I suspect that it might be your use of time padding and frequency tapering when generating the true and background data that caused them.

[LiZHENGYui]

Dear Professor, @ar4

Thank you again for your prompt and patient reply. Inspired by your feedback, I experimented with various parameters of the scalar function and their impact on the final image quality, and I made some interesting observations that I’d like to share.

Following your suggestion, I focused on four parameters: pml_freq, pml_width, freq_taper_frac, and time_pad_frac. As expected, pml_freq should be set to the source center frequency. However, somewhat counterintuitively, when I increased pml_width to 80, the result actually became worse than with the default value. The attached figure shows this effect, which I found quite surprising.

Another interesting phenomenon appeared when I changed the imaging area. When I set the full imaging domain to 4 cm × 4 cm, the interfaces between different velocities were not well represented. But when I extended the imaging area to 4 cm × 8 cm, the behavior changed: the interfaces in the image became more slender, and when I restricted the displayed region to the vertical coverage of the array elements (keeping the element spacing and element count the same in both cases), the interface was better reconstructed on both sides in the 4 cm × 8 cm setting.

Anyway, I have achieved the effect I wanted. DeepWave has proven to be a truly remarkable tool! Thank you again for your generosity in open-sourcing it and for your patient guidance and responses over the days. Your dedication is deeply inspiring and deserves the respect of all researchers.

If there is ever anything I can help you with, please feel free to let me know—I would be very happy to contribute in any way I can.

Wishing you good health, continued success in your work, and many more outstanding achievements!

[ar4]

Dear @LiZHENGYui

Thank you for such a thorough and helpful analysis. Your findings are very interesting. Even though I know that seismic inversion/imaging can often behave in unexpected ways, the effect you observed from increasing the PML width is indeed surprising.

I noticed that the PML seems to be less effective with the source wavelet that you used compared to a Ricker wavelet. I wonder whether the "phantom layer" problem that you encountered when the PML width was 80 would still occur even when a Ricker wavelet was used.