I use TR method to reconstruct the simulated data for layered media cases. However, the results seem exist a lot of artifacts. Here is the code

clearvars;

% =========================================================================

% SIMULATION

% =========================================================================

% create the computational grid

PML_size = 20; % size of the PML in grid points

Nx = 256 - 2 * PML_size; % number of grid points in the x direction

Ny = 256 - 2 * PML_size; % number of grid points in the y direction

dx = 0.1e-3; % grid point spacing in the x direction [m]

dy = 0.1e-3; % grid point spacing in the y direction [m]

kgrid = kWaveGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium

medium.sound_speed = 1500*ones(Nx,Ny);

medium.sound_speed(Nx/2:3*Nx/4,:) = 2900;

medium.density = 1000*ones(Nx,Ny);

medium.density(Nx/2:3*Nx/4,:) = 1900;

% create initial pressure distribution using makeDisc

disc_magnitude = 5; % [Pa]

disc_x_pos = 200; % [grid points]

disc_y_pos = 140; % [grid points]

disc_radius = 1; % [grid points]

disc_2 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_x_pos = 200; % [grid points]

disc_y_pos = 110; % [grid points]

disc_radius = 1; % [grid points]

disc_1 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

% smooth the initial pressure distribution and restore the magnitude

p0 = smooth(kgrid, disc_1 + disc_2, true);

% assign to the source structure

source.p0 = p0;

% define a binary line sensor

sensor.mask = zeros(Nx, Ny);

sensor.mask(1, :) = 1;

% create the time array

kgrid.makeTime(medium.sound_speed);

% set the input arguements: force the PML to be outside the computational

% grid; switch off p0 smoothing within kspaceFirstOrder2D

input_args = {'PMLInside', false, 'PMLSize', PML_size, 'Smooth', false, 'PlotPML', false};

% run the simulation

sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% reset the initial pressure

source.p0 = 0;

% assign the time reversal data

sensor.time_reversal_boundary_data = sensor_data;

% run the time reversal reconstruction

p0_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% add first order compensation for only recording over a half plane

p0_recon = 2 * p0_recon;