做地震资料处理的人十有八九都跟Ricker小波打过交道。它对称、零相位、只有一个主频参数是合成地震记录、正演模拟和反演测试里最常用的子波之一。可很多人用的时候只关心主瓣长什么样真到需要把旁瓣的数学表达写清楚、把旁瓣比算准的时候反而容易卡壳——教科书上通常写一句旁瓣幅度约为主瓣的44%但旁瓣极值到底出现在哪个时间、峰值旁瓣比怎么推出来、能量占比有多少、数值代码怎么写很少有资料一次讲透。这篇文章就从Ricker小波的时域数学表达出发把旁瓣的解析推导、量化指标和数值实现完整走一遍。适合做地震正演、子波分析、合成记录调试、以及想搞清楚Ricker小波旁瓣计算公式的同行参考。1. Ricker小波的核心数学表达1.1 时间域定义与无量纲化Ricker子波的时间域最常用形式是$$ R(t) \left(1 - 2\pi^2 f_0^2 t^2\right) e^{-\pi^2 f_0^2 t^2} $$其中 $f_0$ 是主频单位Hz$t$ 是时间单位s。这个表达式已经做了幅度归一化当 $t0$ 时 $R(0)1$所以主瓣峰值永远是1和主频 $f_0$ 没有关系。括号里的 $\left(1 - 2\pi^2 f_0^2 t^2\right)$ 相当于主瓣的基座指数项 $e^{-\pi^2 f_0^2 t^2}$ 是高斯衰减包络两者相乘就得到了中间鼓起、两侧凹陷再收敛的经典Ricker形态。这里有个特别容易踩的坑$\pi$ 的位置。有的文献把公式写成 $R(t) \left(1 - 2(\pi f_0 t)^2\right) e^{-(\pi f_0 t)^2}$跟上面完全等价如果用角频率 $\omega_0 2\pi f_0$ 表示又会变成 $R(t) \left(1 - \frac{1}{2}\omega_0^2 t^2\right) e^{-\frac{1}{4}\omega_0^2 t^2}$。形式五花八门但展开后都是同一个函数。抄公式的时候如果只抄了括号里、漏了指数里的 $\pi$波形和旁瓣比会差得非常远。为了方便推导我习惯先做一个无量纲变换令 $x \pi f_0 t$于是$$ R(x) (1 - 2x^2) e^{-x^2} $$这个变换非常关键它把主频 $f_0$ 从公式里请了出去剩下一个只跟 $x$ 有关的纯函数。旁瓣位置、旁瓣比、零交叉点本质上都是 $x$ 坐标上的常数最后再用 $t x/(\pi f_0)$ 映射回时间坐标即可。理解了这一步后面所有解析结果都会变得清爽很多。1.2 频率域表达与峰值频率对时域的 $R(t)$ 做傅里叶变换可以得到频域表达$$ R(f) \frac{1}{\sqrt{\pi}} \cdot \frac{f^2}{f_0^3} \cdot e^{-f^2 / f_0^2} $$不同文献的归一化约定不一样可能会差一个常数系数但函数形态是固定的$f^2$ 乘以高斯衰减。这个频谱有几个很重要的性质。首先对 $f$ 求导并令导数为零可以得到频谱峰值出现在 $f f_0$。也就是说$f_0$ 既是时域表达式的参数也是频谱峰值频率这就是主频二字的含义。第二个性质是 $R(0)0$说明Ricker子波是零均值信号不产生直流分量这在去均值、频谱分析的时候非常方便。频域表达式还有一个实用价值如果实际资料里提取的子波频谱形态不符合 $f^2 e^{-f^2/f_0^2}$那它大概率不是纯Ricker子波而是经过了滤波整形或与其它子波褶积。做子波分析时我通常先把频谱形状画出来跟这个理论曲线叠在一起看偏差一目了然。1.3 为什么旁瓣值得单独算旁瓣在Ricker子波里不是缺陷而是有限带宽信号不可避免的时域振荡。任何带限信号都不可能在时间域有限支撑高斯衰减把主瓣约束得越窄两侧的振荡就越明显。旁瓣直接影响地震记录的分辨率解释旁瓣幅度大了会被误判成薄层的反射多个相邻反射的旁瓣互相叠加也会改变薄层调谐振幅的形态。做合成记录时旁瓣造成的假反射尤其讨厌。用生活化的比喻来说Ricker小波就像敲了一口钟主瓣是那一声最响的撞击声旁瓣是之后的余响只是这个余响的方向是反的。单看撞击声很干净可余响没算清楚整个声音的听感就变了。对地震信号而言主瓣决定能不能分辨旁瓣决定会不会看错所以做子波分析和正演模拟时对旁瓣的控制往往比只看主瓣更讲究。2. 旁瓣的解析推导一阶导数定位置二阶导数验极值2.1 无量纲形式下的极值方程用无量纲形式 $R(x) (1 - 2x^2) e^{-x^2}$ 求导过程非常干净$$ \frac{dR}{dx} (4x^3 - 6x) e^{-x^2} 2x(2x^2 - 3) e^{-x^2} $$令一阶导数为零得到三个驻点$x 0$ 和 $x \pm \sqrt{\frac{3}{2}}$。$x0$ 对应主瓣峰值这没什么好说的$x \pm \sqrt{3/2} \approx \pm 1.2247$ 就是两个对称旁瓣的极值位置。需要注意旁瓣的极值是负的所以从函数形态上属于局部极小值而不是极大值。为了确认极值性质再求二阶导$$ \frac{d^2R}{dx^2} (-8x^4 24x^2 - 6) e^{-x^2} $$在 $x0$ 处$R(0) -6 0$是局部极大在 $x \sqrt{3/2}$ 处$R 12 0$确实是局部极小。主瓣和旁瓣的极值性质从数学上分得很清楚这也是为什么旁瓣在波形图上表现为两个对称的负峰。把 $x$ 映射回时间坐标旁瓣极值出现在$$ t_s \pm \frac{\sqrt{3/2}}{\pi f_0} \approx \pm \frac{0.3897}{f_0} $$举例来说$f_0 30\text{Hz}$ 时旁瓣极值大约在 $\pm 13\text{ms}$ 的位置$f_0 60\text{Hz}$ 时大约在 $\pm 6.5\text{ms}$ 的位置。主频越高旁瓣离主瓣越近但不要以为近就等于小相对幅度并不会变。2.2 旁瓣峰值与旁瓣比的解析值把 $x_s \sqrt{3/2}$ 代回原函数$$ R(x_s) \left(1 - 2 \times \frac{3}{2}\right) e^{-3/2} -2e^{-3/2} \approx -0.44626 $$旁瓣峰值幅度是主瓣峰值的 $0.4463$ 倍极性为负。这个 $-2e^{-3/2}$ 是一个干净漂亮的常数跟主频 $f_0$ 完全没有关系。峰值旁瓣比Peak Sidelobe Level RatioPSLR的定义是$$ PSLR 20 \log_{10} \frac{|R(x_s)|}{|R(0)|} 20 \log_{10} \left(2e^{-3/2}\right) $$把数值代进去算一遍$\log_{10}(2e^{-1.5}) \log_{10}2 - 1.5 \times \log_{10}e \approx 0.3010 - 0.6514 -0.3504$再乘以 20 得到 $-7.008\text{dB}$。所以大家常说的旁瓣比约 -7dB严格说其实是 $-7.01\text{dB}$。这个结果同样与主频无关是Ricker子波的自相似特性之一。这里要提醒一句有些资料用功率比表示旁瓣比也就是 $10\log_{10}(\text{幅度比}^2)$算出来和 $20\log_{10}(\text{幅度比})$ 完全相同。真正容易出问题的是把主瓣峰值和旁瓣峰值的符号搞混——旁瓣是负的取幅度比时一定要加绝对值。2.3 零交叉点与主旁瓣分界令 $R(x) 0$解 $\left(1 - 2x^2\right) e^{-x^2} 0$得到 $x \pm 1/\sqrt{2} \approx \pm 0.7071$。时间坐标上$$ t_z \pm \frac{1}{\sqrt{2}\pi f_0} \approx \pm \frac{0.2251}{f_0} $$主瓣宽度也就是两个零交叉点之间的距离为 $\sqrt{2}/(\pi f_0) \approx 0.4502/f_0$。旁瓣极值位置在 $\pm 0.3897/f_0$比零交叉点大约晚了 $0.1646/f_0$ 秒。这个先过零、再旁瓣峰的顺序是判断一个波形是不是纯Ricker的重要特征。把关键特征整理成一张表方便直接用特征量$x$ 坐标时间坐标$t x/\pi f_0$幅度主瓣峰值$0$$0$$1$零交叉点$\pm 1/\sqrt{2}$$\pm 0.2251/f_0$$0$旁瓣极值$\pm \sqrt{3/2}$$\pm 0.3897/f_0$$-0.4463$峰值旁瓣比——$-7.01\text{dB}$这张表在实际调试子波时非常有用。比如你在某套资料里提取了一个30Hz的子波发现旁瓣极值不在13ms附近而是出现在20ms那就说明这个子波不是纯Ricker可能被滤波整形过。3. 旁瓣比与能量占比两个量化指标一起看3.1 数值脚本验证峰值旁瓣比解析公式推完了最好再用数值方法验证一遍防止推导过程出错。这里给出一段可以直接跑的Python代码用密集采样生成Ricker子波再用scipy找局部极小值和解析结果对照。import numpy as np from scipy.signal import argrelextrema def ricker(t, f0): x np.pi * f0 * t return (1 - 2 * x**2) * np.exp(-x**2) f0 30.0 fs 2000.0 # 采样率 2kHz足够密 t np.arange(-0.4, 0.4, 1.0 / fs) w ricker(t, f0) # 找到主瓣峰值 peak_idx np.argmax(w) peak_val w[peak_idx] # 找局部极小值Ricker旁瓣是负瓣 min_idx argrelextrema(w, np.less, order5)[0] # 排除主瓣区域内的数值抖动只保留离主瓣较远的旁瓣极值 side_idx [i for i in min_idx if abs(t[i]) 0.5 / f0] left_idx side_idx[0] right_idx side_idx[-1] for idx in (left_idx, right_idx): print(f旁瓣极值 t{t[idx]*1000:.2f} ms, 幅度{w[idx]:.4f}) print(f数值 PSLR {20*np.log10(abs(w[left_idx]) / abs(peak_val)):.3f} dB) print(f解析 PSLR {20*np.log10(2*np.exp(-1.5)):.3f} dB)实测输出大约是这样旁瓣极值 t-13.00 ms, 幅度-0.4463 旁瓣极值 t 13.00 ms, 幅度-0.4463 数值 PSLR -7.008 dB 解析 PSLR -7.008 dB数值和解析完全吻合。用 $order5$ 是为了避免数值噪声造成伪极值采样率越高越稳。3.2 旁瓣能量占比接近三成不可忽略幅度比之外旁瓣能量占比也是值得关注的一个指标。它的计算方法是把主瓣范围之外$|t| t_z$的能量除以总能量dt 1.0 / fs energy_total np.sum(w**2) * dt mask_side np.abs(t) 1.0 / (np.sqrt(2) * np.pi * f0) energy_side np.sum(w[mask_side]**2) * dt print(f旁瓣能量占比 {energy_side / energy_total * 100:.1f}%)算下来大约是 $29.5%$。这个结果很反直觉旁瓣峰值幅度只有主瓣的 $44.6%$可两个旁瓣持续的时间长累计能量竟然接近总能量的三成。在日常做合成记录时这意味着旁瓣不是可以随手忽略的小尾巴在强反射旁边旁瓣叠加产生的假同相轴幅度完全可能达到可识别水平。为什么主瓣那么突出能量占比却不高因为Ricker子波的主瓣虽然幅度大但持续时间短旁瓣幅度衰减快可拖尾很长平方之后的时间积分不可小觑。实际资料里的子波通常还会被加时窗截断如果截断位置太靠近主瓣能量占比会进一步变化。3.3 主频改变时旁瓣表现保持不变由 $x \pi f_0 t$ 可知主频 $f_0$ 只改变时间尺度不改变函数形态。所以主频变化时旁瓣极值时间 $t_s \pm 0.3897/f_0$随主频升高而线性缩短零交叉点时间 $t_z \pm 0.2251/f_0$同样线性缩短峰值旁瓣比固定为 $-7.01\text{dB}$不随主频变化旁瓣能量占比固定约 $29.5%$不随主频变化。这组性质说明Ricker子波是自相似的。做子波压缩时提高主频只是把整个波形挤到更窄的时间范围旁瓣的相对幅度一点都不会变小。如果目标是压制旁瓣单纯提高主频没有用要考虑整形滤波或者换用旁瓣更小的子波类型。这个认知在子波设计阶段非常重要能帮你少走很多弯路。4. 数值实现从解析公式到可复用代码4.1 Python完整脚本解析与数值一并输出上面的片段比较散我把它们整合成一个完整的、带绘图的分析脚本。日常做子波分析时我一般直接复制这个脚本改主频几秒钟就能看到结果。import numpy as np import matplotlib.pyplot as plt from scipy.signal import argrelextrema def ricker(t, f0): x np.pi * f0 * t return (1 - 2 * x**2) * np.exp(-x**2) def ricker_sidelobe_analysis(f0, fs2000.0, t_range0.4): t np.arange(-t_range, t_range, 1.0 / fs) w ricker(t, f0) peak_idx np.argmax(w) peak_val w[peak_idx] min_idx argrelextrema(w, np.less, order5)[0] side_idx [i for i in min_idx if abs(t[i]) 0.5 / f0] left, right side_idx[0], side_idx[-1] # 解析参考值 t_s np.sqrt(3.0 / 2.0) / (np.pi * f0) a_s -2.0 * np.exp(-1.5) psr_db 20.0 * np.log10(abs(a_s)) print(f主频 f0 {f0} Hz) print(f数值左旁瓣: t{t[left]*1000:.3f} ms, A{w[left]:.5f}) print(f数值右旁瓣: t{t[right]*1000:.3f} ms, A{w[right]:.5f}) print(f解析旁瓣: t{t_s*1000:.3f} ms, A{a_s:.5f}) print(f峰值旁瓣比: {20*np.log10(abs(w[left])/peak_val):.3f} dB (解析 {psr_db:.3f} dB)) # 能量占比 dt 1.0 / fs t_z 1.0 / (np.sqrt(2) * np.pi * f0) energy_total np.sum(w**2) * dt energy_side np.sum(w[np.abs(t) t_z]**2) * dt print(f旁瓣能量占比: {energy_side/energy_total*100:.2f}%) # 绘图 plt.figure(figsize(10, 4)) plt.plot(t * 1000, w, lw1.8) plt.axhline(0, colorgray, lw0.8) plt.axvline(t_s * 1000, colorred, ls--, lw1.0) plt.axvline(-t_s * 1000, colorred, ls--, lw1.0) plt.xlabel(Time (ms)) plt.ylabel(Amplitude) plt.title(fRicker Wavelet f0{f0} Hz) plt.grid(alpha0.3) plt.show() if __name__ __main__: ricker_sidelobe_analysis(30.0)运行后可以同时得到旁瓣极值时间、幅度、峰值旁瓣比、能量占比和波形图。其中解析和数值的输出放在一起对比能直观发现代码有没有低级错误。4.2 MATLAB版最小编程量求旁瓣如果你平时用MATLAB下面这段足够用了。核心思路是用gradient求数值梯度再用diff(sign(...))找局部极小值。f0 30; dt 0.0005; t (-0.4:dt:0.4).; x pi * f0 * t; w (1 - 2 * x.^2) .* exp(-x.^2); [peak, pidx] max(w); dwd gradient(w, dt); cross find(diff(sign(dwd)) 2) 1; % 局部极小值 % 筛掉主瓣附近的点保留旁瓣 side cross(abs(t(cross)) 0.5 / f0); [s1, k1] min(w(side)); s2 max(w(side)); fprintf(左/右旁瓣幅度: %.4f, %.4f\n, s1, s2); fprintf(数值 PSLR: %.3f dB\n, 20*log10(abs(s1)/peak)); fprintf(解析 PSLR: %.3f dB\n, 20*log10(2*exp(-1.5))); tz 1 / (sqrt(2) * pi * f0); E_total sum(w.^2) * dt; E_side sum(w(abs(t) tz).^2) * dt; fprintf(旁瓣能量占比: %.2f%%\n, E_side / E_total * 100); plot(t*1000, w, LineWidth, 1.6); hold on; plot(t(pidx)*1000, peak, ro); plot(t(side(k1))*1000, s1, rv); grid on; xlabel(Time (ms)); ylabel(Amplitude);MATLAB的diff(sign(gradient))在局部极小处会出现从 -1 到 1 的跳变差值为2所以用 2就能定位。如果波形本身有噪声建议先做平滑滤波再求梯度否则会找到一堆伪极值。4.3 计算精度与边界条件控制数值计算Ricker旁瓣最常见的误差来源有三个采样率不足、时间窗截断、噪声引起的伪极值。采样率方面经验法则是采样间隔至少小于主瓣宽度的一半工程上我习惯取 $dt \le 1/(10f_0)$。以30Hz为例$dt \le 3.3\text{ms}$ 就能看到旁瓣但为了极值精度取 $0.5\text{ms}$ 更稳。时间窗方面如果窗长太短旁瓣尾部被截断能量占比会被低估太短甚至可能截到旁瓣峰值本身。建议单边时间窗至少取到 $4/f_0$ 以上让指数项衰减到足够小。以30Hz为例$t$ 从 $-0.4\text{s}$ 到 $0.4\text{s}$ 完全够用尾部幅度已经接近零。数值求导方面梯度方法对采样率和噪声都比较敏感。如果子波是从实际地震记录里提取的不可避免带有噪声最好先对波形做带限滤波或者干脆用解析公式做最小二乘拟合直接从拟合参数里读出旁瓣极值位置和幅度。这比在带噪数据上找极值稳定得多。5. 常见问题与排查技巧实录5.1 为什么我算出来的旁瓣比不是 -7dB这是被问得最多的一个问题。旁瓣比不是 -7dB通常不是公式问题而是数值实现问题。第一采样率太低旁瓣峰值恰好落在两个采样点之间数值上记录的幅度偏小第二时间窗太短旁瓣还没完全展开就被截断了第三找极值时把边界效应或数值噪声当成了旁瓣。排查时按顺序做三件事先把采样率提高到 $f_s \ge 10f_0$ 以上再把时间窗加宽到至少 $\pm 4/f_0$最后把波形尾部打印出来确认它已经衰减到接近零。如果做完这三步数值旁瓣比依然是 -7dB 左右那基本就是对的。5.2 主频的标定方式不同结果差很多主频这个词在不同软件里定义不一样。Ricker公式里的 $f_0$ 是频谱峰值频率也就是频谱上幅度最大的那个频率但有些工具统计的是占优频率频谱重心有些则是视频率。如果你从频谱上量出峰值频率 30Hz却用 $R(t) (1-2\pi^2 f_0^2 t^2)e^{-\pi^2 f_0^2t^2}$ 去生成子波那没问题可如果你用频谱重心频率30Hz去代替公式里的 $f_0$波形就会偏胖或偏瘦旁瓣位置和主瓣宽度全对不上。所以做子波对比时一定要先确认你手里的主频是按哪种方式定义的。最稳妥的方法是用频谱峰值标定 $f_0$再用谱峰位置校验生成出来的子波频谱。5.3 实际提取的子波旁瓣位置对不上怎么办实际资料里提取的子波几乎不可能是纯Ricker。滤波、吸收衰减、叠加过程都会改变子波形态旁瓣可能变得更宽、更矮甚至不对称。判断一个子波偏离纯Ricker多远我推荐一个简单指标把主瓣峰值到旁瓣峰值的时间间隔与主瓣峰值到零交叉点的时间间隔做比值。对于纯Ricker这个比值是$$ \frac{t_s}{t_z} \frac{\sqrt{3/2}}{1/\sqrt{2}} \sqrt{3} \approx 1.732 $$如果你的实际子波这个比值明显偏离1.732比如变成了1.5以下或者2.0以上那它大概率已经不是纯Ricker了。这个比值不随主频变化非常稳定比单纯看波形图更可靠。5.4 旁瓣能不能再压一压Ricker的旁瓣比固定 -7dB能量占比约三成这确实让不少做高分辨率处理的人头疼。如果想要更小的旁瓣可以考虑换用Butterworth子波——它的幅度谱更方时域旁瓣通常更小也可以通过最小平方整形把子波向目标旁瓣特性整形。但要注意压制旁瓣通常以展宽主瓣为代价时间分辨率和旁瓣幅度之间存在取舍没有免费的午餐。我的建议是先搞清楚你的应用场景到底是看薄层还是躲假象再决定要不要牺牲分辨率去压旁瓣。6. 写在最后一点实操习惯我自己的习惯是每设计一套观测系统或做一批合成记录都会先把子波旁瓣的极值时间和幅度写进备注。别小看这一步它能在后面帮你快速分辨强反射旁边的影子反射到底是不是旁瓣引起的。比如某层反射旁边 13ms 处出现一个负极性同相轴而你的激发子波正好是30Hz Ricker那这个同相轴十有八九是旁瓣而不是真实的地层响应。最后再分享一个压箱底的小技巧拿到一个提取子波后先别急着看频谱直接量主瓣峰值到旁瓣峰值的时间间隔和主瓣峰值到零交叉点的时间间隔算一下比值是不是 $\sqrt{3}$。如果是基本可以放心把它当Ricker用如果不是赶紧检查子波提取流程或者考虑用更合适的子波模型。这个技巧花不了半分钟却能在子波分析的一开始就帮你筛掉很多坑。