简介基于Radon正反变换的多次波去除Matlab源码包专为地震勘探数据处理人员与地球物理专业学生设计能有效压制地震记录中呈抛物线形状的多次波干扰。资源共8个文件其中6个Matlab脚本涵盖频率域Radon正变换、反变换、去多次波主程序与演示案例另2个SU格式地震数据文件提供合成记录便于直接运行测试整个RAR压缩包仅218KB轻量高效适合快速学习与复现。已有755人学习过该资源表明该实现具有较好的参考价值。通过阅读脚本注释并运行demo读者能深入理解Radon变换的正反演原理学习将地震数据变换到τ-p域进行分析的方法再通过逆变换恢复有效信号从而掌握去除多次波的标准流程同时可将代码迁移至实际工区根据数据特点调整参数是一种实用性很强的处理方案。1. 为什么压制多次波要先做一次radon正变换而不是直接去速度滤波在 CMP 道集里一次波和多次波都沿双曲线轨迹走差别只在等效速度更慢、曲率更大。很多处理同学第一反应是拉速度谱、上速度滤波结果要么把一次波削平要么远偏移距残留一大截。radon 正反变换的思路是换一个域把道集里的每条同相轴按它自己的剩余曲率“聚焦”到截距时间–曲率平面上的一个点多次波和一次波自然分成两团然后再用反变换把要减掉的那一团重建回时间–偏移距域从原始数据里减掉。这套方法在 matlab 里几十行代码就能搭出原型特别适合处理员、研究多次波压制算法的研究生以及想在正式处理流程前快速评估效果的地球物理工程师。2. radon变换原理讲解从倾斜叠加开始到抛物radon的选型依据2.1 线性radon沿着斜直线求和把同相轴变成一个点radon 变换在地震资料里第一次出名是因为 τ–p 变换也叫线性 radon 或倾斜叠加。原理可以用一句话说清楚取一个斜率 p单位 s/m在整张道集上沿着“截距时间 τ p × 偏移距 x”这条直线对振幅求和。如果道集里恰有一组同相轴沿这条直线走求和结果会明显变强反之不同相位的振幅互相抵消。把所有可能的 p 和 τ 扫一遍就得到一张 radon 谱。从数学上看线性 radon 正变换是m(τ, p) Σ_x d(t τ p x, x)这里的 p 本质上是慢度也就是走时对偏移距的变化率。每个 p 值对应一条直线斜率这个变换相当于把二维数据 (t, x) 重排到 (τ, p) 域。注意反射波在 CMP 道集里是双曲线不同偏移距处的局部斜率是变化的所以一条双曲线同相轴不会聚焦成一个点而会拉成一条能量带。这是线性 radon 对反射资料不太友好的地方但它计算简单、物理含义清楚适合理解后续抛物 radon 的来历。2.2 抛物radon为什么NMO后的多次波残留恰好是抛物线CMP 道集里反射走时满足双曲线关系t² t₀² x² / v²当动校正速度取一次波速度时一次波被拉平多次波等效速度低校正量不足剩余时差在远偏移距展开。把双曲线在 t₀ 附近做泰勒展开t ≈ t₀ x² / (2 v² t₀)忽略高阶项后剩余时差正比于 x²于是抛物 radon 的正变换写成m(τ, q) Σ_x d(t τ q x², x)这里的 q 是曲率参数量纲为 s/m²。对于双曲线同相轴展开后 q ≈ 1 / (2 v² t₀)。一次波在 NMO 后被拉平剩余时差接近零q 趋近于 0多次波等效速度 v_multi 更小q_multi 明显大于 0所以在 (τ, q) 域里一次波和多次波落在不同的能量区一次波压在 q≈0 附近多次波偏向大 q 一侧。这就是抛物 radon 适合去多次波的根本原因。相比线性 radon它把双曲线残差压缩到一个参数 q 上聚焦效果更好相比双曲 radon它不需要在速度维和零偏移距时间维做双参数扫描计算量小一个量级数值稳定性也更好。生产流程里最常见的就是对 NMO 后的道集做抛物 radon。2.3 双曲与线性radon的取舍视角双曲 radon 沿 t sqrt(τ² x²/v²) 求和参数是速度 v它不依赖 NMO对中远偏移距的走时描述更精确适合大偏移距、大剩余时差的资料。代价是速度采样点数和时间点数相乘反演矩阵规模爆炸而且低速多次波与高速深层一次波在速度谱上常常混叠要加很强的阻尼才压得住。线性 radon 则更适合平面波类问题比如 VSP 上行波分离、面波压制用在 CMP 反射道集上会因为双曲线不聚焦而损失分辨率。所以我的选型习惯很简单常规 CMP 道集去多次波先做 NMO再上抛物 radon远偏移距残差明显偏离抛物线时先做各向异性校正或偏移距加权而不是直接换双曲 radon。2.4 用matlab先画一张radon谱判断可分离性在写正式流程之前先在 matlab 里画一张抛物 radon 谱确认一次波和多次波是否真的分得开。% 假设 D_nmo 是已做NMO的道集维度 nt x nxt 是时间轴x 是偏移距 dt t(2) - t(1); q 0 : 0.1e-7 : 8e-7; % 曲率扫描范围单位 s/m^2 tau 0 : dt : t(end); % 截距时间扫描范围与道集对齐 R zeros(numel(tau), numel(q)); for iq 1:numel(q) for ix 1:numel(x) % 每个偏移距沿抛物线平移后累加 t_slice tau q(iq) * x(ix)^2; R(:, iq) R(:, iq) interp1(t, D_nmo(:, ix), t_slice, linear, 0); end end figure; imagesc(q, tau, R); axis xy; colorbar; xlabel(曲率 q (s/m^2)); ylabel(截距时间 tau (s)); title(抛物 radon 谱);这段代码的核心是interp1做逐道时移把每条抛物路径上的振幅累加到一个 (τ, q) 点上。参数上q 的范围建议先按理论值估算对一次波速度 2000 m/s、t₀ 为 0.6 s 的同相轴q ≈ 1/(2×2000²×0.6) ≈ 2.1e-7 s/m²多次波等效速度 1000 m/s 时 q ≈ 8.3e-7所以扫到 8e-7 已经能把两边都覆盖。q 步长取 0.1e-7在最大偏移距 1000 m 下对应时差步长约 0.1 ms足以分辨常规 20–30 Hz 子波。如果谱上两条能量团在 q 方向分得开后面切除边界就好画得多。3. 基于radon正反变换的多次波去除NMO、切除与相减的流程设计3.1 流程全貌从原始CMP道集到干净一次波基于 radon 正反变换的去多次波整个流程可以拆成九个步骤抽取一个 CMP 道集或者沿测线滑动取一组 CMP。做速度分析得到一次波速度场。用一次波速度做 NMO 校正把一次波拉平。对校正后的道集做抛物 radon 正变换得到 (τ, q) 域的 radon 谱。在 radon 谱上确定多次波聚焦带通常是大 q 一侧。设计切除边界把多次波带或者一次波带置零。对保留部分做反变换重建多次波模型或一次波模型。把多次波模型做振幅、相位匹配后从原始 NMO 道集中减去。对输出结果做反 NMO恢复到正常时间轴。第 8 步是最容易翻车的环节。radon 正反变换重建出来的模型只保证了运动学位置正确子波振幅和相位受插值、切除边界影响已经偏离原始数据。直接拿模型去减浅层容易出现“减过头”深层则可能残留同相轴尾巴。常见做法是用最小二乘匹配一个标量或短算子把模型匹配到数据之后再减。3.2 策略一切除多次波带重建模型再相减这是生产流程里最常用的做法。正变换得到 radon 谱后把 q 大于某个阈值的区域保留q 小于阈值的区域置零反变换输出的就是多次波模型。记这个模型为 D_multi原始 NMO 道集为 D_nmo去多次波结果近似为D_demulti D_nmo - α · D_multi其中 α 是最小二乘匹配系数。这个策略的好处是最终输出仍然是原始道集的形态振幅随偏移距的变化保留得相对完整。浅层也不容易出问题直接重建一次波在 τ 很小的地方截断反变换会形成振铃把浅层一次波整个吃掉而先建多次波模型再从原数据减浅层能量的损失要小得多。缺点同样明显多次波和一次波在 radon 域并非完全不重叠。当两者截距时间接近、曲率差异不大时q 方向上会有重叠区切除边界无论怎么画都会误伤一侧。所以切除边界一般不是一条直线而是按 (τ, q) 画一条曲线边界外做渐变衰减而不是硬性置零。3.3 策略二直接在radon域重建一次波另一种做法是切除多次波带后对剩余的一次波区域直接反变换得到一次波道集作为输出。思路直观代码量也小适合快速验证分离效果。但风险很直接第一radon 反变换对切除边界非常敏感边界两侧的截断会在 (t, x) 域产生旁瓣振铃一次波波形被改造第二NMO 校正和反 NMO 两次插值叠加子波频谱会发生畸变后续做 AVO、谱比分析时会被误判。所以我一般在代码原型阶段用它看个趋势真正要交出去的处理结果还是走策略一。策略做法优点风险适用场景策略一切除多次波带重建模型再相减保幅好、浅层稳定重叠区漏删或误删生产流程首选策略二保留一次波带直接反变换输出代码少、直观截断振铃、子波畸变原型验证、分离测试4. matlab骨架实现radon正反变换与多次波去除的代码与参数4.1 合成一个含多次波的CMP道集并做NMO写代码之前先造一组能控制住答案的数据。下面的 matlab 脚本生成 0–1000 m 偏移距、道间距 50 m 的 CMP 道集一次波速度 2000 m/st₀ 0.6 s多次波等效速度 1000 m/st₀ 1.1 s。% 合成含一次波和多次波的 CMP 道集 dt 0.002; % 采样间隔 2ms nt 1001; % 时间采样点数总时长 2s t (0:nt-1) * dt; x 0 : 50 : 1000; % 偏移距 0~1000m间距 50m nx numel(x); fm 30; % 30Hz 雷克子波 tw -0.05 : dt : 0.05; w (1 - 2*(pi*fm*tw).^2) .* exp(-(pi*fm*tw).^2); v1 2000; t01 0.6; % 一次波 v2 1000; t02 1.1; % 多次波等效速度更低 D zeros(nt, nx); for ix 1:nx t1 sqrt(t01.^2 (x(ix)/v1).^2); t2 sqrt(t02.^2 (x(ix)/v2).^2); D(:, ix) interp1(tw t1, w, t, linear, 0); D(:, ix) D(:, ix) 0.8 * interp1(tw t2, w, t, linear, 0); end % 用一次波速度做 NMO 校正 D_nmo zeros(nt, nx); for ix 1:nx t_nmo sqrt(t.^2 - (x(ix)/v1).^2); valid isreal(t_nmo) (t_nmo 0); t_nmo(~valid) 0; D_nmo(:, ix) interp1(t, D(:, ix), t_nmo, linear, 0); D_nmo(~valid, ix) 0; end逻辑说明interp1的作用是把子波移动到双曲线走时位置interp1(twt1, w, t, linear, 0)是一套标准做法t 轴上的每个采样点都能从以 t1 为中心的子波上取值。NMO 部分t_nmo sqrt(t² - (x/v1)²)把原始时间映射回零偏移距时间远偏移距浅层因为 t x/v1 会得到复数或负数先用valid标记再置零避免interp1越界。4.2 抛物radon正变换函数这是一个教学版实现三层循环逻辑最直观也最容易和公式对应。function R radon_p_forward(D_nmo, t, x, q, tau) % 抛物 radon 正变换投影式 % 输入 % D_nmo : nt x nx已做 NMO 的道集 % t : 时间轴 % x : 偏移距轴 % q : 曲率扫描向量 % tau : 截距时间扫描向量 % 输出 % R : radon 域谱维度 numel(tau) x numel(q) nt numel(t); ntau numel(tau); nq numel(q); R zeros(ntau, nq); for it 1:ntau for iq 1:nq s 0; for ix 1:numel(x) tt tau(it) q(iq) * x(ix)^2; if tt t(1) || tt t(end) continue; end s s interp1(t, D_nmo(:, ix), tt, linear, 0); end R(it, iq) s; end end end逻辑说明内层循环对每个偏移距计算目标时间 tt τ q·x²线性插值取振幅并累加这就是沿着抛物线路径的离散叠加。外层两个循环分别扫 τ 和 q。参数上tau 的范围要和道集时间对齐q 的范围按 2.4 节的理论公式估算步长用最大偏移距换算的时差增量来控制让q_step * x_max²明显小于子波主周期的一半。这个三层循环在道集规模大时非常慢生产上一般把这个算子写成稀疏矩阵 L用共轭梯度解最小二乘问题。教学版保留三重循环是为了让你能一步步断点调试看每个样本是怎么累加进去的。4.3 反变换与多次波模型重建反变换采用伴随反投影式把 radon 域每个点的值沿同一条抛物线路径分配回 (t, x) 平面。function D_model radon_p_inverse(R, t, x, q, tau) % 抛物 radon 反变换伴随反投影式 % 输入 R 为 radon 域谱输出为重建的模型道集 nt numel(t); nx numel(x); ntau numel(tau); nq numel(q); D_model zeros(nt, nx); for it 1:ntau for iq 1:nq val R(it, iq); if val 0, continue; end for ix 1:nx tt tau(it) q(iq) * x(ix)^2; if tt t(1) || tt t(end) continue; end idx round((tt - t(1)) / (t(2) - t(1))) 1; D_model(idx, ix) D_model(idx, ix) val; end end end end逻辑说明正变换是“取出来累加”反变换是“放回去累加”两者互为伴随关系不需要解矩阵。val 0的提前跳过在切除之后能省掉大量无效循环。这里用round做最近邻分配是为了让代码短、运行快真要做精细保幅应该用线性插值把val按比例分给相邻两个时间点。接着是完整的切除、重建、相减主脚本% 设置扫描参数并正变换 q 0 : 0.1e-7 : 8e-7; tau (0 : dt : t(end)); R radon_p_forward(D_nmo, t, x, q, tau); % 切除只保留多次波带q 大于阈值 q_cut 3.0e-7; R_multi R; R_multi(:, q q_cut) 0; % 反变换得到多次波模型 D_multi radon_p_inverse(R_multi, t, x, q, tau); % 最小二乘振幅匹配后相减 a sum(sum(D_nmo .* D_multi)) / sum(sum(D_multi .^ 2)); D_demulti D_nmo - a * D_multi; % 反 NMO 回原始时间轴 D_out zeros(nt, nx); for ix 1:nx t_orig sqrt(t.^2 (x(ix)/v1).^2); D_out(:, ix) interp1(t, D_demulti(:, ix), t_orig, linear, 0); end参数说明里最关键是q_cut。它不能盲目取固定值要看 2.4 节那张 radon 谱上两条能量带的分界点。本组合成数据里一次波在 q≈2.1e-7 附近多次波在 q≈8e-7 附近取 3e-7 能把多次波带基本框住。当现场数据存在随机噪声时radon 谱上多次波能量会扩散到更低 qq_cut 往往要降到理论值以下以牺牲少量一次波为代价换取多次波残余不冒头。振幅匹配系数a是标量版本实际生产会用多道最小二乘匹配滤波器我这里用标量是因为演示数据里多次波与一次波正交性较好标量已经够用。4.4 radon参数设置速查表参数作用经验值备注q 扫描范围决定能分离多大剩余时差按 q≈1/(2v²t₀) 估算多次波速度取一次波速度的 0.5–0.7 倍估算q 步长决定曲率分辨率保证 q_step · x_max² 主周期一半0.1e-7 在 1000 m 偏移距下约对应 0.1 mstau 范围覆盖目标层位时间区间与道集时间对齐浅层建议从 0 开始留半个子波长度阻尼 λ控制最小二乘反演的稳定性主对角线均值的 0.01–0.1 倍仅在矩阵解法中用教学版不需要q_cut切除边界先看 radon 谱再定边界处用渐变锥形比硬切除更稳匹配算子模型到数据的振幅/相位修正标量或 3–5 点 FIR 滤波器标量适合高信噪比低信噪比用多道匹配5. 避坑radon正反变换去多次波的5个常见问题与排查5.1 现象radon谱出现“微笑”假象切除边界难定radon 谱上能量团不是圆点而是带弧度的弧形斑常被大家叫“微笑”假象。原因是正变换的插值精度不够线性插值在子波主频较高时会产生频散同相轴越倾斜能量越往下弯。解决先把采样率降低试算如果“微笑”变淡说明是插值问题。把interp1的linear换成spline或者用 sinc 插值能明显改善聚焦。另一个来源是偏移距道间距太大空间假频混进 radon 域这时只能靠道内插或者缩小 q 扫描范围来缓解。5.2 现象反变换后浅层一次波能量被“吃掉”切除多次波带之后反变换输出的一次波模型在浅层往往比原始数据弱很多浅层同相轴像被洗掉一样。原因浅层一次波对应的 τ 很小正变换路径大部分落在时间轴的一小段里radon 域能够利用的采样点数很少反变换叠加回去时能量天生不足。这属于抛物 radon 对浅层覆盖次数的限制不是参数没调好。解决把处理的起始时间从 0 改成目标层顶部以上 50–100 ms给浅层留出足够的 radon 域冗余或者对浅层单独开一个更宽的 tau 扫描窗口再做分窗处理。不要在浅层图和深层共用同一套 q、tau 参数。5.3 现象空道、坏道、偏移距分布不均导致条带噪声实际道集里经常有空道、死道和废道。直接把这些道放进正变换循环radon 谱上会出现垂直于 q 轴的条带噪声反变换后变成横贯道集的斜条纹。原因三维循环里的continue只是跳过该道但道缺失造成的采样不均匀在 radon 域表现为高波数假频。这是空间欠采样问题不是公式问题。解决正变换前先把空道置零但记录索引在正变换循环里跳过它们最稳妥的办法是把所有道先做一次规则化插值使偏移距重新排布成等间距。坏道如果集中在远偏移距宁可整段切除也不要让它们进反演否则影响范围是全道集。5.4 现象未做NMO直接上抛物radon一次波和多次波叠在一起有的人图省事直接对原始未校正的 CMP 道集做抛物 radon结果 radon 谱上所有能量挤成一条连续带根本看不出两个可分离的团块。原因抛物 radon 的假设是剩余时差呈抛物线而原始双曲线的时间随偏移距变化并不满足 t τ qx² 的形式。不做 NMO双曲线的各个偏移距成分会映射到不同 q能量被拉散。解决严格按 3.1 的流程先做 NMO 再正变换。如果因为速度不准导致 NMO 后一次波仍有残余时差宁可把残余时差压到最小再上 radon。速度误差超过 5% 时radon 域聚焦质量会明显下降先去重拾速度而不是调 radon 参数。5.5 现象q_cut边界太硬反变换后出现振铃切除边界从有效值直接跳到零反变换输出的模型在时间轴上会出现高频振铃像在被切除区域周围长了一圈旁瓣。原因radon 域的硬截断等价于在曲率方向乘了一个矩形窗对应到时间域是 sinc 函数的旁瓣。这是傅里叶变换性质决定的绕不开。解决在 q_cut 两侧各设一个过渡带例如边界前后各 0.3e-7 的范围内用余弦渐减加权而不是直接乘 0/1。过渡带的宽度按子波主周期的三分之一到一半来试。还有一个技巧对 radon 域做一次平滑后再反变换振铃也会减弱代价是曲率分辨率略降。6. 用合成数据验证你的radon去多次波流程残差与保幅性的三个检查流程跑通后别急着上真实资料先用合成数据把结果量化。我习惯在调用去多次波函数之后紧跟三个检查。第一个检查是残差道集。把去多次波前后的 radon 谱放在一起对比看大 q 区域的能量是否真的被压下去。如果多次波区域还有明显亮点说明 q_cut 取大了或者匹配算子的相位数不够需要把边界内移。第二个检查是保幅性。把去多次波结果和干净的一次波合成数据做差残差里如果只剩子波拖尾和插值噪声说明保幅过关如果残差里有一次波的波形阴影说明切除边界已经侵入了一次波能量。这里我会计算一次波窗口内的相关系数目标要大于 0.95低于 0.9 基本不能用于后续 AVO 分析。第三个检查是 AVO 响应。沿某个水平同相轴提取振幅随偏移距的变化曲线和理论曲线比对。radon 去多次波最容易犯的错是远偏移距振幅被整体压低因为多次波模型在那里的能量最强匹配相减时会拉着一次波一起减。发现远偏移距振幅掉头向下时优先减小匹配算子长度而不是调 q_cut。进阶方向上我后来会把单次全局 radon 换成时窗分窗处理把道集按 200–400 ms 的时间窗切开每个窗口独立扫描 tau 和 q浅层用小 q 范围深层用大 q 范围。这样能缓解 5.2 节说的浅层能量覆盖不足问题代价是边界处的拼接要留 50% 重叠区再用线性渐变融合。矩阵解法也是值得升级的方向把正变换算子写成稀疏矩阵后用 LSQR 求解谱的聚焦度比教学版投影法高不少处理大偏移距数据时尤其明显。你把这套代码用 matlab 的profile跑一遍会发现三层循环里interp1是最大瓶颈改成内置interp1q或者先做矩阵重排再插值速度能快出两三个量级。希望帮到你。本文还有配套的精品资源点击获取