简介面向机械、材料科学与航空航天等领域的疲劳分析需求这份资源通过MATLAB实现四点法雨流计数用于从非规则的应力或应变时间序列中提取代表性循环载荷数据为后续疲劳寿命评估奠定基础。压缩包内仅含1个m文件整体约1KB代码精简且可直接运行适合工程技术人员、科研人员及学习疲劳分析的学生参考和使用MATLAB环境也让算法调试与结果可视化更为便利。目前已有1742人学习下载足见其在结构耐久性分析中的实用价值。脚本覆盖极值点识别、相邻循环匹配、半径与中心点计算、归一化处理、重复循环合并及按半径排序等完整流程可帮助读者快速掌握四点法的实现细节并将其嵌入自定义数据处理流程高效完成载荷谱转换与寿命估算。1. 雨流计数为什么绕不开四点法拿到一段实测的桥梁应变或风机叶片载荷谱直接按峰谷值统计循环会得到一堆互相嵌套的假循环小幅波动被记成独立载荷大幅加载又被拆分得七零八落寿命估算结果保守得没有参考价值。雨流计数把时间序列重写成一个个闭合迟滞环而四点法用四个连续极值点判断哪些峰谷能构成完整循环不需要预先设定窗口幅值非常适合处理机载记录的长序列。这份 Rainflow.m 正是用四点法把应力谱转成“半径-中心点”循环对为后续 Miner 损伤计算提供直接输入。适合做结构疲劳、机械件寿命预测的工程师也适合需要把雨流算法移植到嵌入式或在线监测平台的人。2. 四点法雨流计数的循环判据与归一化逻辑2.1 极值序列先去掉所有非转折点四点法的输入不应该直接是原始时间序列而是一个严格的峰谷交替极值序列。原始数据中相邻采样点可能连续上升或连续下降这些单调段不产生循环候选如果全送入算法会产生大量无意义的“重叠”判断。常见的做法是先做差分对序列 signal 计算一阶差分 d当相邻差分乘积 d(i)*d(i1)0说明方向在 i1 处翻转该点就是一个局部极大值或极小值。% 提取转折点得到峰谷交替序列 d diff(signal); turn_idx find(d(1:end-1) .* d(2:end) 0) 1; ext signal([1; turn_idx(:); length(signal)]); % 去掉首尾与相邻极值重复的点 keep [true; diff(ext(:)) ~ 0]; ext ext(keep);这段代码里d(1:end-1).*d(2:end)0是核心判断相邻两个差分方向相反则它们中间的那个原始点必然是转折点。turn_idx保存的是原序列索引加 1 是因为 diff 后索引从 2 开始对应原序列。首尾点虽然不一定真的是局部极值但雨流计数通常把端点也视为极值候选这样可以完整覆盖第一个和最后一个半循环。keep逻辑向量用于删除首尾与第二个点相等的情况保证后续四点判断时不会出现连续两个相同点。提取后的ext序列是后续所有工作的基础。注意这里的提取是“等值点友好”的但如果原始信号里存在持续平台比如传感器饱和需要先做压缩这一步放到后面专门讲。对大多数正常载荷谱这样提取出来的极值个数会远小于原始数据长度雨流计数的计算量因此大大降低。2.2 四点判据是循环封闭性的最小判据得到峰谷交替序列后要判断哪些峰谷能够组成一个完整循环。四点法取连续的四个极值点 e1、e2、e3、e4判断由中间两点构成的小区间是否完全被两端点构成的大区间包含。包含的数学条件是min(e2,e3) min(e1,e4)且max(e2,e3) max(e1,e4)这个条件成立时说明 e2 与 e3 之间发生了一次完整的加载-卸载过程且这个过程的幅值没有超出 e1 与 e4 的包围范围。从材料应力-应变回线角度看e2-e3 对应一个可独立闭合的迟滞环可以提取出来不破坏更大循环的完整性。如果条件不成立说明中间两点的波动只是整体趋势的一部分不能单独计作循环需要向后滑动一个点继续判断。为什么不用三点三点只能判断相邻两个极值点是否构成拐点无法区分“小循环嵌套在大循环内”和“整体上升中的锯齿”。雨流计数的核心目标就是把嵌套的小循环逐层剥掉只保留封闭回线四点判据是能实现这一目标的最短窗口。标准雨流算法中的三峰谷法和四点法在多数工况下结果一致但四点法在流式数据输入时更易于实现因为它只需要维护一个长度为 4 的滑动窗口不需要回溯整段历史。2.3 循环半径、中心点与归一化基准每个被提取的循环都有一对峰谷值 high 和 low其中 high 是循环中的极大值low 是极小值。疲劳分析通常用两个量描述循环半径cyclic half-amplitude和中心点middle value。半径等于峰谷差的一半代表循环的应力幅度中心点等于峰谷平均值代表平均应力水平。在相同的幅值下平均应力越高疲劳损伤通常越大所以这两个参数必须同时保留。参数公式物理含义半径 R(high - low) / 2循环幅值的一半决定损伤权重中心点 M(high low) / 2平均应力影响平均应力修正归一化半径 rR / R_max无量纲幅值便于跨工况比较归一化的常见做法是把半径除以全局最大半径使所有循环的半径落在 0 到 1 之间中心点不归一化保留原始单位。为什么中心点不归一化因为平均应力的绝对数值和材料特性、工况零点有关归一化后反而丢失了物理意义。Rainflow.m中归一化基准默认取全部循环半径的最大值这样输出的第一列始终在 0~1 之间后续做不同载荷谱的分布叠加时不需要再处理量纲差异。归一化后的循环对还需要进行去重和排序这部分放在第 4 章但要注意归一化得到的是浮点数不能直接作为分组键必须先离散化比如四舍五入到固定分辨率否则每个循环都可能是唯一值根本合并不了。3. 基于 MATLAB 的 Rainflow.m 实现与核心代码拆解3.1 函数签名与输出矩阵设计Rainflow.m 应设计成一个独立函数而不是脚本这样能在不同载荷数据上反复调用。函数输入为原始时间序列load可选参数norm_base是归一化基准半径输出为一个 N×3 的矩阵cycles每一行对应一个合并后的循环类型三列分别是归一化半径、中心点、出现次数。这个输出格式非常紧凑方便直接喂给 Miner 损伤计算或绘制载荷谱散点图。输入/输出变量说明输入load应力或应变时间序列必须为列向量输入norm_base归一化基准半径留空时自动取最大半径输出cycles(:,1)归一化后的循环半径输出cycles(:,2)循环中心点保持原始单位输出cycles(:,3)相同半径和中心点组合出现的次数输出norm_base归一化基准半径供还原原始幅值使用设计上最重要的点是把“提取循环”和“循环统计”分开。提取循环过程中只记录半径和中心点不即时合并这是因为真实载荷谱中大量循环的半径和中心点非常接近但并非完全相等如果边提取边合并阈值设置会直接影响循环提取结果。先记录再统一合并允许你事后用不同分辨率重做统计而不需要重新跑算法。3.2 极值提取与四点主循环的完整代码下面是对应于第 2 章判据的可运行实现。它保留了四点法的核心回退逻辑同时做了一份简洁的归一化和分组输出。function [cycles, norm_base] rainflow_point4(load, norm_base) % RAINFLOW_POINT4 四点法雨流计数 % 输入 load: 应力或应变时间序列列向量 % 输入 norm_base: 可选归一化基准半径默认取全局最大值 % 输出 cycles: Nx3 矩阵[归一化半径, 中心点, 出现次数] % 输出 norm_base: 实际使用的归一化基准半径 % 1. 提取峰谷交替极值点 d diff(load); turn_idx find(d(1:end-1) .* d(2:end) 0) 1; ext load([1; turn_idx(:); length(load)]); keep [true; diff(ext(:)) ~ 0]; ext ext(keep); % 2. 四点法主循环 amp []; % 半径未归一化 center []; % 中心点 k 1; while length(ext) - k 3 e1 ext(k); e2 ext(k1); e3 ext(k2); e4 ext(k3); % 判断中间两点构成的小区间是否被两端大区间完全包含 if min(e2,e3) min(e1,e4) max(e2,e3) max(e1,e4) amp(end1,1) (max(e2,e3) - min(e2,e3)) / 2; center(end1,1) (e2 e3) / 2; % 提取后删除中间两个点并回退一步重新检查 ext(k1:k2) []; k max(1, k-1); else k k 1; end end % 3. 半径归一化 if nargin 2 || isempty(norm_base) norm_base max(amp); end if norm_base 0 r amp / norm_base; else r amp; end % 4. 按半径和中心点分桶去重 res 1e-3; rq round(r / res) * res; cq round(center / res) * res; [~, ia, ic] unique([rq, cq], rows); count accumarray(ic, 1); cycles sortrows([rq(ia), cq(ia), count], 1); end3.3 主循环中的回退逻辑与判据参数主循环中变量ext(k1:k2) []是四点法最关键的步骤。当 e2-e3 被判定为一个有效循环后这两个点就不能再参与后续循环的配对直接删除。删除后原来的 e1 与后面新露出的点组成新的四点窗口之前因为 e2-e3 的存在而被遮挡的循环结构可能显现出来所以索引必须回退。这里用k max(1, k-1)确保 k 不能小于 1否则下一次访问 ext(0) 会报错。如果 k 本来等于 1回退后仍然从第一个点开始重新检查。判据中的min(e2,e3) min(e1,e4)与max(e2,e3) max(e1,e4)是完整包含关系。实测时如果数据经过滤波直接用严格浮点比较通常没问题但如果是嵌入式数据或含有噪声建议在最外层先做一次轻微平滑否则极值点提取会捕获大量由噪声产生的虚假峰谷四点判据会把它们当真实循环提取出来。这点在第 4 章会给出处理参数。还应注意amp和center在 while 循环中按行追加预分配可以忽略。循环结束后ext中剩下的点无法再构成四点窗口这些残差点是趋势项或半循环严格来说也应该计入疲劳损伤但这里的实现按完整循环处理残差半循环的权重处理在第 4 章末尾说明。3.4 归一化基准和去重精度对后续分析的影响函数中的res 1e-3是分组分辨率半径和中心点都四舍五入到小数点后三位。归一化半径本身在 0~1 之间0.001 的分辨率等价于把幅值分布切成 1000 个桶中心点如果原始量级是几百 MPa0.001 的分辨率过于精细几乎不会合并任何循环。所以实际使用时要根据中心点量级调大res比如设成 0.01 或 0.1。更稳妥的做法是单独给半径和中心点分别设置分辨率例如半径用 1e-3中心点用 1避免一个参数把另一个参数的合并效果抵消。max(amp)作为归一化基准有一个隐含假设最大半径的循环一定存在且包含在提取结果中。对极端载荷谱如果最大半径出现在残差半循环中max(amp)会偏小导致所有归一化半径偏大。遇到这种情况可以显式传入一个物理上更合理的基准例如材料屈服幅值或设计载荷的允许幅值。这样得到的归一化半径才能在不同材料、不同工况间横向比较而不是单纯依赖数据极值。4. 极值点预处理、重复循环合并与输出排序的工程处理4.1 等值平台和单点毛刺的过滤真实载荷数据里传感器在长时间保持同一数值时会产生平台比如停车等待或恒速运行。平台的差分值为 0diff为 0 时d(i)*d(i1)也为 0既不会判为转折点也不会对极值序列产生贡献但平台两端会出现两个方向相反的转折导致算法在平台边界提取出两个紧挨着的极值点形成幅值接近 0 的假循环。处理方式是在进入雨流计数前先压缩连续等值点% 压缩连续相等值只保留每个平台的首个点 change_idx find(diff(signal) ~ 0); compressed signal([1; change_idx 1]); if compressed(end) ~ signal(end) compressed(end1) signal(end); end如果change_idx为空说明整个序列是常数此时不存在任何循环直接返回空矩阵即可。压缩后的序列仍保留首尾值但中间平台段不再贡献多余极值。要注意的是压缩处理必须在差分提取之前否则平台中间的等值点仍会被当成普通样本进入ext。单点毛刺指某一个采样点异常偏离相邻点比如应变片受电磁干扰。毛刺会在序列中制造一正一负两个转折雨流计数会把它识别为一个极小循环。一般用与采样频率匹配的中值滤波处理MATLAB 中medfilt1(signal, 3)可以去除单点脉冲但会让真实的快速峰值变得圆滑所以只对明显噪声段使用。更保守的做法是在提取极值后加入最小幅值阈值如果某相邻峰谷差小于传感器分辨率或材料疲劳极限对应幅值则直接丢弃该循环。这里给出一个极值后筛选示例min_amp 0.1; % 根据实际载荷单位设定 valid (max(ext(1:end-1), ext(2:end)) - min(ext(1:end-1), ext(2:end))) / 2 min_amp; ext ext([true; valid(:)]);这是极值相邻配对筛选仅用于剔除明显噪声循环。4.2 重复循环合并的两类策略原始雨流计数得到的半径和中心点几乎不会完全相等因为数值浮点误差和测量噪声都会让同一工况的循环出现微小差异。重复循环合并直接影响后续损伤计算的分组精度。常用的策略有三种实际中按资源场景选择。策略实现思路优点缺点适用场景取整分桶四舍五入到固定分辨率后 unique简单、快边界截断可能把相近循环拆开快速预筛容差聚类按半径和中心点容差逐个合并保留分布连续性复杂度高需要设定容差高精度损伤评估分位数分桶按分位数划分区间分组数量可控区间边界依赖样本分布多工况统计第 3 章代码中的res取整法属于第一种优点是代码最少但要注意取整会让处于两个桶边界的循环被强行分开。如果在疲劳寿命评估中需要更细致的循环分布推荐用容差聚类。下面是一个兼顾半径和中心点的合并实现% 按半径排序后合并半径和中心点均落在容差内的相邻循环 tol_r 0.002; % 半径容差归一化单位 tol_c 0.5; % 中心点容差按原始单位设定 sorted_rows sortrows(cycles, 1); merged []; i 1; while i size(sorted_rows, 1) j i; while j 1 size(sorted_rows, 1) ... abs(sorted_rows(j1,1) - sorted_rows(i,1)) tol_r ... abs(sorted_rows(j1,2) - sorted_rows(i,2)) tol_c j j 1; end block sorted_rows(i:j, :); merged(end1, :) [mean(block(:,1)), mean(block(:,2)), sum(block(:,3))]; i j 1; end这里的逻辑是先把循环按半径升序排列然后依次把半径差小于tol_r、中心点差小于tol_c的相邻行合并成一组组内半径和中心点取均值次数取和。容差的选择取决于数据量归一化半径 0.002 对应幅值 0.2% 的变化一般载荷谱完全够用中心点容差需要看载荷单位如果应力单位是 MPa0.5 MPa 的合并精度会对平均应力产生过于精细的分组导致每个组合仍有大量零散次数实际使用时根据 S-N 曲线对平均应力的敏感度调整。4.3 排序输出与后续接口的对接合并后的cycles通常按半径排序。排序的主要作用是让后续损伤计算可以从最大循环开始累加同时便于绘制幅值累积频次曲线。MATLAB 中sortrows(merged, 1)即可按第一列升序排列。如果想按降序使用sortrows(-merged(:,1))并手动拼接或者用sort(..., descend)。排序后矩阵的每一行已经是可以直接交给accumarray或histogram的统计量。如果下一步要导入 Python 或数据库建议将归一化半径恢复为原始半径后导出original_radius merged(:,1) * norm_base; export_matrix [original_radius, merged(:,2), merged(:,3)]; writematrix(export_matrix, rainflow_result.csv);这里的norm_base必须是主循环中实际使用的归一化基准也就是调用rainflow_point4时获得的第二个输出。导出时保留归一化半径和原始半径两列更稳妥因为后面做不同载荷谱对比时归一化半径用于横向比较原始半径用于损伤计算。4.4 残差半循环与标准雨流实现的差异四点法主循环结束后ext中剩余的点无法再满足四点判据这些点构成的开放回线就是残差半循环。标准雨流算法会把剩余相邻极值点对当作半循环每个半循环在 Miner 损伤中只计一半权重。这份Rainflow.m代码直接将剩余点丢弃这在长载荷谱中影响不大因为绝大部分循环已在主循环提取残差通常只占总累计幅值的几个百分点。但如果研究对象是低周疲劳或载荷循环次数少残差占比会显著上升必须补上半循环处理。补半循环的常见做法是对剩余的ext序列分别从第 1 个点和第 2 个点开始相邻配对计算半径和中心点并将出现次数按 0.5 记录。在损伤计算中半循环的寿命 N 按完整循环公式计算但累计损伤时乘以 0.5。注意半循环之间也可能存在嵌套关系简单相邻配对会损失精度工程上推荐使用三点法或四点法结合残差回线组合算法这里不展开。5. 从雨流矩阵到 Miner 累积损伤评估的快速验证技巧5.1 用循环矩阵快速计算线性累积损伤雨流循环提取出来以后最常见的落地方式就是 Miner 线性累积损伤。对每个归一化半径 r先还原成原始半径 R_orig r * R0其中 R0 是归一化基准。循环幅值峰谷差S 2 * R_orig。Basquin 公式给出一个材料在某应力幅 S 下的寿命 N N0 * (S0 / S)^k其中 S0 是参考幅值N0 是对应寿命k 是材料 S-N 曲线的斜率。% 材料参数示例k5N01e7 对应 S0200 MPa k 5; N0 1e7; S0 200; R0 norm_base; % 归一化时使用的基准半径 damage 0; for i 1:size(cycles, 1) r cycles(i, 1); count cycles(i, 3); S_amplitude 2 * r * R0; % 循环峰谷差 N_i N0 * (S0 / S_amplitude)^k; damage damage count / N_i; end fprintf(累积损伤 D %.4f\n, damage);S0和k必须从材料手册或试验数据中取这里只是示例。雨流计数和 Miner 假设是配套用的雨流把不规则载荷拆成若干个等幅循环Miner 再把每个循环的损伤线性叠加忽略加载顺序。如果damage接近或超过 1说明结构在该载荷谱下大概率发生疲劳破坏。5.2 用幅值-均值散点图验证循环分布是否异常在跑完整段计算之前先用散点图看循环分布是一种高效的验证手段。将中心点作为横轴、原始半径作为纵轴并用次数作为点大小scatter(cycles(:,2), cycles(:,1) * R0, 10 cycles(:,3), cycles(:,3), filled); xlabel(中心点 / MPa); ylabel(半径 / MPa); colorbar;正常载荷谱的雨流结果应该呈现沿中心点方向连续分布的条带且在低幅值区域点密度最高。如果发现某个孤立点上出现异常大半径或中心点跳变多半是原始信号中的毛刺未被完全过滤或者四点判据把趋势项当成循环。此时可以回到第 4 章的中值滤波和极值筛选步骤重新检查最小幅值阈值的选取。另一个可验证的量是循环总次数雨流提取的循环次数不会超过原始数据点数的一半如果超过说明极值提取有重复点。5.3 快速对比多工况载荷谱的入口当需要对多段实测数据分别做雨流计数时把归一化半径的分布矩阵做核密度估计可以快速看出哪个工况更“残酷”。例如用ksdensity(cycles(:,1))得到归一化半径概率密度分布右尾越厚说明大循环占比越高。中心点均值变化则反映平均应力偏移配合 Goodman 修正能判断是否需要调整安全系数。这个方法不需要额外的疲劳软件MATLAB 自带脚本就能完成很适合在项目早期筛选典型工况。本文还有配套的精品资源点击获取
