1. 环境时间序列数据获取与预处理

季节性分解与趋势评估的质量,几乎在数据进入分解算法之前就已经被决定了。无论使用 STL、X-11 还是回归加季节虚拟变量,模型只能忠实地反映输入序列的结构:如果时间戳错位、缺失值被随意填补、异常峰值未被识别,分解出的"季节项"很可能只是采样制度或设备故障的投影。本节围绕环境监测数据的来源、清洗、缺失值处理、异常值检测与时间序列对齐五个环节展开,目标是产出一条时间连续、量纲一致、异常可解释的分析用序列。

一、数据来源与结构特征

环境时间序列通常来自三类源:空气质量自动监测站(SO₂、NO₂、PM₁₀、PM₂.₅、O₃、CO,常见为小时均值)、水质自动站或手工采样(COD、氨氮、总磷、溶解氧、pH、浊度,小时或日均)、气象观测(气温、相对湿度、风速风向、气压、降水,小时或日值)。三者的时间粒度、缺测机制和量纲差异明显:气象数据通常完整度高、连续性最好;空气质量数据受设备校准和停电影响会出现成段缺失;水质数据在手工采样情形下往往是周或月尺度,且采样时刻不固定。

因此第一步不是清洗,而是建立数据字典:记录每个字段的单位、时间粒度、统计方式(小时均值/日均值/瞬时值)、监测点位编号与经纬度、数据来源与获取时间。这一步看似行政性工作,实际决定了后续能否正确对齐——例如"日均值"是按自然日 24 个小时值平均,还是按有效小时数不少于 18 个计算,不同口径会直接改变序列的方差结构。

二、时间序列对齐:先统一时间轴

多源数据合并前必须统一到同一时间轴。推荐做法是:以分析目标粒度为基准(如日均),将所有序列重采样到该粒度,并显式声明聚合规则。小时转日用时,对浓度类指标取算术平均,对降水取日累计,对风向不可直接平均(应转为 u/v 分量后再合成)。对齐时需注意三个易错点:时区统一(全部转为本地标准时间或 UTC 并注明)、夏令时导致的重复/缺失小时、以及跨年边界上"第 53 周"的处理。

下面的图示展示了一条典型的多源数据从接入到形成分析序列的处理链,以及每个环节的关键判断点:

多源环境监测数据预处理链 空气质量自动站 SO₂/NO₂/PM₂.₅/O₃ · 小时 水质监测 COD/氨氮/总磷 · 日或周 气象观测 气温/湿度/风速/降水 · 小时 统一时间轴与重采样 时区统一 · 小时→日均 · 聚合规则显式声明 数据清洗与异常值检测 物理范围校验 · 3σ/IQR · 设备故障标记 缺失值处理 短缺口线性插值 · 长缺口标记不填补 分析用序列(含质量标记) 进入季节性分解与趋势评估 判断:异常还是真实峰值? 对照气象过程与运维记录 判断:缺口长度阈值 ≤3 点插值;>3 点保留缺失

三、数据清洗与异常值检测

清洗分两层。第一层是物理与逻辑校验:浓度不得为负,PM₂.₅ 不应大于 PM₁₀,相对湿度应在 0–100% 之间,风速为 0 时风向无意义。这类规则可以直接用条件表达式批量标记。第二层是统计异常检测,常用两种互补方法:

  • 3σ 准则:适用于近似正态的序列,将偏离均值超过 3 倍标准差的点标为可疑。缺点是标准差本身受极端值影响,且环境数据常有明显右偏(如 PM₁₀),此时可先做对数变换。
  • IQR 箱线准则:以 Q1−1.5×IQR 与 Q3+1.5×IQR 为界,对偏态分布更稳健,适合作为首选筛查。

关键原则是:检测出的异常点不等于错误点。一次沙尘过程会让 PM₁₀ 出现远超 3σ 的峰值,这是真实的物理信号,若直接删除会削弱季节项与趋势项的估计。正确做法是结合气象记录、站点运维日志与上下游站点比对,将异常分为"设备故障型"(应置为缺失)与"真实极端型"(应保留并可在稳健分解中降权)。建议在数据表中增加 qc_flag 字段,取值如 0=正常、1=可疑、2=缺失、3=极端真实值,使后续每一步都可追溯。

四、缺失值处理策略

缺失机制决定处理方式。若缺失完全随机(如偶发通信中断),短缺口可用线性插值或前向填充;若缺失与季节或工况相关(如冬季某设备停机),简单插值会引入偏差。实践中可按缺口长度分级:连续缺失不超过 3 个时间点,用线性插值;3–24 个点,可用同季节同时段的多年均值或 STL 初步拟合后的残差补全;超过 24 个点(日均序列约一个月),建议保留缺失并在分解时使用能处理缺失的算法(如 STL 的 na.action 设置),而不是强行填补。任何填补都应在数据表中留下标记,避免"填补值"与"实测值"混为一谈。

对季节性分解而言,缺失造成的最大危害不是信息损失,而是周期性断裂:若某年某月整段缺失并被均值填补,该月的季节指数会被系统性拉平,进而污染整条季节曲线。因此长缺口宁可保留,也不要填。

五、一个可复现的处理示例

假设我们拿到某站点 2021–2023 年的 PM₂.₅ 小时数据与同期气象数据,目标是构造日均序列用于 STL 分解。可按以下顺序执行(示例参数为合理假设,非实测结果):

1. 读取原始表,统一时间列为本地标准时间,去除重复时间戳(保留最新修订版本)
2. 物理校验:PM2.5 < 0 或 > 1000 μg/m³ 标记 qc_flag=1
3. 小时→日均:有效小时数 ≥ 18 才计算日均,否则记为缺失
4. 异常检测:对 log(PM2.5) 计算 IQR 边界,越界点标记 qc_flag=1
5. 人工复核:对照沙尘/重污染过程与运维记录,将设备故障点改为 qc_flag=2
6. 缺失处理:连续缺口 ≤3 天线性插值;>3 天保留缺失
7. 输出:date, pm25, qc_flag, temp, rh, wind_speed 七列宽表

完成后的序列应做一次完整性体检:统计各年有效天数、最长连续缺失段、异常标记比例。若某年有效天数低于 300 天,或最长缺失段超过 15 天,应在分析报告中明确说明,因为这会影响年度趋势项的可信度。经过这样处理的序列,才适合进入下一阶段的季节性分解——此时分解结果中的季节项反映的是真实的环境周期,而不是采样缺陷的残留。