气象站异常检测:基于图信号处理与时间序列分析的Python实现
简介一套面向计算机、信号处理方向课程设计的气象站异常检测系统源码包基于Python实现通过图模型对气象站空间关系建模结合纬度差与时间序列历史差异识别异常并融合两类结果提升准确率。压缩包共5个文件涵盖Python主程序、MAT格式数据文件、项目说明文档以及两篇PDF数字信号处理课程大作业报告和信号处理图论前沿论文整体大小约2.92MB压缩包结构精简目录层级清晰便于快速定位与阅读。资料已有70人学习读者可从中获取完整可运行代码、实验数据、算法说明与报告范本既能直接复现异常检测流程也可为课设答辩或图信号处理研究提供参考并可在现有融合算法基础上做进一步改进与扩展。适合需要完成气象监测类项目或学习空间/时间序列融合算法的本专科生及初级开发者。1. 气象站异常检测系统这份 Python 源码到底在检测什么气象站的日平均气温数据表面上是一列数字实际上隐藏着两个维度的问题哪些站点在空间上“不合群”哪些站点在时间上“不对劲”。这份基于 Python 的气象站异常检测系统源码核心思路不是简单地拿阈值卡数据而是用图信号处理Graph Signal ProcessingGSP把气象站之间的空间关系建构成一张图再结合时间序列分析做双重判定。对于做课程设计、数字信号处理方向研究或者刚接触异常检测的开发者来说这套代码的价值在于它把“空间时间”两个维度的检测逻辑完整串了一遍而且配有 MATLAB 数据文件和可以直接运行的 main.py不是那种只有片段、跑不起来的半成品。我拿到这份资源后的第一感觉是它比常见的纯统计阈值方案多了一层空间建模的视角。气象站不是孤立存在的相邻站点的气温存在空间相关性如果一个站点的数据与它的空间邻居差异过大那大概率是设备故障或记录错误。下面我从源码实际实现的角度把这条检测链路拆开讲清楚。2. 空间分析用图信号处理给气象站“编网”算局部变异量2.1 为什么用图模型而不是简单阈值传统的气象数据异常检测最常见的是单站阈值法某个站点的气温超出历史极值范围就标记为异常。这个方法实现简单但有两个硬伤第一它完全忽略了站点之间的空间关系一个站点单独看可能没超阈值但相对于周围站点它就是明显偏移的第二气象数据本身有空间平滑性相邻站点的气温通常是接近的这个先验信息阈值法根本用不上。图模型的做法是把问题换个角度描述把每个气象站当成图的一个节点站点之间的连线边表示它们存在空间关联边的权重体现关联强度。这样一来“某个站点是否异常”就变成了“这个节点上的信号值气温与它的邻居节点是否一致”。一致性好说明该站点数据可信一致性差说明它很可能出了问题。这种建模方式把空间关系显式写进了检测逻辑里比事后拿经纬度做距离筛选要系统得多。2.2 邻接矩阵与拉普拉斯矩阵的构建这份源码的空间分析模块第一步是根据气象站的纬度信息构建邻接矩阵。核心逻辑是两个气象站的纬度差越小它们之间的空间距离越近边的权重就应该越大。这里可以直接用高斯核函数来定义权重这在图信号处理里是主流做法。import numpy as np from scipy.spatial.distance import pdist, squareform def build_adjacency_matrix(lats, sigma1.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, axis1) 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_size30, threshold_std3.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, alpha0.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_scorealpha是融合权重取 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 matplotlibmain.py运行完成后终端会打印检测到的异常气象站编号列表同时项目目录下会生成一张可视化图片把异常站点在图上标出来。我把源码完整读了一遍其实主流程的核心代码量不大真正的难点全在参数选择和边界情况处理上。比如空间分析里的sigma、时间序列里的window_size和threshold_std、融合阶段的alpha这几个参数每一个都直接影响结果。4.3 结果输出与参数调节源码默认输出的是一份异常站点列表但你要真拿去交作业或者落地光有编号不够最好能把异常信息结构化导出。我在复现时直接在main.py尾部加了一段导出逻辑把每个异常站点的编号、纬度、异常分数、异常类型空间异常/时间异常/融合异常写进 CSV 文件。import csv def export_results(station_ids, scores, output_pathanomaly_results.csv): 导出异常检测结果到 CSV 文件 with open(output_path, w, newline, encodingutf-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格式这种格式底层是 HDF5scipy的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 时每个站点几乎都只跟自己最近的一两个邻居相关局部变异量变得非常敏感一条正常的冷锋过境就会被标记成大面积异常。从那以后我每次调完参数都会强制自己跑一遍“极端值测试”确认参数在合理范围内的响应是平滑的。这份源码的价值不在于它能直接给你一个完美部署的异常检测服务而在于它把空间图建模、时间序列分析和融合判决这条链路完整地示范了一遍——把这条链路吃透换数据、换场景、换要素都是水到渠成的事。希望这份拆解能帮你少走一些我走过的弯路。本文还有配套的精品资源点击获取