|
|
@@ -0,0 +1,579 @@
|
|
|
+#!/usr/bin/env python3
|
|
|
+"""Pseudo-online local-feature detection with past-only baselines."""
|
|
|
+from __future__ import annotations
|
|
|
+
|
|
|
+import csv
|
|
|
+import argparse
|
|
|
+import math
|
|
|
+import sys
|
|
|
+from collections import defaultdict, deque
|
|
|
+from datetime import date, datetime, timedelta
|
|
|
+from pathlib import Path
|
|
|
+
|
|
|
+import numpy as np
|
|
|
+
|
|
|
+BACKEND = Path(__file__).resolve().parents[1]
|
|
|
+sys.path.insert(0, str(BACKEND))
|
|
|
+
|
|
|
+from indicator_features.extract_from_db import ( # noqa: E402
|
|
|
+ CYLINDERS, SIDES, connect, fetch_cycles, fetch_file_rows, point_name, sample_per_day,
|
|
|
+)
|
|
|
+from indicator_features.local_evolution_analysis import ( # noqa: E402
|
|
|
+ FAMILY, METRICS, METRIC_CN, normalized_bins, robust_scale,
|
|
|
+)
|
|
|
+
|
|
|
+BAD = {
|
|
|
+ "4": {("二", "盖侧"), ("四", "盖侧"), ("六", "盖侧")},
|
|
|
+ "6": {("二", "盖侧"), ("四", "盖侧"), ("六", "盖侧")},
|
|
|
+ "7": {("四", "盖侧"), ("六", "盖侧")},
|
|
|
+}
|
|
|
+FAILURE_TEMPLATES = {
|
|
|
+ "7": {
|
|
|
+ ("一缸盖侧",240,1),("一缸盖侧",250,1),("一缸盖侧",260,1),
|
|
|
+ ("二缸盖侧",140,-1),("三缸轴侧",160,1),("三缸轴侧",170,1),
|
|
|
+ ("三缸盖侧",20,1),("三缸盖侧",30,1),("三缸盖侧",310,-1),("三缸盖侧",320,-1),
|
|
|
+ },
|
|
|
+ "4": {("四缸轴侧",60,-1)},
|
|
|
+ "6": {
|
|
|
+ ("三缸轴侧",160,1),("三缸轴侧",170,1),("三缸轴侧",120,1),("三缸轴侧",130,1),
|
|
|
+ ("一缸轴侧",290,1),("一缸盖侧",220,-1),("三缸轴侧",50,1),
|
|
|
+ },
|
|
|
+}
|
|
|
+EVENTS = [
|
|
|
+ {"id": "7号_2025维修过程", "unit": "7", "event": "2025-08-08",
|
|
|
+ "baseline": ("2025-04-29", "2025-06-14"), "detect": ("2025-07-25", "2025-07-28"),
|
|
|
+ "note": "事件附近缺测,仅验证7月底运行段状态"},
|
|
|
+ {"id": "7号_2026-05-21", "unit": "7", "event": "2026-05-21",
|
|
|
+ "baseline": ("2026-02-10", "2026-04-05"), "detect": ("2026-05-06", "2026-05-22"), "note": ""},
|
|
|
+ {"id": "4号_2026-06-10", "unit": "4", "event": "2026-06-10",
|
|
|
+ "baseline": ("2026-04-21", "2026-05-31"), "detect": ("2026-06-06", "2026-06-11"), "note": ""},
|
|
|
+ {"id": "6号_2026-04-22_26", "unit": "6", "event": "2026-04-22",
|
|
|
+ "baseline": ("2026-04-12", "2026-04-15"), "detect": ("2026-04-15", "2026-04-27"),
|
|
|
+ "note": "历史不足,仅3个基线日期"},
|
|
|
+]
|
|
|
+UNITS = ("4", "6", "7", "8", "9")
|
|
|
+MERGE_GAP = 2
|
|
|
+HEALTH_SEGMENTS = {
|
|
|
+ "7": [("2026-03-12", "2026-04-05")],
|
|
|
+ "8": [("2026-04-18", "2026-05-07")],
|
|
|
+ "9": [("2026-04-01", "2026-04-19")],
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+def calendar_segments(days):
|
|
|
+ segments=[]
|
|
|
+ for day in sorted(days):
|
|
|
+ if segments and (day-segments[-1][1]).days <= MERGE_GAP+1:
|
|
|
+ segments[-1][1]=day
|
|
|
+ else:
|
|
|
+ segments.append([day,day])
|
|
|
+ return [(start,end) for start,end in segments]
|
|
|
+
|
|
|
+
|
|
|
+def date_ranges(event):
|
|
|
+ start,end=event["detect"]
|
|
|
+ cursor=date.fromisoformat(start);stop=date.fromisoformat(end)
|
|
|
+ while cursor<stop:
|
|
|
+ yield cursor.isoformat()
|
|
|
+ cursor+=timedelta(days=1)
|
|
|
+
|
|
|
+
|
|
|
+def process_unit(connection, unit, start, end, per_day_baseline=3):
|
|
|
+ points={}
|
|
|
+ running_days=set()
|
|
|
+ for cylinder in CYLINDERS:
|
|
|
+ for side in SIDES:
|
|
|
+ if (cylinder,side) in BAD.get(unit,set()):continue
|
|
|
+ files=fetch_file_rows(connection,point_name(unit,cylinder,side),start,end)
|
|
|
+ running_days.update(row["sample_time"].date() for row in files if row.get("rpm") and row["rpm"]>0)
|
|
|
+ baseline_files=sample_per_day(files,per_day_baseline)
|
|
|
+ cycles=fetch_cycles(connection,[int(row["id"]) for row in baseline_files],100) if baseline_files else {}
|
|
|
+ baseline_records=[]
|
|
|
+ for row in baseline_files:
|
|
|
+ values=stats_for_cycle(cycles.get(int(row["id"]),np.array([])))
|
|
|
+ if values is None:continue
|
|
|
+ time=row["sample_time"] if isinstance(row["sample_time"],datetime) else datetime.fromisoformat(str(row["sample_time"]))
|
|
|
+ baseline_records.append({"time":time,"day":time.date().isoformat(),"values":values})
|
|
|
+ if baseline_records:
|
|
|
+ points[f"{cylinder}缸{side}"]={"baseline_records":baseline_records,"detection_records":[]}
|
|
|
+ return points,calendar_segments(running_days)
|
|
|
+
|
|
|
+
|
|
|
+def stats_for_cycle(signal):
|
|
|
+ bins = normalized_bins(signal)
|
|
|
+ if bins is None:
|
|
|
+ return None
|
|
|
+ return {(item["bin_start"], metric): item[metric] for item in bins for metric in METRICS}
|
|
|
+
|
|
|
+
|
|
|
+def point_rows(connection, event, cylinder, side, baseline_per_day=3):
|
|
|
+ point = point_name(event["unit"], cylinder, side)
|
|
|
+ baseline_files = sample_per_day(
|
|
|
+ fetch_file_rows(connection, point, event["baseline"][0], event["baseline"][1]), baseline_per_day
|
|
|
+ )
|
|
|
+ detection_files = fetch_file_rows(connection, point, event["detect"][0], event["detect"][1])
|
|
|
+ cycles = fetch_cycles(connection, [int(row["id"]) for row in baseline_files], 100) if baseline_files else {}
|
|
|
+ baseline, detection_meta = [], []
|
|
|
+ baseline_ids = {int(row["id"]) for row in baseline_files}
|
|
|
+ for row in baseline_files + detection_files:
|
|
|
+ file_id = int(row["id"])
|
|
|
+ time = row["sample_time"] if isinstance(row["sample_time"], datetime) else datetime.fromisoformat(str(row["sample_time"]))
|
|
|
+ if file_id in baseline_ids:
|
|
|
+ values = stats_for_cycle(cycles.get(file_id, np.array([])))
|
|
|
+ if values is None: continue
|
|
|
+ baseline.append({"time":time,"day":time.date().isoformat(),"slot":time.hour//6,"values":values})
|
|
|
+ else:
|
|
|
+ detection_meta.append({"id":file_id,"time":time,"day":time.date().isoformat(),"slot":time.hour//6})
|
|
|
+ return baseline, detection_meta
|
|
|
+
|
|
|
+
|
|
|
+def build_baseline(records):
|
|
|
+ per_key_day = defaultdict(list)
|
|
|
+ per_key_cycles = defaultdict(list)
|
|
|
+ for record in records:
|
|
|
+ for key, value in record["values"].items():
|
|
|
+ per_key_day[(key, record["day"])].append(value)
|
|
|
+ per_key_cycles[key].append(value)
|
|
|
+ daily = defaultdict(list)
|
|
|
+ for (key, _), values in per_key_day.items():
|
|
|
+ daily[key].append(float(np.median(values)))
|
|
|
+ provisional = {}
|
|
|
+ for key in [(angle, metric) for angle in range(0, 360, 10) for metric in METRICS]:
|
|
|
+ days = daily.get(key, [])
|
|
|
+ cycles = per_key_cycles.get(key, [])
|
|
|
+ provisional[key] = {
|
|
|
+ "center": float(np.median(days)) if days else float("nan"),
|
|
|
+ "scale": max(robust_scale(days), robust_scale(cycles)),
|
|
|
+ "days": len(days), "cycles": len(cycles),
|
|
|
+ }
|
|
|
+ for metric in METRICS:
|
|
|
+ scales = [provisional[(angle, metric)]["scale"] for angle in range(0, 360, 10)
|
|
|
+ if provisional[(angle, metric)]["scale"] > 1e-12]
|
|
|
+ floor = .2 * float(np.median(scales)) if scales else 1e-9
|
|
|
+ for angle in range(0, 360, 10):
|
|
|
+ provisional[(angle, metric)]["scale"] = max(provisional[(angle, metric)]["scale"], floor)
|
|
|
+ return provisional
|
|
|
+
|
|
|
+
|
|
|
+def aggregate(records, group_field):
|
|
|
+ grouped = defaultdict(list)
|
|
|
+ for record in records:
|
|
|
+ group = record[group_field] if group_field != "day_slot" else f"{record['day']}#{record['slot']}"
|
|
|
+ for key, value in record["values"].items():
|
|
|
+ grouped[(group, key)].append(value)
|
|
|
+ return {(group, key): float(np.median(values)) for (group, key), values in grouped.items()}
|
|
|
+
|
|
|
+
|
|
|
+def merge_bands(alerts):
|
|
|
+ groups = defaultdict(list)
|
|
|
+ for row in alerts:
|
|
|
+ groups[(row["事件"], row["日期"], row["测点"], row["等级"])].append(row)
|
|
|
+ output = []
|
|
|
+ for key, rows in groups.items():
|
|
|
+ rows.sort(key=lambda row: int(row["角度箱"]))
|
|
|
+ runs = []
|
|
|
+ for row in rows:
|
|
|
+ angle = int(row["角度箱"])
|
|
|
+ if runs and angle == int(runs[-1][-1]["角度箱"]) + 10:
|
|
|
+ runs[-1].append(row)
|
|
|
+ else:
|
|
|
+ runs.append([row])
|
|
|
+ if len(runs) > 1 and int(runs[0][0]["角度箱"]) == 0 and int(runs[-1][-1]["角度箱"]) == 350:
|
|
|
+ runs[0] = runs[-1] + runs[0]; runs.pop()
|
|
|
+ for run in runs:
|
|
|
+ starts = [int(row["角度箱"]) for row in run]
|
|
|
+ band = f"{starts[0]}-360°及0-{starts[-1]+10}°" if starts[0] > starts[-1] else f"{starts[0]}-{starts[-1]+10}°"
|
|
|
+ best = max(run, key=lambda row: float(row["综合RMS"]))
|
|
|
+ output.append({
|
|
|
+ "事件": key[0], "日期": key[1], "测点": key[2], "等级": key[3], "角度带": band,
|
|
|
+ "角度箱数": len(run), "最大支持特征数": max(int(row["支持特征数"]) for row in run),
|
|
|
+ "最大特征族数": max(int(row["支持特征族数"]) for row in run),
|
|
|
+ "最大幅度平均Z": max(float(row["幅度平均Z"]) for row in run),
|
|
|
+ "最大综合RMS": max(float(row["综合RMS"]) for row in run),
|
|
|
+ "最大异常子窗数": max(int(row["异常6小时子窗数"]) for row in run),
|
|
|
+ "代表特征": best["特征摘要"],
|
|
|
+ })
|
|
|
+ return output
|
|
|
+
|
|
|
+
|
|
|
+def run(out_dir: Path, selected_ids=None, slow_channel=False):
|
|
|
+ all_detail, all_bands, validations = [], [], []
|
|
|
+ reuse_dir=Path(str(out_dir)+"_reuse")
|
|
|
+ reuse_dir.mkdir(parents=True,exist_ok=True)
|
|
|
+ connection = connect()
|
|
|
+ try:
|
|
|
+ for event in EVENTS:
|
|
|
+ if selected_ids and event["id"] not in selected_ids:
|
|
|
+ continue
|
|
|
+ point_data = {}
|
|
|
+ for cylinder in CYLINDERS:
|
|
|
+ for side in SIDES:
|
|
|
+ if (cylinder, side) in BAD[event["unit"]]:
|
|
|
+ continue
|
|
|
+ baseline, detection_meta = point_rows(connection, event, cylinder, side)
|
|
|
+ if baseline and detection_meta:
|
|
|
+ point_data[f"{cylinder}缸{side}"]={"baseline":build_baseline(baseline),"baseline_records":baseline,
|
|
|
+ "detection_meta":detection_meta}
|
|
|
+ early_baselines={}
|
|
|
+ for point,data in point_data.items():
|
|
|
+ early_days=sorted({record["day"] for record in data["baseline_records"]})[:3]
|
|
|
+ early_records=[record for record in data["baseline_records"] if record["day"] in early_days]
|
|
|
+ early_baselines[point]=build_baseline(early_records) if early_records else {}
|
|
|
+ detection_days=sorted({record["day"] for data in point_data.values() for record in data["detection_meta"]})
|
|
|
+ history=defaultdict(deque)
|
|
|
+ slow_history=defaultdict(deque)
|
|
|
+ event_alerts=[]
|
|
|
+ for day in detection_days:
|
|
|
+ raw_by_point={}
|
|
|
+ slot_z=defaultdict(dict)
|
|
|
+ for point,data in point_data.items():
|
|
|
+ meta=[record for record in data["detection_meta"] if record["day"]==day]
|
|
|
+ ids=[record["id"] for record in meta]
|
|
|
+ cycles=fetch_cycles(connection,ids,100) if ids else {}
|
|
|
+ records=[]
|
|
|
+ for record in meta:
|
|
|
+ values=stats_for_cycle(cycles.get(record["id"],np.array([])))
|
|
|
+ if values is not None:
|
|
|
+ records.append({**record,"values":values})
|
|
|
+ if not records:continue
|
|
|
+ window_z={}
|
|
|
+ for angle in range(0,360,10):
|
|
|
+ for metric in METRICS:
|
|
|
+ key=(angle,metric);base=data["baseline"][key]
|
|
|
+ if not np.isfinite(base["center"]):continue
|
|
|
+ vals=[rec["values"][key] for rec in records if key in rec["values"]]
|
|
|
+ if not vals:continue
|
|
|
+ window_z[key]=(float(np.median(vals))-base["center"])/base["scale"]
|
|
|
+ for slot in range(4):
|
|
|
+ slot_vals=[rec["values"][key] for rec in records if rec["slot"]==slot and key in rec["values"]]
|
|
|
+ if slot_vals:slot_z[(slot,key)][point]= (float(np.median(slot_vals))-base["center"])/base["scale"]
|
|
|
+ raw_by_point[point]=window_z
|
|
|
+ common={}
|
|
|
+ for (slot,key),values in slot_z.items():common[(slot,key)]=float(np.median(list(values.values())))
|
|
|
+ common24={}
|
|
|
+ for key in {(a,m) for point in raw_by_point for a,m in raw_by_point[point]}:
|
|
|
+ vals=[raw_by_point[point][key] for point in raw_by_point if key in raw_by_point[point]]
|
|
|
+ common24[key]=float(np.median(vals)) if vals else 0.0
|
|
|
+ by_candidate=defaultdict(list)
|
|
|
+ for point,window_z in raw_by_point.items():
|
|
|
+ for (angle,metric),z in window_z.items():
|
|
|
+ row={"point":point,"angle":angle,"metric":metric,"z":z,"local":z-common24[(angle,metric)],"base":point_data[point]["baseline"][(angle,metric)]}
|
|
|
+ by_candidate[(point,angle)].append(row)
|
|
|
+ abnormal_bins_by_point=defaultdict(set)
|
|
|
+ for (point,angle),rows in by_candidate.items():
|
|
|
+ if any(abs(row["local"])>=6 for row in rows):
|
|
|
+ abnormal_bins_by_point[point].add(angle)
|
|
|
+ channel_fault={point for point,bins in abnormal_bins_by_point.items() if len(bins)>=18}
|
|
|
+ for (point,angle),rows in sorted(by_candidate.items(),key=lambda item:(item[0][0],item[0][1])):
|
|
|
+ supported = [row for row in rows if abs(row["local"]) >= 2]
|
|
|
+ families = {FAMILY[row["metric"]] for row in supported}
|
|
|
+ amplitudes = [abs(row["local"]) for row in supported]
|
|
|
+ mean_amp = float(np.mean(amplitudes)) if amplitudes else 0.0
|
|
|
+ rms = float(np.sqrt(np.mean(np.square(amplitudes)))) if amplitudes else 0.0
|
|
|
+ abnormal_slots=0
|
|
|
+ for slot in range(4):
|
|
|
+ slot_values=[]
|
|
|
+ for row in rows:
|
|
|
+ slot_key=(slot,(angle,row["metric"]))
|
|
|
+ if point in slot_z.get(slot_key,{}):
|
|
|
+ slot_values.append(abs(slot_z[slot_key][point]-common.get(slot_key,0.0)))
|
|
|
+ if sum(value>=2 for value in slot_values)>=2:abnormal_slots+=1
|
|
|
+ candidate = len(families) >= 2 and mean_amp >= 2 and rms >= 3 and abnormal_slots >= 2
|
|
|
+ hist = history[(point, angle)]; hist.append(candidate)
|
|
|
+ if len(hist) > 5: hist.popleft()
|
|
|
+ persistent = len(hist) >= 3 and sum(hist) >= 3
|
|
|
+
|
|
|
+ slow_rows=[]
|
|
|
+ for row in rows:
|
|
|
+ metric=row["metric"];key=(angle,metric)
|
|
|
+ early=early_baselines[point].get(key)
|
|
|
+ if not early or not np.isfinite(early["center"]):continue
|
|
|
+ current_value=row["z"]*row["base"]["scale"]+row["base"]["center"]
|
|
|
+ scale=max(early["scale"],row["base"]["scale"]*.2)
|
|
|
+ slow_z=(current_value-early["center"])/scale
|
|
|
+ slow_rows.append({**row,"slow_z":slow_z})
|
|
|
+ slow_public_values=defaultdict(list)
|
|
|
+ for other_point,other_rows in by_candidate.items():
|
|
|
+ other_angle=other_point[1]
|
|
|
+ if other_angle!=angle:continue
|
|
|
+ for other_row in other_rows:
|
|
|
+ other_early=early_baselines[other_point[0]].get((angle,other_row["metric"]))
|
|
|
+ if not other_early or not np.isfinite(other_early["center"]):continue
|
|
|
+ value=other_row["z"]*other_row["base"]["scale"]+other_row["base"]["center"]
|
|
|
+ scale=max(other_early["scale"],other_row["base"]["scale"]*.2)
|
|
|
+ slow_public_values[other_row["metric"]].append((value-other_early["center"])/scale)
|
|
|
+ slow_supported_rows=[]
|
|
|
+ for row in slow_rows:
|
|
|
+ public=float(np.median(slow_public_values[row["metric"]])) if slow_public_values[row["metric"]] else 0.0
|
|
|
+ row["slow_local"]=row["slow_z"]-public
|
|
|
+ if abs(row["slow_local"])>=2:slow_supported_rows.append(row)
|
|
|
+ slow_families={FAMILY[row["metric"]] for row in slow_supported_rows}
|
|
|
+ slow_amplitudes=[abs(row["slow_local"]) for row in slow_supported_rows]
|
|
|
+ slow_mean=float(np.mean(slow_amplitudes)) if slow_amplitudes else 0.0
|
|
|
+ slow_rms=float(np.sqrt(np.mean(np.square(slow_amplitudes)))) if slow_amplitudes else 0.0
|
|
|
+ slow_candidate=len(slow_families)>=2 and slow_mean>=2 and slow_rms>=3
|
|
|
+ slow_hist=slow_history[(point,angle)];slow_hist.append(slow_candidate)
|
|
|
+ if len(slow_hist)>5:slow_hist.popleft()
|
|
|
+ slow_persistent=len(slow_hist)>=3 and sum(slow_hist)>=3
|
|
|
+ raw_supported=supported or slow_supported_rows
|
|
|
+ direction_values=[row["local"] for row in supported] or [row["slow_local"] for row in slow_supported_rows]
|
|
|
+ direction=1 if direction_values and np.median(direction_values)>=0 else -1
|
|
|
+ signature_match=any(tp==point and ta==angle and td==direction
|
|
|
+ for tp,ta,td in FAILURE_TEMPLATES.get(event["unit"],set()))
|
|
|
+ candidate_level=(candidate or slow_candidate)
|
|
|
+ level="通道质控" if point in channel_fault else "持续预警" if ((persistent and candidate) or (slow_persistent and slow_candidate)) else "候选" if candidate_level else "观察" if raw_supported else "正常"
|
|
|
+ if level in ("候选","持续预警") and not signature_match:level="未分类异常"
|
|
|
+ base_days=min(row["base"]["days"] for row in rows)
|
|
|
+ detail={
|
|
|
+ "事件":event["id"],"机组":event["unit"],"事件日期":event["event"],"日期":day,
|
|
|
+ "测点":point,"角度箱":angle,"等级":level,
|
|
|
+ "支持特征数":len(supported),"支持特征族数":len(families),
|
|
|
+ "幅度平均Z":round(mean_amp,4),"幅度中位Z":round(float(np.median(amplitudes)) if amplitudes else 0,4),
|
|
|
+ "综合RMS":round(rms,4),"异常6小时子窗数":abnormal_slots,
|
|
|
+ "最近5窗候选数":sum(hist),"基线有效日数":base_days,
|
|
|
+ "慢通道支持特征数":len(slow_supported_rows),"慢通道特征族数":len(slow_families),
|
|
|
+ "慢通道幅度平均Z":round(slow_mean,4),"慢通道综合RMS":round(slow_rms,4),
|
|
|
+ "慢通道最近5窗候选数":sum(slow_hist),"慢通道持续":int(slow_persistent),
|
|
|
+ "历史事件签名匹配":int(signature_match),
|
|
|
+ "基线置信":"低" if base_days<7 else "中" if base_days<20 else "高",
|
|
|
+ "特征摘要":";".join(f"{METRIC_CN[row['metric']]}{'↑' if row['local']>0 else '↓'}{abs(row['local']):.1f}Z" for row in sorted(supported,key=lambda x:-abs(x['local']))),
|
|
|
+ "慢通道特征摘要":";".join(f"{METRIC_CN[row['metric']]}{'↑' if row['slow_local']>0 else '↓'}{abs(row['slow_local']):.1f}Z" for row in sorted(slow_supported_rows,key=lambda x:-abs(x['slow_local']))),
|
|
|
+ "备注":event["note"],
|
|
|
+ }
|
|
|
+ all_detail.append(detail)
|
|
|
+ if level not in ("正常","通道质控"):event_alerts.append(detail)
|
|
|
+ bands = merge_bands(event_alerts); all_bands.extend(bands)
|
|
|
+ event_date = datetime.fromisoformat(event["event"]).date()
|
|
|
+ grouped_bands = defaultdict(list)
|
|
|
+ for row in bands: grouped_bands[(row["测点"], row["角度带"])].append(row)
|
|
|
+ for (point, band), rows in grouped_bands.items():
|
|
|
+ rows.sort(key=lambda row: row["日期"])
|
|
|
+ candidates = [row for row in rows if row["等级"] in ("候选", "持续预警", "未分类异常")]
|
|
|
+ persistent = [row for row in rows if row["等级"] == "持续预警"]
|
|
|
+ if not candidates: continue
|
|
|
+ first = datetime.fromisoformat(candidates[0]["日期"]).date()
|
|
|
+ validations.append({
|
|
|
+ "事件": event["id"], "事件日期": event["event"], "测点": point, "角度带": band,
|
|
|
+ "首次候选日期": str(first), "提前日数": (event_date-first).days,
|
|
|
+ "首次持续预警日期": persistent[0]["日期"] if persistent else "",
|
|
|
+ "候选窗口数": len(candidates), "持续预警窗口数": len(persistent),
|
|
|
+ "未分类异常窗口数": sum(row["等级"] == "未分类异常" for row in candidates),
|
|
|
+ "最大幅度平均Z": max(float(row["最大幅度平均Z"]) for row in candidates),
|
|
|
+ "最大综合RMS": max(float(row["最大综合RMS"]) for row in candidates),
|
|
|
+ "备注": event["note"],
|
|
|
+ })
|
|
|
+ if all_detail and any(row["事件"]==event["id"] for row in all_detail):
|
|
|
+ event_detail=[row for row in all_detail if row["事件"]==event["id"]]
|
|
|
+ event_band=[row for row in all_bands if row["事件"]==event["id"]]
|
|
|
+ event_validation=[row for row in validations if row["事件"]==event["id"]]
|
|
|
+ for filename,records in (("伪在线_窗口检测明细.csv",event_detail),
|
|
|
+ ("伪在线_角度带告警.csv",event_band),
|
|
|
+ ("伪在线_事件验证.csv",event_validation)):
|
|
|
+ if records:
|
|
|
+ with (reuse_dir/filename).open("w",encoding="utf-8-sig",newline="") as handle:
|
|
|
+ writer=csv.DictWriter(handle,fieldnames=list(records[0]));writer.writeheader();writer.writerows(records)
|
|
|
+ print(event["id"], "测点", len(point_data), "明细", len([r for r in all_detail if r['事件']==event['id']]),
|
|
|
+ "角度带", len(bands), flush=True)
|
|
|
+ finally:
|
|
|
+ connection.close()
|
|
|
+
|
|
|
+ out_dir.mkdir(parents=True, exist_ok=True)
|
|
|
+ def write(name, rows):
|
|
|
+ if not rows: return
|
|
|
+ with (out_dir/name).open("w", encoding="utf-8-sig", newline="") as handle:
|
|
|
+ writer=csv.DictWriter(handle,fieldnames=list(rows[0]));writer.writeheader();writer.writerows(rows)
|
|
|
+ write("伪在线_窗口检测明细.csv", all_detail)
|
|
|
+ write("伪在线_角度带告警.csv", all_bands)
|
|
|
+ write("伪在线_事件验证.csv", validations)
|
|
|
+
|
|
|
+ false_rows=[]
|
|
|
+ for event in EVENTS:
|
|
|
+ rows=[row for row in all_bands if row["事件"]==event["id"]]
|
|
|
+ total_days=len({row["日期"] for row in all_detail if row["事件"]==event["id"]})
|
|
|
+ for point in sorted({row["测点"] for row in rows}):
|
|
|
+ selected=[row for row in rows if row["测点"]==point]
|
|
|
+ false_rows.append({"事件":event["id"],"测点":point,"检测有效日数":total_days,
|
|
|
+ "观察或以上角度带窗口":len(selected),"候选或以上窗口":sum(row["等级"] in ('候选','持续预警') for row in selected),
|
|
|
+ "持续预警窗口":sum(row["等级"]=='持续预警' for row in selected)})
|
|
|
+ write("伪在线_误报统计.csv",false_rows)
|
|
|
+ return all_detail, all_bands, validations
|
|
|
+
|
|
|
+
|
|
|
+def run_health_segments(out_dir: Path, units=("7", "8", "9")):
|
|
|
+ """Replay unlabelled controls; triggers are not confirmed false positives."""
|
|
|
+ detail_rows=[]
|
|
|
+ connection=None
|
|
|
+ cache_dir=out_dir/"unlabelled_control_features_v1"
|
|
|
+ cache_dir.mkdir(parents=True,exist_ok=True)
|
|
|
+ keys=[(angle,metric) for angle in range(0,360,10) for metric in METRICS]
|
|
|
+ try:
|
|
|
+ for unit in units:
|
|
|
+ unit_detail_start=len(detail_rows)
|
|
|
+ for segment_start,segment_end in HEALTH_SEGMENTS.get(unit,[]):
|
|
|
+ point_data={}
|
|
|
+ for cylinder in CYLINDERS:
|
|
|
+ for side in SIDES:
|
|
|
+ if (cylinder,side) in BAD.get(unit,set()) or (unit=="8" and (cylinder,side)==("六","盖侧")):continue
|
|
|
+ cache_path=cache_dir/f"{unit}_{cylinder}_{side}_{segment_start}_{segment_end}.npz"
|
|
|
+ baseline=[];daily=[];slots=[];cycle_count=0
|
|
|
+ if cache_path.exists():
|
|
|
+ with np.load(cache_path,allow_pickle=False) as cached:
|
|
|
+ if list(cached["metrics"])!=list(METRICS):
|
|
|
+ raise ValueError(f"Feature schema mismatch: {cache_path}")
|
|
|
+ for kind,records in (("baseline",baseline),("daily",daily),("slots",slots)):
|
|
|
+ for day,slot,values in zip(cached[f"{kind}_days"],cached[f"{kind}_slots"],cached[f"{kind}_values"]):
|
|
|
+ records.append({"day":str(day),"slot":int(slot),"values":{key:float(value) for key,value in zip(keys,values) if np.isfinite(value)}})
|
|
|
+ cycle_count=int(cached["cycle_count"].item())
|
|
|
+ else:
|
|
|
+ if connection is None:connection=connect()
|
|
|
+ files=fetch_file_rows(connection,point_name(unit,cylinder,side),segment_start,segment_end)
|
|
|
+ by_day=defaultdict(list)
|
|
|
+ for row in files:
|
|
|
+ time=row["sample_time"] if isinstance(row["sample_time"],datetime) else datetime.fromisoformat(str(row["sample_time"]))
|
|
|
+ by_day[time.date().isoformat()].append({**row,"sample_time":time})
|
|
|
+ valid_days=0
|
|
|
+ for day,day_files in sorted(by_day.items()):
|
|
|
+ records=[]
|
|
|
+ day_files.sort(key=lambda row:(row["sample_time"],int(row["id"])))
|
|
|
+ for offset in range(0,len(day_files),100):
|
|
|
+ batch=day_files[offset:offset+100]
|
|
|
+ cycles=fetch_cycles(connection,[int(row["id"]) for row in batch],100)
|
|
|
+ for row in batch:
|
|
|
+ values=stats_for_cycle(cycles.get(int(row["id"]),np.array([])))
|
|
|
+ if values is None:continue
|
|
|
+ values={key:value for key,value in values.items() if key in keys and np.isfinite(value)}
|
|
|
+ if not values:continue
|
|
|
+ time=row["sample_time"]
|
|
|
+ records.append({"sample_time":time,"day":day,"slot":time.hour//6,"values":values})
|
|
|
+ del cycles
|
|
|
+ if not records:continue
|
|
|
+ valid_days+=1
|
|
|
+ if valid_days<=3:
|
|
|
+ baseline.extend(sample_per_day(records,3))
|
|
|
+ continue
|
|
|
+ cycle_count+=len(records)
|
|
|
+ daily_values=aggregate(records,"day")
|
|
|
+ daily.append({"day":day,"slot":-1,"values":{key:value for (group,key),value in daily_values.items()}})
|
|
|
+ slot_values=aggregate(records,"day_slot")
|
|
|
+ for slot in sorted({record["slot"] for record in records}):
|
|
|
+ slots.append({"day":day,"slot":slot,"values":{key:value for (group,key),value in slot_values.items() if group==f"{day}#{slot}"}})
|
|
|
+ payload={"metrics":np.asarray(METRICS),"cycle_count":np.asarray(cycle_count)}
|
|
|
+ for kind,records in (("baseline",baseline),("daily",daily),("slots",slots)):
|
|
|
+ payload[f"{kind}_days"]=np.asarray([record["day"] for record in records],dtype="U10")
|
|
|
+ payload[f"{kind}_slots"]=np.asarray([record["slot"] for record in records],dtype=int)
|
|
|
+ payload[f"{kind}_values"]=np.asarray([[record["values"].get(key,np.nan) for key in keys] for record in records],dtype=float).reshape(-1,len(keys))
|
|
|
+ temporary=cache_path.with_suffix(".tmp.npz")
|
|
|
+ np.savez_compressed(temporary,**payload)
|
|
|
+ temporary.replace(cache_path)
|
|
|
+ if baseline and daily:
|
|
|
+ point_data[f"{cylinder}缸{side}"]={"baseline":build_baseline(baseline),"daily":daily,"slots":slots,"cycle_count":cycle_count}
|
|
|
+ if not point_data:continue
|
|
|
+ z_rows=[];slot_z=defaultdict(dict)
|
|
|
+ for point,data in point_data.items():
|
|
|
+ for record in data["daily"]:
|
|
|
+ for key,value in record["values"].items():
|
|
|
+ base=data["baseline"][key]
|
|
|
+ if not np.isfinite(base["center"]):continue
|
|
|
+ z=(value-base["center"])/base["scale"]
|
|
|
+ z_rows.append({"point":point,"day":record["day"],"angle":key[0],"metric":key[1],"z":z,"baseline":base})
|
|
|
+ for record in data["slots"]:
|
|
|
+ for key,value in record["values"].items():
|
|
|
+ base=data["baseline"][key]
|
|
|
+ if not np.isfinite(base["center"]):continue
|
|
|
+ slot_z[(record["day"],record["slot"],key)][point]=(value-base["center"])/base["scale"]
|
|
|
+ public=defaultdict(list)
|
|
|
+ for row in z_rows:public[(row["day"],row["angle"],row["metric"])].append(row["z"])
|
|
|
+ public={key:float(np.median(values)) for key,values in public.items()}
|
|
|
+ slot_public={key:float(np.median(list(values.values()))) for key,values in slot_z.items()}
|
|
|
+ grouped=defaultdict(list)
|
|
|
+ for row in z_rows:
|
|
|
+ row["local"]=row["z"]-public[(row["day"],row["angle"],row["metric"])]
|
|
|
+ grouped[(row["point"],row["day"],row["angle"])].append(row)
|
|
|
+ history=defaultdict(deque)
|
|
|
+ for (point,day,angle),rows in sorted(grouped.items(),key=lambda item:(item[0][1],item[0][0],item[0][2])):
|
|
|
+ supported=[row for row in rows if abs(row["local"])>=2]
|
|
|
+ families={FAMILY[row["metric"]] for row in supported}
|
|
|
+ amps=[abs(row["local"]) for row in supported]
|
|
|
+ avg=float(np.mean(amps)) if amps else 0.0
|
|
|
+ rms=float(np.sqrt(np.mean(np.square(amps)))) if amps else 0.0
|
|
|
+ abnormal_slots=0
|
|
|
+ for slot in range(4):
|
|
|
+ slot_vals=[]
|
|
|
+ for row in rows:
|
|
|
+ key=(day,slot,(angle,row["metric"]))
|
|
|
+ if point in slot_z.get(key,{}):slot_vals.append(abs(slot_z[key][point]-slot_public[key]))
|
|
|
+ if sum(v>=2 for v in slot_vals)>=2:abnormal_slots+=1
|
|
|
+ candidate=len(families)>=2 and avg>=2 and rms>=3 and abnormal_slots>=2
|
|
|
+ hist=history[(point,angle)];hist.append(candidate)
|
|
|
+ if len(hist)>5:hist.popleft()
|
|
|
+ persistent=len(hist)>=3 and sum(hist)>=3
|
|
|
+ level="持续预警" if persistent and candidate else "候选" if candidate else "观察" if supported else "正常"
|
|
|
+ base_days=min(row["baseline"]["days"] for row in rows)
|
|
|
+ detail_rows.append({
|
|
|
+ "机组":unit,"运行段":f"{segment_start}~{segment_end}","日期":day,"测点":point,
|
|
|
+ "样本标签":"unlabelled control",
|
|
|
+ "备注":"2026-04-18距2026-04-01维修17天,不保证健康" if unit=="8" else "未确认健康,触发不能直接认定为误报",
|
|
|
+ "角度箱":angle,"等级":level,"支持特征数":len(supported),"支持特征族数":len(families),
|
|
|
+ "幅度平均Z":round(avg,4),"幅度中位Z":round(float(np.median(amps)) if amps else 0,4),
|
|
|
+ "综合RMS":round(rms,4),"异常6小时子窗数":abnormal_slots,"最近5窗候选数":sum(hist),
|
|
|
+ "基线有效日数":base_days,"基线置信":"低" if base_days<7 else "中" if base_days<20 else "高",
|
|
|
+ "特征摘要":";".join(f"{METRIC_CN[r['metric']]}{'↑' if r['local']>0 else '↓'}{abs(r['local']):.1f}Z" for r in sorted(supported,key=lambda x:-abs(x['local']))),
|
|
|
+ })
|
|
|
+ print(f"unlabelled control {unit}号 {segment_start}~{segment_end}: 测点{len(point_data)},全量检测周期{sum(x['cycle_count'] for x in point_data.values()):,}",flush=True)
|
|
|
+ unit_rows=detail_rows[unit_detail_start:]
|
|
|
+ if unit_rows:
|
|
|
+ out_dir.mkdir(parents=True,exist_ok=True)
|
|
|
+ with (out_dir/f"伪在线_未标注对照触发明细_{unit}号.csv").open("w",encoding="utf-8-sig",newline="") as handle:
|
|
|
+ writer=csv.DictWriter(handle,fieldnames=list(unit_rows[0]));writer.writeheader();writer.writerows(unit_rows)
|
|
|
+ finally:
|
|
|
+ if connection is not None:connection.close()
|
|
|
+ out_dir.mkdir(parents=True,exist_ok=True)
|
|
|
+ with (out_dir/"伪在线_未标注对照触发明细.csv").open("w",encoding="utf-8-sig",newline="") as handle:
|
|
|
+ writer=csv.DictWriter(handle,fieldnames=list(detail_rows[0]) if detail_rows else ["机组","样本标签"]);writer.writeheader();writer.writerows(detail_rows)
|
|
|
+ summary=[]
|
|
|
+ by_unit=defaultdict(list)
|
|
|
+ for row in detail_rows:by_unit[row["机组"]].append(row)
|
|
|
+ for unit,rows in by_unit.items():
|
|
|
+ days=len({row["日期"] for row in rows});points=len({row["测点"] for row in rows})
|
|
|
+ summary.append({"机组":unit,"样本标签":"unlabelled control","对照检测有效日数":days,"测点数":points,"窗口点箱记录":len(rows),
|
|
|
+ "观察记录":sum(row["等级"]=="观察" for row in rows),
|
|
|
+ "候选记录":sum(row["等级"]=="候选" for row in rows),
|
|
|
+ "持续预警记录":sum(row["等级"]=="持续预警" for row in rows),
|
|
|
+ "候选及以上日测点箱数":len({(row["日期"],row["测点"],row["角度箱"]) for row in rows if row["等级"] in ("候选","持续预警")}),
|
|
|
+ "持续预警日测点箱数":len({(row["日期"],row["测点"],row["角度箱"]) for row in rows if row["等级"]=="持续预警"})})
|
|
|
+ with (out_dir/"伪在线_未标注对照触发统计.csv").open("w",encoding="utf-8-sig",newline="") as handle:
|
|
|
+ writer=csv.DictWriter(handle,fieldnames=list(summary[0]) if summary else ["机组","样本标签"]);writer.writeheader();writer.writerows(summary)
|
|
|
+ return detail_rows,summary
|
|
|
+
|
|
|
+
|
|
|
+def report(out_dir, detail, bands, validations):
|
|
|
+ lines=["# 伪在线故障检测验证", "", "## 方法", "",
|
|
|
+ "- 历史基线每天粗采最多3周期;24小时检测窗口处理全部有效周期。",
|
|
|
+ "- 周期平台归一化、测点自身历史稳健Z标准化、单机公共变化扣除。",
|
|
|
+ "- 检测不使用维修后数据;维修日期仅用于计算提前量。", ""]
|
|
|
+ for event in EVENTS:
|
|
|
+ rows=[row for row in validations if row["事件"]==event["id"]]
|
|
|
+ rows.sort(key=lambda row:(-int(row["持续预警窗口数"]),-float(row["最大综合RMS"])))
|
|
|
+ lines += [f"## {event['id']}", "", f"- {event['note'] or '按事件前数据回放'}。", f"- 形成候选的测点角度带:{len(rows)}。"]
|
|
|
+ for i,row in enumerate(rows[:20],1):
|
|
|
+ lines.append(f"{i}. {row['测点']} {row['角度带']}:首次候选{row['首次候选日期']},提前{row['提前日数']}天,持续预警{row['持续预警窗口数']}窗,最大RMS{float(row['最大综合RMS']):.1f}。")
|
|
|
+ lines.append("")
|
|
|
+ lines += ["## 使用限制", "", "- 当前按有数据的自然日形成0~24小时窗口;下一版可改为每小时滚动24小时。",
|
|
|
+ "- 6号基线只有3个日期,所有结果应降级解释。", "- 7号2025事件附近缺测,只能验证运行段状态变化,不能证明事件提前量。"]
|
|
|
+ (out_dir/"伪在线检测结论.md").write_text("\n".join(lines)+"\n",encoding="utf-8")
|
|
|
+
|
|
|
+
|
|
|
+if __name__ == "__main__":
|
|
|
+ parser=argparse.ArgumentParser(description="局部特征伪在线回放")
|
|
|
+ parser.add_argument("--out-dir",default="cache/indicator/伪在线检测")
|
|
|
+ parser.add_argument("--event",action="append",choices=[event["id"] for event in EVENTS])
|
|
|
+ parser.add_argument("--healthy",action="store_true",help="回放无已知故障的健康运行段")
|
|
|
+ parser.add_argument("--unit",action="append",choices=UNITS,help="限制机组,可重复传入")
|
|
|
+ args=parser.parse_args()
|
|
|
+ output=Path(args.out_dir)
|
|
|
+ if args.healthy:
|
|
|
+ run_health_segments(output,tuple(args.unit or ("7","8","9")))
|
|
|
+ else:
|
|
|
+ detail,bands,validations=run(output,set(args.event or []))
|
|
|
+ report(output,detail,bands,validations)
|
|
|
+ print("结果 ->",output)
|