Kaynağa Gözat

功图:可用性扫描(往前找数据段)、等熵效率、长期偏离检测与事件观察修正

18922397810 3 hafta önce
ebeveyn
işleme
3f4e0cb3aa

+ 112 - 0
backend/indicator_features/availability.py

@@ -0,0 +1,112 @@
+#!/usr/bin/env python3
+"""Scan where indicator-diagram cycles actually exist (read-only).
+
+For each point it samples a few files per day *without* an rpm filter (rpm is
+unreliable for wave_file) and checks whether a valid first cycle exists in
+wave_sample_one with a real pressure span.  This finds the last data segment
+before a fault, which may be weeks before the repair record.
+"""
+from __future__ import annotations
+
+import argparse
+import sys
+from collections import defaultdict
+from pathlib import Path
+
+import numpy as np
+import pymysql
+
+BACKEND = Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(BACKEND))
+
+from app.config import settings  # noqa: E402
+
+CYLINDERS = ("一", "二", "三", "四", "五", "六")
+SIDES = ("盖侧", "轴侧")
+MIN_SPAN = 0.3
+
+
+def connect():
+    return pymysql.connect(
+        host=settings.db_host, port=settings.db_port, user=settings.db_user,
+        password=settings.db_password, database=settings.db_name, charset="utf8mb4",
+        cursorclass=pymysql.cursors.DictCursor, connect_timeout=8, read_timeout=300,
+    )
+
+
+def point_name(unit, cylinder, side):
+    return f"{unit}号机组{cylinder}缸压力{side}"
+
+
+def sample_per_day(rows, per_day):
+    by_day = defaultdict(list)
+    for row in rows:
+        by_day[row["sample_time"].date()].append(row)
+    picked = []
+    for day in sorted(by_day):
+        day_rows = by_day[day]
+        if len(day_rows) <= per_day:
+            picked.extend(day_rows)
+        else:
+            idx = np.linspace(0, len(day_rows) - 1, per_day).round().astype(int)
+            picked.extend(day_rows[int(i)] for i in idx)
+    return picked
+
+
+def spans(connection, file_ids):
+    """Median pressure span of wave_sample_one per file_id."""
+    result = {}
+    with connection.cursor() as cursor:
+        for start in range(0, len(file_ids), 200):
+            chunk = file_ids[start:start + 200]
+            ph = ", ".join(["%s"] * len(chunk))
+            cursor.execute(
+                f"SELECT wave_file_id, COUNT(*) n, MIN(signal_value) mn, MAX(signal_value) mx "
+                f"FROM wave_sample_one WHERE wave_file_id IN ({ph}) GROUP BY wave_file_id",
+                tuple(chunk),
+            )
+            for row in cursor.fetchall():
+                mn, mx = row["mn"], row["mx"]
+                result[int(row["wave_file_id"])] = float(mx) - float(mn) if mn is not None and mx is not None else 0.0
+    return result
+
+
+def main():
+    parser = argparse.ArgumentParser(description="扫描功图周期实际存在的日期")
+    parser.add_argument("--unit", required=True, choices=("7", "8", "9"))
+    parser.add_argument("--cylinders", nargs="+", default=list(CYLINDERS), choices=CYLINDERS)
+    parser.add_argument("--sides", nargs="+", default=["盖侧"], choices=SIDES)
+    parser.add_argument("--start", required=True)
+    parser.add_argument("--end", required=True)
+    parser.add_argument("--per-day", type=int, default=2)
+    args = parser.parse_args()
+
+    connection = connect()
+    try:
+        for cylinder in args.cylinders:
+            for side in args.sides:
+                point = point_name(args.unit, cylinder, side)
+                with connection.cursor() as cursor:
+                    cursor.execute(
+                        "SELECT id, sample_time FROM wave_file WHERE point_name=%s AND measurement_type='压力' "
+                        "AND sample_time >= %s AND sample_time < %s ORDER BY sample_time",
+                        (point, args.start, args.end),
+                    )
+                    files = list(cursor.fetchall())
+                picked = sample_per_day(files, args.per_day)
+                span_by_id = spans(connection, [int(r["id"]) for r in picked])
+                by_day = defaultdict(list)
+                for row in picked:
+                    by_day[row["sample_time"].date()].append(span_by_id.get(int(row["id"]), 0.0))
+                valid_days = [d for d, v in sorted(by_day.items()) if np.median(v) >= MIN_SPAN and len(span_by_id)]
+                all_days = sorted(by_day)
+                print(f"{point}: 文件 {len(files)},抽样日 {len(all_days)},其中有效日 {len(valid_days)}")
+                if valid_days:
+                    print(f"   有效数据段: {valid_days[0]} ~ {valid_days[-1]}  共 {len(valid_days)} 天")
+                    print("   有效日:", ", ".join(str(d) for d in valid_days))
+    finally:
+        connection.close()
+
+
+if __name__ == "__main__":
+    main()

+ 1 - 1
backend/indicator_features/extract_from_db.py

@@ -52,7 +52,7 @@ def fetch_file_rows(connection, point: str, start: str, end: str) -> list[dict]:
             SELECT id, sample_time, rpm
             FROM wave_file
             WHERE point_name = %s AND measurement_type = '压力'
-              AND rpm > 0 AND sample_time >= %s AND sample_time < %s
+              AND sample_time >= %s AND sample_time < %s
             ORDER BY sample_time
             """,
             (point, start, end),

+ 1 - 1
backend/indicator_features/indicator_features.py

@@ -134,7 +134,7 @@ def cycle_features(signal: np.ndarray, cylinder: str, smooth_window: int = 21) -
         )
     else:
         ideal = None
-    efficiency = (indicated_work / ideal) if ideal and ideal > 1e-9 else None
+    efficiency = (ideal / indicated_work) if ideal and indicated_work and indicated_work > 1e-9 else None
 
     # Phase boundaries by pressure crossing (aligned).
     high_mask = p >= threshold

+ 104 - 0
backend/indicator_features/trend.py

@@ -0,0 +1,104 @@
+#!/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:
+            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 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}")
+
+
+if __name__ == "__main__":
+    main()

+ 51 - 34
backend/indicator_features/观察记录.md

@@ -1,47 +1,64 @@
 # 功图找气阀/活塞故障 —— 观察记录
 
-> 只读分析。样本为 `wave_file` + `wave_sample_one`(首周期,DB 实测约 1545 点)。故障锚点来自 `7,8,9号压缩机故障汇总.xlsx`;异常段只认 `cache/pks_scanner/evolution/`。
+> 只读分析。样本 `wave_file` + `wave_sample_one`(首周期,约 1545 点)。故障锚点:`7,8,9号压缩机故障汇总.xlsx`;异常段对齐 `cache/pks_scanner/evolution/`。
+> 目录 `backend/indicator_features/` 为**新增独立**目录,只读数据库、只 `import app.config`,不影响任何现有算法/程序。
 
-## 1. 功图数据方法(确认)
+## 1. 数据方法(确认)
 
-- 采样:25.6 kHz、键相 8 脉冲/转、≈996 RPM → 单周期≈60 ms、≈1540 点;15 分钟一个文件。
-- `wave_sample_one` 现在每文件就是**一个首周期**;`wave_file.cycle_start/cycle_end` 已废弃。
-- 级配:**一级=2/4/6、二级=3/5、三级=1**;压力逐级升高(一级 Pd≈1.4~2.2、二级≈3.2~4.4、三级≈5.0~5.1 MPa)。
+- 25.6 kHz、8 脉冲/转、≈996 RPM → 单周期≈60 ms;15 分钟一个文件。
+- `wave_sample_one` = 每文件**一个首周期**;`cycle_start/end` 废弃。
+- 级配:**一级=2/4/6、二级=3/5、三级=1**(Pd 约 1.4~2.2 / 3.2~4.4 / 5.0~5.1 MPa)。
 
-## 2. 已完成改进(对应任务 1/2/3)
+## 2. 可用性(往前找,直到有周期的日子)
 
-1. **稳健平台检测**:用 Otsu 双峰分割替代分位数,稳定求 `P_s`(吸气平台)、`P_d`(排气平台)、`pressure_ratio`。
-2. **物理相位对齐**:在整转内搜索使 PV 环面积最大的曲柄偏移,得到可用的 `indicated_work`、等熵理想功与 `efficiency`、多变指数 `n_comp`。
-3. **归一化相位形状/对照**(`analyze.py`):
-   - 同级缸间偏差(同侧、同刻相对中位数);
-   - 同缸双腔差(盖侧 − 轴侧,含归一化相位剖面 L2 距离);
-   - 本缸自身近期基线 z(沿时间趋势)。
-   - 通道 QC:按“相对同级期望排气压力”自动标记 dead/weak/high 并排除。
+`availability.py` 扫出“真正存在完整功图周期”的日子;`rpm` 不作为运行判据(不可靠)。
 
-## 3. 7 号机一级 2026-05(故障前两周,任务 4 的最佳样本)
+| 测点 | 有效段(有周期) |
+|---|---|
+| 8 号机三缸盖侧 | 2026-03-16~20、04-19~05-06、05-31~06-18 |
+| 7 号机二缸盖侧 | 2025-06-13、07-25~27、11-18~20、2026-02、03-13~04-04、05-07~21、06-01~06 |
+| 9 号机一缸盖侧 | 2026-01-30~02-09、03-04~13、03-22~27、04-02~18、05-22~31、06-15~22 |
 
-- 通道标记:**六缸盖侧 weak、四缸盖侧 dead**(已排除对照)。
-- 同级(轴侧)偏差稳定:二/四/六 p_discharge 偏差 −0.5% / +1.8% / −6.6%,pressure_ratio −6.3% / +1.9% / +3.7%。
-- 本缸自身趋势:前段 z 与后段 z 基本在 ±0.7 内,**无持续漂移**。
-- 结论:**故障前无功图趋势**,与该事件“进气阀组油压管突发崩开”的性质一致(突发而非渐进)。
+**据此修正后可用的事件**(往前找最后一段):
 
-## 4. 事件可用性(任务 4)
+| 事件 | 故障前最后有功图段 |
+|---|---|
+| 7 号机 2025.8.8 气阀 | 2025-07-25~27 ✅ |
+| 8 号机 2026.3.29 / 4.1 气阀 | 2026-03-16~20 ✅ |
+| 7 号机 2026.5.21 气阀 | 2026-05-07~21 ✅(连续) |
+| 8 号机 2025.12.6 活塞 | 无功图(2025-06 后未采集)❌ |
+| 9 号机 2026.1.18 活塞 | 1 缸 2026-01-30 才有(故障后)❌ |
+| 7 号机 2025.1.10 活塞 | 无功图 ❌ |
 
-| 类型 | 事件 | 故障前功图 |
-|---|---|---|
-| 气阀 | 7 号机 2026.5.21 | ✅ 2026-05-06~20 连续运行(但属突发) |
-| 气阀 | 8 号机 2026.3.29 / 4.1 | ❌ 8 号机 2026-03-20 后停机 |
-| 气阀 | 8 号机 2025.2.17 | ❌ 8 号机 2025-03 才有数据 |
-| 气阀 | 7 号机 2025.8.8 / 2025.10.7 | ❌ 无 2025-08 / 2025-10 数据 |
-| 活塞 | 8 号机 2025.12.6 | ❌ 故障前已停机 |
-| 活塞 | 9 号机 2026.1.18(三级=1缸) | ❌ 1 缸数据 2026-01-30 才有(故障后) |
-| 活塞 | 7 号机 2025.1.10 | ❌ 早于数据范围 |
+## 3. 特征(任务 1/2 完成)
 
-**结论:当前功图数据里没有一个“故障前持续运行 + 渐进退化”的气阀/活塞事件**,因此暂无法证明或否定功图趋势。唯一可用事件是突发型。
+- 稳健平台:Otsu 双峰求 `P_s / P_d / pressure_ratio`。
+- 相位对齐:整转内取 PV 环面积最大的曲柄偏移 → `indicated_work`、`ideal_work_isentropic`、**`efficiency`(等熵效率,实测 0.68~0.76,合理)**、`n_comp`。
+- `analyze.py`:同级缸间、同缸双腔(含相位剖面 L2)、本缸自身 z,含通道 QC(dead/weak/high 自动排除)。
+- `trend.py`:按天自身基线 z,并与 evolution 压力/振动异常段对齐(任务 3)。
 
-## 5. 下一步
+## 4. 关键观察
 
-1. 等新数据积累(2026-06 后有 4/6 号机与 7/8/9 号机的密集数据),找渐进型事件;
-2. 校准等熵指数 k(当前 efficiency≈1.2,偏高),把“压力转换效率”做成可解释指标;
-3. 用 `n_comp`、相位剖面 L2、双腔不对称做长期基线,做“偏离检测”而非分类;
-4. 与 `cache/pks_scanner/evolution/` 的压力/振动异常段做时间对齐,多证据互验。
+### 4.1 8 号机 2026 春季(含 3.29/4.1 气阀)
+- 3-16~20 一级排气压力比 4/5 月低约 10 个 MAD(z≈−10),03-18 落在 `压力/轻微` 段内。
+- 但**同级缸间对照**显示:二/四盖侧仅 ∓4%、各缸与盖/轴两腔**同步移动** → 这是**工况/负荷级变化,不是单缸故障**。
+- 说明:必须用同级/双腔对照消工况;单看自身 z 会误报。
+
+### 4.2 7 号机 2026-05 气阀前两周
+- 通道:六缸盖侧 weak、四缸盖侧 dead。
+- 同级轴侧偏差稳定(Pd ±0.5~6.6%、压比 ±2~6%);无趋势(该事件为管道突发)。
+
+### 4.3 问题
+- `n_comp` 在通道间不稳定(部分为 0/负),暂不宜作为主特征。
+- 坏/弱通道常见,必须 QC。
+
+## 5. 结论
+
+- 抽取链路、稳健平台、等熵效率、三种归一化、长期偏离检测均已就绪。
+- 现有数据里:可用事件都是**气阀**且多为突发;**尚未发现单缸级、故障前可持续观测的功图趋势**。
+- 8 号机 3 月的大偏离是工况级,说明对照归一化是必需的。
+
+## 6. 下一步
+1. 等 2026-06 后新事件(4/6 号机与 7/8/9 密集段);
+2. 稳定 `n_comp`(限制拟合区间/最小样本);
+3. 用“同级 + 双腔”偏差做长期基线与偏离报警,与 evolution 段多证据互验。