简介面向地震勘探资料处理学习者与地球物理专业师生的MATLAB偏移成像代码合集聚焦校正反射点位置、消除波传播造成的成像畸变适合具备一定MATLAB基础、希望从脚本层面理解偏移算法实现的中高级读者。压缩包共24个文件全部为.m源码体积约50KB内容涵盖频率-波数域偏移、有限差分偏移、分裂步偏移、基尔霍夫偏移以及时间-深度转换、锥形窗、数据预处理等辅助模块结构紧凑便于按算法逐项研读。已有2015人学习下载说明其在课程实验与小规模数据验证中具有一定参考价值。读者可据此对照波动方程偏移的基本流程理解频域变换、波场延拓与成像条件在代码中的落地方式并借助不同脚本比较各类方法的适用条件为后续使用实际地震数据开展偏移实验、调试参数与可视化分析提供可复用的起点。1. 地震波偏移成像在 MATLAB 里跑通多数人卡在一个很具体的画面叠加剖面上还算连续的反射同相轴一偏移就变成一圈圈弧形假象深部糊成一片灰。偏移成像要干的事很朴素——按速度模型把被绕射和倾斜地层打散的能量重新放回它真正的空间位置。手上已有观测记录或正演数据、会写 MATLAB 循环但没系统做过波动方程成像的人可以照着搭出一套二维声波逆时偏移的最小闭环顺带知道哪几个参数一动就翻车。2. 有限差分逆时偏移的最小可跑闭环从速度模型到成像剖面在 MATLAB 里做地震波偏移成像核心从来不是「会不会写 FFT」而是能不能把四个环节串成一个自洽的循环正演震源波场、记录地表数据、把数据倒着注回模型、在相同时刻把两个波场相乘累加。这一段把可运行的最小版本直接写出来网格约定、量纲、存储策略、成像条件全部定死后面几章再往上加精度和存储优化。2.1 网格约定与速度模型先把坐标和量纲确定下来MATLAB 的矩阵索引是「行、列」地震波场则习惯「深度、水平距离」。两者能对上行索引对应深度 z 方向向下为正列索引对应水平 x 方向。一旦这个约定定下来就别再改我自己血的教训是中途换了一次行列约定反传波场和正传波场对不上成像结果靠近对称位置排查了大半天才找到。量纲统一用 SI速度 m/s网格间距 m时间 s频率 Hz。别混用 km 和 m也别让速度写 2.0 而网格写 5这类错误不会报错只会让成像结果整体偏移。速度模型对偏移来说必须是平滑的。偏移只负责把能量归位不负责反演速度把速度模型建得过于精细反而容易把绕射能量打散成像剖面反而更乱。常见做法是用一个背景速度加几个高斯型异常体既可控又能验证算法% 网格与速度模型行深度列距离单位统一 m nx 401; nz 201; dx 5; dz 5; % 网格间距 5 m x (0:nx-1)*dx; z (0:nz-1)*dz; [Z,X] ndgrid(z,x); % 生成 nz x nx 的坐标网格 % 背景 1800 m/s叠加一个高速高斯异常体 v 1800 700*exp(-((X-850).^2/(2*150^2) (Z-420).^2/(2*90^2))); v(1:24,:) 1500; % 顶部低速层模拟近地表速度场用高斯函数而不是矩形块原因很实际矩形块的四个角是速度的阶跃面正演时会在角点产生强烈的绕射成像剖面上会出现四条从角点发散出去的弧干扰判断。用平滑的异常体绕射能量集中在异常体本身看起来更接近真实地质体的响应。150和90分别是横向、纵向的高斯半径控制异常体的胖瘦调这两个数就能做出不同尺度的目标体。2.2 二维声波方程的 4 阶差分算子与单步更新函数二维常密度声波方程写成u_tt v²(u_xx u_zz)时间方向用 2 阶中心差分空间方向用 4 阶中心差分。空间阶数选 4 阶而不是 2 阶是精度和成本之间的折中2 阶差分在同样的网格间距下网格频散非常明显剖面上会出现一串沿传播方向拖尾的「波纹」4 阶差分的算子系数是(-1/12, 4/3, -5/2, 4/3, -1/12)同样网格下频散误差降低一个数量级以上代价只是每个点多重两个点在 MATLAB 里就是多写两行索引。function un step2d(u, up, v, dt, dx, dz, damp) % 二维声波方程单步推进时间 2 阶、空间 4 阶 % u : 当前时刻波场 [nz x nx] % up : 上一时刻波场 [nz x nx] % v : 速度模型 [nz x nx]m/s % damp : 每步乘的吸收系数 [nz x nx]内部为 1 c (v*dt).^2; L zeros(size(u)); L(3:end-2,3:end-2) ... (-1/12*(u(3:end-2,1:end-4) u(3:end-2,5:end)) ... 4/3 *(u(3:end-2,2:end-3) u(3:end-2,4:end-1)) ... - 5/2 * u(3:end-2,3:end-2)) / dx^2 ... (-1/12*(u(1:end-4,3:end-2) u(5:end,3:end-2)) ... 4/3 *(u(2:end-3,3:end-2) u(4:end-1,3:end-2)) ... - 5/2 * u(3:end-2,3:end-2)) / dz^2; un (2*u - up c.*L) .* damp;索引3:end-2是内点范围边界两圈留给吸收层不做差分计算。un计算出来之后马上乘以damp把边界附近的能量逐步衰减掉——这一步放在最后是为了让差分算子始终作用在未经衰减的场值上避免吸收系数污染内点。这个函数每调用一次推进一个时间步是整个程序里最热的代码。MATLAB 里做这类循环有几个提升速度的做法值得早点用上把u、up、v、damp全部预先转成singleL不要每次zeros重新分配改成一个在外层声明好的持久变量把c (v*dt).^2提到循环外面算好因为它不随时间变。2.3 正传记录与反传成像把记录倒着注回去逆时偏移的物理图像很直白正传时震源波场从震源向外扩散遇到界面产生反射反射被检波器记录下来反传时把记录当成震源从检波点位置按时间倒序注入波场会沿着原来的路径反向传播回去在反射界面处「聚焦」。把正传波场和反传波场在同一时刻、同一位置相乘再累加聚焦处乘积大其他位置乘积小成像条件就成立。观测系统设置如下。把震源和检波点放在同一条水平线上相当于一条海底电缆这样能避开自由表面产生的强烈多次波干扰把注意力集中在成像本身。参数取值说明震源深度isz30150 m位于吸收层内部避免源点被衰减震源横向isx2011000 m模型中心检波点深度rz30与震源同线检波点范围rxs25:4:377间距 20 m共 89 道主频f012 Hz决定分辨率也决定网格够不够记录长度nt4000 步约 1.96 s正传阶段把每一步的波场抽稀存储抽稀间隔dec 4。这一步是后面第 4 章存储优化的引子先按最省事的方式实现把成像闭环跑通。% ---- 正传记录数据并抽稀保存波场 ---- dec 4; % 每 4 步存一次正传波场 ns floor(nt/dec); Sf zeros(nz, nx, ns, single); rec zeros(nt, nrec, single); u1 zeros(nz,nx,single); u2 zeros(nz,nx,single); for it 1:nt un step2d(u1, u2, v, dt, dx, dz, damp); u2 u1; u1 un; u1(isz,isx) u1(isz,isx) src(it); % 点源注入 rec(it,:) u1(rz, rxs); if mod(it,dec) 0 Sf(:,:,it/dec) u1; end endu2 u1; u1 un;这两行的顺序不能反。step2d需要「当前」和「上一时刻」两个场推进完之后旧的当前场变成上一时刻场新算出来的场成为当前场。写反了程序不报错但波场会原地不动或者发散属于典型的静默失败。源注入直接加src(it)不乘dt²/(dx·dz)。严格来说这个常数因子是量纲需要的但成像条件和最终的归一化会把它消掉为了代码可读性先略去。如果要做振幅保真的反演类工作这个因子就必须补回去。反传阶段把记录倒序注入同时在抽稀的时间点上做互相关累加% ---- 反传记录倒序注入遇抽稀点做互相关 ---- ur1 zeros(nz,nx,single); ur2 zeros(nz,nx,single); img zeros(nz,nx); for it nt:-1:1 unr step2d(ur1, ur2, v, dt, dx, dz, damp); ur2 ur1; ur1 unr; ur1(rz, rxs) ur1(rz, rxs) rec(it,:); if mod(it,dec) 0 img img Sf(:,:,it/dec) .* ur1; % 零延迟互相关 end endmod(it,dec)0在倒序循环里同样成立it/dec从nt/dec递减到 1正好和正传时的存储顺序一一对应。这里只用了 1/4 的时刻做互相关成像项少了四分之三但因为互相关函数在时间上是平滑的少几个采样点带来的误差远小于视觉可辨的程度。真要较真把dec改成 1 再跑一遍对比会发现差异主要出现在深部弱能量区。2.4 照明补偿与 Laplace 滤波让深部能量能看直接算出来的img几乎不能看。两个原因一是几何扩散让深部能量比浅部弱好几个数量级画出来只有顶部一条亮线二是零延迟互相关成像条件本身会产生强烈的低频背景表现为剖面上大片的灰雾和界面处的「画弧」。照明补偿解决第一个问题。把正传波场的能量在时间上累加得到每个点的照明强度再用它去除成像结果illum sum(Sf.^2, 3); % 正传照明强度 img img ./ (illum 0.02*max(illum(:))); % 加小量防止除零分母里的0.02*max(illum(:))是稳定项取值太小吃不消深部极弱的照明除法会把噪声放大成雪花取得太大又退化成不做补偿。0.02是经验值速度横向变化剧烈时可以往上调到0.05。注意照明补偿只能压平能量不能凭空造出信息照明为零的地方补偿之后还是一片噪声这是方法的边界。第二个问题用 Laplace 滤波处理本质是高通把成像结果减去它自己的大尺度平滑版本留下界面附近的尖锐部分。k ones(15)/225; % 15x15 盒式平滑核 imgl conv2(conv2(img, k, same), k, same); imgf img - imgl; % 高通后的成像结果 imgf imgf / max p a hrefhttps://download.csdn.net/download/u014775170/8849083 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p
简介面向地震勘探资料处理学习者与地球物理专业师生的MATLAB偏移成像代码合集聚焦校正反射点位置、消除波传播造成的成像畸变适合具备一定MATLAB基础、希望从脚本层面理解偏移算法实现的中高级读者。压缩包共24个文件全部为.m源码体积约50KB内容涵盖频率-波数域偏移、有限差分偏移、分裂步偏移、基尔霍夫偏移以及时间-深度转换、锥形窗、数据预处理等辅助模块结构紧凑便于按算法逐项研读。已有2015人学习下载说明其在课程实验与小规模数据验证中具有一定参考价值。读者可据此对照波动方程偏移的基本流程理解频域变换、波场延拓与成像条件在代码中的落地方式并借助不同脚本比较各类方法的适用条件为后续使用实际地震数据开展偏移实验、调试参数与可视化分析提供可复用的起点。1. 地震波偏移成像在 MATLAB 里跑通多数人卡在一个很具体的画面叠加剖面上还算连续的反射同相轴一偏移就变成一圈圈弧形假象深部糊成一片灰。偏移成像要干的事很朴素——按速度模型把被绕射和倾斜地层打散的能量重新放回它真正的空间位置。手上已有观测记录或正演数据、会写 MATLAB 循环但没系统做过波动方程成像的人可以照着搭出一套二维声波逆时偏移的最小闭环顺带知道哪几个参数一动就翻车。2. 有限差分逆时偏移的最小可跑闭环从速度模型到成像剖面在 MATLAB 里做地震波偏移成像核心从来不是「会不会写 FFT」而是能不能把四个环节串成一个自洽的循环正演震源波场、记录地表数据、把数据倒着注回模型、在相同时刻把两个波场相乘累加。这一段把可运行的最小版本直接写出来网格约定、量纲、存储策略、成像条件全部定死后面几章再往上加精度和存储优化。2.1 网格约定与速度模型先把坐标和量纲确定下来MATLAB 的矩阵索引是「行、列」地震波场则习惯「深度、水平距离」。两者能对上行索引对应深度 z 方向向下为正列索引对应水平 x 方向。一旦这个约定定下来就别再改我自己血的教训是中途换了一次行列约定反传波场和正传波场对不上成像结果靠近对称位置排查了大半天才找到。量纲统一用 SI速度 m/s网格间距 m时间 s频率 Hz。别混用 km 和 m也别让速度写 2.0 而网格写 5这类错误不会报错只会让成像结果整体偏移。速度模型对偏移来说必须是平滑的。偏移只负责把能量归位不负责反演速度把速度模型建得过于精细反而容易把绕射能量打散成像剖面反而更乱。常见做法是用一个背景速度加几个高斯型异常体既可控又能验证算法% 网格与速度模型行深度列距离单位统一 m nx 401; nz 201; dx 5; dz 5; % 网格间距 5 m x (0:nx-1)*dx; z (0:nz-1)*dz; [Z,X] ndgrid(z,x); % 生成 nz x nx 的坐标网格 % 背景 1800 m/s叠加一个高速高斯异常体 v 1800 700*exp(-((X-850).^2/(2*150^2) (Z-420).^2/(2*90^2))); v(1:24,:) 1500; % 顶部低速层模拟近地表速度场用高斯函数而不是矩形块原因很实际矩形块的四个角是速度的阶跃面正演时会在角点产生强烈的绕射成像剖面上会出现四条从角点发散出去的弧干扰判断。用平滑的异常体绕射能量集中在异常体本身看起来更接近真实地质体的响应。150和90分别是横向、纵向的高斯半径控制异常体的胖瘦调这两个数就能做出不同尺度的目标体。2.2 二维声波方程的 4 阶差分算子与单步更新函数二维常密度声波方程写成u_tt v²(u_xx u_zz)时间方向用 2 阶中心差分空间方向用 4 阶中心差分。空间阶数选 4 阶而不是 2 阶是精度和成本之间的折中2 阶差分在同样的网格间距下网格频散非常明显剖面上会出现一串沿传播方向拖尾的「波纹」4 阶差分的算子系数是(-1/12, 4/3, -5/2, 4/3, -1/12)同样网格下频散误差降低一个数量级以上代价只是每个点多重两个点在 MATLAB 里就是多写两行索引。function un step2d(u, up, v, dt, dx, dz, damp) % 二维声波方程单步推进时间 2 阶、空间 4 阶 % u : 当前时刻波场 [nz x nx] % up : 上一时刻波场 [nz x nx] % v : 速度模型 [nz x nx]m/s % damp : 每步乘的吸收系数 [nz x nx]内部为 1 c (v*dt).^2; L zeros(size(u)); L(3:end-2,3:end-2) ... (-1/12*(u(3:end-2,1:end-4) u(3:end-2,5:end)) ... 4/3 *(u(3:end-2,2:end-3) u(3:end-2,4:end-1)) ... - 5/2 * u(3:end-2,3:end-2)) / dx^2 ... (-1/12*(u(1:end-4,3:end-2) u(5:end,3:end-2)) ... 4/3 *(u(2:end-3,3:end-2) u(4:end-1,3:end-2)) ... - 5/2 * u(3:end-2,3:end-2)) / dz^2; un (2*u - up c.*L) .* damp;索引3:end-2是内点范围边界两圈留给吸收层不做差分计算。un计算出来之后马上乘以damp把边界附近的能量逐步衰减掉——这一步放在最后是为了让差分算子始终作用在未经衰减的场值上避免吸收系数污染内点。这个函数每调用一次推进一个时间步是整个程序里最热的代码。MATLAB 里做这类循环有几个提升速度的做法值得早点用上把u、up、v、damp全部预先转成singleL不要每次zeros重新分配改成一个在外层声明好的持久变量把c (v*dt).^2提到循环外面算好因为它不随时间变。2.3 正传记录与反传成像把记录倒着注回去逆时偏移的物理图像很直白正传时震源波场从震源向外扩散遇到界面产生反射反射被检波器记录下来反传时把记录当成震源从检波点位置按时间倒序注入波场会沿着原来的路径反向传播回去在反射界面处「聚焦」。把正传波场和反传波场在同一时刻、同一位置相乘再累加聚焦处乘积大其他位置乘积小成像条件就成立。观测系统设置如下。把震源和检波点放在同一条水平线上相当于一条海底电缆这样能避开自由表面产生的强烈多次波干扰把注意力集中在成像本身。参数取值说明震源深度isz30150 m位于吸收层内部避免源点被衰减震源横向isx2011000 m模型中心检波点深度rz30与震源同线检波点范围rxs25:4:377间距 20 m共 89 道主频f012 Hz决定分辨率也决定网格够不够记录长度nt4000 步约 1.96 s正传阶段把每一步的波场抽稀存储抽稀间隔dec 4。这一步是后面第 4 章存储优化的引子先按最省事的方式实现把成像闭环跑通。% ---- 正传记录数据并抽稀保存波场 ---- dec 4; % 每 4 步存一次正传波场 ns floor(nt/dec); Sf zeros(nz, nx, ns, single); rec zeros(nt, nrec, single); u1 zeros(nz,nx,single); u2 zeros(nz,nx,single); for it 1:nt un step2d(u1, u2, v, dt, dx, dz, damp); u2 u1; u1 un; u1(isz,isx) u1(isz,isx) src(it); % 点源注入 rec(it,:) u1(rz, rxs); if mod(it,dec) 0 Sf(:,:,it/dec) u1; end endu2 u1; u1 un;这两行的顺序不能反。step2d需要「当前」和「上一时刻」两个场推进完之后旧的当前场变成上一时刻场新算出来的场成为当前场。写反了程序不报错但波场会原地不动或者发散属于典型的静默失败。源注入直接加src(it)不乘dt²/(dx·dz)。严格来说这个常数因子是量纲需要的但成像条件和最终的归一化会把它消掉为了代码可读性先略去。如果要做振幅保真的反演类工作这个因子就必须补回去。反传阶段把记录倒序注入同时在抽稀的时间点上做互相关累加% ---- 反传记录倒序注入遇抽稀点做互相关 ---- ur1 zeros(nz,nx,single); ur2 zeros(nz,nx,single); img zeros(nz,nx); for it nt:-1:1 unr step2d(ur1, ur2, v, dt, dx, dz, damp); ur2 ur1; ur1 unr; ur1(rz, rxs) ur1(rz, rxs) rec(it,:); if mod(it,dec) 0 img img Sf(:,:,it/dec) .* ur1; % 零延迟互相关 end endmod(it,dec)0在倒序循环里同样成立it/dec从nt/dec递减到 1正好和正传时的存储顺序一一对应。这里只用了 1/4 的时刻做互相关成像项少了四分之三但因为互相关函数在时间上是平滑的少几个采样点带来的误差远小于视觉可辨的程度。真要较真把dec改成 1 再跑一遍对比会发现差异主要出现在深部弱能量区。2.4 照明补偿与 Laplace 滤波让深部能量能看直接算出来的img几乎不能看。两个原因一是几何扩散让深部能量比浅部弱好几个数量级画出来只有顶部一条亮线二是零延迟互相关成像条件本身会产生强烈的低频背景表现为剖面上大片的灰雾和界面处的「画弧」。照明补偿解决第一个问题。把正传波场的能量在时间上累加得到每个点的照明强度再用它去除成像结果illum sum(Sf.^2, 3); % 正传照明强度 img img ./ (illum 0.02*max(illum(:))); % 加小量防止除零分母里的0.02*max(illum(:))是稳定项取值太小吃不消深部极弱的照明除法会把噪声放大成雪花取得太大又退化成不做补偿。0.02是经验值速度横向变化剧烈时可以往上调到0.05。注意照明补偿只能压平能量不能凭空造出信息照明为零的地方补偿之后还是一片噪声这是方法的边界。第二个问题用 Laplace 滤波处理本质是高通把成像结果减去它自己的大尺度平滑版本留下界面附近的尖锐部分。k ones(15)/225; % 15x15 盒式平滑核 imgl conv2(conv2(img, k, same), k, same); imgf img - imgl; % 高通后的成像结果 imgf imgf / max p a hrefhttps://download.csdn.net/download/u014775170/8849083 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p
