简介在航空航天、桥梁建筑、机械装备等领域结构常承受交变随机载荷疲劳失效是其主要破坏形式之一而雨流计数法正是提取疲劳载荷谱的关键技术。这套Matlab代码包实现了三点法和四点法两种计数算法文件内包含主程序cycle_counting_3.m、cycle_counting_4.m以及辅助函数fun.m三点法用三个连续数据点判定一个循环四点法通过增加峰谷检查识别更完整的载荷循环两者互补供用户对照学习。压缩包共12个文件以m源码为核心配以4个Excel载荷数据文件load_F*.xls作为验证输入另有docx格式的计算方法说明与数据校验文档以及txt版使用说明整体仅261KB轻量易用。目前已有149人学习适合材料疲劳分析、结构应力应变领域的学生和工程师可通过阅读文档、运行两个算法文件深入理解雨流计数的计算流程也可替换为自定义载荷数据完成循环计数和损伤评估是一份即取即用的工程参考资料。1. 雨流计数法把乱糟糟的载荷时间历程变成损伤账单做结构疲劳的工程师手里拿到一段载荷时间历程十有八九要先过一遍雨流计数法再交给 Palmgren-Miner 线性累积损伤去估寿命。这个方法入门门槛不高用 MATLAB 写一个能跑的实现只要几十行但计数规则、半循环处理和参数设置里的细节能直接影响寿命评估结果。雨流计数法做的事可以概括成一句话把一段任意波形的应力或应变时间序列拆成若干个完整的循环幅值、均值、次数丢掉时间顺序只留下对损伤统计有意义的“账单”。它不关心载荷从哪里来只关心材料的迟滞回线在哪里闭合。适合做疲劳分析、结构健康监测、台架试验数据处理也适合正在学材料力学和信号处理的同学拿来练手。2. 雨流计数法计数逻辑峰值谷值、三点法与四点法2.1 为什么计数对象是迟滞回线而不是峰值材料在循环载荷下的应力应变曲线会形成迟滞回线一个完整的载荷循环对应一条闭合的滞回环。疲劳损伤的本质是这些滞回环反复作用下微裂纹的萌生和扩展。如果只统计载荷峰值会丢失“从哪来到哪去”的信息同一个峰值配上不同的谷值造成的损伤差异很大。雨流计数法的核心思想就是把随机载荷里所有能形成闭合滞回环的路径识别出来。想象雨水从峰顶和谷顶往下流遇到更低的谷或更高的峰就停止一条“雨流路径”就是一个半循环或全循环。这个比喻帮助理解但写成程序时通常用三点法或四点法做数值判定。四点法更直观也是 ASTM E1049 标准中推荐的实现思路本文后面的 MATLAB 代码采用四点法。2.2 极值序列压缩去掉所有中间点雨流计数不关心每个采样点的具体数值只关心局部极值。一个波峰和一个波谷之间经过了多少个点对循环提取没有影响。因此第一步是把原始时间序列压缩成峰谷交替的极值序列。压缩规则非常简单相邻两个差分符号发生变化的点就是极值点符号不变的点全部删掉。function ext extract_extrema(x) % EXTRACT_EXTREMA 抽出载荷时间历程中的全部峰谷值 % x: N x 1 或 1 x N 的载荷向量 % ext: 只包含局部极大、极小值的行向量 x x(:); n numel(x); ext x(1); for k 2:n-1 d1 x(k) - x(k-1); d2 x(k1) - x(k); if d1 * d2 0 ext(end1) x(k); end end ext(end1) x(n); % 删除相邻重复值避免平台段产生误判 keep [true, diff(ext) ~ 0]; ext ext(keep); end这段代码的关键是d1 * d2 0这个判据。载荷从上升到下降、或从下降到上升的转折点两个差分乘积必然为负反之同向变化时乘积为正说明该点不是极值点丢弃。首尾两点无条件保留因为它们是半循环的起点或终点。最后的diff(ext) ~ 0处理等值平台比如长时间保持同一个载荷值的情况避免后续四点法把相等的点当成有效峰谷。2.3 四点法配对规则与循环提取压缩后的极值序列是峰谷交替的相邻两点之间的距离代表一个“半循环跨度”。四点法每次取连续的四个极值点 A、B、C、D判断 B 到 C 这一段是否被 A-B 和 C-D 两段完全夹住。如果|B-C| |A-B|且|B-C| |C-D|说明 B-C 是一个被夹在中间的完整循环把它取出来记录幅值和均值然后删掉 B、C 两个点重新组合序列继续判断。如果不满足就删掉 A窗口向右移动一个点继续找。这个逻辑模拟的是中间那个小循环没有碰到更大的“屋顶”或“地板”它在更大的迟滞回线内部闭合先被提取出来。删除 B、C 之后剩下的序列继续参与更大循环的匹配保证嵌套的小循环先计数、大循环后计数。2.4 半循环的三种去向四点法循环结束后极值序列里往往会剩下 1 到 3 个点。这些点不能成对提取形成半循环。半循环不是噪声它们同样消耗疲劳寿命处理方式决定了寿命评估的保守程度。处理方式做法适用场景镜像拼接法把残余极值序列翻转后接到原序列尾部重新跑一次计数多数半循环会配成完整的周期性载荷、载荷谱较长时推荐系数折减法半循环按 0.5 次全循环计入线性累积损伤快速估算、损伤量级判断直接丢弃忽略全部半循环不推荐会明显低估损伤实际工程中我一般会把残余半循环按 0.5 次计入损伤同时在代码里保留半循环的幅值和均值这样后续如果需要做闭合处理数据还在不用重新数一遍。3. 用 MATLAB 实现雨流计数核心函数与最小可运行代码3.1 输入数据规范与预处理写给雨流计数的输入理论上只需要一列数值向量不需要采样时间。需要注意的是输入序列应该已经去除了趋势项和异常跳变。实测数据里常见的毛刺、掉线导致的尖峰会让极值提取阶段产生大量假峰谷直接扭曲后面的循环统计。一般先做一次中值滤波或滑动平均滤掉高频噪声后再进入极值提取。MATLAB 里加载数据通常来自.mat文件或 Excel读进来后统一转成列向量。如果数据量大到几百万点以上建议用single精度存储载荷值雨流计数对精度要求不高double白白浪费一倍内存。3.2 四点法主循环rainflow_countfunction [cycles, half_cycles] rainflow_count(ext, gate) % RAINFLOW_COUNT 四点法雨流循环提取 % ext: extract_extrema 输出的峰谷序列 % gate: 幅值门限小于门限的全循环不记入输出 % cycles: 全循环表每行 [幅值, 均值] % half_cycles: 剩余半循环表每行 [幅值, 均值] if nargin 2 gate 0; end ext ext(:); cycles zeros(0, 2); half_cycles zeros(0, 2); pts ext; while numel(pts) 4 A pts(1); B pts(2); C pts(3); D pts(4); if abs(B - C) abs(A - B) abs(B - C) abs(C - D) amp abs(B - C) / 2; mean_val (B C) / 2; if amp gate cycles(end1, :) [amp, mean_val]; end pts(2:3) []; % 取出 B-C 段删除这两个极值点 else pts(1) []; % 当前四点不满足配对丢弃最左点 end end for k 1:numel(pts)-1 half_cycles(end1, :) [abs(pts(k1) - pts(k)) / 2, ... (pts(k) pts(k1)) / 2]; end endpts(2:3) []是整段逻辑的核心。删除 B、C 之后原本的 A、D 变成新的相邻极值它们之间可能还有更大的循环需要匹配所以不能让窗口直接跳到后面。这里不额外做“回退”操作因为删除后序列重组下一轮while自动从新的起点开始检查等价于回退了两步。pts(1) []则是在不满足配对条件时把窗口整体右移一个点继续寻找新的四点组合。gate参数控制的是输出过滤而不是删除极值点。也就是说小于门限的循环仍然会被取出并从序列中删除只是不记录到cycles里。这样做的意义在于小循环不参与大循环的配对但它本身也不应该进入损伤统计。3.3 完整调用示例% 原始载荷序列包含嵌套的小循环 x [0, 10, 2, 8, -4]; % 第一步压缩成峰谷序列 ext extract_extrema(x); % 第二步雨流计数 [cycles, half_cycles] rainflow_count(ext, 0); disp(全循环:); disp(cycles); disp(半循环:); disp(half_cycles);这个例子的极值序列本身就是[0, 10, 2, 8, -4]。四点法第一步取 A0, B10, C2, D8|B-C|8|A-B|10|C-D|6条件810成立但86不成立所以删掉 A。第二步取 A10, B2, C8, D-4|B-C|6|A-B|8|C-D|12两个条件都成立提取全循环幅值为(8-2)/23均值为(82)/25。删除 B、C 后剩下[10, -4]记一个半循环幅值 7均值 3。跑出来的结果应该是全循环: 3 5 半循环: 7 3这个简单的输入输出对用来验证代码正确性非常方便手算和程序输出完全一致再接真实数据。3.4 输出结果的单位与方向雨流计数输出的幅值严格来说是载荷差的一半。如果输入是应变输出的就是应变幅如果输入是应力输出的就是应力幅。mean_val是 B、C 两点的代数平均值代表这个循环的均值载荷。疲劳分析中均值对损伤影响很大同样的幅值叠加在拉均值上比叠加在压均值上更危险所以均值不能丢弃。有一点容易踩坑如果输入的是应力范围即直接采集到的峰值减谷值那输出的amp就是应力范围的一半。后续查 S-N 曲线时要确认曲线横坐标是幅值还是范围单位搞错寿命能差出几个数量级。我习惯在代码注释里写明输入的类型避免半年后回来看代码时还要猜。4. 雨流计数法的工程化参数门限、内存与批量处理4.1 幅值门限 gate 怎么定实测载荷里永远有大量小幅值抖动不一定是真实损伤更多是测量噪声。如果统统计入损伤寿命评估会偏保守如果门限设得太大又把真实的小损伤循环滤掉了。通常的做法是先看载荷幅值分布取最大幅值的 1% 到 10% 作为门限具体数值取决于传感器的信噪比和结构的疲劳敏感度。load_data load(channel_1.mat).channel_data; max_amp max(abs(load_data)); gate 0.02 * max_amp; % 按最大幅值的2%过滤 ext extract_extrema(load_data); [cycles, half_cycles] rainflow_count(ext, gate);这一段代码里gate的选择依据是最大幅值而不是均值或标准差因为雨流计数关心的是局部循环的跨度用全局统计量做门限容易误删大均值循环。实际项目中我会先跑一遍不带门限的计数画幅值直方图观察底部噪声平台再反过来确定门限。这个流程比拍脑袋设一个百分比更可靠。参数推荐做法常见误区gate最大幅值的 1%~10%设太大丢失真实损伤循环数据类型single用 double 存储大数据内存翻倍首尾处理保留半循环并计入损伤直接丢弃半循环寿命偏乐观循环统计幅值、均值、次数三元组只存幅值均值信息丢失4.2 大数据量的内存与性能处理一段 24 小时、采样率 100 Hz 的载荷数据约 864 万点用 double 存储需要 64 MB看起来不大但如果是 32 通道同时采集就是 2 GB直接爆内存。雨流计数不关心时间戳所有通道可以先各自压缩成极值序列再计数。实测载荷的极值点通常只占原始采样点的 5% 到 10%压缩这一步能省掉九成内存。data single(load(mat_file.mat).data); % 转 single 省一半内存 ext extract_extrema(data(:)); [cycles, half] rainflow_count(ext, gate);如果数据量大到连single都撑不住就分块读取。注意分块不能简单地把每块独立计数因为块与块交界处的极值会丢失前后关系导致循环计数不完整。常见做法是按通道整列读取或者对单一大文件用matfile对象做部分读取再手工拼接极值序列。拼接的规则很简单把上一块最后一个极值和下一块第一个极值做一次合并如果两者之间没有新的极值就直接接上。4.3 批量处理多通道载荷台架试验常常一次记录十几个通道的信号逐通道复制粘贴计数代码显然不合理。写一个循环遍历通道文件把结果汇总到一个结构体里。files dir(channel_*.mat); all_results struct(); k 1; for f files ch_data load(fullfile(f.folder, f.name)).channel; ch_data single(ch_data(:)); ext extract_extrema(ch_data); max_amp max(abs(ch_data)); [all_results(k).cycles, all_results(k).half_cycles] ... rainflow_count(ext, 0.02 * max_amp); all_results(k).channel_name f.name; k k 1; end每个通道的门限按各自的最大幅值独立计算避免一个通道的强载荷把另一个弱通道的真实循环全部滤掉。保存结果时我一般会同时存原始数据的文件名和门限值方便追溯。dir返回的结构体在for循环里按列遍历所以上面用了f files转置成行这是 MATLAB 新手容易卡住的一个点。4.4 结果可视化幅值均值三维直方图计数完成后直接看矩阵很难发现规律画一张三维直方图能快速了解载荷的分布形态。figure; histogram2(cycles(:,1), cycles(:,2), 30, DisplayStyle, bar3); xlabel(幅值); ylabel(均值); zlabel(循环次数); view(45, 30);histogram2的第一个参数是幅值第二个是均值第三个是分箱数量。bar3样式把每个箱子画成一根柱子适合观察循环次数的集中区域。如果发现循环全部集中在小幅值区域说明结构长期处于低应力工作状态疲劳寿命主要受少数大循环控制如果大幅值区域也有明显分布则载荷谱中存在典型的重载工况。画图前最好先统计一下数值范围把异常大的离群点剔除否则坐标轴会被拉得很长三分之二的图都是空的。5. 从雨流矩阵算疲劳寿命并验证代码正确性5.1 用 Palmgren-Miner 线性累积损伤把循环变成寿命雨流计数输出的是循环清单寿命评估还要靠累积损伤法则。最简单的形式是 Palmgren-Miner每个循环的损伤是 1/NN 由 S-N 曲线给出把所有循环的损伤相加得到总损伤总损伤达到 1 时认为发生疲劳破坏。function D rainflow_miner_damage(cycles, half_cycles, C, m) % RAINFLOW_MINER_DAMAGE 基于 Basquin 公式的线性累积损伤 % S-N 曲线形式: N C * S^(-m) % cycles: N x 2 矩阵 [幅值, 均值] % half_cycles: M x 2 矩阵 [幅值, 均值] D 0; for i 1:size(cycles, 1) S cycles(i, 1); N C * S^(-m); D D 1 / N; end for i 1:size(half_cycles, 1) S half_cycles(i, 1); N C * S^(-m); D D 0.5 / N; end endC和m是 S-N 曲线的材料常数不同材料差异巨大必须从材料手册或试验获取。半循环按 0.5 计入这是一个偏保守的处理方式工程上可接受。均值在这里没有参与计算因为 Basquin 公式不带均值修正如果均值变化大需要改用 Goodman 或 Gerber 修正把循环的均值折算到等效零均值幅值上。这里的S用的是幅值如果 S-N 曲线横坐标是应力范围就要传入2 * cycles(i, 1)这是最容易搞错的一步。5.2 手工验证用例验证雨流计数代码用真实载荷反而难判断对错。推荐用手工可推的极值序列。下面这组是经过精心设计的最小用例ext_test [0, 10, 2, 8, -4]; [cycles, half] rainflow_count(ext_test, 0);手推过程第一次四点检查(0, 10, 2, 8)不满足配对条件第二次检查(10, 2, 8, -4)满足条件提取循环幅值 3均值 5剩余[10, -4]构成半循环幅值 7均值 3。如果代码输出的数值与此一致说明四点法和删除重组的逻辑基本正确。在此基础上再测试一段真实载荷观察循环数量级是否合理就能发现多数实现上的边界问题。5.3 边界条件与报错排查雨流计数代码报错主要集中在极值提取和空数组处理两个地方。extract_extrema对长度小于 3 的序列会产生空结果或单点结果调用前先加长度判断。另外如果载荷全程单调上升或下降没有峰谷交替四点法循环一次都不会执行cycles为空但half_cycles仍然会输出一段半循环这是正常的不要当成 bug。常见的一个隐蔽错误极值序列里出现相邻相等点比如[1, 1, -2, 3]会导致d1 * d2 0判断失效。extract_extrema里的diff(ext) ~ 0去重就是专门防这个问题的。另一个容易忽略的是gate设负值amp gate恒成立等于没有过滤写参数校验时建议assert(gate 0)。调试时把pts数组在while循环里打印出来对照手推过程一步一步看比反复猜测哪里出错快得多。本文还有配套的精品资源点击获取