简介这份资源面向地震勘探数据处理的学习者与研究人员聚焦基于Radon正反变换去除抛物线型多次波这一经典难题。包内共8个文件以6个m脚本和2个su地震数据文件为主压缩包约218KBm文件承担正反Radon变换、频率域实现与segy读取等核心运算su文件则提供合成单炮记录及含多次波记录用于对比验证。资源从Radon变换的积分定义与投影原理讲起逐步过渡到抛物线型多次波的识别、投影值滤波置零以及反Radon变换的数值恢复并配有可直接运行的Matlab演示脚本便于读者复现完整处理流程。目前已有755人学习下载适合希望理解Radon变换原理并动手实践多次波压制的地震资料处理入门者与进阶者参考。1. Radon 变换去多次波为什么它至今仍是地震资料处理的第一道防线做过海上地震资料处理的人都有一个共识多次波不去掉后面的速度分析、偏移成像全是白费功夫。而 Radon 变换去多次波几乎是每个处理流程里绕不开的一步。它的核心逻辑不复杂——把时空域t-x的地震道集变换到 Radon 域τ-p 或 τ-q利用一次波和多次波在 Radon 域中不同的曲率分布把它们分离开再反变换回时空域达到压制多次波的目的。听起来简单但实际调参时动校正拉伸、截断效应、假频这些坑一个比一个深。这篇笔记从 Radon 正反变换的数学原理讲起用 MATLAB 把线性 Radon 和高分辨率 Radon 的完整流程跑通最后给出我在实际资料上踩过的坑和参数建议。适合正在做多次波压制、需要快速复现 Radon 流程的地震处理从业者。2. Radon 变换的数学骨架从连续积分到离散矩阵2.1 线性 Radon 变换的正变换公式推导线性 Radon 变换也称倾斜叠加的正变换定义为沿直线积分u(τ, p) ∫ d(t τ p·x, x) dx其中d(t, x)是时空域地震数据u(τ, p)是 Radon 域结果p是射线参数斜率τ是截距时间。物理含义很直观对每一个斜率p沿着该斜率对应的直线把振幅累加累加结果就是 Radon 域中该(τ, p)点的值。反变换则是把 Radon 域的每个点按对应斜率撒回时空域d(t, x) ∫ u(τ t - p·x, p) dp连续形式看着简洁但实际数据是离散采样的必须转成矩阵形式才能编程实现。把时空域数据拉成列向量d长度Nt × NxRadon 域拉成列向量m长度Nτ × Np正变换写成d L · m其中L是正变换算子矩阵元素L[i,j]表示第j个 Radon 域点对第i个时空域点的贡献权重。反变换就是求解这个线性方程组。2.2 离散化时的三个关键参数离散化不是随便切一刀就行有三个参数直接决定变换质量参数含义经验取值影响Npp 方向采样点数81~161太少则斜率分辨率不够多次波和一次波在 Radon 域重叠dpp 采样间隔由pmax/Np决定太大会产生空间假频太小则计算量暴增pmax最大斜率由最大偏移距和最深目标层速度估算必须覆盖多次波的最大斜率否则截断pmax的估算公式pmax xmax / (vmin · sqrt(xmax² 4·zmax²))其中xmax是最大偏移距zmax是目标层最大深度vmin是最小层速度。这个公式我一般会算完之后再乘 1.2 倍留余量。2.3 为什么反变换不能直接求逆理论上m L⁻¹ · d就能反变换但L通常是高度病态的——条件数动辄 10⁶ 以上。直接求逆的结果就是噪声被无限放大Radon 域里稍微有一点误差反变换回来就是满屏雪花。所以实际做法是求最小二乘解并加正则化项m (LᵀL λI)⁻¹ · Lᵀ · dλ是阻尼因子控制解的平滑程度。λ太大则分辨率下降一次波和多次波分不开λ太小则数值不稳定。我一般从λ 0.01 × trace(LᵀL)/M开始试M是 Radon 域点数。3. MATLAB 实现从正变换到最小二乘反变换3.1 构建 Radon 正变换算子矩阵function L build_radon_operator(Nt, Nx, dt, dx, Np, pmax) % 构建线性 Radon 正变换算子矩阵 % 输入 % Nt - 时间采样点数 % Nx - 道数 % dt - 时间采样间隔 (s) % dx - 道间距 (m) % Np - p 方向采样点数 % pmax - 最大斜率 (s/m) % 输出 % L - 正变换矩阵大小 (Nt*Nx) x (Ntau*Np) Ntau Nt; % 截距时间采样点数一般与 Nt 相同 dp 2 * pmax / (Np - 1); p_axis linspace(-pmax, pmax, Np); tau_axis (0:Ntau-1) * dt; L zeros(Nt * Nx, Ntau * Np); for ip 1:Np p p_axis(ip); for ix 1:Nx x (ix - 1) * dx; for itau 1:Ntau tau tau_axis(itau); t tau p * x; % 沿直线的时间 it round(t / dt) 1; % 最近邻插值到时间网格 if it 1 it Nt row (ix - 1) * Nt it; col (ip - 1) * Ntau itau; L(row, col) 1; end end end end end这段代码的核心逻辑是对每一个(τ, p)组合计算它对应到每一道的到达时间t τ p·x然后把该时间采样点的权重设为 1。这里用的是最近邻插值简单但精度有限。实际生产中我会换成线性插值或 sinc 插值代价是矩阵更稠密、计算更慢。参数说明Np建议取 121 或 161pmax按 2.2 节公式估算后乘 1.2。Ntau一般等于Nt但如果数据时间长度很大可以适当截断因为远偏移距的多次波截距时间不会超过最大记录时间。3.2 阻尼最小二乘反变换function m radon_inverse(L, d, lambda) % 阻尼最小二乘求解 Radon 域模型 % 输入 % L - 正变换矩阵 % d - 时空域数据列向量 % lambda - 阻尼因子 % 输出 % m - Radon 域模型列向量 [M, N] size(L); % 构造法方程矩阵 A L * L lambda * eye(N); b L * d; % 求解 m A \ b; end这里直接用 MATLAB 的左除求解法方程。对于中等规模数据Nt1000, Nx120, Np121矩阵大小约(120000 × 121000)内存占用约 116 GB普通工作站扛不住。实际做法有两种一是分时间窗口处理每次只处理 200~300 个时间采样二是用共轭梯度法迭代求解不显式存储LᵀL。lambda的选取我一般用试错法先取lambda 0.001 * trace(L*L)/N看 Radon 域聚焦效果如果一次波和多次波还是混在一起就减小如果出现明显噪声条带就增大。3.3 完整流程脚本与合成数据验证% 生成合成数据一个一次波 一个多次波 Nt 512; Nx 96; dt 0.004; dx 25; t (0:Nt-1) * dt; x (0:Nx-1) * dx; d zeros(Nt, Nx); % 一次波斜率 0.0002 s/m截距时间 0.3 s for ix 1:Nx t0 0.3 0.0002 * x(ix); it round(t0/dt) 1; if it Nt d(it, ix) 1; end end % 多次波斜率 0.0004 s/m截距时间 0.5 s for ix 1:Nx t0 0.5 0.0004 * x(ix); it round(t0/dt) 1; if it Nt d(it, ix) d(it, ix) 0.7; end end % 构建算子并反变换 Np 121; pmax 0.0008; L build_radon_operator(Nt, Nx, dt, dx, Np, pmax); d_vec d(:); lambda 0.01; m_vec radon_inverse(L, d_vec, lambda); m reshape(m_vec, Nt, Np); % 显示 Radon 域结果 figure; imagesc(linspace(-pmax, pmax, Np), t, m); xlabel(p (s/m)); ylabel(tau (s)); title(Radon 域结果); colorbar;跑完这段代码你会在 Radon 域看到两个明显的能量团一个在p0.0002附近一个在p0.0004附近。接下来在 Radon 域里把p较大的区域多次波置零或衰减再反变换回时空域就完成了多次波压制。反变换回时空域直接用d_filtered L * m_filtered其中m_filtered是切除后的 Radon 域模型。注意这里不需要再求逆因为正变换算子L是已知的。4. 高分辨率 Radon 变换为什么你的多次波总是去不干净4.1 最小二乘 Radon 的分辨率瓶颈上一节的最小二乘解有一个致命问题它在 Radon 域产生严重的能量拖尾。一个理想的一次波在 Radon 域应该是一个点但最小二乘解会把它抹成一个椭圆形的能量团多次波的能量团也变大两者一重叠切除时就分不干净。这就是很多人反映“多次波去不干净一次波还伤了”的根本原因。分辨率瓶颈来自两个方面一是有限孔径效应地震数据的偏移距范围有限Radon 变换的积分区间不完整二是最小二乘解本身倾向于能量最小化而不是聚焦最大化。4.2 迭代重加权与稀疏约束高分辨率 Radon 变换的核心思想是引入稀疏约束让 Radon 域的能量尽可能集中。常用做法是迭代重加权最小二乘IRLSfunction m radon_highres(L, d, lambda, niter) % 迭代重加权高分辨率 Radon 变换 % 输入 % L - 正变换矩阵 % d - 时空域数据列向量 % lambda - 初始阻尼因子 % niter - 迭代次数 % 输出 % m - Radon 域模型 [~, N] size(L); m zeros(N, 1); epsilon 1e-4; % 防止除零 for iter 1:niter % 构建加权矩阵 W diag(1 ./ (abs(m) epsilon)); % 加权法方程 A L * L lambda * W; b L * d; m A \ b; end end迭代次数一般取 5~10 次。第一次迭代时m初始化为零向量W是单位矩阵的缩放等价于普通最小二乘。后续迭代中m中能量大的位置权重变小能量小的位置权重变大迫使解向稀疏方向收敛。lambda在这里的作用和之前不同它控制稀疏性的强弱。lambda越大解越稀疏Radon 域聚焦越好但可能丢失弱能量lambda越小解越接近最小二乘分辨率下降。我一般从lambda 0.1开始试根据 Radon 域聚焦效果调整。4.3 切除策略在 Radon 域怎么切才不伤一次波切除是 Radon 去多次波的最后一步也是最容易翻车的一步。常见做法有三种硬切除设定一个p阈值大于阈值的区域全部置零。简单粗暴但阈值附近的一次波和多次波会同时被切掉或同时保留。自适应切除对每个τ时间点根据该时间点一次波和多次波的能量比动态确定切除边界。实现时可以用滑动窗口计算能量比再设定一个比例阈值。减去法不直接切除而是把 Radon 域中多次波区域反变换回时空域得到多次波模型再从原始数据中减去。这样做的好处是避免切除边界处的能量突变减少假象。我一般用减去法虽然多一步反变换但对一次波的保护更好。具体做法是在 Radon 域中把一次波区域置零保留多次波区域反变换得到多次波模型d_mult然后d_primary d - d_mult。5. 避坑指南Radon 去多次波最常见的五个翻车现场5.1 假频导致 Radon 域出现虚假能量团现象Radon 域中出现多个对称的能量团切除后时空域数据出现规律性噪声。原因dp采样间隔太大不满足 Nyquist 条件。Radon 变换的假频条件与时空域采样有关dp必须满足dp 1 / (2 · fmax · xmax)其中fmax是数据最高频率xmax是最大偏移距。解决减小dp即增大Np。如果计算量不允许可以先对数据做低通滤波降低fmax。5.2 动校正拉伸导致远偏移距波形畸变现象远偏移距道在 Radon 域中能量发散切除后远偏移距多次波残留严重。原因动校正拉伸使远偏移距子波变宽Radon 变换的直线假设不再成立。解决在 Radon 变换前先做部分动校正PMO把双曲线拉平到接近直线再做 Radon 变换。或者采用抛物线 Radon 变换它对动校正后的数据适应性更好。5.3 截距时间范围不够导致多次波能量折叠现象Radon 域中多次波能量出现在τ边界处反变换后多次波压制不彻底。原因Ntau取值太小远偏移距多次波的截距时间超出了 Radon 域的τ范围能量被折叠到边界。解决Ntau至少取Nt pmax · xmax / dt确保所有可能的截距时间都被覆盖。5.4 阻尼因子选取不当导致一次波损伤现象多次波去掉了但一次波同相轴也变弱了分辨率下降。原因lambda太大Radon 域解过于平滑一次波和多次波的能量团重叠切除时误伤一次波。解决减小lambda或改用高分辨率 Radon 变换。判断标准是看 Radon 域中一次波和多次波的能量团是否分离清晰。5.5 分窗口处理时窗口重叠不足导致边界假象现象分时间窗口处理后窗口边界处出现水平条带假象。原因窗口之间没有重叠或者重叠区域没有做加权过渡。解决窗口重叠至少 50%重叠区域用余弦加权或线性加权过渡。我一般用 50% 重叠加汉宁窗加权。6. 从合成到实际资料Radon 去多次波的验证方法与调参习惯合成数据跑通只是第一步实际资料上的验证才是关键。我一般按这个顺序做第一步用合成数据验证算子正确性。构建一个已知的一次波和多次波模型跑完 Radon 正反变换后检查反变换结果与原始数据的残差。如果残差能量大于原始数据能量的 1%说明算子构建或求解有问题。第二步用实际资料的近偏移距道集做参数扫描。近偏移距道集多次波能量弱一次波清晰适合用来判断 Radon 变换是否损伤了一次波。扫描Np、pmax、lambda三个参数每个参数取 3~5 个值看 Radon 域聚焦效果和反变换后一次波同相轴的连续性。第三步用远偏移距道集验证多次波压制效果。远偏移距道集多次波能量强是检验切除策略的关键。重点看多次波同相轴是否被压制到噪声水平同时一次波是否保持连续。第四步用速度谱验证。多次波压制前后各做一次速度分析如果多次波能量在速度谱上明显减弱说明压制有效。注意不要过度压制否则速度谱上会出现一次波能量也被削弱的情况。调参时我有一个习惯每次只改一个参数改完立刻看 Radon 域和时空域的双重效果。同时改多个参数出了问题根本不知道是哪个参数导致的。另外lambda的调整幅度不要超过 2 倍否则解的性质会突变。还有一个容易被忽略的点Radon 变换去多次波的效果与数据本身的信噪比密切相关。如果原始数据信噪比很低Radon 域中一次波和多次波的能量团都会被噪声淹没这时候应该先做去噪再做 Radon 变换。我见过不少人跳过去噪直接上 Radon结果多次波没去掉反而把噪声变成了规律性假象后悔药都没得吃。最后说一个验证技巧把压制掉的多次波模型单独反变换回时空域显示出来。如果多次波模型里混入了明显的一次波同相轴说明切除策略有问题需要调整。这个步骤相当于给 Radon 去多次波装了一个黑匣子能直观看到哪些能量被去掉了哪些被保留了。希望帮到你。本文还有配套的精品资源点击获取
