MATLAB实现滑动窗口信息熵:声发射信号损伤演化分析
发布时间:2026/9/28 14:28:01来源:尧图网络
做声发射实验的同学十有八九都遇到过这种情况传感器采回来几十万甚至上百万个点波形看着密密麻麻但想从里面读出“信号变复杂了没有”“损伤发展到哪一步了”单凭肉眼根本看不出来。把整段信号直接算一个信息熵出来的只是一个孤零零的数字中途的变化全被平均掉了。这时候滑动窗口法就派上用场了——它能把信号切成一段段小窗口逐段计算信息熵得到一条随时间变化的熵值曲线配合MATLAB把曲线和信号图一起保存下来整套流程跑完整个损伤过程的“复杂度演化”就一目了然。这篇内容面向的是正在做声发射信号处理、结构健康监测或者材料损伤在线检测的朋友也适合刚接触信息熵概念的初学者。我会把滑动窗口法背后的设计逻辑讲清楚再给出完整的MATLAB代码、参数建议和踩坑记录照着改改就能用到自己的数据上。1. 声发射信号的信息熵到底反映了什么先说一个最基础的问题信息熵为什么能用来分析声发射信号。声发射Acoustic Emission简称AE是材料在变形、裂纹扩展、断裂等过程中瞬间释放弹性波的现象。传感器采集到的AE信号有几个明显特点非平稳、突发性、成分复杂。损伤越活跃信号包含的“信息量”往往越大频率成分越丰富波形形态越杂乱。香农信息熵就是衡量这种“杂乱程度”的指标。对离散信号把信号的幅值分布看成一组概率H -Σ p_i · log2(p_i)其中p_i是幅值落在第i个分箱中的概率。如果信号幅值分布非常集中比如全是零附近的小抖动熵值就低说明信号“没什么花样”如果幅值分布很均匀、各种大小幅值都出现熵值就高说明信号复杂、成分杂。单位是bit表示平均编码每个样本需要多少个二进制位。把信息熵用到AE信号上一个很直观的意义是它可以用一个标量刻画“声发射活动是否剧烈”。研究岩石破裂、金属疲劳、复合材料损伤时很多文献里都出现过类似的现象——临近破坏时熵值会出现明显抬升因为裂纹活动加剧、突发事件增多信号变得更加不可预测。不过这里有一个坑如果计算整段信号的全局熵不同时刻的损伤演化信息全被打包在一起结果既看不出“什么时候开始活跃”也看不出“破坏前有什么征兆”。尤其对持续数分钟的加载实验整段信号熵可能稳定在某个值附近但中间其实已经经历了多次明显的应力降。所以滑动窗口法的核心价值不是换一种算法而是换一个“时间尺度”。它把信号切成许多小切片按时间顺序计算每个切片的熵值把一维的时间序列“翻译”成一维的熵值演化曲线。这条曲线能直观展示信号复杂度什么时候升高、什么时候回落、什么时候持续攀升——这些对于判断损伤状态非常重要。2. 滑动窗口法的设计思路与原理2.1 为什么不用固定的一整段信号用一个生活化的例子类比你要分析一个人一整天的运动状态最原始也是信息量最低的方法是统计他今天“动没动过”——只要动过就算。更好的做法是每隔十分钟记录一次步数、心率这样就能看出他什么时候在跑步、什么时候在休息、什么时候心跳异常。整段信号计算全局熵就相当于第一招滑动窗口法相当于第二招。滑动窗口法的基本过程是设定一个窗口长度从信号起点开始截取第一个窗口计算该窗口的信息熵然后让窗口向前滑动一定步长继续截取、计算直到窗口移动到信号末尾。最终得到一条以“窗口中心时间”为横轴的熵值序列。这个过程有几个关键点需要提前想清楚第一窗口长度。窗口太短里面包含的AE事件太少幅值分布统计不稳定熵值会像噪声一样来回跳窗口太长时间分辨率太低短暂但关键的信号突变会被抹平。窗口长度的选择应该和信号本身的时间尺度匹配。比如采样率1 MHz、AE事件持续时间通常只有几百微秒到几毫秒窗口长度取1024到4096个采样点是比较常见的起点。第二步长。步长决定了相邻窗口的重叠程度。步长等于窗口长度相邻窗口不重叠计算量最小但曲线会呈锯齿状且可能漏掉跨窗口边界的特征步长取窗口长度的1/4相邻窗口有75%重叠曲线平滑很多代价是计算量增加。如果信号不算太长我一般取windowLen/4。第三直方图分箱数。信息熵计算依赖幅值分布概率而概率估计依赖怎么分箱。这是一个容易忽略但影响很大的细节。分箱数太少分辨不出幅值细节所有窗口的熵值都会偏低分箱数太多每个窗口里的样本数不够大量分箱概率为0统计波动剧烈。窗口长度是1024点的时候分箱数我取64到128之间比较稳。2.2 直接算熵值的两种常见做法对比这里多说一句计算AE信号熵值其实有两条路线一是直接对原始波形做幅值直方图统计属于“信号熵”二是先对信号做小波分解或经验模态分解再计算各尺度的熵值属于“尺度熵”。滑动窗口法这两条路线都能用但对大多数常规AE监测场景直接对原始波形做分箱统计已经足够实现简单、计算速度快、物理意义也清楚。如果后续想做更细致的区分可以在相同滑动窗口框架下把小波能量熵、样本熵都算一遍。不过这一步不是必需的先把普通的香农熵跑通形成套路再升级也不迟。3. MATLAB完整实现与代码解读3.1 数据读取与预处理我的代码以单通道AE信号为例假设数据已经导出为CSV文件或MAT文件。如果只有一列幅值数据、没有时间列需要根据采样率自己构造时间轴。如果是多通道数据建议先用索引把需要的通道取出来再做一次去直流偏置处理。AE传感器采集到的信号往往带有直流分量尤其是经过电荷放大器后容易出现基线偏移。直流分量对幅值分布影响很大会让直方图的主要概率集中在某个非零位置干扰熵值计算。所以读完数据后第一时间做去均值处理这一步很关键。% 1. 读取数据 % 方式A从CSV读取第一列为时间第二列为信号幅值单位可自行转换 data readmatrix(AE_data.csv); timeVec data(:,1); sig data(:,2); % 方式B如果直接从MAT文件读取 % load(AE_data.mat); % sig AE_signal; % 请按实际变量名修改 % fs 1000000; % 采样率 1 MHz % timeVec (0:length(sig)-1)/fs; % 构造时间轴 % 2. 去直流偏置 sig sig - mean(sig); % 3. 可选去除明显异常尖峰 % sigma std(sig); % sig(abs(sig) 6*sigma) 0; % 慎用仅在噪声毛刺明显时开启读数据时建议先做一次简单检查length(sig)、min(sig)、max(sig)看看幅值量级是否正常。有些采集软件会把幅值存储为int16读进MATLAB后变成长整型计算没问题但画图时纵轴刻度可能需要乘一个校准系数。幅值的绝对大小不影响熵值归一化概率后和幅值单位无关所以这里不用太纠结单位。3.2 信息熵计算核心函数下一步写一个独立的函数来计算某段信号的信息熵。这个函数只需要两个输入信号片段和分箱个数输出该段的熵值。写成函数的好处是主循环里调用方便之后换成小波能谱熵也只需要改这个函数内部。常见的信息熵实现方法是对信号画直方图统计落在每个区间内的样本数再归一化成概率。我在公式里加了一点保护逻辑——如果某段信号全为常数概率分布只有一个分箱非零熵为零如果某个分箱的概率为零直接跳过避免log2(0)产生NaN。function H shannonEntropy(seg, numBins) % 功能计算一段信号的信息熵 % 输入seg - 信号片段列向量numBins - 直方图分箱数 % 输出H - 香农信息熵单位为bit seg seg(:); seg seg(~isnan(seg)); % 去除NaN if isempty(seg) || numBins 2 H 0; return; end % 等间距分箱覆盖该窗口信号的幅值范围 edges linspace(min(seg), max(seg), numBins 1); counts histcounts(seg, edges); % 统计每个分箱的样本数 p counts / sum(counts); % 归一化为概率 p(p 0) []; % 去掉零概率项 H -sum(p .* log2(p)); % 香农熵 end这里有一个选型上的讲究为什么用linspace(min(seg), max(seg), numBins1)而不用MATLAB默认的histogram自动分箱自动分箱比如默认的auto会根据数据范围自适应调整分箱位置不同窗口的分箱边界还不一样。滑动窗口法要求各窗口之间“统计口径一致”如果每个窗口的分箱边界都变熵值之间就不可比了。固定numBins、固定等间距边界虽然不能保证每个窗口都分箱最优但保证了可比性这是滑动窗口法的核心原则。3.3 滑动窗口主循环接下来是主程序。窗口长度windowLen、步长stepSize、分箱数numBins这三个参数建议集中放在脚本开头方便后续调整。循环里不断取出新的信号片段调用熵函数同时记录窗口中心时间。主循环的时间复杂度是O(N/M×W)N为总长度M为步长W为窗口长度。如果信号很长且步长很小循环次数可能达到几十万次需要考虑后面第4节里的向量化方案。这里先用清晰明了的循环版本便于理解和调试。% 4. 滑动窗口参数设置 windowLen 2048; % 窗口长度根据AE事件持续时间调整 stepSize windowLen / 4; % 步长默认窗口长度四分之一 numBins 128; % 直方图分箱数 n length(sig); if windowLen n error(窗口长度超过信号总长度); end % 计算所有窗口的起始索引 startIdx 1:stepSize:(n - windowLen 1); numWindows length(startIdx); entropySeries zeros(1, numWindows); timeCenter zeros(1, numWindows); for k 1:numWindows idx startIdx(k):(startIdx(k) windowLen - 1); seg sig(idx); entropySeries(k) shannonEntropy(seg, numBins); timeCenter(k) mean(timeVec(idx)); % 取窗口中心的时间作为横坐标 end窗口中心时间的计算用mean(timeVec(idx))而不是timeVec(startIdx(k))。原因是当重叠比较大时窗口起点变化很小曲线会显得“堆积”在起点附近而用中心时间能更准确地反映熵值对应的时刻。如果你更在意“当前窗口结束时刻的熵”也可以取timeVec(idx(end))两种口径在绘图时不会改变曲线形状只影响横轴平移。3.4 图片绘制与保存从显示到输出信息熵曲线本身没多少可调整的空间图漂亮与否更多取决于绘图参数设置。我要重点说说图片保存因为这是很多人在MATLAB里最后一步栽跟头的地方——代码跑完图像文件要么空白、要么中文乱码、要么分辨率不够。先画曲线再设置字体、坐标范围最后用print输出。% 5. 绘图 fig figure(Color, w); set(fig, Position, [100 100 1100 420]); % 宽图适合展示演化曲线 plot(timeCenter, entropySeries, LineWidth, 1.5, Color, [0.15 0.35 0.65]); xlabel(时间 (s)); ylabel(信息熵 (bit)); title(声发射信号滑动窗口信息熵曲线); grid(on); set(gca, FontName, SimHei, FontSize, 11); set(gca, Box, on); % 6. 保存图片 outDir output; if ~exist(outDir, dir) mkdir(outDir); end saveas(fig, fullfile(outDir, AE_entropy_curve.fig)); % 保存FIG源文件 print(fig, fullfile(outDir, AE_entropy_curve.png), -dpng, -r300); % 高分辨率PNG使用print而不是saveas直接保存图片是因为saveas在保存PNG时默认分辨率往往不够字体也可能走样。print配合-r300可以指定300 DPI印刷和论文插图都够了。如果投期刊要求矢量图将保存格式改为-depsc或-dpdf出来的就是矢量格式放大也不会模糊。还有一个细节如果把fig创建时的Visible设为on窗口会在屏幕上闪一下。服务器或者批量处理场景下可以改成figure(Visible,off)计算完成后直接保存文件屏幕不显示速度也更快。3.5 参数怎么定一张速查表这三组参数的组合没有绝对标准但可以根据场景快速定初值。我做了个速查表按常见采样率1 MHz、AE事件持续1-2 ms来给参考值参数初值建议说明windowLen2048-8192点约2-8 ms越长越平滑越短越灵敏stepSize0.25×windowLen重叠率高曲线平滑numBins64-256窗口长时取大值窗口短时取小值采样率按实际记录时间轴换算用绘图分辨率300 DPI论文插图推荐如果是连续非突发的平稳信号比如匀速摩擦声窗口可以适当缩短比如1024点甚至512点熵曲线能体现出平稳背景下的微小波动。如果是离散事件很稀疏的信号比如单轴压缩下的大破裂建议窗口尽量包含1到2个完整AE事件否则很多窗口里没有任何事件熵值会长时间处于同一低水平看起来像一条“平线”信息量不够。4. 实际调试与常见问题排查4.1 熵值波动异常先查这三个地方我在调试中遇到的熵值问题大部分不是算法本身错了而是参数或数据细节没处理好。这里列三个最典型的排查方向第一个方向是分箱数与窗口长度的匹配。如果窗口长度只有512点却强行用256个分箱平均每个分箱只有2个样本统计噪声极大熵曲线毛刺会非常多。解决方法是让“每个分箱平均样本数”大于5即numBins windowLen/5。第二个方向是数据里的尖峰毛刺。AE信号如果采集时受到强电磁干扰会出现超出正常幅值数倍的尖峰。这些尖峰一旦进入窗口会大幅改变直方图的分布范围少数分箱占据极大权重导致熵值异常降低甚至变成“死水一潭”。处理方式是在预处理阶段用中位绝对偏差或幅值阈值识别并剔除尖峰。第三个方向是直流偏置没有去干净。如果窗口内存在缓慢变化的直流趋势幅值分布会出现多峰现象熵值虚高。对AE信号建议在滑窗前做一次高速滤波比如100 Hz高通而不是仅做减均值。毕竟减均值只能消除常数基线不能消除缓慢漂移。4.2 MATLAB中文注释与图片字体乱码这个问题在中文用户群体中特别常见值得单独说。MATLAB的编辑器默认编码在部分版本里是GBK注释里写了中文后到另一台机器或者换版本打开就乱码。图片里的中文乱码则是另一个问题系统没有对应字体或者字体名写错。对代码乱码有几种治本方案在MATLAB偏好设置里把文件编码改为UTF-8。不同版本路径不同一般位于“主页 → 预设 → 常规 → 编码”选择“UTF-8”。脚本开头不加任何中文注释图片标签里的中文改为英文“Entropy (bit)”。用外部编辑器按UTF-8编码写代码再在MATLAB里追加slCharacterEncoding(UTF-8)仅对Simulink模型有效脚本文件主要在预设中设置。对图片乱码注意MATLAB中文字体名不能随便写。常见可用的字体名有SimHei黑体、SimSun(宋体)、Microsoft YaHei微软雅黑。很多示例代码喜欢写FontName, 宋体在某些版本里也能显示但换到没有中文环境字体包的机器就会乱码。我一般统一用SimHei兼容性最好。如果最后要投英文期刊干脆全部用英文标注彻底避开字体问题。4.3 保存的图片内容异常全白、只显示一半、坐标轴丢失一个容易忽略的问题是如果不是用figure的句柄而是直接调用plotMATLAB会因为当前图窗被其他图形覆盖而把内容画到错误的对象上。稳妥做法是每个绘图段都显式使用fig figure(...)保存时也把fig作为第一个参数传给print。还有一个常见操作plot之后手动调整了坐标范围结果放大后发现曲线被截断。如果保存前调用了axis tight范围会缩到数据边界上但首尾的突变可能被裁掉。我建议先观察熵值范围再手动设置ylim给上下各留5%的余量曲线看起来更舒服。比如熵值范围是1.2到5.8设置成ylim([1 6])。如果保存的PNG在系统看图软件里显示“全白”检查print的格式字符串是否写对。print(fig, output.png, -dpng, -r300)中的-dpng不能写成-dpng -r300中间缺少空格也不能把输出文件名放在格式选项之后。文件名后跟格式选项是合法用法但多数人容易混淆建议保持“文件名紧跟print格式选项跟在文件后”这样直观的顺序。4.4 长信号下循环慢两种应对方案一个100万点的信号窗口2048点、步长512点大概会产生1952个窗口每个窗口调一次histcounts速度还能接受。但如果是10个通道、每个通道500万点循环可能要跑到十几万次MATLAB脚本循环慢的问题就暴露出来了。第一种方案是把循环换成arrayfun或纯粹的循环优化在进入循环前把信号切片缓存到矩阵里用buffer函数做滑窗矩阵化。buffer的用法是buf buffer(sig, windowLen, windowLen-stepSize)直接把整个信号切分成多列窗口最后一列不足时自动补零。切片完成后用矩阵运算一次性计算直方图速度提升显著。代价是内存占用上升大信号时会爆内存。第二种方案是先用MATLAB的parfor并行循环。如果电脑有多核把内层for改成parfor再设置parpool滑动窗口这种每个窗口独立计算的场景天然适合并行。但要注意如果电脑内存不大parfor会把数据复制到每个worker内存占用成倍增加反而可能更慢。我个人在超长信号上的做法是混合先取前1%的样本确定合理的windowLen然后用buffer分块、块内向量化计算最后再用parfor扫完全部数据。这个过程有一定编码量但跑起来很舒服。5. 实操记录一套岩石破裂实验数据的完整处理案例用一套我处理过的加载实验数据来说明整个流程。实验对象是花岗岩试样单轴压缩加载速率每秒0.5 kN同时用声发射仪记录波形。数据规模大约800万点采样率1 MHz时长约80秒。先做预处理。读取原始CSV后第一件事就是绘制整段信号的幅值时间图检查有没有明显异常段。当时的原始数据在58秒附近出现了几簇幅值突然减小又恢复的毛刺排查后发现是采集卡增益自动切换造成的不是真实的声发射事件。如果带着这些毛刺算熵58秒附近会有几个异常偏高的熵值尖峰误导判断。预处理完成后我设定了窗口长度4096点步长1024点分箱数128。这样每个窗口时长约4毫秒足够覆盖大多数岩石微破裂事件而且100 ms时间分辨率对压缩试验的损伤演化追踪也够用。计算出的熵值曲线整体趋势非常清楚加载初期0-20秒试样内部只有少量原生裂隙压密AE事件稀疏熵值平稳在1.8 bit附近波动。中期20-50秒微破裂事件逐渐增多熵值出现多次锯齿状上升每次上升对应一次局部能量释放。临近破坏50-78秒事件活动显著增强熵值整体抬升到3.5 bit以上且回落变少曲线形态从“锯齿”变成“台阶”。这组结果让我确定了滑动窗口熵值对破裂前兆的指示作用比单纯统计事件数要灵敏得多。在论文里我把原始的AE能量、振铃计数和滑动窗口信息熵放在同一张子图里对比熵值曲线的转折点明显领先于能量曲线的峰值。审稿人后来也专门确认了这一点因为这说明信息熵能反映损伤状态的“质变”而AE能量只能反映“量变”。6. 写给自己也给读者的几点体会写到这基本流程都说完了。最后分享几个我踩过几次坑之后形成的习惯。第一每次保存图片时顺手把对应的计算参数保存成一个文本文件。窗口长度、步长、分箱数、信号文件名、处理日期这些信息看起来不起眼但在你三个月后回看实验数据时如果没有记录就只能重新猜参数。我一般会在输出目录里自动生成一个params.txt用几条fopen加fprintf就能实现。第二对熵值曲线不要只看绝对数值。不同分箱数、不同窗口长度下熵值绝对大小没有可比性比较的是它的相对变化趋势。同一批数据里只要保证参数一致不同工况之间的曲线对比才是有效的。跨实验对比时建议先用统一的规范化流程相同滤波、相同分箱、相同窗口重新计算然后再比较。第三信息熵做损伤预警不是万能的。它在突发型AE事件明显的场景下效果突出但如果信号呈连续稳态背景比如水流噪声、机械振动噪声熵值曲线的辨识度会下降。遇到这类信号一个方向是做带通滤波把声发射频段以外的噪声去掉另一个方向是改用变分模态分解加小波熵逐尺度观察异常变化。如果你手上已经有采集好的AE数据不妨先按文中的代码跑一遍看看熵值曲线有没有帮助你理清损伤过程的时间段。这套方法改数据、改参数的成本都很低但得到的曲线会成为分析报告里非常有力的一张图。后续如果想扩展可以把滑动窗口信息熵与AE事件定位结合起来在损伤累计图上标注熵值跃升的时刻定位精度会更直观。
网站建设高端定制企业官网