SHPB实验数据处理全流程：MATLAB自动化实现方案

简介：面向材料动态力学性能研究与军事工程领域，这份PDF文献围绕霍普金森杆（SHPB）实验数据处理程序展开，内容涵盖实验基本原理、入射波/反射波/透射波分离难点及基于VC++的程序设计思路，适合需要了解SHPB数据处理方法或开发相关工具的研究人员与工程技术人员。资源为单份PDF文档，共1个文件，压缩包约176KB，轻量易用。目前已有438人学习/下载。文中详细介绍了一维假定与应力应变均匀假定，并给出三波公式推导，还展示了数据输入、数据处理、波形处理、曲线显示等模块的设计框架，读者可据此快速理解完整的数据处理流程与程序实现方案。 做SHPB实验的人都知道，整套装置搭起来不算太难，真正的功夫全在数据处理上。从示波器里导出的原始波形，到最终能写进论文里的应力-应变曲线，中间隔着滤波、对时、截取、计算一大堆步骤，每一步都可能埋着坑。我最初做霍普金森杆实验时，完全靠Origin手工处理，一条曲线折腾半小时是常事，遇到波形稍微难看一点，整个人都要崩溃。后来咬咬牙，把这套数据处理程序完整写了出来，从读取CSV到输出应力-应变曲线全自动跑通，这才算是把实验效率真正提了上来。这篇博文就完整梳理一下这套程序的设计思路和实现细节，希望能给同样被数据处理折磨的朋友一点参考。

1. 整体方案设计与思路拆解

1.1 核心需求与难点分析

霍普金森杆（SHPB）实验数据处理程序，本质上解决的是“把原始电压信号变成材料动态力学性能参数”这个问题。实验时，撞击杆以一定速度撞击入射杆，产生的应力波经过入射杆、试件、透射杆，由贴在杆上的应变片记录下三组波形信号：入射波、反射波和透射波。我们需要做的，就是对这些信号进行平滑、对时、截取，再代进经典的三波法公式，算出声学应变率、应力和应变，最终绘制出动态应力-应变曲线。

这个过程中有几个绕不开的难点：

• 原始信号噪声大，尤其是应变片信号，基线漂移和毛刺非常严重，不处理直接算会得到一堆荒谬的数字。
• 三组波形的时间轴必须严格对齐，入射波起跳点和透射波起跳点之间的时间差，直接决定了应力平衡的判定，差一个采样点结果就完全不一样。
• 特征点的选取是个精细活，起跳点定在哪、平台段取多宽、截断位置怎么选，这些参数直接决定最终曲线的形态。

我当时写程序时给自己定了一条原则：能自动判断的绝不手动干预，但必须保留手动微调的人工接口。这样既保证了批量处理的效率，又兼顾了单条数据的精细修正。

1.2 程序模块划分与数据流设计

整套程序采用模块化结构，按数据流方向拆分为五个核心模块：数据读取模块、信号预处理模块、特征点识别模块、三波法计算模块和结果输出模块。这种划分方式的好处是，每个模块的输入输出都是明确的中间数据，调试时可以单独验证每一段的正确性。

程序的主控流程大概是：读取原始CSV文件，解析出时间轴和入射波、反射波、透射波三组电压信号，然后依次进行滤波、对时、特征点识别、截取，再代入计算公式得到各项力学参数，最后自动绘制波形图、应力平衡检查图和应力-应变曲线，并把关键参数导出到Excel或者文本文件。

数据流的设计上有一个关键细节，就是所有中间结果都保留在内存里的同时，也允许用户选择性地导出中间数据。我吃过亏，有时候程序跑完发现结果不对，但中间数据已经丢了，只能从头再来。所以实时保存中间结果这个习惯，真是拿踩坑换来的经验。

1.3 工具选型：为什么用MATLAB而不是Python或Origin

选型的时候，我在MATLAB、Python和Origin之间犹豫了很久。Python胜在免费开源，数据科学库丰富，numpy和scipy做信号处理也完全够用。Origin做图表漂亮，交互方便，是实验室很多人的首选。但最终我还是选了MATLAB，原因有三点：

第一，MATLAB的Signal Processing Toolbox在信号处理函数上封装程度高，filtfilt做零相位滤波、findpeaks找峰值、ginput做手动选点，都是开箱即用，开发和调试效率远高于Python需要自己造轮子的方式。

第二，MATLAB在高校力学实验室中仍然是主流工具，很多已发表的SHPB数据处理论文和代码都是MATLAB写的，参考和验证起来非常方便。遇到问题的时候，可以快速找到现成的代码片段对照排查。

第三，MATLAB的App Designer做交互界面门槛不高，可以比较快地给程序包一个图形界面，方便课题组其他不会写代码的同事实操。虽然我自己平时更习惯直接跑脚本，但考虑到程序的可用性，界面还是值得做的。

当然，如果你的实验室已经有Python的完整环境，也没必要非得迁移到MATLAB。核心的算法逻辑和我在下文讲解的流程是通用的，换成Python也就是换一套语法的问题。

2. 核心细节解析与实操要点

2.1 原始信号的平滑与去噪处理

SHPB的原始信号噪声主要来自三个方面：应变片电桥的固有噪声、电磁干扰以及采集系统自身的量化噪声。噪声的幅值正常情况下相对信号本身要小很多，但如果实验时屏蔽没做好或者杆件对中不佳，噪声水平会显著抬升，直接影响后续特征点识别的准确性。

滤波方案上，我测试过滑动平均、移动中值、巴特沃斯低通滤波三种方式，最终采用的是“滑动平均+巴特沃斯零相位低通滤波”的组合策略。滑动平均的窗口宽度取采样点数的0.5%到1%，大小为5到11个点为佳。太大的窗口会把真实的波形前沿抹平，导致起跳点位置偏移。

我实际使用的滤波代码片段如下：

function data_filtered = shpb_filter(raw_data, fs, cutoff_freq) % 滑动平均去高频毛刺，窗口大小为5个点 data_smooth = movmean(raw_data, 5); % 巴特沃斯低通滤波，截止频率根据采样率设定 [b, a] = butter(4, cutoff_freq / (fs / 2), 'low'); % 使用filtfilt做零相位滤波，避免相位偏移影响起跳点判定 data_filtered = filtfilt(b, a, data_smooth); end

这里有个关键点：千万不要用filter函数替代filtfilt。filter是有相位延迟的，处理反射波和透射波这种需要精确对时的信号，相位偏移会让起跳点的位置系统性地偏移几个采样点，最终算出的应力-应变曲线会有明显的偏差。filtfilt虽然计算量稍大一点，但零相位特性让波形中每个特征点的位置都保持准确，这在SHPB数据处理中是必须的。

2.2 时间轴对齐与零点修正

原始数据里，三组信号的时间轴虽然是同一块采集卡采出来的，理论上是天然对齐的，但实际操作中还是会遇到需要修正的地方。最常见的情况是实验前应变片电桥没有完全平衡，导致零漂现象——信号在入射波到来之前就已经偏离零点，有时甚至稳定在一个直流偏置电平上。

零点修正的常规做法是取触发信号到来之前那段“安静的”前基线数据，计算其均值作为基线值，然后用整条信号减去这个基线值。这段前基线一般取触发前200到500微秒的数据，既足够长能稳住统计均值，又不会太靠前以至于引入其他干扰。

时间轴对齐方面，还有一个容易忽略的操作是检查采样率设置。不同采集卡的标称采样率与实际采样率之间可能存在千分之几的误差，对于SHPB这种微秒量级的事件来说，这个误差积累下来会让入射波和透射波的对时产生几个微秒的偏差。严谨的做法是用标准方波信号对采集卡做一次时间标定，或者在每次实验后用示波器输出的脉冲信号校验一下时间轴。

2.3 入射波、反射波和透射波的特征点识别

特征点识别是程序里最核心、也最容易出问题的地方。入射波和反射波的起跳点对应应力波到达入射杆应变片位置的时刻，透射波的起跳点对应应力波穿过试件到达透射杆应变片位置的时刻。三个起跳点的时间关系承载着物理信息：入射波和反射波起跳点的时间差，就是应力波从入射杆应变片位置传到试件端面再反射回来的时间；入射波和透射波起跳点的时间差，反映应力波穿过试件的时间加上两段杆上应变片位置到试件端面的传播时间差。

自动识别起跳点的算法我用的是“阈值+斜率双重判定法”。核心思路是：在信号找到最小值（对于入射波取极大值的负向，因为压杆用压电或半导体应变片测压缩波通常输出为负向信号）后，往回搜索第一个斜率变化超过设定阈值的点，作为候选起跳点。具体步骤如下：

1. 找到信号的极值点（入射波取最小值，因为压缩波通常是负脉冲）。
2. 从极值点向前搜索，对信号求一阶差分，找到差分绝对值最大的点，作为最陡峭的上升/下降沿位置。
3. 从最陡峭位置继续向前，找到信号首次超过噪声带（即基线均值±n倍标准差）的点，定为起跳点。

这个逻辑在处理干净波形时相当可靠，但遇到波形前沿比较缓的情况会偏保守。所以我保留了手动修正的接口，用户可以在界面上直接点击调整起跳点的位置。程序会把自动识别的结果画出来，人工确认后进入下一环节。

2.4 三波法计算公式与参数标定

SHPB实验的理论基础是一维应力波理论，前提假设是杆中应力波无弥散、试件内部应力应变均匀。数据处理的核心公式是三波法：

• 应变率：(\dot{\varepsilon}(t) = \frac{2C_0}{L_s} \varepsilon_r(t))
• 应变：(\varepsilon(t) = \frac{2C_0}{L_s} \int_0^t \varepsilon_r(\tau) d\tau)
• 应力：(\sigma(t) = \frac{E A}{A_s} \varepsilon_t(t))

其中，C0是杆中的弹性波速，Ls是试件的初始长度，E是杆材的弹性模量，A是杆的横截面积，As是试件的横截面积，εr(t)是反射波应变信号，εt(t)是透射波应变信号。

这里需要特别注意的是，公式中用的εr和εt是杆表面应变片测得的应变，不是电压信号。所以程序里必须有一个“电压→应变”的换算环节，换算系数由应变片的灵敏系数、电桥接法和增益共同决定。很多初学者在这个环节出问题，用的参数不对，算出的应变和应力差了整整一个数量级。

应变换算公式为：(\varepsilon = \frac{2U_0}{K \cdot E_{bridge} \cdot G})，其中U0是采集到的电压，K是应变片灵敏系数（通常为2.0左右），E_bridge是桥压，G是放大器增益。如果是半桥接法，系数会不同，一定要根据实际电路选对公式。

杆的弹性模量和波速也是需要实际标定的参数，不能直接查表。标准做法是用已知长度的杆做一次撞击实验，测量应力波在杆中来回传播的时间，通过(C0 = 2L / \Delta t)算出实际波速。这样做的好处是把应变片粘贴位置误差、杆长误差等系统误差都包含在标定系数里，比单纯查手册要可靠得多。

3. 实操过程与核心环节实现

3.1 程序主界面与数据读取

我写的这套程序，主界面用MATLAB App Designer搭建，左边是文件列表和参数设置面板，右侧是波形显示区和结果输出区。运行程序后，第一步是选择实验数据文件夹，程序会自动读取文件夹内所有CSV格式的原始数据文件，并按照文件名中的编号进行排序。

数据读取模块的核心代码大致如下：

function data = read_shpb_data(filepath) opts = delimitedTextImportOptions('NumVariables', 4); opts.DataLines = [start_row, end_row]; opts.Delimiter = ','; opts.VariableNames = ['time', 'incident', 'reflected', 'transmitted']; data = readtable(filepath, opts); % 提取时间轴和信号，注意单位统一为秒和伏 t = data.time; incident_raw = data.incident; reflected_raw = data.reflected; transmitted_raw = data.transmitted; end

CSV文件的格式各个实验室可能不一样，有的导出四列数据，有的只导出三列，还有的带表头带说明文字。我建议在程序里做一个“列映射”的小功能，让用户手动指定哪一列是时间、哪一列是入射波、反射波和透射波，避免因为格式差异反复改代码。

读取完数据之后，程序会立刻画出原始波形图，让用户直观检查信号质量。这一步非常关键，如果原始信号明显饱和或者断线，就根本不值得继续处理，直接标记为无效数据跳过，能节省大量时间。

3.2 批量自动化处理与关键参数配置

批量处理是这套程序提升效率的核心能力。一次SHPB实验中，同一批材料往往要测3到5种应变率，每种应变率下做3到5根试件，算下来一次实验跑完少说也有15到25条数据。手工处理的话，半天时间就搭进去了，用程序批量跑，配制好参数后几分钟就能全部算完。

批量处理的关键参数包括：

• 杆材参数：弹性模量、密度、杆径、波速。
• 试件参数：初始长度、直径。
• 采集参数：采样率、桥压、增益、应变片灵敏系数。
• 算法参数：滤波截止频率、起跳点搜索阈值、信号截断位置。

这些参数全部放在一个配置文件中，程序启动时自动加载。配置文件用简单的键值对格式，方便直接文本编辑，也方便保存多套配置，比如“铝杆配置”“钢杆配置”“高温实验配置”，切实验条件时一键切换，不用重新输入。

批量处理时，每一条数据的处理结果会立即显示出来，包括计算出的峰值应力、峰值应变率、平均应变率和应力平衡系数。如果某个数值超出合理的物理范围，程序会在界面上用醒目的颜色标出，方便人工复核。

3.3 波形截取与积分算法实现

波形截取的目的是把反射波和透射波中有效的信号段分离出来，用于后续计算。截取位置的选取有一定的灵活性，但基本逻辑是固定的：

• 反射波的起点是入射波起跳点，终点是应力波在杆中来回反射多轮后信号衰减到接近基线处。
• 透射波的起点是透射波起跳点，终点是试件破坏或信号衰减到基线处。
• 为了避免末尾基线部分的随机波动影响积分结果，通常会适当保留一段尾部基线，在积分结束前再将基线扣除。

积分算法方面，我采用梯形法数值积分。MATLAB自带的trapz函数就可以完成这个工作，但直接对整个信号积分有个问题：如果尾部基线的均值不为零，积分结果会有一个线性漂移的误差积累。所以我在积分前先计算尾部基线的均值，然后从反射波信号中扣除这个值，再做积分，这样可以明显改善应变曲线的末端回零情况。

以下是计算应变、应力、应变率的核心代码：

function [strain_rate, strain, stress] = three_wave_method(c0, E, A_rod, L_s, A_s, reflected, transmitted, dt) % 应变率与反射波幅度成正比 strain_rate = -2 * c0 / L_s * reflected; % 应变是应变率对时间的积分 strain = cumtrapz(strain_rate) * dt; % 应力与透射波幅度成正比 stress = E * A_rod / A_s * transmitted; end

负号的处理需要特别说明一下：压缩波在应变片上的响应是负向信号，而力学中通常约定压应力和压应变为正，所以公式里要乘一个负号把符号扳回来。这个细节不处理好，画出来的应力-应变曲线整个倒过来，还会让人一头雾水。

3.4 动态平衡校核与数据有效性检验

做完三波法计算后，程序会同步做一个动态平衡校核。SHPB数据的有效性，很大程度上依赖于试件内部是否达到了应力平衡状态。简单来说，就是在试件两端分别测得的力（入射杆端的力由入射波加反射波得到，透射杆端的力由透射波得到），在试件加载的大部分时间内应该基本相等。

校验公式为：入射力(F_1(t) = E A (\varepsilon_i(t) + \varepsilon_r(t)))，透射力(F_2(t) = E A \varepsilon_t(t))。如果F1和F2两条曲线基本重合，说明试件受力均匀，数据有效；如果在加载上升段就分离很开，说明试件还没达到平衡就破坏了，这组数据需要谨慎使用或废弃。

我程序里会自动计算一个“应力平衡系数”：定义为整个加载段内F1和F2之差的绝对值与F1平均值的比值，比值小于5%就认为平衡良好，5%到10%之间需要手动判定，大于10%直接给出警告。这个检验逻辑做进程序后，我处理数据时心里踏实多了，不用再靠肉眼去对比两条曲线。

4. 常见问题与排查技巧实录

4.1 原始信号饱和与应变片粘贴问题

遇到的第一个高频问题是信号饱和。应变片测量范围有限，如果撞击速度太高，应变片输出超出采集卡的量程，波形就会出现削顶现象。削顶后的波形看起来平台很漂亮，但实际上是假的，算出来的应力偏低。

排查方法是看原始波形是否出现明显平顶。如果饱和发生了，只能重新做实验，降低撞击速度，或者减小应变片的灵敏系数。程序里我给原始波形绘制加了一个“削顶检测”功能，自动检查信号最大值是否接近采集卡满量程的一定比例，超过这个比例直接提醒用户。

另一个常见的物理层问题是应变片粘贴质量不佳，导致信号中出现大量高频振荡，尤其是透射波信号特别弱的情况。透射波弱通常意味着试件太软或者截面积太小，应力波大部分被反射回去了，这种情况的原始信号信噪比很低，靠滤波也救不回来。经验上是试件和杆的横截面积比值最好控制在0.6到1.0之间，太小的比值会让透射信号弱到无法使用。

4.2 起跳点误判的典型情况与修正策略

自动起跳点识别在两种情况下最容易出问题：

一种是波形前沿有“台阶”，也就是信号起跳后先缓慢上升一段，再突然加快。这种台阶通常是应变片粘贴不均匀或者撞击杆端面不平造成的。算法在选择最陡峭上升点时，可能选中后面那段，导致起跳点整体后移。修正策略是在算法里增加一个“局部极值检查”，在找到最陡峭点后，继续向前搜索，看是否存在一个明显斜率骤减的位置，如果有则把起跳点前移到那里。

另一种情况是信号噪声过大，导致算法在噪声带内就触发了阈值。针对这个问题，在阈值设定时采用自适应策略，先计算前基线段的统计特性，然后以基线均值加10倍标准差作为起跳阈值。10倍标准差这个经验值比固定阈值靠谱得多，能适应不同实验室的噪声底，不必每次手动改。

如果在自动识别后还需要手动微调，界面上提供了光标选点功能，直接点击目标位置即可更新起跳点。程序会自动记录手动修正的坐标值，便于事后追溯哪条数据被人为干预过。

4.3 计算结果的合理性校验与常见错误

程序跑完并不代表万事大吉，结果合理性校验是最后一道关卡。我一共设了三个快速的物理合理性检查：

第一，应变率曲线的形状是否符合预期。对于均质材料，应变率曲线通常在一个平台附近波动，如果应变率峰值是平均值的两三倍甚至更高，说明试件发生了明显的非均匀变形，这组数据的代表性存疑。

第二，应力-应变曲线的弹性段斜率是否合理。虽然SHPB是动态实验，但加载初始阶段材料仍处于线弹性状态，应力-应变曲线的斜率应当接近材料的弹性模量，或者至少是同一个数量级。如果差了两三个数量级，基本可以肯定是换算参数出了问题。

第三，试件最终应变是否与实物吻合。实验做完后通常可以观察到试件的变形甚至破碎，估算的最终应变应该和肉眼观察的结果相互印证。如果程序算出的应变为20%，试件却几乎没变形，那就说明积分环节存在问题，比如截取段选得太短或者基线扣除没做好。

常见错误汇总成表格如下：

4.4 数据导出与绘图的细节优化

最后一步是结果输出。程序自动绘制的图包括：原始波形图、滤波后波形图、应力平衡检查图和最终的应力-应变曲线图。每种图都按照论文发表的基本要求设置了字体大小、线条粗细和图例位置，基本做到“一键出图”，不用再去Origin里二次修整。

数据导出方面，程序会把每条试件的关键参数（峰值应力、峰值应变率、平均应变率、最大应变、应力平衡系数）汇总成一个表格导出到Excel，同时把完整的三波法计算结果（时间、应变率、应力、应变）保存成CSV文件，方便后续用其他软件做深入分析。

绘图和导出这两个环节，我的建议是格式统一化。实验室里大家风格迥异，有的喜欢红色曲线，有的要黑色实线，有的图例放右下角有的放左上角。把出图样式统一在程序设置里，能避免大量无谓的返工。

5. 几个实操心得与后续扩展方向

用这套程序跑了不下几百组数据之后，有几个体会特别深。第一个是：数据处理程序的价值，不只是省时间，更在于它强制你梳理清楚每一个处理步骤的物理含义。以前手工处理的时候，很多环节是“凭感觉来”，但写程序逼着你把每个参数、每个判断准则都明确下来。比如“起跳点在哪里”，以前鼠标点一下就行，写程序时你就得想清楚：什么样的点算起跳点？判断依据是什么？不同的判断方式对结果影响有多大？这个过程本身对实验理解的提升，比程序跑出来的结果更值钱。

第二个体会是：程序里一定要保留人工干预的接口，但尽量让干预变得“昂贵”。我的做法是，手动修正的每条数据，界面上都会显示一个红色的“M”标记，表示这条数据的某些参数被人工调整过。这能提醒自己，也对数据的可追溯性有好处。如果后续发现某条数据的结果明显异常，可以先检查它是不是被人工干预过，排除操作失误的因素。

第三个体会是关于批量处理的，我个人建议不要完全盲目地批量跑。最稳妥的做法是先手动处理一条参考数据，确认所有参数设置正确、结果合理，再启动批量模式。我见过不少同事图省事，拿到数据直接批量跑，结果某一步参数设置错了，一整批数据全部报废，又得重新处理。先验证再批量化，这个习惯真得养成。

后续的扩展方向上，我目前正在做两个改进。一个是把Python深度学习的方法引入起跳点自动识别，看看能不能用更鲁棒的方式处理噪声很大或者波形形状怪异的数据。另一个是增加试件的高温/低温环境参数修正模块，因为温度对杆材波速和试件本构都有影响，目前SHPB的高温数据处理大多还要靠人工修正，程序化难度比较大，但很值得做。如果你经常做SHPB实验，建议也把这套流程掌握起来，不一定非得完全用我的方案，但数据处理逻辑和里面的坑，应该能帮你少走不少弯路。

本文还有配套的精品资源，点击获取