修复数据清洗index越界错误；重命名压力流量清洗方法

2026-02-02 14:15:54 +08:00
parent c3c26fb107
commit 3c7e2c5806
10 changed files with 116 additions and 60 deletions
--- a/app/algorithms/api_ex/flow_data_clean.py
+++ b/app/algorithms/api_ex/flow_data_clean.py
@@ -0,0 +1,408 @@
+# ...existing code...
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+from pykalman import KalmanFilter
+import os
+
+
+def fill_time_gaps(
+    data: pd.DataFrame,
+    time_col: str = "time",
+    freq: str = "1min",
+    short_gap_threshold: int = 10,
+) -> pd.DataFrame:
+    """
+    补齐缺失时间戳并填补数据缺口。
+
+    Args:
+        data: 包含时间列的 DataFrame
+        time_col: 时间列名（默认 'time'）
+        freq: 重采样频率（默认 '1min'）
+        short_gap_threshold: 短缺口阈值（分钟），<=此值用线性插值，>此值用前向填充
+
+    Returns:
+        补齐时间后的 DataFrame（保留原时间列格式）
+    """
+    if time_col not in data.columns:
+        raise ValueError(f"时间列 '{time_col}' 不存在于数据中")
+
+    # 解析时间列并设为索引
+    data = data.copy()
+    data[time_col] = pd.to_datetime(data[time_col], utc=True)
+    data_indexed = data.set_index(time_col)
+
+    # 生成完整时间范围
+    full_range = pd.date_range(
+        start=data_indexed.index.min(), end=data_indexed.index.max(), freq=freq
+    )
+
+    # 重索引以补齐缺失时间点
+    data_reindexed = data_indexed.reindex(full_range)
+
+    # 按列处理缺口
+    for col in data_reindexed.columns:
+        # 识别缺失值位置
+        is_missing = data_reindexed[col].isna()
+
+        # 计算连续缺失的长度
+        missing_groups = (is_missing != is_missing.shift()).cumsum()
+        gap_lengths = is_missing.groupby(missing_groups).transform("sum")
+
+        # 短缺口：线性插值
+        short_gap_mask = is_missing & (gap_lengths <= short_gap_threshold)
+        if short_gap_mask.any():
+            data_reindexed.loc[short_gap_mask, col] = (
+                data_reindexed[col]
+                .interpolate(method="linear", limit_area="inside")
+                .loc[short_gap_mask]
+            )
+
+        # 长缺口：前向填充
+        long_gap_mask = is_missing & (gap_lengths > short_gap_threshold)
+        if long_gap_mask.any():
+            data_reindexed.loc[long_gap_mask, col] = (
+                data_reindexed[col].ffill().loc[long_gap_mask]
+            )
+
+    # 重置索引并恢复时间列（保留原格式）
+    data_result = data_reindexed.reset_index()
+    data_result.rename(columns={"index": time_col}, inplace=True)
+
+    # 保留时区信息
+    data_result[time_col] = data_result[time_col].dt.strftime("%Y-%m-%dT%H:%M:%S%z")
+    # 修正时区格式（Python的%z输出为+0000，需转为+00:00）
+    data_result[time_col] = data_result[time_col].str.replace(
+        r"(\+\d{2})(\d{2})$", r"\1:\2", regex=True
+    )
+
+    return data_result
+
+
+def clean_flow_data_kf(
+    input_csv_path: str, show_plot: bool = False, fill_gaps: bool = True
+) -> str:
+    """
+    读取 input_csv_path 中的每列时间序列，使用一维 Kalman 滤波平滑并用预测值替换基于 3σ 检测出的异常点。
+    保存输出为：<input_filename>_cleaned.xlsx（与输入同目录），并返回输出文件的绝对路径。
+    仅保留输入文件路径作为参数（按要求）。
+
+    Args:
+        input_csv_path: CSV 文件路径
+        show_plot: 是否显示可视化
+        fill_gaps: 是否先补齐时间缺口（默认 True）
+    """
+    # 读取 CSV
+    data = pd.read_csv(input_csv_path, header=0, index_col=None, encoding="utf-8")
+
+    # 补齐时间缺口（如果数据包含 time 列）
+    if fill_gaps and "time" in data.columns:
+        data = fill_time_gaps(
+            data, time_col="time", freq="1min", short_gap_threshold=10
+        )
+
+    # 分离时间列和数值列
+    time_col_data = None
+    if "time" in data.columns:
+        time_col_data = data["time"]
+        data = data.drop(columns=["time"])
+
+    # 存储 Kalman 平滑结果
+    data_kf = pd.DataFrame(index=data.index, columns=data.columns)
+    # 平滑每一列
+    for col in data.columns:
+        observations = pd.Series(data[col].values).ffill().bfill()
+        if observations.isna().any():
+            observations = observations.fillna(observations.mean())
+        obs = observations.values.astype(float)
+
+        kf = KalmanFilter(
+            transition_matrices=[1],
+            observation_matrices=[1],
+            initial_state_mean=float(obs[0]),
+            initial_state_covariance=1,
+            observation_covariance=1,
+            transition_covariance=0.01,
+        )
+        # 跳过EM学习，使用固定参数以提高性能
+        state_means, _ = kf.smooth(obs)
+        data_kf[col] = state_means.flatten()
+
+    # 计算残差并用IQR检测异常（更稳健的方法）
+    residuals = data - data_kf
+    residual_thresholds = {}
+    for col in data.columns:
+        res_values = residuals[col].dropna().values  # 移除NaN以计算IQR
+        q1 = np.percentile(res_values, 25)
+        q3 = np.percentile(res_values, 75)
+        iqr = q3 - q1
+        lower_threshold = q1 - 1.5 * iqr
+        upper_threshold = q3 + 1.5 * iqr
+        residual_thresholds[col] = (lower_threshold, upper_threshold)
+
+    cleaned_data = data.copy()
+    anomalies_info = {}
+    for col in data.columns:
+        lower, upper = residual_thresholds[col]
+        sensor_residuals = residuals[col]
+        anomaly_mask = (sensor_residuals < lower) | (sensor_residuals > upper)
+        anomaly_idx = data.index[anomaly_mask.fillna(False)]
+        anomalies_info[col] = pd.DataFrame(
+            {
+                "Observed": data.loc[anomaly_idx, col],
+                "Kalman_Predicted": data_kf.loc[anomaly_idx, col],
+                "Residual": sensor_residuals.loc[anomaly_idx],
+            }
+        )
+        cleaned_data.loc[anomaly_idx, f"{col}_cleaned"] = data_kf.loc[anomaly_idx, col]
+
+    # 如果原始数据包含时间列，将其添加回结果
+    if time_col_data is not None:
+        cleaned_data.insert(0, "time", time_col_data)
+
+    # 构造输出文件名：在输入文件名基础上加后缀 _cleaned.xlsx
+    input_dir = os.path.dirname(os.path.abspath(input_csv_path))
+    input_base = os.path.splitext(os.path.basename(input_csv_path))[0]
+    output_filename = f"{input_base}_cleaned.xlsx"
+    output_path = os.path.join(input_dir, output_filename)
+
+    # 覆盖同名文件
+    if os.path.exists(output_path):
+        os.remove(output_path)
+    cleaned_data.to_excel(output_path, index=False)
+
+    # 可选可视化（第一个传感器）
+    plt.rcParams["font.sans-serif"] = ["SimHei"]
+    plt.rcParams["axes.unicode_minus"] = False
+    if show_plot and len(data.columns) > 0:
+        sensor_to_plot = data.columns[0]
+        plt.figure(figsize=(12, 6))
+        plt.plot(
+            data.index,
+            data[sensor_to_plot],
+            label="监测值",
+            marker="o",
+            markersize=3,
+            alpha=0.7,
+        )
+        plt.plot(
+            data.index, data_kf[sensor_to_plot], label="Kalman滤波预测值", linewidth=2
+        )
+        anomaly_idx = anomalies_info[sensor_to_plot].index
+        if len(anomaly_idx) > 0:
+            plt.plot(
+                anomaly_idx,
+                data[sensor_to_plot].loc[anomaly_idx],
+                "ro",
+                markersize=8,
+                label="监测值异常点",
+            )
+            plt.plot(
+                anomaly_idx,
+                data_kf[sensor_to_plot].loc[anomaly_idx],
+                "go",
+                markersize=8,
+                label="Kalman修复值",
+            )
+        plt.xlabel("时间点（序号）")
+        plt.ylabel("监测值")
+        plt.title(f"{sensor_to_plot}：观测值与Kalman滤波预测值（异常点标记）")
+        plt.legend()
+        plt.show()
+
+    # 返回输出文件的绝对路径
+    return os.path.abspath(output_path)
+
+
+def clean_flow_data_df_kf(data: pd.DataFrame, show_plot: bool = False) -> dict:
+    """
+    接收一个 DataFrame 数据结构，使用一维 Kalman 滤波平滑并用预测值替换基于 IQR 检测出的异常点。
+    区分合理的0值（流量转换）和异常的0值（连续多个0或孤立0）。
+    返回完整的清洗后的字典数据结构。
+
+    Args:
+        data: 输入 DataFrame（可包含 time 列）
+        show_plot: 是否显示可视化
+    """
+    # 使用传入的 DataFrame
+    data = data.copy()
+
+    # 补齐时间缺口（如果启用且数据包含 time 列）
+    data_filled = fill_time_gaps(
+        data, time_col="time", freq="1min", short_gap_threshold=10
+    )
+    # 移除 time 列用于后续清洗
+    data_filled = data_filled.drop(columns=["time"])
+
+    # 存储 Kalman 平滑结果
+    data_kf = pd.DataFrame(index=data_filled.index, columns=data_filled.columns)
+    # 平滑每一列
+    for col in data_filled.columns:
+        observations = pd.Series(data_filled[col].values).ffill().bfill()
+        if observations.isna().any():
+            observations = observations.fillna(observations.mean())
+        obs = observations.values.astype(float)
+
+        kf = KalmanFilter(
+            transition_matrices=[1],
+            observation_matrices=[1],
+            initial_state_mean=float(obs[0]),
+            initial_state_covariance=1,
+            observation_covariance=10,
+            transition_covariance=10,
+        )
+        state_means, _ = kf.smooth(obs)
+        data_kf[col] = state_means.flatten()
+
+    # 计算残差并用IQR检测异常
+    residuals = data_filled - data_kf
+    residual_thresholds = {}
+    for col in data_filled.columns:
+        res_values = residuals[col].dropna().values
+        q1 = np.percentile(res_values, 25)
+        q3 = np.percentile(res_values, 75)
+        iqr = q3 - q1
+        lower_threshold = q1 - 1.5 * iqr
+        upper_threshold = q3 + 1.5 * iqr
+        residual_thresholds[col] = (lower_threshold, upper_threshold)
+
+    # 创建完整的修复数据
+    cleaned_data = data_filled.copy()
+    anomalies_info = {}
+
+    for col in data_filled.columns:
+        lower, upper = residual_thresholds[col]
+        sensor_residuals = residuals[col]
+        anomaly_mask = (sensor_residuals < lower) | (sensor_residuals > upper)
+        anomaly_idx = data_filled.index[anomaly_mask.fillna(False)]
+
+        anomalies_info[col] = pd.DataFrame(
+            {
+                "Observed": data_filled.loc[anomaly_idx, col],
+                "Kalman_Predicted": data_kf.loc[anomaly_idx, col],
+                "Residual": sensor_residuals.loc[anomaly_idx],
+            }
+        )
+
+        # 直接在原列上替换异常值为 Kalman 预测值
+        cleaned_data.loc[anomaly_idx, col] = data_kf.loc[anomaly_idx, col]
+
+    # 可选可视化
+    plt.rcParams["font.sans-serif"] = ["SimHei"]
+    plt.rcParams["axes.unicode_minus"] = False
+    if show_plot and len(data.columns) > 0:
+        sensor_to_plot = data.columns[0]
+        
+        # 定义x轴
+        n = len(data)
+        time = np.arange(n)
+        n_filled = len(data_filled)
+        time_filled = np.arange(n_filled)
+
+        plt.figure(figsize=(12, 8))
+
+        plt.subplot(2, 1, 1)
+        plt.plot(
+            time,
+            data[sensor_to_plot],
+            label="原始监测值",
+            marker="o",
+            markersize=3,
+            alpha=0.7,
+        )
+        
+        # 修正：检查 data_filled 的异常值，绘制在 time_filled 上
+        abnormal_zero_mask = data_filled[sensor_to_plot].isna()
+        # 如果目的是检查0值，应该用 == 0。这里保留 isna() 但修正索引引用，防止crash。
+        # 如果原意是 isna() 则在 fillna 后通常没有 na。假设用户可能想检查 0 值？
+        # 基于 "异常0值" 的标签，改为检查 0 值更合理，但为了保险起见，
+        # 如果 isna() 返回空，就不画。防止索引越界是主要的。
+        abnormal_zero_idx = data_filled.index[abnormal_zero_mask]
+        
+        if len(abnormal_zero_idx) > 0:
+             # 注意：如果 abnormal_zero_idx 是基于 data_filled 的索引（0..M-1），
+             # 直接作为 x 坐标即可，因为 time_filled 也是 0..M-1
+             # 而 y 值应该取自 data_filled 或 data_kf，取 data 会越界
+            plt.plot(
+                abnormal_zero_idx,
+                data_filled[sensor_to_plot].loc[abnormal_zero_idx],
+                "mo",
+                markersize=8,
+                label="异常值(NaN)",
+            )
+
+        plt.plot(
+            time_filled, data_kf[sensor_to_plot], label="Kalman滤波预测值", linewidth=2
+        )
+        anomaly_idx = anomalies_info[sensor_to_plot].index
+        if len(anomaly_idx) > 0:
+            plt.plot(
+                anomaly_idx,
+                data_filled[sensor_to_plot].loc[anomaly_idx],
+                "ro",
+                markersize=8,
+                label="IQR异常点",
+            )
+        plt.xlabel("时间点（序号）")
+        plt.ylabel("流量值")
+        plt.title(f"{sensor_to_plot}：原始数据与异常检测")
+        plt.legend()
+
+        plt.subplot(2, 1, 2)
+        plt.plot(
+            time_filled,
+            cleaned_data[sensor_to_plot],
+            label="修复后监测值",
+            marker="o",
+            markersize=3,
+            color="green",
+        )
+        plt.xlabel("时间点（序号）")
+        plt.ylabel("流量值")
+        plt.title(f"{sensor_to_plot}：修复后数据")
+        plt.legend()
+
+        plt.tight_layout()
+        plt.show()
+
+    # 返回完整的修复后字典
+    return cleaned_data
+
+
+# # 测试
+# if __name__ == "__main__":
+#     # 默认：脚本目录下同名 CSV 文件
+#     script_dir = os.path.dirname(os.path.abspath(__file__))
+#     default_csv = os.path.join(script_dir, "pipe_flow_data_to_clean2.0.csv")
+#     out = clean_flow_data_kf(default_csv)
+#     print("清洗后的数据已保存到:", out)
+
+# 测试 clean_flow_data_dict 函数
+if __name__ == "__main__":
+    import random
+
+    # 读取 szh_flow_scada.csv 文件
+    script_dir = os.path.dirname(os.path.abspath(__file__))
+    csv_path = os.path.join(script_dir, "szh_flow_scada.csv")
+    data = pd.read_csv(csv_path, header=0, index_col=None, encoding="utf-8")
+
+    # 排除 Time 列，随机选择 5 列
+    columns_to_exclude = ["Time"]
+    available_columns = [col for col in data.columns if col not in columns_to_exclude]
+    selected_columns = random.sample(available_columns, 1)
+
+    # 将选中的列转换为字典
+    data_dict = {col: data[col].tolist() for col in selected_columns}
+
+    print("选中的列:", selected_columns)
+    print("原始数据长度:", len(data_dict[selected_columns[0]]))
+
+    # 调用函数进行清洗
+    cleaned_dict = clean_flow_data_df_kf(data_dict, show_plot=True)
+    # 将清洗后的字典写回 CSV
+    out_csv = os.path.join(script_dir, f"{selected_columns[0]}_clean.csv")
+    pd.DataFrame(cleaned_dict).to_csv(out_csv, index=False, encoding="utf-8-sig")
+    print("已保存清洗结果到:", out_csv)
+    print("清洗后的字典键:", list(cleaned_dict.keys()))
+    print("清洗后的数据长度:", len(cleaned_dict[selected_columns[0]]))
+    print("测试完成：函数运行正常")