企业官网建设流程全解析

简介：面向机械、材料科学与航空航天等领域的疲劳分析需求，这份资源通过MATLAB实现四点法雨流计数，用于从非规则的应力或应变时间序列中提取代表性循环载荷数据，为后续疲劳寿命评估奠定基础。压缩包内仅含1个m文件，整体约1KB，代码精简且可直接运行，适合工程技术人员、科研人员及学习疲劳分析的学生参考和使用；MATLAB环境也让算法调试与结果可视化更为便利。目前已有1742人学习下载，足见其在结构耐久性分析中的实用价值。脚本覆盖极值点识别、相邻循环匹配、半径与中心点计算、归一化处理、重复循环合并及按半径排序等完整流程，可帮助读者快速掌握四点法的实现细节，并将其嵌入自定义数据处理流程，高效完成载荷谱转换与寿命估算。

1. 雨流计数为什么绕不开四点法

拿到一段实测的桥梁应变或风机叶片载荷谱，直接按峰谷值统计循环会得到一堆互相嵌套的假循环：小幅波动被记成独立载荷，大幅加载又被拆分得七零八落，寿命估算结果保守得没有参考价值。雨流计数把时间序列重写成一个个闭合迟滞环，而四点法用四个连续极值点判断哪些峰谷能构成完整循环，不需要预先设定窗口幅值，非常适合处理机载记录的长序列。这份 Rainflow.m 正是用四点法把应力谱转成“半径-中心点”循环对，为后续 Miner 损伤计算提供直接输入。适合做结构疲劳、机械件寿命预测的工程师，也适合需要把雨流算法移植到嵌入式或在线监测平台的人。

2. 四点法雨流计数的循环判据与归一化逻辑

2.1 极值序列：先去掉所有非转折点

四点法的输入不应该直接是原始时间序列，而是一个严格的峰谷交替极值序列。原始数据中相邻采样点可能连续上升或连续下降，这些单调段不产生循环候选，如果全送入算法，会产生大量无意义的“重叠”判断。常见的做法是先做差分：对序列 signal 计算一阶差分 d，当相邻差分乘积 d(i)*d(i+1)<0，说明方向在 i+1 处翻转，该点就是一个局部极大值或极小值。

% 提取转折点，得到峰谷交替序列 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）。半径等于峰谷差的一半，代表循环的应力幅度；中心点等于峰谷平均值，代表平均应力水平。在相同的幅值下，平均应力越高，疲劳损伤通常越大，所以这两个参数必须同时保留。

归一化的常见做法是把半径除以全局最大半径，使所有循环的半径落在 0 到 1 之间，中心点不归一化，保留原始单位。为什么中心点不归一化？因为平均应力的绝对数值和材料特性、工况零点有关，归一化后反而丢失了物理意义。Rainflow.m中归一化基准默认取全部循环半径的最大值，这样输出的第一列始终在 0~1 之间，后续做不同载荷谱的分布叠加时，不需要再处理量纲差异。

归一化后的循环对还需要进行去重和排序，这部分放在第 4 章，但要注意：归一化得到的是浮点数，不能直接作为分组键，必须先离散化比如四舍五入到固定分辨率，否则每个循环都可能是唯一值，根本合并不了。

3. 基于 MATLAB 的 Rainflow.m 实现与核心代码拆解

3.1 函数签名与输出矩阵设计

Rainflow.m 应设计成一个独立函数，而不是脚本，这样能在不同载荷数据上反复调用。函数输入为原始时间序列load，可选参数norm_base是归一化基准半径；输出为一个 N×3 的矩阵cycles，每一行对应一个合并后的循环类型，三列分别是归一化半径、中心点、出现次数。这个输出格式非常紧凑，方便直接喂给 Miner 损伤计算或绘制载荷谱散点图。

设计上最重要的点是把“提取循环”和“循环统计”分开。提取循环过程中只记录半径和中心点，不即时合并，这是因为真实载荷谱中大量循环的半径和中心点非常接近但并非完全相等，如果边提取边合并，阈值设置会直接影响循环提取结果。先记录再统一合并，允许你事后用不同分辨率重做统计，而不需要重新跑算法。

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(k+1); e3 = ext(k+2); e4 = ext(k+3); % 判断中间两点构成的小区间是否被两端大区间完全包含 if min(e2,e3) >= min(e1,e4) && max(e2,e3) <= max(e1,e4) amp(end+1,1) = (max(e2,e3) - min(e2,e3)) / 2; center(end+1,1) = (e2 + e3) / 2; % 提取后删除中间两个点，并回退一步重新检查 ext(k+1:k+2) = []; 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); end

3.3 主循环中的回退逻辑与判据参数

主循环中变量ext(k+1:k+2) = []是四点法最关键的步骤。当 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 个桶；中心点如果原始量级是几百 MPa，0.001 的分辨率过于精细，几乎不会合并任何循环。所以实际使用时要根据中心点量级调大res，比如设成 0.01 或 0.1。更稳妥的做法是单独给半径和中心点分别设置分辨率，例如半径用 1e-3，中心点用 1，避免一个参数把另一个参数的合并效果抵消。

max(amp)作为归一化基准有一个隐含假设：最大半径的循环一定存在且包含在提取结果中。对极端载荷谱，如果最大半径出现在残差半循环中，max(amp)会偏小，导致所有归一化半径偏大。遇到这种情况，可以显式传入一个物理上更合理的基准，例如材料屈服幅值或设计载荷的允许幅值。这样得到的归一化半径才能在不同材料、不同工况间横向比较，而不是单纯依赖数据极值。

4. 极值点预处理、重复循环合并与输出排序的工程处理

4.1 等值平台和单点毛刺的过滤

真实载荷数据里，传感器在长时间保持同一数值时会产生平台，比如停车等待或恒速运行。平台的差分值为 0，diff为 0 时d(i)*d(i+1)也为 0，既不会判为转折点，也不会对极值序列产生贡献，但平台两端会出现两个方向相反的转折，导致算法在平台边界提取出两个紧挨着的极值点，形成幅值接近 0 的假循环。

处理方式是在进入雨流计数前先压缩连续等值点：

% 压缩连续相等值，只保留每个平台的首个点 change_idx = find(diff(signal) ~= 0); compressed = signal([1; change_idx + 1]); if compressed(end) ~= signal(end) compressed(end+1) = 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 重复循环合并的两类策略

原始雨流计数得到的半径和中心点几乎不会完全相等，因为数值浮点误差和测量噪声都会让同一工况的循环出现微小差异。重复循环合并直接影响后续损伤计算的分组精度。常用的策略有三种，实际中按资源场景选择。

第 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(j+1,1) - sorted_rows(i,1)) < tol_r && ... abs(sorted_rows(j+1,2) - sorted_rows(i,2)) < tol_c j = j + 1; end block = sorted_rows(i:j, :); merged(end+1, :) = [mean(block(:,1)), mean(block(:,2)), sum(block(:,3))]; i = j + 1; end

这里的逻辑是先把循环按半径升序排列，然后依次把半径差小于tol_r、中心点差小于tol_c的相邻行合并成一组，组内半径和中心点取均值，次数取和。容差的选择取决于数据量：归一化半径 0.002 对应幅值 0.2% 的变化，一般载荷谱完全够用；中心点容差需要看载荷单位，如果应力单位是 MPa，0.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 曲线的斜率。

% 材料参数示例：k=5，N0=1e7 对应 S0=200 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 自带脚本就能完成，很适合在项目早期筛选典型工况。

本文还有配套的精品资源，点击获取