guanlan-src/app_ontology/app_ontology_guanlan/sop/vib_impact.py

# -*- coding: utf-8 -*-
"""冲击轴 (时域六窗层) 分型 → 六枚举定级. 库层单一实现: windcms 插件与 fusion 脚本都调这里 (落库≠接上 的教训: 生产路径必须走同一函数).

口径 (2026-08-23 如东 M5 实证, 见 skill vibration-spectral §标量指标的劣化分型 / B 类):
- 基线段 = 末窗之前; 末段 = 最后一个窗 (取"末两窗"会把 8/7 阶跃稀释成间歇尖峰, 实测逮);
- 分型由 sop.discriminators.vib_scalar_typology 给: 持续恶化 / 间歇尖峰 / 稳定高位 / 正常;
- 定级: 持续恶化→候选·新发 (B1 阶跃); 稳定高位→候选 (B2); 间歇尖峰 且 剔共模后高值日≥3→候选 (B3 反复发作), 否则→参考; 正常→正常;
- 共模日 (台风/电网等外部激励): 当日超自身阈 (max(4×自身基线, 15)) 的台数 ≥ max(6, 15% 台数) ⇒ 该日不计入个体高值日;
  全部高值日都落在共模日 ⇒ 个体证据为零 → 参考 (不作设备判).
★以上阈值 (8×/2.5×/15 m/s²/共模台数) 由如东主轴承 Peak 标定, 跨测点/跨场用前须复标 [暂行 n=1 场].
"""
import numpy as np
import pandas as pd

LEVEL_BY_TYPE = {'持续恶化': '候选·新发', '稳定高位': '候选', '正常': '正常'}
VALID_SENSORS = ('Main_bearing_front', 'Main_bearing_rear')   # 阈值标定范围; 其他测点产出但标 _caliber 未标定


def impact_axis_typology(scalars, turbine, sensor, meas='Peak', event_windows=None, min_late_days=3):
    """scalars: 标量长表 (列 turbine, window, sensor_name, meas_name, trigger_time, scalar_value), 含全场 (fleet 参照与共模日都要全场).
    event_windows: 场级外部事件窗 [(start, end, label), ...] (台风等, 须有独立于 CMS 的确认: 公开气象/运行记录; 如东见 findings TYPHOON_COMMON_MODE_20260819),
        窗内日与共模探测日一样不计入个体高值日 (skill: 海上先剔天气; 共模探测器靠"≥6 台超阈"对台风边缘日 (2-4 台超阈) 不敏感, 2026-08-23 实测 25#/29# 被误升 B3).
    min_late_days: 末窗有效天数低于此 → 末段不可判 (A1 盲区嫌疑), 只按历史 B3 定级, 禁判"已归位" (2026-08-23 实测 20# 末窗 2 天被判归位→参考, 实为带病进盲区).
    返回 dict(type, level, detail, raw, late_start)."""
    from .discriminators import vib_scalar_typology
    s = scalars
    q = s[(s.sensor_name == sensor) & (s.meas_name == meas)].copy()
    if q.empty:
        return dict(type='INSUFFICIENT', level='INSUFFICIENT', detail='无标量', raw={}, late_start=None)
    q['day'] = pd.to_datetime(q.trigger_time).dt.normalize()
    mine = q[q.turbine == turbine]
    if mine.empty:
        return dict(type='INSUFFICIENT', level='INSUFFICIENT', detail=f'{turbine} 无 {sensor}/{meas} 标量', raw={}, late_start=None)
    daily = mine.groupby('day').scalar_value.agg(['max', 'median', 'count'])
    wins = sorted(q.window.unique())
    late_start = q[q.window == wins[-1]].trigger_time.min().normalize()
    fleet_med = float(q.groupby(['turbine', 'day']).scalar_value.median().groupby('turbine').median().median())
    r = vib_scalar_typology(daily, baseline_end=late_start, late_start=late_start, fleet_daily_median=fleet_med)
    typ = r.get('type', 'INSUFFICIENT')
    level = LEVEL_BY_TYPE.get(typ, 'INSUFFICIENT')   # 间歇尖峰 的级别在剔共模/事件窗之后再定 (见下), 此处先占位
    nums = {k: (round(float(v), 3) if isinstance(v, (int, float)) and np.isfinite(v) else v) for k, v in r.items() if k != '_caliber'}
    base = float(r.get('baseline_max', np.nan)) if r.get('baseline_max') is not None else np.nan
    # 共模日: 当日超自身阈值的台数 ≥ max(6, 15% 台数)
    dm = q.groupby(['turbine', 'day']).scalar_value.max().unstack('turbine')
    bases = dm[dm.index < late_start].median() if (dm.index < late_start).any() else dm.median()
    thr_t = np.maximum(4 * bases, 15.0)
    exceed = (dm >= thr_t).sum(axis=1)
    n_turb = int(dm.notna().any().sum())
    common_days = set(exceed[exceed >= max(6, int(0.15 * n_turb))].index)
    ev_days = set()
    for w in (event_windows or []):
        a, b = pd.Timestamp(w[0]).normalize(), pd.Timestamp(w[1]).normalize()
        ev_days |= set(pd.date_range(a, b))
    common_days |= ev_days
    hi_days = daily.index[daily['max'] >= max(4 * base, 15.0)] if np.isfinite(base) else []
    n_high_raw = int(len(hi_days))
    n_high_all = int(sum(1 for d_ in hi_days if d_ not in common_days))
    n_common = n_high_raw - n_high_all
    late_days = daily[daily.index >= late_start]
    n_late = int(len(late_days)); n_late_ev = int(sum(1 for d_ in late_days.index if d_ in ev_days))
    hist = ''
    # 全场大风日注释 (非闸): rms_200 全场逐日中位 ≥1.1× 全期中位 的日子; 台风日 ≥1.4×, 普通大风日 1.1–1.3× (如东 07-23/24, 08-07)
    windy_days = set()
    try:
        r2 = scalars[(scalars.sensor_name == sensor) & (scalars.meas_name == 'rms_200')].copy()
        if len(r2):
            r2['day'] = pd.to_datetime(r2.trigger_time).dt.normalize()
            fr2 = r2.groupby('day').scalar_value.median(); windy_days = set(fr2[fr2 >= 1.1 * float(fr2.median())].index)
    except Exception:
        windy_days = set()
    n_windy = int(sum(1 for d_ in hi_days if d_ not in common_days and d_ in windy_days))
    if typ == '间歇尖峰':
        level = '候选' if n_high_all >= 3 else '参考'   # B3
        if n_high_all >= 3 and n_windy == n_high_all:
            level = '参考'
            hist += f'; 剔窗后 {n_high_all} 个高值日**全部**为全场大风日 (rms_200 全场中位 ≥1.1×) 且逐日中位未抬 → 尖峰与风激励相关, 按 参考 记基线·下窗跟踪 (若非大风日再发 → 候选)'
        elif n_windy:
            hist += f'; 剔窗后高值日中 {n_windy} 天为全场大风日 (rms_200 全场中位 ≥1.1×), 尖峰与风激励相关, 候选须下窗复核'
    # B3 反复发作: ≥3 个**剔窗后**高值日 (2026-08-23 修: 原按剔窗前计数, 25#/29# 被台风边缘日误升)
    # ── B1 阶跃抬高 (skill B 类第一条): 逐日中位 分段比 ≥1.6 且 全场同期对照 ≈1.0; 分段点 = 各窗起始日 (末窗除外), 两段各 ≥7 天 ──
    step = _b1_step(q, turbine, daily, wins, ev_days=ev_days)
    if step and step.get('gap_shift'):
        g = step['gap_shift']
        if g.get('reason') == 'gap':
            hist += f"; 跨缺口抬升 {g['ratio']}× (分段 {g['split']} 前有 {g['gap_days']} 天无数据, fleet {g['fleet_ratio']}×): 时点不可定, 不作 B1 事件"
        else:
            hist += f"; 末段局部抬升 {g['ratio']}× @ {g['split']} 但剔事件窗后 {'阶跃后' if g.get('reason') == 'short_post' else '阶跃前'}仅 {g['post_days'] if g.get('reason') == 'short_post' else g['pre_days']} 天 (<{7}), 不足以判 B1 → 待下窗确认"
        nums.update(gap_shift_ratio=g['ratio'], gap_shift_split=g['split'], gap_days=g['gap_days'], gap_shift_reason=g.get('reason'))
    if step and step.get('ratio') and typ in ('正常', '间歇尖峰', '稳定高位'):
        # 阶跃发生在末窗之前时, "末窗 vs 之前"的基线已被阶跃后数据污染 (16# 7/15 阶跃实测漏判) → 以 B1 为准
        if step['still_high']:
            typ, level = '阶跃抬高', '候选·新发'
            if step.get('sparse_tail'):
                # 2026-08-24 用户裁: 回落未定 (末7日样本稀) 不挂报警级 → 记基线观察, 下窗再抬即候选·新发
                # (16# 实案: 末窗中位/基线 0.93/0.86 已在基线下, 仅因样本稀不能确认回落 — "维持候选·新发"过强)
                level = '候选·记基线'
                hist += f"; ⚠ 阶跃后末 7 个有数据日仅 {step['tail_records']} 条记录, 回落未定 → 记基线观察 (下窗再抬即候选·新发)"
        else:
            hist += f"; 曾阶跃 (分段比 {step['ratio']}× @ {step['split']}, fleet {step['fleet_ratio']}×) 但已回落 (末 7 个有数据日 {step['last7_ratio']}× < 回落阈 {1 + 0.5 * (step['ratio'] - 1):.2f}×, {step['tail_records']} 条) → 参考 (记基线, 下窗再抬即候选·新发)"
            level = max(level, '参考', key=['INSUFFICIENT', '正常', '参考', '候选·记基线', '候选', '候选·新发', '准定论·预警'].index)
        nums.update(b1_split=step['split'], b1_ratio=step['ratio'], b1_fleet_ratio=step['fleet_ratio'], b1_post_days=step['post_days'],
                    b1_last7_ratio=step.get('last7_ratio'), b1_tail_records=step.get('tail_records'))
    if n_late < min_late_days:
        # 末窗稀疏 = A1 盲区嫌疑: 末段不可判, 禁"已归位"; 级别只看历史 B3
        if n_high_all >= 3:
            typ, level = '反复发作(末段不可判)', '候选'
            hist = f'; ⚠ 末窗仅 {n_late} 天有记录 (<{min_late_days}, A1 盲区嫌疑) ⇒ 末段不可判、禁判已归位; 历史剔共模后高值日 {n_high_all} 天 ≥3 ⇒ B3 反复发作记候选 + A1 (恢复采集后复判)'
        else:
            typ, level = 'INSUFFICIENT', 'INSUFFICIENT'
            hist = f'; ⚠ 末窗仅 {n_late} 天有记录 (<{min_late_days}, A1 盲区嫌疑) ⇒ 末段不可判; 历史剔共模后高值日 {n_high_all} 天 <3 ⇒ 不可判 (非正常)'
    elif typ in ('间歇尖峰', '持续恶化') and n_high_raw and n_high_all == 0:
        hist = f'; ⚠ {n_high_raw} 个高值日全部落在全场共模日 (台风/外部激励), 个体证据为零 → 不作设备判'
        typ, level = '正常', '参考'
    elif typ == '正常' and n_high_all >= 3:
        level = '参考'
        hist = f'; 历史反复发作 {n_high_all} 天 (末窗前, 已剔共模日 {n_common} 天), 末窗已归位 → 记基线, 查 A1 盲区与外部事件窗 (B3)'
    elif n_common:
        hist = f'; 已剔全场共模/事件窗日 {n_common} 天'
    if typ == '持续恶化' and n_late_ev:
        hist += f'; ⚠ 末窗 {n_late} 天中 {n_late_ev} 天落在外部事件窗 (台风), 末段中位倍数含外部放大, 不作损伤演进速率 (起病日须在窗外才算新发)'
    cal = '' if sensor in VALID_SENSORS else '; ⚠ 阈值按主轴承 Peak 标定, 本测点未标定 (只作参考)'
    detail = (f'末段 = 末窗 {wins[-1]} 起 {late_start.date()}, 逐日 {len(daily)} 天; ' + ', '.join(f'{k}={v}' for k, v in nums.items() if k != 'type')
              + f', n_high_days_all={n_high_all}' + ('; 高值日≥3 按 B3 反复发作记候选 (需剔除台风/外部事件窗后复核)' if typ == '间歇尖峰' and n_high_all >= 3 else '') + hist + cal)
    return dict(type=typ, level=level, detail=detail, raw=dict(nums, n_high_days_all=n_high_all, n_common_days=n_common, n_late_days=n_late, n_late_days_in_event=n_late_ev, n_windy_high_days=n_windy), late_start=str(late_start.date()))


def _b1_step(q, turbine, daily, wins, ratio_min=1.6, fleet_max=1.2, min_days=7, max_gap_days=5, local_days=30, ev_days=None):
    """B1 阶跃 (skill B 类第一条): 逐日试分段点; 本台与全场的 pre/post 逐日中位都只用分段点前后各 local_days 天的**局部窗口**、并剔除外部事件窗日
    (全史 pre 会把远期低值混进基线: 如东 23# 1 月窗拉低基线造出 2.28× 假阶跃; 台风日抬 post). 成立: 比≥ratio_min ∧ 全场同分段比≤fleet_max ∧
    剔窗后两侧各≥min_days 天 ∧ 分段点前相邻缺口≤max_gap_days (否则是跨缺口 level shift, 时点不可定 → 只记 gap_shift 注释).
    still_high = 末窗 (剔事件窗日, 不足 2 天则全用) 逐日中位 / 阶跃前局部中位 ≥ ratio_min."""
    ev = set(ev_days or [])
    dd = daily[~daily.index.isin(ev)] if ev else daily
    if len(dd) < 2 * min_days:
        return None
    fm = q.groupby(['turbine', 'day']).scalar_value.median().unstack('turbine')
    fleet_daily = fm.median(axis=1)
    fleet_daily = fleet_daily[~fleet_daily.index.isin(ev)] if ev else fleet_daily
    last_start = q[q.window == wins[-1]].trigger_time.min().normalize()
    last_all = daily[daily.index >= last_start]['median']
    last = last_all[~last_all.index.isin(ev)] if ev else last_all
    if len(last) < 2:
        last = last_all
    idx = daily.index
    best, gap_shift = None, None
    for k in range(1, len(idx)):
        split = idx[k]
        lo, hi = split - pd.Timedelta(days=local_days), split + pd.Timedelta(days=local_days)
        pre = dd['median'][(dd.index >= lo) & (dd.index < split)]
        post = dd['median'][(dd.index >= split) & (dd.index < hi)]
        p0, p1 = float(np.nanmedian(pre)) if len(pre) else np.nan, float(np.nanmedian(post)) if len(post) else np.nan
        if not (np.isfinite(p0) and np.isfinite(p1)) or p0 <= 1e-9:
            continue
        r = p1 / p0
        if r < ratio_min:
            continue
        f_pre = fleet_daily[(fleet_daily.index >= lo) & (fleet_daily.index < split)]; f_post = fleet_daily[(fleet_daily.index >= split) & (fleet_daily.index < hi)]
        f0, f1 = float(np.nanmedian(f_pre)) if len(f_pre) else np.nan, float(np.nanmedian(f_post)) if len(f_post) else np.nan
        fr = (f1 / f0) if (np.isfinite(f0) and f0 > 1e-9 and np.isfinite(f1)) else np.nan
        if not (np.isfinite(fr) and fr <= fleet_max):
            continue
        gap = (split - idx[k - 1]).days
        lr = float(np.nanmedian(last)) / p0 if len(last) else np.nan
        # 回落判据: 最后 7 个有数据日 (含事件窗日) 的中位 / 阶跃前局部基线; 回落 = 让回超过阶跃幅度一半 (阈 1+0.5(r-1));
        # 这 7 天记录 <5 条 → 样本稀, 回落未定 (16# 8/9-8/11 每日 1-2 条实测), 维持阶跃判; 36# 8/1-8/7 隆起后 8/8-8/11 回到基线 (每日多条) = 真回落
        tail = daily.iloc[-7:]
        l7 = float(np.nanmedian(tail['median'])) / p0
        n_tail = int(tail['count'].sum()) if 'count' in tail else 99
        thr_rec = 1.0 + 0.5 * (r - 1.0)
        sparse_tail = n_tail < 5
        rec = dict(split=str(split.date()), ratio=round(r, 2), fleet_ratio=round(fr, 2), pre_days=int(len(pre)), post_days=int(len(post)), gap_days=int(gap),
                   last_ratio=(round(lr, 2) if np.isfinite(lr) else None), last7_ratio=round(l7, 2), tail_records=n_tail, sparse_tail=sparse_tail,
                   still_high=bool(sparse_tail or l7 >= thr_rec))
        if gap > max_gap_days or len(pre) < min_days or len(post) < min_days:
            rec['reason'] = 'gap' if gap > max_gap_days else ('short_post' if len(post) < min_days else 'short_pre')
            if gap_shift is None or r > gap_shift['ratio']:
                gap_shift = rec
        elif best is None or r > best['ratio']:
            best = rec
    if best:
        best['gap_shift'] = gap_shift
    return best or (dict(gap_shift=gap_shift) if gap_shift else None)