From db6032bd841cc5c2ac83b110e17147c5240aae38 Mon Sep 17 00:00:00 2001 From: Jiang Date: Fri, 17 Jul 2026 16:31:46 +0800 Subject: [PATCH] feat(burst-detection): update scada analysis flow --- app/api/v1/endpoints/burst_detection.py | 10 + app/infra/db/timescaledb/internal_queries.py | 32 ++ .../db/timescaledb/repositories/scada.py | 21 + app/services/burst_detection.py | 392 ++++++++++++++++-- tests/unit/test_burst_detection_service.py | 160 +++++++ 5 files changed, 586 insertions(+), 29 deletions(-) create mode 100644 tests/unit/test_burst_detection_service.py diff --git a/app/api/v1/endpoints/burst_detection.py b/app/api/v1/endpoints/burst_detection.py index 37d7849..e1cb237 100644 --- a/app/api/v1/endpoints/burst_detection.py +++ b/app/api/v1/endpoints/burst_detection.py @@ -30,6 +30,16 @@ class BurstDetectionRequest(BaseModel): points_per_day: int = Field(1440, description="每天的数据点数") mu: int = Field(100, description="异常值检测的参数") iforest_params: dict[str, Any] | None = Field(None, description="隔离森林算法参数") + target_time: datetime | None = Field( + None, + description="目标侦测时刻；为空时自动使用最近一个完整的监测时刻", + ) + sampling_interval_minutes: int | None = Field( + None, + ge=1, + le=1440, + description="采样间隔（分钟）；为空时根据压力 SCADA 传输频率自动推断", + ) scada_start: datetime | None = Field(None, description="SCADA数据起始时间") scada_end: datetime | None = Field(None, description="SCADA数据结束时间") sensor_nodes: list[str] | None = Field(None, description="传感器节点列表") diff --git a/app/infra/db/timescaledb/internal_queries.py b/app/infra/db/timescaledb/internal_queries.py index 16d45fe..5126cf8 100644 --- a/app/infra/db/timescaledb/internal_queries.py +++ b/app/infra/db/timescaledb/internal_queries.py @@ -169,6 +169,38 @@ class InternalQueries: else: raise + @staticmethod + def query_latest_scada_time( + device_ids: List[str], + before_time: str | datetime | None = None, + db_name: str = None, + max_retries: int = 3, + ) -> datetime | None: + """Return the latest SCADA timestamp for the selected devices.""" + before_dt = ( + parse_utc_time(before_time, field_name="before_time") + if before_time is not None + else None + ) + for attempt in range(max_retries): + try: + conn_string = ( + get_timescaledb_pgconn_string(db_name=db_name) + if db_name + else get_timescaledb_pgconn_string() + ) + with psycopg.Connection.connect(conn_string) as conn: + return ScadaRepository.get_latest_scada_time_sync( + conn, + device_ids, + before_dt, + ) + except Exception: + if attempt < max_retries - 1: + time.sleep(1) + else: + raise + @staticmethod def query_realtime_simulation_by_ids_timerange( element_ids: List[str], diff --git a/app/infra/db/timescaledb/repositories/scada.py b/app/infra/db/timescaledb/repositories/scada.py index d5a6348..b28dfea 100644 --- a/app/infra/db/timescaledb/repositories/scada.py +++ b/app/infra/db/timescaledb/repositories/scada.py @@ -54,6 +54,27 @@ class ScadaRepository: ) return cur.fetchall() + @staticmethod + def get_latest_scada_time_sync( + conn: Connection, + device_ids: List[str], + before_time: datetime | None = None, + ) -> datetime | None: + with conn.cursor(row_factory=dict_row) as cur: + if before_time is None: + cur.execute( + "SELECT max(time) AS time FROM scada.scada_data WHERE device_id = ANY(%s)", + (device_ids,), + ) + else: + cur.execute( + "SELECT max(time) AS time FROM scada.scada_data " + "WHERE device_id = ANY(%s) AND time <= %s", + (device_ids, before_time), + ) + row = cur.fetchone() + return row["time"] if row else None + @staticmethod async def get_scada_field_by_id_time_range( conn: AsyncConnection, diff --git a/app/services/burst_detection.py b/app/services/burst_detection.py index 9665934..e9d35a4 100644 --- a/app/services/burst_detection.py +++ b/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: diff --git a/tests/unit/test_burst_detection_service.py b/tests/unit/test_burst_detection_service.py new file mode 100644 index 0000000..d901d90 --- /dev/null +++ b/tests/unit/test_burst_detection_service.py @@ -0,0 +1,160 @@ +from __future__ import annotations + +from datetime import datetime, timedelta, timezone + +import numpy as np +import pandas as pd + +from app.services import burst_detection + + +TARGET = datetime(2026, 6, 20, 5, 30, tzinfo=timezone.utc) + + +def _complete_records(*, target: datetime, offset: float = 0.0) -> list[dict]: + start = target - timedelta(days=15) + timedelta(minutes=15) + return [ + { + "time": (start + timedelta(minutes=15 * index)).isoformat(), + "value": float(index % 96) + offset, + } + for index in range(15 * 96) + ] + + +def test_build_target_pressure_aligns_timestamps_and_excludes_incomplete_sensor( + monkeypatch, +): + nodes = [f"J{index}" for index in range(6)] + mapping = {node: f"D{index}" for index, node in enumerate(nodes)} + scada_data = { + query_id: _complete_records(target=TARGET, offset=float(index)) + for index, query_id in enumerate(mapping.values()) + } + scada_data["D5"] = scada_data["D5"][:-1] + + monkeypatch.setattr( + burst_detection, + "_get_pressure_sensor_mapping", + lambda _network: mapping, + ) + monkeypatch.setattr( + burst_detection.InternalQueries, + "query_latest_scada_time", + lambda **_kwargs: TARGET, + ) + monkeypatch.setattr( + burst_detection.InternalQueries, + "query_scada_by_ids_timerange", + lambda **_kwargs: scada_data, + ) + + frame, resolved_target, excluded = ( + burst_detection._build_target_pressure_from_scada( + network="test", + sensor_nodes=nodes, + requested_target_time=TARGET, + sampling_interval_minutes=15, + points_per_day=96, + ) + ) + + assert resolved_target == TARGET + assert frame.shape == (1440, 5) + assert frame.index[-1].to_pydatetime() == TARGET + assert excluded == [ + {"sensor_node": "J5", "reason": "missing_or_invalid_samples"} + ] + + +def test_target_mode_uses_fixed_parameters_and_only_classifies_target(monkeypatch): + index = pd.date_range( + start=TARGET - timedelta(days=15) + timedelta(minutes=15), + end=TARGET, + freq="15min", + ) + values = np.tile(np.arange(96, dtype=float), 15) + frame = pd.DataFrame( + {f"J{sensor}": values + sensor for sensor in range(5)}, + index=index, + ) + monkeypatch.setattr( + burst_detection, + "_get_pressure_sensor_nodes", + lambda _network: list(frame.columns), + ) + monkeypatch.setattr( + burst_detection, + "_build_target_pressure_from_scada", + lambda **_kwargs: (frame, TARGET, []), + ) + + payload = burst_detection.run_burst_detection( + network="test", + username="tester", + sampling_interval_minutes=15, + ) + + assert payload["target_time"] == TARGET.isoformat() + assert payload["sample_count"] == 1440 + assert payload["points_per_day"] == 96 + assert payload["algorithm_params"]["mu"] == 1 + assert payload["summary"]["score_threshold"] == -0.04 + assert [row["Role"] for row in payload["rows"]].count("target") == 1 + assert all(not row["IsBurst"] for row in payload["rows"][:-1]) + assert payload["reference_window"] == { + "start": (TARGET - timedelta(days=14)).isoformat(), + "end": (TARGET - timedelta(days=1)).isoformat(), + "day_count": 14, + } + + +def test_sampling_interval_uses_scada_frequency_and_can_be_overridden(monkeypatch): + monkeypatch.setattr( + burst_detection, + "get_all_scada_info", + lambda _network: [ + { + "type": "pressure", + "associated_element_id": "J1", + "transmission_frequency": "0:15:00", + }, + { + "type": "pressure", + "associated_element_id": "J2", + "transmission_frequency": "0:15:00", + }, + ], + ) + + assert ( + burst_detection._resolve_sampling_interval_minutes( + network="test", + sensor_nodes=["J1", "J2"], + requested_interval=None, + ) + == 15 + ) + assert ( + burst_detection._resolve_sampling_interval_minutes( + network="test", + sensor_nodes=["J1", "J2"], + requested_interval=30, + ) + == 30 + ) + + +def test_target_threshold_is_applied_only_to_latest_row(): + result = pd.DataFrame( + { + "Day": [1, 2, 3], + "Score": [-0.3, -0.2, -0.04], + "Prediction": [-1, -1, 1], + "IsBurst": [True, True, False], + } + ) + + rows = burst_detection._serialize_result_rows(result, target_only=True) + + assert [row["IsBurst"] for row in rows] == [False, False, True]