★ FEATURED ARTICLE 第一次接触声发射波形是在给金属试件做拉伸实验的时候。传感器贴在试件表面材料一有裂纹扩展主机上就蹦出一个像是被狠狠衰减过的振荡脉冲——信号持续时间很短、幅值一下拉起来又指数往下掉。当时老师傅说这叫声发射事件是一个瞬态能量释放过程。后来自己用MATLAB模拟单次声发射事件才发现这东西用几分钟就能写出来但要让波形“像回事”还得给频谱找个物理依据这也是这篇以Planck谱做加权包裹的主线思路。这篇文章适合三类人刚接触声发射检测、想快速生成仿真波形做算法验证的工程师做无损检测课程设计或论文复现的学生以及那些手头有真实AE数据、想理解“为什么波形长这样”的信号处理爱好者。我会从最经典的衰减正弦模型讲起再引入Planck谱做多频叠加最后给出可直接跑的MATLAB代码和一堆实际踩过的坑。1. 把“声发射”和“Planck”放在一起到底能模拟什么1.1 声发射波形为什么是衰减振荡声发射简称AE指材料内部因局部应力集中、裂纹扩展、位错运动或摩擦等原因突然释放应变能并产生弹性波的现象。这个波传到材料表面被压电传感器拾取就成了我们看到的AE波形。单次声发射事件在时域上有两个非常明显的特征一是持续时间极短通常只有几十微秒到几毫秒二是幅值急剧上升后按指数规律衰减。为什么是这种形态因为能量释放是瞬时的波在结构内部来回反射每一次反射都会损失一部分能量传感器本身也有阻尼特性所以接收到的信号自然就是“一个快速起步、拖着一条尾巴”的衰减振荡。打个比方你拨一根琴弦声音不是突然消失的而是越来越弱直到听不见。AE波形就是那根弦只不过它的频率更高、衰减更快而且往往不只一个频率在响。1.2 Planck在这里扮演的角色Planck也就是普朗克常数、普朗克辐射定律那个Planck。很多人第一反应是黑体辐射跟声发射有什么关系关系不在于“热辐射”而在于谱加权思想。普朗克定律描述了热平衡状态下不同频率的能量分布密度。如果把材料内部一次声发射事件看成大量晶格振动的集体释放类比成一个个“声子”的发射过程那么不同频率分量携带的能量强弱就可以借鉴Planck分布来做加权。说得直白点我不再让所有频率的音量一样大而是按照某个谱形决定哪个频率贡献多、哪个贡献少。这里必须说明一个工程上的澄清声发射频段的频率只有几十到几百千赫对应的声子能量h·f非常小远小于k_B·T所以普朗克公式在这个频段实际上退化成瑞利-琼斯形式谱形近似按频率平方上升。这意味着如果你严格用真实普朗克常数去算得到的是一个高频更强的包络而不是经典黑体曲线那个“先升后降”的钟形峰。这并不妨碍我们把它当做一个建模工具反而给了我们一个非常灵活的频域加权函数。本文就是基于这个思路来做的。2. 建模思路与数学表达2.1 经典单事件模型绝大多数声发射波形模拟用的都是下面这个经典表达式x(t) A·exp(-α·(t - t0))·sin(2π·f0·(t - t0))·u(t - t0)其中A信号幅值可以理解为声发射事件的能量强度α衰减系数单位1/sα越大波形衰减越快f0中心频率单位Hz对应传感器谐振频率或AE源主频t0事件到达时刻也就是波形在时间轴上“蹦出来”的那一瞬u(t-t0)单位阶跃函数保证 t 小于 t0 时波形为0避免出现“负时间”信号。这个模型的好处是简单、可解释性强。你调大α波形尾巴变短调高f0波形变密改A改变能量。很多商用AE仿真软件里的标准数据库本质上也是在用这个模型配合不同参数生成的。但它的局限也很明显只有一个频率。真实声发射信号往往是多模态的材料边界反射、传感器响应、不同传播路径叠加在一起波形频谱是宽带的。想更贴近实际就得把多个频率分量合起来。2.2 Planck谱加权怎么做既然要多个频率就得给每个频率分配一个权重。这里我用Planck分布作为权重函数写成频率域的能量密度形式B(f) (8π·h·f³ / v³) / (exp(h·f / (k_B·Tem)) - 1)参数含义如下h普朗克常数6.62607015e-34 J·sk_B玻尔兹曼常数1.380649e-23 J/KTem等效温度注意这里不是材料实际温度而是描述AE源能量释放强度的等效参数v介质声速钢中纵波大约5900 m/sf频率。实际写代码时我一般先定义一段频率网格比如50 kHz到500 kHz每10 kHz取一个点然后按这个公式算出每个频率的权重再归一化。归一化后B(f)的形状就决定了哪些频率成分在合成波形里更突出。2.3 从单频到多频的合成策略合成思路不复杂把经典模型里的单个正弦项替换成多个正弦项的加权叠加。具体做法是对频率网格里的每一个频率fi生成一个衰减正弦信号幅度乘上对应的Planck权重B(fi)然后全部累加。写成公式就是x(t) Σ_i B(fi)·exp(-α·(t - t0))·sin(2π·fi·(t - t0))·u(t - t0)为什么这样做因为真实AE事件在频域里就是一个连续谱不同频段的能量在衰减过程中并不是独立的而是同时在结构里传播。叠加出来的波形会更接近传感器实际拾取到的信号看起来更“毛糙”、更真实而不是一个干净到不自然的单频正弦。我在实际项目中经常用这种多频叠加的信号来当合成AE数据用于测试阈值触发算法、到达时差定位算法和机器学习分类模型。因为真实数据量不够时这种物理上有依据的仿真数据能极大地扩充样本多样性。3. MATLAB实操从参数设置到出图3.1 环境准备与基础参数设置MATLAB版本没什么特殊要求R2019b以后都能跑不需要额外工具箱只用最基础的语言和绘图函数。建议先清一下工作区避免历史变量干扰。clear; clc; close all; % 采样参数 fs 5e6; % 采样率 5 MHz T_dur 1e-3; % 信号时长 1 ms N round(T_dur * fs); % 采样点数 t (0:N-1) / fs; % 时间轴 % 声发射源参数 f0 150e3; % 中心频率 150 kHz alpha 2.0e4; % 衰减系数 2e4 1/s A0 1; % 幅值 t0 2e-4; % 事件到达时刻 200 us这里有几个关键点。采样率5 MHz对150 kHz中心频率来说完全够用满足了奈奎斯特条件也为FFT观察频谱留足了带宽。信号时长1 ms在5 MHz采样率下就是5000个点绘图、运算都很舒服。3.2 方法A经典衰减正弦生成% 生成经典单频AE波形 idx t t0; x_classic zeros(size(t)); x_classic(idx) A0 * exp(-alpha * (t(idx) - t0)) .* sin(2 * pi * f0 * (t(idx) - t0));这一小段代码是整个仿真的核心逻辑。先用t t0制作掩码idx确保事件到达前的样本都是0到达后用指数项exp(-alpha·dt)控制衰减用sin(2π·f0·dt)产生振荡。我一开始写的时候偷懒直接在整个时间轴上算exp(-alpha·(t-t0))结果t小于t0的部分指数项变成正的叠加出莫名其妙的振荡。后来才学乖了一律先做逻辑掩码再对子集运算。这一点在后面的Planck叠加里同样重要。3.3 方法BPlanck谱加权多频叠加% 声速与物理常数 v 5900; % 钢中纵波声速 5900 m/s h 6.62607015e-34; % Planck 常数 kB 1.380649e-23; % Boltzmann 常数 Tem 850; % 等效温度 % 频率网格与Planck权重 fc (50:10:500) * 1e3; % 50 kHz ~ 500 kHz Bw (8 * pi * h * fc.^3 / v^3) ./ (exp(h * fc / (kB * Tem)) - 1); Bw Bw / max(Bw); % 归一化权重 % 生成Planck加权的多频AE波形 x_planck zeros(size(t)); for k 1:length(fc) tmp zeros(size(t)); tmp(idx) Bw(k) * exp(-alpha * (t(idx) - t0)) .* sin(2 * pi * fc(k) * (t(idx) - t0)); x_planck x_planck tmp; end x_planck x_planck / max(abs(x_planck));这段代码最值得注意的是分母里的exp(h·f / (kB·Tem)) - 1。因为h·f比kB·Tem小好几个数量级这个指数项非常接近1分母很小整体计算不会溢出但确实能看出低频和高频权重的差异。你可以试着把Tem分别改成300、850、5000再画Bw曲线对比会发现整体形状几乎不变只是幅值成比例变化。原因就是前面提到的经典极限。想看到Planck分布明显的“尖峰”等效温度要拉到10^6 K量级但那已经不是AE频段的真实物理了工程上不建议这么调。3.4 绘图与频谱对比% 时域波形对比 figure(Color, w, Position, [100, 100, 1000, 750]); subplot(2, 2, 1); plot(t * 1e3, x_classic, b, LineWidth, 1.2); grid on; xlim([0, 1]); xlabel(时间 (ms)); ylabel(幅值); title(经典单频衰减AE波形); subplot(2, 2, 2); plot(t * 1e3, x_planck, r, LineWidth, 1.2); grid on; xlim([0, 1]); xlabel(时间 (ms)); ylabel(幅值); title(Planck加权多频AE波形); % 频谱对比 X1 abs(fft(x_classic)); X2 abs(fft(x_planck)); f_ax (0:N-1) / N * fs; subplot(2, 2, 3); plot(f_ax(1:N/2) / 1e3, X1(1:N/2), b, LineWidth, 1.2); hold on; plot(f_ax(1:N/2) / 1e3, X2(1:N/2), r, LineWidth, 1.2); grid on; xlim([0, 600]); xlabel(频率 (kHz)); ylabel(幅度谱); legend({经典单频, Planck多频}, Location, northeast); title(FFT频谱对比); % Planck权重曲线 subplot(2, 2, 4); plot(fc / 1e3, Bw, k, LineWidth, 1.5); grid on; xlabel(频率 (kHz)); ylabel(归一化权重); title(Planck谱权重曲线);跑完这段代码你会看到两个关键结果。第一经典单频波形频谱只在一根谱线上有明显能量而Planck多频波形频谱铺开了一个宽带区域看起来更接近真实AE信号。第二Planck权重曲线在50~500 kHz范围内整体是随频率上升的所以合成波形里高频成分占比更大波形振荡得更细碎。如果觉得高频太强可以把频率网格上限降低或者用更低的等效温度搭配一个自定义传感器传递函数来整形。4. 参数敏感性、物理意义与工程建议4.1 四个关键参数怎么调参数取值范围建议对波形的影响物理意义fs采样率2.5~10 MHz决定波形时间分辨率太小会让高频成分混叠数据采集卡采样率f0中心频率50~500 kHz决定振荡疏密程度越高峰值频率越高AE源主频或传感器谐振频率alpha衰减系数1e3~1e5衰减越快波形尾巴越短能量越集中结构阻尼和传播路径损耗t0到达时刻通常几十到几百微秒决定窗口内事件出现的位置AE事件发生并传到传感器的时间Tem等效温度300~5000 K调幅值经典极限下对谱形状影响很小AE源能量释放强度的等效描述实际调试时我通常先用alpha2e4跑出一版波形看尾巴长度合不合适再调采样率和Grid。alpha太小波形拖得很长触发算法会把多个事件混在一起alpha太大则衰减太快幅值很快就变成噪声水平定位算法反倒找不到峰值。4.2 我在调参时踩过的坑第一个坑是给t0留的静默段太短。如果t0设在10 us以内你在时域图里几乎看不到信号起始的“台阶”整条波形看起来像是从0秒就开始振荡不利于展示事件到达特征。我一般把t0放在信号总时长的20%左右比如1 ms的信号放在200 us这样既有清晰的触发沿也保留了前置静默段后面用来测试阈值检测很顺手。第二个坑是频率网格步长选太大。如果你只用50 kHz和500 kHz两个频率点叠加波形会周期性地“打架”出现明显的拍频现象。网格步长最好在5~10 kHz以内叠加出来才连贯。计算量完全不用担心5000个采样点乘几十个频率循环一次也就毫秒级。第三个坑是FFT画谱时只看前半段却忘了频率分辨率。频率分辨率是1/T_dur1 ms信号对应1 kHz分辨率这对150 kHz量级的主频来说完全够用但如果你要分辨靠得很近的两个峰值就得延长信号时长。补零只能让谱线更平滑并不能真正提高分辨率这个细节很多人会忽略。4.3 这个模型怎么用到真实工程里真实AE采集系统里信号经过传感器、前置放大器、带通滤波器之后波形已经不是源信号的原始样貌。我建议在仿真波形后面串联一个简单的传感器传递函数比如用一个谐振频率150 kHz、品质因数Q10的二阶带通滤波器去卷积合成信号波形会立刻“失真”出压电传感器的味道。另一个工程实践是拿这个仿真波形测试声发射阈值触发。给x_planck加上-20 dB左右的随机噪声设置一个固定阈值统计触发时间和真实t0之间的偏差就能评估不同信噪比下事件检测的可靠性。这比直接在真实数据上反复实验要快得多也能提前发现参数设置的问题。5. 常见问题速查与排查技巧5.1 画出来的波形为什么没有衰减尾巴如果你看到波形从头到尾幅度差不多先检查alpha是不是太小比如设成了2e2而不是2e4。再检查是不是把exp(-alpha·(t-t0))写成了exp(alpha·(t-t0))符号反了尾巴不但不衰减还会越来越大。最后检查时间间隔dt到底是秒还是毫秒t的单位不一致会导致alpha看起来失效。5.2 FFT频谱为什么出现奇怪的周期成分最常见的来源是信号末尾被硬截断。时域波形在1 ms结束时如果幅值还没衰减到接近零FFT就会把截断处当成一个阶跃频谱里出现旁瓣振荡。处理办法有两个一是把T_dur加长到3~5 ms让衰减更充分二是给时域信号加一个Hann窗再去做FFT。注意加窗会改变幅值做归一化时要用窗函数能量修正。5.3 hilbert包络不光滑怎么办用hilbert提取包络时如果信号里混入高频噪声包络会毛刺明显。不要直接修改原始波形改用滑动RMS窗口来画包络窗口长度取10~50个振荡周期。比如150 kHz信号周期约6.7 usRMS窗口取100 us左右画出来的能量包络就非常平滑也方便和朋友解释“AE事件能量随时间衰减”这个物理量。5.4 多频叠加后波形幅值忽大忽小这是相位干涉的正常现象不是代码错误。不同频率分量在某个时间点同相叠加幅度就高反相抵消幅度就低。想要稳定幅值可以像代码里那样最后做一次max归一化或者在每个频率分量上随机一个初始相位多次实验取平均包络。我在做批量生成训练数据时会给每个分量加随机相位获得的样本多样性会好很多。6. 最后分享一个我常用的经验和习惯做单次AE波形模拟时我很少只出“一张完美图交差”。我会把参数定义成一个结构体或者参数表alpha、f0、Tem、频率网格全部集中管理然后一口气生成几十组不同参数的波形批量保存成MAT文件。后续做分类模型训练、阈值算法验证、定位精度测试想用哪组调哪组不用反复改代码重跑。另外一个小习惯每次改完参数先画Planck权重曲线看一眼再画时域波形。权重曲线能直接告诉你哪些频段在主导信号如果发现合成波形和预期不符十有八九是权重形状不对而不是时域生成代码有问题。这种“先调谱再调时域”的顺序能帮你少走很多弯路。这个仿真模型后续还可以继续扩展比如加入多个声发射事件、模拟传感器阵列接收到的多通道信号、给每个通道加不同的传播时延和衰减系数就能从“单波形模拟”升级成“整场声发射事件定位仿真”。到那一步你手里这份Planck加权波形的价值会被放大很多倍。
阅读完成 · 觉得有帮助？
