4. 雨流计数循环的滞回环闭合判定与可视化

滞回环闭合的物理含义与判定条件

雨流计数把应力时程中的一次“加载—卸载”往返识别为一个完整循环。在应力-应变平面上，这个往返对应一条闭合的滞回环：从某个起点出发，沿加载路径到达峰值，再沿卸载路径回到与起点同一应力水平的点，形成闭合。闭合判定的核心是应力水平相等，而不是时间点相同——雨流计数中的“闭合点”通常出现在起点之后的某个峰谷序列位置，其应力值与起点相等（在离散采样下允许容差）。

设峰谷序列为 \(P_0, P_1, \dots, P_n\)（应力值），雨流计数中一个完整循环由四个特征点构成：起点 \(P_i\)、峰值 \(P_{i+1}\)、谷值 \(P_{i+2}\)、闭合点 \(P_{i+3}\)，且满足闭合条件：

循环幅值 \(S_a = (P_{i+1} − P_{i+2})/2\)，均值 \(S_m = (P_{i+1} + P_{i+2})/2\)。注意：幅值和均值由峰谷两点决定，而闭合点只用于确认循环边界，不参与幅值计算。这是本案例验证的关键——闭合判定正确不等于幅值均值提取正确，两者需分别校验。

实现：记录起点、闭合点与滞回环路径

下面用四点法实现雨流计数，并在识别出完整循环时记录起点索引、闭合点索引、峰值、谷值以及滞回环的路径点序列。代码可直接运行，输入为一段教学示例峰谷序列（模拟海洋结构浪致应力波动）。

import numpy as np

def rainflow_with_hysteresis(peaks, eps=0.5):
    """
    peaks: 峰谷序列（应力值，MPa），已去除相邻等值点
    eps:   闭合容差（MPa）
    返回: 循环列表，每个元素含 start_idx, close_idx, peak, valley, Sa, Sm, path
    """
    stack = []          # 栈中存 (索引, 应力值)
    cycles = []
    for idx, val in enumerate(peaks):
        stack.append((idx, val))
        # 四点法：检查栈顶四个点是否构成闭合循环
        while len(stack) >= 4:
            i0, v0 = stack[-4]
            i1, v1 = stack[-3]
            i2, v2 = stack[-2]
            i3, v3 = stack[-1]
            # 闭合判定：起点与闭合点应力相等（容差内）
            if abs(v3 - v0) <= eps:
                # 峰谷交替：v1 为峰、v2 为谷（或反之）
                peak, valley = (v1, v2) if v1 > v2 else (v2, v1)
                Sa = (peak - valley) / 2.0
                Sm = (peak + valley) / 2.0
                # 滞回环路径：起点→峰→谷→闭合点
                path = [(i0, v0), (i1, v1), (i2, v2), (i3, v3)]
                cycles.append({
                    'start_idx': i0, 'close_idx': i3,
                    'peak': peak, 'valley': valley,
                    'Sa': Sa, 'Sm': Sm, 'path': path
                })
                # 弹出中间两点，保留起点与闭合点（闭合点成为新起点）
                stack.pop()   # 弹出 i3
                stack.pop()   # 弹出 i2
                stack.pop()   # 弹出 i1
                # 此时栈顶为 i0，i3 已作为新起点压回
                stack.append((i3, v3))
            else:
                break
    return cycles

# 教学示例峰谷序列（MPa），模拟浪致应力波动
peaks = np.array([0.0, 80.0, -40.0, 60.0, -20.0, 100.0, -60.0, 40.0, -80.0, 20.0, 0.0])
cycles = rainflow_with_hysteresis(peaks, eps=0.5)

for c in cycles:
    print(f"起点#{c['start_idx']}→闭合#{c['close_idx']}  "
          f"峰={c['peak']:.1f} 谷={c['valley']:.1f}  "
          f"Sa={c['Sa']:.1f} Sm={c['Sm']:.1f}  "
          f"路径={c['path']}")

运行后应得到若干完整循环。以示例序列为例，可观察到：起点 0.0 与闭合点 0.0 构成一个循环，峰 80、谷 -40，Sa=60、Sm=20；随后闭合点 0.0 作为新起点继续参与后续计数。若把 eps 改为 0.1，某些近似闭合的循环将不再被识别，说明容差直接影响循环数量——这正是海洋结构实测数据中需要谨慎设定容差的原因。

滞回环可视化与闭合判定校验

下图把识别出的循环在应力-应变平面上绘制为滞回环，并标注起点、峰值、谷值、闭合点。应力-应变关系采用简化的线性硬化模型（弹性模量 E、塑性斜率 H 为教学假设值），仅用于示意滞回环形状，不表示真实材料本构。

闭合判定与幅值均值提取的分别校验

可视化之后需要做两项独立校验，避免“图看起来闭合了”就认为结果正确：

• 闭合判定校验：对每个识别出的循环，检查 |P_close − P_start| ≤ ε 是否成立，并统计被容差拒绝的候选循环数量。若拒绝数量异常多，说明 ε 过小或峰谷序列存在漂移。
• 幅值均值校验：用峰谷两点独立重算 Sa、Sm，与记录值比对；同时检查循环总数与峰谷序列长度的关系（完整循环数 ≤ (n−1)/2），防止重复计数或漏计。

# 校验：闭合条件与幅值均值一致性
ok_close, ok_amp = 0, 0
for c in cycles:
    p_start = peaks[c['start_idx']]
    p_close = peaks[c['close_idx']]
    if abs(p_close - p_start) <= 0.5:
        ok_close += 1
    Sa_check = (c['peak'] - c['valley']) / 2.0
    Sm_check = (c['peak'] + c['valley']) / 2.0
    if abs(Sa_check - c['Sa']) < 1e-9 and abs(Sm_check - c['Sm']) < 1e-9:
        ok_amp += 1

print(f"闭合判定通过: {ok_close}/{len(cycles)}")
print(f"幅值均值一致: {ok_amp}/{len(cycles)}")
print(f"循环总数: {len(cycles)}，峰谷点数: {len(peaks)}")