简介：这是一份基于MATLAB的Lamb波频散曲线计算与绘制资料，面向从事结构健康监测、无损检测及相关研究的工程师与科研人员。内容围绕薄板中Lamb波的传播特性，涵盖相速度与群速度的频散求解思路，并通过可视化图形直观展示波速随频率的变化规律，有助于理解频散机制并为缺陷检测、材料均匀性评估提供参考。压缩包共4个文件，包含1个m脚本、2个fig图形文件以及1个docx说明文档，整体约431KB。m文件承担频散计算与绘图核心逻辑，fig文件分别对应群速度和相速度曲线结果，docx则提供步骤说明与代码示例，便于对照学习。目前已有1014人学习下载，适合需要快速上手Lamb波频散分析、或希望借助现成代码开展仿真验证的初学者与工程技术人员。
1. 结构健康监测里，Lamb波频散曲线是一张比实验更早存在的地图
同样是一块 1 mm 厚的铝板，换能器频率从 0.5 MHz 提到 2 MHz，测到的 Lamb 波波形可能从“一个脉冲”变成“一串拖尾”，这不是仪器抖动，而是频散：每个频率分量以不同速度传播。工程上要区分 S0、A0、A1、S1，要选激励频率，要对长距离传播做补偿，全都要回到同一张图——Lamb波频散曲线。很多同行手里都有一份名为“Lamb波频散曲线.rar”的程序包，里面通常是某个平台的求解脚本或者现成数据表，但只会在界面里点“计算”是不够的。这篇文章把频散方程的来源、根扫描实现、实测信号对照和频散补偿一次讲透，适合超声无损检测、SHM 传感器布置和导波仿真方向的人读。
2. 从瑞利-兰姆方程到S/A模态：频散曲线不是查表查出来的
频散曲线在物理上由瑞利-兰姆（Rayleigh-Lamb）方程决定。平板上下表面自由，位移和应力场同时满足波动方程和零牵引力边界条件，最后得到的是一个关于频率和波数的超越方程，而不是某个显式解析式。这也是频散曲线麻烦的来源：方程不可直接求逆，只能在给定材料参数条件下数值找根。
2.1 S模态和A模态的分野来自板厚方向位移分布的对称性
取板厚为 d，半厚 h = d/2，纵波声速 cL，横波声速 cT，角频率 ω = 2πf，波数 k = ω/cp。按照平面应变假设，解可以按厚度方向位移分布分成两组：对称模态（S0、S1、S2……）对应面内位移对称、离面位移反对称；反对称模态（A0、A1、A2……）则反过来。代入上下表面应力为零的边界条件，行列式为零，得到下面两个实频散方程：
S（对称模态）: (q^2 - k^2)^2 * cos(p h) * sin(q h) + 4 k^2 p q * sin(p h) * cos(q h) = 0 A（反对称模态）: 4 k^2 p q * cos(p h) * sin(q h) + (q^2 - k^2)^2 * sin(p h) * cos(q h) = 0其中 p = sqrt((ω/cL)^2 - k^2)，q = sqrt((ω/cT)^2 - k^2)。注意这里 p、q 并没有规定必须为实数。低频段 S0 的相速度低于纵波声速，此时 p 本身是虚数；高频段模态趋近瑞利波速度，p、q 都可能变成虚数。很多初学者在这直接踩坑：直接在实数域开方，把 p^2 < 0 的区间全部丢掉，结果低频 S0 和高频 A0 全都画不出来。正确做法是用复数开方，最后统一判断方程函数本身的符号变化。
从方程结构还能看出一个特点：频率、波数、板厚总是以 ω、k、h 的组合出现，因此把横坐标写成频率与板厚的乘积 fd，相速度和群速度曲线对同一种材料就是唯一的。这就是为什么到处都强调用 MHz·mm 而不是只用 MHz 来标定 Lamb 波频散曲线。换板厚，曲线横坐标缩放即可，不需要重新测量材料声速。
2.2 相速度cp和群速度cg在频散曲线里各管一件事
相速度 cp = ω/k 描述的是某个频率分量的波前传播速度，群速度 cg = dω/dk 描述的是包络能量传播速度。无损检测里大家更关心后者，因为“信号到靶点的时间”和缺陷定位用的是群速度；传感器设计、波束角度计算则更依赖相速度。频散曲线的本质矛盾在于：cp 与 cg 并不相等，而且两者随频率变化的规律完全不同。低频 A0 的 cp 随频率升高而增大，cg 却在一个频率点附近达到峰值后再下降；S0 低频段两条速度都相对平稳，但到高频区和 A1、S1 一旦耦合，速度曲线会剧烈起伏。理解这张图，等于同时理解了“哪个频率段信号形状变化小”和“哪个频率段速度对厚度、温度最敏感”这两个工程问题。
3. 用Python解Lamb波频散曲线：从根扫描到可出图的最小脚本
求解瑞利-兰姆方程没有“一步到位”的数值库，常见做法是自己写根扫描。网络上流传的“Lamb波频散曲线.rar”里，多数实现也是这个思路，区别只在扫描轴选取、括号算法和分支归类。下面给出一个去掉界面包装、能在几分钟内跑通全过程的最小求解脚本。
3.1 在k域里求ω，而不是在f域里求cp
一个常见的错误做法是固定频率 f，扫描相速度 cp，找方程的符号变化。这样在截止频率附近会很难处理：A1、S1 等高阶模态在截止点处 cp 趋于无穷，频散曲线几乎是垂直的，步长稍微大一点就会漏根。我一般改成固定波数 k，在 ω 方向上找根。这样每个模态都是从截止频率出发、随 k 平滑延伸的一条曲线，根的位置更好追踪，分组也更自然。
另一个注意点是求根函数不能只在实数范围计算。p、q 为虚数时，方程函数可能变成纯虚数，用brentq前需要把“真正变号的量”取出来。下面的实现把复数结果按实部/虚部选择，保证 p、q 任意虚实组合下都有稳定符号。
3.2 完整求解代码与参数说明
import numpy as np from scipy.optimize import brentq import matplotlib.pyplot as plt # 材料参数：6061铝合金常用值，单位 m/s cL = 6.38e3 cT = 3.12e3 d = 1.0e-3 # 板厚 1 mm，换成其它厚度只需改这里 def dispersion_val(k, w, sym): """瑞利-兰姆方程的函数值。 sym=True 对应 S 模态，sym=False 对应 A 模态。 返回值是可用于符号判断的实数。 """ w = np.complex128(w) p = np.sqrt((w / cL) ** 2 - k * k) # 允许 p、q 为复数 q = np.sqrt((w / cT) ** 2 - k * k) ph = p * d / 2.0 qh = q * d / 2.0 k2 = k * k D = (q * q - k2) ** 2 B = 4.0 * k2 * p * q if sym: F = D * np.cos(ph) * np.sin(qh) + B * np.sin(ph) * np.cos(qh) else: F = B * np.cos(ph) * np.sin(qh) + D * np.sin(ph) * np.cos(qh) # p、q同实或同虚时F取实部或虚部，括号扫描时才有稳定符号变化 if np.abs(F.imag) < 1e-12: return F.real return F.imag f_max = 2.0e6 # 最高求解频率 2 MHz w_max = 2.0 * np.pi * f_max k_grid = np.linspace(1.0, 2.0e4, 400) # 波数网格，单位 rad/m w_grid = np.linspace(1.0, w_max, 500) # 频率网格，单位 rad/s roots_by_k = [] # 每个 k 下按频率升序存放各模态根 for k in k_grid: row = {'S': [], 'A': []} for tag, sym in (('S', True), ('A', False)): F = np.array([dispersion_val(k, float(w), sym) for w in w_grid]) for i in range(len(F) - 1): if not np.isfinite(F[i]) or not np.isfinite(F[i + 1]): continue if F[i] * F[i + 1] < 0: root = brentq( lambda ww, kk=k, ss=sym: dispersion_val(kk, ww, ss), w_grid[i], w_grid[i + 1] ) row[tag].append(root) roots_by_k.append(row) # 按分支序号画图：0--A0/S0，1--A1/S1，以此类推 fig, ax = plt.subplots(figsize=(8, 4.5)) for tag, color in (('S', '#1f77b4'), ('A', '#d62728')): for branch in range(6): cp_list, fd_list = [], [] for i, k in enumerate(k_grid): if branch < len(roots_by_k[i][tag]): w = roots_by_k[i][tag][branch] cp_list.append(w / k) fd_list.append((w / (2.0 * np.pi)) * d) if len(cp_list) > 2: cp_arr = np.array(cp_list) / 1e3 # m/s -> km/s fd_arr = np.array(fd_list) * 1e3 # m·Hz -> MHz·mm ax.plot(fd_arr, cp_arr, color=color, lw=0.9, label=tag if branch == 0 else None) ax.set_xlabel('fd (MHz·mm)') ax.set_ylabel('相速度 cp (km/s)') ax.set_xlim(0, 2.0) ax.set_ylim(0, 12) ax.legend() plt.tight_layout() plt.savefig('lamb_dispersion_cp.png', dpi=150)代码逻辑分三段：dispersion_val计算方程函数值；外层循环对每个波数 k 沿频率轴扫描，检测相邻两点符号变化后用brentq精确取根；最后按分支序号把同一模态连成曲线。扫描中每个 k 下得到的根按频率升序排列，顺序天然对应 S0、S1、S2……或 A0、A1、A2……但这种排序只在根未发生相交的频段有效，画高阶模态交叉区域时还要人工核对。
3.3 材料声速、板厚和网格步长怎么给
| 参数 | 推荐取值 | 说明 |
|---|---|---|
| cL | 铝取 6.38 km/s，钢取 5.90 km/s | 优先用同批次试块实测值 |
| cT | 铝取 3.12 km/s，钢取 3.20 km/s | 可用横波直探头测，也可由剪切模量和密度换算 |
| d | 按真实板厚 | 只影响 fd 横轴，不改变模态个数 |
| k_grid 点数 | 300~800 | 太疏会漏掉高阶模态靠近截止处的细碎根 |
| w_grid 点数 | 400~600 | 太疏会在模态曲线接近垂直时漏括号 |
w_grid 的上限 f_max 决定了能画到多少阶模态。比如 2 mm 铝板画到 4 MHz·mm，网格取 500 点能稳定出到 A2、S2；再往高频率走，模态间距变小，建议把 w_grid 提到 800 点以上，否则高阶模态之间会出现“断线”。这里给的材料声速只是典型值，实际构件涂层、各向异性、温度变化都会让频散曲线偏移，严谨的工程计算应以实测声速为准。
4. 实测Lamb波信号与频散曲线对照：模态识别与激励频率选择
算出频散曲线只是第一步。做检测时真正的问题是：手头激励出的这一串波形里哪一段是 A0，哪一段是 S0？频散曲线必须能和实测数据对上，否则曲线就只是一张好看的图。
4.1 从曲线反推激励频率和模态组合
频散曲线可以直接指导激励频率选择。比如薄板缺陷检测常用 A0，因为其离面位移大、对表面和近表面缺陷敏感，但低频 A0 的群速度随频率变化剧烈，传播 20 cm 后脉冲会展宽；S0 低频段频散小，但波长短、对表面质量要求高。更合理的做法是先根据频散曲线画出“fd-群速度”图，选群速度曲线比较平缓的频率段作为激励中心频率。一般我会这样选：A0 避开 cg 峰值附近，选在峰值右侧下降较缓的频段；S0 选 fd ≤ 1 MHz·mm 的低频区，此时能量集中且波形接近无频散。
| 检测目的 | 推荐模态 | 推荐 fd 范围 | 主要理由 |
|---|---|---|---|
| 表面缺陷、涂层剥离 | A0 | 1~2 MHz·mm | 离面位移大，导向性好 |
| 长距离焊缝扫查 | S0 | 0.5~1 MHz·mm | 频散弱，回波容易解释 |
| 厚度减薄检测 | S1/A1 附近 | 靠近截止频率 | 高阶模态对厚度变化敏感 |
| 缺陷成像算法验证 | A0 单一模态 | 先选窄带激励 | 减少多模态叠加干扰 |
表中“单一模态”四个字很关键。Lamb 波实验做不好，多数不是因为曲线算错，而是激励信号带宽太宽，同时在板内激发出 A0、S0，甚至 A1。激励信号一般用汉宁窗调制的 5~10 周期正弦脉冲，带宽由周期数控制，周期数越多带宽越窄，模态越纯。
4.2 用二维FFT验证实测信号里的模态成分
要确认实测波形里到底有哪些模态，最有效的办法是做二维傅里叶变换。沿传播方向等间距布置一列测点，记录每个位置的时间信号，得到一个“空间-时间”矩阵 S。对 S 做二维 FFT 后，幅值峰值所在位置就是实测频散关系。
# S: (M, N)，M为沿传播方向的测点序号，N为时间采样点 S = np.array(scan_signals) S_hat = np.fft.fftshift(np.fft.fft2(S)) dx = 1.5e-3 # 测点间距，单位 m dt = 2.5e-8 # 采样时间间隔，单位 s M, N = S.shape k_axis = np.fft.fftshift(np.fft.fftfreq(M, dx)) * 2.0 * np.pi # 波数 rad/m f_axis = np.fft.fftshift(np.fft.fftfreq(N, dt)) # 频率 Hz extent = [f_axis.min() / 1e6, f_axis.max() / 1e6, k_axis.min(), k_axis.max()] plt.imshow(np.abs(S_hat), aspect='auto', extent=extent, cmap='inferno') plt.xlabel('频率 (MHz)') plt.ylabel('波数 k (rad/m)')二维 FFT 里每一条亮线对应一种 Lamb 波模态。理论频散曲线可以按 k = 2πf / cp 换算到同一张图上叠加，峰值落在理论线上，就证明该模态确实被激励出来了。实际操作中测点间距 dx 必须小于最短波长的二分之一，也就是成像采样定理；波长可以从频散曲线上读，取所用频带内最小波长，dx 至少比它小一半。
二维 FFT 还有副作用：它能直接区分传播方向上的正负行波，反射波和端面回波会显示在负波数一侧。做模态识别时只看正波数半平面，能避免把反射波误判成第二模态。
5. Lamb波频散补偿与“按图反查”的两个实用技巧
曲线算出来、模态也对上了，下一步要解决的是“频散让波形变宽”的问题。如果已经知道传播距离 L，并且从频散曲线上读到了波数与角频率的关系 k(ω)，就可以在频域做相位补偿，把所有频率成分的相位拉回到同一起点。
def compensate_dispersion(signal, fre, k_fun, L): """按照频散关系 k(ω) 对传播距离 L 做相位补偿。 k_fun: 由频散曲线插值得到的函数，输入频率返回波数 """ spec = np.fft.fft(signal) spec = spec * np.exp(-1j * k_fun(fre) * L) return np.fft.ifft(spec).real # 频率轴需要和fft结果对齐 N = len(signal) dt = 2.5e-8 fre = np.fft.fftfreq(N, dt) # k_fun 可以由前面求根结果插值得到，例如针对A0模态 # from_scipy_interpolate: k_fun = interp1d(f, k, bounds_error=False, fill_value=0) corr = compensate_dispersion(raw_a0, fre, k_fun, 0.3)补偿后的波形中，原本被拉伸包络会重新变窄，这有利于精确读取飞行时间。要注意补偿函数只对单一模态有效，如果信号里同时混着 A0 和 S0，要先通过窄带激励或二维 FFT 把目标模态分离出来，再各自做补偿。
最后一个实用技巧是注意频散曲线数据表的分支对齐。从网上下载的“Lamb波频散曲线.rar”解压后，常见格式是四到六列文本，分别存 fd、cp、cg 和模态标识；有些脚本把 S0、A0 分文件存，有些则混在一个数据矩阵里。使用前先画散点图看分支是否连续，凡是出现跳变的行，基本是根扫描时把两个模态接在了一起。对于工程应用，只保留自己关心频段内的连续分支，再对 k(ω) 做插值，否则频散补偿会把错误的相位修正叠加到信号上。
频散补偿做完后，可以进一步用“误差包络法”验证：把补偿后的 A0 波包中心与激励脉冲对齐，对比两端的包络宽度。如果补偿正确，波包宽度应接近原始激励脉冲宽度，偏差通常能压到半个周期以内，这也是频散曲线从理论表格变成现场判据的最终一步。
本文还有配套的精品资源，点击获取
