企业官网建设流程全解析

简介：一份关于频率—波数域频散曲线提取方法及程序设计的学术论文，面向地震勘探、工程物探及瑞利波数据处理方向的研究人员与专业技术人员。资源包共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=ω/k=f/κ（κ 是循环空间频率，单位 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 相速度与波数的换算：v=f/κ 这条公式的适用边界

用 numpy 的 fftfreq 生成坐标轴时，时间轴的单位是 Hz，空间轴的单位是 cycles/m，两者直接相处 v=f/κ 就得到相速度，不需要再乘 2π。这个简洁关系成立的前提是：介质水平分层、排列方向与传播方向一致、检波点在直线上等间隔布设。野外如果测线沿地形起伏布设，或者排列中间有几道被障碍物挪位，波数轴的几何关系会被破坏，谱峰位置偏移，提取出的速度值就不准。常见做法是在预处理里先做道头校正或重排，把偏移量折回规则网格。

另一个隐含条件是远场平面波近似。面波在近炮检距处包含近场成分，波前面不是平面，二维傅里叶变换的"平面波分解"语义会打折扣；工程实践里一般把炮检距最近的道控制在 5～10 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_start=0.0, win_end=None): """对原始炮记录做预处理. data: (nt, nx) 二维数组, 第一维时间, 第二维道号 dt: 时间采样间隔(秒) win_start, win_end: 保留的时间窗(秒), 一般圈面波到达段 """ nt, nx = data.shape # 1) 去常数和线性趋势, 消除直流漂移与超低频扰动 data = detrend(data, axis=0, type='constant') data = detrend(data, axis=0, type='linear') # 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, axis=0)) + 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，这样后面 v=f/κ 直接换速度。

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, d=dt)) # 空间频率轴: 范围[-1/(2dx), 1/(2dx)], 单位cycles/m kappa = np.fft.fftshift(np.fft.fftfreq(nx, d=dx)) 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]，其中负半轴对应反方向传播的能量，拾取时可以只取 kappa>0 的半边。

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 一般取谱上信噪比高的频带，浅层面波通常给 3～80 Hz；vmin/vmax 先按工区经验给宽范围，比如 120～1200 m/s，跑完看结果再收紧，防止拾取线在低频端跳跃。对高频段，vmin 还要受空间混叠限制，具体见第 4.1 节。输出用 np.savetxt() 存成三列文本即可，常见做法是把结果写成 (f, v) 两列，方便直接导入反演程序。

4. 参数怎么设：采样间隔、道间距和排列长度对提取结果的硬约束

f-k 域频散提取不是"数据丢进去就出结果"的黑匣子，谱的质量由几个观测参数直接决定。这些参数多数在野外采集时就已固定，室内只能通过截窗、补零来补救。下面这张总结表先给出关系，再逐条展开。

4.1 道间距 dx 决定最大可分辨波数，限制的是最小相速度

空间采样同样有奈奎斯特限：波数最高只能分辨到 κ_nyq=1/(2dx) cycles/m，超过这个值的能量会折叠回低波数区域，表现成高频段谱峰突然"兜回来"形成假分支。把 κ_nyq 代进相速度公式，得到可分辨的最小相速度 v_min=2f·dx。以 dx=2m、频率 50Hz 为例，v_min=200m/s；如果浅层表土剪切波速度只有 150m/s，那 50Hz 以上的基阶相速度低于 200m/s，必然发生混叠。

这也解释了一个常见现象：某些工区用 2m 道距采的面波记录，高频段提出来速度值一路走高，看起来像速度倒转，其实是混叠造成的假象。室内处理能做的有效手段是：把采集道距不满足 v_min 约束的高频部分直接截掉，或者按相邻两道求和合并成 4m 等效道距（提高道距会加重混叠，必须配合低通滤波），再在拾取结果里删除对应频段。真正要从根源上解决，只能在采集设计阶段按目标最浅层速度反算道间距。

4.2 记录长度和时窗决定频率分辨率

频率分辨率由记录时长决定：Δf=1/T。2 秒记录对应 0.5Hz 分辨率，对 5Hz 以下的低频段来说，谱峰的半宽度可能占掉好几十个百分比，低频端频散曲线看起来就像被低通滤波过一样。野外面波记录经常只记录了 1～2 秒，低频信息不是被仪器滤掉，而是被时间窗截断了。

室内处理常见的补救办法是零填充：把时窗截取后的数据在时间方向补零到 4 倍长度，再做二维 FFT。补零不增加真实信息，但能让频率轴加密，谱峰位置更平滑，对自动拾取有帮助。需要注意的是，如果截取面波到时窗本身只有 0.5 秒，补零后频率分辨率仍由 0.5 秒决定，只是插值点变密，真正有效的低频上限大约在 2Hz 左右。此外，时间采样率 dt 决定 f_max=1/(2dt)，工程仪器普遍用 0.5ms 或 1ms 采样，对 100Hz 以内的面波完全够用，一般不是瓶颈。

4.3 速度搜索范围与模态识别策略

速度窗的设置直接影响拾取结果在哪条模态上落脚。常见做法是先用宽窗跑一遍，得到一条"毛刺很多"的曲线，再根据曲线的大致走向把窗收紧：比如基阶曲线在 5～60Hz 内从 500m/s 降到 250m/s，就设 vmin=200、vmax=600，重新提取一遍。收紧后谱峰不会跳到高阶分支，曲线连续性明显变好。

如果谱上明显存在两条靠近的脊线，可以用两套搜索窗分别提取基阶和二阶模态。还有一种更稳的方法：对振幅谱做对数压缩后，在每个频点搜索局部极大值（比如取前三个峰），再按"相邻频率峰值连续性"原则做路径追踪，而不是每频率独立取全局最大。这个思路在模态混叠时很管用，但代码量比全局搜索多出一截；实际工程里多数数据用速度窗+全局峰就够了，只在工区存在明显高速夹层时才需要多模态追踪。

5. f-k 频散提取的常见翻车与排查：五条真实踩坑记录

5.1 现象：f-k 谱被水平亮条纹盖住，脊线看不清

谱图上出现贯穿所有波数的水平亮条纹，本质是时间方向上的强能量截断或周期性扰动：要么是记录存在大幅度直流漂移，要么是面波时窗外还有强噪声，要么是某道放大器饱和产生方波状记录。时间窗边缘的硬截断也会在谱上形成垂直于频率轴的旁瓣条纹。

排查顺序：先看原始道集的 wiggle 图，确认是否有坏道和饱和道；再做 detrend 和时窗平滑，把窗边缘的汉宁长度从 5% 加到 15%；最后逐道检查 RMS 异常大的道，直接置零或剔除。这三个步骤做完，条带通常消失。还有一种情况是 50Hz 工频干扰，会在固定频率位置形成亮点，和我们要的脊线叠加，处理时用陷波滤波单独切除。

5.2 现象：频率越高能量峰越向右歪，或者某段频率后峰位突然跳回低波数

这是空间混叠的典型表现。波数超过奈奎斯特限后折叠回低波数侧，谱峰看起来"猝不及防"地出现在不该出现的位置。判断方法很简单：把提取结果里的 (f, v) 和 v_min=2f·dx 画在同一张图上，凡是落到该曲线之下（速度更低）的点，基本可以判定为混叠产物。

解决方法是把超过混叠界限的频率点从结果里删掉。如果非要保住高频段，只能从采集端改：缩小道间距 dx 或采用不等距排列。有些工区用地表低速带测出浅层速度只有 130m/s，用 2m 道距在 40Hz 以上就没法采了，这类数据要直接说明高频截止频率，而不是硬提。

5.3 现象：低频段峰值左右乱跳，提取曲线不连续

低频段波数小，谱峰离开原点不远，空间孔径 L=nx·dx 如果不够大，波数分辨率 Δκ=1/L 就比峰位间距还粗，这时相邻几个频率的峰一会落在 0.01cycles/m，一会落在 0.03cycles/m，折线很难看。另一个原因是低频能量本身弱，时窗截断后低频段信噪比下降。

解决思路分两步：先补零到 2 倍或 4 倍长度，让谱在波数方向更平滑；再做 3～5 点频率方向中值滤波，把孤立跳点压掉。如果低频段必须要到 3Hz，而排列只有 24 道×2m=48m，物理上就不足以分辨这个波数，滤波只是让曲线好看，信息其实已经丢了。这种情况不如主动把低频截止提到 5Hz，保证提出来的是可靠数据。

5.4 现象：提取出的频散曲线低速段一直压在某个固定速度上

常发生在 330～350m/s 附近，这条"平台"多半是空气波。空气波在时距图上是一条斜率极陡的强线性同相轴，进入 f-k 域后成为一条独立亮线，如果它对应的相速度落在你设的搜索窗内，峰值搜索很容易被它吸走。

排查方法：把原始记录按 340m/s 画一条理论走时线，看空气波是否覆盖了面波时窗。处理上有三层手段：先时窗切除空气波到达之前的样点；再把搜索窗 vmin 抬到 400m/s 以上；如果还残留，就在 f-k 谱上用扇形滤波把低速区抹平再提取。三层都做完，平台就会消失。

5.5 现象：基阶和高阶模态混成一片，峰值搜索来回切换

当覆盖层和半空间速度差较大，一阶模态在某个频段能量反超基阶，这时谱峰搜索会在相邻几个频率点之间来回跳跃，曲线上表现为一小段突然"抬升"又"落回"。用全局最大搜索基本无法避免。

稳妥处理是改成"多峰+连续性追踪"：对每个频率取幅度前两个局部极大值，分别作为候选，然后用上一频率的峰位预测本频率峰位，选距离最近的那个候选。代码上可以在 extract_dispersion 里加一个 prev_k 参数，引入 15% 的速度变化容忍度。这个逻辑能有效钉住一条模态，代价是参数需要根据谱图微调，一般工区跑两遍就能定下合适阈值。

6. 用合成记录给提取结果上户口：一套五分钟能跑完的验证流程

6.1 先造一份已知答案的合成记录

在接触野外数据之前，我用一套很轻的检查流程给提取程序做"体检"：先构造一条理论频散曲线，再反向合成 t-x 记录，跑完整提取流程后对比误差。理论速度可以给解析式，比如 v(f)=320+180·exp(−f/20)，表示低频端接近深层高速约 500m/s，高频端接近浅层低速约 320m/s。合成记录时每个频率成分按对应波数在空间上旋转相位，叠加后就是一张带频散特征的炮集。

def synthetic_shot(nx=24, dx=2.0, nt=1024, dt=0.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), axis=1) 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% 以上，先检查该频段是否接近空间混叠边界，再检查搜索窗口是否卡在高阶分支上。这套合成验证是我每换一个工区、每改一次参数就重跑一遍的保留项目，它能把"提取程序有没有写对"和"野外数据能不能提出来"两个问题分开——程序的问题不该甩锅给野外数据。

我的习惯是：野外数据正式处理前先跑合成，确认无误后才上手；处理完成后，把提取结果和相邻钻孔实测的波速分层对一遍速度趋势，高频段曲线能否对上浅层低速带，是最后一道照妖镜。希望这套思路对你也有用。

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