请用matlab写一个逆时偏移成像的例子

好的，下面是一个简单的逆时偏移成像的 matlab 代码示例：

%% 设置模拟参数
clear all; close all;
v = 2000;           % 声速
h = 100;            % 地表深度
dx = 10;            % 空间采样间隔
dt = 0.001;         % 时间采样间隔
tmax = 1;           % 最大时间
x = -500:dx:500;    % x方向采样点坐标
z = 0:dx:h;         % z方向采样点坐标
nt = round(tmax/dt);% 时间采样点数
nx = length(x);     % x方向采样点数
nz = length(z);     % z方向采样点数

%% 生成初始波场
f0 = 100;           % 主频
t = (0:nt-1)*dt;    % 时间坐标
src = sin(2*pi*f0*t);% 源波形
src = [src zeros(1, nx*nt-length(src))];% 补零
src = reshape(src, nt, nx);% 转换为矩阵形式
p = zeros(nz, nx, nt);% 初始化波场矩阵
p(:,:,1) = src;% 将源波形赋值给初始时刻波场矩阵

%% 逆时偏移成像
for it = 1:nt-1
    % 计算横向梯度
    px = diff(p(:,:,it), 1, 2)/dx;
    px(:,end+1) = 0; % 补足梯度矩阵最后一列
    % 计算纵向梯度
    pz = diff(p(:,:,it), 1, 1)/dx;
    pz(end+1,:) = 0; % 补足梯度矩阵最后一行
    % 横向和纵向梯度加权平方和
    grad2 = px.^2 + pz.^2;
    % 逆时偏移成像
    for ix = 1:nx
        for iz = 1:nz
            % 计算当前成像点在当前时刻的时间
            tcur = sqrt((x(ix)-x).^2 + (z(iz)-z).^2)/v;
            % 累积逆时偏移成像因子
            p(iz,ix,it+1) = p(iz,ix,it+1) + sum(src(:,ix).*src(:,ix).*grad2(iz,:)./(4*pi*(t(1:nt-1)+tcur(1:nt-1)).^3))*dt^2;
        end
    end
end

%% 显示逆时偏移成像结果
figure;
imagesc(x, z, squeeze(p(:,:,nt)));
xlabel('距离 (m)');
ylabel('深度 (m)');
title('逆时偏移成像结果');

代码中先生成一个简单的二维声波模型并产生一个初始波场。然后在每个时间步骤中计算横向和纵向梯度，再计算梯度平方和，并累积每个成像点的逆时偏移成像因子。最后显示逆时偏移成像结果。

注意：该代码仅为一个简单的示例，实际应用中需要根据具体问题进行参数调整和算法优化。