Files
TJWaterServerBinary/app/services/burst_detection.py
T
jiang db6032bd84
Server CI/CD / docker-image (push) Successful in 24s
Server CI/CD / deploy-fallback-log (push) Has been cancelled
feat(burst-detection): update scada analysis flow
2026-07-17 16:31:46 +08:00

765 lines
26 KiB
Python

from __future__ import annotations
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
from app.infra.db.timescaledb.internal_queries import InternalQueries
from app.services.scheme_management import (
query_burst_detection_scheme_detail,
query_burst_detection_schemes,
scheme_name_exists,
store_scheme_info,
)
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,
username: str,
observed_pressure_data: (
pd.DataFrame
| dict[str, list[Any]]
| list[dict[str, Any]]
| list[list[Any]]
| None
) = None,
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,
scheme_name: str | None = None,
data_source: str = "monitoring",
simulation_scheme_name: str | None = None,
simulation_scheme_type: str | None = None,
) -> dict[str, Any]:
"""
运行爆管侦测服务入口。
调用方式三选一:
- 不传数据时间窗,自动侦测最近完整时刻;可用 `target_time` 回放历史时刻
- 直接传 `observed_pressure_data`
- 或传 `scada_start/scada_end` 让后端自动查询 SCADA 压力数据
`observed_pressure_data` 支持格式:
- `pd.DataFrame`
行表示时间点,列表示传感器;列名应为传感器/节点 ID。
- `dict[str, list[Any]]`
键为传感器/节点 ID,值为按时间顺序排列的压力序列。
例如:`{"J1": [101.2, 101.0], "J2": [99.8, 99.7]}`。
- `list[dict[str, Any]]`
每个元素代表一个时间点的多传感器观测。
例如:`[{"J1": 101.2, "J2": 99.8}, {"J1": 101.0, "J2": 99.7}]`。
- `list[list[Any]]`
二维数组式 JSON,格式为 `(时间点数, 传感器数)`。
这是最接近原始 `burst_detector` 示例代码的调用方式。
数据约束:
- 统一要求“行=时间点,列=传感器”。
- 总样本点数必须能被 `points_per_day` 整除。
- 至少要有 2 天数据,即 `sample_count >= 2 * points_per_day`。
- 若传入 `sensor_nodes`,输入数据必须包含这些列;SCADA 模式下也会只按这些节点取数。
"""
if not network:
raise ValueError("network is required.")
selected_sensor_nodes = (
list(dict.fromkeys([node for node in (sensor_nodes or []) if node]))
if sensor_nodes
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
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
else _get_pressure_sensor_nodes(network)
)
if data_source == "simulation":
if not simulation_scheme_name:
raise ValueError("模拟方案模式必须提供 simulation_scheme_name。")
observed_df = _build_observed_pressure_from_simulation(
network=network,
sensor_nodes=scada_sensor_nodes,
scada_start=scada_start,
scada_end=scada_end,
simulation_scheme_name=simulation_scheme_name,
simulation_scheme_type=simulation_scheme_type,
)
observed_input = observed_df
observed_source = "simulation_scheme_timerange"
else:
observed_df = _build_observed_pressure_from_scada(
network=network,
sensor_nodes=scada_sensor_nodes,
scada_start=scada_start,
scada_end=scada_end,
)
observed_input = observed_df
observed_source = "backend_timerange"
else:
if observed_pressure_data is None:
raise ValueError(
"未提供 observed_pressure_data,且未提供 scada_start/scada_end。"
)
observed_input = observed_pressure_data
observed_source = "request_payload"
detector = BurstDetector(
mu=mu,
points_per_day=points_per_day,
iforest_params=iforest_params,
)
result_df = detector.run_detection(
observed_input,
sensor_nodes=selected_sensor_nodes,
)
resolved_sensor_nodes = list(result_df.attrs.get("sensor_nodes", []))
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,
"observed_source": observed_source,
"sample_count": int(result_df.attrs.get("sample_count", 0)),
"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": 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"
payload["simulation_scheme"] = {
"name": simulation_scheme_name,
"type": simulation_scheme_type,
}
else:
payload["data_source"] = "monitoring"
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(),
}
if scheme_name:
_store_burst_detection_scheme(
network=network,
scheme_name=scheme_name,
username=username,
payload=payload,
mu=mu,
points_per_day=points_per_day,
iforest_params=detector.iforest_params,
)
payload["scheme_name"] = scheme_name
return payload
def _build_observed_pressure_from_simulation(
*,
network: str,
sensor_nodes: list[str],
scada_start: datetime | str | None,
scada_end: datetime | str | None,
simulation_scheme_name: str | None,
simulation_scheme_type: str | None = None,
) -> pd.DataFrame:
if scada_start is None or scada_end is None:
raise ValueError("使用模拟方案查询时必须同时提供 scada_start 与 scada_end。")
start_dt = _to_datetime(scada_start)
end_dt = _to_datetime(scada_end)
if start_dt >= end_dt:
raise ValueError("SCADA 时间窗非法:scada_start 必须早于 scada_end。")
# Reuse burst_location logic partially here or call internal queries directly
# For burst detection, we need time series for all sensor nodes.
# Check for missing nodes in simulation result if needed, but InternalQueries handles some of it.
# We assume sensor_nodes are valid pressure nodes.
scheme_type = simulation_scheme_type or "burst_analysis"
simulation_data = InternalQueries.query_scheme_simulation_by_ids_timerange(
db_name=network,
scheme_type=scheme_type,
scheme_name=simulation_scheme_name,
element_ids=sensor_nodes,
start_time=start_dt.isoformat(),
end_time=end_dt.isoformat(),
element_type="node",
field="pressure",
)
# simulation_data is {sensor_id: [{time, value}, ...]}
# Convert to DataFrame: index=time, columns=sensor_ids
data_dict = {}
timestamps = set()
for sensor_id in sensor_nodes:
records = simulation_data.get(sensor_id, [])
if not records:
continue
# Convert records to Series with time index
ts_values = []
ts_index = []
for r in records:
if r.get("value") is not None:
ts_values.append(float(r["value"]))
ts_index.append(pd.to_datetime(r["time"]))
if ts_values:
s = pd.Series(ts_values, index=ts_index)
data_dict[sensor_id] = s
timestamps.update(ts_index)
if not data_dict:
raise ValueError("指定时间窗内未查询到模拟压力数据。")
observation_df = pd.DataFrame(data_dict)
# Handle missing timestamps if any (though simulation usually has uniform steps)
# Forward fill or interpolate might be needed if steps differ, but typically they align.
observation_df = observation_df.sort_index()
# Fill NaN if any missing points for some sensors
observation_df = observation_df.fillna(method="ffill").fillna(method="bfill")
if observation_df.empty:
raise ValueError("模拟压力数据无法构建观测矩阵。")
return observation_df
def list_burst_detection_schemes(
network: str,
query_date: datetime | str | None = None,
) -> list[dict[str, Any]]:
parsed_date = extract_date(query_date, field_name="query_date") if query_date is not None else None
return query_burst_detection_schemes(
name=network,
network=network,
query_date=parsed_date,
)
def get_burst_detection_scheme_detail(network: str, scheme_name: str) -> dict[str, Any]:
result = query_burst_detection_scheme_detail(network, scheme_name)
if not result:
raise ValueError(f"未找到爆管侦测方案: {scheme_name}")
return result
def _store_burst_detection_scheme(
*,
network: str,
scheme_name: str,
username: str,
payload: dict[str, Any],
mu: int,
points_per_day: int,
iforest_params: dict[str, Any],
) -> None:
if scheme_name_exists(network, scheme_name):
raise ValueError(f"方案名称已存在: {scheme_name}")
now_iso = utc_now().isoformat()
scheme_detail = {
"network": network,
"sensor_nodes": payload.get("sensor_nodes", []),
"observed_source": payload.get("observed_source"),
"scada_window": payload.get("scada_window"),
"algorithm_params": {
"mu": mu,
"points_per_day": points_per_day,
"iforest_params": iforest_params,
},
"result_summary": payload.get("summary", {}),
"result_payload": payload,
}
store_scheme_info(
name=network,
scheme_name=scheme_name,
scheme_type="burst_detection",
username=username,
scheme_start_time=now_iso,
scheme_detail=scheme_detail,
)
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]] = []
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": 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,
*,
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("爆管侦测结果为空。")
score_series = result_df["Score"]
most_anomalous_index = int(score_series.idxmin())
latest_row = rows[-1]
anomaly_days = [row["Day"] for row in rows if row["IsBurst"]]
summary = {
"burst_detected": bool(latest_row["IsBurst"]),
"latest_day": latest_row,
"most_anomalous_day": int(result_df.iloc[most_anomalous_index]["Day"]),
"anomaly_days": anomaly_days,
"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]]:
feature_matrix = result_df.attrs.get("high_freq_features")
sensor_nodes = list(result_df.attrs.get("sensor_nodes", []))
if feature_matrix is None or len(sensor_nodes) == 0:
return []
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,
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, 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):
if str(item.get("type", "")).lower() != "pressure":
continue
node_id = item.get("associated_element_id")
if isinstance(node_id, str) and node_id:
sensor_nodes.append(node_id)
sensor_nodes = list(dict.fromkeys(sensor_nodes))
if not sensor_nodes:
raise ValueError("未找到压力传感器对应节点(scada_info.type=pressure)。")
return sensor_nodes
def _build_observed_pressure_from_scada(
*,
network: str,
sensor_nodes: list[str],
scada_start: datetime | str | None,
scada_end: datetime | str | None,
) -> pd.DataFrame:
if scada_start is None or scada_end is None:
raise ValueError("使用后端 SCADA 查询时必须同时提供 scada_start 与 scada_end。")
start_dt = _to_datetime(scada_start)
end_dt = _to_datetime(scada_end)
if start_dt >= end_dt:
raise ValueError("SCADA 时间窗非法:scada_start 必须早于 scada_end。")
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:
preview = ", ".join(missing_nodes[:10])
raise ValueError(f"未找到可用于压力观测的 SCADA api_query_id: {preview}")
query_ids = [node_query_id[node_id] for node_id in sensor_nodes]
scada_data = InternalQueries.query_scada_by_ids_timerange(
db_name=network,
device_ids=query_ids,
start_time=start_dt.isoformat(),
end_time=end_dt.isoformat(),
)
available_lengths = [
len(scada_data.get(query_id, []))
for query_id in query_ids
if len(scada_data.get(query_id, [])) > 0
]
if not available_lengths:
raise ValueError("指定时间窗内未查询到压力 SCADA 数据。")
min_len = min(available_lengths)
observation_df = pd.DataFrame()
for node_id in sensor_nodes:
query_id = node_query_id[node_id]
records = scada_data.get(query_id, [])[:min_len]
if len(records) < min_len:
continue
observation_df[node_id] = [float(item["value"]) for item in records]
if observation_df.empty:
raise ValueError("SCADA 压力数据无法构建观测矩阵。")
return observation_df
def _to_datetime(value: datetime | str) -> datetime:
return parse_utc_time(value)