最近在复现一篇FRP复合材料加载试验的声发射数据分析图卡在一个很基础但特别容易绕晕的词上——“强度分析”。论文里用到了两个缩写HI和Sr图表看着很高级一条曲线平缓另一条在某处突然抬头。但等我翻到方法部分只有短短一段公式和一句“按ASTM标准计算”细节全靠猜。折腾了两天试了好几种写法总算是把曲线和原文对上号了。这篇博客就是把我摸索出来的完整过程写下来从指标的物理意义、Matlab实现、到绘图复现每一步都给出可以直接跑的代码和避坑提醒。无论是做复合材料损伤检测、金属裂纹监测还是纯想搞懂声发射强度分析的新手都能按这份记录复现出自己的HI和Sr图。1. 先搞清楚HI和Sr到底在算什么很多教材把声发射信号处理讲得玄乎其实落到实际上强度分析的核心思路非常直白声发射事件不是均匀发生的强度突变往往对应损伤模式的改变。1.1 信号强度为什么比计数可靠声发射系统通常会记录好几个参数到达时间、峰值幅度单位通常是dB AE、能量、持续时间、上升时间。论文里如果出现“signal strength”指的往往是信号强度这个概念单位一般是伏·秒或皮欧·米也可以理解成信号包络的积分面积。一个常见误区是拿“事件计数”当主要指标。计数的问题是一次大裂纹扩展产生的信号和一百次微小基体开裂产生的信号计数值可能差别很大但对材料损伤的实际贡献完全不同。信号强度保留了信号的“大小”信息所以强度分析更能反映真实损伤。简单说信号强度序列就是按时间顺序排列的S_ii1,2,...,N后面的所有计算都围绕这个序列展开。1.2 HI历史指数找累积曲线上的“抬升点”HI这个名字有不同的展开写法常见的是Histogram Index或Historical Index中文翻译成历史指数。它干的事情是考察截止到当前时刻累积信号强度曲线的“局部斜率”是否比历史平均更高。公式在ASTM标准里的定义大概是[ H(t) \frac{N}{N - k} \times \frac{\sum_{ik1}^{N} S_i}{\sum_{i1}^{N} S_i} ]其中N是截止当前时刻的事件总数S_i是第i个事件按到达时间排列的信号强度值k是偏移量取值跟N有关。这个公式的意思是把前面k个事件剔除后剩下这部分信号强度占整体累积的比例再乘以一个缩放系数。如果损伤进入剧烈阶段尾部信号强度占比突然变大HI就会显著上升。k的取值每个标准版本略有出入我用的惯例如下当前事件数 N偏移 kN ≤ 50051 ≤ N ≤ 20025201 ≤ N ≤ 50075N 500125需要说明的是ASTM标准本身和后来各路论文里的实现不完全一样做复现时第一件事是看你参考的论文怎么定义的k。后面代码里我会把k做成一个可配置参数方便切换。1.3 Sr严重率盯着最大信号强度的尾巴SrSeverity严重率的定义比HI简单它衡量的是当前时刻之前最大的J个信号强度的平均值[ S_r \frac{1}{J} \sum_{j1}^{J} S_{max,j} ]J在不同文献里有取10的也有取事件总数一定比例的。常见的简化实现是取J min(10, N)也就是最多10个最大信号强度求平均。如果当前事件数太少就取全部。从物理上讲HI反映的是累积历史中“增量异常”而Sr反映的是“最强事件有多强”。二者配合使用既能找到剧烈活动的起点也能定性判断活动猛烈程度。很多标准里会用Sr作为损伤严重度的直接判据单看数值大小决定是否达到临界状态。2. 数据准备光有采集文件还不够2.1 从采集系统导出什么样的数据声发射采集系统比如Physical Acoustics的PAC系统、Vallen等一般都能导出ASCII文本或CSV文件。列名可能五花八门但关键的几列通常包括时间或到达序号Time/Segment、峰值幅度Amplitude、信号强度Signal Strength。有个很容易踩的坑很多系统导出的是带行号的自由格式文本表头在第二行列之间用制表符或逗号混用。直接readmatrix可能报错建议先看一眼文件头几行再决定导入方式。我建议的第一步不是急着写计算函数而是把数据读进来确认列对得上。以最典型的三列CSV为例% 读取声发射导出数据 filename AE_data.csv; opts detectImportOptions(filename); data readmatrix(filename, opts); % 假设列顺序为序号/时间、幅度、信号强度、其他 t data(:, 1); % 时间或序号 amp data(:, 2); % 幅度单位通常是dB AE S data(:, 3); % 信号强度单位需确认如果文件里有表头readmatrix默认会把它跳过去但如果表头有多行就会把第二行当数据读进来。稳妥一点的做法是raw readcell(filename, FileType, text); % 找到表头行的行号再用readmatrix指定起始行 headerRow find(contains(raw(:,1), Time), 1); data readmatrix(filename, NumHeaderLines, headerRow);这一步看似无关紧要但后面计算HI时N是从1开始的自然数如果混入表头字符串整段代码都会崩溃所以值得多花10秒检查。2.2 信号强度单位陷阱dB和线性别搞混声发射系统的“幅度”通常用dB AE表示但“信号强度”一般已经是线性值。个别系统导出时会把信号强度也转成dB此时需要对数反变换才能用于强度分析。拿常见公式来说如果系统手册写着“Signal Strength (dB) 20*log10(线性值)”那线性信号强度就应该这样恢复linear_S 10 .^ (S_dB / 20);但这里有一个需要特别小心的问题声发射系统的dB基准不一定是同一个电压基准。每个厂家对0 dB的定义可能不同比如1 μV或1 μVpp所以如果你要复现一篇论文最好从论文的实验部分找到采集设备型号再按手册确认。我自己就遇到过因为没注意基准导致Sr曲线整体偏移约30 dB对照不上原图的情况。2.3 数据清洗剔掉机械噪声和电脉冲实测声发射信号里除了材料损伤产生的真实AE事件还有夹具摩擦、电磁干扰、液压系统噪声等。这些伪事件如果不剔掉会直接影响信号强度序列的统计尤其会污染Sr——最强的几个事件里只要混入一个电脉冲Sr就会跳得很高。常见清洗手段按幅度过滤选择只保留幅度高于背景噪声比如40 dB AE的事件按频率过滤保留主频在典型AE频段100 kHz - 1000 kHz的事件按时段过滤有些机械噪声在特定加载阶段集中出现可以手动划定有效时间窗。% 按幅度阈值过滤阈值需要结合系统噪声底看 threshold 45; % dB valid amp threshold; t t(valid); amp amp(valid); S S(valid);注意这里的过滤阈值不能乱拍最好先画一个幅度-时间散点图看看背景噪声水平。我一般直接找加载前那段纯环境噪声的最大幅度再往上加3~6 dB作为事件触发阈值。3. Matlab核心函数实现H与Sr的计算理论清楚了写代码就不难。下面直接给出两个独立的函数输入输出都做了清晰的说明。3.1 computeHI函数这里我按ASTM惯例实现k值做成可配置参数function HI computeHI(S, varargin) % 计算声发射强度分析中的历史指数 Hi % 输入 % S - 信号强度列向量按时间/到达顺序排列 % 可选参数 % kMode - astm(默认不变) 或自定义数值 % 输出 % HI - 与S等长的历史指数序列 % % 说明HI反映了累积信号强度曲线的局部抬升程度 % 具体定义参考ASTM声发射材料评价相关标准。 S S(:); N length(S); CSS cumsum(S); % 累积信号强度 HI ones(N, 1); % 默认最小值为1 % 解析可选参数 kMode astm; if ~isempty(varargin) kMode varargin{1}; end for t 1:N n t; % 当前事件数 % 根据当前事件数确定 k if isnumeric(kMode) k min(kMode, n); % 自定义固定k else % 默认按ASTM分段取值 if n 50 k 0; elseif n 200 k 25; elseif n 500 k 75; else k 125; end end % 防止k过大的边界情况 if k n HI(t) 1; continue; end % 核心公式 numerator CSS(t) - CSS(max(k,0)); denominator CSS(t); HI(t) (n / (n - k)) * (numerator / denominator); end end这里有一个在代码里体现得不明显、但很容易被忽略的边界问题当k0时“去掉前k个事件”等于不去除任何事件CSS(max(k,0))取CSS(0)时在Matlab里会索引越界所以我用了max(k,0)但CSS在Matlab里索引从1开始所以更严谨的写法是if k 0 prefix CSS(k); else prefix 0; end numerator CSS(t) - prefix;用max(k,0)并不能真正解决问题因为CSS(0)仍然非法。让我修正这个细节实际代码里应该用上面这个if判断。3.2 computeSr函数function Sr computeSr(S, varargin) % 计算声发射强度分析中的严重率 Sr % 输入 % S - 信号强度列向量 % J - 参与平均的事件数量默认min(10,N) % 输出 % Sr - 与S等长的严重率序列 % % 说明Sr取前J个最大事件信号强度的平均值可以作为损伤严重度的相对判据。 S S(:); N length(S); Sr zeros(N, 1); % 默认J min(10, N) J min(10, N); if ~isempty(varargin) J min(varargin{1}, N); if J 0 error(J必须为正整数); end end for t 1:N if t J Sr(t) mean(S(1:t)); else % 当前时刻所有数据中的前J个最大值 sortedVals maxk(S(1:t), J); Sr(t) mean(sortedVals); end end end这个实现是逐步滑动的复杂度O(N^2 log J)对N几千没问题。但如果你的导出文件里有几万甚至几十万个事件这个写法会非常慢。我后面会单独说优化方案。3.3 效率优化滑动窗口和数据量过大的处理声发射试验一跑几个小时事件数轻松破万。N20000时上面的computeSr每次循环都要对S(1:t)做一次maxk运行时间会让人怀疑人生。优化思路是这样的当N特别大时不需要每一步都精确计算基于全历史最大值的Sr。可以这样做方案一降采样。只对关键时间节点计算Sr比如每隔2~5个事件计算一次再用interp1做线性插值。画图时人眼根本看不出差别。方案二利用累积最大值堆结构。Matlab没有现成的最大堆但可以用accumarray和排序替代。方案三等时间间隔分箱。把时间轴切成500~1000个区间每个区间内只保留最大信号强度之后在区间级别上计算Sr。实用性最强的是方案一实现非常简单% 降采样计算Sr每step个点计算一次之后插值 step 5; idx 1:step:N; Sr_sampled arrayfun((x) mean(maxk(S(1:x), min(10,x))), idx); Sr_full interp1(idx, Sr_sampled, 1:N, linear, extrap);这样速度提升了step倍精度损失在可接受范围。实际复现论文图时如果原图本身是连续曲线降采样后插值根本看不出来。4. 绘图复现从计算值到论文配图4.1 基础绘图HI和Sr随事件序号/时间的变化计算完HI和Sr后第一张图通常是横轴为时间或事件序号纵轴为HI和Sr。一种常用做法是把两个指标放在双y轴图里因为它们的量纲完全不同HI是接近1~20的比值量级Sr可能是数千到数万。% 假设已经计算出HI和Sr且t为时间向量 figure(Color, w, Position, [100 100 800 420]); yyaxis left plot(t, HI, b-, LineWidth, 1.2); ylabel(Hi (历史指数), FontName, Helvetica, FontSize, 11); yyaxis right plot(t, Sr, r-, LineWidth, 1.2); ylabel(Sr (严重率), FontName, Helvetica, FontSize, 11); xlabel(试验时间 (s), FontName, Helvetica, FontSize, 11); set(gca, FontName, Helvetica, FontSize, 10); grid on;画这种图我会加一句提醒双y轴虽然直观但在正式论文里容易引起误解尤其当左右两个轴的刻度范围没有对齐时曲线在视觉上的“突变点”可能被放大或缩小。如果审稿人较真建议分上下两个子图共用x轴figure(Color, w, Position, [100 100 800 680]); subplot(2,1,1); plot(t, HI, b-, LineWidth, 1.2); ylabel(Hi); grid on; subplot(2,1,2); plot(t, Sr, r-, LineWidth, 1.2); xlabel(试验时间 (s)); ylabel(Sr); grid on; % 上下子图共用x轴保证对齐 linkaxes(findobj(gcf, Type, axes), x);linkaxes这行是关键不加上下两个子图的x轴范围可能不一致肉眼对比时会误判两个指标变化的时间对应关系。4.2 累积信号强度CSS曲线作为辅助参考复现论文图时我通常会把累积信号强度同时画出来因为它和HI/Sr有直接关系但CSS是对数增长的直接画在同一个坐标系里会挤压其他曲线。更好的方式是双y轴布局左边放CSS对数刻度右边放HI或Sr。这样一眼就能看出CSS斜率变化跟HI抬升的关系。figure(Color, w, Position, [100 100 800 420]); CSS cumsum(S); yyaxis left semilogy(t, CSS, k-, LineWidth, 1.0); ylabel(累积信号强度 (对数刻度)); yyaxis right plot(t, HI, b-, LineWidth, 1.2); ylabel(Hi); xlabel(试验时间 (s)); grid on;做这个图的好处是CSS垂直台阶跳变的位置通常对应HI的局部峰值两者互相印证能排除偶然的个别大信号导致的误判。4.3 多通道数据循环绘图与图例实际声发射试验往往不只一个传感器四个通道是比较常见的。把四组HI/Sr画在一起时最省事的方式是循环但我用过很多次之后有个教训刷参数时容易把通道间的颜色搞混。配色上建议用Matlab自带的parula或lines不要每轮循环里用默认顺序色循环第二次起会重复可以提前定义好颜色矩阵colors lines(channelNum); % 取通道数不同颜色 figure(Color, w, Position, [100 100 900 420]); hold on; for ch 1:channelNum % ch_data 读取对应通道数据 HI_ch computeHI(ch_data); plot(t_ch, HI_ch, Color, colors(ch, :), LineWidth, 1.2, ... DisplayName, sprintf(CH%d, ch)); end hold off; legend(Location, northwest); xlabel(试验时间 (s)); ylabel(Hi); grid on;legend的显示名称用DisplayName属性来设置比用legend({CH1,CH2,...})更灵活尤其当循环里有跳过某些通道的逻辑时不会出现图例比曲线多的情况。4.4 导出高清图论文投稿级别的输出设置MATLAB导出图片到Word或论文时容易遇到两个问题字体不对和分辨率不够。我最常用的导出方式是exportgraphics能直接控制分辨率exportgraphics(gcf, HI_Sr_analysis.png, Resolution, 600);如果投期刊还需要矢量图建议输出为PDF或EPS格式exportgraphics(gcf, HI_Sr_analysis.pdf, ContentType, vector);字体方面中文字体在很多期刊模板里会出问题。如果图表里只保留英文标签用Helvetica或Arial即可。如果必须写中文最好先把图导出成PDF再在排版软件里补上文字比在Matlab里硬调字体省事得多。5. 复现过程中最容易翻车的几个细节5.1 HI的前端异常为什么第一段曲线接近一条平线HI在事件数小于等于50时k0公式退化为n/n×(CSS/CSS)1所以曲线前端一定是一条1.0的平线。这不是bug而是定义使然——当样本量太少时历史指数没有统计意义。但有些论文里的HI图没有这截平台因为他们画图的横轴是从50或某个阈值之后开始的。复现时如果原图的横轴起点不是0要注意检查对方是否做了截断。如果你觉得平台太长不好看可以设置绘图时只显示HI 1的数据idx_plot HI 1; plot(t(idx_plot), HI(idx_plot), b-);但这种处理要如实说明否则会掩盖试验初期的真实信息。5.2 Sr对个别异常事件高度敏感Sr取的是前J个最大信号强度的平均值这带来一个实际麻烦只要出现一个超强噪声事件Sr会在很长一段时间内维持在高位因为它始终被排在前J个最大值里一直到后面出现更多的大事件把它挤出去。如果你复现的图里Sr曲线出现明显的“阶梯跳水”大多不是算法错误而是某个异常事件被更后面的事件从top J中挤掉了。此时做数据清洗时要有意识地检查最大信号强度的那几个事件来源% 找出信号强度最大的10个事件 [~, idxTop] maxk(S, 10); % 打印这些事件的到达时间和幅度方便回溯 table(idxTop, t(idxTop), amp(idxTop), S(idxTop));从工程角度说这是声发射强度分析的一个固有特性也是为什么Sr要与HI配合使用HI看整体趋势Sr看极端事件。5.3 时间字段对齐加载阶段切分有些试验在加载过程中分成多个阶段比如加-保-卸-加声发射事件不是均匀分布在时间轴上。如果你直接按“事件序号”作为横轴计算HI会丢失时间信息但如果你按“时间”作为横轴需要注意事件密度不均匀导致曲线前段稀疏、后段密集。我建议的默认做法是先画“事件序号-时间”图判断时间轴上的事件分布是否均匀。如果存在明显的加载停顿要按加载阶段分别计算HI和Sr而不是把整个试验数据一口气跑完。% 简单分组按时间间隔中是否存在超过threshold_s的空窗来切分数据段 gaps diff(t); gapIdx find(gaps 5); % 超过5秒无事件视为阶段分隔 segStart [1; gapIdx 1]; segEnd [gapIdx; length(t)]; % 对每段分别调用computeHI和computeSr如果不做分段保载阶段的稀疏事件会把整体曲线压得很平遗漏真正的损伤演化信息。5.4 与论文原图对照时注意横轴坐标到底是时间还是归一化百分比这个问题非常实在。很多论文的声发射强度分析图横轴写的是“加载时间”或“加载百分比”但他们的数据实际是按事件序号画的只是把事件序号按比例换算成了时间或百分比。如果直接按你自己的时间轴画曲线形状相同但横轴刻度会对不上。对策是对照原图数据点的稀疏程度和曲线上的突变节点确认映射关系。具体做法可以在两条曲线的突变点位置取几个时间坐标做线性回归判断是否存在恒定的缩放因子。% 假设突变点处你的时间坐标是 t1, t2, t3 % 原图横轴坐标是 x1, x2, x3 % 拟合 x a * t b p polyfit([t1 t2 t3], [x1 x2 x3], 1); disp([缩放因子 a , num2str(p(1)), , 偏移 b , num2str(p(2))]);如果线性拟合残差很小说明你的时间轴与原图就是线性关系画图时直接用这个缩放关系转换即可。6. 一个完整案例FRP拉伸试件的强度分析6.1 数据情况和目标我用一组玻纤增强复合材料板的拉伸加载AE数据做演示。采集系统是PCI-2PAC公司事件记录包括幅度、能量、持续时间、信号强度。加载过程大约300秒事件总数约4200个中间有两次5秒的保载换挡整体事件密度在加载后期显著上升。复现目标是画出“信号强度分析图”包括CSS、HI和Sr三条曲线并判断损伤大致在哪个时间段进入加速阶段。6.2 逐步执行结果先按第2节的方法导入数据做45 dB的幅度阈值过滤剩约3917个事件。计算HI和SrHI computeHI(S_valid); Sr computeSr(S_valid, 10); % 取前10个最大值平均然后画出三合一图CSS对数、HI、Sr双y轴加CSS。从结果图上看CSS曲线在前180秒近似线性增长属于典型的基体微裂纹稳定扩展期从约190秒开始CSS斜率突然变陡同时HI从接近1的数值在10秒内跳到3.5左右说明有较大的损伤事件集中出现Sr在约195秒达到第一个峰值对应一次较大的开裂事件之后下降但没有回到峰值前水平在240秒后HI出现第二次抬升Sr继续波动上升说明损伤进入更剧烈的层板脱粘阶段。这些现象和试件加载后期刚度下降的节点基本吻合说明强度分析指标能有效捕捉到单一事件计数看不到的损伤转折。6.3 阈值设置的参考用HI和Sr做损伤严重度判断时不同材料的阈值差异很大。对我测过的FRP材料经验值大概是HI稳定在1.0~1.2之间时可视为正常运行状态HI超过2.0通常对应可见的基体开裂聚集HI上升到4.0以上并伴随Sr大幅升高往往接近宏观分层或纤维断裂。但这些数值只能作为相对参考不能直接套到金属材料或者其他工艺状态的复合材料上。不同采集系统、不同阈值、不同探头增益下绝对数值都会变化。正确打开方式是先用同一批次中无损伤或低损伤试件建立基线再拿基线数据的HI/Sr统计分布做判据。我在实际项目中也是这样做的比直接搬文献里的阈值可靠得多。6.4 从指标到结论的完整逻辑最后总结一下我在这个项目里建立的分析链路后续遇到类似声发射数据都可以按这个顺序来快速浏览原始数据剔除明显噪声段计算CSS判断全局增长趋势计算HI定位累积曲线异常抬升的起点计算Sr确认这些异常点是否伴随强事件把HI和Sr的突变时间与加载过程中的力学参数位移、载荷、刚度对照确认物理意义生成报告图时优先选择上下子图或双y轴图把CSS和HI/Sr画在一起提供完整证据链。这样得到的分析结果既不会因为单看计数指标而漏掉损伤转折也不会因为个别异常事件而误判严重度比单纯读峰值幅度或累加能量都更能反映结构损伤的真实演化过程。 
长河君 · 2026/9/9 15:28:52 · 阅读 1.2k 声发射强度分析HI与Sr计算及Matlab复现指南
最近在复现一篇FRP复合材料加载试验的声发射数据分析图卡在一个很基础但特别容易绕晕的词上——“强度分析”。论文里用到了两个缩写HI和Sr图表看着很高级一条曲线平缓另一条在某处突然抬头。但等我翻到方法部分只有短短一段公式和一句“按ASTM标准计算”细节全靠猜。折腾了两天试了好几种写法总算是把曲线和原文对上号了。这篇博客就是把我摸索出来的完整过程写下来从指标的物理意义、Matlab实现、到绘图复现每一步都给出可以直接跑的代码和避坑提醒。无论是做复合材料损伤检测、金属裂纹监测还是纯想搞懂声发射强度分析的新手都能按这份记录复现出自己的HI和Sr图。1. 先搞清楚HI和Sr到底在算什么很多教材把声发射信号处理讲得玄乎其实落到实际上强度分析的核心思路非常直白声发射事件不是均匀发生的强度突变往往对应损伤模式的改变。1.1 信号强度为什么比计数可靠声发射系统通常会记录好几个参数到达时间、峰值幅度单位通常是dB AE、能量、持续时间、上升时间。论文里如果出现“signal strength”指的往往是信号强度这个概念单位一般是伏·秒或皮欧·米也可以理解成信号包络的积分面积。一个常见误区是拿“事件计数”当主要指标。计数的问题是一次大裂纹扩展产生的信号和一百次微小基体开裂产生的信号计数值可能差别很大但对材料损伤的实际贡献完全不同。信号强度保留了信号的“大小”信息所以强度分析更能反映真实损伤。简单说信号强度序列就是按时间顺序排列的S_ii1,2,...,N后面的所有计算都围绕这个序列展开。1.2 HI历史指数找累积曲线上的“抬升点”HI这个名字有不同的展开写法常见的是Histogram Index或Historical Index中文翻译成历史指数。它干的事情是考察截止到当前时刻累积信号强度曲线的“局部斜率”是否比历史平均更高。公式在ASTM标准里的定义大概是[ H(t) \frac{N}{N - k} \times \frac{\sum_{ik1}^{N} S_i}{\sum_{i1}^{N} S_i} ]其中N是截止当前时刻的事件总数S_i是第i个事件按到达时间排列的信号强度值k是偏移量取值跟N有关。这个公式的意思是把前面k个事件剔除后剩下这部分信号强度占整体累积的比例再乘以一个缩放系数。如果损伤进入剧烈阶段尾部信号强度占比突然变大HI就会显著上升。k的取值每个标准版本略有出入我用的惯例如下当前事件数 N偏移 kN ≤ 50051 ≤ N ≤ 20025201 ≤ N ≤ 50075N 500125需要说明的是ASTM标准本身和后来各路论文里的实现不完全一样做复现时第一件事是看你参考的论文怎么定义的k。后面代码里我会把k做成一个可配置参数方便切换。1.3 Sr严重率盯着最大信号强度的尾巴SrSeverity严重率的定义比HI简单它衡量的是当前时刻之前最大的J个信号强度的平均值[ S_r \frac{1}{J} \sum_{j1}^{J} S_{max,j} ]J在不同文献里有取10的也有取事件总数一定比例的。常见的简化实现是取J min(10, N)也就是最多10个最大信号强度求平均。如果当前事件数太少就取全部。从物理上讲HI反映的是累积历史中“增量异常”而Sr反映的是“最强事件有多强”。二者配合使用既能找到剧烈活动的起点也能定性判断活动猛烈程度。很多标准里会用Sr作为损伤严重度的直接判据单看数值大小决定是否达到临界状态。2. 数据准备光有采集文件还不够2.1 从采集系统导出什么样的数据声发射采集系统比如Physical Acoustics的PAC系统、Vallen等一般都能导出ASCII文本或CSV文件。列名可能五花八门但关键的几列通常包括时间或到达序号Time/Segment、峰值幅度Amplitude、信号强度Signal Strength。有个很容易踩的坑很多系统导出的是带行号的自由格式文本表头在第二行列之间用制表符或逗号混用。直接readmatrix可能报错建议先看一眼文件头几行再决定导入方式。我建议的第一步不是急着写计算函数而是把数据读进来确认列对得上。以最典型的三列CSV为例% 读取声发射导出数据 filename AE_data.csv; opts detectImportOptions(filename); data readmatrix(filename, opts); % 假设列顺序为序号/时间、幅度、信号强度、其他 t data(:, 1); % 时间或序号 amp data(:, 2); % 幅度单位通常是dB AE S data(:, 3); % 信号强度单位需确认如果文件里有表头readmatrix默认会把它跳过去但如果表头有多行就会把第二行当数据读进来。稳妥一点的做法是raw readcell(filename, FileType, text); % 找到表头行的行号再用readmatrix指定起始行 headerRow find(contains(raw(:,1), Time), 1); data readmatrix(filename, NumHeaderLines, headerRow);这一步看似无关紧要但后面计算HI时N是从1开始的自然数如果混入表头字符串整段代码都会崩溃所以值得多花10秒检查。2.2 信号强度单位陷阱dB和线性别搞混声发射系统的“幅度”通常用dB AE表示但“信号强度”一般已经是线性值。个别系统导出时会把信号强度也转成dB此时需要对数反变换才能用于强度分析。拿常见公式来说如果系统手册写着“Signal Strength (dB) 20*log10(线性值)”那线性信号强度就应该这样恢复linear_S 10 .^ (S_dB / 20);但这里有一个需要特别小心的问题声发射系统的dB基准不一定是同一个电压基准。每个厂家对0 dB的定义可能不同比如1 μV或1 μVpp所以如果你要复现一篇论文最好从论文的实验部分找到采集设备型号再按手册确认。我自己就遇到过因为没注意基准导致Sr曲线整体偏移约30 dB对照不上原图的情况。2.3 数据清洗剔掉机械噪声和电脉冲实测声发射信号里除了材料损伤产生的真实AE事件还有夹具摩擦、电磁干扰、液压系统噪声等。这些伪事件如果不剔掉会直接影响信号强度序列的统计尤其会污染Sr——最强的几个事件里只要混入一个电脉冲Sr就会跳得很高。常见清洗手段按幅度过滤选择只保留幅度高于背景噪声比如40 dB AE的事件按频率过滤保留主频在典型AE频段100 kHz - 1000 kHz的事件按时段过滤有些机械噪声在特定加载阶段集中出现可以手动划定有效时间窗。% 按幅度阈值过滤阈值需要结合系统噪声底看 threshold 45; % dB valid amp threshold; t t(valid); amp amp(valid); S S(valid);注意这里的过滤阈值不能乱拍最好先画一个幅度-时间散点图看看背景噪声水平。我一般直接找加载前那段纯环境噪声的最大幅度再往上加3~6 dB作为事件触发阈值。3. Matlab核心函数实现H与Sr的计算理论清楚了写代码就不难。下面直接给出两个独立的函数输入输出都做了清晰的说明。3.1 computeHI函数这里我按ASTM惯例实现k值做成可配置参数function HI computeHI(S, varargin) % 计算声发射强度分析中的历史指数 Hi % 输入 % S - 信号强度列向量按时间/到达顺序排列 % 可选参数 % kMode - astm(默认不变) 或自定义数值 % 输出 % HI - 与S等长的历史指数序列 % % 说明HI反映了累积信号强度曲线的局部抬升程度 % 具体定义参考ASTM声发射材料评价相关标准。 S S(:); N length(S); CSS cumsum(S); % 累积信号强度 HI ones(N, 1); % 默认最小值为1 % 解析可选参数 kMode astm; if ~isempty(varargin) kMode varargin{1}; end for t 1:N n t; % 当前事件数 % 根据当前事件数确定 k if isnumeric(kMode) k min(kMode, n); % 自定义固定k else % 默认按ASTM分段取值 if n 50 k 0; elseif n 200 k 25; elseif n 500 k 75; else k 125; end end % 防止k过大的边界情况 if k n HI(t) 1; continue; end % 核心公式 numerator CSS(t) - CSS(max(k,0)); denominator CSS(t); HI(t) (n / (n - k)) * (numerator / denominator); end end这里有一个在代码里体现得不明显、但很容易被忽略的边界问题当k0时“去掉前k个事件”等于不去除任何事件CSS(max(k,0))取CSS(0)时在Matlab里会索引越界所以我用了max(k,0)但CSS在Matlab里索引从1开始所以更严谨的写法是if k 0 prefix CSS(k); else prefix 0; end numerator CSS(t) - prefix;用max(k,0)并不能真正解决问题因为CSS(0)仍然非法。让我修正这个细节实际代码里应该用上面这个if判断。3.2 computeSr函数function Sr computeSr(S, varargin) % 计算声发射强度分析中的严重率 Sr % 输入 % S - 信号强度列向量 % J - 参与平均的事件数量默认min(10,N) % 输出 % Sr - 与S等长的严重率序列 % % 说明Sr取前J个最大事件信号强度的平均值可以作为损伤严重度的相对判据。 S S(:); N length(S); Sr zeros(N, 1); % 默认J min(10, N) J min(10, N); if ~isempty(varargin) J min(varargin{1}, N); if J 0 error(J必须为正整数); end end for t 1:N if t J Sr(t) mean(S(1:t)); else % 当前时刻所有数据中的前J个最大值 sortedVals maxk(S(1:t), J); Sr(t) mean(sortedVals); end end end这个实现是逐步滑动的复杂度O(N^2 log J)对N几千没问题。但如果你的导出文件里有几万甚至几十万个事件这个写法会非常慢。我后面会单独说优化方案。3.3 效率优化滑动窗口和数据量过大的处理声发射试验一跑几个小时事件数轻松破万。N20000时上面的computeSr每次循环都要对S(1:t)做一次maxk运行时间会让人怀疑人生。优化思路是这样的当N特别大时不需要每一步都精确计算基于全历史最大值的Sr。可以这样做方案一降采样。只对关键时间节点计算Sr比如每隔2~5个事件计算一次再用interp1做线性插值。画图时人眼根本看不出差别。方案二利用累积最大值堆结构。Matlab没有现成的最大堆但可以用accumarray和排序替代。方案三等时间间隔分箱。把时间轴切成500~1000个区间每个区间内只保留最大信号强度之后在区间级别上计算Sr。实用性最强的是方案一实现非常简单% 降采样计算Sr每step个点计算一次之后插值 step 5; idx 1:step:N; Sr_sampled arrayfun((x) mean(maxk(S(1:x), min(10,x))), idx); Sr_full interp1(idx, Sr_sampled, 1:N, linear, extrap);这样速度提升了step倍精度损失在可接受范围。实际复现论文图时如果原图本身是连续曲线降采样后插值根本看不出来。4. 绘图复现从计算值到论文配图4.1 基础绘图HI和Sr随事件序号/时间的变化计算完HI和Sr后第一张图通常是横轴为时间或事件序号纵轴为HI和Sr。一种常用做法是把两个指标放在双y轴图里因为它们的量纲完全不同HI是接近1~20的比值量级Sr可能是数千到数万。% 假设已经计算出HI和Sr且t为时间向量 figure(Color, w, Position, [100 100 800 420]); yyaxis left plot(t, HI, b-, LineWidth, 1.2); ylabel(Hi (历史指数), FontName, Helvetica, FontSize, 11); yyaxis right plot(t, Sr, r-, LineWidth, 1.2); ylabel(Sr (严重率), FontName, Helvetica, FontSize, 11); xlabel(试验时间 (s), FontName, Helvetica, FontSize, 11); set(gca, FontName, Helvetica, FontSize, 10); grid on;画这种图我会加一句提醒双y轴虽然直观但在正式论文里容易引起误解尤其当左右两个轴的刻度范围没有对齐时曲线在视觉上的“突变点”可能被放大或缩小。如果审稿人较真建议分上下两个子图共用x轴figure(Color, w, Position, [100 100 800 680]); subplot(2,1,1); plot(t, HI, b-, LineWidth, 1.2); ylabel(Hi); grid on; subplot(2,1,2); plot(t, Sr, r-, LineWidth, 1.2); xlabel(试验时间 (s)); ylabel(Sr); grid on; % 上下子图共用x轴保证对齐 linkaxes(findobj(gcf, Type, axes), x);linkaxes这行是关键不加上下两个子图的x轴范围可能不一致肉眼对比时会误判两个指标变化的时间对应关系。4.2 累积信号强度CSS曲线作为辅助参考复现论文图时我通常会把累积信号强度同时画出来因为它和HI/Sr有直接关系但CSS是对数增长的直接画在同一个坐标系里会挤压其他曲线。更好的方式是双y轴布局左边放CSS对数刻度右边放HI或Sr。这样一眼就能看出CSS斜率变化跟HI抬升的关系。figure(Color, w, Position, [100 100 800 420]); CSS cumsum(S); yyaxis left semilogy(t, CSS, k-, LineWidth, 1.0); ylabel(累积信号强度 (对数刻度)); yyaxis right plot(t, HI, b-, LineWidth, 1.2); ylabel(Hi); xlabel(试验时间 (s)); grid on;做这个图的好处是CSS垂直台阶跳变的位置通常对应HI的局部峰值两者互相印证能排除偶然的个别大信号导致的误判。4.3 多通道数据循环绘图与图例实际声发射试验往往不只一个传感器四个通道是比较常见的。把四组HI/Sr画在一起时最省事的方式是循环但我用过很多次之后有个教训刷参数时容易把通道间的颜色搞混。配色上建议用Matlab自带的parula或lines不要每轮循环里用默认顺序色循环第二次起会重复可以提前定义好颜色矩阵colors lines(channelNum); % 取通道数不同颜色 figure(Color, w, Position, [100 100 900 420]); hold on; for ch 1:channelNum % ch_data 读取对应通道数据 HI_ch computeHI(ch_data); plot(t_ch, HI_ch, Color, colors(ch, :), LineWidth, 1.2, ... DisplayName, sprintf(CH%d, ch)); end hold off; legend(Location, northwest); xlabel(试验时间 (s)); ylabel(Hi); grid on;legend的显示名称用DisplayName属性来设置比用legend({CH1,CH2,...})更灵活尤其当循环里有跳过某些通道的逻辑时不会出现图例比曲线多的情况。4.4 导出高清图论文投稿级别的输出设置MATLAB导出图片到Word或论文时容易遇到两个问题字体不对和分辨率不够。我最常用的导出方式是exportgraphics能直接控制分辨率exportgraphics(gcf, HI_Sr_analysis.png, Resolution, 600);如果投期刊还需要矢量图建议输出为PDF或EPS格式exportgraphics(gcf, HI_Sr_analysis.pdf, ContentType, vector);字体方面中文字体在很多期刊模板里会出问题。如果图表里只保留英文标签用Helvetica或Arial即可。如果必须写中文最好先把图导出成PDF再在排版软件里补上文字比在Matlab里硬调字体省事得多。5. 复现过程中最容易翻车的几个细节5.1 HI的前端异常为什么第一段曲线接近一条平线HI在事件数小于等于50时k0公式退化为n/n×(CSS/CSS)1所以曲线前端一定是一条1.0的平线。这不是bug而是定义使然——当样本量太少时历史指数没有统计意义。但有些论文里的HI图没有这截平台因为他们画图的横轴是从50或某个阈值之后开始的。复现时如果原图的横轴起点不是0要注意检查对方是否做了截断。如果你觉得平台太长不好看可以设置绘图时只显示HI 1的数据idx_plot HI 1; plot(t(idx_plot), HI(idx_plot), b-);但这种处理要如实说明否则会掩盖试验初期的真实信息。5.2 Sr对个别异常事件高度敏感Sr取的是前J个最大信号强度的平均值这带来一个实际麻烦只要出现一个超强噪声事件Sr会在很长一段时间内维持在高位因为它始终被排在前J个最大值里一直到后面出现更多的大事件把它挤出去。如果你复现的图里Sr曲线出现明显的“阶梯跳水”大多不是算法错误而是某个异常事件被更后面的事件从top J中挤掉了。此时做数据清洗时要有意识地检查最大信号强度的那几个事件来源% 找出信号强度最大的10个事件 [~, idxTop] maxk(S, 10); % 打印这些事件的到达时间和幅度方便回溯 table(idxTop, t(idxTop), amp(idxTop), S(idxTop));从工程角度说这是声发射强度分析的一个固有特性也是为什么Sr要与HI配合使用HI看整体趋势Sr看极端事件。5.3 时间字段对齐加载阶段切分有些试验在加载过程中分成多个阶段比如加-保-卸-加声发射事件不是均匀分布在时间轴上。如果你直接按“事件序号”作为横轴计算HI会丢失时间信息但如果你按“时间”作为横轴需要注意事件密度不均匀导致曲线前段稀疏、后段密集。我建议的默认做法是先画“事件序号-时间”图判断时间轴上的事件分布是否均匀。如果存在明显的加载停顿要按加载阶段分别计算HI和Sr而不是把整个试验数据一口气跑完。% 简单分组按时间间隔中是否存在超过threshold_s的空窗来切分数据段 gaps diff(t); gapIdx find(gaps 5); % 超过5秒无事件视为阶段分隔 segStart [1; gapIdx 1]; segEnd [gapIdx; length(t)]; % 对每段分别调用computeHI和computeSr如果不做分段保载阶段的稀疏事件会把整体曲线压得很平遗漏真正的损伤演化信息。5.4 与论文原图对照时注意横轴坐标到底是时间还是归一化百分比这个问题非常实在。很多论文的声发射强度分析图横轴写的是“加载时间”或“加载百分比”但他们的数据实际是按事件序号画的只是把事件序号按比例换算成了时间或百分比。如果直接按你自己的时间轴画曲线形状相同但横轴刻度会对不上。对策是对照原图数据点的稀疏程度和曲线上的突变节点确认映射关系。具体做法可以在两条曲线的突变点位置取几个时间坐标做线性回归判断是否存在恒定的缩放因子。% 假设突变点处你的时间坐标是 t1, t2, t3 % 原图横轴坐标是 x1, x2, x3 % 拟合 x a * t b p polyfit([t1 t2 t3], [x1 x2 x3], 1); disp([缩放因子 a , num2str(p(1)), , 偏移 b , num2str(p(2))]);如果线性拟合残差很小说明你的时间轴与原图就是线性关系画图时直接用这个缩放关系转换即可。6. 一个完整案例FRP拉伸试件的强度分析6.1 数据情况和目标我用一组玻纤增强复合材料板的拉伸加载AE数据做演示。采集系统是PCI-2PAC公司事件记录包括幅度、能量、持续时间、信号强度。加载过程大约300秒事件总数约4200个中间有两次5秒的保载换挡整体事件密度在加载后期显著上升。复现目标是画出“信号强度分析图”包括CSS、HI和Sr三条曲线并判断损伤大致在哪个时间段进入加速阶段。6.2 逐步执行结果先按第2节的方法导入数据做45 dB的幅度阈值过滤剩约3917个事件。计算HI和SrHI computeHI(S_valid); Sr computeSr(S_valid, 10); % 取前10个最大值平均然后画出三合一图CSS对数、HI、Sr双y轴加CSS。从结果图上看CSS曲线在前180秒近似线性增长属于典型的基体微裂纹稳定扩展期从约190秒开始CSS斜率突然变陡同时HI从接近1的数值在10秒内跳到3.5左右说明有较大的损伤事件集中出现Sr在约195秒达到第一个峰值对应一次较大的开裂事件之后下降但没有回到峰值前水平在240秒后HI出现第二次抬升Sr继续波动上升说明损伤进入更剧烈的层板脱粘阶段。这些现象和试件加载后期刚度下降的节点基本吻合说明强度分析指标能有效捕捉到单一事件计数看不到的损伤转折。6.3 阈值设置的参考用HI和Sr做损伤严重度判断时不同材料的阈值差异很大。对我测过的FRP材料经验值大概是HI稳定在1.0~1.2之间时可视为正常运行状态HI超过2.0通常对应可见的基体开裂聚集HI上升到4.0以上并伴随Sr大幅升高往往接近宏观分层或纤维断裂。但这些数值只能作为相对参考不能直接套到金属材料或者其他工艺状态的复合材料上。不同采集系统、不同阈值、不同探头增益下绝对数值都会变化。正确打开方式是先用同一批次中无损伤或低损伤试件建立基线再拿基线数据的HI/Sr统计分布做判据。我在实际项目中也是这样做的比直接搬文献里的阈值可靠得多。6.4 从指标到结论的完整逻辑最后总结一下我在这个项目里建立的分析链路后续遇到类似声发射数据都可以按这个顺序来快速浏览原始数据剔除明显噪声段计算CSS判断全局增长趋势计算HI定位累积曲线异常抬升的起点计算Sr确认这些异常点是否伴随强事件把HI和Sr的突变时间与加载过程中的力学参数位移、载荷、刚度对照确认物理意义生成报告图时优先选择上下子图或双y轴图把CSS和HI/Sr画在一起提供完整证据链。这样得到的分析结果既不会因为单看计数指标而漏掉损伤转折也不会因为个别异常事件而误判严重度比单纯读峰值幅度或累加能量都更能反映结构损伤的真实演化过程。