3. 监测数据异常值识别与标记

任务与数据约定

本节针对某站点一段逐小时 PM2.5 浓度序列，完成三件事：用物理量程约束先剔除不可能值，再用滑动统计（滑动中位数 + MAD）识别局部离群点，最后输出带标记列的明细与异常记录清单，并比较不同阈值策略下的标记数量。示例数据为教学构造，用于演示判断逻辑，不代表真实监测结果。

约定：PM2.5 单位 μg/m³，物理有效量程取 [0, 1000]；滑动窗口取中心对齐的 24 点（约 1 天），最小有效点数为 12；稳健 z 分数定义为 z = 0.6745 × (x − 滑动中位数) / MAD，其中 MAD 为窗口内 |x − 中位数| 的中位数。阈值策略分别取 |z| > 3.5、|z| > 3.0、|z| > 2.5 三档进行对比。

实现代码

下面代码自包含：构造一段含缺失、量程越界与局部尖峰的示例序列，完成量程过滤、滑动稳健统计、标记与清单输出，并统计三档阈值下的标记数量。运行环境为 Python 3 + pandas + numpy。

import numpy as np
import pandas as pd

# 1) 构造教学示例序列：72 小时，含缺失、越界值与局部尖峰
rng = np.random.default_rng(7)
idx = pd.date_range("2024-03-01 00:00", periods=72, freq="h")
base = 45 + 12 * np.sin(np.arange(72) / 24 * 2 * np.pi)
pm25 = base + rng.normal(0, 4, 72)

pm25[10] = np.nan          # 缺失
pm25[20] = -3.0            # 负值，量程越界
pm25[35] = 1500.0          # 超上限，量程越界
pm25[50] = 180.0           # 局部尖峰，量程内但统计离群

df = pd.DataFrame({"timestamp": idx, "pm25": pm25})

# 2) 物理量程约束
LOW, HIGH = 0.0, 1000.0
df["range_bad"] = df["pm25"].notna() & ((df["pm25"] < LOW) | (df["pm25"] > HIGH))

# 3) 滑动中位数与 MAD（中心对齐，窗口 24，最少 12 个有效点）
WIN, MINP = 24, 12
roll = df["pm25"].rolling(window=WIN, center=True, min_periods=MINP)
med = roll.median()
mad = roll.apply(lambda s: np.median(np.abs(s - np.median(s))), raw=True)

# 4) 稳健 z 分数；MAD 为 0 时置为 NaN，避免除零
z = 0.6745 * (df["pm25"] - med) / mad.replace(0, np.nan)
df["roll_median"] = med
df["mad"] = mad
df["robust_z"] = z

# 5) 三档阈值标记与计数
for t in (3.5, 3.0, 2.5):
    df[f"outlier_z{t}"] = z.abs() > t

df["is_outlier"] = df["range_bad"] | df[f"outlier_z3.0"]
df["reason"] = np.select(
    [df["range_bad"], df["outlier_z3.0"]],
    ["range", "robust_z"],
    default=""
)

# 6) 输出标记明细与异常清单
cols = ["timestamp", "pm25", "roll_median", "mad", "robust_z", "is_outlier", "reason"]
marked = df[cols]
outliers = df.loc[df["is_outlier"], cols].sort_values("timestamp")

print("各阈值策略标记数量：")
for t in (3.5, 3.0, 2.5):
    print(f"  |z| > {t}: {int(df[f'outlier_z{t}'].sum())}")
print(f"  量程越界: {int(df['range_bad'].sum())}")
print(f"  合并标记(is_outlier): {int(df['is_outlier'].sum())}")
print("\n异常记录清单：")
print(outliers.to_string(index=False))

结果判读与阈值选择

示例序列中，量程越界点有 2 个（-3.0 与 1500.0），它们不参与滑动统计的有效点计数，但会被 range_bad 直接标记。局部尖峰 180.0 位于量程内，需依靠稳健 z 分数识别：由于窗口内中位数约在 45 附近、MAD 较小，该点的 |z| 通常远大于 3.5，因此三档阈值都会把它标出。阈值越宽松，被标记的点越多，其中会混入正常波动；阈值越严格，越可能漏掉幅度较小的真实异常。

• 先量程、后统计：量程越界值若参与滑动统计，会污染中位数与 MAD，因此应先剔除或置为缺失再计算。
• MAD 为 0 的处理：当窗口内超过一半的点取值相同（如长时间恒定低浓度），MAD 为 0，z 分数无定义，应置为 NaN 而非判为异常。
• 阈值选择依据：以 |z| > 3.0 作为起点，检查异常清单中是否包含已知的仪器故障、传输错误或真实污染过程；若清单中大量为正常日变化峰值，说明阈值偏松，应上调。
• 标记而非删除：保留原始值与 reason 列，便于后续插补或复核，避免不可逆的数据损失。