feat(burst-detection): update scada analysis flow

app/services/burst_detection.py

@@ -1,8 +1,10 @@
 from __future__ import annotations
 
-from datetime import datetime
+from collections import Counter
+from datetime import datetime, timedelta
 from typing import Any
 
+import numpy as np
 import pandas as pd
 
 from app.algorithms.burst_detection.burst_detector import BurstDetector
@@ -17,6 +19,15 @@ from app.services.tjnetwork import get_all_scada_info
 from app.services.time_api import extract_date, parse_utc_time, utc_now
 
 
+TARGET_DAY_COUNT = 15
+DEFAULT_SAMPLE_INTERVAL_MINUTES = 15
+TARGET_MU = 1
+TARGET_N_ESTIMATORS = 50
+TARGET_RANDOM_STATE = 42
+TARGET_SCORE_THRESHOLD = -0.04
+MIN_COMPLETE_SENSORS = 5
+
+
 def run_burst_detection(
     *,
     network: str,
@@ -31,6 +42,8 @@ def run_burst_detection(
     points_per_day: int = 1440,
     mu: int = 100,
     iforest_params: dict[str, Any] | None = None,
+    target_time: datetime | str | None = None,
+    sampling_interval_minutes: int | None = None,
     scada_start: datetime | str | None = None,
     scada_end: datetime | str | None = None,
     sensor_nodes: list[str] | None = None,
@@ -42,7 +55,8 @@ def run_burst_detection(
     """
     运行爆管侦测服务入口。
 
-    调用方式二选一：
+    调用方式三选一：
+    - 不传数据时间窗，自动侦测最近完整时刻；可用 `target_time` 回放历史时刻
     - 直接传 `observed_pressure_data`
     - 或传 `scada_start/scada_end` 让后端自动查询 SCADA 压力数据
 
@@ -74,8 +88,65 @@ def run_burst_detection(
         else None
     )
     use_scada_source = scada_start is not None or scada_end is not None
+    use_target_mode = (
+        observed_pressure_data is None
+        and not use_scada_source
+        and data_source != "simulation"
+    ) or target_time is not None
 
-    if use_scada_source:
+    resolved_target_time: datetime | None = None
+    requested_target_time: datetime | None = None
+    excluded_sensors: list[dict[str, str]] = []
+    daily_times: list[datetime] | None = None
+    resolved_sampling_interval_minutes: int | None = None
+
+    if use_target_mode:
+        if observed_pressure_data is not None or use_scada_source:
+            raise ValueError(
+                "target_time 不能与 observed_pressure_data 或 scada_start/scada_end 同时使用。"
+            )
+        scada_sensor_nodes = (
+            selected_sensor_nodes
+            if selected_sensor_nodes is not None
+            else _get_pressure_sensor_nodes(network)
+        )
+        requested_target_time = (
+            _to_datetime(target_time) if target_time is not None else None
+        )
+        resolved_sampling_interval_minutes = _resolve_sampling_interval_minutes(
+            network=network,
+            sensor_nodes=scada_sensor_nodes,
+            requested_interval=sampling_interval_minutes,
+        )
+        target_points_per_day = 1440 // resolved_sampling_interval_minutes
+        (
+            observed_input,
+            resolved_target_time,
+            excluded_sensors,
+        ) = _build_target_pressure_from_scada(
+            network=network,
+            sensor_nodes=scada_sensor_nodes,
+            requested_target_time=requested_target_time,
+            sampling_interval_minutes=resolved_sampling_interval_minutes,
+            points_per_day=target_points_per_day,
+        )
+        selected_sensor_nodes = list(observed_input.columns)
+        observed_source = (
+            "latest_monitoring" if target_time is None else "historical_monitoring"
+        )
+        points_per_day = target_points_per_day
+        mu = TARGET_MU
+        iforest_params = {
+            "n_estimators": TARGET_N_ESTIMATORS,
+            "random_state": TARGET_RANDOM_STATE,
+            "contamination": "auto",
+        }
+        daily_times = [
+            resolved_target_time - timedelta(days=offset)
+            for offset in range(TARGET_DAY_COUNT - 1, -1, -1)
+        ]
+
+    elif use_scada_source:
         scada_sensor_nodes = (
             selected_sensor_nodes
             if selected_sensor_nodes is not None
@@ -121,7 +192,16 @@ def run_burst_detection(
         sensor_nodes=selected_sensor_nodes,
     )
     resolved_sensor_nodes = list(result_df.attrs.get("sensor_nodes", []))
-    rows = _serialize_result_rows(result_df)
+    rows = _serialize_result_rows(
+        result_df,
+        daily_times=daily_times,
+        target_only=use_target_mode,
+    )
+    summary = _build_detection_summary(
+        result_df,
+        daily_times=daily_times,
+        target_only=use_target_mode,
+    )
     payload: dict[str, Any] = {
         "network": network,
         "sensor_nodes": resolved_sensor_nodes,
@@ -130,7 +210,17 @@ def run_burst_detection(
         "points_per_day": int(result_df.attrs.get("points_per_day", points_per_day)),
         "day_count": int(result_df.attrs.get("day_count", len(result_df))),
         "rows": rows,
-        "summary": _build_detection_summary(result_df),
+        "summary": summary,
+        "algorithm_params": {
+            "mu": mu,
+            "points_per_day": points_per_day,
+            "iforest_params": detector.iforest_params,
+            **(
+                {"score_threshold": TARGET_SCORE_THRESHOLD}
+                if use_target_mode
+                else {}
+            ),
+        },
     }
     if data_source == "simulation":
         payload["data_source"] = "simulation"
@@ -141,7 +231,50 @@ def run_burst_detection(
     else:
         payload["data_source"] = "monitoring"
 
-    if use_scada_source:
+    if (
+        use_target_mode
+        and resolved_target_time is not None
+        and resolved_sampling_interval_minutes is not None
+    ):
+        sample_start = resolved_target_time - timedelta(
+            days=TARGET_DAY_COUNT,
+            minutes=-resolved_sampling_interval_minutes,
+        )
+        payload.update(
+            {
+                "requested_target_time": (
+                    requested_target_time.isoformat()
+                    if requested_target_time is not None
+                    else None
+                ),
+                "target_time": resolved_target_time.isoformat(),
+                "reference_window": {
+                    "start": (resolved_target_time - timedelta(days=14)).isoformat(),
+                    "end": (resolved_target_time - timedelta(days=1)).isoformat(),
+                    "day_count": 14,
+                },
+                "sampling_interval_minutes": resolved_sampling_interval_minutes,
+                "daily_scores": [
+                    {
+                        "timestamp": row["Timestamp"],
+                        "role": row["Role"],
+                        "score": row["Score"],
+                        "raw_prediction": row["Prediction"],
+                    }
+                    for row in rows
+                ],
+                "data_quality": {
+                    "included_sensors": resolved_sensor_nodes,
+                    "excluded_sensors": excluded_sensors,
+                    "minimum_required_sensors": MIN_COMPLETE_SENSORS,
+                },
+                "scada_window": {
+                    "start": sample_start.isoformat(),
+                    "end": resolved_target_time.isoformat(),
+                },
+            }
+        )
+    elif use_scada_source:
         payload["scada_window"] = {
             "start": _to_datetime(scada_start).isoformat(),
             "end": _to_datetime(scada_end).isoformat(),
@@ -294,22 +427,49 @@ def _store_burst_detection_scheme(
     )
 
 
-def _serialize_result_rows(result_df: pd.DataFrame) -> list[dict[str, Any]]:
+def _serialize_result_rows(
+    result_df: pd.DataFrame,
+    *,
+    daily_times: list[datetime] | None = None,
+    target_only: bool = False,
+) -> list[dict[str, Any]]:
     rows: list[dict[str, Any]] = []
-    for row in result_df.to_dict(orient="records"):
+    raw_rows = result_df.to_dict(orient="records")
+    for index, row in enumerate(raw_rows):
+        is_target = index == len(raw_rows) - 1
+        is_burst = bool(row["IsBurst"])
+        if target_only:
+            is_burst = is_target and float(row["Score"]) <= TARGET_SCORE_THRESHOLD
         rows.append(
             {
                 "Day": int(row["Day"]),
                 "Score": float(row["Score"]),
                 "Prediction": int(row["Prediction"]),
-                "IsBurst": bool(row["IsBurst"]),
+                "IsBurst": is_burst,
+                **(
+                    {
+                        "Timestamp": daily_times[index].isoformat(),
+                        "Role": "target" if is_target else "reference",
+                    }
+                    if daily_times is not None
+                    else {}
+                ),
             }
         )
     return rows
 
 
-def _build_detection_summary(result_df: pd.DataFrame) -> dict[str, Any]:
-    rows = _serialize_result_rows(result_df)
+def _build_detection_summary(
+    result_df: pd.DataFrame,
+    *,
+    daily_times: list[datetime] | None = None,
+    target_only: bool = False,
+) -> dict[str, Any]:
+    rows = _serialize_result_rows(
+        result_df,
+        daily_times=daily_times,
+        target_only=target_only,
+    )
     if not rows:
         raise ValueError("爆管侦测结果为空。")
 
@@ -318,7 +478,7 @@ def _build_detection_summary(result_df: pd.DataFrame) -> dict[str, Any]:
     latest_row = rows[-1]
     anomaly_days = [row["Day"] for row in rows if row["IsBurst"]]
 
-    return {
+    summary = {
         "burst_detected": bool(latest_row["IsBurst"]),
         "latest_day": latest_row,
         "most_anomalous_day": int(result_df.iloc[most_anomalous_index]["Day"]),
@@ -326,6 +486,18 @@ def _build_detection_summary(result_df: pd.DataFrame) -> dict[str, Any]:
         "anomaly_day_count": len(anomaly_days),
         "latest_sensor_rankings": _build_latest_sensor_rankings(result_df),
     }
+    if target_only:
+        target_score = float(latest_row["Score"])
+        summary.update(
+            {
+                "target_score": target_score,
+                "score_threshold": TARGET_SCORE_THRESHOLD,
+                "target_rank": int(result_df["Score"].rank(method="min").iloc[-1]),
+                "target_time": latest_row.get("Timestamp"),
+                "reference_day_count": TARGET_DAY_COUNT - 1,
+            }
+        )
+    return summary
 
 
 def _build_latest_sensor_rankings(result_df: pd.DataFrame) -> list[dict[str, Any]]:
@@ -334,20 +506,194 @@ def _build_latest_sensor_rankings(result_df: pd.DataFrame) -> list[dict[str, Any
     if feature_matrix is None or len(sensor_nodes) == 0:
         return []
 
-    latest_values = feature_matrix[-1]
+    latest_values = np.asarray(feature_matrix[-1], dtype=float)
+    history = np.asarray(feature_matrix[:-1], dtype=float)
+    history_means = history.mean(axis=0)
+    history_stds = history.std(axis=0)
+    safe_stds = np.where(history_stds > 1e-9, history_stds, 1e-9)
+    deviations = (latest_values - history_means) / safe_stds
     ranking = sorted(
-        zip(sensor_nodes, latest_values, strict=False),
-        key=lambda item: item[1],
+        zip(
+            sensor_nodes,
+            latest_values,
+            history_means,
+            history_stds,
+            deviations,
+            strict=False,
+        ),
+        key=lambda item: item[4],
     )
     return [
         {
             "sensor_node": sensor_id,
             "latest_high_frequency_value": float(value),
+            "historical_mean": float(history_mean),
+            "historical_std": float(history_std),
+            "standardized_deviation": float(deviation),
         }
-        for sensor_id, value in ranking[: min(10, len(ranking))]
+        for sensor_id, value, history_mean, history_std, deviation in ranking[
+            : min(10, len(ranking))
+        ]
     ]
 
 
+def _build_target_pressure_from_scada(
+    *,
+    network: str,
+    sensor_nodes: list[str],
+    requested_target_time: datetime | None,
+    sampling_interval_minutes: int,
+    points_per_day: int,
+) -> tuple[pd.DataFrame, datetime, list[dict[str, str]]]:
+    node_query_id = _get_pressure_sensor_mapping(network)
+    mapped_nodes = [node for node in sensor_nodes if node in node_query_id]
+    excluded_without_mapping = [
+        {"sensor_node": node, "reason": "missing_api_query_id"}
+        for node in sensor_nodes
+        if node not in node_query_id
+    ]
+    if len(mapped_nodes) < MIN_COMPLETE_SENSORS:
+        raise ValueError(
+            f"可查询的压力测点少于 {MIN_COMPLETE_SENSORS} 个，无法执行爆管侦测。"
+        )
+
+    query_ids = [node_query_id[node] for node in mapped_nodes]
+    candidate_before = requested_target_time
+    last_excluded: list[dict[str, str]] = excluded_without_mapping
+
+    for _ in range(4):
+        resolved_target = InternalQueries.query_latest_scada_time(
+            db_name=network,
+            device_ids=query_ids,
+            before_time=candidate_before,
+        )
+        if resolved_target is None:
+            break
+
+        sample_start = resolved_target - timedelta(
+            days=TARGET_DAY_COUNT,
+            minutes=-sampling_interval_minutes,
+        )
+        expected_index = pd.date_range(
+            start=sample_start,
+            end=resolved_target,
+            freq=f"{sampling_interval_minutes}min",
+        )
+        scada_data = InternalQueries.query_scada_by_ids_timerange(
+            db_name=network,
+            device_ids=query_ids,
+            start_time=sample_start,
+            end_time=resolved_target,
+        )
+
+        complete_columns: dict[str, pd.Series] = {}
+        excluded = list(excluded_without_mapping)
+        for node_id in mapped_nodes:
+            query_id = node_query_id[node_id]
+            records = scada_data.get(query_id, [])
+            if not records:
+                excluded.append({"sensor_node": node_id, "reason": "no_data"})
+                continue
+
+            record_frame = pd.DataFrame.from_records(records)
+            record_frame["time"] = pd.to_datetime(record_frame["time"], utc=True)
+            record_frame["value"] = pd.to_numeric(
+                record_frame["value"], errors="coerce"
+            )
+            series = (
+                record_frame.drop_duplicates(subset="time", keep="last")
+                .set_index("time")["value"]
+                .reindex(expected_index)
+            )
+            if len(series) != TARGET_DAY_COUNT * points_per_day:
+                excluded.append(
+                    {"sensor_node": node_id, "reason": "unexpected_sample_count"}
+                )
+                continue
+            if series.isna().any():
+                excluded.append(
+                    {"sensor_node": node_id, "reason": "missing_or_invalid_samples"}
+                )
+                continue
+            complete_columns[node_id] = series
+
+        if len(complete_columns) >= MIN_COMPLETE_SENSORS:
+            observation_df = pd.DataFrame(complete_columns, index=expected_index)
+            return observation_df, resolved_target, excluded
+
+        last_excluded = excluded
+        candidate_before = resolved_target - timedelta(microseconds=1)
+
+    excluded_preview = ", ".join(
+        item["sensor_node"] for item in last_excluded[:10]
+    )
+    raise ValueError(
+        f"最近数据中完整压力测点少于 {MIN_COMPLETE_SENSORS} 个；"
+        f"请检查 15 天数据完整性。排除测点: {excluded_preview or '无'}"
+    )
+
+
+def _resolve_sampling_interval_minutes(
+    *,
+    network: str,
+    sensor_nodes: list[str],
+    requested_interval: int | None,
+) -> int:
+    if requested_interval is not None:
+        interval = int(requested_interval)
+    else:
+        selected_nodes = set(sensor_nodes)
+        inferred_intervals = [
+            parsed
+            for item in get_all_scada_info(network)
+            if str(item.get("type", "")).lower() == "pressure"
+            and str(item.get("associated_element_id", "")) in selected_nodes
+            and (
+                parsed := _parse_sampling_interval_minutes(
+                    item.get("transmission_frequency")
+                )
+            )
+            is not None
+        ]
+        interval = (
+            Counter(inferred_intervals).most_common(1)[0][0]
+            if inferred_intervals
+            else DEFAULT_SAMPLE_INTERVAL_MINUTES
+        )
+
+    if interval <= 0 or 1440 % interval != 0:
+        raise ValueError("采样间隔必须是能整除 1440 分钟的正整数。")
+    return interval
+
+
+def _parse_sampling_interval_minutes(value: Any) -> int | None:
+    if value is None:
+        return None
+    if isinstance(value, (int, float)):
+        minutes = float(value)
+    else:
+        try:
+            minutes = pd.to_timedelta(str(value)).total_seconds() / 60
+        except (TypeError, ValueError):
+            return None
+    rounded = round(minutes)
+    if minutes <= 0 or abs(minutes - rounded) > 1e-6:
+        return None
+    return int(rounded)
+
+
+def _get_pressure_sensor_mapping(network: str) -> dict[str, str]:
+    node_query_id: dict[str, str] = {}
+    for item in get_all_scada_info(network):
+        if str(item.get("type", "")).lower() != "pressure":
+            continue
+        node_id = item.get("associated_element_id")
+        query_id = item.get("api_query_id")
+        if node_id and query_id is not None:
+            node_query_id[str(node_id)] = str(query_id)
+    return node_query_id
+
+
 def _get_pressure_sensor_nodes(network: str) -> list[str]:
     sensor_nodes: list[str] = []
     for item in get_all_scada_info(network):
@@ -377,19 +723,7 @@ def _build_observed_pressure_from_scada(
     if start_dt >= end_dt:
         raise ValueError("SCADA 时间窗非法：scada_start 必须早于 scada_end。")
 
-    node_query_id: dict[str, str] = {}
-    for item in get_all_scada_info(network):
-        if str(item.get("type", "")).lower() != "pressure":
-            continue
-        node_id = item.get("associated_element_id")
-        query_id = item.get("api_query_id")
-        if (
-            isinstance(node_id, str)
-            and node_id
-            and isinstance(query_id, str)
-            and query_id
-        ):
-            node_query_id[node_id] = query_id
+    node_query_id = _get_pressure_sensor_mapping(network)
 
     missing_nodes = [node_id for node_id in sensor_nodes if node_id not in node_query_id]
     if missing_nodes: