功图特征抽取脚本与观察记录（只读、抽样，先找特征）

backend/indicator_features/README.md

@@ -0,0 +1,33 @@
+# 功图特征抽取（indicator_features）
+
+从 `wave_file` + `wave_sample_one`（每文件一个**首周期**）抽取单周期功图特征，用于找气阀/活塞故障特征。**只读、抽样、不建模**。
+
+级配（现场确认）：**一级=2/4/6 缸、二级=3/5 缸、三级=1 缸**；气体从 2/4/6 进，1 缸压缩后出，压力逐级升高。
+
+## 用法
+
+```bash
+cd backend
+python indicator_features/extract_from_db.py \
+  --unit 7 --cylinders 二 四 六 --sides 盖侧 轴侧 \
+  --start 2026-05-01 --end 2026-05-22 --per-day 3 \
+  --output cache/indicator/7_1stage_valve_202605.csv
+```
+
+- `--per-day N`：每个测点每天最多取 N 个文件（均匀取样），控制拉取量。
+- 输出每周期一行的 CSV：吸/排气平台、压比、相位占比、斜率、拐角、指示功等。
+
+## 特征（当前）
+
+- 平台：`p_suction`、`p_discharge`、`p_span`、`pressure_ratio`
+- 相位：`suction_frac`、`discharge_frac`
+- 线型：`compression_slope`、`expansion_slope`、`corner_sharpness`
+- 功：`indicated_work`、`imep`（**注意**：需物理相位锚定，当前不可靠）
+
+## 已知问题
+
+1. **坏通道**：如 8 号机六缸盖侧、7 号机六缸盖侧/四缸盖侧在部分时段压力≈0，必须识别并排除。
+2. **指示功/效率**：单周期起始角与活塞相位未锚定，`∮P dV` 目前不可用；需用键相/相位对齐后再算。
+3. **数据覆盖**：功图只在设备运行日存在，很多故障发生在停机附近，故障前数据可能不足。
+
+见 `观察记录.md`。

backend/indicator_features/extract_from_db.py

@@ -0,0 +1,165 @@
+#!/usr/bin/env python3
+"""Extract per-cycle indicator-diagram features from the wave database (read-only).
+
+Fetches a *small, sampled* set of running cycles (wave_sample_one holds one first
+cycle per wave_file) for the chosen cylinder points and time range, computes the
+phase-level features in indicator_features.py, and writes one CSV row per cycle.
+
+Sampling keeps the query volume low: at most --per-day files per point per day,
+picked evenly across that day's running files.
+"""
+from __future__ import annotations
+
+import argparse
+import csv
+import sys
+from collections import defaultdict
+from datetime import datetime
+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
+from indicator_features.indicator_features import (  # noqa: E402
+    CYLINDER_STAGE,
+    cycle_features,
+)
+
+CYLINDERS = ("一", "二", "三", "四", "五", "六")
+SIDES = ("盖侧", "轴侧")
+
+
+def connect() -> pymysql.connections.Connection:
+    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: str, cylinder: str, side: str) -> str:
+    return f"{unit}号机组{cylinder}缸压力{side}"
+
+
+def fetch_file_rows(connection, point: str, start: str, end: str) -> list[dict]:
+    with connection.cursor() as cursor:
+        cursor.execute(
+            """
+            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
+            ORDER BY sample_time
+            """,
+            (point, start, end),
+        )
+        return list(cursor.fetchall())
+
+
+def sample_per_day(rows: list[dict], per_day: int) -> list[dict]:
+    if per_day <= 0:
+        return rows
+    by_day: dict = defaultdict(list)
+    for row in rows:
+        by_day[row["sample_time"].date()].append(row)
+    picked: list[dict] = []
+    for day in sorted(by_day):
+        day_rows = by_day[day]
+        if len(day_rows) <= per_day:
+            picked.extend(day_rows)
+            continue
+        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 fetch_cycles(connection, file_ids: list[int], batch: int) -> dict[int, np.ndarray]:
+    cycles: dict[int, list[tuple[int, float]]] = defaultdict(list)
+    with connection.cursor() as cursor:
+        for start in range(0, len(file_ids), batch):
+            chunk = file_ids[start:start + batch]
+            placeholders = ", ".join(["%s"] * len(chunk))
+            cursor.execute(
+                f"""
+                SELECT wave_file_id, sample_index, CAST(signal_value AS FLOAT) AS sig
+                FROM wave_sample_one
+                WHERE wave_file_id IN ({placeholders})
+                ORDER BY wave_file_id, sample_index
+                """,
+                tuple(chunk),
+            )
+            for row in cursor.fetchall():
+                if row["sig"] is None:
+                    continue
+                cycles[int(row["wave_file_id"])].append((int(row["sample_index"]), float(row["sig"])))
+    result = {}
+    for file_id, rows in cycles.items():
+        rows.sort()
+        result[file_id] = np.asarray([value for _, value in rows], dtype=float)
+    return result
+
+
+def main() -> None:
+    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=list(SIDES), choices=SIDES)
+    parser.add_argument("--start", required=True, help="YYYY-MM-DD")
+    parser.add_argument("--end", required=True, help="YYYY-MM-DD（不含）")
+    parser.add_argument("--per-day", type=int, default=3, help="每测点每天最多取几个文件")
+    parser.add_argument("--batch", type=int, default=200, help="IN 查询批大小")
+    parser.add_argument("--output", required=True)
+    args = parser.parse_args()
+
+    connection = connect()
+    rows_out: list[dict] = []
+    try:
+        for cylinder in args.cylinders:
+            for side in args.sides:
+                point = point_name(args.unit, cylinder, side)
+                files = fetch_file_rows(connection, point, args.start, args.end)
+                picked = sample_per_day(files, args.per_day)
+                print(f"{point}: 运行文件 {len(files)}，抽样 {len(picked)}", flush=True)
+                if not picked:
+                    continue
+                id_to_time = {int(r["id"]): r for r in picked}
+                cycles = fetch_cycles(connection, list(id_to_time), args.batch)
+                for file_id, signal in cycles.items():
+                    feats = cycle_features(signal, cylinder)
+                    if not feats:
+                        continue
+                    meta = id_to_time[file_id]
+                    rows_out.append({
+                        "unit": args.unit,
+                        "cylinder": cylinder,
+                        "stage": CYLINDER_STAGE[cylinder],
+                        "side": side,
+                        "point_name": point,
+                        "wave_file_id": file_id,
+                        "sample_time": str(meta["sample_time"]),
+                        "rpm": float(meta["rpm"] or 0.0),
+                        **feats,
+                    })
+    finally:
+        connection.close()
+
+    if not rows_out:
+        print("没有特征行，退出")
+        return
+    rows_out.sort(key=lambda r: (r["cylinder"], r["side"], r["sample_time"]))
+    output = Path(args.output)
+    output.parent.mkdir(parents=True, exist_ok=True)
+    with output.open("w", encoding="utf-8-sig", newline="") as handle:
+        writer = csv.DictWriter(handle, fieldnames=list(rows_out[0].keys()))
+        writer.writeheader()
+        writer.writerows(rows_out)
+    print(f"特征行 {len(rows_out):,} -> {output}")
+    print(f"时间：{datetime.now():%Y-%m-%d %H:%M:%S}")
+
+
+if __name__ == "__main__":
+    main()

backend/indicator_features/indicator_features.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+"""Phase-level features from one compressor indicator diagram cycle.
+
+One stored cycle (wave_sample_one) is the pressure of a single cylinder chamber
+over one crankshaft revolution (~60 ms at ~996 RPM).  Absolute pressure levels
+depend on the process, so these features are meant to be compared:
+  * against the same cylinder's own recent baseline (trend), and
+  * against the peer cylinders of the same compression stage at the same time,
+  * and between the cover/shaft chambers of the same cylinder.
+
+Compression staging (confirmed with the site):
+  一级(1st) = 2/4/6 缸,  二级(2nd) = 3/5 缸,  三级(3rd) = 1 缸.
+No model is trained here; the goal is to expose interpretable features.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+CYLINDER_BORE_MM = {"一": 360.0, "二": 490.0, "三": 390.0, "四": 490.0, "五": 390.0, "六": 490.0}
+CYLINDER_STAGE = {"一": "三级", "二": "一级", "三": "二级", "四": "一级", "五": "二级", "六": "一级"}
+PISTON_STROKE_MM = 148.0
+CONNECTING_ROD_LENGTH_MM = 460.0
+CLEARANCE_VOLUME_L_BY_BORE = {490.0: 1.51, 390.0: 0.74, 360.0: 0.62}
+
+
+def smooth(values: np.ndarray, window: int) -> np.ndarray:
+    window = max(1, int(window))
+    if window <= 1 or len(values) < window:
+        return values
+    kernel = np.ones(window, dtype=float) / window
+    return np.convolve(values, kernel, mode="same")
+
+
+def volume_from_angle(count: int, bore_mm: float) -> np.ndarray:
+    """Chamber volume (L) over one revolution, from crank-slider geometry."""
+    crank_radius = PISTON_STROKE_MM / 2.0
+    clearance = CLEARANCE_VOLUME_L_BY_BORE.get(bore_mm)
+    if clearance is None:
+        clearance = 1.0
+    theta = np.linspace(0.0, 2.0 * np.pi, count, endpoint=False)
+    travel = (
+        crank_radius * (1.0 - np.cos(theta))
+        + CONNECTING_ROD_LENGTH_MM
+        - np.sqrt(CONNECTING_ROD_LENGTH_MM**2 - (crank_radius * np.sin(theta)) ** 2)
+    )
+    piston_area = np.pi * (bore_mm / 2.0) ** 2
+    return clearance + piston_area * travel / 1_000_000.0
+
+
+def _plateau_levels(p: np.ndarray) -> tuple[float, float]:
+    """Suction (low) and discharge (high) plateau levels via robust quantiles."""
+    low = np.quantile(p, 0.12)
+    high = np.quantile(p, 0.88)
+    near_low = p[p <= low + 0.15 * (high - low)]
+    near_high = p[p >= high - 0.15 * (high - low)]
+    suction = float(np.median(near_low)) if near_low.size else float(np.min(p))
+    discharge = float(np.median(near_high)) if near_high.size else float(np.max(p))
+    return suction, discharge
+
+
+def cycle_features(signal: np.ndarray, cylinder: str, smooth_window: int = 21) -> dict:
+    signal = np.asarray(signal, dtype=float)
+    signal = signal[np.isfinite(signal)]
+    count = len(signal)
+    if count < 50:
+        return {}
+    p = smooth(signal, smooth_window)
+    p = np.where(np.isfinite(p), p, signal)
+    suction, discharge = _plateau_levels(p)
+    span = discharge - suction
+    threshold = suction + 0.5 * span
+    high = p >= threshold
+    dp = np.diff(p)
+    rising = dp[dp > 0]
+    falling = dp[dp < 0]
+    # Compression/expansion line steepness (per sample).
+    compression_slope = float(np.quantile(rising, 0.9)) if rising.size else 0.0
+    expansion_slope = float(-np.quantile(falling, 0.1)) if falling.size else 0.0
+    # Plateau occupancy: fraction of the cycle spent near high/low plateau.
+    near = 0.08 * span + 1e-9
+    suction_frac = float(np.mean(np.abs(p - suction) <= near))
+    discharge_frac = float(np.mean(np.abs(p - discharge) <= near))
+    # Indicated work over the revolution (net area) and mean effective pressure.
+    volume = volume_from_angle(count, CYLINDER_BORE_MM.get(cylinder, 490.0))
+    dvolume = np.diff(volume, append=volume[0])
+    work = float(-np.sum(p * dvolume))  # net work done on the gas (positive)
+    sweep = float(np.max(volume) - np.min(volume))
+    imep = work / sweep if sweep else 0.0
+    # Corner sharpness: mean |curvature| of pressure where it crosses the mid band.
+    d2 = np.diff(p, 2)
+    band = np.abs(p[1:-1] - threshold) <= 0.1 * span
+    corner = float(np.mean(np.abs(d2[band]))) if band.any() else 0.0
+    return {
+        "n_samples": count,
+        "p_min": float(np.min(p)),
+        "p_max": float(np.max(p)),
+        "p_mean": float(np.mean(p)),
+        "p_suction": suction,
+        "p_discharge": discharge,
+        "p_span": span,
+        "pressure_ratio": (discharge / suction) if suction else None,
+        "pressure_rise": span,
+        "suction_frac": suction_frac,
+        "discharge_frac": discharge_frac,
+        "compression_slope": compression_slope,
+        "expansion_slope": expansion_slope,
+        "indicated_work": work,
+        "imep": imep,
+        "corner_sharpness": corner,
+    }

backend/indicator_features/观察记录.md

@@ -0,0 +1,57 @@
+# 功图找气阀/活塞故障 —— 观察记录
+
+> 只读分析，样本为 `wave_file` + `wave_sample_one`（首周期）。故障锚点来自 `7,8,9号压缩机故障汇总.xlsx`。
+
+## 1. 功图数据方法（确认）
+
+- 采样：25.6 kHz、键相 8 脉冲/转、≈996 RPM → 单周期≈60 ms、≈1540 点。
+- `wave_sample_one` 现在每个文件就是**一个首周期**（DB 实测 1545 点）。
+- `wave_file.cycle_start/cycle_end` 已废弃，不使用。
+- 只认 `cache/pks_scanner/evolution/` 作为异常段来源。
+
+## 2. 数据覆盖（按月，rpm>0 运行文件）
+
+| 测点 | 运行覆盖 |
+|---|---|
+| 8 号机 全部测点 | 2025-03~05、2026-03~06 |
+| 7 号机 一级(2/4/6)、四级… | 2025-04~07、2025-11、2026-02~06（2025-08/10 无） |
+| 9 号机 一级(2缸) | 2025-03、06、09、10、2026-01~06 |
+| 7 号机 一/六缸、9 号机 一/五/六缸 | 2026-01 起才有 |
+
+**含义**：多数 2025 年气阀/活塞事件没有功图数据；可用的候选事件：
+
+- **7 号机 2026.5.21 气阀**：2026-05-06~05-20 连续密集运行（最优，故障前两周）
+- **8 号机 2026.3.29 / 4.1 气阀**：运行仅到 2026-03-20，随后停机
+- **9 号机 2026.1.18 活塞（三级=1 缸）**：1 缸数据从 2026-01-30 起（故障后才运行）
+
+## 3. 已验证的观察
+
+### 3.1 7 号机 一级 2026-05 故障前两周
+特征按天均值基本平稳，无显著趋势：
+
+| 缸·侧 | p_discharge 走势 | 压比 前段→末段 |
+|---|---|---|
+| 二·盖 | 1.79→1.91→1.85 | 2.93→2.99 (+2.1%) |
+| 二·轴 | 1.69→1.89→1.81 | 2.69→2.69 (−0.2%) |
+| 四·轴 | 1.77→1.92→1.82 | 2.91→3.01 (+3.5%) |
+| 六·轴 | 1.62→1.76→1.68 | 3.05→3.10 (+1.5%) |
+
+- 该事件是“进气阀组油压管崩开”，属**突发**，**故障前期无趋势符合预期**。
+- 相邻时刻功图几乎重合；跨天变化在 ±3% 内 → 与“压力变化微弱”的观察一致。
+
+### 3.2 问题
+- **坏通道**：8 号机六缸盖侧 Pd≈0.02；7 号机六缸盖侧 Pd≈0.8、四缸盖侧 NaN。做缸间对照前必须剔除。
+- **指示功/IMEP 不可用**：`∮P dV` 依赖周期起始与活塞相位的对齐，当前未锚定，出现负值/爆量。
+
+## 4. 下一步（找特征）
+
+1. **稳健平台检测**：用压力直方图双峰/众数定 P_s、P_d，替代分位数，稳住压比。
+2. **物理相位锚定**：用键相/触发把周期对齐到活塞相位，使 `∮P dV`、压力转换效率可用。
+3. **归一化形状特征**：每周期按相位重采样 P(θ) 并做“同级缸间 / 同缸双腔 / 本缸近期基线”对照。
+4. **只挑渐进型故障**：现有最佳是 7 号机 2026-05（突发），需再找带“故障前持续运行”的渐进事件验证。
+5. 坏通道自动剔除 + 每测点近期基线。
+
+## 5. 结论（当前）
+
+- 抽取链路可用，平台类特征稳定（同侧同级可比、双腔可比）。
+- 但在目前可用的故障前数据上**尚未找到明显趋势**；需先完善相位/形状特征，并找合适的渐进型事件。