基于高阶累积量与模拟退火的地震子波相位恢复方法

简介：面向地震资料数字处理中的子波提取需求，这份资源提供基于模拟退火的高阶累积量子波提取方法的全套MATLAB源代码。算法将高阶累积量作为目标函数，利用模拟退火策略在解空间内进行全局寻优，避开局部极值，适用于地震子波估计、信号反演等科研与工程场景；同时模块划分明确，既适合新手快速上手，也方便有经验者针对特定资料进行二次开发。压缩包为rar格式，共六个文件，均为m脚本或函数文件，整体仅四KB左右，代码量轻量紧凑，便于快速查阅。内容涵盖扰动生成、模拟退火迭代、地震记录合成、误差计算等核心模块，按功能拆分为独立文件，便于按步骤阅读、修改与调试；通过调整扰动幅度与退火参数，可直观对比不同设定下的子波提取结果。源码已经过测试校正，确认可正常运行，目前已有237人学习下载，可作为地震勘探、地球物理专业算法研究与课程设计的实用参考。

1. 高阶累积量子波提取：把相位从地震道里抠出来

地震子波估计的难点从来不在振幅谱，而在相位谱。自相关函数把相位信息丢得干干净净，传统反褶积只能靠最小相位假设硬补，可陆上采集、井旁道和做过地表一致性处理的记录里混合相位子波非常普遍，假设一旦失守，后续的波阻抗反演、薄层厚度解释全部跟着偏。用高阶累积量做子波提取，思路上是绕开二阶统计量，直接从四阶累积量里同时重建振幅谱和相位谱；而模拟退火负责求解这个天然多峰的非凸目标函数。这套 MATLAB 工程包含从合成记录生成、加噪、四阶累积量估计、模拟退火反演到误差评估的完整链路，适合做反褶积、井震标定、子波整形，或者给深度学习模型制作训练集的人——只要手里有一段反射界面清晰的叠后道，流程就能直接跑起来。

2. 高阶累积量为什么能把相位找回来

2.1 欠定反演与相位丢失的根源

地震记录用离散褶积模型描述：

x(t) = w(t) * r(t) + n(t)

观测到的只有 x，子波 w 和反射系数 r 都未知，这是个本质欠定问题。维纳滤波和预测反褶积的做法是人为指定 w 为最小相位，把欠定问题变成定解问题，代价是相位谱被强行绑定到振幅谱上。一旦真实子波不是最小相位，反褶积输出会残留明显的旁瓣和时移。自相关函数是二阶矩，功率谱只含幅度信息，相位谱在计算过程中被积掉了。因此凡是只依赖功率谱的方法，理论上都无法恢复混合相位子波。高阶累积量统计的是波形在不同时延位置上的联合分布，相位结构被保留在这些高阶矩里，这为混合相位子波估计提供了信息基础。

2.2 四阶累积量切片的定义与高斯抑制性

零均值平稳随机过程 x(t) 的四阶累积量定义为：

C4x(τ1, τ2, τ3) = E{x(t)x(t+τ1)x(t+τ2)x(t+τ3)} − E{x(t)x(t+τ1)}·E{x(t+τ2)x(t+τ3)} − E{x(t)x(t+τ2)}·E{x(t+τ1)x(t+τ3)} − E{x(t)x(t+τ3)}·E{x(t+τ1)x(t+τ2)}

高斯过程的概率密度完全由一阶和二阶矩决定，因此任意高于二阶的累积量恒等于零。把这个性质放到子波提取场景里：地震噪声通常被建模为高斯或近似高斯过程，它的四阶累积量贡献为零；而反射系数序列是稀疏尖峰式的超高斯分布，四阶累积量携带大量有效信息。这意味着目标函数天然具备对高斯噪声的结构性免疫，不需要预先估计噪声方差再做加权。这个特性和信噪比描述的是两回事——后者依赖能量比例，前者依赖分布类型，这正是高阶统计量方法在低信噪比下仍然可用的原因。

2.3 目标函数：相关系数型的累积量拟合

在反射系数 r 为独立同分布超高斯白噪声的假设下，线性系统高阶统计量输入输出关系给出了一个非常紧凑的公式：

C4x(τ1, τ2, τ3) = γ4 · Σ_t w(t)·w(t−τ1)·w(t−τ2)·w(t−τ3)

其中 γ4 = E[r⁴] − 3(E[r²])² 是反射系数的四阶累积量。取 τ3=0 的二维切片，即可将待估计子波直接代入上式生成理论累积量 C4w，再与观测记录估计出的 C4x 比较。最直接的方式是最小二乘：

J = ‖C4x − λ·C4w‖²

但这个形式需要额外估计尺度因子 λ，因为褶积模型里子波振幅和反射系数振幅存在固有的尺度模糊，搜索过程会在 λ 方向上来回震荡。更稳的做法是相关系数型目标函数：

J(w) = 1 − ρ², ρ = (C4xᵀ·C4w) / (‖C4x‖·‖C4w‖)

尺度 λ 被相关系数天然归一化掉，目标函数只关心波形的形状匹配。项目里的 wave4st.m 对应的就是四阶累积量估计模块，核心实现如下：

function C4 = cum4slice(x, maxlag) % 固定 tau3=0 的四阶累积量二维切片 % 输入 x 为地震记录列向量，maxlag 为最大滞后 x = x(:) - mean(x); % 零均值化 N = length(x); M = maxlag; C4 = zeros(2*M+1, 2*M+1); % 索引覆盖 -M..M for tau1 = -M:M for tau2 = -M:M % 有效区间：保证 t, t+tau1, t+tau2 都在 1..N 内 lo = max([1, 1-tau1, 1-tau2]); hi = min([N, N-tau1, N-tau2]); if hi >= lo + 4 % 样本数过少时跳过 C4(tau1+M+1, tau2+M+1) = ... mean( x(lo:hi) .* x(lo+tau1:hi+tau1) ... .* x(lo+tau2:hi+tau2) .* x(lo:hi) ); end end end end

这段代码用时间平均代替统计期望，是实际工程中最常见的累积量估计方式。循环边界 lo/hi 的作用是防止索引越界，同时保证四个序列片段长度一致；hi >= lo + 4是经验性的最小样本数保护，样本太少时该滞后点的累积量估计方差会非常大，直接置零反而更有利于后续反演稳定。

有了观测数据的 C4x，目标函数计算如下：

function J = objective(w, C4target) % 相关系数型目标函数 M = (size(C4target,1)-1)/2; L = length(w); C4w = zeros(size(C4target)); for tau1 = -M:M for tau2 = -M:M lo = max([1, 1-tau1, 1-tau2]); hi = min([L, L-tau1, L-tau2]); if hi >= lo C4w(tau1+M+1, tau2+M+1) = ... sum(w(lo:hi) .* w(lo+tau1:hi+tau1) ... .* w(lo+tau2:hi+tau2) .* w(lo:hi)); end end end a = C4target(:); b = C4w(:); rho = (a' * b) / (norm(a) * norm(b) + eps); J = 1 - rho^2; end

理论累积量 C4w 中省略了 γ4 和振幅尺度，因为它们只影响整体缩放，而相关系数对尺度不敏感。这个目标函数是子波采样值的四次多项式，天然多峰，且存在周期性时移模糊和尺度模糊，梯度下降法很容易卡在局部极值——这正是引入模拟退火的原因。

3. 模拟退火：逃逸累积量多峰曲面的实现

3.1 Metropolis 准则与温度控制

模拟退火的核心是 Metropolis 接受准则：

P(accept) = 1, 当 ΔJ ≤ 0 P(accept) = exp(−ΔJ / T), 当 ΔJ > 0

温度 T 高时，能量上升的劣化解也有较高概率被接受，搜索可以在目标函数曲面上大范围跳跃；温度逐步降低后，接受劣解的概率收缩，算法逐渐退化为局部精细搜索。温度下降通常采用指数退火策略：

T(k) = α · T(k−1)

α 取 0.9~0.98，温度从 T0 指数衰减到终止温度 Tmin。α 越接近 1，退火越慢，搜索越充分，但计算量也越大。

3.2 主循环：完整反演函数实现

项目中的 monituihuo.m 是模拟退火主程序，核心框架如下：

function [wbest, Jbest, trace] = sa_inversion(C4target, params) % 模拟退火反演地震子波 w = randn(params.wLen, 1); % 随机初值 w = w - mean(w); % 去除直流分量，子波不允许有 DC J = objective(w, C4target); wbest = w; Jbest = J; T = params.T0; wstd = std(w); step = params.sigma0 * wstd; % 初始扰动步长 trace = zeros(300, 2); % 温度与接受率记录 k = 0; while T > params.Tmin acc = 0; % 本轮接受次数 for i = 1:params.L % 马尔可夫链长度 wNew = perturb(w, step, T); JNew = objective(wNew, C4target); dJ = JNew - J; if dJ < 0 || rand < exp(-dJ / T) w = wNew; J = JNew; acc = acc + 1; if J < Jbest Jbest = J; wbest = w; end end end k = k + 1; trace(k, :) = [T, acc / params.L]; T = params.alpha * T; step = params.sigma0 * sqrt(T / params.T0); % 扰动随温度收缩 end trace = trace(1:k, :); end

每次扰动后计算一次目标函数，接受后更新当前解；最优解单独保存，防止后期温度过低时错过更好的状态。step 随 √(T/T0) 收缩，使搜索从粗扫逐渐过渡到精细调整，这是比固定步长更稳定的做法。

3.3 扰动算子的三个选择

扰动算子的设计直接影响搜索效率。我一般会做三种扰动的组合：

function wNew = perturb(w, step, T) % 交替使用 Cauchy 重尾扰动与高斯局部扰动 if rand < 0.5 % Cauchy 扰动：重尾分布，偶发大步长跳跃 delta = step * tan(pi * (rand(size(w)) - 0.5)); else % 高斯扰动：局部小步长精细调整 delta = step * randn(size(w)); end wNew = w + delta; wNew = wNew - mean(wNew); % 防止直流漂移 end

Cauchy 分布的重尾特性使算法有概率产生大步长跳跃，配合温度控制可以更快逃离局部极小；高斯扰动则负责在收敛阶段做精细搜索。两者交替使用，比单一扰动策略覆盖更全面的搜索尺度。

3.4 参数表与接受率监控

运行后观察 trace 中温度与接受率的变化趋势：

figure; yyaxis left; plot(trace(:,1)); ylabel('温度 T'); yyaxis right; plot(trace(:,2)); ylabel('接受率');

正常情况接受率应从 0.7 左右逐步降到 0.1 以下。如果某轮接受率长期高于 0.9，说明温度降得太慢或扰动步长太小；如果一开始就低于 0.1，说明 T0 过低或 sigma0 过大，算法根本没有进行有效搜索。

4. 合成记录全链路验证与抗噪表现

4.1 构造带真解的混合相位子波

验证模拟退火反演效果，必须构造一个带真解的测试场景。用零相位 Ricker 子波串联全通滤波器，生成混合相位子波：

N = 2000; dt = 0.001; f0 = 30; tc = (-0.1:dt:0.1)'; w0 = (1 - 2*(pi*f0*tc).^2) .* exp(-(pi*f0*tc).^2); % Ricker a1 = 0.55; a2 = 0.35; % 全通滤波器参数 w1 = filter([a1 1], [1 a1], w0); % 一阶全通 w2 = filter([a2 1], [1 a2], w1); % 二级级联 w_true = w2 / norm(w2); % 归一化便于误差比较

一阶全通滤波器 H(z) = (a + z⁻¹)/(1 + a·z⁻¹) 的幅频响应恒为 1，只改变相位谱，因此 w_true 和 w0 的振幅谱完全一致，相位谱不再是最小相位。关键是这里保留了真解，可以直接定量评估反演误差。

4.2 合成反射系数、褶积与加噪

r = zeros(N, 1); % 稀疏反射系数 p = randperm(N, round(N * 0.06)); % 6% 非零 r(p) = randn(length(p), 1); r(p) = r(p) / std(r); x = conv(w_true, r, 'same'); % 褶积 snr = 15; sigma_n = std(x) * 10^(-snr/20); xn = x + sigma_n * randn(N, 1); % 加高斯噪声

反射系数序列只有约 6% 的非零样本，分布呈现明显的超高斯特征——这是高阶累积量方法能起作用的前提条件。加噪使用固定 SNR 方式，噪声方差由信号标准差换算得到。

4.3 反演命令与收敛过程解读

整个流程对应项目文件的分工：sawave.m/sawavee.m 负责生成子波，synthesis.m 完成合成记录，Disturbance.m 完成加噪，wave4st.m 完成四阶累积量估计，monituihuo.m 完成模拟退火反演，wucha.m 完成误差评估。实验调用顺序如下：

M = 20; % 累积量滞后范围 C4target = cum4slice(xn, M); params = struct('T0', 1.5, 'alpha', 0.96, 'Tmin', 1e-3, ... 'L', 120, 'wLen', 21, 'sigma0', 0.12); [w_est, J_best, trace] = sa_inversion(C4target, params);

子波长度取 21 个样点，覆盖主频 30 Hz 时约 60 ms 的时间窗；累积量滞后 M=20 与子波长度匹配，保证目标函数有足够的约束信息。

4.4 抗噪性能对比表

反演结束后用 wucha.m 中的归一化均方误差和相关系数评估质量，对齐后再计算相位残差。在子波长度 25、记录长度 2000、反射系数稀疏度 6% 的条件下，单次运行的典型结果如下：

SNR 降到 10 dB 以下时，四阶累积量的估计方差迅速增大，目标函数曲面上出现大量虚假极值，反演结果的旁瓣明显增多。这进一步说明：高阶累积量对高斯噪声有结构性免疫，但免疫的前提是累积量估计本身要足够准。

4.5 失效边界与常见坑

反射系数接近高斯分布时，四阶累积量趋近于零，目标函数失去有效梯度信息，反演结果基本不可用。记录长度太短（小于 500 个样点）时，累积量估计方差过大，解决办法是分段估计后取平均。子波长度取得过长，反演参数维度增加，模拟退火搜索效率大幅下降，建议先做短窗反演确定主波形，再逐步加长细化。

5. 多链校验、反褶积自检与初值策略

5.1 时移对齐与多链稳定性校验

模拟退火是单链随机过程，单次运行存在漏掉最优盆地的风险。常见做法是跑 4 条独立链，比较各链终解的 J 值和波形。比较之前必须先做时移对齐，否则累积量目标函数的平移模糊会让 NMSE 虚高：

function weAlign = alignWavelet(we, wt) % 以 wt 为参考对 we 做时移对齐 [c, lags] = xcorr(we, wt, 'coeff'); [~, i] = max(c); sh = lags(i); if sh >= 0 weAlign = [zeros(sh,1); we]; else weAlign = we(-sh+1:end); end weAlign = weAlign(1:max(length(we), length(wt))); % 截齐 end

对齐后用 xcorr 的最大相关系数作为波形相似度的判据。如果多条链的 J 值接近但对齐后波形差异明显，说明目标函数存在严重的时移或尺度模糊，需要对累积量切片做能量归一化后再跑一轮。

5.2 反褶积自检

恢复子波是否可靠，最终要看它能否把地震记录反褶积成尖脉冲。频域除法加规则化是最常用的自检方式：

wPad = zeros(size(xn)); wPad(1:length(w_est)) = w_est; Xf = fft(xn); Wf = fft(wPad); reg = 0.05 * max(abs(Wf)); % 规则化系数 S = Xf .* conj(Wf) ./ (abs(Wf).^2 + reg^2); spike = ifft(S);

reg的作用是抑制子波谱零点附近的噪声放大，取值范围通常取最大谱值的 0.01~0.1 倍。检查 spike 序列的主瓣宽度和旁瓣电平：旁瓣电平超过主瓣 30% 时，说明反演子波的相位谱仍然存在偏差。

5.3 混合初值与降采样加速

初值选择对模拟退火的收敛速度和最终精度影响显著。我一般会准备四组初值：随机高斯序列、截断 sinc 函数、倒序的随机序列、常数序列，分别代表不同相位特征，各自跑一遍后取 J 值最小的解。这样即使某条链掉进局部极小，其他链仍有概率找到更优盆地。提速方面，可以先把记录降采样到原采样率的 1/4，在低分辨率数据上快速锁定子波大致形状，再用全分辨率数据做第二轮精化，两步合计耗时通常只有直接全分辨率反演的一半。

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