企业官网建设流程全解析

简介：多道面波分析是地震勘探中研究浅地表结构的重要手段，相移法频散曲线提取是该领域常用的一类技术流程。面向地球物理专业学生与科研人员，这套MATLAB程序用于从多道面波记录中提取相速度-频率关系曲线，进而推断地层速度结构与物性参数。压缩包体积仅119KB，共包含2个文件：一个MATLAB主程序脚本与一个示例地震数据文件，用户可直接运行并对照输出结果理解每一步处理逻辑。目前已有774人学习/下载。程序实现涵盖数据预处理、快速傅里叶变换相位提取、相邻道相位差计算、频散曲线构建等关键环节；示例数据让初学者无需额外采集资料即可上手，同时也可将脚本改造后用于实际勘探数据的批量处理。借助MATLAB的fft、unwrap等内置函数，程序兼顾教学演示与基础科研需求，对理解面波频散特性和开展浅层速度结构反演具有实用参考价值。

1. 多道面波分析相移法：频散曲线提取的MATLAB实现核心

多道面波分析（MASW）里，相移法（Phase Shift）是提取频散曲线最稳的一条路。很多做岩土勘察、场地波速测试的同行第一次跑通频散提取时，遇到的第一个坎不是数学，而是“为什么我的频散谱一团黑，看不出能量团”。这个问题的答案往往不在算法，而在数据组织和参数设定。这篇笔记围绕一套完整的MATLAB相移法程序，把从原始炮集到频散曲线输出的全过程拆开讲：它解决的是“从几十道面波记录里稳定地抽出相速度-频率关系”这个具体问题，适合正在写论文、做场地调查或者自己搭处理流程的物探和地震工程从业者。先说一个反直觉的结论：相移法并不需要你手动估计视速度范围，它靠二维傅里叶变换把“速度扫描”这件事变成了频谱峰值搜索，所以它的成败更多取决于道间距、时间窗和频率截断。

2. 把炮集变成频散谱：相移法原理与最小可复现流程

2.1 相移法为什么比 τ-p 变换更抗噪

频散曲线提取有两条主流路线：一条是 τ-p 变换，把时间-偏移距记录变换到慢度域，再沿慢度轴找峰值；另一条就是相移法。相移法的核心思想是利用面波在频率域中的相位信息。对原始记录做二维傅里叶变换后，面波信号的能量会在频率-波数域中沿着某条脊线分布，这条脊线的位置对应不同频率下的相速度。你不必预先给定速度搜索网格，因为傅里叶变换天然把“不同速度”映射成了“不同斜率”的能量分布，谱峰就是频散曲线。

这条路的抗噪优势来自累加效应。多道记录在波数域里是对所有道做空间傅里叶变换得到的，随机噪声在空间上不相关，叠加后被压低，而面波信号在道上连续传播，相位一致性在波数谱上会形成尖锐的峰值。相比之下，τ-p 变换对道间时差更敏感，道间距不均匀或者浅层低速层发育时容易产生假峰。相移法对道数不足的容忍度也更好，12 道、24 道都能出可用的频散谱，只是分辨率有差别。

2.2 最小可复现流程：从原始炮集到频散谱的五个步骤

拿到一份 MATLAB 相移法程序，我习惯先不碰它的主函数，而是先把它读数据的入口弄清楚。常见做法是程序里包含一个read_seismic.m或直接在main.m里用load读 .mat 或 SEG-Y 格式。下面是这套流程的骨架，对应程序中phase_shift.m的核心逻辑：

function [f, c, D] = phase_shift(data, dt, dx) % 相移法频散曲线提取主函数 % data: 二维数组，行是时间采样点，列是接收道 % dt : 时间采样间隔，单位秒 % dx : 道间距，单位米 [n, ntrace] = size(data); nfft1 = 2^nextpow2(n); % 时间方向补零到2的幂 nfft2 = 2^nextpow2(ntrace); % 空间方向补零到2的幂 % 二维傅里叶变换到频率-波数域 FK = fft2(data, nfft1, nfft2); FK = fftshift(FK, 2); % 只对空间轴平移，让零波数在中心 % 频率轴和波数轴 f = (0:nfft1-1) / (nfft1 * dt); k = (-nfft2/2 : nfft2/2-1) / (nfft2 * dx); % 相速度网格：c = 2*pi*f / k % 只取正频率、正波数象限 fpos = f(1:floor(nfft1/2)); kpos = k(k > 0); [F, K] = meshgrid(fpos, kpos); C = 2 * pi * F ./ K; % 对每个频率沿相速度轴搜索峰值 % D: 频散谱矩阵，行是相速度，列是频率 amp = abs(FK(1:floor(nfft1/2), k > 0)); amp = amp ./ max(amp(:)); % 归一化 % 峰值搜索 nf = length(fpos); nc = 200; % 相速度采样数，可调 cmin = 100; cmax = 600; % 相速度搜索范围，根据工区调整 c = linspace(cmin, cmax, nc); D = zeros(nf, nc); for i = 1:nf % 把波数域振幅映射到相速度域 for j = 1:nc kk = 2 * pi * fpos(i) / c(j); [~, kidx] = min(abs(kpos - kk)); D(i, j) = amp(i, kidx); end end D = D'; end

这段代码做了四件事：补零、二维傅里叶变换、频率-波数坐标构建、相速度域重采样。关键参数是nfft1和nfft2，它们决定频散谱的频率分辨率和波数分辨率。时间方向补零能加密频率轴，空间方向补零能加密波数轴，但补零不会增加真实信息，只是让谱峰看起来更平滑。cmin和cmax设的是相速度搜索范围，这个范围如果偏离实际值太远，峰值会被截在边界上，出来的曲线会在两端贴边。nc是相速度轴采样点数，设 200 已经够用，再多只是让图像更密，不会提高精度。

2.3 数据组织的隐藏要求

程序的第二个入口是数据排列方式。很多人第一次跑通程序后画出的频散谱是“斜线”而非“能量团”，原因通常是数据矩阵的维度搞反了。data要求是时间 × 道数，也就是每一列是一道地震记录。如果你的原始数据是道数 × 时间，一定要先转置。另一个容易忽略的点是dx必须和原始排列的道间距一致。如果野外采集用的是 2 米道距，室内输入写成了 1 米，所有相速度都会被高估一倍，频散曲线整体向上飘。

3. 从频谱峰值到频散曲线：核心函数的逐段拆解

3.1 峰值搜索与频散曲线自动拾取

频散谱矩阵 D 生成之后，下一步是从谱峰里提取频散曲线。程序里通常有两个拾取策略：一是沿频率轴逐列找振幅最大点，二是用局部峰值检测避开多峰干扰。逐列最大值的实现最直接，但缺点是碰到能量分叉或高阶模式时，会突然跳到另一支能量团上，曲线出现“阶梯状跳跃”。

function [f_out, c_out] = pick_dispersion(f, c, D, fmin, fmax) % 自动拾取频散曲线 % f: 频率轴（Hz），c: 相速度轴（m/s），D: 频散谱矩阵 % fmin, fmax: 拾取频率范围 f_idx = find(f >= fmin & f <= fmax); f_out = f(f_idx); c_out = zeros(size(f_out)); for i = 1:length(f_out) % 当前频率列 col = D(:, f_idx(i)); % 平滑后找最大值，避免毛刺 col_s = smooth(col, 7); [~, idx] = max(col_s); c_out(i) = c(idx); end % 中值滤波去掉跳点 c_out = medfilt1(c_out, 3); % 剔除与邻近点速度差过大的异常值 for i = 2:length(c_out)-1 if abs(c_out(i) - c_out(i-1)) > 0.2 * c_out(i-1) c_out(i) = nan; end end f_out(isnan(c_out)) = []; c_out(isnan(c_out)) = []; end

这段代码用smooth对振幅谱做七点平滑，再用medfilt1做中值滤波，两道防线处理高频抖动。最后一层是速度突变检测，相邻频率点相速度变化超过 20% 就标记为异常值并剔除。任何一个单步都可能把有效信号误杀，但三步组合起来基本能把人工拾取的质量追平。使用时要根据数据情况调节阈值：如果频散曲线本身很陡（比如浅层强速度梯度），20% 的突变判定可能把真实拐点当异常删掉，这种情况建议放宽到 30%。

3.2 相速度搜索范围设定的门道

cmin和cmax不是随便填的。相速度下限通常由场地表层介质决定，土层剪切波速低于 100 m/s 的情况很常见，所以下限设 100 是安全的。上限要看最深目标层，如果关心到 30 米深度，波速可能到 800 m/s，上限设 1000 比较合理。搜索范围过大的代价是频散谱被压缩，低幅值信号在显示时几乎不可见。我一般会先跑一次宽范围，看到能量团大致分布后，再收窄到能量团周围，这样谱图分辨率更高，峰值拾取也更稳定。

3.3 时间窗截断对低频端的影响

程序里通常还会有一个切除直达波和面波后到达的折射波的时窗操作。面波数据里最常见的干扰是折射波和反射波晚至能量，它们会在地震记录的后半段形成强振幅，进入傅里叶变换后污染低频段的波数谱。常见做法是先做底波切除，再做时间方向加窗。加窗函数用汉宁窗或余弦窗，窗长要覆盖面波的主要能量到达范围。

% 时间窗截断示例 tmax_wave = 1.2; % 面波主要能量到达时刻，单位秒 ncut = round(tmax_wave / dt); window = ones(n, 1); window(ncut:end) = 0; window = conv(window, hann(51), 'same'); % 平滑过渡 data_w = data .* window;

ncut这个值决定了低频端能保留多少信号。时窗截短，频率分辨率下降，频散谱的低频端会变糊；时窗截长，晚至干扰进入，高频端可能出假峰。没有一个万能值，需要对着地震记录看一眼面波波列的持续时间再定。

4. 参数设置与边界条件：为什么两组参数会得出两张谱

4.1 频率轴尺度与 FFT 补零的权衡

MATLAB 的fft2在时间和空间两个方向都做快速傅里叶变换，补零操作直接决定频散谱的纵横尺寸。时间方向补零到nfft1，相当于在时域数据末尾加了零值，这会加密频率轴采样，但不会改变真实频率分辨率。真实分辨率由原始记录长度决定：df = 1 / (n * dt)。如果原始数据只有 1 秒长，采样率 1000 Hz，那么频率分辨率就是 1 Hz。这意味着在 5 Hz 和 6 Hz 之间的能量团无法被分开，你看到“低频端糊成一团”不是程序bug，而是数据长度本身的物理限制。

空间方向补零的效应类似。dk = 1 / (ntrace * dx)，如果只有 24 道、道距 2 米，波数分辨率是 1/48 per meter。两个相邻模式在波数域里靠得足够近时，空间补零能让视觉上分开一些，但真实可区分度还是受道数和排列长度约束。这是相移法里最容易被误解的部分：补零是对谱的插值，不是对数据的扩展。

4.2 频散谱显示的增益与裁剪

程序画图部分通常会有clim或vmax参数控制色彩范围。默认画法是把谱归一化后做imagesc，这时能量团周围的小幅值背景噪声也会被显示出来。野外观测数据里，背景噪声可以占到最大振幅的 20%~30%，不裁剪就画，整个谱面会布满噪点，峰值反而被淹没。一般做法是取振幅的 70%~80% 分位数作为显示上限，让低于阈值的部分显示为深色背景，能量团才“浮出来”。这属于显示层面的操作，不会改变计算结果，但会直接影响你对频散曲线可靠性的判断。

4.3 道间距与排列长度的真实约束

MASW 现场采集里，道间距的选择决定了频散曲线的频率上限。空间采样定理要求排列长度内至少有两个采样点覆盖一个波长，所以最高可提取频率大致是fmax = vmin / (2 * dx)。用 1 米道间距理论上能提取到更高频率，但实际受地表衰减和震源频谱限制，高频端信号本来就弱。dx输入不对，整个频散谱的横轴（相速度）会出现系统性偏移，这个错误排查起来很隐蔽，因为谱的形状看起来完全正常。

5. 相移法避坑指南：三个让结果变形的隐蔽问题

5.1 阶段式相位折叠导致的高频振荡

现象：频散曲线在 20 Hz 以上出现锯齿状抖动，不是平滑下降。原因：相移法计算相速度时用的公式是c = 2πf / k，当频率升高、波数变大，波数轴上的采样间隔dk不够密时，会出现相位折叠。相位折叠本质是空间采样不足导致的波数混叠，对应的最高无混叠频率由道间距决定。解决：先检查dx是否和道间距一致，然后确认目标频率上限是否超过vmin/(2*dx)。如果超过了，要么减小道间距，要么在拾取阶段把高频段的曲线做中值滤波。程序里已经写了medfilt1，但注意它默认窗口为 3，频率点很密时窗口太小，压不住振荡，可以改成 5 或 7。

5.2 高阶模式能量盖过基阶模式

现象：频散谱上同时出现两条能量团，自动拾取在某个频率点突然跳到另一条上，曲线出现不连续台阶。原因：高阶模式（如 Love 波或高阶 Rayleigh 波）在某些频段的能量可能强于基阶，自动最大值搜索会锁定高阶模式。解决：把拾取函数里的max改成局部峰值搜索，限定在基阶模式所在的波数范围内。另一种做法是先看谱图，人工确定基阶模式的大致频率-速度范围，然后用这个范围约束自动拾取。程序里如果在pick_dispersion中加入了cmin和cmax输入参数，就可以直接限定搜索区域。高阶模式不是错误数据，它反而携带浅层信息，但如果你要的是基阶曲线，就必须在拾取阶段把它剔除。

5.3 时间窗太短导致频散谱低频端变形

现象：频散曲线低频端（3~8 Hz）出现不正常的“上翘”或者整体偏离实测值。原因：时间窗截取过短，频域分辨率下降，低频段能量被平滑摊开，峰值位置偏移。时间窗对应的频率分辨率是df = 1/T，T 是窗长。窗长 1 秒时低频分辨率 1 Hz，8 Hz 附近的峰值会被 7.5~8.5 Hz 的邻频能量“拉偏”。解决：放宽时间窗，至少保证窗长是目标最低频率的 2~3 倍周期。如果你关心 5 Hz 的信号，窗长至少要有 0.4~0.6 秒的完整面波信号。如果原始记录里面波波列本身就短，低频信息就先天不足，再处理也补不回来，这时候只能换低频震源或者加长排列。

6. 频散曲线能不能信：正演检验与成果输出技巧

频散曲线提取完了，下一步不是直接拿去反演，而是先做一次“能不能信”的检验。我的习惯是抽出曲线后用频散公式正演回理论相速度曲线做对比。MATLAB 里可以用简化的两层或多层模型跑一个快速正演，把计算出的理论频散曲线和拾取曲线叠在一张图里。偏差在 5% 以内说明提取流程没问题；偏差如果集中在某个频段，就要回头检查频谱峰值搜索是否在那一频段被噪声带偏。这个方法比肉眼盯着频散谱看半天有效得多。

成果输出方面，程序一般会保存c_out和f_out到.mat文件或直接写.txt。我建议输出两列文本，第一列是频率，第二列是相速度，分隔符用空格或逗号都可以，后续用 Python 或 GMT 画图都方便。如果是批量处理，把所有测点的频散曲线写进同一个文件，每两条曲线之间空一行，方便后期反演软件读取。

% 输出频散曲线到文本文件 out = [f_out(:), c_out(:)]; writematrix(out, 'dispersion_curve.txt', 'Delimiter', 'space');

最后再提一个平时容易忽略的操作：把提取出的频散曲线和原始频散谱同时画出来。曲线叠在谱上，如果曲线边缘贴合能量团中心，说明拾取正确；如果曲线在能量团上方或下方飘着，说明拾取的频率-速度映射出了问题，最常见的原因就是在波数插值那一步弄错了F和K的网格方向。从那以后我每次跑完相移法，都强制自己把拾取曲线和频散谱叠图检查一遍，确保曲线不是画在了噪声背景里，而是稳地压在能量脊上。这一步花不了两分钟，但能挡掉绝大多数“看着像模像样、实际全是错的”结果。希望帮到你。

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