简介一份关于频率—波数域频散曲线提取方法及程序设计的学术论文面向地震勘探、工程物探及瑞利波数据处理方向的研究人员与专业技术人员。资源包共1个文件为PDF格式大小约365KB目前已有551人学习下载。内容系统论述了基于二维傅氏变换与半波长理论的频散曲线提取原理并给出基于Delphi7.0平台的程序设计思路与实现过程涉及f-k能量谱分析、速度—深度域曲线绘制等关键环节可帮助读者理解多道瑞利波数据处理的完整技术路径。该方法利用瑞利波频散特性与传播速度同介质物理性质的关联可应用于地层划分、确定地基夯实效果、划分软弱地层埋深与范围、评价基岩完整性等工程场景尤其对需要掌握f-k法程序实现的工程师具有直接参考价值。1. 频率—波数域频散曲线提取从时距记录到相速度数据之间隔着什么做过浅层地震面波勘探的人基本都经历过这种场面野外采回来的炮集在时距图上就是一团连续抖动的扫帚眼睛很难从里面直接读出地层速度的变化。可一旦把整张记录做一次二维傅里叶变换从时间—空间域转到频率—波数域大家都叫 f-k 域情况立刻变得清爽——面波能量会沿着一条或者几条弯曲的脊线集中分布脊线的几何形态就是频散关系的直接体现。频率—波数域频散曲线提取要做的事就是把这些脊线逐频点拾取出来换算成频率—相速度数据对交给后面的剪切波速度反演。这篇技术笔记面向的读者是正在做工程物探、场地波速测试、路面或隧道检测的工程师和学生你们手里有炮集或面波记录想不靠商业软件自己把频散曲线提出来那需要先弄清原理、参数再躲开那些让人翻车的坑。2. f-k域频散谱的物理含义面波能量为什么会在谱上连成一条脊线2.1 二维傅里叶变换把时距关系翻译成频率—波数关系一个沿水平方向传播的简谐面波可以写成 d(t,x)A·exp[i(ωt−kx)]ω2πf 是角频率k 是角波数相速度 vω/kf/κκ 是循环空间频率单位 cycles/m与角波数相差 2π 倍。对于水平层状介质瑞利波存在频散不同频率对应不同相速度因此每个频率 f 下的面波都有固定的波数 k(f) 2πf/v(f)。检波器以道间距 dx 等间隔采样记录按采样率 dt 的时间网格排列整个炮集就是一组简谐分量的叠加而二维傅里叶变换恰好能把叠加打散把能量按频率和波数重新归位。二维傅里叶变换的离散形式写成 F(f,κ)Σ_t Σ_x d(t,x)·exp[−i2π(ftκx)]。入射波场里一个沿固定速度传播的平面波变换后会变成一条过原点的亮直线如果速度随频率变化这条亮线就弯曲成脊线。在 f-k 振幅谱上脊线的峰值位置就是该频率下主要模态的波数。这里要提醒一句实际程序里用的是快速傅里叶变换做两个维度的变换默认把零频、零波数放在数组角上所以变换完必须做 fftshift把原点挪到谱中心再可视化或拾取否则峰值搜索的坐标全部错位。从程序设计实践的角度看把变换、取模、坐标轴生成写成独立函数比在一个循环里堆几层 for 要稳得多后面第 3 章给的代码就是按这个思路拆的。理解这一步的关键不在于背公式而在于建立直觉时距图里斜率越陡的线性同相轴变换后越靠近零波数频率越高、波数越大对应的相速度就越低。基阶面波在高频段能量集中在浅层低速区所以在谱上会看到脊线随频率升高而向右上方大波数方向弯曲。2.2 相速度与波数的换算vf/κ 这条公式的适用边界用 numpy 的 fftfreq 生成坐标轴时时间轴的单位是 Hz空间轴的单位是 cycles/m两者直接相处 vf/κ 就得到相速度不需要再乘 2π。这个简洁关系成立的前提是介质水平分层、排列方向与传播方向一致、检波点在直线上等间隔布设。野外如果测线沿地形起伏布设或者排列中间有几道被障碍物挪位波数轴的几何关系会被破坏谱峰位置偏移提取出的速度值就不准。常见做法是在预处理里先做道头校正或重排把偏移量折回规则网格。另一个隐含条件是远场平面波近似。面波在近炮检距处包含近场成分波前面不是平面二维傅里叶变换的平面波分解语义会打折扣工程实践里一般把炮检距最近的道控制在 510 m 以外把前几道切除或者在后续拾取时只参考中层道的能量分布。还要注意波数轴的正负正波数对应沿排列正方向传播的能量负波数对应反向传播。单边激发时反方向能量主要来自端射和散射拾取时只取正波数半边通常更干净。2.3 基阶和高阶模态在 f-k 谱上的位置规律同一个地层模型下基阶模态的相速度最低而一阶、二阶等高阶模态相速度更高。换算到 f-k 谱上基阶脊线更靠近大波数侧右上方高阶脊线则依次偏向小波数侧左下方。振幅方面基阶能量往往最强但在高频段高阶模态能量不一定弱尤其是覆盖层与半空间速度反差大的时候高阶分支甚至会成为局部最大峰值。这就引出一个频散提取里最典型的判断问题单纯做全局峰值搜索很容易在某个频点跳到高阶分支上导致频散曲线断裂成不连续的折线。我在实际处理中会先看整张谱的 grayscale 图确认哪条脊线从低频到高频是连续的再决定拾取哪条分支如果两条分支靠得近就分别限制速度搜索范围一条一条提。后面的第 3.3 节代码也是按这个策略设计的先给搜索窗再找峰而不是搜完全谱再猜模态。3. 程序设计用一段可运行的 Python 把 f-k 频散提取流程串起来这一章给一套能在本地直接跑通的最小实现。数据形式假设为你已经整理好的二维数组 data形状是 (nt, nx)第一维是时间采样第二维是空间道号dt 是时间采样间隔dx 是道间距。如果手里是 SEG-Y 或 MiniSEED常见做法先用 ObsPy 或 segyio 读进来再转成 numpy 数组后面的流程完全不用改。3.1 数据读入与预处理去直流、时窗、道均衡一个都不能少先看预处理函数。很多炮集里近道振幅远大于远道直接做傅里叶变换会让谱上的能量峰被拉宽波数分辨率变差所以要按道做 RMS 归一化面波时窗外若有折射波或声波到达会给谱里叠加强干扰所以用汉宁窗把记录两端削掉。代码里我顺手把时间轴的线性趋势也去掉了。import numpy as np from scipy.signal import detrend def preprocess_record(data, dt, win_start0.0, win_endNone): 对原始炮记录做预处理. data: (nt, nx) 二维数组, 第一维时间, 第二维道号 dt: 时间采样间隔(秒) win_start, win_end: 保留的时间窗(秒), 一般圈面波到达段 nt, nx data.shape # 1) 去常数和线性趋势, 消除直流漂移与超低频扰动 data detrend(data, axis0, typeconstant) data detrend(data, axis0, typelinear) # 2) 时间窗截取: 超出面波到达范围的能量直接清零 if win_end is None: win_end (nt - 1) * dt i0 int(max(win_start / dt, 0)) i1 int(min(win_end / dt, nt)) data[:i0, :] 0.0 data[i1:, :] 0.0 # 3) 汉宁窗: 压低截断旁瓣, 代价是频率分辨率略降 taper np.hanning(i1 - i0) data[i0:i1, :] * taper[:, None] # 4) 道均衡: 每道按RMS归一, 防止近道能量压过远道 rms np.sqrt(np.mean(data ** 2, axis0)) 1e-12 data data / rms[None, :] return data逻辑说明detrend 按时间轴做两次去趋势第一次去掉常数直流第二次去掉线性漂移这两步对手持锤击或落重这种能量不均衡的数据特别有效。时间窗 i0、i1 的圈定很关键——窗太宽会把折射波和声波放进谱里窗太窄则低频成分被削掉。常见做法是先画一张 wiggle 图把面波锥明显到达的区间读出来填进 win_start 和 win_end后续再微调。参数说明dt 要和数据匹配若读数据时已经做了抽稀这里必须用抽稀后的实际采样间隔rms 归一化会改变振幅绝对值但对峰值搜索没有影响因为每个频率的谱峰比较是相对强弱汉宁窗会让主瓣变宽对能量本来就弱的低频段影响更明显所以后面补零时通常把时窗截取后的记录再延拓一倍。3.2 计算 f-k 谱二维 FFT 与频率波数轴的生成预处理之后进入变换。这里用一个独立函数计算 f-k 谱把坐标轴和振幅谱一起返回。注意 fftfreq 生成的是循环频率空间轴的单位是 cycles/m这样后面 vf/κ 直接换速度。def compute_fk_spectrum(data, dt, dx): 二维傅里叶变换并返回频率轴、空间频率轴和振幅谱. data: 预处理后 (nt, nx) 数组 dt: 时间采样间隔(秒); dx: 道间距(米) nt, nx data.shape # 二维FFT, 第0轴是时间, 第1轴是空间 fk np.fft.fft2(data, axes(0, 1)) fk np.fft.fftshift(fk) # 零频移到谱中心 # 频率轴: 范围[-1/(2dt), 1/(2dt)], 单位Hz freqs np.fft.fftshift(np.fft.fftfreq(nt, ddt)) # 空间频率轴: 范围[-1/(2dx), 1/(2dx)], 单位cycles/m kappa np.fft.fftshift(np.fft.fftfreq(nx, ddx)) amp np.abs(fk) # 振幅谱 return freqs, kappa, amp逻辑说明fft2 一次完成两个维度的傅里叶变换axes(0,1) 明确告诉 numpy 第 0 轴是时间、第 1 轴是空间避免二维数组转置后轴顺序混掉。fftshift 在这里做了两次一次是谱的移中一次是坐标轴的移中两者顺序必须一致否则谱峰和坐标错位。amp 取模后单位是时域振幅累积值对搜索无影响如果谱的动态范围太大可以在搜索前做 log10 压缩。参数说明dx 的单位直接决定波数轴范围和最终相速度值填错会整体偏移如果野外道间距并不均匀比如坏道导致间隔是 2m 和 4m 交替代码仍按等间隔处理谱上会出现虚假的旁瓣峰值稳妥做法是先把排列重抽样到等间隔道网格。fftfreq 返回的是双半边坐标已经按 fftshift 对齐到 [-0.5/dt, 0.5/dt]其中负半轴对应反方向传播的能量拾取时可以只取 kappa0 的半边。3.3 峰值搜索与频散曲线输出从谱峰到频率-相速度表得到振幅谱后核心拾取函数对每个频率切一刀在允许的波数范围内找最大峰值并用抛物线插值把峰位修细。为了不让全局最大值乱跳到高阶模态搜索之前先用一个速度窗把范围卡死。def extract_dispersion(freqs, kappa, amp, fmin, fmax, vmin, vmax): 逐频率搜索谱峰并换算相速度. 返回数组, 每行三项: 频率(Hz), 相速度(m/s), 空间频率(cycles/m) dk kappa[1] - kappa[0] results [] for idx_f, f in enumerate(freqs): if f fmin or f fmax: continue # 当前频率的一维切片 slice_amp amp[idx_f, :].copy() # 速度窗换算成波数搜索区间: v f/kappa kappa f/v k_lo f / vmax # 高速对应小波数 k_hi f / vmin # 低速对应大波数 mask (kappa k_lo) (kappa k_hi) slice_amp[~mask] 0.0 # 正波数半边, 排除反向传播能量 slice_amp[kappa 0] 0.0 idx int(np.argmax(slice_amp)) if slice_amp[idx] 0.0: continue # 抛物线插值: 用相邻三个样点拟合亚波数峰值位置 if 1 idx len(kappa) - 2: a slice_amp[idx - 1] b slice_amp[idx] c slice_amp[idx 1] denom a - 2.0 * b c if abs(denom) 1e-12: delta 0.5 * (a - c) / denom if abs(delta) 1.0: idx_frac idx delta else: idx_frac float(idx) else: idx_frac float(idx) k_peak kappa[int(np.floor(idx_frac))] (idx_frac - np.floor(idx_frac)) * dk else: k_peak kappa[idx] if k_peak 0: results.append((f, f / k_peak, k_peak)) return np.array(results)逻辑说明核心是 mask 那段——把用户给的速度范围换算成波数上下界再把范围之外的振幅置零。这样即使谱上有声波、空气波等强能量只要不在速度窗内就不会被搜到。正波数半边处理把反向传播的能量排除这对单边激发数据很关键。抛物线插值是对克制的改进离散谱峰落在某两个采样点之间时三点拟合出的亚波数位置比整数网格更接近真实峰位在高频段、谱峰较陡时能明显减少速度抖动。参数说明fmin/fmax 一般取谱上信噪比高的频带浅层面波通常给 380 Hzvmin/vmax 先按工区经验给宽范围比如 1201200 m/s跑完看结果再收紧防止拾取线在低频端跳跃。对高频段vmin 还要受空间混叠限制具体见第 4.1 节。输出用 np.savetxt() 存成三列文本即可常见做法是把结果写成 (f, v) 两列方便直接导入反演程序。4. 参数怎么设采样间隔、道间距和排列长度对提取结果的硬约束f-k 域频散提取不是数据丢进去就出结果的黑匣子谱的质量由几个观测参数直接决定。这些参数多数在野外采集时就已固定室内只能通过截窗、补零来补救。下面这张总结表先给出关系再逐条展开。参数符号对 f-k 谱的影响硬约束道间距dx决定最大可分辨波数相速度不能低于 2 f dx排列长度Lnx·dx决定波数分辨率低频需要大 L时间采样率dt决定最大可分辨频率f_max≤1/(2dt)记录时长Tnt·dt决定频率分辨率Δf1/T4.1 道间距 dx 决定最大可分辨波数限制的是最小相速度空间采样同样有奈奎斯特限波数最高只能分辨到 κ_nyq1/(2dx) cycles/m超过这个值的能量会折叠回低波数区域表现成高频段谱峰突然兜回来形成假分支。把 κ_nyq 代进相速度公式得到可分辨的最小相速度 v_min2f·dx。以 dx2m、频率 50Hz 为例v_min200m/s如果浅层表土剪切波速度只有 150m/s那 50Hz 以上的基阶相速度低于 200m/s必然发生混叠。这也解释了一个常见现象某些工区用 2m 道距采的面波记录高频段提出来速度值一路走高看起来像速度倒转其实是混叠造成的假象。室内处理能做的有效手段是把采集道距不满足 v_min 约束的高频部分直接截掉或者按相邻两道求和合并成 4m 等效道距提高道距会加重混叠必须配合低通滤波再在拾取结果里删除对应频段。真正要从根源上解决只能在采集设计阶段按目标最浅层速度反算道间距。4.2 记录长度和时窗决定频率分辨率频率分辨率由记录时长决定Δf1/T。2 秒记录对应 0.5Hz 分辨率对 5Hz 以下的低频段来说谱峰的半宽度可能占掉好几十个百分比低频端频散曲线看起来就像被低通滤波过一样。野外面波记录经常只记录了 12 秒低频信息不是被仪器滤掉而是被时间窗截断了。室内处理常见的补救办法是零填充把时窗截取后的数据在时间方向补零到 4 倍长度再做二维 FFT。补零不增加真实信息但能让频率轴加密谱峰位置更平滑对自动拾取有帮助。需要注意的是如果截取面波到时窗本身只有 0.5 秒补零后频率分辨率仍由 0.5 秒决定只是插值点变密真正有效的低频上限大约在 2Hz 左右。此外时间采样率 dt 决定 f_max1/(2dt)工程仪器普遍用 0.5ms 或 1ms 采样对 100Hz 以内的面波完全够用一般不是瓶颈。4.3 速度搜索范围与模态识别策略速度窗的设置直接影响拾取结果在哪条模态上落脚。常见做法是先用宽窗跑一遍得到一条毛刺很多的曲线再根据曲线的大致走向把窗收紧比如基阶曲线在 560Hz 内从 500m/s 降到 250m/s就设 vmin200、vmax600重新提取一遍。收紧后谱峰不会跳到高阶分支曲线连续性明显变好。如果谱上明显存在两条靠近的脊线可以用两套搜索窗分别提取基阶和二阶模态。还有一种更稳的方法对振幅谱做对数压缩后在每个频点搜索局部极大值比如取前三个峰再按相邻频率峰值连续性原则做路径追踪而不是每频率独立取全局最大。这个思路在模态混叠时很管用但代码量比全局搜索多出一截实际工程里多数数据用速度窗全局峰就够了只在工区存在明显高速夹层时才需要多模态追踪。5. f-k 频散提取的常见翻车与排查五条真实踩坑记录5.1 现象f-k 谱被水平亮条纹盖住脊线看不清谱图上出现贯穿所有波数的水平亮条纹本质是时间方向上的强能量截断或周期性扰动要么是记录存在大幅度直流漂移要么是面波时窗外还有强噪声要么是某道放大器饱和产生方波状记录。时间窗边缘的硬截断也会在谱上形成垂直于频率轴的旁瓣条纹。排查顺序先看原始道集的 wiggle 图确认是否有坏道和饱和道再做 detrend 和时窗平滑把窗边缘的汉宁长度从 5% 加到 15%最后逐道检查 RMS 异常大的道直接置零或剔除。这三个步骤做完条带通常消失。还有一种情况是 50Hz 工频干扰会在固定频率位置形成亮点和我们要的脊线叠加处理时用陷波滤波单独切除。5.2 现象频率越高能量峰越向右歪或者某段频率后峰位突然跳回低波数这是空间混叠的典型表现。波数超过奈奎斯特限后折叠回低波数侧谱峰看起来猝不及防地出现在不该出现的位置。判断方法很简单把提取结果里的 (f, v) 和 v_min2f·dx 画在同一张图上凡是落到该曲线之下速度更低的点基本可以判定为混叠产物。解决方法是把超过混叠界限的频率点从结果里删掉。如果非要保住高频段只能从采集端改缩小道间距 dx 或采用不等距排列。有些工区用地表低速带测出浅层速度只有 130m/s用 2m 道距在 40Hz 以上就没法采了这类数据要直接说明高频截止频率而不是硬提。5.3 现象低频段峰值左右乱跳提取曲线不连续低频段波数小谱峰离开原点不远空间孔径 Lnx·dx 如果不够大波数分辨率 Δκ1/L 就比峰位间距还粗这时相邻几个频率的峰一会落在 0.01cycles/m一会落在 0.03cycles/m折线很难看。另一个原因是低频能量本身弱时窗截断后低频段信噪比下降。解决思路分两步先补零到 2 倍或 4 倍长度让谱在波数方向更平滑再做 35 点频率方向中值滤波把孤立跳点压掉。如果低频段必须要到 3Hz而排列只有 24 道×2m48m物理上就不足以分辨这个波数滤波只是让曲线好看信息其实已经丢了。这种情况不如主动把低频截止提到 5Hz保证提出来的是可靠数据。5.4 现象提取出的频散曲线低速段一直压在某个固定速度上常发生在 330350m/s 附近这条平台多半是空气波。空气波在时距图上是一条斜率极陡的强线性同相轴进入 f-k 域后成为一条独立亮线如果它对应的相速度落在你设的搜索窗内峰值搜索很容易被它吸走。排查方法把原始记录按 340m/s 画一条理论走时线看空气波是否覆盖了面波时窗。处理上有三层手段先时窗切除空气波到达之前的样点再把搜索窗 vmin 抬到 400m/s 以上如果还残留就在 f-k 谱上用扇形滤波把低速区抹平再提取。三层都做完平台就会消失。5.5 现象基阶和高阶模态混成一片峰值搜索来回切换当覆盖层和半空间速度差较大一阶模态在某个频段能量反超基阶这时谱峰搜索会在相邻几个频率点之间来回跳跃曲线上表现为一小段突然抬升又落回。用全局最大搜索基本无法避免。稳妥处理是改成多峰连续性追踪对每个频率取幅度前两个局部极大值分别作为候选然后用上一频率的峰位预测本频率峰位选距离最近的那个候选。代码上可以在 extract_dispersion 里加一个 prev_k 参数引入 15% 的速度变化容忍度。这个逻辑能有效钉住一条模态代价是参数需要根据谱图微调一般工区跑两遍就能定下合适阈值。6. 用合成记录给提取结果上户口一套五分钟能跑完的验证流程6.1 先造一份已知答案的合成记录在接触野外数据之前我用一套很轻的检查流程给提取程序做体检先构造一条理论频散曲线再反向合成 t-x 记录跑完整提取流程后对比误差。理论速度可以给解析式比如 v(f)320180·exp(−f/20)表示低频端接近深层高速约 500m/s高频端接近浅层低速约 320m/s。合成记录时每个频率成分按对应波数在空间上旋转相位叠加后就是一张带频散特征的炮集。def synthetic_shot(nx24, dx2.0, nt1024, dt0.001): 按理论频散关系合成无噪声面波记录. t np.arange(nt) * dt x np.arange(nx) * dx freqs np.arange(2, 80, 0.5) data np.zeros((nt, nx)) for j, xj in enumerate(x): v 320 180 * np.exp(-freqs / 20) # 理论相速度 kappa freqs / v # 空间频率 phase 2 * np.pi * (np.outer(t, freqs) - xj * kappa) data[:, j] np.sum(np.cos(phase), axis1) return data这段代码故意不加噪声、不加窗目的是看提取流程本身有没有系统偏差。跑 preprocess → compute_fk → extract_dispersion 三步把结果和理论曲线画在同一张图上。6.2 把提取曲线和理论曲线直接对比误差落在哪一段最值得关注对比时用 np.interp 把提取曲线插值到理论曲线相同频率点逐点算相对误差。我会把误差按低频段、中频段、高频段分段打印而不是只给一个均值均值容易被中间表现好的频段掩盖问题。f_test np.arange(5, 50, 5) v_theory 320 180 * np.exp(-f_test / 20) v_extract np.interp(f_test, disp[:, 0], disp[:, 1]) err_pct (v_extract - v_theory) / v_theory * 100 for f, e in zip(f_test, err_pct): print(f{f:5.1f} Hz 误差 {e:6.2f}%)正常情况低频段误差应小于 1%高频段因波数分辨率限制可能出现 3%5% 的抖动。如果某段误差达到 10% 以上先检查该频段是否接近空间混叠边界再检查搜索窗口是否卡在高阶分支上。这套合成验证是我每换一个工区、每改一次参数就重跑一遍的保留项目它能把提取程序有没有写对和野外数据能不能提出来两个问题分开——程序的问题不该甩锅给野外数据。我的习惯是野外数据正式处理前先跑合成确认无误后才上手处理完成后把提取结果和相邻钻孔实测的波速分层对一遍速度趋势高频段曲线能否对上浅层低速带是最后一道照妖镜。希望这套思路对你也有用。本文还有配套的精品资源点击获取