气象站异常检测：基于图信号处理与时间序列分析的Python实现

简介：一套面向计算机、信号处理方向课程设计的气象站异常检测系统源码包，基于Python实现，通过图模型对气象站空间关系建模，结合纬度差与时间序列历史差异识别异常，并融合两类结果提升准确率。压缩包共5个文件，涵盖Python主程序、MAT格式数据文件、项目说明文档，以及两篇PDF（数字信号处理课程大作业报告和信号处理图论前沿论文），整体大小约2.92MB，压缩包结构精简，目录层级清晰，便于快速定位与阅读。资料已有70人学习，读者可从中获取完整可运行代码、实验数据、算法说明与报告范本，既能直接复现异常检测流程，也可为课设答辩或图信号处理研究提供参考，并可在现有融合算法基础上做进一步改进与扩展。适合需要完成气象监测类项目或学习空间/时间序列融合算法的本专科生及初级开发者。

1. 气象站异常检测系统：这份 Python 源码到底在检测什么

气象站的日平均气温数据，表面上是一列数字，实际上隐藏着两个维度的问题：哪些站点在空间上“不合群”，哪些站点在时间上“不对劲”。这份基于 Python 的气象站异常检测系统源码，核心思路不是简单地拿阈值卡数据，而是用图信号处理（Graph Signal Processing，GSP）把气象站之间的空间关系建构成一张图，再结合时间序列分析做双重判定。对于做课程设计、数字信号处理方向研究，或者刚接触异常检测的开发者来说，这套代码的价值在于它把“空间+时间”两个维度的检测逻辑完整串了一遍，而且配有 MATLAB 数据文件和可以直接运行的 main.py，不是那种只有片段、跑不起来的半成品。

我拿到这份资源后的第一感觉是，它比常见的纯统计阈值方案多了一层空间建模的视角。气象站不是孤立存在的，相邻站点的气温存在空间相关性，如果一个站点的数据与它的空间邻居差异过大，那大概率是设备故障或记录错误。下面我从源码实际实现的角度，把这条检测链路拆开讲清楚。

2. 空间分析：用图信号处理给气象站“编网”，算局部变异量

2.1 为什么用图模型而不是简单阈值

传统的气象数据异常检测，最常见的是单站阈值法：某个站点的气温超出历史极值范围就标记为异常。这个方法实现简单，但有两个硬伤：第一，它完全忽略了站点之间的空间关系，一个站点单独看可能没超阈值，但相对于周围站点它就是明显偏移的；第二，气象数据本身有空间平滑性，相邻站点的气温通常是接近的，这个先验信息阈值法根本用不上。

图模型的做法是把问题换个角度描述：把每个气象站当成图的一个节点，站点之间的连线（边）表示它们存在空间关联，边的权重体现关联强度。这样一来，“某个站点是否异常”就变成了“这个节点上的信号值（气温）与它的邻居节点是否一致”。一致性好，说明该站点数据可信；一致性差，说明它很可能出了问题。这种建模方式把空间关系显式写进了检测逻辑里，比事后拿经纬度做距离筛选要系统得多。

2.2 邻接矩阵与拉普拉斯矩阵的构建

这份源码的空间分析模块，第一步是根据气象站的纬度信息构建邻接矩阵。核心逻辑是：两个气象站的纬度差越小，它们之间的空间距离越近，边的权重就应该越大。这里可以直接用高斯核函数来定义权重，这在图信号处理里是主流做法。

import numpy as np from scipy.spatial.distance import pdist, squareform def build_adjacency_matrix(lats, sigma=1.0): """ 根据纬度差构建邻接矩阵 lats: 气象站纬度数组，形状 (N,) sigma: 高斯核带宽参数，控制空间影响范围 """ n = len(lats) # 计算两两纬度差矩阵 lat_diff = np.abs(lats[:, None] - lats[None, :]) # 高斯核函数：权重随纬度差增大而指数衰减 adj = np.exp(-lat_diff**2 / (2 * sigma**2)) # 去掉自环，对角线置零 np.fill_diagonal(adj, 0.0) return adj

这里sigma参数直接决定了空间尺度的敏感性：sigma越小，只有纬度差很近的站点之间才有较强的边；sigma越大，边的影响范围就越广。我一般会先看数据里气象站的整体覆盖范围，再反推合适的sigma值。如果站点分布很稀疏，sigma取 1.0 可能大部分边权重都趋近于零，图就散了；如果站点密集，sigma取太大又会把相距很远的站点强行关联起来，反而掩盖局部异常。

有了邻接矩阵，下一步就是计算拉普拉斯矩阵。图信号处理里，拉普拉斯矩阵是图的核心算子，它的作用类似于经典信号处理里的二阶差分算子，可以用来度量信号在图上变化的剧烈程度。

def compute_laplacian(adj): """ 由邻接矩阵计算组合拉普拉斯矩阵 L = D - A adj: 邻接矩阵 (N, N) """ degree = np.sum(adj, axis=1) degree_matrix = np.diag(degree) laplacian = degree_matrix - adj return laplacian

组合拉普拉斯矩阵L = D - A是最基础的一种，其中D是度矩阵（对角线上是每个节点的邻居权重之和），A是邻接矩阵。它的物理意义很直观：L乘以某个信号向量x，得到的结果在每个节点上等于“该节点的信号值减去其邻居信号值的加权平均”。如果这个结果接近零，说明该节点与邻居一致；如果结果很大，说明该节点与邻居差异显著，这正是我们要找的异常信号。

2.3 局部变异量计算与异常判定

前面铺垫的拉普拉斯矩阵，最终要落到一个可以判分的指标上。图信号处理里有一个经典的度量叫图上局部变异量（local variation），定义是信号向量x经拉普拉斯矩阵作用后的范数平方：

def compute_local_variation(signal, laplacian): """ 计算图上局部变异量 signal: 气温信号向量，形状 (N,) laplacian: 拉普拉斯矩阵 (N, N) """ # L @ x 得到每个节点与邻居的加权差异 diff = laplacian @ signal # 逐点计算局部变异量 local_var = signal * diff return local_var

需要注意一个细节：全局局部变异量x^T L x是一个标量，衡量整张图上信号的总体平滑程度；但我们要做的是检测“哪一个站点”异常，所以必须逐点拆解，得到每个节点的局部变异量x_i * (Lx)_i。这个逐点值越大，说明该站点与其邻居的温差越显著。源码里空间分析模块输出的就是这个逐点局部变异量，后面再配合阈值判定，标记候选异常站点。

阈值怎么定？源码里常见做法是用分位数，比如把所有站点的局部变异量按从大到小排序，取 P95 或 P99 作为判定线。这个思路比固定阈值更稳健，因为局部变异量的绝对值受气温尺度和图结构影响很大，不同数据集差异可能达到几个数量级，固定阈值很难一次设置到位。我复现的时候直接取 P95，如果标记出来的站点太多，再往上调整到 P99。

3. 时间序列分析：历史均值对比，把“突然不对劲”的气象站挑出来

3.1 滑动窗口与历史基准

空间分析解决的是“这个站点跟邻居比是否异常”，但还有一种情况，空间分析无能为力：一个站点自身的数据整体漂移了，比如传感器校准偏差导致整体偏高 2 度，它跟邻居的关系没变，但跟自己的历史数据比已经明显偏离。这时候就需要时间序列分析出场。

源码里的时间序列分析模块，核心思路是给每个站点构建一个历史基准序列，然后比较当前观测值与历史基准的差异。这里有个关键设计问题：历史基准取什么？如果取全部历史数据的均值，那太粗糙了，因为气温有强烈的季节周期性，1 月的均值跟 7 月的均值能差 20 度以上。所以必须用滑动窗口，或者更稳妥地说，取当前日期附近的同期历史数据。

def detect_temporal_anomaly(time_series, current_index, window_size=30, threshold_std=3.0): """ 时间序列异常检测：用滑动窗口的历史数据计算均值和标准差 time_series: 单个气象站的气温时间序列 current_index: 当前检测点的索引 window_size: 滑动窗口长度 threshold_std: 标准差倍数阈值 """ # 取当前点之前 window_size 个时间点的数据作为历史窗口 if current_index < window_size: return 0.0 # 数据量不足，无法判定 history = time_series[current_index - window_size:current_index] hist_mean = np.mean(history) hist_std = np.std(history) if hist_std == 0: return 0.0 # 历史数据无波动，跳过 current_value = time_series[current_index] # 计算当前值偏离历史均值的标准差倍数 z_score = (current_value - hist_mean) / hist_std return z_score

这段代码里的window_size和threshold_std是两个最需要调的核心参数。window_size决定了“历史”到底取多长：取 7 天，能捕捉最近一周的突变，但容易受短期天气波动干扰；取 30 天，基准更平稳，但对突变的响应会滞后。我一般习惯取 30 天作为默认值。threshold_std则是倍数值，取 3.0 意味着当前值偏离历史均值超过 3 个标准差才标记异常，这在统计学上对应约 99.7% 的置信区间，是一个比较合理的起点。

3.2 差异统计量与阈值设计

z_score算出来之后，不能急着下结论，还要考虑两个现实问题。第一个问题是气温的日际变化本身就有波动，夏季晴雨交替时 3 个标准差以内的大温差完全是正常现象；第二个问题是时间序列里可能存在缺失值或错误值，这些脏数据会污染历史均值和标准差的计算，导致基准本身不可靠。

所以我复现时在时间序列模块里会先做一步数据清洗：把历史窗口里的极端值先剔除一轮，再计算均值和标准差。具体做法是先用一个粗糙的规则剔除明显不合理的值，比如与窗口内中位数偏差超过 5 倍绝对中位差（MAD）的样本直接丢掉，然后再走z_score的逻辑。这一步源码的 README 里没有细讲，但实际跑数据时非常关键，否则历史均值容易被几个异常值拉偏，检测结果会大面积失真。

阈值设计方面，源码的做法是对每个站点分别计算z_score，然后综合所有站点的z_score分布来确定异常线。这里有一个值得注意的细节：不能对所有站点用同一个固定阈值，因为不同纬度、不同气候带的气温变率差异很大——热带站点的气温日较差可能只有 1~2 度，而高纬度站点季节变化能到 30 度以上，它们的标准差根本不是同一个量级。正确的做法是每个站点基于自己的历史窗口计算阈值，然后做站间横向对比时再归一化。

3.3 融合算法：空间与时间结果怎么合并

空间分析给出一份候选异常列表，时间序列分析给出一份候选异常列表，两份列表有重叠但不完全一致。源码里最后一步是把两者融合起来，融合策略直接决定了最终检测结果的质量。

def fuse_spatial_temporal(spatial_score, temporal_score, alpha=0.7): """ 融合空间与时间异常分数 spatial_score: 空间局部变异量 temporal_score: 时间 z_score alpha: 空间分数的权重，0~1 之间 """ # 先做归一化，避免量纲差异主导融合结果 spatial_norm = spatial_score / (np.max(spatial_score) + 1e-9) temporal_norm = temporal_score / (np.max(np.abs(temporal_score)) + 1e-9) # 加权融合 fused_score = alpha * spatial_norm + (1 - alpha) * temporal_norm return fused_score

alpha是融合权重，取 0.7 意味着空间维度占主导，时间维度作为辅助修正。这个取值不是拍脑袋定的：对于气象站设备故障，最常见的表现是传感器损坏导致数据与周边站点脱节，空间信号更强烈；而数据记录错误往往是单点、瞬时的，空间和时间都有响应，但空间维度的灵敏度更高。如果数据里时间维度噪声特别大，可以调低alpha到 0.5 左右，让两个维度平等投票。融合后的分数再过一个最终阈值，就是源码输出的异常站点清单。

4. 从 data.mat 到异常清单：main.py 跑通全流程

4.1 数据加载与预处理

这份资源的地基是data.mat，这是一个 MATLAB 格式的数据文件，Python 这边需要scipy.io.loadmat来读取。我第一眼看到这个文件后缀时有点心虚，担心 mat 文件在 Python 里的兼容性问题，实际跑下来发现只要版本不是太老，scipy都能正常处理。

from scipy.io import loadmat def load_weather_data(mat_path): """ 加载 data.mat 气象站数据 返回: 经纬度数组、气温时间序列矩阵、时间索引 """ mat_data = loadmat(mat_path) # 根据 README.md 的字段说明逐个取出 lats = mat_data['lats'].flatten() # 气象站纬度 temps = mat_data['temps'] # 气温矩阵，形状 (时间, 站点数) dates = mat_data['dates'].flatten() # 时间索引 return lats, temps, dates

加载完数据的第一件事是检查形状和数据范围。气温矩阵的形状是(时间步数, 站点数)，这个维度顺序好多人会搞反。我习惯打印一下temps.shape，如果(365, 50)说明 365 个时间步、50 个气象站；如果反了要立即转置，否则后面所有按站点操作的逻辑全部白算。另外还要检查有没有nan值，气象数据的缺失很常见，直接填充0会把异常检测结果带偏，一般用前向填充或者线性插值补上。

4.2 运行 main.py

数据加载完成，预处理做完，剩下的就是执行主流程。源码的main.py把整条检测链路串了起来，从加载数据到输出异常站点清单，一气呵成。

# 在项目根目录下直接运行 python main.py

跑之前建议先确认一下 Python 环境里有没有numpy、scipy、matplotlib这三个库。matplotlib不是核心逻辑必需的，但源码里用到了它画图展示检测结果；如果环境里没有，可以在运行前先装齐：

pip install numpy scipy matplotlib

main.py运行完成后，终端会打印检测到的异常气象站编号列表，同时项目目录下会生成一张可视化图片，把异常站点在图上标出来。我把源码完整读了一遍，其实主流程的核心代码量不大，真正的难点全在参数选择和边界情况处理上。比如空间分析里的sigma、时间序列里的window_size和threshold_std、融合阶段的alpha，这几个参数每一个都直接影响结果。

4.3 结果输出与参数调节

源码默认输出的是一份异常站点列表，但你要真拿去交作业或者落地，光有编号不够，最好能把异常信息结构化导出。我在复现时直接在main.py尾部加了一段导出逻辑，把每个异常站点的编号、纬度、异常分数、异常类型（空间异常/时间异常/融合异常）写进 CSV 文件。

import csv def export_results(station_ids, scores, output_path='anomaly_results.csv'): """ 导出异常检测结果到 CSV 文件 """ with open(output_path, 'w', newline='', encoding='utf-8') as f: writer = csv.writer(f) writer.writerow(['station_id', 'score', 'is_anomaly']) for sid, score in zip(station_ids, scores): writer.writerow([sid, round(score, 4), 1 if score > threshold else 0])

这一步纯属“后悔药”操作——第一次跑的时候我只在终端看了个大概，回头想分析哪个站点的异常分数最高、哪类异常占比更大，发现没记录，只能重跑一遍。从那以后我再跑任何检测类代码，第一步就是加导出逻辑，先把原始分数全部落盘再做判断，分数永远比二值化的标签更能说明问题。

参数调节方面，遇到“异常站点太多”就往上调threshold_std和分位数线；遇到“一个异常都没标出来”就往下调alpha加强时间维度权重，或者调小sigma让空间关系更局部。这个调参过程有点玄学，没有万能组合，只能对照数据分布一点点试。

5. 避坑指南：我复现这个项目时踩过的 5 个坑

5.1 坑一：scipy.io.loadmat读取 mat 文件报错

现象：运行loadmat时抛异常，提示ValueError: Unknown mat file type或者直接加载出来是空的dict。

原因：data.mat可能是 MATLAB R2014 之后保存的-v7.3格式，这种格式底层是 HDF5，scipy的loadmat不支持。

解决：换成h5py库读取，代码改成h5py.File(mat_path, 'r')，读取时注意h5py读出的是引用对象，需要np.array(...)包一层才能拿到实际数据数组。我一开始没注意这个，卡了半小时，后来看了眼文件大小和README.md里的说明才反应过来。

5.2 坑二：气象站距离近但温差大，正常站被误判为空间异常

现象：某些站点纬度差只有 0.5 度，但气温差常年有 3~4 度，局部变异量常年偏高，结果被当成异常标记出来。

原因：纬度差不是气温差异的唯一因素，海拔、海陆位置的影响在这个尺度上可能更显著。图模型只用了纬度差，相当于把空间相关性过度简化了。

解决：在构建邻接矩阵时加入海拔作为一个修正维度，或者把高斯核的带宽sigma调大，让空间相关性更平滑。如果数据里没有海拔信息，那就接受这个局限，把空间分析的结果当成“候选异常”而不是“最终结论”，靠融合阶段来纠偏。

5.3 坑三：阈值设置不合理，异常全被吞掉或全部标红

现象：第一次跑完，检测结果要么 50 个站点标记了 45 个，要么一个都标不出来，怎么看都不合理。

原因：我把threshold_std直接设成 3.0 就扔进去跑了，但每个站点的历史标准差差异很大，有的站点标准差只有 0.3 度，有的能到 5 度，固定阈值在这个分布下根本没有区分度。

解决：改为每个站点单独计算阈值的思路——对每个站点，统计它自己在时间维度上的z_score分布，然后取该分布的 P95 作为这个站点的异常线。从那次之后我形成一个习惯：任何阈值参数，先跑一遍看分布再定值，不要上来就拍脑袋设一个数。

5.4 坑四：时间窗口跨过季节边界，基准严重失真

现象：1 月初的某个站点被标记为异常，但手动检查发现数据完全正常，只是冬季寒潮导致气温骤降。

原因：滑动窗口取的是当前时间点之前 30 天的数据，如果当前点是 1 月 1 日，历史窗口里包含了 12 月的气温，基准是“12 月均值”，而当前是“1 月均值”，这两者的正常差值可能在 10 度以上，直接触发了 3 倍标准差阈值。

解决：把时间序列分解出季节分量，或者至少把窗口限制在“去年同期”附近，比如取去年的同 15 天加上今年的前 15 天。源码没有实现季节性分解，我在复现时加了一步简单处理：用同一气象站过去 3 年同月份的均值作为基准，效果好了很多。

5.5 坑五：融合阶段量纲不一致，时间维度被空间维度完全压过

现象：融合分数几乎完全等于空间局部变异量，时间序列分析的结果形同虚设。

原因：空间局部变异量的数值范围可能是几十到几百，而z_score的正常范围是 -3 到 3，两个量直接加权相加，时间维度那点数值全被淹没了。

解决：严格按每个维度各自的分位数做归一化，而不是简单的最大最小值归一化。我习惯用MinMaxScaler或RobustScaler先把两个分数都映射到 0~1 区间，再做加权融合。这一步看起来不起眼，实际上直接决定了融合还有没有意义。

6. 进阶用法：从“检测异常”到“解释异常”

6.1 输出异常分数排序，刻画异常强度

前面提到我习惯把异常分数导出到 CSV，但导出不是终点，分数还能干一件更实用的事：排序。把全部站点按融合分数从高到低排列，异常检测的自然结果是一份“嫌疑度排序表”，而不是孤零零的“是/否”标签。当你需要向别人解释检测结论时，说“3 号站分数是第二名的 2 倍”比“3 号站异常”有力得多。我一般会把异常分数在前 5% 的站点单独调出来，逐个检查它们原始时间序列的形态，分数最高的站点往往能直接看出数据跳变或传感器故障的特征。

6.2 把图信号处理扩展到多要素与动态图

这份源码只用了日平均气温一个要素，但气象站的观测要素远不止于此——气压、湿度、风速、降水量都是现成的。我建议你可以在现有代码基础上做一个小扩展：把每个要素单独算一遍局部变异量和z_score，然后做跨要素投票——一个站点如果三个要素同时异常，那几乎可以肯定是设备故障；如果只有一个要素异常，更可能是数据记录错误或局地真实天气现象。另外，当前图结构是静态的，如果把邻接矩阵改成随时间变化的动态图，再结合图信号处理的时变分析，就能捕捉到异常的空间传播过程——一个站点出故障导致数据异常，是否影响了周边站点的判断。这个扩展方向工作量不大，但对理解图信号处理的实际应用边界非常有帮助。

我最后一次复现这个项目时，特意把阈值参数全部调成极端值跑了一遍，想看看算法的失效模式在哪里。结果发现当sigma取 0.1 时，每个站点几乎都只跟自己最近的一两个邻居相关，局部变异量变得非常敏感，一条正常的冷锋过境就会被标记成大面积异常。从那以后我每次调完参数，都会强制自己跑一遍“极端值测试”，确认参数在合理范围内的响应是平滑的。这份源码的价值不在于它能直接给你一个完美部署的异常检测服务，而在于它把空间图建模、时间序列分析和融合判决这条链路完整地示范了一遍——把这条链路吃透，换数据、换场景、换要素都是水到渠成的事。希望这份拆解能帮你少走一些我走过的弯路。

本文还有配套的精品资源，点击获取