做声发射检测或者材料损伤监测的朋友一定都遇到过这个尴尬的环节算法写好了传感器也接上了但手上没有真实缺陷信号连最基本的波动特性验证都无从下手。这时候用 MATLAB 模拟并绘制一个单个声发射波形就成了整个信号处理项目里的第一张多米诺骨牌。本文以 Planck 模型为例带你把从公式到波形的全过程跑一遍不仅给出可复现的代码还会讲清楚每一步为什么这么设置参数、哪些地方最容易翻车。适合刚接触声发射AE信号模拟的科研初学者也适合工程现场需要快速造一组测试样本的开发者。1. 为什么先说清楚“模拟单个声发射波形”这件事1.1 声发射波形的特征突发型瞬态信号先达成一个共识单个声发射事件比如一次微裂纹扩展、一次纤维断裂、一次摩擦滑动到达传感器时产生的波形不是持续的周期振动而是一个典型的瞬态脉冲。它有极短的上升前沿、幅度迅速爬到峰值随后又经过一个相对缓慢的衰减尾巴整个波形持续时间通常只有几百微秒到几毫秒。频域上则体现为宽带成分叠加但能量主要集中在传感器谐振频率附近的窄带范围内。这种“突然爆发又快速衰减”的形状决定了我们在模拟时不能拿一个简单正弦波敷衍了事。很多刚入门的朋友直接用sin(2*pi*f0*t)生成一段固定频率的正弦然后加个白噪声就当声发射波形用这其实是把“声发射”和“超声波连续波”混为一谈了。真实 AE 波形的识别特征恰恰在包络的不对称性和衰减趋势上。这也是我为什么选用 Planck 模型的核心原因——它能天然地生成一个“快上升、缓下降”的非对称包络和真实突发型 AE 信号在形态上非常贴近。1.2 模拟单个波形能做什么不能做什么模拟单个波形这件事定位上要清楚它不是用来替代真实实验的。它的价值主要体现在四个场景。第一算法开发阶段的基准信号。无论是做阈值触发、AR 模型特征提取还是基于机器学习的 AE 事件分类你都需要一个能调节“信噪比、幅值、主频、衰减速度”的干净样本用来验证算法逻辑是否跑得通。第二传感器与采集系统的功能测试。把模拟波形导出为任意波形发生器可以播放的文件在线缆、放大器和采集卡链路里过一遍可以快速确认系统增益和触发设置是否合理。第三教学演示与论文示意图。一个标注清晰的模拟波形远比一段噪声淹没的现场截屏更适合讲清楚 AE 特征参数的定义。第四为更复杂的传播仿真提供激励源。比如你在做有限元声传播模拟需要把“等效源信号”作为输入这时候单个波形模拟就是整个仿真的上游环节。但也要强调模拟波形不能用来验证诸如“动态裂纹扩展机制”“材料非线性行为”“多模态频散特性”这类深层次物理问题——那些需要精细的数值建模或真实实验标定。明白这个边界下面写代码时才不会过度追求复杂而失去实用性。2. Planck 波形模型原理理解与参数直观感受2.1 Planck 函数形式从黑体辐射到信号包络听到 Planck很多人第一反应是普朗克黑体辐射公式。没错我们确实是借用了那个著名的函数形式但物理含义完全不同。Primitive Planck 分布的形状是单峰非对称的——低频段快速爬升高频段拖尾下降。这种数学特征被信号处理领域借鉴过来用来构造声发射波形的包络$$h(t) \frac{A \cdot t^p}{e^{q t} - 1}$$其中p 控制上升沿的“锐利程度”q 控制衰减段的下降速度A 是幅度系数。t 为时间秒。当 t 很小时分母按泰勒展开可近似为 qt所以包络近似正比于 t 的 p-1 次方p 大于 1 时波形从零开始平缓爬升而不是瞬间跳变。这就天然避开了指数衰减包络在 t0 处“陡峭跳变”的伪迹更接近真实传感器响应的有限上升时间。需要说明的是直接把黑体辐射里的波长 λ 替换成时间 t 不是唯一做法也可以把 Planck 函数用在频域上构造频谱包络再通过逆傅里叶变换回到时域。但根据我的实操经验时域直接构造包络再与载波相乘参数更直观、调试更方便也是目前实践中最常用的方案。下面正文中的“Planck 模型”指的就是这种时域包络构造方法。2.2 参数与波形外观的对应关系为了让参数调起来有方向我给你列一个直观对照表这是我自己反复调试总结出来的经验值不一定绝对精确但方向可靠参数增大时的效果减小时的效果典型取值范围相对于 1 MHz 采样率p上升沿变缓峰值时刻后移上升沿变陡波形更像尖脉冲1.5 ~ 4q衰减加快波形总时长缩短衰减变慢拖尾变长2×10^4 ~ 8×10^4A整体幅值线性放大整体幅值缩小按需要归一化f0载波频率升高周期变密载波频率降低周期变疏100 kHz ~ 1 MHz从表里能看出p 和 q 分别主导了“波形的上升”和“波形的下降”两者相对独立。这就是 Planck 包络比单一指数衰减包络更有优势的地方你可以单独调节上升沿的形状而不影响衰减段的大致趋势。实际调试时我习惯先把 q 定下来让波形总长度符合目标再微调 p 来匹配上升时间。2.3 和其他包络模型的横向对比模拟 AE 波形常见的包络还有几种纯指数衰减包络、Gaussian 包络、Hann 窗包络以及将指数衰减和上升延迟结合的“延迟指数”包络。纯指数衰减包络实现最简单形式为 e^(-αt)但它在 t0 处幅值最大缺少上升沿和真实传感器冲击响应差得有点远。Gaussian 包络形状对称做正弦脉冲还行但真实 AE 波形普遍不对称“慢衰减”特征体现不出来。Hann 窗包络是限长周期信号的好选择适用于持续振铃型 AE但同样没有“物理衰减”语义。Planck 包络刚好处于“指数衰减”和“对称窗”之间的折中位置它既包含非对称的物理瞬态特征又不会因为上升沿缺失导致波形过于失真。而且它是光滑可导的做后续求导、包络分析时不会引入数值跳变。正因为这些特性我个人的经验是如果你模拟的是“突发型声发射事件”Planck 模型是性价比最高的选择如果你模拟的是“连续型声发射”如摩擦噪声那就应该换用随机噪声叠加窄带带通滤波的方案不要强行套 Planck。3. MATLAB 代码实现从公式到可落地波形3.1 环境与基础参数设定MATLAB 入门门槛低这里不需要任何额外工具箱——只用基础 MATLAB 即可甚至 Octave 也能跑。先定义采样参数和物理参数。采样率的选择是个关键点声发射传感器的主频通常在 100 kHz 到 1 MHz 之间根据奈奎斯特采样定理采样率至少是最高频率分量的两倍。但如果后续要做特征提取我建议采样率要高一些至少达到载波频率的 6~10 倍这样波形绘制出来也更平滑。% --- 基础参数 --- fs 1e6; % 采样率 1 MHz duration 1.5e-3; % 波形时长 1.5 ms t (0:round(fs*duration)-1) / fs; % 时间向量单位 s n length(t); % --- Planck 包络参数 --- p 2.5; % 上升沿指数系数 q 4e4; % 衰减系数 A 1; % 峰值幅度后面再归一化 % --- 载波参数 --- f0 150e3; % 主频 150 kHz phase 0; % 初始相位实际工程中如果你的传感器谐振频率是 150 kHz采集系统采样率一般会设置在 1 MHz 以上上面的参数就是一个很常见的组合。要注意的是round(fs*duration)这一步确保时间向量长度是整数否则后面索引、绘图都会出问题。3.2 生成 Planck 包络与载波合成接下来是核心代码。按 Planck 函数形式生成包络再和正弦载波相乘。这里有一个很容易踩的坑当 t 为 0 时分母exp(q*t)-1等于 0直接计算会得到 NaN导致整条波形崩溃。解决方法有两个一个是给分母加一个极小数eps另一个是从 t 略大于 0 的点开始计算。我习惯加eps这样能保持向量长度一致。% --- 构造 Planck 包络 --- denominator exp(q * t) - 1 eps; % 加 eps 防除零 envelope A * t.^p ./ denominator; % --- 构造载波 --- carrier sin(2 * pi * f0 * t phase); % --- 合成 AE 波形 --- ae_wave envelope .* carrier; % --- 归一化到 [-1, 1] --- ae_wave ae_wave / max(abs(ae_wave));注意几个细节。第一点运算不能少.^和./是对向量逐元素操作少了那个点就是矩阵运算直接报错或者结果完全不对。第二归一化放在合成之后、加噪之前这样后续加噪时幅值基准是明确的。第三合成出的ae_wave理论上应该像“一个带上铃响的脉冲”前面几个周期振幅迅速增大中间的峰值周期最明显之后振幅按指数规律缩小。如果画出来不是这个形状多半是 p 和 q 的量级没配对后面第 5 节我会专门讲调参规律。3.3 让波形更真实加噪、归一化、幅值微调干净波形有了但直接用干净波形去验证算法“过于乐观”。实际采集时前置放大器噪声、电磁干扰都会以带状噪声叠加到信号上。所以我通常会在模拟脚本里预留一个信噪比控制项。MATLAB 中可以使用awgn函数但这个函数属于通信工具箱很多环境不一定装。为了保险我习惯自己写一个加噪语句% --- 可选加入高斯白噪声设置信噪比 SNR --- snr_db 15; % 信噪比 15 dB noise_power var(ae_wave) / (10^(snr_db/10)); noise sqrt(noise_power) * randn(size(ae_wave)); ae_wave_noisy ae_wave noise;这里用var(ae_wave)算信号功率再根据目标 SNR 反推噪声功率。randn生成标准高斯白噪声乘以标准差后即得到指定功率的噪声序列。一个容易被忽略的点真实 AE 采集系统中的噪声背景往往不是白噪声而是带限噪声主要能量集中在传感器通带附近。如果项目要求高仿真度建议再加一个 Butterworth 带通滤波器对噪声整形。基础版本先用白噪声算法验证阶段的差异不大。最后幅值微调要看你的后续用途。如果模拟波形是给阈值触发算法用的可以把峰值幅值设置为传感器输出电压对应的伏特值比如 10 mV如果只是特征提取算法的输入直接归一化处理反而更稳定。4. 结果可视化与特征提取验证4.1 时域绘制与 AE 特征参数标注波形生成完画图是第一件要做的事。我推荐用两张子图第一张画完整时域波形第二张放大局部区间看细节。画图时顺手把 AE 信号最常用的几个特征参数标出来——峰值幅值、上升时间、持续时间、振铃计数。这些参数既是后续特征工程的基础也让读者一眼看出波形形态是否符合真实 AE 信号。figure; subplot(2,1,1); t_ms t * 1e3; % 转为毫秒 plot(t_ms, ae_wave_noisy, b); xlabel(时间 (ms)); ylabel(幅值 (V)); title(模拟声发射波形Planck 包络 正弦载波); xlim([0 duration*1e3]); grid on; subplot(2,1,2); zoom_range t 0.1e-3 t 0.4e-3; plot(t_ms(zoom_range), ae_wave_noisy(zoom_range), b); xlabel(时间 (ms)); ylabel(幅值 (V)); title(局部放大上升沿与主周期细节); grid on;上升时间在真实 AE 信号里普遍很短典型值在几微秒到几十微秒之间。如果你的 Planck 模型把上升沿拖到了几百微秒那看起来就不像突发型声发射更像一个慢变的冲击响应。后面第 5 节会给出如何按上升时间反推 p 参数的经验公式。4.2 频域验证FFT 与主频检查画完时域一定要看一眼频域。原因有两个一是检查主频是否真的落在 f0 附近二是检查 Planck 包络引入的频谱展宽是否处于合理范围。FFT 实现不过几行代码但坐标轴转成物理单位这一步经常有人做错。我给出经过验证的写法% --- 频域分析 --- N_fft 2^nextpow2(n); % 2 的幂次 FFT 点数 X fft(ae_wave_noisy, N_fft); X_mag abs(X(1:N_fft/21)); % 取正频率部分 f_axis (0:N_fft/2) * fs / N_fft / 1e3; % 转 kHz figure; plot(f_axis, X_mag); xlabel(频率 (kHz)); ylabel(幅值); title(模拟声发射波形的幅频特性); xlim([0 fs/2/1e3]); grid on;正常结果应该在 150 kHz 附近出现一个明显的谱峰同时在峰周围有一定带宽展宽——这是因为 Planck 包络的时域形状给载波调制出了边带频谱不再是单根谱线而是有一定宽度的峰。如果谱峰出现在 0 Hz 或者低频段异常偏高多半是波形没有去均值直流分量泄漏出来了。解决方法是先减去均值再 FFT即X fft(ae_wave_noisy - mean(ae_wave_noisy), N_fft);。另外要注意窗函数对频谱的泄露影响也不小。矩形窗的旁瓣比较高如果你发现谱图旁瓣和主瓣几乎一样高可以在 FFT 前加一个 Hann 窗。不过加了窗之后幅值会变小想读回准确幅值需要乘以 2 的修正系数。实测经验是做算法验证时先用矩形窗就是不加窗读主频位置再用 Hann 窗做精细谱分析两个结果对照更稳妥。4.3 包络提取匹配测试模拟波形做出来之后我还习惯做一道“自闭门测试”用常见的包络提取算法去还原包络验证提取结果是否和已知的 Planck 包络对得上。这不光是数学自检更是给后续特征提取算法打底。最简单的方法是 Hilbert 变换取包络MATLAB 里就是一行analytic hilbert(ae_wave_noisy); envelope_extracted abs(analytic);把envelope_extracted和原始的envelope画在同一张图上如果你没有加过量噪声两条曲线应该高度重合。注意 Hilbert 变换对“带限信号”效果最好而你的 AE 波形里载波频率 f0 一定远小于采样率 fs所以满足带限条件。如果没做归一化而且噪声设得很高比如 SNR 低于 5 dB包络提取会偏离原始包络较远这在真实场景中也一样说明你的模拟参数踩在了信噪比极限附近算法必须要用更稳健的包络估计方法。4.4 常见的不合理参数组合与调整经验把我在调试中遇过的“看起来奇怪”的波形总结一下给你当避雷指南。常用问题一整条波形几乎只有一个周期衰减快得像被一刀切。这通常是 q 取值过大导致的包络宽度远小于载波周期以致包络还没来得及爬升载波已经振完好几轮。解决办法是减小 q让包络持续宽度至少覆盖 5~10 个载波周期。经验公式持续时间粗略估计为 3 / q 到 6 / q 秒你可以反过来根据期望持续时间和载波频率选 q。常用问题二波形上升沿太陡峰值在首个周期就出现跟指数衰减包络没有本质区别。这通常意味着 p 取值过小。把 p 提到 2 以上上升段就会有一个明显的“渐入”过程。如果你希望上升时间在 20 微秒左右可以先设 p3再观察峰值位置是否在 20 微秒附近之后微调。常用问题三加噪声后波形看起来还是太干净不像实测信号。真实传感器波形往往伴随低频背景漂移。可以加一个低幅值、低频率比如 1 kHz 以下的慢变分量模拟基线漂移。这是一个很实用的工程技巧在算法验证时能暴露高通滤波不彻底的问题。5. 参数调整心得与常见坑5.1 Planck 参数的调节规律为了让你调参时不靠瞎猜我给出一个更定量的思路。假设你希望波形的峰值出现在 t_peak那么上节公式h(t) A t^p / (e^{qt} - 1)对 t 求导等于 0可以得到峰值时间点的近似关系。在 t 较小、e^{qt} - 1 ≈ qt的近似下退化出一个简化公式t_peak ≈ p / q。这个近似在 p 不太大时相当好用。举个例子如果我想让峰值在 60 微秒附近fs 为 1 MHz那么可以设 p3、q5×10^4t_peak ≈ 3 / 5e4 6e-5 秒 60 微秒刚好满足。然后在这个基础上微调 p 和 q 来适配上升沿细节和拖尾长度。这个方法帮我在项目里省了大量试探时间你可以先记住 p/q 的比值关系再去细化。再补充一个针对“振铃计数”的调整技巧振铃计数定义为包络超过阈值电平的次数。Pack 参数不变的情况下阈值电平设得越高超过阈值的次数就越少提取出来的特征差异很大。我在做特征一致性验证时会固定阈值电平为峰值幅值的 10% 或 20%并把这个设定写进脚本注释里防止后面忘了。5.2 采样率、时间窗、载波频率的匹配这三者的匹配关系是模拟波形工程化中最重要的部分。采样率决定时间分辨率时间窗决定能容纳多少载波周期载波频率决定波形振荡的快慢。三者的约束关系如下。载波频率 f0 不能太高。f0 必须小于 fs/2 是硬约束但实际操作中建议 f0 小于 fs/5。如果 f0150kHz、fs1MHz余量充足如果 f0800kHz、fs1MHz虽然满足奈奎斯特定理但模拟波形每个载波周期只有不到 2 个采样点绘制出来锯齿感明显特征提取精度也会下降。时间窗长度要和包络持续时间匹配。如果 duration 设得太短包络衰减段被硬截断频谱上会多出截断泄漏。我的习惯是把 duration 设为衰减到峰值 1% 所需时间的 1.2~1.5 倍。这可以用 q 粗估总时长约 6/q 到 10/q 秒。载波初始相位。在单波形模拟里phase0意味着波形从正弦零点开始。如果你希望波形从峰值开始让每次事件都能从最大幅值触发检测就把相位设成pi/2。两种设定都有意义但注意要和时间触发逻辑对应上。另外MATLAB 跑这类模拟几乎不耗时但对后续大规模参数扫描时间窗也不能无限拉长。数据点越多FFT 和 Hilbert 的计算量越大批量跑几千组参数时耗时差就很明显了。5.3 从单波形到批量参数扫描的脚本技巧单个波形验证完通常你还需要扫参数得到一批样本比如不同衰减程度q 变化、不同主频f0 变化的模拟 AE 波形用于后续分类网络的数据增强。这时候写循环要注意两点第一每次循环都要重新生成时间向量和随机噪声种子第二仿真的结果要打标签存档方便后续训练模型直接读取。我给一个精简的批量生成骨架% 参数列表 q_list [3e4, 5e4, 7e4]; f0_list [100e3, 150e3, 200e3]; n_files 10; % 每个组合生成 10 个样本 for idx_q 1:length(q_list) for idx_f 1:length(f0_list) for k 1:n_files % 更新参数 q q_list(idx_q); f0 f0_list(idx_f); p 2.5; duration 8 / q; % 动态设定时间窗 t (0:round(fs*duration)-1) / fs; denom exp(q*t) - 1 eps; envelope t.^p ./ denom; carrier sin(2*pi*f0*t); ae envelope .* carrier; ae ae / max(abs(ae)); % 加噪与保存 snr_db 10 5*randn; % 随机 SNR ae_noisy add_noise(ae, snr_db); save(sprintf(synthetic_ae_q%df0%dk%d.mat, idx_q, idx_f, k), ... ae_noisy, fs, f0, q, snr_db); end end end其中add_noise是前面第 3.3 节的加噪函数封装。批量生成的脚本价值在于可复现性所以每条数据的物理参数、信噪比都要随 mat 文件保存下来。训练分类器时这些标签就是你的 ground truth。5.4 波形导出与后续工程接口模拟做完最后一步通常是导出。两个常见场景一是写论文或报告需要高分辨率图二是接硬件测试需要波形文件。导出图片推荐exportgraphics这是 MATLAB R2020a 之后内置的高质量导出函数比print更清晰exportgraphics(gcf, ae_waveform.png, Resolution, 600);导出波形数据给任意波形发生器播放可用audiowrite或writematrix。如果后续要导入 Python 做深度学习建议直接存 CSVwritematrix([t(:), ae_wave_noisy(:)], ae_waveform.csv);CSV 两侧跨语言读取最省事不用纠结 mat 文件版本兼容性。需要注意的是真实声发射采集系统一般有自己的专有格式模拟波形导入时可能需要重采样或者幅值变换导出前先看清楚硬件接口要求不要导出完才发现采样率不匹配。6. 写在最后的实操体会这套 Planck 波形模拟方法已经在我的多个小项目里反复打磨过。做得次数越多我越觉得模拟信号这件事难点从来不在代码本身而在于你是否理解每个参数在物理世界中的对应物。p 对应能量释放的陡峭程度q 对应传播路径对信号的衰减吸收f0 对应传感器谐振特性——把这三个参数理解透你产出的就不只是一条数学上好看的曲线而是一条可以拿来和真实信号做对照、可以用来验证算法的工程样本。最后分享一个小技巧你可以在调参脚本里加一个和真实 AE 波形叠加的对比视图就是直接把现场采集的一段波形调进来和 Planck 模拟波形对齐画在一起。第一次看到二者形态接近时那种“模拟和现实碰上了”的可靠感是纯粹用代码堆出来的波形给不了的。从这个意义上说模拟单个声发射波形不是终点它是你后续做检测算法、做传播特性研究、做传感器布点验证的一块稳固的地基。