简介一套基于四点法的雨流计数MATLAB实现聚焦疲劳分析中不规则应力/应变时间序列的循环载荷提取。适用于材料、机械、航空航天等领域的工程师和科研人员尤其适合需要处理实测载荷谱并计算疲劳寿命的初学者参考。压缩包仅包含1个m脚本文件体积仅1KB便于直接阅读与二次修改。代码覆盖极值点识别、循环匹配、半径与中心点计算、归一化处理以及重复循环去除等关键步骤最终输出可排序的归一化循环对方便不同工况载荷谱之间的横向比较理解后既能掌握四点法雨流计数的完整实现逻辑也可进一步嵌入更完整的疲劳寿命评估流程。目前已有1742人学习下载整体适合作为MATLAB实现四点法雨流计数的入门样例与参考脚本。1. 雨流计数的四点法实现从载荷谱里抽出闭合滞回环一批随机载荷时程放到你面前第一条直觉通常是直接统计峰谷和幅值分布。真正落到疲劳寿命评估时这个直觉会立刻失灵同一个应力幅值如果出现在封闭滞回环里它对损伤的贡献是一整圈的如果只是半程爬坡贡献就要减半。雨流计数要做的就是把一条连续的载荷时间序列拆成一个个完整的闭合循环而四点法实现是这套逻辑里收敛最快、也最容易写出问题的一种写法。它用四个相邻转折点比较三段行程一旦中间两点的行程小于两边就判定这一段闭合并从序列中抽走。下面按四点法实现这条线展开先讲判据和计数流程再给一份能直接跑出循环记录的最小 Python 实现接着把阈值、峰谷容差和残波处理这三个参数讲透最后用正弦叠加噪声的自检序列验证计数结果。这套东西适用于做疲劳耐久、结构强度或载荷谱处理的人也适合想理解 Rainflow 内部到底在算什么的人。2. 四点法判据相邻四点与三段行程的闭合条件2.1 四点法判据Δ1、Δ2、Δ3 的大小关系四点法每次从转折点序列里取连续的四个点 x1、x2、x3、x4计算三段行程差Δ1 |x2 - x1|Δ2 |x3 - x2|Δ3 |x4 - x3|当 Δ2 ≤ Δ1 且 Δ2 ≤ Δ3 时认为 x2 到 x3 这一段构成了一个可以抽出的闭合循环候选。这个条件的物理含义很直白中间段是三个相邻段里最短的一段它两侧的行程都比它长说明载荷在 x2→x3 的反向段里已经足够形成一个完整的滞回环。注意中间点 x2 还必须是峰或谷也就是 (x2 x1 且 x2 x3) 或者 (x2 x1 且 x2 x3)。如果 x1、x2、x3 三点单调递增或递减那 x2 只是爬坡途中的一个点d2 ≤ d1 也会成立但这不是一个闭合循环抽出来会浪费后续配对机会。这个极值限定是四点法实现里最容易漏掉的条件漏掉的后果是在单调长序列中间误抽大量小循环。满足条件后循环的幅值取 Δ2均值取 (x2 x3) / 2。为什么不取 x1 和 x4因为闭合滞回环里 x1 和 x4 并不一定严格相同取 x2、x3 代表的是这一段回程的真实应力水平x1 和 x4 只是用来限定闭合条件的边界。把 x2、x3 从序列里抽走后剩下的序列继续保持原有顺序下一次循环从抽走位置的前一个点开始重新组合。2.2 三点法与四点法判断粒度与误判风险雨流计数还有一个常见写法是三点法只在三个点里判断相邻两段行程。三点法的判定条件是 |x2 - x1| ≤ |x3 - x2|当当前点 x2 相对 x1 的反向行程足够短时就把 x1→x2 作为半个候选段处理。三点法在大量工程实现里也能用但它比较的是两段行程遇到连续同向走一段再回头时很容易把尚未闭合的半循环提前当成完整循环记录下来。四点法多取了一个边界点等于多等了一段数据再下结论对趋势斜率和慢变分量的容忍度更高误抽率明显低于三点法。对比项三点法四点法参与判断的转折点数34行程段数23闭合判定后段 ≤ 前段中段 ≤ 前段且中段 ≤ 后段对单调趋势的误判较高较低典型应用快速在线计数离线载荷谱、疲劳寿命评估四点法并不是三点法的简单加长。三点法把计数粒度压在两个点之间四点法则把判定粒度压在三个点之间要让中间段成为局部最小值这天然过滤掉了单向漂移带来的伪循环。实际处理长时程信号时四点法比三点法多出的这个边界条件能显著减少小环的重复计数。2.3 四点法的计数流程与伪代码完整的四点法雨流计数流程可以拆成五步。第一步把原始时程压缩成转折点序列去掉单调段里中间冗余点。第二步从序列头开始取四个连续点。第三步检查中间点是否为峰或谷以及中间段是否同时短于两侧行程。第四步满足条件则抽取 x2、x3 为一个循环记录并从序列中删除这两个点扫描位置回退一位不满足则前进一位。第五步序列剩余不足四点时停止剩余的点作为残波保留。seq extrema[:] # extrema 为预处理后的转折点序列 cycles [] i 1 # 从第二个点开始找候选段 while i len(seq) - 2: x1, x2, x3, x4 seq[i-1], seq[i], seq[i1], seq[i2] d1 abs(x2 - x1) d2 abs(x3 - x2) d3 abs(x4 - x3) is_peak (x2 x1) and (x2 x3) is_valley (x2 x1) and (x2 x3) if (is_peak or is_valley) and d2 d1 and d2 d3: cycles.append(((x2 x3) / 2.0, d2)) del seq[i:i2] # 抽走 x2 和 x3 i max(0, i - 1) # 回退一位重新拼段 else: i 1这段伪代码的逻辑核心在del seq[i:i2]和i max(0, i - 1)。抽走两个点后原本隔开的 x1 与 x4 变成相邻点它们和更前面的点可能组成新的闭合循环所以扫描位置必须回退一位重新判断。如果回退到 0下一轮就从头开始拼四个点。i从 1 而不是 0 开始是为了让seq[i]始终作为候选循环的中间段起点。这里记录的d2可以进一步乘上材料常数换算成应力幅均值则用于后续的平均应力修正。3. 四点法雨流计数的 Python 最小实现预处理、主循环与输出3.1 峰谷预处理把原始时程压缩成转折点原始载荷时程通常是高采样率的连续曲线甚至带有大量轻微抖动。直接丢进四点法循环会产生两类问题一是单调上升或下降段里的中间点会被反复判断白白消耗计算二是小幅抖动会制造大量幅值很小的伪循环把统计结果淹没。所以计数前必须先做峰谷压缩只保留方向真正发生反转的转折点。def to_extrema(series, tol1e-9): 把原始载荷时程压缩成转折点序列。 tol 为峰谷容差相邻两个采样点差值小于 tol 时视为无效抖动。 ext [] for v in series: if ext: d v - ext[-1] if abs(d) tol: continue if len(ext) 2 and (ext[-1] - ext[-2]) * d 0: # 与上一段方向相同替换末点保留更大行程 ext[-1] v else: ext.append(v) else: ext.append(v) return ext函数里最关键的是方向判断(ext[-1] - ext[-2]) * d 0表示新点与最近一段行程同向这时用新点替换掉ext[-1]保留更远的行程末端异向则说明方向反转新点作为新的转折点入栈。这样处理完序列里相邻两点一定是交替升降的。tol参数用来屏蔽采集噪声它比较的是相邻采样的绝对值差不是累计位移差。对采样率 1kHz 的应变信号tol通常设在应变分辨率的两到三倍具体值取决于传感器量程和 ADC 位数。3.2 四点法主循环与可运行代码预处理之后进入四点法主循环。下面这份代码可以直接复制运行输入一串载荷数组输出闭合循环列表和剩余残波序列。def rainflow_four_point(extrema, threshold0.0): 四点法雨流计数主循环。 threshold 为循环幅值阈值小于阈值的循环不记录 但仍然会从序列中抽出以保证后续拼接正确。 seq list(extrema) if len(seq) 4: return [], seq cycles [] i 1 while i len(seq) - 2: x1, x2, x3, x4 seq[i-1], seq[i], seq[i1], seq[i2] d1 abs(x2 - x1) d2 abs(x3 - x2) d3 abs(x4 - x3) is_peak (x2 x1) and (x2 x3) is_valley (x2 x1) and (x2 x3) if (is_peak or is_valley) and d2 d1 and d2 d3: if d2 threshold: mean (x2 x3) / 2.0 cycles.append((mean, d2)) del seq[i:i2] i max(0, i - 1) else: i 1 return cycles, seq if __name__ __main__: series [0, 2, 5, 3, 8, 1, 6, 4, 9, 2] ext to_extrema(series) cycles, residual rainflow_four_point(ext) print(转折点序列:, ext) print(闭合循环 (均值, 幅值):, cycles) print(残波:, residual)主循环里threshold只控制是否记录不控制是否抽出。很多第一次写雨流计数的工程师会把小于阈值的循环直接留在序列里这样做的坏处是这些小段不闭合会在后续拼接时改变相邻点的配对关系最终得到的循环幅值整体偏离真实值。正确的做法是先抽除、再判断是否记录代码里del seq[i:i2]始终执行if d2 threshold只包住cycles.append。3.3 运行结果与循环记录的三个字段运行上面代码输出大致如下转折点序列: [0, 5, 3, 8, 1, 6, 4, 9, 2] 闭合循环 (均值, 幅值): [(4.0, 2.0), (3.5, 5.0)] 残波: [0, 6, 4, 2]第一条记录(4.0, 2.0)来自 5→3 这一段均值 4、幅值 2第二条记录(3.5, 5.0)来自 8→1 这一段。残波[0, 6, 4, 2]是序列剩余不足四点时留下的它不再闭合需要后续单独处理。每一条循环记录都包含均值、幅值两个核心字段必要时还要加上计数时刻、起始索引或载荷方向用于后处理中计算等效损伤、平均应力修正和载荷谱编制。关于循环记录的字段经验上建议保留原始索引或时间戳。雨流计数输出的循环不再携带时间顺序信息但如果要把它映射回某一时段的载荷特征例如晚间大跨度桥梁的涡振、风电机组的阵风工况索引能帮你快速定位到原始时程的对应区间否则后处理阶段只能按统计分布做事后归因。4. 雨流计数参数阈值、峰谷容差与残波循环的取舍4.1 循环幅值阈值过滤低幅振动的正确位置循环幅值阈值threshold是最直观的参数但它的放置位置比数值本身更重要。很多人在拿到数据后先用滑动平均滤波把高频小振幅信号抹平再计数。这个做法会改变峰值的真实高度因为滑动平均是双向平均会把原本的尖峰压低也会把陡峭的谷底抬高最终得到的循环幅值整体偏小。比较稳妥的顺序是原始时程 → 峰谷压缩 → 四点法抽取 → 阈值过滤。先让四点法把所有闭合循环都找出来再按幅值阈值决定哪些进入损伤统计这样峰值高度没有被提前改变。阈值大小的取法要结合疲劳曲线的截断范围。如果材料的 S-N 曲线在低应力段有无限寿命水平低于该应力幅的循环对损伤没有任何贡献阈值就可以设在这个水平对应的载荷值附近。实际工程里我一般会把阈值放在最大循环幅值的 1% 到 5% 之间再根据低幅循环的数量占比微调。如果低幅循环数量占比超过 90%说明峰谷压缩里的tol太小了应该先去处理抖动而不是靠阈值硬切。4.2 峰谷容差抖动的清洗发生在计数前峰谷容差tol控制的是转折点序列的生成粒度。它和循环阈值有两个本质区别第一tol作用于相邻采样点不作用于循环幅值第二tol改变了四点法看到的数据形态而threshold只改变统计输出。tol设置过小采集信号里的量化噪声会残留在转折点序列中四点法会把每个微小起伏都配对成循环循环总数可能膨胀几十倍tol设置过大真实的结构载荷峰值会被压缩掉特别是那些斜率较缓的长周期大振幅循环转折点还没来得及记录就已被后续同向点替换最终输出的幅值分布缺少大环。# 统计原始时程的相邻采样差值分布 import numpy as np series np.array(series) diffs np.abs(np.diff(series)) print(差值分位数:, np.percentile(diffs, [50, 90, 99]))把相邻采样差值的中位数和 99% 分位数打出来tol取在 50% 分位数到 99% 分位数之间比较合理。如果 99% 分位数和最大值相差一个数量级说明数据里有少量大幅跳变此时tol取中位数的两倍足够如果 99% 分位数已经接近最大值说明信号本身信噪比很差需要回到采集端解决单纯加大tol会连真实载荷一起丢掉。4.3 残波循环闭合滞回环与半循环的统计口径四点法主循环结束后剩余的残波是另一个容易被忽略的参数点。残波本身没有闭合按照 Rainflow 的原始定义它不应该被记成一个完整的循环。但许多商用软件会在残波的首尾复制一段再跑一次把能闭合的段全部抽出来剩下实在闭合不了的段按半循环计入损伤。这里要区分两种口径如果按闭合循环计残波只能算半循环损伤贡献按一半计如果做残波重排把残波首尾相接后再跑一次四点法得到的闭合循环可以按完整循环计入。经验法则是残波点数占原始转折点数的比例超过 20% 时必须做残波重排否则会丢掉相当一部分损伤贡献低于 10% 时直接按半循环处理即可重排带来的增量很小。残波点数占比处理方式循环计数影响 10%半循环计入损伤可忽略10%~20%半循环计入并标注保守性小 20%残波首尾相接重排一次显著“首尾相接”不是简单地把残波序列复制一遍。常见的做法是去掉残波的最后一个点再把残波接到自己后面形成“原残波 原残波去掉末点”的组合序列。这么做是为了避免首尾相接处出现连续相等等值点影响四点法的极值判断。5. 用正弦叠加噪声自检四点法验证、残波重排与结果修正5.1 自检序列正弦波叠加噪声后的循环个数期望写完四点法不能直接拿真实载荷谱去跑因为真实数据你根本不知道正确答案。先构造一个理论结果确定的序列来验证实现一个完整的正弦周期 0→A→0→-A→0理论上包含一个幅值 A 的整循环和两个半循环两个完整周期应该得到两个整循环。如果主循环输出和理论值不一致优先检查极值限定条件和回退逻辑。import math t [i * 0.01 for i in range(201)] # 两个完整周期200 个采样步 sig [5.0 * math.sin(2 * math.pi * x) for x in t] cycles, residual rainflow_four_point(to_extrema(sig, tol1e-9)) semi_cycles len(residual) // 2 print(闭合循环数:, len(cycles), 残波半循环数:, semi_cycles)正弦信号经过峰谷压缩后只剩首、峰值、谷值、末点这几个转折点四点法会把第一段 5→-5 抽成循环剩余残波再首尾相接后又能抽出一个循环最终结果通常是一个完整循环加两个半循环。如果你跑出来的循环数低于这个预期检查to_extrema是否把首尾点都保留住了如果循环数偏高说明主循环把残波里的半循环也当成了整循环。5.2 残波重排一次法把未闭合的尾巴桥接起来验证残波处理时可以用一个最简单的手算序列0→10→0→10→0。这个序列在直观上是“一个完整往返循环外加一次半程爬坡”四点法主循环会先抽出 10→0 这一段循环记录为均值 5、幅值 10残波剩下 [0, 10, 0]。此时把残波去掉末点再拼一次得到 [0, 10, 0, 10]继续运行四点法又可以得到一个循环。def rerun_residual(residual): 残波首尾相接去掉最后一个点后重跑一次四点法。 if len(residual) 4: return residual new_seq residual residual[:-1] extra_cycles, _ rainflow_four_point(new_seq) return extra_cycles重跑后得到的循环需要合并进主循环结果里。注意重跑出来的循环幅值分布和主循环不一定连续因为残波本身是主循环抽取后留下的“骨头”它的长度往往比较短重跑时四点法的判定窗口也会少一些。合并时建议单独列出“主循环”和“残波重排循环”两组记录而不是直接混在一起这样后期编制载荷谱时能清楚看到主计数和残波桥接各自的贡献。验证的最终标准是看总循环数与理论预期的一致性原始时程里整循环的个数加上残波重排后新抽出的整循环个数应等于原始曲线实际包含的闭合滞回环数量。对 0→10→0→10→0 这个序列正确答案是两个整循环对应两次完整的 0→10→0 往返如果你得到的是三个循环那说明残波首尾相接时没有去掉末点把首尾重复点算成了额外的拐点。本文还有配套的精品资源点击获取