做疲劳分析的人对“雨流计数法”这个词应该都不陌生。我做结构耐久性评估那阵子天天跟随机载荷谱打交道被问得最多的就是这段乱七八糟的实测载荷到底能折合成多少个疲劳循环这个问题不搞清楚后面S-N曲线、Miner累积损伤全是空中楼阁。我自己的解决方案是用MATLAB手写一套雨流计数工具分别实现了三点法和四点法这篇就完整复盘一下这两套算法的思路、代码实现和实际使用中踩过的坑。这个内容适合机械结构、车辆、航空、风电等领域做疲劳分析的工程师也适合做载荷谱数据处理的研究生。你不需要有很高深的MATLAB水平只要会基本的数组操作就能看懂。我会把算法逻辑拆开讲清楚再给出可以直接复现代码最后分享几个我用下来觉得特别值钱的注意事项。1. 先搞明白雨流计数法到底在数什么很多人一上来就找代码其实代码反而最简单真正容易卡住的是理解不完备。我先用大白话把原理讲透。1.1 疲劳循环的本质迟滞回线与损伤对应关系材料在循环载荷作用下应力-应变关系会画出一个闭合的迟滞回线。每画出一个完整的封闭回线材料就经历了一次完整的疲劳损伤事件对应一个循环。工程上最头疼的问题是实测的载荷时间历程往往是乱的今天大峰套小谷明天小谷里藏大包根本没法直接数出来究竟有多少个完整回线。雨流计数的核心思想就是把这个杂乱无章的时域波形通过“雨滴从屋顶流下”的隐喻拆解成一系列封闭的应力循环和半循环。为什么叫雨流想象一个宝塔状屋顶剖面雨滴从内侧往下流遇到比起始点更低的屋檐就会滴落。这个比喻对应到载荷谱上就是从峰值出发追踪载荷的路径直到你追到的点低于起始点才停止每一个完整路径代表一个循环。这个方法的厉害之处在于它是从局部应力应变行为出发的跟材料实际承受的损伤机制是对应的所以提取出来的循环数据可以直接喂给疲劳损伤公式用。1.2 为什么这项技术在工程里这么重要你去看一份真实的随机载荷谱比如风电叶片根部弯矩、汽车悬架弹簧的垂向力、飞机机翼的过载谱它们几乎不可能像教科书里那种等幅正弦波一样规律。如果拿这个波形直接跟S-N曲线对比你会发现完全无法判断“这算几个循环”。雨流计数法就是用来解决这个问题的标准化工具。它把不规则的载荷-时间历程转化为一系列等效的恒幅循环输出每个循环的幅值、均值和发生次数。这样你就可以用Miner线性累积损伤公式把每个循环的损伤量累加起来得到总损伤进而估算疲劳寿命。在工程产品开发流程里这是一道绕不开的工序。载荷谱采集完之后的下一步几乎必然是雨流计数。如果这一层处理得不严谨后面所有仿真和试验都是在沙滩上盖楼。1.3 为什么偏偏用MATLAB实现市面上确实有各种雨流计数工具箱也有商业软件自带模块但我在实际工作中还是更喜欢自己用MATLAB写一套。原因有三个。第一MATLAB处理数组和循环逻辑非常直观雨流计数的核心操作其实就是数组元素的比较、删除、重排用MATLAB写出来代码非常短调试也方便。第二我需要把计数结果直接跟载荷谱数据、有限元结果、寿命计算公式无缝衔接用MATLAB做这些数据流转是最顺手的不用来回倒格式。第三雨流计数算法本身并不复杂尤其是三点法和四点法逻辑清晰MATLAB实现一遍之后你能完全掌控每一步细节不会被工具箱的黑盒子困住。这一点在写论文或者审疲劳报告的时候特别重要别人问起来你能解释清楚。2. 动手之前把原始载荷变成干干净净的峰谷序列我在写算法之前走了不少弯路一开始直接把原始采集数据扔进计数函数里结果惨不忍睹。后来才明白雨流计数必须建立在预处理过的峰谷序列上。2.1 预处理第一步去平台、去毛刺实测数据最常见的两个问题是零漂和噪声。零漂导致载荷数据整体往上漂或者往下飘如果不处理计数出来的均值会失真。噪声会导致载荷曲线出现大量细小的锯齿这些锯齿会被误识别成循环最终放大损伤。第一步工作是去除平台段也就是把连续相等的点只保留一个。为什么要这么干因为雨流算法在判断峰谷时依赖差分的符号变化如果中间有一段水平线段符号变化会被破坏产生错误的转折点。我在代码里的处理比较直接用diff判断相邻两个点的差值是否约等于0小于阈值就跳过。第二步是去毛刺。通常我会先用一个很小的去幅阈值把小于设定阈值的波动过滤掉。举个例子如果整个载荷范围是正负50千牛那我可能会把0.5千牛以下的微小波动先抹平。这一步一定要格外小心阈值设大了会把真实存在的小载荷循环也删掉导致损伤被低估。我的经验是先画图看一遍载荷谱的波动尺度再定阈值不要拍脑袋。2.2 提取峰谷只有转折点才是雨流的骨架去完平台和毛刺后数据中间可能还有大量单调上升或下降的中间点。这些点对雨流计数毫无意义真正的关键是那些转折点峰值和谷值。提取峰谷的标准做法是检查每个点与前后两个点的关系如果中间点同时大于前后两个点或者同时小于前后两个点那就是一个转折点。用MATLAB写就是判断两个差分值的乘积是否小于零。这一步得到的基本就是雨流计数的“骨架”。需要注意序列的首尾两个点必须保留。它们不是转折点但它们是整个载荷历程的边界对残余循环的计算很重要。你说不清首尾这两个点在未来会跟哪个内部点闭合所以先留着交给计数算法去处理。2.3 重排起点从最大峰或最小谷开始预处理完峰谷序列之后还需要做一次重排找到整个序列中绝对值最大的点把它作为新的起点然后把后续数据接上再从原序列开头接到这个最大点之前。我用的判断标准是max(abs(peaks))也就是比较最大峰值和最小谷值的绝对大小。为什么必须这样重排雨流计数法在计数循环时本质上要求从整个载荷历程的最高点或最低点开始。你想一下如果起点不是极值那最高点和最低点之间的那个大循环就会被切到序列的两端一边一半导致一个大循环被拆成两半计数结果全乱套了。这一步是三点法和四点法共同的前提条件千万不能省。3. 三点法实现麻雀虽小五脏俱全三点法是最容易理解的雨流计数实现方式。它的代码短、逻辑直观适合刚接触雨流计数的朋友用来建立直觉。3.1 三点法的判断逻辑三点法的基本思路是依次取峰谷序列中连续的三个点判断这三个点能否形成一个封闭循环。假设三个点分别为a、b、c也就是序列中连续的三段路径。首先判断中间这段路径的方向是不是和前面那段相反也就是(b-a)乘以(c-b)小于零说明这里出现了一个转折是峰或者谷。如果方向没有反转那么就不构成一个潜在的闭合路径继续往下走。判断完方向后还需要比较两段路径的纵向跨度。如果从a到b的跨度小于等于从b到c的跨度那么a到b这一段就可以看作是一个已经被完整包裹的内层循环可以提取出来。提取的依据是在雨流的规则里小的循环会被大的循环包在里面当外层的跨度足够大时里层那一段已经构成了完整的闭合回线。一旦满足提取条件我们就记录一个完整循环循环幅值就是ab段纵向距离的一半均值也就是a和b的平均值。然后从峰谷序列中删掉b点和它前面的a点让新的邻居重新组成三点窗口继续判断。这个逻辑写出来非常简单但它有一个天然的局限窗口内的判断是基于局部信息。对于那种波形极其复杂、大小循环嵌套紧密的数据三点法容易留下较多残余循环并且在某些边界情况下和严格的雨流规则不完全一致。3.2 最简版MATLAB代码我把三点法实现封装成一个函数输入是一维载荷时间序列输出是一个循环统计表格。代码里我特意加了注释方便你自己调试。function [cycles, residual] rainflow_three_point(data, thresh) % rainflow_three_point - 三点法雨流计数 % 输入 % data : 一维载荷时间序列行向量或列向量 % thresh : 去小幅波动的阈值可选默认0 % 输出 % cycles : n行3列的矩阵每行为 [幅值, 均值, 循环类型] % 循环类型取值为1(完整循环)或0.5(半循环) % residual : 剩余峰谷序列未计成完整循环的部分 if nargin 2 || isempty(thresh) thresh 0; end % --- 1. 基本预处理 --- data data(~isnan(data)); data data(:); % 统一成列向量 % --- 2. 提取峰谷 --- peaks data(1); for i 2:length(data)-1 if (data(i) data(i-1) data(i) data(i1)) || ... (data(i) data(i-1) data(i) data(i1)) peaks(end1, 1) data(i); %#okAGROW end end peaks(end1, 1) data(end); % --- 3. 去除相邻等值点和微小波动 --- j 1; for i 2:length(peaks) if abs(peaks(i) - peaks(j)) thresh j j 1; peaks(j) peaks(i); end end peaks peaks(1:j); % --- 4. 重排从绝对值最大的点开始 --- [~, idx_max] max(abs(peaks)); peaks [peaks(idx_max:end); peaks(1:idx_max)]; % --- 5. 三点法计数 --- v peaks; cycles []; while length(v) 3 a v(1); b v(2); c v(3); % 判断 a-b 和 b-c 方向是否相反存在转折 if (b - a) * (c - b) 0 % 如果 ab 跨度 bc 跨度则 a-b 构成完整循环 if abs(b - a) abs(c - b) amp abs(b - a) / 2; mean_val (a b) / 2; cycles(end1, :) [amp, mean_val, 1]; %#okAGROW % 删除 a 和 b窗口回退 v(1:2) []; else % 否则右移一个点再看下一组三连点 v(1) []; end else % 没有转折右移 v(1) []; end end % --- 6. 残余序列按半循环处理 --- residual v; for i 1:length(v)-1 amp abs(v(i1) - v(i)) / 2; mean_val (v(i) v(i1)) / 2; cycles(end1, :) [amp, mean_val, 0.5]; %#okAGROW end % --- 7. 汇总成表格形式并排序 --- if isempty(cycles) cycles zeros(0, 3); else cycles sortrows(cycles, 1); end end这段代码应该是可以直接跑通的。你只需要把载荷序列传进去比如这样的形式data [0 2 -1 3 -2 1 0]; [cycles, residual] rainflow_three_point(data, 0); disp(cycles);输出结果的每一行代表一个被识别出的循环第一列是半幅值第二列是均值第三列是循环类型。注意我这里幅值用的是半幅值也就是循环范围的一半。如果你习惯用整个应力范围做寿命计算要记得乘2。3.3 三点法实测输出与局限用上面的示例数据跑一遍你会发现算法只用了三步就完成了主要计数剩余序列变短了输出的完整循环数量也不算多运行速度非常快。但是我在实际使用中很快发现了三点法的问题它太依赖局部三点之间的跨度关系对波峰波谷嵌套的复杂工况不够灵敏。什么意思呢当一个大型载荷循环内部有许多次生循环时三点法往往会漏掉个别次生循环或者把一些本该继续延展的循环提前封口导致计数结果和商业软件存在偏差。更加严格的说法是三点法其实更适合作为教学演示或者用于数据量非常大、只需要快速估算的场景。如果做正式的项目报告我开始偏向使用四点法。4. 四点法实现工程中更推荐的版本四点法是业界公认比三点法更稳定的雨流计数实现。它多引入一个点看起来只是窗口多了一个元素实际上这个改变让判断闭合循环时多了一个约束条件从而更接近雨流定义。4.1 为什么多一个点就更稳四点法一次取连续四个峰谷点a、b、c、d。它要判断的是中间b到c这一段是否闭合依据是bc跨度不大于ab跨度同时也不大于cd跨度。如果成立说明中间这段不仅被起点一侧包住也被终点一侧包住完全处于一个更大的轨迹内部这时候提取b到c这个循环是稳妥的。多出来的那个d点实际上起到了“外部范围确认”的作用。三点法在判断ab能构成循环时只看了后面的c跨度没考虑再往后还会有可能更大的外框。四点法则把b到c这一段放进了前后两个更大范围的上下文中避免了提前封口。这个特性在处理大循环套小循环的信号时优势特别明显。用一句话概括三点法看“我是不是比你小”四点法看“我前后都是大的所以我完整闭合了”。工程上四点法输出的循环更接近材料真实迟滞回线。4.2 完整MATLAB代码function [cycles, residual] rainflow_four_point(data, thresh) % rainflow_four_point - 四点法雨流计数 % 输入 % data : 一维载荷时间序列行向量或列向量 % thresh : 小幅波动的阈值默认0 % 输出 % cycles : n行3列的矩阵每行为 [幅值, 均值, 循环类型] % residual : 剩余峰谷序列 if nargin 2 || isempty(thresh) thresh 0; end % --- 1. 预处理 --- data data(~isnan(data)); data data(:); % --- 2. 提取峰谷 --- peaks data(1); for i 2:length(data)-1 if (data(i) data(i-1) data(i) data(i1)) || ... (data(i) data(i-1) data(i) data(i1)) peaks(end1, 1) data(i); %#okAGROW end end peaks(end1, 1) data(end); % --- 3. 去等值点和微小波动 --- j 1; for i 2:length(peaks) if abs(peaks(i) - peaks(j)) thresh j j 1; peaks(j) peaks(i); end end peaks peaks(1:j); % --- 4. 重排从绝对值最大的点开始 --- [~, idx_max] max(abs(peaks)); peaks [peaks(idx_max:end); peaks(1:idx_max)]; % --- 5. 四点法主循环 --- v peaks; cycles []; while length(v) 4 a v(1); b v(2); c v(3); d v(4); % 判断bc段是否可以提取为完整循环 if abs(b - c) abs(a - b) abs(b - c) abs(c - d) amp abs(b - c) / 2; mean_val (b c) / 2; cycles(end1, :) [amp, mean_val, 1]; %#okAGROW % 删除 b 和 c窗口回退 v(2:3) []; else % 右移一个点继续检查 v(1) []; end end % --- 6. 残余序列按半循环处理 --- residual v; for i 1:length(v)-1 amp abs(v(i1) - v(i)) / 2; mean_val (v(i) v(i1)) / 2; cycles(end1, :) [amp, mean_val, 0.5]; %#okAGROW end if isempty(cycles) cycles zeros(0, 3); else cycles sortrows(cycles, 1); end end这里最关键的地方在主循环的条件判断abs(b-c) abs(a-b) abs(b-c) abs(c-d)。实践中我发现如果数据已经重排并且预处理干净四点法的收敛速度通常比三点法更快因为每次提取循环时删掉两个点序列长度减少得更快。4.3 三点法和四点法结果怎么对齐在实际项目中我经常需要对比不同计数方法带来的寿命估算差异。同一个载荷谱分别用三点法和四点法跑结果并排放在一起你会发现三点法输出的循环数量可能更多但很多循环幅值偏小四点法输出的循环更“整”更接近手画迟滞回线的感觉。为什么会有这种差异因为三点法把一些实际未闭合的半循环提前算成了完整循环或者把若干小循环拆分得更碎。四点法由于多了外框约束小循环只有在确认被包住之后才会提取。我的建议是报告里写清楚用的是什么算法。两种方法都是允许的但你不能这个项目用三点法下个项目用四点法然后寿命结果直接做横向对比。如果真要对比建议至少用四点法作为基准因为它在复杂载荷谱下的稳定性更好。5. 实操中躲不开的坑与排查技巧代码写完只是第一步真正用起来你会发现坑全在后头。我把自己踩过的问题和排查经验整理一下希望能帮你少走弯路。5.1 数据边界问题处理不好会漏记半循环雨流计数最容易被忽略的是首尾数据边界。整个载荷时间序列不可能刚好跑完一个完整的闭合循环最后一定会有残余路径。这段残余路径对应的是半循环半循环对疲劳损伤也有贡献。我一开始的处理方法特别粗暴直接把残余序列倒过来再跑一遍结果发现输出的循环数量和手算对不上。后来才想明白正确的做法是把残余序列按照相邻转折点的差值记录成半循环也就是我上面代码里第6步的处理方式。这一步如果漏掉总的循环累计损伤会偏低在某些极端工况下偏低幅度能超过10%。5.2 等值点和浮点误差小心让判断条件失效实测采集的载荷数据经常会出现相邻两个值完全相等的情况。如果直接用(data(i) data(i-1) data(i) data(i1))去判断转折点遇到平台段就会判断失败损失真实转折点。我的做法是先用阈值把等值点压掉同时在比较中允许浮点误差。MATLAB的浮点运算有时候会给你带来莫名的差异比如两个本应相等的浮点数计算后可能差出1e-15导致判断条件不成立。这种问题不太好查最快的办法是对数据整体做一次中心化或者归一化把数值范围压到比较舒适的区间比如正负1浮点误差问题会缓解很多。5.3 如何验证计数结果是对的写完算法之后我强烈建议你用一个已知答案的简单载荷谱去验证。比如取一个典型的正弦变幅载荷谱手动画迟滞回线数出有哪些闭合循环然后跟程序输出对比。我自己的验证方法有三个。第一画迟滞回线对照图把载荷时程画出来同时把程序提取的每个循环用矩形或者平行四边形叠加画上去肉眼看是否能对应上。第二循环总数守恒验证。原始峰谷序列如果看作相邻点差值的总和那么程序提取出的所有循环和半循环的范围总和应该和原始序列的总变程保持一致。我写过一个辅助函数专门校验这个守恒关系一旦对不上说明计数逻辑里有bug。第三和商业软件做一次对标。拿同一段数据丢到成熟工具里跑一遍对比循环幅值分布直方图和累计疲劳损伤。偏差在百分之几以内说明实现正确偏差一大就需要回头查边界处理。这里我整理了一个常见问题排查表基本涵盖了我这些年遇到的主要问题。现象可能原因排查方法提取的循环数量明显偏多数据中的噪声没有被有效滤除先画载荷时程图观察锯齿密度适当调大去幅阈值结果与商业软件对不上起点重排规则不一致确认是否从最大绝对值点开始重排幅值整体偏小等值点没有被压缩检查峰谷序列中是否存在相邻相等点半循环统计异常残余序列处理方式不同手动从剩余序列末尾倒推一遍确认每条半循环是否正确浮点判断时好时坏数据数量级跨度过大对载荷统一做归一化处理后再计数5.4 性能优化处理超长载荷谱时的建议实际工程载荷谱数据量往往很大动辄几十万个点甚至上百万个点。上面这种基于数组删除和重建的算法在数据量特别大的时候会拖慢速度因为每次删元素MATLAB都要重新分配内存。我的建议是先用峰谷提取把数据压缩到只剩转折点通常能把数据量减少一个数量级然后计数循环里的判断次数也会大幅降低。如果数据还是特别大可以考虑把数组索引操作改成while循环加左右指针的方式避免频繁删除重建数组。不过对于大多数项目和论文场景峰谷压缩后再跑四点法性能完全够用。这段路径我是一步步踩出来的。最开始只会套公式后来自己推演迟滞回线再看标准里的计数步骤才真正把三点法和四点法吃透。我自己现在的主力工具是四点法版本日常快速预扫一遍载荷谱用三点法辅助对比。最后再分享一个小技巧如果你想快速检验自己手上的载荷谱到底有多少个有效循环先用峰谷提取压缩数据再用四点法跑速度会快得让你怀疑人生。