简介这份资源面向地震勘探数据处理方向的学习者与研究人员聚焦抛物线型多次波干扰的压制问题提供基于Radon正反变换的完整Matlab实现方案。Radon变换能将地震记录投影到不同方向使抛物线型多次波在变换域中呈现可分离的特征再通过滤波或置零处理最后经反Radon变换恢复至空间域从而提升地震资料的解析度。资源包共8个文件包含6个m脚本与2个su地震数据文件压缩包约218KB脚本覆盖正反Radon变换、频率域实现、segy读取及演示流程su文件则提供合成单炮记录与含多次波记录用于对比验证。目前已有755人学习下载适合希望理解Radon变换原理并动手实践多次波去除的读者参考可借助现成脚本快速复现处理流程观察去除前后效果差异为后续算法调优与资料处理打下基础。1. Radon 变换去多次波为什么它至今仍是地震资料处理里的硬通货做过海上地震资料处理的人都有一个共同体会多次波这东西你不管它它就在叠加剖面上糊成一片把真实的一次波反射结构盖得严严实实。而基于 radon 正反变换的多次波去除几乎是每个处理流程里绕不开的一步。它的核心思路其实不复杂——把时空域t-x 域的地震道集变换到 Radon 域τ-p 域或 τ-q 域在变换域里一次波和多次波因为动校正速度不同会被分离开压制掉多次波能量后再反变换回 t-x 域。Matlab 在这个环节里扮演的角色是绝大多数人做算法验证和参数扫描的首选工具矩阵操作天然适合 Radon 变换的离散实现可视化方便你实时看切除效果而且从算法原型到生产代码的迁移路径也清晰。这篇文章面向的是已经接触过地震资料处理基础、想自己动手把 Radon 去多次波跑通并调好参数的人——不管你是刚入行的处理员还是想从零搭一套验证流程的开发者下面从原理到代码到踩坑我按实际做的顺序讲一遍。2. Radon 变换的数学骨架与 Matlab 离散实现2.1 从连续 Radon 变换到离散线性算子连续形式的 Radon 变换本质上就是沿某条曲线做线积分。在地震数据处理里我们用的是线性 Radon 变换τ-p 变换和抛物 Radon 变换τ-q 变换两种。线性 Radon 沿直线积分抛物 Radon 沿抛物线积分。为什么是这两种因为反射地震的时距曲线在动校正后近似为抛物线而线性 Radon 对线性同相轴的聚焦更好。连续形式写出来是这样线性 Radon 正变换m(τ, p) ∫ d(t τ p·x, x) dx抛物 Radon 正变换m(τ, q) ∫ d(t τ q·x², x) dx其中 d(t, x) 是时空域数据m(τ, p) 或 m(τ, q) 是 Radon 域模型x 是偏移距p 是射线参数慢度q 是曲率参数。离散化之后问题变成求解一个线性方程组d L · md 是时空域数据向量所有道所有时间采样排成一列m 是 Radon 域模型向量L 是正变换算子矩阵。L 的每一个元素代表某条 Radon 曲线经过某个 (t, x) 采样点的权重。反变换就是求 m L⁻¹ · d。但 L 通常不是方阵而且条件数很大直接求逆会炸。所以实际用的是最小二乘或者迭代求解。2.2 Matlab 里怎么构建 L 矩阵这是整个流程里最核心的一步。L 矩阵构建得好不好直接决定后面反变换能不能收敛、分辨率够不够。function L build_radon_operator(t, x, p, mode) % BUILD_RADON_OPERATOR 构建 Radon 变换算子矩阵 % 输入: % t - 时间采样向量 (nt x 1) % x - 偏移距向量 (nx x 1) % p - 射线参数/曲率向量 (np x 1) % mode - linear 或 parabolic % 输出: % L - 正变换算子矩阵 (nt*nx x nt*np) nt length(t); nx length(x); np length(p); % 初始化稀疏矩阵用稀疏存储省内存 rows []; cols []; vals []; for ip 1:np for it 1:nt % 当前 Radon 域网格点对应的列索引 col_idx (ip - 1) * nt it; for ix 1:nx % 根据模式计算该偏移距上的到达时间 if strcmp(mode, linear) t_arrival t(it) p(ip) * x(ix); else % parabolic t_arrival t(it) p(ip) * x(ix)^2; end % 线性插值找到 t_arrival 在时间轴上的位置 if t_arrival t(1) t_arrival t(end) % 定位相邻采样点 idx_float (t_arrival - t(1)) / (t(2) - t(1)) 1; idx_low floor(idx_float); idx_high ceil(idx_float); if idx_low idx_high % 正好落在采样点上 row_idx (ix - 1) * nt idx_low; rows [rows; row_idx]; cols [cols; col_idx]; vals [vals; 1.0]; else % 线性插值权重 w_high idx_float - idx_low; w_low 1 - w_high; row_low (ix - 1) * nt idx_low; row_high (ix - 1) * nt idx_high; rows [rows; row_low; row_high]; cols [cols; col_idx; col_idx]; vals [vals; w_low; w_high]; end end end end end % 组装稀疏矩阵 L sparse(rows, cols, vals, nt*nx, nt*np); end这段代码的逻辑说明外层双重循环遍历 Radon 域的所有网格点 (τ, p)对每个网格点再遍历所有偏移距计算该 Radon 参数对应的到达时间然后用线性插值把能量映射到时间轴的相邻采样点上。最终组装成一个稀疏矩阵。参数说明t的时间采样间隔决定了插值精度一般用 2ms 或 4msp的范围要根据实际数据的最大慢度来定太小了分不开太大了浪费计算量mode选 linear 还是 parabolic 取决于你的数据动校正后的同相轴形态。注意这个双重循环在 nt 和 np 较大时非常慢。实际使用中建议把内层循环向量化或者用 Matlab 的sparse函数直接批量构建索引。我一般会预分配 rows/cols/vals 数组而不是动态增长。2.3 反变换求解为什么不能直接求逆有了 L 之后反变换就是解Lm d。但 L 的条件数通常在 10⁶ 以上直接m L \ d会得到一堆震荡解。常见做法有三种第一种是阻尼最小二乘Tikhonov 正则化m (LᵀL λI)⁻¹ Lᵀd。λ 是阻尼因子控制解的平滑程度。λ 太大分辨率低太小噪声放大。第二种是共轭梯度迭代CGLS适合大规模稀疏矩阵不用显式求逆。第三种是频域实现把每个频率切片单独做 Radon 变换计算量小很多也是生产中最常用的。function m solve_radon_inverse(L, d, lambda, max_iter) % SOLVE_RADON_INVERSE 阻尼最小二乘求解 Radon 反变换 % 输入: % L - 正变换算子 (稀疏矩阵) % d - 时空域数据向量 % lambda - 阻尼因子典型值 0.01 ~ 1.0 % max_iter - 最大迭代次数 % 输出: % m - Radon 域模型向量 % 构建法方程矩阵 A L * L lambda * speye(size(L, 2)); b L * d; % 用共轭梯度法求解 A*m b m zeros(size(L, 2), 1); r b - A * m; p r; rsold r * r; for iter 1:max_iter Ap A * p; alpha rsold / (p * Ap); m m alpha * p; r r - alpha * Ap; rsnew r * r; % 收敛判断 if sqrt(rsnew) 1e-6 fprintf(CG 在第 %d 次迭代收敛\n, iter); break; end p r (rsnew / rsold) * p; rsold rsnew; end end逻辑说明这里没有直接对 L 求逆而是构建法方程 LᵀL λI然后用共轭梯度迭代求解。这样做的好处是避免了显式求逆的数值不稳定同时稀疏矩阵乘法效率高。参数说明lambda是关键参数我一般从 0.1 开始试看 Radon 域的能量聚焦程度再调max_iter通常 50 到 200 次就够了再多收益很小。3. 从合成记录到实际数据完整去多次波流程3.1 用合成数据验证正反变换的正确性在拿实际数据开刀之前先用合成记录验证你的算子和求解器对不对。这一步很多人跳过结果到了实际数据上出了问题根本分不清是算法错了还是数据太差。% 生成合成地震记录一个一次波 一个多次波 clear; clc; nt 512; % 时间采样点数 dt 0.002; % 采样间隔 2ms t (0:nt-1) * dt; nx 60; % 道数 x (0:nx-1) * 25; % 偏移距 0 到 1475m % 一次波t0 0.3s动校正速度 2000m/s t0_primary 0.3; v_primary 2000; t_primary sqrt(t0_primary^2 (x / v_primary).^2); % 多次波t0 0.6s视速度更低 t0_multiple 0.6; v_multiple 1400; t_multiple sqrt(t0_multiple^2 (x / v_multiple).^2); % 生成数据 d zeros(nt, nx); for ix 1:nx % 一次波Ricker 子波 d(:, ix) d(:, ix) ricker_wavelet(t - t_primary(ix), 0.01); % 多次波幅度 0.6 倍 d(:, ix) d(:, ix) 0.6 * ricker_wavelet(t - t_multiple(ix), 0.01); end % 定义 Ricker 子波函数 function w ricker_wavelet(t, f_dom) % t: 时间向量相对于子波中心 % f_dom: 主频 w (1 - 2 * (pi * f_dom * t).^2) .* exp(-(pi * f_dom * t).^2); end这段代码生成了一个包含一次波和多次波的合成道集。一次波的动校正速度是 2000m/s多次波视速度 1400m/s两者在 t-x 域里是交叉的但在 Radon 域里会落到不同的 p 值上。3.2 Radon 域切除怎么切、切多少变换到 Radon 域之后一次波和多次波会分布在不同的 p 值区间。切除就是在 Radon 域里把多次波对应的能量置零或者衰减。% 接上面的合成数据 % 构建 Radon 算子 p linspace(-0.0005, 0.0005, 100); % 慢度范围 L build_radon_operator(t, x, p, linear); % 正变换时空域 - Radon 域 d_vec d(:); % 展平 m solve_radon_inverse(L, d_vec, 0.05, 100); m_2d reshape(m, nt, length(p)); % 在 Radon 域做切除 % 一次波在 p 接近 0 的区域多次波在 |p| 较大的区域 p_cut 0.0002; % 切除阈值需要根据实际数据调整 mask ones(size(m_2d)); for ip 1:length(p) if abs(p(ip)) p_cut mask(:, ip) 0; % 压制大 p 值区域的能量 end end m_cut m_2d .* mask; % 反变换Radon 域 - 时空域 d_cut_vec L * m_cut(:); d_cut reshape(d_cut_vec, nt, nx); % 对比去多次波前后的道集 figure; subplot(1,3,1); imagesc(x, t, d); title(原始道集); xlabel(偏移距(m)); ylabel(时间(s)); subplot(1,3,2); imagesc(p, t, m_2d); title(Radon 域); xlabel(p (s/m)); ylabel(τ (s)); subplot(1,3,3); imagesc(x, t, d_cut); title(去多次波后); xlabel(偏移距(m)); ylabel(时间(s));逻辑说明先把时空域数据展平成一维向量求解 Radon 反变换得到 Radon 域模型然后在 Radon 域里根据 p 值做切除最后用正变换算子把切除后的模型映射回时空域。参数说明p_cut是最关键的参数。太小了会伤到一次波太大了多次波压不干净。实际操作中我会先看 Radon 域的能量分布图找到一次波和多次波之间的分界 p 值然后在这个值附近做几组扫描测试。3.3 参数扫描p 范围、阻尼因子、切除阈值怎么定这三个参数是耦合的不能单独调。我一般按这个顺序来先定 p 范围。看数据的最大偏移距和最深目的层时间估算最大慢度 p_max ≈ x_max / (v_min · t_max)。p 的采样数一般取 80 到 150太少了 Radon 域分辨率不够太多了计算量爆炸。再定阻尼因子 λ。从 0.01 开始每次乘 2 往上试。判断标准是 Radon 域的能量是否聚焦成团——如果散成一片说明 λ 太小如果糊成一坨说明 λ 太大。最后定切除阈值。这个必须结合 Radon 域的能量分布来定。我通常会把切除前后的道集并排显示反复对比确保一次波没被伤到。% 参数扫描示例 lambda_list [0.01, 0.05, 0.1, 0.5]; p_cut_list [0.0001, 0.0002, 0.0003]; for il 1:length(lambda_list) for ic 1:length(p_cut_list) m_test solve_radon_inverse(L, d_vec, lambda_list(il), 100); m_test_2d reshape(m_test, nt, length(p)); % 应用切除 mask_test ones(size(m_test_2d)); for ip 1:length(p) if abs(p(ip)) p_cut_list(ic) mask_test(:, ip) 0; end end d_test reshape(L * (m_test_2d .* mask_test)(:), nt, nx); % 计算残余能量比 residual_ratio norm(d_test(:) - d(:)) / norm(d(:)); fprintf(lambda%.2f, p_cut%.4f, 残余比%.3f\n, ... lambda_list(il), p_cut_list(ic), residual_ratio); end end这段扫描代码会输出不同参数组合下的残余能量比。残余比不是越小越好——太小了可能把一次波也去掉了。理想情况是残余比在 0.3 到 0.5 之间说明多次波被压制了但一次波保留完好。4. 避坑指南Radon 去多次波最容易翻车的五个地方4.1 坑一空间假频导致 Radon 域出现能量拖尾现象Radon 域里一次波和多次波的能量不是聚焦成团而是拖成一条一条的尾巴切除时根本分不清哪个是哪个。原因偏移距采样太稀空间假频严重。Radon 变换对空间采样率非常敏感当偏移距间隔大于最小视波长的一半时变换域就会出现假频能量。解决在 Radon 变换之前先做道插值把偏移距加密到原来的 2 到 4 倍。或者改用高分辨率 Radon 变换比如用稀疏约束代替 L2 正则它对假频的容忍度更高。4.2 坑二切除太狠伤了一次波现象去多次波后的叠加剖面上浅层反射的振幅明显变弱甚至出现空白。原因一次波和多次波在 Radon 域里并不是完全分离的在近偏移距区域两者有重叠。切除阈值设得太小把重叠区域的一次波能量也一起干掉了。解决不要用硬切除置零改用软切除加权衰减。在重叠区域用平滑的过渡带让一次波的能量逐渐衰减而不是突然消失。% 软切除用余弦过渡带代替硬置零 p_cut_inner 0.00015; % 完全保留的 p 范围 p_cut_outer 0.00025; % 完全压制的 p 范围 mask_soft ones(size(m_2d)); for ip 1:length(p) abs_p abs(p(ip)); if abs_p p_cut_inner mask_soft(:, ip) 1; elseif abs_p p_cut_outer mask_soft(:, ip) 0; else % 余弦过渡 ratio (abs_p - p_cut_inner) / (p_cut_outer - p_cut_inner); mask_soft(:, ip) 0.5 * (1 cos(pi * ratio)); end end4.3 坑三反变换后振幅不一致现象去多次波后的道集和原始道集对比一次波的振幅对不上有的道强有的道弱。原因Radon 正反变换不是保幅变换。正变换算子 L 和反变换算子 L⁻¹ 的乘积不是单位矩阵而是一个近似。在 Radon 域做切除后这个近似误差会被放大。解决在反变换之后做一次振幅补偿用原始道集和去多次波后道集在一次波窗口内的振幅比作为补偿因子。或者用保幅的 Radon 变换实现在算子构建时考虑几何扩散补偿。4.4 坑四计算量太大跑不动现象nt1000、nx200、np150 的数据构建 L 矩阵就要好几分钟求解更是要等半小时以上。原因L 矩阵的大小是 (nt·nx) × (nt·np)上面的例子就是 200000 × 150000即使稀疏存储也占很大内存。而且 CG 迭代每次都要做两次稀疏矩阵乘法。解决改用频域实现。对每个频率切片单独做 Radon 变换矩阵大小降到 nx × np计算量降低几个数量级。这也是生产中的标准做法。% 频域 Radon 变换的核心思路 % 对每个频率 f构建 L_f (nx x np)求解 m_f L_f \ d_f % 然后对每个频率反变换回时间域 freqs linspace(0, 125, nt/21); % 频率向量 D_freq fft(d, [], 1); % 时间方向 FFT for ifreq 1:length(freqs) f freqs(ifreq); % 构建频率域算子L_f(ix, ip) exp(-1i*2*pi*f*p(ip)*x(ix)) L_f exp(-1i * 2 * pi * f * x * p); % 求解 d_f D_freq(ifreq, :); m_f (L_f * L_f 0.01 * eye(length(p))) \ (L_f * d_f); % 切除 % ... (在频率域做切除) % 反变换 D_cut(ifreq, :) (L_f * m_f); end d_cut real(ifft([D_cut; conj(D_cut(end-1:-1:2, :))], [], 1));4.5 坑五多次波周期估计错误导致切除位置偏移现象Radon 域里明明看到多次波能量了但切除之后叠加剖面上多次波还在。原因多次波在 Radon 域的位置和它的周期有关。如果动校正速度用错了多次波在 Radon 域的位置会偏移你按理论值切的其实不是多次波。解决不要死算理论值直接在 Radon 域里看能量团的实际位置。把 Radon 域的能量图放大找到多次波能量团的中心 p 值用这个实测值来定切除阈值。另外动校正速度的精度也很关键速度误差超过 5% 就会导致 Radon 域聚焦变差。5. 高分辨率 Radon 变换从 L2 到 L1 的进阶技巧前面讲的都是基于 L2 范数的最小二乘它的缺点是 Radon 域的能量聚焦不够好一次波和多次波的分界模糊。如果你对去多次波的效果要求更高可以试试高分辨率 Radon 变换——用 L1 范数或者混合范数代替 L2。核心思路是把求解问题从 min ||Lm - d||₂² λ||m||₂² 改成 min ||Lm - d||₂² λ||m||₁。L1 正则化倾向于产生稀疏解也就是说 Radon 域的能量会更集中一次波和多次波的分界会更清晰。Matlab 里没有现成的 L1 正则化稀疏求解器但可以用迭代重加权最小二乘IRLS来近似。思路是每次迭代时根据当前解的幅度调整权重幅度大的地方权重小幅度小的地方权重大逼着解往稀疏方向走。function m solve_radon_l1(L, d, lambda, max_outer, max_inner) % SOLVE_RADON_L1 用 IRLS 求解 L1 正则化 Radon 反变换 % 输入: % L - 正变换算子 % d - 数据向量 % lambda - 正则化参数 % max_outer - IRLS 外迭代次数 % max_inner - 内层 CG 迭代次数 % 输出: % m - Radon 域模型 m zeros(size(L, 2), 1); epsilon 1e-6; % 防止除零 for outer 1:max_outer % 根据当前解计算权重 W diag(1 ./ sqrt(abs(m).^2 epsilon)); % 加权最小二乘min ||Lm - d||^2 lambda * m * W * m A L * L lambda * W; b L * d; % 内层 CG 求解 m cg_solve(A, b, m, max_inner); % 输出当前稀疏度 sparsity sum(abs(m) 0.01 * max(abs(m))) / length(m); fprintf(外迭代 %d: 稀疏度 %.3f\n, outer, sparsity); end end function x cg_solve(A, b, x0, max_iter) % 标准共轭梯度求解 x x0; r b - A * x; p r; rsold r * r; for iter 1:max_iter Ap A * p; alpha rsold / (p * Ap); x x alpha * p; r r - alpha * Ap; rsnew r * r; if sqrt(rsnew) 1e-8 break; end p r (rsnew / rsold) * p; rsold rsnew; end end这段代码的关键在权重矩阵 W 的更新。每次外迭代时W 的对角元素是 1/|m_i|幅度大的 m_i 权重小幅度小的 m_i 权重大。这样在最小化过程中小幅度成分会被进一步压制最终解趋向稀疏。参数说明lambda控制稀疏程度越大越稀疏但可能丢失弱能量max_outer一般 5 到 10 次就收敛了epsilon是数值稳定项取 1e-6 左右。实际用的时候我会先用 L2 求解跑一遍看 Radon 域的能量分布确定大致的 p 范围。然后用 L1 求解再跑一遍对比两者的聚焦效果。如果 L1 的 Radon 域能量团明显更集中切除后的道集上多次波压制更干净且一次波没损伤那就用 L1 的结果。如果 L1 把弱一次波也压掉了说明 lambda 太大往回退。还有一个实用技巧不要对整个 Radon 域做 L1 正则化只对多次波所在的 p 区间做。一次波区域的解保持 L2 求解的结果这样既提高了多次波区域的聚焦度又不会损伤一次波。实现上就是构建一个对角权重矩阵在多次波区域加大正则化权重在一次波区域用很小的权重。最后说一个我自己的习惯每次调完参数一定把去多次波前后的道集和叠加剖面并排显示反复对比。不要只看叠加剖面道集上的细节更重要——多次波在近偏移距的残余、一次波在远偏移距的振幅变化这些在叠加剖面上看不出来。我见过太多人只看叠加剖面觉得效果不错结果道集上一看近偏移距的多次波根本没压干净。这个行业里没有后悔药参数调错了重新跑一遍就是几个小时所以每一步都留好中间结果方便回退对比。希望帮到你。本文还有配套的精品资源点击获取
