简介：这是一套基于MATLAB的被动源面波频散曲线反演程序，主要面向地球物理、地震学与地质工程领域的研究人员和学习者。程序利用天然或人为地震信号提取面波频散信息，进而反演地下分层结构与物性参数，相比主动源方法成本更低、环境限制更少，但数据处理与非线性反演求解也更具挑战。压缩包内包含十七个M脚本文件和一个Markdown说明文档，总大小约二十四KB，涵盖了数据预处理、频散曲线计算、多种反演算法（如Muller方法、快速Haskell算法、介质模型构建等）以及结果可视化模块，目录结构清晰，便于修改和二次开发。目前已有七百一十五人下载学习。借助这套程序，读者可以理解Rayleigh波频散曲线与层状介质模型之间的关系，掌握从信号到反演结果的完整流程，并通过运行示例脚本快速上手，适合作为科研实验、教学演示或工程评估的实用参考。

1. 被动源面波频散曲线反演：这份 MATLAB 程序包解决了什么

接触过主动源面波勘探的同行都知道，敲击或振源激发虽然可控，但遇到城市环境、深厚覆盖层或者需要探测 50 米以下结构时，主动源往往力不从心。被动源面波（环境噪声面波）不需要人工震源，直接利用车辆、工厂、风等背景噪声激发的面波，低频成分丰富，探测深度天然占优——但代价是反演不再是回车就能出结果：频散曲线怎么提、初始模型怎么给、反演算法会不会陷进局部极值，每一步都可能让结果变成「看起来合理，实际不能用」的废图。这份基于 MATLAB 的被动源面波频散曲线反演程序，把从频散曲线到横波速度剖面的完整流程封装成了可改、可跑、可查的脚本集合，适合正在做场地波速测试、覆盖层调查或微动勘探试验的从业者，也适合课题组里刚接手面波数据处理的新手。它不是黑匣子，每个关键参数都摆在台面上，下面我从原理到踩坑逐个拆。

2. 反演先懂原理：从频散曲线到横波速度结构的映射关系

2.1 面波频散：频率、相速度与介质深度之间的耦合逻辑

在层状介质中，瑞利面波的传播速度不是常数，而是与频率有关的函数。高频成分波长短，主要携带浅层信息；低频成分波长长，能穿透更深。这个物理事实直接决定了反演的本质：从一条实测频散曲线出发，反向推算出各层的横波速度、层厚度和纵波速度与密度参数，使其正演出的理论频散曲线尽可能逼近实测曲线。

被动源面波反演的输入一般有两种情况：一种是直接从台阵记录提取出频散曲线（频率-相速度对），输入给反演程序；另一种是输入原始时间序列，程序内部先做频谱分析和频散提取。这套程序采用的是第一种模式——把频散提取和反演解耦开，好处很明显：你可以用任意一种方法（SPAC、FK 或时频分析）先提取频散曲线，单独验证曲线质量，再进入反演环节，避免把提取阶段的误差传导到反演结果里。

频散曲线的质量直接决定反演可信度。实测曲线如果频点稀疏、相邻频点的相速度跳变超过 5%～8%，反演算法会把这种抖动解释成层速度的剧烈变化，产出虚假分层，所以拿到程序后第一件事不是跑反演，而是审视输入的频散曲线是否平滑可用。

2.2 反演参数与模型分层：横波速度为什么是主角

层状模型的参数包括每层的横波速度 Vs、纵波速度 Vp、密度 ρ 和层厚度 h。对于浅层工程勘探，横波速度是场地分类和抗震设计最敏感的指标，纵波速度和密度对频散曲线的敏感度相对低，所以反演中一般固定 Vp 和 ρ，只让 Vs 与 h 参与搜索。程序中 Vp 通常由经验公式从 Vs 推算：

Vp = Vs * sqrt((1-ν) / (0.5-ν))

其中 ν 是泊松比。注意这个公式中 ν 取值对 Vp 影响很大，饱和土层取 0.45～0.48，干燥砂土取 0.25～0.35，取值错误会让理论频散曲线高频段明显偏离实测，出现「反演拟合很好但波速整体偏大」的假象。

层数和层厚的设定遵循从粗到细的原则。覆盖层较厚时，浅层厚度取得过小会需要大量层数支撑，计算量成倍增加；浅层厚度取太大又会让高频信息无法约束。常见的做法是最浅层厚度不超过最短波长对应探测深度的三分之一，层厚度按等比数列递增，相邻层厚度比值在 1.2～1.5 之间。

2.3 正演与反演的闭环：目标函数与搜索策略

反演的本质是最小化实测频散曲线与理论频散曲线的残差。目标函数一般写成：

Φ = Σ [ (V_obs(fi) - V_cal(fi, m)) / σi ]²

m 是待反演的模型参数向量（各层 Vs 与 h），σi 是该频点的误差估计。如果 σi 偏小，说明该频点权重高，反演会更努力拟合它；相反则放松约束。实践中常见的问题是所有频点设成同一 σ，导致高频段噪声被过度拟合。

搜索策略上，这类程序大多采用群体智能算法，遗传算法和粒子群算法是主流。粒子群实现简单、收敛快，但容易早熟；遗传算法全局搜索能力强但参数（种群数量、交叉率、变异率）需要调。程序默认采用粒子群算法时，建议把种群数设在 60～100 之间，迭代次数 200～500 次，惯性权重从 0.9 线性降到 0.4——这个设置对大多数覆盖层场地都能兼顾收敛速度与全局搜索能力。

3. 程序结构拆解：从数据读入到剖面输出的完整链路

3.1 主脚本与核心函数模块划分

拿到程序包后先别急着双击 main，我一般用 MATLAB 直接浏览整个目录结构，把文件分成四个功能群，再按调用顺序走一遍代码。正常的调用次序是：读入频散曲线文件 → 设定模型层数与参数边界 → 调用正演函数 → 运行反演算法 → 绘图输出。

主脚本通常是一个将上述步骤串起来的编排层，它不承载算法实现，只负责参数配置与流程控制。如果主脚本里混杂了大段矩阵运算，反而不利于排查问题。程序包内正演函数和反演函数的接口设计决定了你能否替换算法，好的接口会统一成 function [residual, fitted] = inverseCore(para, dispData)，para 是模型参数向量，dispData 是实测频散曲线结构体。

3.2 正演函数：理论频散曲线怎么算

正演是反演的内核，它必须足够快才能在多次迭代中反复调用。层状介质中瑞利波频散曲线的计算常用互易定理或快速标量传递矩阵法，核心代码逻辑如下：

function [phaseV, groupV] = forwardModel(model, freq) % 输入: model 为结构体, 含每层厚度 h, 横波速度 vs, 纵波速度 vp, 密度 rho % freq 为频率向量 (Hz) % 输出: phaseV 相速度, groupV 群速度 nLayer = length(model.h); phaseV = zeros(size(freq)); for iF = 1:length(freq) w = 2*pi*freq(iF); kn = w / 200; % 初始猜测波数 % 使用二分法求频散方程的根: 波数满足瑞利特征行列式 D(k)=0 [pV, err] = bisectionRoot(@(k) rayleighDet(k, w, model), ...); end end

代码里最耗时的部分是瑞利特征行列式的求根，正演函数效率决定了整个反演耗时。如果单次正演耗时超过 0.5 秒，500 次迭代 × 80 个粒子的计算量会达到数小时，这是不可接受的。遇到这种情况，先检查是否有矩阵重复计算没有缓存，再考虑减少层数。

3.3 反演引擎：粒子群参数与种群初始化

反演引擎负责生成初始种群并按目标函数迭代更新，粒子群部分的核心逻辑通常如下：

for iter = 1:maxIter w = 0.9 - (0.9 - 0.4) * iter / maxIter; for i = 1:nPop % 速度更新公式: velocity(i,:) = w * velocity(i,:) + ... c1 * rand * (pbest(i,:) - pos(i,:)) + ... c2 * rand * (gbest(:) - pos(i,:)); % 越界处理: 参数超出边界则拉回并反向飞 pos(i,:) = pos(i,:) + velocity(i,:); pos(i,:) = max(lo, min(hi, pos(i,:))); % 目标函数计算 [res, ~] = forwardModel(struct('h', hset, 'vs', pos(i,1:nLayer)), freq); end end

注意速度更新公式中 c1、c2（自我认知系数与社会认知系数）的取值通常都设 1.5 左右，两者之和大于 4 时会出现震荡不收敛。边界约束采用拉回策略，速度快但可能造成粒子在边界堆积，更稳的办法是反弹策略（超出边界时速度取反）。如果反演出的模型参数大量落在边界上，说明边界设置过窄或模型层数不适合当前频散曲线。

3.4 结果输出与绘图工具

输出部分一般包含三个图：实测频散曲线与理论频散曲线对比图、横波速度随深度变化剖面图、进化曲线（目标函数值随迭代下降）。看进化曲线能快速判断是否收敛：曲线在 30 代内就平坦，说明种群多样性丢失，早熟风险大；曲线持续下降且最后趋于平稳，则收敛状态正常。剖面图上应标注每一层的速度值与层底深度，方便直接读取覆盖层厚度与等效剪切波速。

4. 运行与参数调试：把自己的数据在程序上跑通

4.1 输入数据准备：格式与预处理

程序读入的频散曲线文件一般要求两列，第一列是频率，第二列是相速度，单位需要统一检查。常见翻车场景是数据文件里混入了中英文符号的逗号、表头文字或者 Nano 标识符，MATLAB 在 textscan 读入时会报错或读入错位数据。

fid = fopen('dispersion_curve.txt', 'r'); data = textscan(fid, '%f %f', 'CommentStyle', '#'); fclose(fid); freq = data{1}; phaseV = data{2}; % 剔除不合理的频点: 相速度小于 0 或大于预设范围的值 validIdx = phaseV > 50 & phaseV < 2000; freq = freq(validIdx); phaseV = phaseV(validIdx);

数据处理中一个容易被忽略的环节是频带裁剪。实测频散曲线在低频端如果信噪比低、相速度误差大，应该手动剔除这些频点，而不是全部交给反演加权。高频端同理，尤其是记录长度有限时，高频信息往往受噪声干扰严重，一般保留信噪比大于 2 的频率段。

4.2 反演参数配置表：关键参数的取值范围与调试建议

参数配置集中在程序头部，下面的表格总结了各参数在典型浅层工程勘探场景下的推荐值，以及偏离推荐值时的表现，建议直接抄作业后再按实测数据微调。

需要特别提醒的是，横波速度边界不能都设成一个固定范围，浅层 Vs 可能在 150～400 m/s，深层可能到 600～1500 m/s，需要根据场地经验逐层设定。某次在覆盖层场地反演时，因为第一层边界设得过高（最低 250 m/s），低频段拟合尚可，但高频段理论曲线明显偏离浅层低速特征，最后把首层下限降到 120 m/s，拟合残差瞬间降了 40%。

4.3 运行结果检验：反演收敛、曲线拟合与地质合理性

反演程序跑完不代表工作结束。拿到剖面后，第一件事是看理论频散曲线与实测数据的拟合残差是否在每个频段都均匀分布——如果在部分频段系统性偏低或偏高，说明模型结构与实际情况有偏差，需要通过增加层数或调整边界解决。

其次，把反演得到的 Vs 剖面与该区域已有钻孔资料或主动源面波结果交叉验证。反演是典型的病态问题，等效性意味着不同模型可以产生几乎一致的频散曲线，所以没有任何约束的反演结果一定要谨慎解读。如果程序支持多初始模型反演，我一般跑 3～5 次并叠加结果，各次结果的速度差异在 10% 以内，说明模型可信度高；超过 20% 则说明约束不足，需要根据地质前提收紧边界。

5. 反演实战避坑指南：六个高频翻车场景与排查思路

5.1 高频段拟合正常但低频段曲线抬不起来

现象：反演完成后，理论频散曲线在高频段与实测贴合很好，低频段却整体偏移或呈锯齿状，频率低于 5 Hz 时相速度明显偏高，且无论怎么调参数都无法下降。

原因：低频信息携带深层介质特征，如果初始模型深层横波速度边界设得太低，反演搜索空间根本不包含真正深层高速结构，再迭代也无法触及。另一种可能是实测频散曲线在低频段本身存在大幅误差，低频端受台阵尺寸限制，提取的相速度已不可靠。

解决：先检查实测曲线低频段的提取质量：用 SPAC 法时台阵的最长半径决定了有效最低频率，一般要求最长半径对应的波长不小于感兴趣最深结构的四倍。然后按场地经验和区域资料重新设定深层 Vs 上限。我遇到一次翻车就是因为把第四层 Vs 上限设成 800 m/s，实际该地区第四层是风化岩，速度 1200 m/s 以上，改完边界后低频段立刻拟合上了。

5.2 反演结果整体平滑但浅层有一串交替高低速薄层

现象：剖面图上浅层 5 米内出现间隔 0.5 米左右的高速度与低速度交替薄层，从频散曲线看拟合残差不大，但地质解释上完全不成立。

原因：反演目标函数只考虑频散残差，没有加入模型平滑约束。高频段对浅层速度变化极其敏感，实测曲线的微小误差（相速度 5 m/s 级别）就能让反演用薄层结构去补偿噪声。这种现象在频点间距小于 0.5 Hz 的高密度频散数据中尤其常见。

解决：程序若支持正则化，在目标函数中加入相邻层速度差惩罚项：

Φ_new = Φ_original + λ * Σ (Vs(i+1) - Vs(i))²

λ 取值从 0.01 开始试，逐步增大到模型不再出现交替薄层。此外限制最小层厚，比如规定最浅层厚度不小于 1 米，也能从物理上排除虚假结构。需要强调，正则化强度不要过度，否则会把真实的速度界面也抹平。

5.3 目标函数值收敛但最终速度剖面每次运行结果差异很大

现象：同一份频散曲线，多次运行时每次跳出略微不同的收敛值，得到的 Vs 剖面在深层显著不同，有的解释为第四层高速层，有的解释为速度过渡带，无法确定哪个可信。

原因：反演的非唯一性没有通过多起点搜索覆盖。粒子群算法单次运行只能探索搜索空间的一部分，尤其搜索空间维度高（4 层以上）、参数相关性强的场景中，不同初始种群极易落进不同的局部极值。

解决：在程序外部增加循环，用不同随机种子跑 5～10 次。

for trial = 1:10 rng(trial * 100 + 7); [optModel, bestRes] = runInversion(dispData, config); finalModels(trial,:) = optModel; end % 统计每次反演结果的层速度均值和标准差, 剔除偏离均值超过两倍标准差的结果

统计后的模型中，各层速度标准差小于 10% 的结果才值得进入解释流程。这算是我个人的一个习惯，从那以后每次做被动源数据处理都强制把多起始反演写进流程，宁可多花一个小时跑批，也不愿意在解释会上被问到「你这个结果唯一吗」时语塞。

5.4 程序运行耗时过长，单次迭代要好几秒

现象：MATLAB 启动后运行主脚本，性能监视器显示正演函数被反复调用，但单次正演计算耗时已超过数秒，按 200 次迭代 × 80 粒子估算，整个反演需要通宵才能出结果。

原因：正演代码中可能包含重复的矩阵求逆或特征值求解。比如在频散方程求根循环里，每次都重复计算各层的传递矩阵，而没有利用同一频率下矩阵的中间参数缓存。另一个常见原因是模型层数设置过多（比如 12 层以上），导致传递矩阵阶数增大和求根耗时的非线性增长。

解决：用 MATLAB Profiler 定位耗时热点函数，优先优化被调用次数最多的矩阵计算。可做三件事：一是把所有不依赖迭代变化的常量矩阵移到循环外一次性计算；二是把频散求根中的高精度容差放宽到 1e-6（从 1e-10 提升）；三是如果层数超过 8 层且频点超过 200 个，考虑是否能用覆盖层等效模型近似。常见做法是把层数控制在 5～7，既抓住主要界面又控制计算量。

5.5 反演报错「矩阵维度不匹配」或「下标越界」

现象：程序运行到中段突然报错，错误信息提示矩阵维度不匹配，或者索引超出数组边界。

原因：大多发生在读入频散数据的长度与初始化模型时的频点数量不一致。少数情况是子函数里固定了数组长度，而实际输入层数与配置的 nLayer 不一致。比如初始化模型参数向量长度为 2*nLayer，但正演函数从 model 结构体中读取的 h 向量长度只有 nLayer-1，程序直接崩溃。

解决：在调用反演前加一层防御性校验：

assert(length(model.h) == nLayer, '层厚度数量与层数不匹配'); assert(length(model.vs) == nLayer, '横波速度数量与层数不匹配');

同时检查频散曲线输入是否包含 NaN 或 Inf，这些值会导致正演函数跳出数值异常错误。把错误定位到数据和模型配置的不一致上，比逐行检查主脚本高效得多。

5.6 反演出的 Vs 剖面与钻孔波速结果系统性偏低

现象：反演剖面整体趋势合理，但各层 Vs 与同一场地实测的钻孔波速相比普遍低 15% 以上，低频段尤其明显。

原因：被动源面波反演中，相速度频散曲线可能是群速度而不是相速度，两者在频散明显的频段差异显著。另一个常见原因是正演代码中 Vp 和密度赋值不合理，如果泊松比设置过低，瑞利波相速度整体会被压低。

解决：先确认输入的频散曲线提取阶段界定的速度类型是相速度而非群速度，观察频散曲线上高频端是否收敛到表层剪切波速附近。然后核对程序中 Vp 和 ρ 的计算公式，特别是泊松比的选择要与实际场地岩土性质匹配。

6. 进阶技巧：给反演加约束条件和稳健性评估

6.1 分层约束与先验信息融合

当有钻孔资料、主动源勘探结果或区域地质资料时，可以把这些先验信息转化为参数边界约束，而不是只依赖目标函数的惩罚项。具体做法是：有钻孔结果的层位深度，把厚度边界收窄到 ±1 米；已知某些层是高速硬层，则将该层 Vs 下限设为已知值的 80%。约束越紧，反演的多解性空间就越小，结果越稳定，代价是如果先验信息本身有误，反演结果也会被带偏。所以在加入强约束前，建议先做一次无约束反演与有约束反演对比，两者结果若偏差超过 15%，说明先验信息可能不可靠。

6.2 正则化参数的自适应调整

传统的固定 λ 需要反复试验才能找到平衡点。更快的方案是先用较大的 λ 快速获得平滑模型，再逐步降低 λ 并用上一轮得到的模型作为下一轮初始模型，这个过程称作冷却法。实际操作中我一般从 λ=1 开始，每轮迭代后检查目标函数残差残差下降率，低于 2% 时就将 λ 减半继续，两到三轮就能定位合适的正则化强度，省去大量手动调参时间。

6.3 灵敏度分析与不确定性输出

对反演结果说服力升级的关键一步是输出灵敏度信息。改变某一层参数（±10%）后重新正演，观察理论频散曲线变化的幅度与该频段实测误差的关系。若某层参数的改变对曲线影响很小（小于实测误差的一半），说明该层参数在当前频带内是不可辨识的，应明确标注给报告阅读者。灵敏度分析代码结构如下：

% 对第 i 层的 vs 做 ±10% 扰动 for iLayer = 1:nLayer modelPert = model; modelPert.vs(iLayer) = model.vs(iLayer) * 1.1; [phaseV_pert, ~] = forwardModel(modelPert, freq); sensitivity(iLayer) = norm(phaseV_pert - phaseV_base) / norm(phaseV_base); end

灵敏度的阈值判断建议设置为：如果所有频率范围内灵敏度系数均值小于 0.05，则标记该层为不可辨识层。在最终的剖面图上，用透明度或虚线表示不可辨识层，这样报告不会给人虚假的精确感。从那以后，我每次交付反演结果时都附上灵敏度表和多次运行稳定性统计，虽然多花时间，但解释风险大幅下降，希望帮到你。

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