简介在机械工程与材料疲劳分析中雨流计数法用于将复杂的应力应变历程转化为可统计的完整循环。这份基于MATLAB实现的资源包面向需要处理载荷谱的工程师、科研人员和相关专业学生解决了从数据预处理、极值识别、循环匹配到雨流矩阵生成的全流程需求。压缩包共包含15个文件其中7个m脚本与函数覆盖算法主体与示例演示辅以C源码、dll动态库、HTML说明文档及可视化图像整体仅112KB精简且便于理解和二次开发。目前已有1423人学习下载是疲劳寿命预测场景中较为实用的参考工具。借助计数函数、特征提取脚本和多个demo使用者可直接导入自己的应力应变数据快速获得等效应力循环与雨流矩阵从而支撑结构疲劳损伤评估与寿命预测。1. 用Matlab跑雨流计数为什么说“关键在参数和残谱”打开rainflow.zip里面通常只有一个rainflow.m和几个辅助文件很多人第一反应就是把载荷数组塞进去然后对着一串输出矩阵发懵。雨流计数Rainflow counting在疲劳载荷谱分析里的地位相当于傅里叶变换在信号处理里面它把一段不规则的时间历程拆成完整的应力循环每个循环用幅值和均值表征再映射到S-N曲线做Miner损伤累积。但实际跑通后你会发现结果行数比肉眼数的载波周期多或者全循环和半循环混在一起不知道该怎么用。这往往不是算法代码坏了而是你没有先校准两件事输入数据的峰谷抽取规则以及残谱到底按0.5还是按1处理。这篇文章给出一个完整、可直接运行的Matlab实现从峰谷抽取、四点法主循环到参数设置和损伤计算全部对齐读完你也能判断网上下载的rainflow.zip是否靠谱。2. 雨流计数原理与Matlab函数实现从峰谷提取到四点法出循环2.1 雨流计数法在疲劳分析中的位置疲劳寿命预估最怕的不是平均应力高而是随机载荷里藏着各种大小不等的闭环。直接数峰谷或者用穿越计数会把内部的非完整循环漏掉或者重复计入。雨流计数把载荷历程当作屋檐雨滴从极值点流下形成一个个完整的滞后环。ASTM E1049-85(2011)是几乎所有雨流实现的参考基准它规定输出至少包含循环的幅值Sa、均值Sm和循环次数。对于全循环计数为1对于残谱中无法闭合的半循环计数为0.5。后面做Miner累积时只需要这三列数据所以实现雨流计数的核心其实就是把峰谷序列按“包含关系”反复裁剪直到剩下不能配对的残余点。2.2 先做峰谷抽取过滤掉无效斜率点很多原始载荷信号是直线段连接起来的直接拿相邻点去做四点法会留下大量无意义的中间点。必须先把信号压缩成峰谷交替的极值序列。常见做法是用diff做一阶差分再对差分结果求二阶差分斜率变号的位置就是极值点。function pv extractPeakValley(signal) % extractPeakValley 从原始载荷时间历程中提取峰谷序列 % signal: 一维列向量载荷时间历程 % pv : 只保留极值点和两端的序列相邻点符号必然交替 signal signal(:); delta diff(signal); % 相邻点差值 if all(delta 0) pv signal([1 end]); return; end signChange find(diff(delta) ~ 0) 1; % 斜率变号的位置 idx unique([1; signChange; length(signal)]); pv signal(idx); end这段代码里diff(delta) ~ 0表示当前点的前一段差值和后一段差值不同号也就是局部极值。unique用来去重防止平台段被重复选中。不过要特别注意如果原始信号存在宽度大于1的平顶diff会出现连续的0这时候diff(delta) ~ 0会漏掉平顶终点。工业数据里这种平台很常见建议在调用前先对连续重复值做一次去重也就是只保留每个平台的首尾两点。2.3 四点法提取循环Matlab核心算法实现峰谷序列准备好之后主循环就可以开始了。这里采用A. Nieslony经典四点法思路用一个缓冲数组R模拟雨流过程每读入一个新峰谷点就检查R末尾四个点是否满足“中间两个点形成的行程被前后两段行程包住”如果满足则提取一个完整循环然后删除中间两个点让修改后的序列继续参与检查。function [cycles, residual] rainflowCounting(pv) % rainflowCounting 四点法雨流计数 % pv : extractPeakValley() 输出的峰谷序列 % cycles : Mx3矩阵依次为 [幅值 Sa, 均值 Sm, 循环次数] % residual : 未配对成完整循环的剩余峰谷点 pv pv(:); n length(pv); if n 4 cycles []; residual pv; return; end R pv(1:2); % 缓冲区先放前两个点 cycles zeros(floor(n/2) 64, 3); cIdx 0; i 3; while i n R(end1) pv(i); % 新点送入缓冲 while length(R) 4 y1 R(end-3); y2 R(end-2); y3 R(end-1); y4 R(end); d1 y2 - y1; d2 y3 - y2; d3 y4 - y3; if d1*d2 0 d2*d3 0 ... abs(d2) abs(d1) abs(d2) abs(d3) cIdx cIdx 1; cycles(cIdx,:) [abs(d2)/2, (y2y3)/2, 1]; R(end-2:end-1) []; % 删除y2,y3让剩余序列继续比较 else break; end end i i 1; end % 残余点按相邻点对拆成半循环 for k 1:length(R)-1 cIdx cIdx 1; cycles(cIdx,:) [abs(R(k1)-R(k))/2, (R(k)R(k1))/2, 0.5]; end cycles(cIdx1:end,:) []; residual R; end代码里d1*d20用来保证y2是极值点d2*d30保证y3也是极值点两者方向相反。abs(d2) abs(d1)和abs(d2) abs(d3)表示y2到y3这段行程不会超过相邻的前段和后段因此在雨流轨迹上它先被“切”出来。循环提取后删除中间两个点相当于这段已经闭合剩下的重新连接继续检查。最后缓冲区里剩下的点无法再形成闭合循环按半循环处理计数为0.5。这里输出的第一列是幅值abs(d2)/2不是峰谷差注意别和某些zip包里的“range”混淆。2.4 输出结果格式与ASTM E1049的对应关系无论是自己写的还是从rainflow.zip里拿到的代码第一步都应该确认输出矩阵的列顺序。以下是我建议的统一约定输出列字段说明第1列幅值Sa(峰-谷)/2即循环的半程第2列均值Sm(峰谷)/2对应平均应力第3列循环计数全循环1半循环0.5ASTM E1049并不强制规定列顺序所以很多zip包会输出[range, mean, count]甚至[start_index, end_index, range, mean]。拿到外部代码后先用一个已知波形跑一遍逐列对比再做后续操作否则Miner损伤算出来会完全对不上。3. 在Matlab中设置雨流计数参数数据预处理、门槛值、残谱处理3.1 输入数据的常见转换从Excel/CSV导入和时间戳处理实际工作时载荷数据往往是csv或xlsx文件第一列时间戳第二列应变或应力。读取后第一步是清洗去掉NaN和重复时间戳。Matlab的readmatrix从R2019a开始可用之前的版本要退回csvread或xlsread。data readmatrix(load_history.csv); % 或者 .xlsx strain data(:,2); % 假设第二列是应力或应变 strain strain(~isnan(strain)); % 删除无效行 % 如果采样率太高可以按时间抽稀 t data(:,1); keepIdx 1:round(length(t)/1e6):length(t); % 保留最多100万点 strain strain(keepIdx);keepIdx的步长根据总点数调整工程上100万点做雨流完全够用。抽稀前最好先对信号做峰值保持或者低通滤波否则会丢掉真实的极值点。直接每隔N点抽取虽然简单但可能切掉局部峰导致后续循环计数偏小。3.2 幅值门槛和均值离散化避免小循环淹没统计载荷信号里总是混着高频噪声噪声产生的微小循环会让cycles矩阵膨胀几十倍而对疲劳损伤几乎没有贡献。设置一个幅值门槛是最直接的过滤方式。门槛值可以用整个信号最大最小值的百分比也可以用绝对应力值。参数推荐设置作用ampThreshold(max-min)*0.01 ~ 0.05去掉噪声小循环meanBinSize分成20~50个区间压缩后续统计矩阵大小residualModehalf 或 closed决定残谱半循环是否闭合过滤代码很简单amp cycles(:,1); meanVal cycles(:,2); cnt cycles(:,3); valid amp ampThreshold; cyclesFiltered cycles(valid, :);要注意的是门槛设太高会把真实的小幅值载荷循环也去掉低估损伤设太低则过滤效果差。最稳妥的做法是先绘制幅值分布直方图看0.5倍噪声以下的峰在哪里把ampThreshold放在噪声峰和真实循环峰之间的谷位上。3.3 残谱的两种处理方式半循环 vs 闭合循环雨流计数结束后缓冲数组里剩下的点无法配对成完整循环这就是残谱。标准做法是把残谱的相邻点对当作半循环计数为0.5。但这种处理在Miner损伤中会带来一个尴尬半循环的两端没有闭合实际损伤可能被低估。另一个常见做法是把残谱首尾相连形成闭合循环这样所有残余贡献都被归入一个完整加载路径。如果要把残谱当成一个闭合循环可以在rainflowCounting返回后手动拼接if ~isempty(residual) length(residual) 2 delta residual(end) - residual(1); closedCycle [abs(delta)/2, (residual(1)residual(end))/2, 1]; cyclesFiltered(end1, :) closedCycle; end大多数工程规范比如风力发电机组载荷计算会把残谱按半循环处理因为首尾闭合会引入一个实际不存在的大循环所以除非你的规范明确要求闭合否则保留半循环计数0.5更安全。3.4 一个完整的demo脚本模板把前面的代码串起来就是一个可以直接运行的测试脚本% demo_rainflow.m t (0:0.01:50); signal 10*sin(2*pi*0.15*t) 4*sin(2*pi*0.61*t 1.2); pv extractPeakValley(signal); [cycles, residual] rainflowCounting(pv); fprintf(总循环数: %d\n, size(cycles,1)); fprintf(最大幅值: %.3f\n, max(cycles(:,1))); histogram(cycles(:,1), 30); xlabel(幅值); ylabel(循环数);脚本先构造一个双正弦叠加信号提取峰谷后做雨流计数最终打印循环总数和最大幅值。如果输出循环数是0说明extractPeakValley返回的序列长度小于4检查原始数据是否全是常量或者清洗过度。4. 实际案例用雨流计数法处理正弦波叠加随机谱并验证计数结果4.1 构造测试信号为了验证函数的正确性先用一个可控信号慢波加快波慢波周期10秒快波周期2.5秒叠加后肉眼能数出大约5个慢波周期但快波会产生附加循环。构造信号代码如下t (0:0.001:50); signal 10*sin(0.2*pi*t) 3*sin(0.8*pi*t 0.5) 0.5*randn(size(t));加入高斯噪声是为了测试ampThreshold的实际效果。这个信号总点数5万直接雨流计数会得到大量噪声循环正好用来观察阈值过滤前后的差异。4.2 运行计数并绘制雨流参数直方图调用我们的函数后用二维直方图看均值-幅值分布pv extractPeakValley(signal); [cycles, residual] rainflowCounting(pv); figure; histogram2(cycles(:,2), cycles(:,1), 30, DisplayStyle,tile); xlabel(均值 Sm); ylabel(幅值 Sa); colorbar;histogram2的第一个输入是均值第二个是幅值每个tile的颜色代表该应力区间内循环数量。从图里能直观看到慢波产生的大幅值循环集中在均值0附近噪声循环则散落在小幅值区域。此时可以先设置ampThreshold 0.5再重新统计过滤后的循环数通常过滤后循环总数会下降一到两个数量级。4.3 验证计数结果循环数和幅值总和为什么对不上很多人在这一步会发问明明信号只有5个慢波周期为什么输出循环数大于5答案在于快波循环和噪声循环也都被计数了另外残谱半循环也会增加计数。为了单独验证算法用一段无噪声正弦波测试testT (0:0.01:20); testSig 5*sin(2*pi*0.1*testT); % 2个完整周期 pvTest extractPeakValley(testSig); [cycTest, ~] rainflowCounting(pvTest); expected 2; % 两个完整正弦周期产生两个全循环 actual sum(cycTest(:,3) 1);说明如下测试信号理论循环数雨流输出2周期正弦波2个全循环2个全循环少量残谱半循环5周期正弦波快波5若干快波循环大于5但快波循环可辨识慢波和快波会形成不同幅值的循环落在直方图的不同区间。如果输出里出现半循环计数0.5这属于残谱不是错误。反过来如果全循环数量大于理论值通常是因为峰值谷值抽取时把平台段重复计数检查extractPeakValley对平顶信号的处理即可。4.4 不同实现之间的差异rainflow.zip中的旧版和新版不一致怎么办网上流传的rainflow.zip版本很多有的输出[range, mean, count]有的输出[mean, range, start, end]还有的是C-MEX编译版本。处理方式不是一个个试而是先构造一个非常简单的闭环序列做基准测试比如[0 1 0 -1 0]。该序列包含一个从-1到1再到-1的闭合滞后环理论上应该输出1个全循环幅值1均值0。用这段基准去跑你下载的代码逐列比对哪个输出列与这些值吻合就能快速定位列顺序和单位。baseline [0 1 0 -1 0]; pvBase extractPeakValley(baseline); [cycBase, resBase] rainflowCounting(pvBase); disp(cycBase);这段脚本输出应该包含一行[1, 0, 1]幅值1均值0计数1。如果实际输出是[2, 0, 1]说明单位是峰谷差而不是幅值如果均值列是0.5说明列顺序不同。用这样10个点就能排查完外部代码再放到真实数据上。5. 进阶应用把雨流计数结果直接做Miner疲劳损伤累积5.1 从雨流结果到S-N曲线Goodman平均应力修正雨流计数输出的均值对应循环的平均应力。同一个幅值在不同均值下对损伤的影响不同通常用Goodman直线修正到等效应力幅。设材料极限强度为Su等效应力幅Sa_eq Sa / (1 - Sm / Su)。Sa cycles(:,1); Sm cycles(:,2); Su 800; % 材料抗拉强度MPa Sa_eq Sa ./ (1 - Sm ./ Su); Sa_eq(Sm Su) inf; % 超过极限的直接判为失效注意Sm有可能是负的也就是压缩平均应力此时Goodman修正后的Sa_eq会变小这是合理的因为压平均应力通常提高疲劳寿命。如果分母出现0或负数说明该循环已静载破坏直接标记为inf或从统计中剔除。5.2 用循环结果矩阵快速计算Miner损伤得到等效应力幅后用Basquin公式估算每个循环的寿命N_f C * Sa_eq^(-m)其中C和m是材料S-N曲线参数。Miner线性累积损伤为所有循环损伤的总和C 1e12; % 材料常数 m 3.0; % S-N曲线斜率 N_f C * Sa_eq.^(-m); D_total sum(cycles(:,3) ./ N_f); % 半循环的0.5也会计入 if D_total 1 fprintf(预测疲劳失效累积损伤 D %.3f\n, D_total); else fprintf(安全累积损伤 D %.4f\n, D_total); end这里cycles(:,3)同时包含全循环和半循环半循环的0.5会直接参与求和不需要额外标记。计算过程中如果Sa_eq中有0N_f会变成无穷损伤为0对应零幅值循环可以提前过滤。整个流程从雨流计数到损伤计算就闭环了。5.3 性能优化与大规模数据的处理技巧长周期载荷信号动辄几百万点提取峰谷后点数通常能压缩到原来的1/10到1/100因此rainflowCounting输入规模并不大。但输出循环矩阵可能仍然很大尤其是噪声未滤除时。优化思路有三个第一优先做峰谷抽取和幅值门槛过滤把循环数降到十万以内第二把幅值和均值离散化到二维直方图后面损伤计算直接对每个bin累计而不是逐行循环第三如果单个文件实在太长把信号分段做雨流分段之间会丢失跨界循环所以拆段时要在首尾各保留一段重叠区域重叠长度至少覆盖一个最大周期。numBins 50; binEdgesAmp linspace(0, max(cycles(:,1)), numBins); binEdgesMean linspace(min(cycles(:,2)), max(cycles(:,2)), numBins); damageMap histcounts2(cycles(:,2), cycles(:,1), ... binEdgesMean, binEdgesAmp);这段代码把循环按均值和幅值投到50×50的bin矩阵中后续Miner损伤计算就可以对每个bin里的循环数统一处理显著加速。如果下载的rainflow.zip里带C-MEX编译版本记得先运行mex -setup配置编译器再编译对应文件否则Matlab会静默回退到较慢的m文件版本。本文还有配套的精品资源点击获取