1. 单通道应力时程的雨流计数基础实现

任务对象与算法约定

本例处理一段等采样间隔的单通道应力时程，采样间隔 Δt = 0.5 s，共 24 个采样点，单位 MPa。雨流计数的输入是离散应力序列，但算法真正操作的是峰谷序列：相邻同向点（连续上升或连续下降中的中间点）不构成转折，必须先压缩掉。本例采用经典的四点法（也称 4 点递归法）实现，其核心判断是：取连续四个转折点 R1、R2、R3、R4，若区间 R2–R3 的幅值不大于 R1–R2 与 R3–R4 的幅值，则 R2–R3 构成一个完整循环，可被提取并从序列中删除，然后回退重新判断。

为便于验证，先明确三个约定：幅值取半幅还是全幅——本例统一用全幅 S = |R3 − R2|，均值 M = (R2 + R3)/2；半循环指序列处理结束后残留的单调段，每个残留段计 0.5 次；循环总数按“完整循环数 + 半循环数 × 0.5”统计。理论上，对 n 个转折点的封闭序列，完整循环数与半循环数之和满足 2×完整循环 + 半循环 = n − 1（每提取一个完整循环消耗 2 个内部点，残留段数为剩余点数减 1）。

示例时程与峰谷压缩

取如下 24 点应力序列（MPa）：

s = [10, 25, 40, 30, 18, 35, 52, 44, 28, 15,
     22, 38, 50, 42, 30, 20, 33, 46, 36, 24,
     16, 29, 41, 31]

压缩规则：保留首点；对每个内部点，若 (s[i]−s[i−1]) 与 (s[i+1]−s[i]) 异号或存在零值，则该点为转折点；末点保留。压缩后得到转折点序列（索引从 0 计）：

turns = [10, 40, 18, 52, 15, 50, 20, 46, 16, 41, 31]
索引    0   2   4   6   8  12  15  17  20  22  23

共 11 个转折点。按公式，若全部处理完，完整循环数 C 与半循环数 H 应满足 2C + H = 11 − 1 = 10。

四点法实现

下面给出可直接运行的 Python 实现。算法用栈保存已确认的转折点，每次压入新点后检查栈顶四个点是否满足提取条件；满足则弹出中间两点作为完整循环，并继续回退检查。

def rainflow_four_point(turns):
    stack = []
    cycles = []          # (幅值S, 均值M, 计数)
    for x in turns:
        stack.append(x)
        while len(stack) >= 4:
            r1, r2, r3, r4 = stack[-4], stack[-3], stack[-2], stack[-1]
            a = abs(r2 - r1)
            b = abs(r3 - r2)
            c = abs(r4 - r3)
            if b <= a and b <= c:
                S = b
                M = (r2 + r3) / 2.0
                cycles.append((S, M, 1))
                del stack[-3:-1]      # 删除 R2、R3
            else:
                break
    # 残留单调段：相邻点构成半循环
    for i in range(len(stack) - 1):
        S = abs(stack[i+1] - stack[i])
        M = (stack[i] + stack[i+1]) / 2.0
        cycles.append((S, M, 0.5))
    return cycles

turns = [10, 40, 18, 52, 15, 50, 20, 46, 16, 41, 31]
res = rainflow_four_point(turns)
full = sum(1 for c in res if c[2] == 1)
half = sum(1 for c in res if c[2] == 0.5)
print("完整循环:", full, "半循环:", half, "总计数:", full + half*0.5)
for S, M, n in res:
    print(f"S={S:5.1f}  M={M:6.1f}  计数={n}")

按上述输入，程序输出的完整循环为 4 个、半循环为 2 个，总计数 4 + 2×0.5 = 5，满足 2×4 + 2 = 10 = 11 − 1，与理论关系一致。提取出的完整循环幅值依次为 22、34、32、30 MPa（对应区间 40→18、52→15、50→20、46→16），残留半循环为 10→40 与 41→31 两段。

结果结构与验证要点

输出应至少包含三列：幅值 S、均值 M、计数（1 或 0.5）。验证分三步：其一，检查 2C + H = n − 1 是否成立；其二，检查所有循环的幅值之和是否等于序列总应力范围的一半（对封闭序列成立，开放序列需按残留段修正）；其三，用同一输入分别以“原始采样序列直接压缩后计数”和“先压缩再计数”两条路径执行，结果应完全一致——这正是峰谷预处理不改变计数结果的体现。

常见实现陷阱

• 未压缩峰谷直接计数：把每个采样点都当作转折点，会得到大量虚假小循环，幅值分布严重失真。必须先做峰谷压缩。
• 比较条件写成严格小于：四点法要求 b ≤ a 且 b ≤ c，等号必须包含，否则等幅循环会被漏掉。
• 删除后未回退：提取一个循环后必须继续用栈顶四点判断，否则会漏掉嵌套循环。
• 半循环重复计数：残留段每段只计 0.5，不能把首尾两段各计 1。

完成上述实现后，把输出与 2C + H = n − 1 的关系对照，即可确认算法在计数守恒层面正确；幅值与均值的具体数值则用于后续的 S-N 曲线损伤累积。