企业官网建设流程全解析

简介：基于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 nx，t 是时间轴，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 正反变换的去多次波，整个流程可以拆成九个步骤：

1. 抽取一个 CMP 道集，或者沿测线滑动取一组 CMP。
2. 做速度分析，得到一次波速度场。
3. 用一次波速度做 NMO 校正，把一次波拉平。
4. 对校正后的道集做抛物 radon 正变换，得到 (τ, q) 域的 radon 谱。
5. 在 radon 谱上确定多次波聚焦带，通常是大 q 一侧。
6. 设计切除边界，把多次波带或者一次波带置零。
7. 对保留部分做反变换，重建多次波模型或一次波模型。
8. 把多次波模型做振幅、相位匹配后，从原始 NMO 道集中减去。
9. 对输出结果做反 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/s，t₀ = 0.6 s，多次波等效速度 1000 m/s，t₀ = 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(tw+t1, 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 谱上多次波能量会扩散到更低 q，q_cut 往往要降到理论值以下，以牺牲少量一次波为代价换取多次波残余不冒头。振幅匹配系数a是标量版本，实际生产会用多道最小二乘匹配滤波器，我这里用标量是因为演示数据里多次波与一次波正交性较好，标量已经够用。

4.4 radon参数设置速查表

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或者先做矩阵重排再插值，速度能快出两三个量级。希望帮到你。

本文还有配套的精品资源，点击获取