#!/usr/bin/env python3 """Long-term indicator-feature deviation vs the scanner's anomaly segments. For each cylinder/side it reduces the feature CSV to daily medians, scores them as robust z against the point's own baseline (median/MAD over its valid days), and marks days that fall inside a rule anomaly segment from cache/pks_scanner/ evolution (pressure 阶段段 / 振动异常段). Read-only. """ from __future__ import annotations import argparse import csv from collections import defaultdict from datetime import datetime from pathlib import Path import numpy as np UNIT_BATCH = {"7": 30, "8": 31, "9": 32} TREND_FEATURES = ("p_discharge", "pressure_ratio", "n_comp", "efficiency") def parse_dt(value): return datetime.strptime(value[:19], "%Y-%m-%d %H:%M:%S") def load_features(path): rows = list(csv.DictReader(open(path, encoding="utf-8-sig"))) for row in rows: row["day"] = row["sample_time"][:10] for name in TREND_FEATURES + ("indicated_work",): try: row[name] = float(row[name]) except (TypeError, ValueError): row[name] = None return rows def load_segments(evolution_dir, unit): """Rule anomaly segments for a unit from the two stage CSVs.""" segments = [] for name, kind in (("润滑油压力小时阶段段.csv", "压力"), ("振动异常段.csv", "振动")): path = Path(evolution_dir) / name if not path.is_file(): continue for row in csv.DictReader(open(path, encoding="utf-8-sig")): if row["机组"] != f"{unit}号机": continue segments.append({ "kind": kind, "level": row.get("压力阶段") or row.get("振动阶段"), "start": parse_dt(row["开始小时"]), "end": parse_dt(row["结束小时"]), }) return segments def mark_segments(day: str, segments): moment = datetime.strptime(day, "%Y-%m-%d") return [s["kind"] + "/" + str(s["level"]) for s in segments if s["start"] <= moment <= s["end"]] def peer_sudden_changes(rows, feature, segments): """Day-to-day jumps of the peer-relative deviation (cancels operating level).""" by_key = defaultdict(dict) for row in rows: if row[feature] is not None: by_key[(row["stage"], row["side"], row["day"])][row["cylinder"]] = row[feature] series = defaultdict(dict) for (_stage, side, day), members in by_key.items(): values = [v for v in members.values() if v is not None] if len(values) < 2: continue med = float(np.median(values)) if abs(med) < 1e-9: continue for cylinder, value in members.items(): if value is not None: series[(cylinder, side)][day] = (value - med) / med print(f"\n== 同级缸间偏差 突发变化（{feature}）==") for (cylinder, side), daily in sorted(series.items()): days = sorted(daily) if len(days) < 3: continue values = np.array([daily[d] for d in days], dtype=float) scale = float(np.median(np.abs(values - np.median(values)))) * 1.4826 or float(np.std(values)) if scale <= 1e-9: continue for i in range(1, len(days)): gap = (datetime.strptime(days[i], "%Y-%m-%d") - datetime.strptime(days[i - 1], "%Y-%m-%d")).days if gap > 3: continue delta = values[i] - values[i - 1] if abs(delta) / scale >= 3.0: flags = mark_segments(days[i], segments) print(f" {cylinder}缸{side}: 偏差 {values[i-1]*100:+.1f}% → {values[i]*100:+.1f}% 于 {days[i]} ({delta/scale:+.1f}×) {' '.join(flags)}") def _fit_drift(days, values): x = np.array([(datetime.strptime(d, "%Y-%m-%d") - datetime.strptime(days[0], "%Y-%m-%d")).days for d in days], dtype=float) y = np.asarray(values, dtype=float) if x[-1] <= 0 or len(y) < 8: return None slope, intercept = np.polyfit(x, y, 1) pred = slope * x + intercept ss_res = float(np.sum((y - pred) ** 2)) ss_tot = float(np.sum((y - y.mean()) ** 2)) r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 0.0 center = float(np.median(y)) rel30 = (slope * 30.0 / center) if abs(center) > 1e-9 else 0.0 return rel30, r2, float(y[0]), float(y[-1]) def slow_trends(rows, features=("p_discharge", "pressure_ratio", "efficiency", "n_comp", "indicated_work")): """Slow, weak drifts: linear fit per (cylinder, side) and per peer-relative series.""" print("\n== 慢趋势：自基线（斜率 %/30天，要求 |%|>=3 且 R2>=0.3） ==") by_series = defaultdict(list) for row in rows: by_series[(row["cylinder"], row["side"])].append(row) for (cylinder, side), series in sorted(by_series.items()): days = sorted({r["day"] for r in series}) for feature in features: daily = {d: np.median([r[feature] for r in series if r["day"] == d and r[feature] is not None]) for d in days} present = [d for d in days if daily.get(d) is not None] fit = _fit_drift(present, [daily[d] for d in present]) if fit and abs(fit[0]) >= 0.03 and fit[1] >= 0.3: print(f" {cylinder}缸{side} {feature:<16} {fit[0]*100:+.1f}%/30天 R2={fit[1]:.2f} ({fit[2]:.3f}->{fit[3]:.3f})") print("\n== 慢趋势：同级 peer 相对（%/30天） ==") by_ps = defaultdict(dict) for row in rows: if row["pressure_ratio"] is not None: by_ps[(row["stage"], row["side"], row["day"])][row["cylinder"]] = row["pressure_ratio"] rel = defaultdict(dict) for (stage, side, day), members in by_ps.items(): clean = {c: v for c, v in members.items() if v is not None} if len(clean) < 2: continue med = float(np.median(list(clean.values()))) if abs(med) < 1e-9: continue for cylinder, value in clean.items(): rel[(cylinder, side)][day] = (value - med) / med for (cylinder, side), daily in sorted(rel.items()): days = sorted(daily) fit = _fit_drift(days, [daily[d] for d in days]) if fit and abs(fit[0]) >= 0.03 and fit[1] >= 0.3: print(f" {cylinder}缸{side} pressure_ratio(peer) {fit[0]*100:+.1f}%/30天 R2={fit[1]:.2f} ({fit[2]*100:+.1f}%->{fit[3]*100:+.1f}%)") def main(): parser = argparse.ArgumentParser(description="功图特征长期偏离检测（对齐 evolution 异常段）") parser.add_argument("csv_path") parser.add_argument("--evolution-dir", default=str(Path(__file__).resolve().parents[2] / "cache" / "pks_scanner" / "evolution")) parser.add_argument("--unit", required=True, choices=("7", "8", "9")) args = parser.parse_args() rows = load_features(args.csv_path) segments = load_segments(args.evolution_dir, args.unit) by_series = defaultdict(list) for row in rows: by_series[(row["cylinder"], row["side"])].append(row) print(f"文件 {args.csv_path} 行 {len(rows)} evolution段 {len(segments)}") for (cylinder, side), series in sorted(by_series.items()): days = sorted({r["day"] for r in series}) print(f"\n== {cylinder}缸{side} ({len(days)} 天)") for feature in TREND_FEATURES: daily = {} for day in days: values = [r[feature] for r in series if r["day"] == day and r[feature] is not None] if values: daily[day] = float(np.median(values)) if len(daily) < 5: continue values = np.array([daily[d] for d in days if d in daily], dtype=float) center = float(np.median(values)) mad = float(np.median(np.abs(values - center))) scale = mad * 1.4826 or float(np.std(values)) if scale <= 1e-9: continue print(f" {feature:<16} 基线 {center:.3f}±{scale:.3f}") for day in days: if day not in daily: continue z = (daily[day] - center) / scale flags = mark_segments(day, segments) if abs(z) >= 2.0 or flags: tag = " ".join(flags) print(f" {day} {daily[day]:>8.3f} z={z:+.2f} {tag}") # Abrupt day-to-day change (sudden step), independent of any trend. present = [d for d in days if d in daily] jumps = [] for i in range(1, len(present)): gap = (datetime.strptime(present[i], "%Y-%m-%d") - datetime.strptime(present[i - 1], "%Y-%m-%d")).days if gap <= 3: jumps.append((present[i], daily[present[i]] - daily[present[i - 1]])) if jumps: day_jump, delta = max(jumps, key=lambda item: abs(item[1])) if abs(delta) / scale >= 3.0: flags = mark_segments(day_jump, segments) tag = " ".join(flags) print(f" [突发] {feature}: Δ={delta:+.3f} ({delta / scale:+.1f}×基线尺度) 于 {day_jump} {tag}") peer_sudden_changes(rows, "p_discharge", segments) slow_trends(rows) if __name__ == "__main__": main()