简介在MATLAB环境下实现的雨流计数法工具包适用于机械工程、材料科学和疲劳分析场景帮助用户从复杂应力应变历史中提取关键循环进而评估结构疲劳寿命。压缩包共15个文件大小仅112KB以.m脚本为主包含可直接调用的雨流计数函数、动态链接库、C语言源文件及HTML说明文档并提供多个演示示例与结果绘图工具便于快速理解算法流程和进行二次开发。该程序通过查找极值、分割上升/下降段并配对循环最终生成雨流矩阵为疲劳损伤计算提供数据基础。已有1423人学习下载适合需要处理实测载荷谱、开展疲劳寿命预测的MATLAB使用者。下载后可对照示例学习完整实现思路也可将程序集成到自己的数据流程中直接获取循环统计结果。1. 为什么疲劳分析把雨流计数当成必修课做疲劳分析的人手里最麻烦的不是有限元算出的应力云图而是那一串动辄几万行的应力时间历程。直接按峰值数循环小幅波动也会被当成损伤算出来的寿命保守到没法用只找最大最小值又把嵌套的小循环全部吞掉结果偏危险。雨流计数法Rainflow Counting通过把任意变幅载荷拆成一组组闭合的应力-应变滞后环同时保留大循环的幅值、均值和小循环的局部损伤贡献再配合Miner线性累积损伤规则估算寿命。这套方法在材料科学、机械工程、风电和汽车疲劳评估里几乎是标准前置步骤。rainflow.zip 里的 MATLAB 雨流计数程序正好覆盖了从原始载荷到雨流矩阵的完整链路尤其适合写论文、做课程设计和工程评估的人拿来直接改。2. 雨流计数原理与rainflow.zip的文件分工2.1 从应力-应变滞后环到“雨滴流过屋檐”的配对规则雨流计数的底层逻辑来自材料在循环载荷下的应力-应变滞后回线。每经历一个完整循环材料内部就留下一个闭合的滞后环这个环的顶点和底点对应的应力幅与均值是计算疲劳损伤的核心参数。真实载荷不会按整齐的正弦波走而是一条上下反复、大小嵌套的曲线所以算法需要先把曲线转成峰谷序列再按“雨流”规则提取循环。标准四峰谷规则是按时间顺序取四个峰谷点如果中间段始终小于等于两端点或者中间段始终大于等于两端点就认为这一段构成了一个完整循环可以直接抠出来。抠出来的循环记录其幅值、均值和次数剩下的不完整段会在下一轮迭代中继续配对。rainflow.zip 里的核心算法就是按照这个规则实现的。% 用一个简单序列观察峰谷提取效果 sig [0; 3; 1; 4; -2; 5; 0]; ext sig2ext(sig); disp(原始序列: ); disp(sig); disp(峰谷序列: ); disp(ext);这里sig2ext.m负责把sig转换成语义上的峰谷序列只保留转折点。原始序列中[3;1;4]这种“中间点既不是最高也不是最低”的转折就是雨流算法要处理的嵌套结构。逻辑说明sig2ext输出给rainflow之后算法才能正确识别出1到3到4这样的小循环而不会把它误判成一个大上升段。参数方面如果数据里有大量抖动应该先滤波否则峰谷提取会产生大量伪转折点。2.2 rainflow.m、rainflow.c、sig2ext.m 各自负责什么zip 包内的文件并不是单一函数而是一整套前后处理工具。我把主要文件的分工整理成一张表方便对照实际代码阅读。文件名作用调用时机sig2ext.m从原始时间序列提取峰谷转点删除等值和中间点在调用rainflow之前运行rainflow.m主计数函数接收峰谷序列并返回雨流矩阵核心入口rainflow.c雨流算法的 C 语言实现用于编译成动态库加速需要高性能或批量处理时使用rainflow.dll编译好的动态库老版本 MATLAB 中可直接被调用作为rainflow.m的底层实现rfhist.m绘制雨流矩阵的三维直方图结果可视化rfmatrix.m生成雨流矩阵的平面图或矩阵数据后处理rfdemo1.m/rfdemo2.m/rfdemo3.m三个不同场景的演示脚本快速上手时直接运行实际使用时流程通常是先对原始应力应变序列调用sig2ext得到峰谷再将其传入rainflow得到矩阵最后用rfhist画图。rainflow.dll是较早的 C 语言编译产物在 64 位 MATLAB 上一般需要重新编译rainflow.c这一点在最后一章会详细讲。2.3 跑通rfdemo1.m从随机载荷到雨流直方图rfdemo1.m 是整个包中最容易入门的脚本。先看它的整体结构再单独执行% rfdemo1.m 核心流程示例 % 生成一段带有噪声的随机载荷 t (0:0.01:20); stress 80 30*sin(0.5*t) 8*sin(2.3*t) 3*randn(size(t)); % 第一步提取峰谷 extr sig2ext(stress); % 第二步雨流计数 rf rainflow(extr); % 第三步绘制雨流直方图 rfhist(rf);执行后能在一个三维图上看到横轴为应力幅、纵轴为均值、竖轴为循环次数的柱状分布。说明一下stress是模拟应变片实测数据包含多个叠加简谐波和随机噪声sig2ext把几十万点压缩成几个转折点rainflow返回的rf矩阵第一列是循环幅值第二列是循环均值第三列是每个区间内的循环次数。运行前需要确认当前目录在解压后的文件夹内且路径里没有中文。老版本脚本里出现过.asv这种自动保存文件不影响运行可以直接忽略。3. 手把手处理自己的载荷谱输入输出、阈值与结果导出3.1 数据预处理连续时间序列如何变成干净的峰谷序列从信号采集器里拿到的数据是等间隔采样点可能包含重复值、平台段和高频噪声。雨流算法对这类数据非常敏感连续重复值会让内部比较逻辑产生误判噪声则会产生大量无关紧要的小循环。所以我的习惯是先做峰谷提取再设置幅值阈值过滤最后才交给雨流计数。下面是一个简化版的sig2ext实现逻辑和包内的思路一致function ext sig2ext(sig) % 确保输入是列向量 if size(sig, 1) 1 sig sig; end n length(sig); keep true(n, 1); % 删除单调段中间点只保留转折点 for i 2 : n - 1 d1 sig(i) - sig(i - 1); d2 sig(i 1) - sig(i); if (d1 * d2 0) || (d1 0) keep(i) false; end end ext sig(keep); end这段代码的关键在于d1 * d2 0当左右两点同向或有一个为零时当前点就不是转折点删掉后能大幅压缩数据。压缩后原始波形中一段缓慢上升的直线变成了首尾两个点但雨流算法只关心幅值和均值不关心中间路径所以这样处理是安全的。实际包里的sig2ext.m可能还做了首尾保留和空值处理原理一致。参数上如果数据量特别大可以先用resample降到合适采样率。3.2 rainflow函数调用约定与常见格式坑rainflow.m 的调用格式在不同版本里差异不大但返回参数的顺序有时会让人踩坑。我拿到这个包的第一件事就是用help rainflow查看说明再运行一次 demo确认返回的是单个雨流矩阵还是多个分开的数组。% 调用方式一只取雨流矩阵 rfm rainflow(extr); % 调用方式二返回幅值、均值、计数数组 [amp, mean_val, counts] rainflow(extr); % 如果用方式二需要转成矩阵再画图 rfm2 [amp, mean_val, counts]; rfhist(rfm2);这里amp和mean_val是每个循环的幅值与均值counts是循环次数。很多老代码直接把三个输出拼成矩阵但新版本可能直接返回矩阵所以需要判断判断。逻辑说明方式二适合后续做概率统计比如按平均应力修正方式一适合直接丢给rfhist画图。如果提示输入数据包含NaN先执行extr extr(~isnan(extr));清理。常见的输入错误和解决办法见下表现象可能原因处理提示维度不匹配输入是行向量而函数需要列向量统一用extr extr(:)雨流矩阵全是 0载荷本身没有完整循环或峰谷提取失败先单独运行sig2ext检查输出长度运行时间过长数据点达到百万级M 脚本迭代慢用编译后的rainflow.c版本三维图空白循环数量太少或数据范围太大调整rfhist的坐标轴范围还有一个隐藏坑雨流计数前不要对数据做“平滑处理”过度因为过度平滑会滤掉真实的小循环使损伤评估偏小。正确的做法是先保留原始信号做一次计数再对比滤波后的结果两者差异可以反映噪声对寿命评估的影响。3.3 从rfhist.m到自定义导出雨流矩阵怎么读rfhist.m 把雨流矩阵映射成三维直方图每个柱子的位置代表箱体中心高度代表循环数。读取这个图时重点关注两个区域高幅值低均值区域通常对应主导损伤低幅值高均值区域对应高平均应力下的轻微循环。不过三维图只能用来定性观察定量计算还是要读取矩阵本身。% 导出雨流计数结果到 CSV rfm rainflow(extr); if size(rfm, 2) 3 T table(rfm(:,1), rfm(:,2), rfm(:,3), ... VariableNames, {Amplitude, MeanStress, Cycles}); writetable(T, rainflow_result.csv); else T table(rfm(:,1), rfm(:,2), ... VariableNames, {Amplitude, Cycles}); writetable(T, rainflow_result.csv); end这段代码先判断返回类型再导出为标准的三列表格。为什么要把均值也导出因为后面的寿命估算需要用 Goodman 或 Gerber 公式修正平均应力。很多初学者只保留幅值就把均值丢了导致后续无法做修正只能做最粗糙的等幅寿命估算。参数说明Amplitude是循环应力半幅MeanStress是均值Cycles是循环次数单位与输入数据一致。CSV 文件可以继续用 Python 的 pandas 读取做进一步的概率统计分析。4. 进阶把雨流矩阵变成疲劳寿命估算和试验载荷谱4.1 Miner线性累积损伤从循环数到寿命的数学桥雨流矩阵本身只是载荷统计结果真正得到寿命还需要结合材料 S-N 曲线。最常见的是 Basquin 公式配合 Miner 线性损伤% 材料 S-N 曲线参数 Sf 300; % 疲劳强度系数单位与载荷数据一致 k 8; % Basquin 指数常见金属材料约 5~10 % 从雨流结果取幅值和循环数 amp rfm(:, 1); mean_s rfm(:, 2); cycles rfm(:, 3); % 平均应力修正Goodman 形式 Su 600; % 材料极限强度 amp_adj amp ./ (1 - mean_s / Su); % 每个循环对应的寿命 N_cycle (Sf ./ amp_adj) .^ k; % Miner 损伤累积 D sum(cycles ./ N_cycle); % 以当前载荷块作为单位 1 的寿命 life_blocks 1 / D; fprintf(累积损伤 D %.4f\n预测寿命 %.2f 个载荷块\n, D, life_blocks);代码里amp_adj是用 Goodman 公式把非零均值循环折算成等效零均值幅值。这个修正非常关键相同幅值下拉压均值的影响可能差出好几倍寿命。D是单位载荷块造成的总损伤小于 1 表示当前载荷块重复life_blocks次后会发生破坏。注意这里没有考虑小循环的疲劳极限截断如果amp_adj低于疲劳极限通常把对应的N_cycle设为无穷大也就是不贡献损伤。实现时直接把这些行的1 ./ N_cycle改为 0 即可。4.2 rfdemo2/3演示的残差循环与载荷重排列逻辑雨流计数会留下一个或几个不闭合的残差段这些段实际上是整个载荷历程里最大的那个主循环。如果直接把残差丢掉损伤会被低估所以 rfdemo2.m 和 rfdemo3.m 里通常演示了将残差单独提取并加入最后结果的方法。残差处理的核心是把峰值和谷值中绝对值最大的两点找出来作为主循环的上下端点再与前面抠出来的循环叠加。我一般会用下面的步骤来处理残差先调用雨流计数得到完整循环的下标然后剩下未匹配的峰谷点按原顺序组成一个新的序列接着把这个序列的首尾连接成一个闭合循环但要注意不能重复计算已经计过的点。工程上更常用的做法是把残差当作一个半循环计入损伤时按半循环折算。实际拆包时rfdemo2.m 里应该能看到类似的逻辑对同一个载荷序列先正序计数一遍再倒序计数一遍最后把两个结果合并避免方向性偏差。% 残差处理示意正序和倒序各计一次 rfm_forward rainflow(extr); rfm_backward rainflow(flipud(extr)); % 把倒序结果的幅值均值保持原方向不变 rfm_backward(:, 1) rfm_backward(:, 1); rfm_backward(:, 2) rfm_backward(:, 2); % 合并计数 rfm_combined [rfm_forward; rfm_backward]; % 按箱体合并相同幅值和均值区间的循环次数这段代码展示了为什么需要双向计数某些非对称载荷历程正序和倒序会留下不同残差取两者的总和会更接近真实的损伤累积。参数flipud只是反转时间方向不改变幅值物理意义。需要提醒的是双向计数结果不能直接相加后马上画图要先按箱体分组聚合否则同一个循环会被计两次。常见做法是用accumarray或uniquetol将非常接近的箱体合并。4.3 批量处理多个应力文件的自动化脚本实际项目中不可能只处理一条载荷谱。风电、汽车路谱数据往往来自几十个通道或几十个工况所以我习惯把雨流计数封装成一个函数再用循环批处理。下面的脚本处理目录下所有 CSV 文件% 批量雨流计数脚本 files dir(load_case_*.csv); results cell(length(files), 1); for i 1:length(files) data readmatrix(files(i).name); sig data(:, 2); % 假设第二列是应力 ext sig2ext(sig); rfm rainflow(ext); results{i} rfm; fprintf(已处理 %s循环数 %d\n, files(i).name, size(rfm, 1)); end % 将多个结果累加为累计雨流矩阵 total_rfm cell2mat(results);批量处理时最容易遇到的是数据量过大导致内存不足。一个 50 万点的信号sig2ext之后变成几千个峰谷rainflow很快就能算完但如果直接把 50 万点原始数据丢进去点对点比较会非常慢。所以这条链路的顺序非常重要先压缩、再计数。代码里的cell2mat拼接后的矩阵可能很大建议再用histcounts2统计成箱体矩阵方便后续损伤计算。5. 排错与边界rainflow.dll、阈值和验证技巧5.1 64位MATLAB下重新编译rainflow.c压缩包里的rainflow.dll是很早以前用 32 位 MATLAB 编译的在 64 位 MATLAB 上直接调用会报“无效的 MEX 文件”错误。常见的解决办法是在 MATLAB 里用mex重新编译mex -setup mex rainflow.c编译成功后会生成rainflow.mexw64或类似文件运行时 MATLAB 会优先调用新生成的 mex 文件。如果编译时报缺少头文件通常是缺少 MEX 配置安装配套的 C 编译器就好。5.2 用正弦波快速验证结果是否正确为了确认算法理解没出错我每次拿到新雨流程序都会先用一个完整正弦周期做验证sig [0; 1; 0; -1; 0]; ext sig2ext(sig); rfm rainflow(ext); disp(rfm);一个完整正弦周期从 0 到 1 再到 0 到 -1 再到 0理论上应该输出 1 个循环幅值 1均值 0。如果多出来循环说明sig2ext把首尾连接部分误判成了另一个循环。这个验证方法尤其适合检查二手代码能快速定位是算法问题还是数据预处理问题。5.3 阈值过滤与等值点的边界雨流计数里最隐蔽的坑是平台段。设备停机、传感器饱和都会产生连续相同数值的点这些点在峰谷提取时如果不特殊处理会产生无数个幅值为 0 的循环。处理方式是在sig2ext内部判断相邻差值是否严格不为零并把平台段压缩为单个点。阈值过滤则放在计数之前用最大幅值的百分比过滤小幅波动才能避免噪声主导损伤结果。最后提醒一句任何雨流程序都只能忠实反映你输入的数据决定寿命评估质量的往往是你对信号的预处理和对残差的理解。本文还有配套的精品资源点击获取