企业官网建设流程全解析

简介：在机械工程与材料疲劳分析中，雨流计数法用于将复杂的应力应变历程转化为可统计的完整循环。这份基于MATLAB实现的资源包，面向需要处理载荷谱的工程师、科研人员和相关专业学生，解决了从数据预处理、极值识别、循环匹配到雨流矩阵生成的全流程需求。压缩包共包含15个文件，其中7个m脚本与函数覆盖算法主体与示例演示，辅以C源码、dll动态库、HTML说明文档及可视化图像，整体仅112KB，精简且便于理解和二次开发。目前已有1423人学习下载，是疲劳寿命预测场景中较为实用的参考工具。借助计数函数、特征提取脚本和多个demo，使用者可直接导入自己的应力应变数据，快速获得等效应力循环与雨流矩阵，从而支撑结构疲劳损伤评估与寿命预测。

1. 用Matlab跑雨流计数，为什么说“关键在参数和残谱”

打开rainflow.zip，里面通常只有一个rainflow.m和几个辅助文件，很多人第一反应就是把载荷数组塞进去，然后对着一串输出矩阵发懵。雨流计数（Rainflow counting）在疲劳载荷谱分析里的地位，相当于傅里叶变换在信号处理里面：它把一段不规则的时间历程拆成完整的应力循环，每个循环用幅值和均值表征，再映射到S-N曲线做Miner损伤累积。但实际跑通后你会发现，结果行数比肉眼数的载波周期多，或者全循环和半循环混在一起不知道该怎么用。这往往不是算法代码坏了，而是你没有先校准两件事：输入数据的峰谷抽取规则，以及残谱到底按0.5还是按1处理。这篇文章给出一个完整、可直接运行的Matlab实现，从峰谷抽取、四点法主循环到参数设置和损伤计算全部对齐，读完你也能判断网上下载的rainflow.zip是否靠谱。

2. 雨流计数原理与Matlab函数实现：从峰谷提取到四点法出循环

2.1 雨流计数法在疲劳分析中的位置

疲劳寿命预估最怕的不是平均应力高，而是随机载荷里藏着各种大小不等的闭环。直接数峰谷或者用穿越计数，会把内部的非完整循环漏掉或者重复计入。雨流计数把载荷历程当作屋檐，雨滴从极值点流下，形成一个个完整的滞后环。ASTM E1049-85(2011)是几乎所有雨流实现的参考基准，它规定输出至少包含循环的幅值Sa、均值Sm和循环次数。对于全循环，计数为1；对于残谱中无法闭合的半循环，计数为0.5。后面做Miner累积时，只需要这三列数据，所以实现雨流计数的核心其实就是把峰谷序列按“包含关系”反复裁剪，直到剩下不能配对的残余点。

2.2 先做峰谷抽取：过滤掉无效斜率点

很多原始载荷信号是直线段连接起来的，直接拿相邻点去做四点法会留下大量无意义的中间点。必须先把信号压缩成峰谷交替的极值序列。常见做法是用diff做一阶差分，再对差分结果求二阶差分，斜率变号的位置就是极值点。

function pv = extractPeakValley(signal) % extractPeakValley 从原始载荷时间历程中提取峰谷序列 % signal: 一维列向量，载荷时间历程 % pv : 只保留极值点和两端的序列，相邻点符号必然交替 signal = signal(:); delta = diff(signal); % 相邻点差值 if all(delta == 0) pv = signal([1 end]); return; end signChange = find(diff(delta) ~= 0) + 1; % 斜率变号的位置 idx = unique([1; signChange; length(signal)]); pv = signal(idx); end

这段代码里，diff(delta) ~= 0表示当前点的前一段差值和后一段差值不同号，也就是局部极值。unique用来去重，防止平台段被重复选中。不过要特别注意，如果原始信号存在宽度大于1的平顶，diff会出现连续的0，这时候diff(delta) ~= 0会漏掉平顶终点。工业数据里这种平台很常见，建议在调用前先对连续重复值做一次去重，也就是只保留每个平台的首尾两点。

2.3 四点法提取循环：Matlab核心算法实现

峰谷序列准备好之后，主循环就可以开始了。这里采用A. Nieslony经典四点法思路，用一个缓冲数组R模拟雨流过程：每读入一个新峰谷点，就检查R末尾四个点是否满足“中间两个点形成的行程被前后两段行程包住”，如果满足则提取一个完整循环，然后删除中间两个点，让修改后的序列继续参与检查。

function [cycles, residual] = rainflowCounting(pv) % rainflowCounting 四点法雨流计数 % pv : extractPeakValley() 输出的峰谷序列 % cycles : Mx3矩阵，依次为 [幅值 Sa, 均值 Sm, 循环次数] % residual : 未配对成完整循环的剩余峰谷点 pv = pv(:); n = length(pv); if n < 4 cycles = []; residual = pv; return; end R = pv(1:2); % 缓冲区，先放前两个点 cycles = zeros(floor(n/2) + 64, 3); cIdx = 0; i = 3; while i <= n R(end+1) = pv(i); % 新点送入缓冲 while length(R) >= 4 y1 = R(end-3); y2 = R(end-2); y3 = R(end-1); y4 = R(end); d1 = y2 - y1; d2 = y3 - y2; d3 = y4 - y3; if d1*d2 < 0 && d2*d3 < 0 && ... abs(d2) <= abs(d1) && abs(d2) <= abs(d3) cIdx = cIdx + 1; cycles(cIdx,:) = [abs(d2)/2, (y2+y3)/2, 1]; R(end-2:end-1) = []; % 删除y2,y3，让剩余序列继续比较 else break; end end i = i + 1; end % 残余点按相邻点对拆成半循环 for k = 1:length(R)-1 cIdx = cIdx + 1; cycles(cIdx,:) = [abs(R(k+1)-R(k))/2, (R(k)+R(k+1))/2, 0.5]; end cycles(cIdx+1:end,:) = []; residual = R; end

代码里d1*d2<0用来保证y2是极值点，d2*d3<0保证y3也是极值点，两者方向相反。abs(d2) <= abs(d1)和abs(d2) <= abs(d3)表示y2到y3这段行程不会超过相邻的前段和后段，因此在雨流轨迹上它先被“切”出来。循环提取后删除中间两个点，相当于这段已经闭合，剩下的重新连接继续检查。最后缓冲区里剩下的点无法再形成闭合循环，按半循环处理，计数为0.5。这里输出的第一列是幅值abs(d2)/2，不是峰谷差，注意别和某些zip包里的“range”混淆。

2.4 输出结果格式与ASTM E1049的对应关系

无论是自己写的还是从rainflow.zip里拿到的代码，第一步都应该确认输出矩阵的列顺序。以下是我建议的统一约定：

ASTM E1049并不强制规定列顺序，所以很多zip包会输出[range, mean, count]甚至[start_index, end_index, range, mean]。拿到外部代码后，先用一个已知波形跑一遍，逐列对比，再做后续操作，否则Miner损伤算出来会完全对不上。

3. 在Matlab中设置雨流计数参数：数据预处理、门槛值、残谱处理

3.1 输入数据的常见转换：从Excel/CSV导入和时间戳处理

实际工作时，载荷数据往往是csv或xlsx文件，第一列时间戳，第二列应变或应力。读取后第一步是清洗，去掉NaN和重复时间戳。Matlab的readmatrix从R2019a开始可用，之前的版本要退回csvread或xlsread。

data = readmatrix('load_history.csv'); % 或者 .xlsx strain = data(:,2); % 假设第二列是应力或应变 strain = strain(~isnan(strain)); % 删除无效行 % 如果采样率太高，可以按时间抽稀 t = data(:,1); keepIdx = 1:round(length(t)/1e6):length(t); % 保留最多100万点 strain = strain(keepIdx);

keepIdx的步长根据总点数调整，工程上100万点做雨流完全够用。抽稀前最好先对信号做峰值保持或者低通滤波，否则会丢掉真实的极值点。直接每隔N点抽取虽然简单，但可能切掉局部峰，导致后续循环计数偏小。

3.2 幅值门槛和均值离散化：避免小循环淹没统计

载荷信号里总是混着高频噪声，噪声产生的微小循环会让cycles矩阵膨胀几十倍，而对疲劳损伤几乎没有贡献。设置一个幅值门槛是最直接的过滤方式。门槛值可以用整个信号最大最小值的百分比，也可以用绝对应力值。

过滤代码很简单：

amp = cycles(:,1); meanVal = cycles(:,2); cnt = cycles(:,3); valid = amp > ampThreshold; cyclesFiltered = cycles(valid, :);

要注意的是，门槛设太高会把真实的小幅值载荷循环也去掉，低估损伤；设太低则过滤效果差。最稳妥的做法是先绘制幅值分布直方图，看0.5倍噪声以下的峰在哪里，把ampThreshold放在噪声峰和真实循环峰之间的谷位上。

3.3 残谱的两种处理方式：半循环 vs 闭合循环

雨流计数结束后，缓冲数组里剩下的点无法配对成完整循环，这就是残谱。标准做法是把残谱的相邻点对当作半循环，计数为0.5。但这种处理在Miner损伤中会带来一个尴尬：半循环的两端没有闭合，实际损伤可能被低估。另一个常见做法是把残谱首尾相连形成闭合循环，这样所有残余贡献都被归入一个完整加载路径。

如果要把残谱当成一个闭合循环，可以在rainflowCounting返回后手动拼接：

if ~isempty(residual) && length(residual) > 2 delta = residual(end) - residual(1); closedCycle = [abs(delta)/2, (residual(1)+residual(end))/2, 1]; cyclesFiltered(end+1, :) = closedCycle; end

大多数工程规范（比如风力发电机组载荷计算）会把残谱按半循环处理，因为首尾闭合会引入一个实际不存在的大循环，所以除非你的规范明确要求闭合，否则保留半循环计数0.5更安全。

3.4 一个完整的demo脚本模板

把前面的代码串起来，就是一个可以直接运行的测试脚本：

% demo_rainflow.m t = (0:0.01:50)'; signal = 10*sin(2*pi*0.15*t) + 4*sin(2*pi*0.61*t + 1.2); pv = extractPeakValley(signal); [cycles, residual] = rainflowCounting(pv); fprintf('总循环数: %d\n', size(cycles,1)); fprintf('最大幅值: %.3f\n', max(cycles(:,1))); histogram(cycles(:,1), 30); xlabel('幅值'); ylabel('循环数');

脚本先构造一个双正弦叠加信号，提取峰谷后做雨流计数，最终打印循环总数和最大幅值。如果输出循环数是0，说明extractPeakValley返回的序列长度小于4，检查原始数据是否全是常量或者清洗过度。

4. 实际案例：用雨流计数法处理正弦波叠加随机谱，并验证计数结果

4.1 构造测试信号

为了验证函数的正确性，先用一个可控信号：慢波加快波，慢波周期10秒，快波周期2.5秒，叠加后肉眼能数出大约5个慢波周期，但快波会产生附加循环。构造信号代码如下：

t = (0:0.001:50)'; signal = 10*sin(0.2*pi*t) + 3*sin(0.8*pi*t + 0.5) + 0.5*randn(size(t));

加入高斯噪声是为了测试ampThreshold的实际效果。这个信号总点数5万，直接雨流计数会得到大量噪声循环，正好用来观察阈值过滤前后的差异。

4.2 运行计数并绘制雨流参数直方图

调用我们的函数后，用二维直方图看均值-幅值分布：

pv = extractPeakValley(signal); [cycles, residual] = rainflowCounting(pv); figure; histogram2(cycles(:,2), cycles(:,1), 30, 'DisplayStyle','tile'); xlabel('均值 Sm'); ylabel('幅值 Sa'); colorbar;

histogram2的第一个输入是均值，第二个是幅值，每个tile的颜色代表该应力区间内循环数量。从图里能直观看到慢波产生的大幅值循环集中在均值0附近，噪声循环则散落在小幅值区域。此时可以先设置ampThreshold = 0.5，再重新统计过滤后的循环数，通常过滤后循环总数会下降一到两个数量级。

4.3 验证计数结果：循环数和幅值总和为什么对不上？

很多人在这一步会发问：明明信号只有5个慢波周期，为什么输出循环数大于5？答案在于快波循环和噪声循环也都被计数了，另外残谱半循环也会增加计数。为了单独验证算法，用一段无噪声正弦波测试：

testT = (0:0.01:20)'; testSig = 5*sin(2*pi*0.1*testT); % 2个完整周期 pvTest = extractPeakValley(testSig); [cycTest, ~] = rainflowCounting(pvTest); expected = 2; % 两个完整正弦周期产生两个全循环 actual = sum(cycTest(:,3) == 1);

说明如下：

慢波和快波会形成不同幅值的循环，落在直方图的不同区间。如果输出里出现半循环计数0.5，这属于残谱，不是错误。反过来，如果全循环数量大于理论值，通常是因为峰值谷值抽取时把平台段重复计数，检查extractPeakValley对平顶信号的处理即可。

4.4 不同实现之间的差异：rainflow.zip中的旧版和新版不一致怎么办

网上流传的rainflow.zip版本很多，有的输出[range, mean, count]，有的输出[mean, range, start, end]，还有的是C-MEX编译版本。处理方式不是一个个试，而是先构造一个非常简单的闭环序列做基准测试，比如[0 1 0 -1 0]。该序列包含一个从-1到1再到-1的闭合滞后环，理论上应该输出1个全循环，幅值1，均值0。用这段基准去跑你下载的代码，逐列比对哪个输出列与这些值吻合，就能快速定位列顺序和单位。

baseline = [0 1 0 -1 0]'; pvBase = extractPeakValley(baseline); [cycBase, resBase] = rainflowCounting(pvBase); disp(cycBase);

这段脚本输出应该包含一行[1, 0, 1]（幅值1，均值0，计数1）。如果实际输出是[2, 0, 1]，说明单位是峰谷差而不是幅值；如果均值列是0.5，说明列顺序不同。用这样10个点就能排查完外部代码，再放到真实数据上。

5. 进阶应用：把雨流计数结果直接做Miner疲劳损伤累积

5.1 从雨流结果到S-N曲线：Goodman平均应力修正

雨流计数输出的均值对应循环的平均应力。同一个幅值在不同均值下对损伤的影响不同，通常用Goodman直线修正到等效应力幅。设材料极限强度为Su，等效应力幅Sa_eq = Sa / (1 - Sm / Su)。

Sa = cycles(:,1); Sm = cycles(:,2); Su = 800; % 材料抗拉强度，MPa Sa_eq = Sa ./ (1 - Sm ./ Su); Sa_eq(Sm >= Su) = inf; % 超过极限的直接判为失效

注意Sm有可能是负的，也就是压缩平均应力，此时Goodman修正后的Sa_eq会变小，这是合理的，因为压平均应力通常提高疲劳寿命。如果分母出现0或负数，说明该循环已静载破坏，直接标记为inf或从统计中剔除。

5.2 用循环结果矩阵快速计算Miner损伤

得到等效应力幅后，用Basquin公式估算每个循环的寿命N_f = C * Sa_eq^(-m)，其中C和m是材料S-N曲线参数。Miner线性累积损伤为所有循环损伤的总和：

C = 1e12; % 材料常数 m = 3.0; % S-N曲线斜率 N_f = C * Sa_eq.^(-m); D_total = sum(cycles(:,3) ./ N_f); % 半循环的0.5也会计入 if D_total >= 1 fprintf('预测疲劳失效，累积损伤 D = %.3f\n', D_total); else fprintf('安全，累积损伤 D = %.4f\n', D_total); end

这里cycles(:,3)同时包含全循环和半循环，半循环的0.5会直接参与求和，不需要额外标记。计算过程中如果Sa_eq中有0，N_f会变成无穷，损伤为0，对应零幅值循环，可以提前过滤。整个流程从雨流计数到损伤计算就闭环了。

5.3 性能优化与大规模数据的处理技巧

长周期载荷信号动辄几百万点，提取峰谷后点数通常能压缩到原来的1/10到1/100，因此rainflowCounting输入规模并不大。但输出循环矩阵可能仍然很大，尤其是噪声未滤除时。优化思路有三个：第一，优先做峰谷抽取和幅值门槛过滤，把循环数降到十万以内；第二，把幅值和均值离散化到二维直方图，后面损伤计算直接对每个bin累计，而不是逐行循环；第三，如果单个文件实在太长，把信号分段做雨流，分段之间会丢失跨界循环，所以拆段时要在首尾各保留一段重叠区域，重叠长度至少覆盖一个最大周期。

numBins = 50; binEdgesAmp = linspace(0, max(cycles(:,1)), numBins); binEdgesMean = linspace(min(cycles(:,2)), max(cycles(:,2)), numBins); damageMap = histcounts2(cycles(:,2), cycles(:,1), ... binEdgesMean, binEdgesAmp);

这段代码把循环按均值和幅值投到50×50的bin矩阵中，后续Miner损伤计算就可以对每个bin里的循环数统一处理，显著加速。如果下载的rainflow.zip里带C-MEX编译版本，记得先运行mex -setup配置编译器，再编译对应文件，否则Matlab会静默回退到较慢的m文件版本。

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