简介频率—波数域频散曲线提取方法及程序设计.pdf 是一份面向地震勘探与工程物探技术人员的地学专业文献聚焦瑞利波数据处理中的 f-k 频散曲线提取技术。文档从瑞利波基本概念入手系统讲解二维傅氏变换将时间-空间域数据转换为频率-波数域能量谱的原理并结合半波长理论阐述速度-深度域频散曲线绘制流程其中对能量谱能量最大点选取、视速度与波长计算等关键步骤也做了说明同时以 Delphi7.0 为开发环境给出程序化实现思路与关键模块设计对开展相关算法研究、程序编写或工程应用具有直接参考价值。资源包内包含 1 个 PDF 文件文件大小约 365KB内容精炼、图文结合便于下载后快速阅读与存档。目前已有 551 人浏览学习适合物探、岩土工程及地震数据处理方向的师生或工程师作为专业指导文献使用。1. 为什么是频率—波数域一张 f-k 谱能省掉半天的排列试验处理主动源面波数据时我一直有个习惯拿到炮集先不急着画时距图直接做一次二维傅里叶变换到频率—波数域。原因很简单时距图里面波和折射波叠在一起肉眼很难分清哪条能量是基阶频散而在 f-k 谱里不同视速度的信号按斜率分开面波频散能量会聚成一条窄窄的亮带。频率—波数域频散曲线提取方法及程序设计做的就是这件事把时间-空间域的炮集变成 f-k 谱再从谱面上自动提取频散曲线供后续反演使用。这套方案适合做 MASW、主动源面波勘探的工程师和研究人员也适合刚接触频散分析、想摆脱手工拾取的人。读完之后你至少能自己写出一套能跑的提取程序并且知道参数该往哪个方向调。2. 从时间-空间波场到 f-k 谱二维傅里叶变换怎么做才不白做2.1 一张地震记录在 f-k 域里长什么样能量条带与速度线的对应地震记录可以看成二维函数 u(x, t)x 是检波点位置t 是时间。对时间做傅里叶变换得到频率 f对空间做傅里叶变换得到波数 k二维变换之后就是 U(k, f)也就是频率—波数域谱。一个沿 x 方向以速度 c 传播的平面简谐波在 f-k 谱上对应一条过原点的直线斜率就是速度。体波、折射波视速度高直线更陡能量靠近波数轴面波视速度低能量靠近频率轴。因为频散波的相速度随频率变化它的能量不会聚成直线而是形成一条弯曲的亮带。频散曲线提取就是从这条亮带的脊线上按频率逐个取相速度。f-k 域最大的价值在于把不同视速度的能量在平面上分开自动拾取比在时间域里追相位转折点可靠得多。2.2 用 Python 把波场变换到 f-k 域预处理、补零与 FFT 的完整代码这里给出一段可以直接落地的实现。Python 程序设计里最常用的就是 numpy 和 scipy二维 FFT 一行就能算完真正影响结果的是变换前的预处理。import numpy as np from scipy import signal def prepare_shot(data, bad_scale10.0): 把原始炮集整理成适合做 2D FFT 的形态。 data: 二维数组形状 (nx, nt)nx 为道数nt 为采样点数。 nx, nt data.shape # 逐道去均值、去趋势避免零频附近直流泄漏 data data - np.mean(data, axis1, keepdimsTrue) data signal.detrend(data, axis1, typelinear) # 坏道判定振幅超过中位数 10 倍的道直接充零 amp np.max(np.abs(data), axis1) bad amp bad_scale * np.median(amp) data[bad, :] 0.0 # 时间维加窗抑制谱泄漏空间维加窗抑制排列两端截断效应 data * np.hanning(nt)[None, :] data * np.hanning(nx)[:, None] return data def fk_amp_spectrum(data, dt, dx, x_pad4, t_pad4): 对预处理后的炮集做二维 FFT返回振幅谱和坐标轴。 dt: 时间采样间隔(s) dx: 道距(m) nx, nt data.shape # 补零到 2 的幂FFT 更快谱面也更平滑 nfx int(2 ** np.ceil(np.log2(nx * x_pad))) nft int(2 ** np.ceil(np.log2(nt * t_pad))) spec np.fft.fft2(data, s(nfx, nft)) spec np.fft.fftshift(spec, axes0) # 只对波数轴做中心化 freq np.fft.fftfreq(nft, dt) # 频率轴单位 Hz k np.fft.fftfreq(nfx, dx) # 波数轴单位 1/m half nfx // 2 return np.abs(spec[half:, :]), freq, k[half:]先说逻辑。prepare_shot 里先做逐道去均值和线性去趋势这一步不能省否则零波数附近会残留一条强能量把低频端曲线往上拉。坏道充零要在加窗之前因为加窗会把坏道异常值扩散到相邻频率。时间维和空间维都用了汉宁窗汉宁窗主瓣略宽但旁瓣低比矩形窗更适合自动拾取。再看 fk_amp_spectrum。fft2 的第一个维度是道数也就是空间维s(nfx, nft) 指定补零后的形状。fftshift 只在波数轴移动把零波数放到中间然后只取后半段对应正波数。这样做的原因是常见单边排列里波只沿一个方向传播负波数半轴基本只有对称的镜像能量。返回的 k 单位是 cycles/m注意这里不是角波数所以后面相速度直接 c f / k不用乘 2π。提示np.fft.fftfreq 生成的是空间频率单位 1/m。如果换成角波数 rad/m相速度换算要写成 c 2πf / k_angular。x_pad 和 t_pad 是补零倍数一般取 4 到 8。补零能加密谱线让峰值拾取更光滑但它不改变真实分辨率真实分辨率只由原始道数和排列长度决定。频率轴里包含负频率后面提取函数会通过 fmin/fmax 把它们排除掉。2.3 f-k 谱分辨率到底被什么锁死三个算式与最小道距判定在做 f-k 分析前先用三个算式估算谱的分辨率和可用范围能避免后面反复调参。第一个是频率分辨率Δf 1 / (nt·dt)也就是记录总时间长度的倒数。1 秒的记录Δf 就是 1 Hz想在 2 Hz 附近分辨两个很近的模式基本没戏。第二个是波数分辨率Δk 1 / (nx·dx)也就是排列长度的倒数。24 道、道距 2 m排列长 48 mΔk 约 0.0208 1/m。第三个是空间混叠条件道距 dx 必须小于目标频率下最低相速度对应波长的一半写成 dx c_min / (2 fmax)。这三个式子直接决定方案能不能用。比如目标最高频率 40 Hz浅层最低相速度 150 m/s算出来道距要小于 1.875 m野外用的 2 m 道距就已经开始混叠了。又比如在 10 Hz 处相速度 200 m/s 对应波数 0.05 1/m250 m/s 对应 0.04 1/m两者波数差只有 0.01小于 24 道排列的 0.0208 1/m常规 f-k 谱分不开这两条能量。这种情况要么加长排列要么后面考虑高分辨率 f-k。3. 频散曲线提取的程序设计从峰值搜索到曲线输出3.1 为什么直接找峰值不行低速能量与体波、混叠的干扰f-k 谱上除了面波频散能量还有很强的直达波、折射波、声波和混叠回卷能量。如果对每个频率直接取全局振幅最大值低频段大概率抓到高视速度的折射波高频段会抓到空间混叠产生的假能量团。程序设计实践里第一件事就是给峰值搜索加约束而不是写一个 np.argmax 就收工。约束分两层第一层是速度窗用 vmin/vmax 把搜索范围限定在面波可能出现的相速度区间第二层是连续性约束下一频率的峰值必须落在上一频率峰值附近否则认为发生跳变。这两层约束缺一不可。3.2 分频带峰值拾取的模块化实现三个子函数与主流程把提取逻辑拆成函数方便后面单独调试和替换。整个程序可以分成三个模块prepare_shot 负责预处理fk_amp_spectrum 负责频谱计算extract_dispersion 负责峰值拾取。核心拾取代码如下。def extract_dispersion(amp, freq, k, fmin, fmax, vmin, vmax, frac0.2): 从 f-k 振幅谱中按频率逐点拾取峰值。 amp: 振幅谱形状 (n_k, n_f) freq: 频率轴单位 Hz k: 正波数轴单位 1/m vmin, vmax: 相速度搜索范围单位 m/s frac: 相邻频率速度变化的最大比例 idx_f np.where((freq fmin) (freq fmax))[0] out [] for jf in idx_f: f freq[jf] # 速度窗换算成波数窗k f / c klow f / vmax # 速度高波数小 khigh f / vmin # 速度低波数大 # 用上一个点的相速度做连续性限制 if len(out) 0: c_prev out[-1][1] klow max(klow, f / (c_prev * (1.0 frac))) khigh min(khigh, f / (c_prev * (1.0 - frac))) ids np.where((k klow) (k khigh))[0] if len(ids) 0: continue jk ids[np.argmax(amp[ids, jf])] out.append((f, f / k[jk])) return np.array(out)主流程调用也很直接amp, freq, k fk_amp_spectrum(data, dt0.001, dx2.0) curves extract_dispersion(amp, freq, k, fmin3.0, fmax40.0, vmin150.0, vmax800.0, frac0.2) if len(curves) 0: # 中值滤波去掉单点毛刺kernel_size 取奇数 curves[:, 1] signal.medfilt(curves[:, 1], kernel_size5) np.savetxt(dispersion_curve.txt, curves, fmt%.4f, headerfreq_hz velocity_m_s, comments)参数说明frac0.2 表示相邻频率点之间的相速度变化不超过 20%。这个值根据地层变化剧烈程度调整地层横向变化大可以放宽到 0.3变化小可以收紧到 0.1。vmin/vmax 一开始可以给宽一点比如 100 到 1000 m/s先看曲线大致落在哪个区间再逐步收紧。如果某个频率在窗口内没有峰值程序会直接跳过后面用中值滤波把缺口填掉一部分。中值滤波核大小 5 是经验值核太大会把真实弯曲细节抹掉太小起不到平滑作用。输出文件是两列文本第一列频率第二列相速度后续反演脚本直接读这个文件。3.3 相速度换算与频散曲线输出把 f-k 索引映射成可反演的 (f, c) 数据相速度换算最容易出错的地方是单位。上面代码里 k 是空间频率单位 1/m所以 c f / k 得到的单位是 m/s。如果从别的程序里拿到的是角波数注意必须用 c 2πf / k_angular。单边排列只取正波数半轴即可如果是中间放炮的双边排列正负波数各有一条能量常见做法是两边分别提取曲线再取平均能抵消一部分震源不对称带来的偏差。输出时建议把表头写好频率列和速度列的单位写清楚因为很多反演程序只按列位置读取不会自动识别单位。f-k 谱本身也可以存一份 npy 或文本矩阵方便后面重新画图、调整速度窗不用每次重新做 FFT。4. 把 f-k 方法用准5 个必调参数与 3 条输入要求4.1 必调参数表道距、排列长度、时窗长度、频率扫描范围、峰值搜索窗这五个参数基本决定了提取结果的可用性也决定了这个方法值不值得在当前数据上花时间。参数经验取值主要影响调整要点道距 dxdx vmin / (2 fmax)空间混叠与波数范围优先满足混叠约束再谈分辨率排列长度 L最低频率波长的 2~4 倍波数分辨率两个模式分不开就加长排列时窗长度 T完整覆盖面波波列频率分辨率截短会让低频峰变大变宽频率范围 fmin/fmax震源主频附近如 3~50 Hz提取范围先看 f-k 谱里亮带边界速度搜索窗 vmin/vmax根据工区浅层速度粗估拾取稳定性宁宽后靠连续性收缩道距的选择没有商量余地。道距太大高频低速面波直接混叠提取出来的曲线在高频端会向上折叠。排列长度决定波数分辨率最低频率的目标波长至少要有 2 个排列长度否则低频端只有两三个波数点根本找不出峰。时窗长度要覆盖整个面波波列但也不能太长把后面噪声也框进来常见做法是先画一张时距图看面波波列的到达时间范围再截取。频率范围不要拍脑袋先把振幅谱画出来看能量亮带在哪个频率区间再设置 fmin/fmax。4.2 炮集质量的第一道关去均值、去趋势与空间加窗的顺序预处理顺序不能乱。先逐道去均值再去趋势然后坏道充零最后才加窗。如果先加窗再去趋势窗函数会把趋势的形状改变去趋势反而去不干净。空间加窗会让排列两端道幅值衰减相当于有效排列长度变短波数分辨率会有轻微损失。为了减少这种影响空间窗不一定用完整汉宁窗可以用 taper 只在两端各衰减 5%~10% 的余弦窗。判断预处理是否合格可以在 f-k 变换前把每道 RMS 振幅画出来如果有某道明显高出周围说明坏道判据没起作用需要调低 bad_scale。注意空间加窗之后排列两端的信息被削弱。如果目标模式能量刚好落在两端道上加窗会把信号压没。遇到这种情况优先检查排列长度是否足够而不是一味加窗。4.3 高阶模式与基阶模式的区分振幅差异和能量连续性怎么配合基阶模式通常是 f-k 谱里振幅最强、速度最低的那条亮带能量从低频连续延伸到高频。高阶模式往往出现在更高频率或者和基阶交叉。最常用的区分方法有两个信号看振幅基阶一般比高阶强看连续性基阶能量条带从头到尾不中断高阶在低频端常常若隐若现。提取时如果想同时提多阶可以把 extract_dispersion 改成多候选模式在每个频率上取前若干个局部极大值再用速度连续性把点连成曲线。这里的关键是速度搜索窗要分阶设置基阶用低速窗高阶把 vmin/vmax 整体抬高否则峰值搜索会一直在强能量附近打转。5. f-k 频散提取常见问题与排查我踩过的 5 个坑这套方法看起来简单真正跑数据时踩坑的地方全在细节里。下面五条都是我用实际炮集调过、返工过的血泪经验。5.1 低频端曲线突然翘起明明是平坦地层也这样现象2~5 Hz 区间提取的相速度忽然从 300 m/s 跳到 800 m/s曲线在低频端整体上翘和地层模型明显矛盾。原因低频段折射波和体波能量很强如果速度窗给得太宽峰值搜索会在低频段抓到高视速度的折射波而不是面波。另一方面记录时窗不够长时低频段频率分辨率差亮带变宽峰值位置容易被旁边能量拉走。解决先把 vmax 压到合理的面波速度上限比如 600 m/s重新提取看曲线是否回落。如果还翘把 fmin 从 2 Hz 提到 4 Hz牺牲低频段换稳定性。再不行在 f-k 变换前做一次带通滤波把 2~5 Hz 以外能量压低。低频端频散信息本身金贵但在信噪比不足时强行保留反演结果会更差。5.2 基阶能量被拉成斜线提取速度整体偏快现象f-k 谱上基阶亮带明显变宽峰值拾取后整条速度偏高而且高频端偏得更明显曲线像是被从左上往右下拉过。原因最常见的是零波数附近直流泄漏。去均值不彻底时每条道残留直流成分会在 k0 附近形成一条强能量带把亮带向低波数方向拉低波数对应高速度。另一个原因是空间加窗过强主瓣变宽峰值从脊线偏到低波数一侧。解决回到预处理逐道做 detrend并检查是否还有常数偏移。空间加窗改成 taper 更小的余弦窗比如 c0.05 的 Tukey 窗只衰减两端 5% 的道。改完后重新画 f-k 谱看亮带宽度是否明显收窄。5.3 f-k 谱面上出现规律性横条纹峰值拾取反复跳动现象谱面上出现平行于波数轴的周期性亮暗条纹像百叶窗一样。自动提取的曲线在几个速度值之间反复跳。原因这种条纹多半来自坏道或强振幅近道。坏道在空间域像一个窄脉冲傅里叶变换后会在波数方向产生周期性旁瓣频率方向上表现为横条纹。另一个常见来源是 50 Hz 工频干扰会在特定频率上形成水平亮线。解决先检查原始道集把坏道充零或插值。50 Hz 干扰用陷波器滤掉。还有一种容易被忽略的情况是各道增益不一致造成空间采样不均匀先做道间振幅归一化再进 f-k。条纹消除后峰值拾取的跳动通常会立刻缓解。5.4 高频端出现折叠回卷的假能量团现象15~30 Hz 处出现一个向低波数方向回卷的假能量团提取曲线在高频端突然掉到很低的速度或跳到别的模式上。原因这是空间混叠本质上是道距太大短波长面波在空间采样时被折叠。表面波波长小于 2dx 时波数超过奈奎斯特波数混叠能量从高波数回卷到低波数正好落在真实基阶能量旁边的位置。解决先检查道距是否满足 dx c_min / (2 fmax)。不满足时有两个选择降低 fmax 到安全范围野外重新采集时缩小道距。道间插值只能让谱面看起来更连续不能创造真实的空间采样信息混叠照样存在。这个坑是物理限制程序怎么改都绕不过去。5.5 自动拾取曲线抖得没法用手工返修一小时现象提取结果在相邻频率之间上蹿下跳曲线像锯齿。中值滤波之后仍然有毛刺手工修起来比从头画还慢。原因单频点峰值信噪比不足是最直接的原因。另一个常见原因是补零太少谱面呈锯齿状峰值位置在网格间来回摆。连续性搜索窗口和频率步长不匹配也会放大抖动。解决把补零倍数提高到 8谱面更平滑。拾取时不要直接用离散最大值改用峰值附近三点的抛物线拟合得到亚网格精度。最后再用 5 点中值滤波。如果抖动仍然严重考虑对谱面做一次沿频率方向的轻平滑内核长度取 5 个频率点但注意不要过度平滑把真实低频细节抹掉。6. 上反演台之前用正演、Capon 与 τ-p 交叉验证6.1 用正演模拟数据验证程序一个四步验证流程程序写完后先用正演数据验证再上真实炮集。常见做法是先给定层状速度模型用反射率法或广义 R/T 系数法正演理论频散曲线再用同一模型合成炮集跑完提取流程后和理论曲线做差。这样能检查参数设置和代码逻辑比直接拿野外数据调参快得多。误差评估只需要一行代码err np.interp(f_theory, curves[:, 0], curves[:, 1]) - c_theory rmse np.sqrt(np.mean(err**2)) / np.mean(c_theory)如果 RMSE 超过 3%先查速度窗和预处理不要急着改反演算法。6.2 高分辨率 f-kCapon 方法什么时候值得上常规 f-k 谱分辨率受排列长度限制两个模式的波数差小于 1/L 时在谱面上就是一团。Capon 方法用数据协方差矩阵做最小方差无失真响应扫描能在短排列下提高波数分辨能力代价是对坏道和噪声非常敏感还需要做对角加载正则化。只有当常规 FFT 谱明显分不开两条能量带时才值得上数据质量太差时一上就翻车。6.3 和 τ-p 变换结果做差一个有用的交叉验证习惯τ-p 变换和 f-k 域在数学上本质是同一信息的不同坐标但实际代码实现里插值方式不同会引入各自的系统偏差。我的习惯是同一炮集分别用 f-k 和 τ-p 提取把两条曲线画在一起统一到相同频率网格后做差。偏差超过 2% 就回头查参数而不是急着判断哪条更准。独立方法互验比盲目相信某一次输出可靠得多。这也是我现在的固定习惯。希望帮到你。本文还有配套的精品资源点击获取
频率-波数域频散曲线提取：f-k谱分析与自动拾取程序设计
简介频率—波数域频散曲线提取方法及程序设计.pdf 是一份面向地震勘探与工程物探技术人员的地学专业文献聚焦瑞利波数据处理中的 f-k 频散曲线提取技术。文档从瑞利波基本概念入手系统讲解二维傅氏变换将时间-空间域数据转换为频率-波数域能量谱的原理并结合半波长理论阐述速度-深度域频散曲线绘制流程其中对能量谱能量最大点选取、视速度与波长计算等关键步骤也做了说明同时以 Delphi7.0 为开发环境给出程序化实现思路与关键模块设计对开展相关算法研究、程序编写或工程应用具有直接参考价值。资源包内包含 1 个 PDF 文件文件大小约 365KB内容精炼、图文结合便于下载后快速阅读与存档。目前已有 551 人浏览学习适合物探、岩土工程及地震数据处理方向的师生或工程师作为专业指导文献使用。1. 为什么是频率—波数域一张 f-k 谱能省掉半天的排列试验处理主动源面波数据时我一直有个习惯拿到炮集先不急着画时距图直接做一次二维傅里叶变换到频率—波数域。原因很简单时距图里面波和折射波叠在一起肉眼很难分清哪条能量是基阶频散而在 f-k 谱里不同视速度的信号按斜率分开面波频散能量会聚成一条窄窄的亮带。频率—波数域频散曲线提取方法及程序设计做的就是这件事把时间-空间域的炮集变成 f-k 谱再从谱面上自动提取频散曲线供后续反演使用。这套方案适合做 MASW、主动源面波勘探的工程师和研究人员也适合刚接触频散分析、想摆脱手工拾取的人。读完之后你至少能自己写出一套能跑的提取程序并且知道参数该往哪个方向调。2. 从时间-空间波场到 f-k 谱二维傅里叶变换怎么做才不白做2.1 一张地震记录在 f-k 域里长什么样能量条带与速度线的对应地震记录可以看成二维函数 u(x, t)x 是检波点位置t 是时间。对时间做傅里叶变换得到频率 f对空间做傅里叶变换得到波数 k二维变换之后就是 U(k, f)也就是频率—波数域谱。一个沿 x 方向以速度 c 传播的平面简谐波在 f-k 谱上对应一条过原点的直线斜率就是速度。体波、折射波视速度高直线更陡能量靠近波数轴面波视速度低能量靠近频率轴。因为频散波的相速度随频率变化它的能量不会聚成直线而是形成一条弯曲的亮带。频散曲线提取就是从这条亮带的脊线上按频率逐个取相速度。f-k 域最大的价值在于把不同视速度的能量在平面上分开自动拾取比在时间域里追相位转折点可靠得多。2.2 用 Python 把波场变换到 f-k 域预处理、补零与 FFT 的完整代码这里给出一段可以直接落地的实现。Python 程序设计里最常用的就是 numpy 和 scipy二维 FFT 一行就能算完真正影响结果的是变换前的预处理。import numpy as np from scipy import signal def prepare_shot(data, bad_scale10.0): 把原始炮集整理成适合做 2D FFT 的形态。 data: 二维数组形状 (nx, nt)nx 为道数nt 为采样点数。 nx, nt data.shape # 逐道去均值、去趋势避免零频附近直流泄漏 data data - np.mean(data, axis1, keepdimsTrue) data signal.detrend(data, axis1, typelinear) # 坏道判定振幅超过中位数 10 倍的道直接充零 amp np.max(np.abs(data), axis1) bad amp bad_scale * np.median(amp) data[bad, :] 0.0 # 时间维加窗抑制谱泄漏空间维加窗抑制排列两端截断效应 data * np.hanning(nt)[None, :] data * np.hanning(nx)[:, None] return data def fk_amp_spectrum(data, dt, dx, x_pad4, t_pad4): 对预处理后的炮集做二维 FFT返回振幅谱和坐标轴。 dt: 时间采样间隔(s) dx: 道距(m) nx, nt data.shape # 补零到 2 的幂FFT 更快谱面也更平滑 nfx int(2 ** np.ceil(np.log2(nx * x_pad))) nft int(2 ** np.ceil(np.log2(nt * t_pad))) spec np.fft.fft2(data, s(nfx, nft)) spec np.fft.fftshift(spec, axes0) # 只对波数轴做中心化 freq np.fft.fftfreq(nft, dt) # 频率轴单位 Hz k np.fft.fftfreq(nfx, dx) # 波数轴单位 1/m half nfx // 2 return np.abs(spec[half:, :]), freq, k[half:]先说逻辑。prepare_shot 里先做逐道去均值和线性去趋势这一步不能省否则零波数附近会残留一条强能量把低频端曲线往上拉。坏道充零要在加窗之前因为加窗会把坏道异常值扩散到相邻频率。时间维和空间维都用了汉宁窗汉宁窗主瓣略宽但旁瓣低比矩形窗更适合自动拾取。再看 fk_amp_spectrum。fft2 的第一个维度是道数也就是空间维s(nfx, nft) 指定补零后的形状。fftshift 只在波数轴移动把零波数放到中间然后只取后半段对应正波数。这样做的原因是常见单边排列里波只沿一个方向传播负波数半轴基本只有对称的镜像能量。返回的 k 单位是 cycles/m注意这里不是角波数所以后面相速度直接 c f / k不用乘 2π。提示np.fft.fftfreq 生成的是空间频率单位 1/m。如果换成角波数 rad/m相速度换算要写成 c 2πf / k_angular。x_pad 和 t_pad 是补零倍数一般取 4 到 8。补零能加密谱线让峰值拾取更光滑但它不改变真实分辨率真实分辨率只由原始道数和排列长度决定。频率轴里包含负频率后面提取函数会通过 fmin/fmax 把它们排除掉。2.3 f-k 谱分辨率到底被什么锁死三个算式与最小道距判定在做 f-k 分析前先用三个算式估算谱的分辨率和可用范围能避免后面反复调参。第一个是频率分辨率Δf 1 / (nt·dt)也就是记录总时间长度的倒数。1 秒的记录Δf 就是 1 Hz想在 2 Hz 附近分辨两个很近的模式基本没戏。第二个是波数分辨率Δk 1 / (nx·dx)也就是排列长度的倒数。24 道、道距 2 m排列长 48 mΔk 约 0.0208 1/m。第三个是空间混叠条件道距 dx 必须小于目标频率下最低相速度对应波长的一半写成 dx c_min / (2 fmax)。这三个式子直接决定方案能不能用。比如目标最高频率 40 Hz浅层最低相速度 150 m/s算出来道距要小于 1.875 m野外用的 2 m 道距就已经开始混叠了。又比如在 10 Hz 处相速度 200 m/s 对应波数 0.05 1/m250 m/s 对应 0.04 1/m两者波数差只有 0.01小于 24 道排列的 0.0208 1/m常规 f-k 谱分不开这两条能量。这种情况要么加长排列要么后面考虑高分辨率 f-k。3. 频散曲线提取的程序设计从峰值搜索到曲线输出3.1 为什么直接找峰值不行低速能量与体波、混叠的干扰f-k 谱上除了面波频散能量还有很强的直达波、折射波、声波和混叠回卷能量。如果对每个频率直接取全局振幅最大值低频段大概率抓到高视速度的折射波高频段会抓到空间混叠产生的假能量团。程序设计实践里第一件事就是给峰值搜索加约束而不是写一个 np.argmax 就收工。约束分两层第一层是速度窗用 vmin/vmax 把搜索范围限定在面波可能出现的相速度区间第二层是连续性约束下一频率的峰值必须落在上一频率峰值附近否则认为发生跳变。这两层约束缺一不可。3.2 分频带峰值拾取的模块化实现三个子函数与主流程把提取逻辑拆成函数方便后面单独调试和替换。整个程序可以分成三个模块prepare_shot 负责预处理fk_amp_spectrum 负责频谱计算extract_dispersion 负责峰值拾取。核心拾取代码如下。def extract_dispersion(amp, freq, k, fmin, fmax, vmin, vmax, frac0.2): 从 f-k 振幅谱中按频率逐点拾取峰值。 amp: 振幅谱形状 (n_k, n_f) freq: 频率轴单位 Hz k: 正波数轴单位 1/m vmin, vmax: 相速度搜索范围单位 m/s frac: 相邻频率速度变化的最大比例 idx_f np.where((freq fmin) (freq fmax))[0] out [] for jf in idx_f: f freq[jf] # 速度窗换算成波数窗k f / c klow f / vmax # 速度高波数小 khigh f / vmin # 速度低波数大 # 用上一个点的相速度做连续性限制 if len(out) 0: c_prev out[-1][1] klow max(klow, f / (c_prev * (1.0 frac))) khigh min(khigh, f / (c_prev * (1.0 - frac))) ids np.where((k klow) (k khigh))[0] if len(ids) 0: continue jk ids[np.argmax(amp[ids, jf])] out.append((f, f / k[jk])) return np.array(out)主流程调用也很直接amp, freq, k fk_amp_spectrum(data, dt0.001, dx2.0) curves extract_dispersion(amp, freq, k, fmin3.0, fmax40.0, vmin150.0, vmax800.0, frac0.2) if len(curves) 0: # 中值滤波去掉单点毛刺kernel_size 取奇数 curves[:, 1] signal.medfilt(curves[:, 1], kernel_size5) np.savetxt(dispersion_curve.txt, curves, fmt%.4f, headerfreq_hz velocity_m_s, comments)参数说明frac0.2 表示相邻频率点之间的相速度变化不超过 20%。这个值根据地层变化剧烈程度调整地层横向变化大可以放宽到 0.3变化小可以收紧到 0.1。vmin/vmax 一开始可以给宽一点比如 100 到 1000 m/s先看曲线大致落在哪个区间再逐步收紧。如果某个频率在窗口内没有峰值程序会直接跳过后面用中值滤波把缺口填掉一部分。中值滤波核大小 5 是经验值核太大会把真实弯曲细节抹掉太小起不到平滑作用。输出文件是两列文本第一列频率第二列相速度后续反演脚本直接读这个文件。3.3 相速度换算与频散曲线输出把 f-k 索引映射成可反演的 (f, c) 数据相速度换算最容易出错的地方是单位。上面代码里 k 是空间频率单位 1/m所以 c f / k 得到的单位是 m/s。如果从别的程序里拿到的是角波数注意必须用 c 2πf / k_angular。单边排列只取正波数半轴即可如果是中间放炮的双边排列正负波数各有一条能量常见做法是两边分别提取曲线再取平均能抵消一部分震源不对称带来的偏差。输出时建议把表头写好频率列和速度列的单位写清楚因为很多反演程序只按列位置读取不会自动识别单位。f-k 谱本身也可以存一份 npy 或文本矩阵方便后面重新画图、调整速度窗不用每次重新做 FFT。4. 把 f-k 方法用准5 个必调参数与 3 条输入要求4.1 必调参数表道距、排列长度、时窗长度、频率扫描范围、峰值搜索窗这五个参数基本决定了提取结果的可用性也决定了这个方法值不值得在当前数据上花时间。参数经验取值主要影响调整要点道距 dxdx vmin / (2 fmax)空间混叠与波数范围优先满足混叠约束再谈分辨率排列长度 L最低频率波长的 2~4 倍波数分辨率两个模式分不开就加长排列时窗长度 T完整覆盖面波波列频率分辨率截短会让低频峰变大变宽频率范围 fmin/fmax震源主频附近如 3~50 Hz提取范围先看 f-k 谱里亮带边界速度搜索窗 vmin/vmax根据工区浅层速度粗估拾取稳定性宁宽后靠连续性收缩道距的选择没有商量余地。道距太大高频低速面波直接混叠提取出来的曲线在高频端会向上折叠。排列长度决定波数分辨率最低频率的目标波长至少要有 2 个排列长度否则低频端只有两三个波数点根本找不出峰。时窗长度要覆盖整个面波波列但也不能太长把后面噪声也框进来常见做法是先画一张时距图看面波波列的到达时间范围再截取。频率范围不要拍脑袋先把振幅谱画出来看能量亮带在哪个频率区间再设置 fmin/fmax。4.2 炮集质量的第一道关去均值、去趋势与空间加窗的顺序预处理顺序不能乱。先逐道去均值再去趋势然后坏道充零最后才加窗。如果先加窗再去趋势窗函数会把趋势的形状改变去趋势反而去不干净。空间加窗会让排列两端道幅值衰减相当于有效排列长度变短波数分辨率会有轻微损失。为了减少这种影响空间窗不一定用完整汉宁窗可以用 taper 只在两端各衰减 5%~10% 的余弦窗。判断预处理是否合格可以在 f-k 变换前把每道 RMS 振幅画出来如果有某道明显高出周围说明坏道判据没起作用需要调低 bad_scale。注意空间加窗之后排列两端的信息被削弱。如果目标模式能量刚好落在两端道上加窗会把信号压没。遇到这种情况优先检查排列长度是否足够而不是一味加窗。4.3 高阶模式与基阶模式的区分振幅差异和能量连续性怎么配合基阶模式通常是 f-k 谱里振幅最强、速度最低的那条亮带能量从低频连续延伸到高频。高阶模式往往出现在更高频率或者和基阶交叉。最常用的区分方法有两个信号看振幅基阶一般比高阶强看连续性基阶能量条带从头到尾不中断高阶在低频端常常若隐若现。提取时如果想同时提多阶可以把 extract_dispersion 改成多候选模式在每个频率上取前若干个局部极大值再用速度连续性把点连成曲线。这里的关键是速度搜索窗要分阶设置基阶用低速窗高阶把 vmin/vmax 整体抬高否则峰值搜索会一直在强能量附近打转。5. f-k 频散提取常见问题与排查我踩过的 5 个坑这套方法看起来简单真正跑数据时踩坑的地方全在细节里。下面五条都是我用实际炮集调过、返工过的血泪经验。5.1 低频端曲线突然翘起明明是平坦地层也这样现象2~5 Hz 区间提取的相速度忽然从 300 m/s 跳到 800 m/s曲线在低频端整体上翘和地层模型明显矛盾。原因低频段折射波和体波能量很强如果速度窗给得太宽峰值搜索会在低频段抓到高视速度的折射波而不是面波。另一方面记录时窗不够长时低频段频率分辨率差亮带变宽峰值位置容易被旁边能量拉走。解决先把 vmax 压到合理的面波速度上限比如 600 m/s重新提取看曲线是否回落。如果还翘把 fmin 从 2 Hz 提到 4 Hz牺牲低频段换稳定性。再不行在 f-k 变换前做一次带通滤波把 2~5 Hz 以外能量压低。低频端频散信息本身金贵但在信噪比不足时强行保留反演结果会更差。5.2 基阶能量被拉成斜线提取速度整体偏快现象f-k 谱上基阶亮带明显变宽峰值拾取后整条速度偏高而且高频端偏得更明显曲线像是被从左上往右下拉过。原因最常见的是零波数附近直流泄漏。去均值不彻底时每条道残留直流成分会在 k0 附近形成一条强能量带把亮带向低波数方向拉低波数对应高速度。另一个原因是空间加窗过强主瓣变宽峰值从脊线偏到低波数一侧。解决回到预处理逐道做 detrend并检查是否还有常数偏移。空间加窗改成 taper 更小的余弦窗比如 c0.05 的 Tukey 窗只衰减两端 5% 的道。改完后重新画 f-k 谱看亮带宽度是否明显收窄。5.3 f-k 谱面上出现规律性横条纹峰值拾取反复跳动现象谱面上出现平行于波数轴的周期性亮暗条纹像百叶窗一样。自动提取的曲线在几个速度值之间反复跳。原因这种条纹多半来自坏道或强振幅近道。坏道在空间域像一个窄脉冲傅里叶变换后会在波数方向产生周期性旁瓣频率方向上表现为横条纹。另一个常见来源是 50 Hz 工频干扰会在特定频率上形成水平亮线。解决先检查原始道集把坏道充零或插值。50 Hz 干扰用陷波器滤掉。还有一种容易被忽略的情况是各道增益不一致造成空间采样不均匀先做道间振幅归一化再进 f-k。条纹消除后峰值拾取的跳动通常会立刻缓解。5.4 高频端出现折叠回卷的假能量团现象15~30 Hz 处出现一个向低波数方向回卷的假能量团提取曲线在高频端突然掉到很低的速度或跳到别的模式上。原因这是空间混叠本质上是道距太大短波长面波在空间采样时被折叠。表面波波长小于 2dx 时波数超过奈奎斯特波数混叠能量从高波数回卷到低波数正好落在真实基阶能量旁边的位置。解决先检查道距是否满足 dx c_min / (2 fmax)。不满足时有两个选择降低 fmax 到安全范围野外重新采集时缩小道距。道间插值只能让谱面看起来更连续不能创造真实的空间采样信息混叠照样存在。这个坑是物理限制程序怎么改都绕不过去。5.5 自动拾取曲线抖得没法用手工返修一小时现象提取结果在相邻频率之间上蹿下跳曲线像锯齿。中值滤波之后仍然有毛刺手工修起来比从头画还慢。原因单频点峰值信噪比不足是最直接的原因。另一个常见原因是补零太少谱面呈锯齿状峰值位置在网格间来回摆。连续性搜索窗口和频率步长不匹配也会放大抖动。解决把补零倍数提高到 8谱面更平滑。拾取时不要直接用离散最大值改用峰值附近三点的抛物线拟合得到亚网格精度。最后再用 5 点中值滤波。如果抖动仍然严重考虑对谱面做一次沿频率方向的轻平滑内核长度取 5 个频率点但注意不要过度平滑把真实低频细节抹掉。6. 上反演台之前用正演、Capon 与 τ-p 交叉验证6.1 用正演模拟数据验证程序一个四步验证流程程序写完后先用正演数据验证再上真实炮集。常见做法是先给定层状速度模型用反射率法或广义 R/T 系数法正演理论频散曲线再用同一模型合成炮集跑完提取流程后和理论曲线做差。这样能检查参数设置和代码逻辑比直接拿野外数据调参快得多。误差评估只需要一行代码err np.interp(f_theory, curves[:, 0], curves[:, 1]) - c_theory rmse np.sqrt(np.mean(err**2)) / np.mean(c_theory)如果 RMSE 超过 3%先查速度窗和预处理不要急着改反演算法。6.2 高分辨率 f-kCapon 方法什么时候值得上常规 f-k 谱分辨率受排列长度限制两个模式的波数差小于 1/L 时在谱面上就是一团。Capon 方法用数据协方差矩阵做最小方差无失真响应扫描能在短排列下提高波数分辨能力代价是对坏道和噪声非常敏感还需要做对角加载正则化。只有当常规 FFT 谱明显分不开两条能量带时才值得上数据质量太差时一上就翻车。6.3 和 τ-p 变换结果做差一个有用的交叉验证习惯τ-p 变换和 f-k 域在数学上本质是同一信息的不同坐标但实际代码实现里插值方式不同会引入各自的系统偏差。我的习惯是同一炮集分别用 f-k 和 τ-p 提取把两条曲线画在一起统一到相同频率网格后做差。偏差超过 2% 就回头查参数而不是急着判断哪条更准。独立方法互验比盲目相信某一次输出可靠得多。这也是我现在的固定习惯。希望帮到你。本文还有配套的精品资源点击获取
