企业官网建设流程全解析

简介：超声Lamb波的频散曲线是无损检测中模式识别和参数优化的关键依据。这份Matlab程序包面向从事结构健康监测、超声检测研究的工程师与研究生，基于Rayleigh2Lamb方程分析无限大板中的频散特性，提供相速度与群速度的数值计算思路，并通过编程绘制铝板中的频散曲线，便于直接对照试验结果。压缩包共9个文件，以8个m脚本为核心，另有1个prj工程文件；m脚本覆盖模态分类、扫频计算、曲线绘图等环节，prj用于统一管理工程配置，整体大小仅9KB，轻量且易于二次开发。已有5325人学习浏览，表明该方法在Lamb波检测教学与工程验证中获得较广泛认可。借助完整代码与清晰的文件分工，读者可快速复现频散曲线，开展模态分析、缺陷检测参数选择，或扩展为其他板结构的求解程序。

1. Lamb波频散曲线：方程好抄，根难求

用MATLAB把Lamb波的频散曲线画出来，真正卡人的不是绘图函数，而是“稳定地求根”。频散曲线是相速度cp或群速度cg对频厚积F的曲线，它直接决定导波检测里选什么频率、按什么波速算定位。做过结构健康监测或者薄板探伤的朋友都有印象：低频段S0、A0两条线干干净净，频厚积一过几个MHz·mm，高阶模式挤在一起，漏根、跳线、单位错乱轮着来。这篇文章面向正要搭导波仿真、或者想用频散曲线设计激励参数的研究生和工程师，按“方程定义—求解实现—群速度与验证—避坑—延伸应用”的顺序，把这个MATLAB实现讲透。

2. 频散方程与参数约定：先弄懂对称和反对称，再谈画图

2.1 S模式与A模式：两组方程，一套求解逻辑

Lamb波在自由板中沿板面方向传播时，按照位移关于板中面的对称性分成两类：对称模式记作Sn，反对称模式记作An。它们各自满足Rayleigh-Lamb方程，形式很像，但求解时谁也不能替代谁。

以相速度cp、角频率ω、波数k = ω/cp定义两个辅助量：

p² = ω²(1/cL² - 1/cp²)
q² = ω²(1/cT² - 1/cp²)

对称模式满足：

tan(qd) / tan(pd) = -4k²pq / (q²-k²)²

反对称模式满足：

tan(qd) / tan(pd) = -(q²-k²)² / (4k²pq)

这里d是半板厚。需要提醒的是，当cp小于cT甚至小于cL时，p或q会变成虚数，tan作用在虚宗量上会退化为tanh。很多初版代码在这里直接报NaN，不是方程写错了，而是没有处理复变量。

下方这张表把变量含义和这套资源里推荐的单位制列清楚，后面所有公式和代码都按这个约定走。

mm/us这个单位制的妙处在于数值好看：钢的cL约5.9、cT约3.2，画图时纵轴直接标注km/s，完全不需要换算。频厚积用MHz·mm，数值上等于km/s量纲，也方便和文献对照。

2.2 材料参数换算：从E、ν、ρ得到cL和cT

多数人手里只有杨氏模量、泊松比和密度，写代码第一步就是把它们换算成纵波和横波波速。公式如下：

cL = sqrt(E·(1-ν) / (ρ·(1+ν)·(1-2ν)))
cT = sqrt(E / (2ρ·(1+ν)))

下面这段MATLAB代码以普通钢板参数为例，输出单位为m/s，实际使用时换算成mm/us：

% 材料参数换算：E[GPa], nu[-], rho[kg/m^3] E = 210; % 钢的杨氏模量 nu = 0.3; % 泊松比 rho = 7800; % 密度 cL = sqrt(E*1e9*(1-nu)/(rho*(1+nu)*(1-2*nu))); % m/s cT = sqrt(E*1e9/(2*rho*(1+nu))); % m/s % 转成 mm/us，方便后续与 MHz、mm 搭配 cL = cL * 1e-3; cT = cT * 1e-3; fprintf('cL = %.3f mm/us, cT = %.3f mm/us\n', cL, cT);

这段代码的逻辑是先按国际单位制算出波速，再做一次单位缩放。之所以要在最后统一转成mm/us，是因为频厚积F = f × h在MHz和mm组合下数值适中，残差函数中的πF再乘上速度倒数后，量级接近1，不容易触发数值溢出。参数方面，如果你换成铝或者复合材料，只需要替换E、ν、ρ三行，后续求解函数完全不用改。

2.3 “半板厚”与频厚积：你以为的厚度可能差了一倍

Rayleigh-Lamb方程里的d是半板厚，不是整板厚度。早期版本如果把总厚度h直接当成d代进去，算出来的频散曲线会在频厚积轴上被压缩一半：S1、A1的截止位置全都对不上文献。这个错误很隐蔽，因为曲线形状看起来是合理的，只有和参考图对比时才发现错位。

更稳妥的做法是从一开始就用频厚积F = f × h作为自变量，把半板厚d从公式里消掉。因为pd = πF·sqrt(1/cL² - 1/cp²)正好是一个只依赖F的表达式，具体推导只需要把ω = 2πf、d = h/2代入pd的定义，约分后F自然出现。后面所有求解代码都不再显式出现d，只接收F、cL、cT三个参数，这样既避开了半板厚陷阱，也让求解器可以在不同板厚之间复用。

3. MATLAB求解实现：扫频、二分、去重、跟踪的完整流程

3.1 核心残差函数：用sin/cos形式避开tan奇点

直接按tan形式写残差函数会遇到两类麻烦。第一，tan在角度接近π/2时发散，符号求根法会把奇点误判为过零点；第二，当cp小于cT时q变成虚数，tan函数在复数域的行为不如sin/cos直观。因此，把方程两边同时乘以cos(pd)cos(qd)，改写成sin/cos组合形式：

对称模式残差：
sin(qd)cos(pd)(q²-k²)² + 4k²pq·cos(qd)sin(pd) = 0

反对称模式残差：
4k²pq·sin(qd)cos(pd) + (q²-k²)²·cos(qd)sin(pd) = 0

进一步把ω因子从p、q、k中提取出来，残差函数只保留F和cp两个变量。这样写的好处是输入参数少，且在F范围0.1到20 MHz·mm内不会因为频率量级产生数值问题。

function res = lambResidualNorm(cp, F, cL, cT, isSym) % lambResidualNorm 归一化残差 % cp : 相速度，mm/us % F : 频厚积，MHz*mm，F = 频率(MHz)*板厚(mm) % cL,cT: 纵波/横波波速，mm/us % isSym: true=对称模式，false=反对称模式 p = sqrt((1/cL)^2 - (1/cp)^2); % 可能为虚数，MATLAB复数运算自动处理 q = sqrt((1/cT)^2 - (1/cp)^2); pd = pi * F * p; qd = pi * F * q; sp = sin(pd); cp2 = cos(pd); sq = sin(qd); cq = cos(qd); q2k2 = (1/cT)^2 - 2*(1/cp)^2; % (q^2-k^2)/omega^2 coeff = 4 * (1/cp)^2 * p * q; % 4*k^2*p*q/omega^4 if isSym res = sq .* cp2 .* (q2k2.^2) + coeff .* cq .* sp; else res = coeff .* sq .* cp2 + (q2k2.^2) .* cq .* sp; end res = real(res); % 复数运算的虚部理论上抵消，浮点残余直接去掉 end

这段代码看起来短，但两个设计点值得说明。第一，p和q保持复数形式，MATLAB的sin和cos原生支持复数参数，当cp < cT时自动退化为双曲函数，不需要手工分情况讨论。第二，残差函数只依赖F而不是频率和厚度分开传入，这是上一节“消去半板厚”的落地实现。参数方面，cp传入的是mm/us数值，如果直接传入m/s，残差函数会完全失效，因为1/cp的量级差了1000倍。

3.2 每个频厚积下的求根器：对数网格扫一遍，符号变化定区间

残差函数准备好之后，下一步就是在给定F下找出所有满足方程的cp。常见做法是在cp搜索域内密集采样，找到残差符号变化的区间，再对每个区间调fzero精确收敛。

搜索域需要覆盖全部物理上可能的根。A0模式在低频时相速度趋近于0，所以下界不能太高；高阶模式在截止频率附近相速度趋于无穷，上界也不能只设到cL。实践中常用范围是0.02cT到4cL。低频加一个对数采样，对A0这种低相速度根更友好。

function rootsCps = solveLambAtF(F, cL, cT, isSym) % solveLambAtF 在给定频厚积下求解某类模式的所有相速度根 cpMin = cT * 0.02; cpMax = cL * 4; N = 6000; % 对数网格：低频段更密，照顾A0的小相速度根 cpGrid = logspace(log10(cpMin), log10(cpMax), N); r = arrayfun(@(cp) lambResidualNorm(cp, F, cL, cT, isSym), cpGrid); % 符号变化的网格段即为疑似根区间 signProd = sign(r(1:end-1)) .* sign(r(2:end)); idxSeg = find(signProd < 0); rootsCps = zeros(1, numel(idxSeg)); for i = 1:numel(idxSeg) a = cpGrid(idxSeg(i)); b = cpGrid(idxSeg(i) + 1); try rootsCps(i) = fzero(@(cp) lambResidualNorm(cp, F, cL, cT, isSym), [a b]); catch rootsCps(i) = NaN; end end % 去掉求解失败的项，并按相对容差合并离得过近的根 rootsCps = rootsCps(~isnan(rootsCps)); if ~isempty(rootsCps) rootsCps = uniquetol(rootsCps, 1e-5, 'DataScale', max(rootsCps)); end end

这段代码的运行逻辑是：先在整个物理区间内找符号变化点，再用fzero把根精确出来。之所以要加try-catch，是因为网格段端点处残差可能出现恰好为0或者接近0的情况，fzero的闭区间要求两端严格异号，遇到等号会报错。参数N=6000对常见的01到20 MHz·mm区间足够，如果扫到高阶模式密集区仍发现漏根，把N提到12000，耗时大约增加一倍但稳定性明显更好。

3.3 模式追踪：不要简单排序，要靠连续性“认领”根

每个频率点求出的根是一组无序数。如果直接把每个F点上的根从小到大排列后连线，低频段还能对上，高频段会出现S1和S2互换位置的情况，画出来的曲线在交叉点附近会突然跳线。

解决思路是对每个模式单独跟踪：从上一个频率点的根出发，在当前频率的所有根里找距离最近的那个，认领为本模式的延续。已经被认领的根要从候选列表里剔除，防止两个模式抢同一个根。

function cpMat = traceModes(Fvec, cL, cT, isSym) % traceModes 沿频率轴追踪每一条模式曲线 nF = numel(Fvec); cpMat = NaN(nF, 20); % 预留20列，实际模式数自动扩展 for iF = 1:nF roots = solveLambAtF(Fvec(iF), cL, cT, isSym); if iF == 1 % 第一个频率点，按大小顺序初始化为S0,S1,... roots = sort(roots); nMode = numel(roots); cpMat(iF, 1:nMode) = roots; else prev = cpMat(iF-1, :); prev = prev(~isnan(prev)); for im = 1:numel(prev) if isempty(roots), break; end [~, jMin] = min(abs(roots - prev(im))); cpMat(iF, im) = roots(jMin); roots(jMin) = NaN; % 防止两个模式认领同一个根 end end end cpMat = cpMat(:, ~all(isnan(cpMat), 1)); end

这段代码的关键在于“认领后置NaN”。如果不做这一步，当两个模式的根距离相近时，后面的模式可能把前面已经选走的根再选一次，造成两个模式共用同一个解。参数方面，列数20对常见钢板厚度已经够用，复合材料或者厚板上的高阶模式超过20个时，可以把这个数字改大。

4. 群速度计算与可视化：从相速度到实际传播速度

4.1 数值差分前先平滑：不然群速度曲线全是毛刺

相速度曲线连续求出来后，群速度不能直接抄公式里的解析表达式，而是要用数值差分。群速度和相速度的关系是cg = dω/dk，其中k = ω/cp。如果直接用相邻两个频率点做差分，任何细微的相速度抖动量级会被放大，画出来的曲线在高频段全是毛刺。

% 假设 cpMat 维度为 nF x nM，Fvec 为对应频厚积列向量 Fvec = Fvec(:); omega = 2 * pi * Fvec; % 与厚度无关的比例关系，画图时足够用 kMat = omega ./ cpMat; % 波数矩阵 cgMat = zeros(size(kMat)); for im = 1:size(kMat, 2) kSm = movmean(kMat(:, im), 11); % 先做11点滑动平均 cgMat(:, im) = gradient(omega) ./ gradient(kSm); end % 两端diff噪声大，裁掉前后5个点再画图 Fplot = Fvec(6:end-5); cgPlot = cgMat(6:end-5, :);

逻辑上先平滑k再差分，比先差分再平滑更稳。gradient在均匀网格上计算中心差分，端点本身就不准，所以裁头裁尾是必要操作。参数方面，11点滑动窗口适合F步长小于0.05的密集扫频，如果F步长较大，窗口要相应缩小，否则会把群速度曲线上真实存在的凹陷也抹平。

4.2 绘图与坐标习惯：横轴MHz·mm，纵轴km/s

频散曲线的标准画法是横轴频厚积F、纵轴速度。相速度和群速度可以画在同一张图里，实线表示相速度、虚线表示群速度，同一条模式用同一颜色，这样一眼就能看出频散强弱。

figure('Color', 'w'); hold on; colors = lines(6); linestyles = {'-', '--'}; for im = 1:size(cpMat, 2) plot(Fvec, cpMat(:, im), 'Color', colors(im, :), 'LineWidth', 1.4); plot(Fplot, cgPlot(:, im), '--', 'Color', colors(im, :), 'LineWidth', 1.2); end xlabel('频厚积 F [MHz\cdotmm]'); ylabel('速度 v [km/s]'); legend('S0 cp', 'S0 cg', 'A0 cp', 'A0 cg');

因为mm/us和km/s数值相等，纵轴直接标成km/s。读图时需要特别注意A0模式：它的相速度在F趋近0时也趋近0，群速度曲线却不是单调的，会先升后降，这是弯曲波典型的频散行为。S0模式低频段接近板波速度，大约5.4 km/s，然后逐渐下降。

4.3 极限验证：低频板速和高频瑞利速

代码跑通后第一件事不是看全貌，而是取两个极限点和公式对比。钢板例子中，S0在F趋近0时应收敛到板波速度：

c_plate = 2·cT·sqrt(1 - (cT/cL)²)

高频时所有模式相速度和群速度都收敛到瑞利波速，工程近似式为：

cR ≈ cT·(0.87 + 1.12·ν) / (1 + ν)

用前面钢的材料参数代入，c_plate约5.44 km/s，cR约2.98 km/s。如果代码画出的S0低频值在5.4附近、高频所有曲线尾部收在3.0附近，说明求解和单位制都没有问题。这两处验证价值很高，因为它们是解析值，不依赖任何其他软件。

5. 避坑清单：漏根、跳线、单位混乱、溢出和搜索域过窄

5.1 漏根：符号变化法在“接触零点”处失效

现象：某一阶模式曲线在某个频厚积区间突然缺了一小段，把网格加密后又恢复了，最初以为是fzero收敛失败，实际上根压根没被检测到。

原因：符号变化判据依赖相邻采样点乘积为负。如果网格恰好落在根上，或者两个根离得太近导致网格段内残差同号，符号变化法就漏检了。

解决：两层保险。第一，在网格采样后额外检查abs(r) < 1e-8的网格点，直接把它们也作为候选根；第二，对F方向步长加密，让相邻F点的相速度变化不超过搜索网格间距，保证漏掉的机会变小。我一般把F范围再分成两段扫，低F段用细网格高F段用粗网格，漏根概率能压到很低。

5.2 模式跳线：根不排序直接连线，曲线互相串

现象：单独求每个F的根再画图，S1和A1在某个频段发生明显交叉，两条曲线像麻花一样拧在一起，但对照文献图谱发现真实曲线并不会这样。

原因：求解器返回的根无序，按大小排序后，高频段不同模式相互逼近时顺序会发生交换，连线就串了。

解决：用上一节traceModes里的“认领”逻辑，每次以上一频率点的根为基准找最近根，并且同一个根不允许被认领两次。这个方案比排序稳健得多，代价只是每帧需要多跑一次min查找，对几千个频率点来说耗时完全可以接受。

5.3 单位混乱：m/s配kHz，算出满屏NaN

现象：所有设置看起来都对，但fzero提示端点函数值不异号，或者画出来的曲线数量级差了1000倍。检查残差函数里1/cp，发现cp如果用的是m/s，而cL、cT用的是mm/us，速度倒数直接差出3个数量级。

原因：公式里所有速度必须同一单位制。mm/us加MHz是最舒服的组合，但很多人习惯先把波速换成m/s，频率用Hz，这时如果没同步转换，方程就废了。

解决：在脚本最开头做一次单位断言，比如检查cL是否在1到20之间，因为常见材料的纵波速度在mm/us制下最多是个位数。单位制统一后，不妨再验证一次S0低频板速值，这一步能排查掉90%以上的单位问题。

5.4 溢出现象：高频厚板时残差函数直接爆掉

现象：板厚达到几十毫米、频率到MHz时，数组函数输出Inf或NaN，符号变化判断完全失效，求根器一个根都返回不了。

原因：sin/cos复数形式中，虚部很大时指数增长，double超过约1e308就溢出。归一化后的残差已经比tan形式好很多，但极端频厚积下依然有风险。

解决：把残差整体除一个尺度因子cosh(imag(pd))·cosh(imag(qd))，让sin/cos项变成有界组合。关键代码就一行：

% 在lambResidualNorm返回前加一句防溢出缩放 scale = cosh(imag(pd)) .* cosh(imag(qd)); res = res ./ scale;

这样处理不会改变零点位置，因为除以的是一个恒正的数。高频厚板场景里，这条路比直接改N有效得多。

5.5 搜索域太窄：cp上限定到cL，高次模式全丢

现象：S0、A0、S1、A1都在，往上S2、A2、S3怎么也出不来，或者到高频末尾突然断掉。

原因：高阶模式在截止频率附近相速度趋于无穷，搜索域上限如果只到cL，它们的根根本没有进入求解区间。

解决：把上界放宽到4倍cL以上，计算量增加不多，边界效果却很显著。低频段还要注意下界，A0在F很小时相速度低于0.1cT，下界用0.02cT在大多数情况下够用，如果扫F低于0.05，建议下界降到0.005cT。

6. 验证与延伸：把频散曲线落到激励设计上

6.1 三种低成本验证手段

逐点验证极限值是第一种手段。除了前面提到的板波速度和瑞利波速，还可以看A0在F趋近0时相速度和群速度同时趋近0，这是弯曲波低频行为。第二种手段是和文献里的标准谱图对比，重点看S1群速度曲线的“凹陷”位置，这个特征类似模式指纹，不同材料参数下凹陷位置变化明显，如果对得上，代码基本没有大问题。第三种手段是切换材料参数，比如换成铝，观察S0低频板速是否等于对应公式值，能检验代码对材料参数的依赖是否正确。

6.2 从频散曲线到激励信号设计

工程上拿到频散曲线后最常用的场景是选择激励频率。做法是看群速度曲线的平坦区，平坦区意味着传播时波形包络畸变小。比如要激发S0模式，就选S0群速度曲线在目标频段内梯度最小的点作为中心频率。另一个常见需求是估计两个模式到达时间差，以便设计时间窗分离波包。到达时间差公式为：

Δt = L·(1/cg1 - 1/cg2)

其中L是传播距离，cg1和cg2是两个模式的群速度。从曲线上读出对应F下的群速度值，直接代入就能算出分离窗口。这套MATLAB实现里附带的扫频和模式追踪脚本可以直接改参数跑通，不需要再从零搭框架。

6.3 一点个人习惯

以前我调一套复合材料的参数时，把搜索域上限设成cL，画出的曲线少了一整支高阶模式，当时盯着屏幕怀疑是求解器写错了，排查了两天才发现是边界问题。从那以后我每次换材料、换板厚，都会先跑一遍极限验证，确认S0低频板速和高频瑞利速都落在合理范围，再去看整条曲线。二维曲线好看是表象，求根稳定才是这套实现真正值钱的地方。希望帮到你。

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