写这篇东西的起因是上周给一批刚入职的同事讲合成记录标定又有人在看完30Hz雷克子波后问我这个负瓣是坏道吗能不能直接把它砍掉这个问题我几乎每年都要回答一遍。Ricker小波也就是雷克子波大概是地震勘探里最常用的零相位理论子波它的数学形式看起来极简频谱也干净但真拿到实际资料里大家操心的往往不是主瓣而是主瓣两侧那对负旁瓣。旁瓣不是噪声不是采集脚印也不是多次波——它写在Ricker小波的数学表达式里是子波本身的一部分。搞清楚它的幅度多大、峰值出现在什么时间、能量占比多少、能不能压掉是做好薄层分辨和反演的基础功课。下面我用自己的推导习惯和工程经验把Ricker小波旁瓣的数学表达与计算过程完整过一遍。1. 为什么要把旁瓣单独拎出来算合成记录里最熟悉的假峰1.1 一个30Hz子波的直观印象标准的时间域Ricker子波表达式长这样[ w(t)\left[1-2(\pi f_m t)^2\right]e^{-(\pi f_m t)^2} ]其中 (f_m) 是主频也就是幅度谱的峰值频率。我最早接触这个公式时觉得它人畜无害直到把它画出来才意识到问题主瓣确实又窄又尖但主瓣两侧各挂着一个负极性的小山包。取 (f_m30) Hz旁瓣峰值时间大约在 (t\pm 0.013) s也就是主瓣两侧各13毫秒的位置相对幅度接近主瓣的44.6%。这个44.6%是个什么概念道间互相关、层位追踪、自动增益的振幅统计都会受到它的干扰。最典型的场景是合成记录一个30Hz的Ricker子波和一个反射系数序列做褶积后旁瓣会形成一系列额外的同相轴它们与真实界面没有一一对应关系却长得和弱反射一模一样。你如果没有事先算过旁瓣位置很容易把这些假峰当小层解释进去。1.2 旁瓣和分辨率的三角关系地震分辨率的经典判据是Ricker提出的可分辨距离两个等幅反射同相轴当时间间距小于子波主瓣宽度时合成波形的峰值位置会偏移、幅度会重组。这个判据有一个隐含前提——我们眼睛盯着的都是主瓣。但实际处理中两个反射相距更远时旁瓣之间的干涉会先发生。举一个我踩过坑的例子模拟一个顶底反射时间差为 (t_s) 的薄层用主频30Hz的Ricker子波做正演。当薄层时间厚度等于旁瓣峰位置13ms量级时顶反射的主瓣和底反射的旁瓣叠在一起波形看起来会出现一个中等强度的中间峰新手很容易把它解释成第三个薄层。这就是为什么我坚持任何涉及薄层厚度的分析开工前都要先把所用子波的旁瓣峰值时间、幅度和极性打印出来贴在屏幕上。否则后面所有解释都是在跟一个自己没看清的子波属性较劲。2. 从时间域表达式手推旁瓣峰值、位置、过零点一个不落2.1 先做变量代换把公式里的常数请出去直接对带 (f_m) 的公式求导虽然也能算但表达式会显得很啰嗦。我习惯先做一个无量纲代换[ x\pi f_m t ]这样Ricker子波就变成[ R(x)(1-2x^2)e^{-x^2} ]这个代换的物理意义很直接(x) 相当于用主频周期把时间轴归一化。所有关于旁瓣相对幅度、相对位置的结论在这个坐标下都与主频无关真正要换算成秒再乘回 (1/(\pi f_m))。这样做的好处是数学推导干净而且所有结果天然具备尺度不变性不会因为子波主频不同就改变相对结论。需要提醒一句如果你在文献里看到 (w(s)\left[1-2(\pi f_m s)^2\right] \exp(-\pi^2 f_m^2 s^2))这个是同一种子波只是把指数里的平方展开了别被不同写法绕晕。归一化方式也分好几种有的按峰值幅度归一有的按能量归一有的按频谱峰值归一。本文统一采用时间域峰值归一的写法这也是地震软件里最常用的约定。2.2 求导找驻点旁瓣峰值和出现位置直接解析解出要求旁瓣峰值我对 (R(x)) 求一阶导[ R(x)(-4x)e^{-x^2}(1-2x^2)(-2x)e^{-x^2} ]整理一下[ R(x)(4x^3-6x)e^{-x^2}2x(2x^2-3)e^{-x^2} ]令 (R(x)0)得到三个驻点[ x_10,\quad x_{2,3}\pm \sqrt{\frac{3}{2}} ]其中 (x_10) 是主瓣极大值另外两个就是旁瓣极值。把 (x^23/2) 代回[ R\left(\pm\sqrt{\frac{3}{2}}\right)\left(1-2\cdot\frac{3}{2}\right)e^{-3/2}-2e^{-3/2} ]数值是[ -2e^{-3/2}\approx -0.44626 ]这就是旁瓣相对幅度的解析结果约为正极性主瓣峰值的44.6%。换算到真实时间轴[ t_s\pm \frac{\sqrt{3/2}}{\pi f_m}\approx \pm \frac{0.3898}{f_m} ]我用 (f_m30) Hz算一下(t_s\approx \pm 0.013) s用 (f_m50) Hz则只有 (±0.0078) s。这个公式我在实际标定时几乎每周都用强烈建议直接背下来。2.3 过零点与主瓣宽度比峰值更容易算错的地方旁瓣的位置解决了还有一个工程上经常要用的量过零点和主瓣宽度。令 (R(x)0)由于指数项不可能为零只能让括号为零[ 1-2x^20 \Rightarrow x\pm \frac{1}{\sqrt{2}}\approx \pm 0.7071 ]换算成时间[ t_0\pm \frac{0.7071}{\pi f_m}\approx \pm \frac{0.2251}{f_m} ]主瓣宽度就是两个过零点之间的时间长度[ T_\text{main}\frac{2\times 0.7071}{\pi f_m}\approx \frac{0.4502}{f_m} ]30Hz子波的主瓣宽度约为15ms50Hz约为9ms。我见过不少人把旁瓣峰值位置和过零点搞混要么把 (t_0) 当旁瓣位置要么把主瓣宽度算成 (1/f_m)。注意 (1/f_m) 是主频对应的周期和主瓣宽度不是一回事。Ricker子波的主瓣宽度大约是主频周期的0.45倍旁瓣峰位置大约在主频周期的0.39倍处这两个数相差很小但物理含义完全不同用的时候要分清楚。2.4 旁瓣极性的物理解释很多人都注意到Ricker子波的旁瓣是负极性的这不是偶然。观察 (R(x)(1-2x^2)e^{-x^2})指数项恒大于零所以整个子波的极性完全由括号里的 (1-2x^2) 决定当 (|x|1/\sqrt{2})括号为正子波为正极性这是主瓣所在区间当 (|x|1/\sqrt{2})括号为负子波立刻翻转为负极性形成旁瓣。换句话说Ricker子波只有一对负旁瓣旁瓣与主瓣相位差180度。这带来一个很实用的判据在做层位追踪时如果看到紧邻强反射之后、时间距离大约等于旁瓣峰值时间处出现一个极性反转的同相轴先别急着解释成岩性变化大概率是子波旁瓣。另外Ricker子波在 (x) 很大时虽然负值但以指数速度衰减并不会像一些滤波器的响应那样产生振荡的多次旁瓣。很多文献说Ricker子波是最小旁瓣的零相位子波严格说应该理解为它只有单一边瓣而不是旁瓣幅度很小。3. 旁瓣不是幅度问题而是能量问题3.1 频谱视角Ricker子波为什么会有旁瓣Ricker子波的幅度谱可以用傅里叶变换得到[ |W(f)| \propto \left(\frac{f}{f_m}\right)^2 e^{-(f/f_m)^2} ]这个谱的特征是低频端按 (f^2) 快速抬升高频端按高斯函数快速衰减峰值正好落在 (ff_m)。在时域里一个带限信号必然伴随时间上的旁瓣振荡——这是傅里叶变换的基本性质你不可能既把频带压得干干净净又让时间波形只保留一个孤零零的主瓣。Ricker子波的旁瓣本质上就是它的高斯型频谱在有限带宽约束下产生的振铃。我用半幅度点估算过它的有效带宽大致在 (0.48f_m) 到 (1.64f_m) 之间总带宽约 (1.16f_m)。这个不对称的频谱形状决定了时域子波左右对称但旁瓣形态相对简单。3.2 旁瓣能量占比的积分算式与数值结果很多人以为旁瓣峰值44.6%能量占比也差不多是40%量级实际一算会吓一跳。用无量纲变量 (x\pi f_m t)子波能量可以写成[ E\int_{-\infty}^{\infty} (1-2x^2)^2 e^{-2x^2},dx ]展开后逐项积分[ E\int_{-\infty}^{\infty} e^{-2x^2}dx -4\int_{-\infty}^{\infty} x^2 e^{-2x^2}dx 4\int_{-\infty}^{\infty} x^4 e^{-2x^2}dx ][ E\sqrt{\frac{\pi}{2}} -4\cdot\frac{\sqrt{\pi}}{4\sqrt{2}} 4\cdot\frac{3\sqrt{\pi}}{16\sqrt{2}} ]前两项正好抵消最后得到[ E\frac{3\sqrt{\pi}}{4\sqrt{2}}\approx 0.9399 ]主瓣区间取 ([-1/\sqrt{2}, 1/\sqrt{2}]) 的积分正常子波只有一对旁瓣所以旁瓣能量占比为[ \eta_\text{side}1-\frac{E_\text{main}}{E_\text{total}}\approx 0.296 ]也就是说旁瓣峰值只有主瓣的44.6%但旁瓣能量占了整个子波能量的近30%。这个峰值低、能量不低的现象很容易被忽视。在反演或者振幅分析里如果你用的是加权降旁瓣算法却只盯着波形峰值看就可能漏掉旁瓣能量的真实影响。3.3 一个常被忽略的事实提高主频不会压低相对旁瓣有一次做高分辨率处理项目有人提出主频提高以后旁瓣自然就小了。这句话对了一半。提高主频确实让旁瓣在时间轴上更靠近主瓣、绝对时间更短但由于我们用的是峰值归一的Ricker子波(x\pi f_m t) 代换之后波形形状完全不变。换句话说不管 (f_m) 是20Hz还是80Hz归一化Ricker子波的旁瓣高度永远是44.6%旁瓣能量占比永远是29.6%。真正改变旁瓣相对幅度的因素是什么是子波频谱的形状而不是位置。如果你需要更低旁瓣的零相位子波就得改谱形比如用更接近矩形谱的Ormsby子波、或者对Ricker频谱做整形。这个结论在后续章节会详细讲。记住一点单靠提高主频来让旁瓣变小在相对幅度意义上是做不到的。4. 旁瓣参数的工程计算解析公式和Python脚本对照着来4.1 用SciPy验证解析结果解析公式虽然已经给出但工程上免不了要批量计算和验证。我一般用Python把结果复核一遍顺便生成子波图形。下面这段代码直接可用里面用的是无量纲自变量 (x)这样算出来的结果和主频无关import numpy as np from scipy.signal import find_peaks from scipy.integrate import quad def ricker_z(x): 峰值归一化的Ricker子波x pi * fm * t return (1.0 - 2.0 * x**2) * np.exp(-x**2) x np.linspace(-5.0, 5.0, 20001) w ricker_z(x) # 找负旁瓣峰值 peaks, _ find_peaks(-w, distance100) print(旁瓣峰值位置 x , x[peaks]) print(旁瓣峰值幅度 , w[peaks]) # 过零点位置 x_zero 1.0 / np.sqrt(2.0) print(过零点位置 x , -x_zero, x_zero) # 能量 E_total, _ quad(lambda s: ricker_z(s)**2, -np.inf, np.inf) E_main, _ quad(lambda s: ricker_z(s)**2, -x_zero, x_zero) print(总能量 , E_total) print(主瓣能量 , E_main) print(旁瓣能量占比 , (E_total - E_main) / E_total)跑出来的结果应该是旁瓣位置 (x\pm1.2247)旁瓣幅度 (-0.4463)过零点 (\pm0.7071)总能量约0.9399旁瓣能量占比约0.296。这些都是无量纲结果。如果你要算某个具体主频下的时间位置直接除一个 (\pi f_m) 就行def ricker_side_time(fm): return np.sqrt(1.5) / (np.pi * fm), -2.0 * np.exp(-1.5) for fm in [10, 20, 30, 50]: t_side, amp ricker_side_time(fm) t_zero 1.0 / (np.sqrt(2.0) * np.pi * fm) print(ffm{fm:3d} Hz: 旁瓣时间 /-{t_side*1000:6.2f} ms, f过零点 /-{t_zero*1000:6.2f} ms, 旁瓣幅度 {amp:.4f})4.2 批量生成不同主频的旁瓣参数表实际项目里我们经常要一口气评估多个主频子波的旁瓣影响。我习惯直接把常用参数打表贴在工位旁边主频 (f_m) (Hz)旁瓣峰时间 (\pm t_s) (ms)过零点 (\pm t_0) (ms)主瓣宽度 (ms)旁瓣相对幅度1038.9822.5145.020.4462019.4911.2522.510.4463012.997.5015.010.446409.745.6311.260.446507.804.509.000.446804.872.815.630.446这张表的价值在于当你在剖面上看到一个强反射旁瓣距离大约等于上表数值时第一反应不是发现新层位而是验证旁瓣位置。比如主频40Hz的资料旁瓣峰大约在10ms处如果某个局部振幅异常正好出现在强轴两侧10ms左右且极性反转那几乎可以断定是旁瓣干涉而不是真实构造。4.3 计算时容易踩的采样和边界坑用程序算Ricker子波参数有几个坑我经常看新人踩第一个坑是采样率不够。离散化时如果时间采样间隔 (dt) 太大旁瓣峰值可能正好落在两个采样点之间导致 (find_peaks) 找到的幅度严重偏小。以30Hz子波为例旁瓣峰在13ms处若采样率只有4ms旁瓣区域还勉强有点若采样率8ms整个旁瓣包络就废了。稳妥做法是保证每个主频周期至少有20个采样点也就是 (dt \le 1/(20 f_m))。第二个坑是截断边界。很多人生成子波时只截 (-1/f_m) 到 (1/f_m)这会把旁瓣后半段砍掉。虽然旁瓣幅度在 (|x|1.5) 以后快速衰减但能量积分对截断位置很敏感。我做能量计算时最低截到 (|x|5)对应30Hz子波约±53ms才敢说积分稳定。第三个坑是归一化方式混用。前面说过峰值归一、能量归一、频谱归一三种写法下旁瓣的绝对数值完全不同。如果你拿别人代码里的能量归一化子波再跟我这里的44.6%峰值对比对不上是正常的。所以在任何文档里只要写出旁瓣幅度必须同时写清归一化定义否则后面所有分析都失真。5. 处理现场怎么跟旁瓣过招削弱、利用以及我的个人习惯5.1 为什么我不建议直接对旁瓣做截窗碰到旁瓣明显最直接的想法是用时间窗把负旁瓣抹掉。这个操作我劝阻过无数次。原因很简单时域截窗等价于频谱卷积你把旁瓣硬砍掉主瓣旁边的频谱会重新生出更复杂的振铃经常是按下葫芦浮起瓢。Ricker子波的旁瓣是频谱带限性的结果不是独立附加的噪声源。正确的思路不是拿走旁瓣而是换一个旁瓣更小的子波或者说对原始数据做整形让等效子波从Ricker朝低旁瓣目标子波靠拢。5.2 几种常见的旁瓣整形思路我试过并且认为有效的方法有三类第一类是频谱整形。既然Ricker子波频谱是 (f^2 e^{-(f/f_m)^2})低频端抬升太慢、高频端衰减太快那就可以通过谱白化、宽频处理把频谱往目标谱形状上靠。目标谱越接近矩形时域旁瓣越低但主瓣也会相应变宽时间分辨率会打折扣。所以处理参数上要在旁瓣最低和主瓣不垮之间折中。第二类是反褶积类方法。最小平方反褶积、预测反褶积都能把子波向尖脉冲压缩但代价是相位关系可能改变零相位子波会变成混合相位旁瓣结构反而变得不对称、更难预测。如果资料信噪比不够高反褶积带来的旁瓣放大可能比Ricker原生的负旁瓣更麻烦。第三类是匹配滤波整形。设计一个目标低旁瓣零相位子波比如适当加宽频带的Ormsby子波然后用匹配滤波把实际子波整形过去。这个办法比较可控适合需要保幅又不想引入太多处理假象的场景。我实际项目里更倾向于第三种因为它允许你在频带和旁瓣之间显式做权衡而不是盲目压。5.3 把旁瓣当刻度尺薄层调谐中的伪峰值识别旁瓣除了被动挨打其实可以主动利用。最实用的一个场景是薄层调谐分析。当地层顶底反射时间差靠近旁瓣峰位置 (t_s) 时顶反射主瓣与底反射旁瓣、底反射主瓣与顶反射旁瓣会形成干涉产生一系列固定时间间隔的假振幅峰值。这个时间间隔恰好可以用前面那张旁瓣参数表预估。我的操作习惯是这样拿到工区主频后先算出 (t_s)再在解释软件里生成一条旁瓣警戒线——以强反射同相轴为基准上下各画两条与旁瓣峰位置对齐的辅助线。凡是落在辅助线附近的短同相轴先不做层位解释而是用谱分解或者正演模拟验证它是不是旁瓣干涉的结果。这套办法帮我在三个项目里避免了把旁瓣当薄层解释的误判。特别是薄层厚度在调谐厚度附近时这种做法尤其重要因为调谐曲线本身就有厚度接近四分之一波长时振幅最大的特征再叠加上旁瓣的相干增强很容易让人误以为发现了高阻抗层。另外旁瓣的反极性也是个免费的信息。如果剖面里一个强轴旁边不远处出现极性反转的弱同相轴且时间距离约等于旁瓣峰位置这通常意味着子波旁瓣而不是真实岩性界面。我写处理成果报告时也会专门加一张所用子波旁瓣特征表告诉解释组哪些同相轴是子波贡献、哪些是真实的。这个细节在验收时经常得到好的反馈。最后再分享一个经验任何新工区我一定会在做首轮叠加前把理论子波旁瓣特征算一遍并把参数表存档。做反演、做薄层分析、做属性提取后面每一步几乎都要回查这张表。很多人把Ricker子波的旁瓣当作不可避免的背景噪声去抱怨其实只要把它的数学表达吃透、把计算流程固化下来旁瓣完全可以从干扰项变成判别工具。你用它的位置和极性去解释资料比单纯想方设法躲开它要来得舒服得多。