国家级地面气象站日值数据集V3.0清洗与特殊值处理实战

1. 这份V3.0数据集里到底装着什么

第一次拿到这批数据的人，十个里有八个会先愣一下——压缩包里没有一张图，没有一句解释，全是按年份摊开的纯文本文件。它不像很多公开数据集那样给你一个现成的CSV加一份字段说明，而是保留了气象业务系统里最原始的样子：一个要素、一年、一个文件。国家级地面气象站基本气象要素日值数据集这个名字里的每个词都有实际含义：国家级站意味着站点是经过长期业务运行、观测环境和仪器配置相对规范的骨干站，而不是后来大量补充的区域自动站；日值意味着每条记录对应一个站一天；基本气象要素则限定了它涵盖的是气温、气压、湿度、风、降水、蒸发、日照、地温这一套常规观测项目。

我把这套数据用在气候统计、农业气象指标计算和能源评估的交叉项目里已经有些年头了，越用越觉得它像一块未经打磨的原材料：信息密度极高，但直接下锅会硌牙。这篇文章想讲清楚的就是从"下载下来一堆TXT"到"能直接拿去做分析的面板数据"，中间到底隔着哪些不写进说明文档、但一定会绊你一脚的东西。

1.1 站点覆盖与时间跨度的真实形态

站点数量上，两三千个国家级站听着不多，但乘以七十余年的时间轴，再乘以逐日的记录频率，行数就到了千万级甚至上亿级别。真正让处理变麻烦的不是行数，而是站点的时间跨度极度不均匀。建站早的站能回溯到上世纪五十年代初，建站晚的站可能到七八十年代才有第一条记录，还有一些站中途撤销、迁址、恢复观测。所以你在做全国平均的时候，如果只是简单地对"当天有记录的站"求平均，得到的曲线里混进了站点构成变化带来的虚假趋势——1955年可能只有几百个站参与平均，2000年有两千多个站参与平均，两者不可比。

使用这套数据时我通常建议先做一张站点可用性矩阵：行是站号，列是年份，格子里填该年有效观测日数。这张矩阵一画出来，哪些站是"全程在场"的、哪些是"半路加入"的、哪些是"中间断了一截又续上"的，一目了然。做长时间序列分析时，用这张矩阵筛出一个"平衡面板"子集，虽然牺牲了站点数量，但换来的是序列的可比性。

时间维度上还有一层容易被忽略：年份文件的起止并不统一。有些年份的文件里包含了上一年末的补录记录，有些年份的12月31日之后还有额外的"年统计"记录混在末尾。我第一次做拼接时没做过滤，结果发现某几年凭空多出几十条"日期为空"的行，一度以为解析代码写错了，后来才确认是文件尾部附加的汇总行。

1.2 要素清单与实际可用性

把要素列出来容易，但每个要素在空间和时间上的可用性差别巨大，这一点比清单本身重要得多。

蒸发和日照这两项是典型的"看着有、实际用不了几个站"的要素。全国有蒸发观测的站长期维持在一个相对小的规模，而日照虽然站多一些，但序列的仪器断点问题比气温严重得多。我见过有人拿全站平均蒸发量画了一条1960—2020年的曲线，前段平滑上升后段急剧下降，其实那段下降几乎全是蒸发皿从20厘米口径换成E601大型蒸发皿造成的，跟气候没太大关系。

1.3 文件组织：一个要素一年一个TXT

目录结构基本是按"数据集标识—要素—年份"三级展开的，文件名里带了要素代号和年份。要素代号是一套内部编码，比如降水、气温、气压各有自己的编号，不同版本的说明文档在代号上可能有细微差别。我的习惯是不硬编码代号到脚本里，而是先扫一遍目录、把实际存在的文件名列成表，用文件名里的中文或英文要素关键词去匹配，再把匹配结果打印出来人工核对一遍。这看起来笨，但比事后发现"气温文件被当成了露点文件"要省事得多。

每个TXT文件的内部结构是固定的：开头若干行是台站信息头（站号、经纬度、海拔、年），后面是逐日记录，字段按固定顺序排列，用逗号分隔。这里有个经常被忽略的细节——同一个要素在不同年份的文件里，字段数量可能不一致。某些年份多了质量控制相关的辅助列，某些年份又少一两列。用固定列名的读法在跨年份合并时就会报错或者错位。

提示：不要用一个固定的列名列表去读所有年份。先写一个探测脚本，把每个文件的列数、首行内容打印出来汇总成表格，看清楚跨年份的字段结构差异，再决定统一映射方案。这一步花二十分钟，能省掉后面两天的排查。

2. 打开文件前的三件事：编码、分隔符、字段顺序

很多人栽在第一步不是技术难度，而是想当然。默认UTF-8读、默认逗号分隔、默认第一行是表头——这三条假设在这套数据上至少错两条。

2.1 编码踩坑与首行表头处理

文件是早期业务系统导出的，编码常见的是GB系（ANSI），直接用UTF-8读会抛UnicodeDecodeError，或者更阴险地读成功但站点名、备注字段变成乱码。乱码如果只出现在非关键字段还无所谓，但如果把乱码带进了分组键，后面做站点聚合时就会出现"同一个站分裂成两个站"的诡异现象。

处理方式很简单：读的时候显式指定encoding='gbk'，遇到个别站报错就用encoding_errors='replace'兜底，同时在读取后统计一次替换字符出现的次数，超过阈值就单独拎出来看。我一般会把每个文件的解码成功率记进日志，跑完一遍批处理再回头检查哪些年份有问题。

首行处理更微妙。有些文件第一行确实就是表头，有些文件第一行是站点信息头（不是列名），真正的列名在第二行或者根本没有列名。最稳的做法是先读前5行做人工判断：

import pathlib sample = next(pathlib.Path("raw").rglob("*.TXT")) with open(sample, encoding="gbk", errors="replace") as f: for i in range(5): print(i, repr(f.readline()[:200]))

看到内容长什么样，再决定header=参数怎么设。这一步不能省，我见过因为header=0把一个站点信息行当成了列名，结果所有列名变成"站号 经纬度 海拔..."，后面所有取列操作都在报 KeyError。

注意：字段分隔符在绝大多数文件里是英文逗号，但要确认没有混用中文逗号或制表符的情况。用repr()打印一行就能看出来，别靠肉眼看文本编辑器——肉眼分辨中英文逗号是会出错的。

2.2 用pandas读第一份文件的完整代码

确定结构之后，第一份文件读通，剩下的就是批量化。下面这段是我常用的模板，重点是显式指定每列的数据类型，避免pandas把站号推断成整数、把含特殊值的要素列推断成浮点后又悄悄变成NaN。

import pandas as pd COLUMNS = [ "station_id", "lat", "lon", "elev", "year", "month", "day", "precip_raw", "precip_qc", ] dtypes = { "station_id": "string", "lat": "float32", "lon": "float32", "elev": "float32", "year": "int16", "month": "int8", "day": "int8", "precip_raw": "int32", # 原始整数编码，先不转物理量 "precip_qc": "string", } df = pd.read_csv( sample, encoding="gbk", sep=",", header=None, # 先不认表头，读进来再按位置切 names=COLUMNS, dtype=dtypes, na_values=[], keep_default_na=False, # 关键：不让pandas自作主张把特殊值当缺失 )

keep_default_na=False这一句是我用血换来的。pandas默认会把空字符串、NA、NULL之类识别成缺失值，而气象数据里的缺失是用具体的整数编码表示的，例如某几个固定数值代表缺测、微量、累计不可用。如果不关掉默认NA识别，这些编码会被部分吃掉，部分保留，最后你在统计"有多少个缺测"时根本对不上账。

2.3 站点元数据表别漏掉

原始数据里虽然每行都带了经纬度、海拔，但那是"记录时的站点信息"还是"制表时的站点信息"，在不同年份并不一致。真正可靠的站点属性来自配套的台站元数据表，它包含建站时间、撤站时间、地址变更、观测场迁移、海拔变动等历史。这张表的价值在于给每一条日记录提供正确的时间上下文——一个站在1998年迁过址，那么1998年前后的记录来自不同的地表环境，直接连成一条序列做趋势分析是错误的。

元数据表通常不大，几百KB到几MB，但处理起来要小心：同一个站号可能出现多行（对应不同时段），需要按时间区间做"区间连接"（interval join），而不是简单的一对一merge。用pandas的话可以先merge_asof或按条件筛选，用DuckDB的ASOF JOIN会更舒服一些。

3. 特殊值才是这套数据真正的门槛

如果说前面两章是入门，这一章是分水岭。能不能把特殊值处理干净，直接决定了你后面所有的统计结果能不能用。

3.1 缺测、微量、静风：32XXX那套编码

这套数据里最反直觉的设计是：某些要素列用大整数编码来承载"非数值信息"。典型的就是三万多这一档的数字，它们不是"很大的降水量"，而是标记位。常见的语义大致如下（具体以随数据发布的说明为准，不同要素的约定略有差异）：

我第一次看到"32700"以为是个异常大值，还做过一次三倍标准差剔除，结果把"今天下了点毛毛雨"全删了。微量降水在干旱区的年降水统计里影响不大，但在"降水日数"这个指标上影响是决定性的——某北方站一年可能有二十多天微量降水，全删掉之后降水日数直接少了两成。

处理微量值的关键是先想清楚你要算什么。算年降水总量，微量记0没问题；算降水日数、算连续无降水日，微量必须记为一个有效的小值，否则会凭空造出并不存在的"长期干旱"。我现在的做法是新增一列标记is_trace，原始编码保留，物理量列里微量记0，需要用到"是否降水"的判断时去读is_trace列而不是读数值。

风向的静风编码是另一个例子。风向采用16方位编码，其中特定取值表示静风。静风不是"缺失"，它是一个有意义的观测结果。很多做风玫瑰图的脚本直接用dropna()把静风扔了，画出来的图会系统性地高估主导风向的频率。正确做法是在做角度统计前把静风单独计为一类，或者用"静风频率"单独出图。

3.2 质量控制码怎么用才不冤枉好数据

伴随大部分要素列的还有一个质量控制码字符。常见取值体系大致是：0表示正确，1表示可疑，2表示错误，8表示缺测，9表示未做质控。这个字段极其有用，但用错了比不用更糟。

最常见的错误是"一刀切只保留0"。这么做看起来很干净，实际会带来两个后果：一是样本量在早期年份骤减，因为早期质控规则和现在不同，很多现在判为"可疑"的值在当年是正常记录；二是引入选择性偏差，某些气象条件下（比如雷暴、强降温）观测值更容易被判为可疑，你把这些删掉相当于把极端事件系统性剔除了，最后算出来的极值偏小、变率偏小。

我的做法分场景：

做气候平均值、趋势分析时，保留0和1，剔除2和8，把9当作"无质控信息"按数值本身的合理性判断——如果落在气候学界限内就保留。

做极端值统计、重现期分析时，对1类值额外做人工核查：看相邻站同一天的记录是否也有类似量级，看前后几天的走势是否连续，能解释的保留，解释不了的标注保留但不纳入拟合。

做数据发布的成品时，才考虑最严格的过滤，因为发布出去的东西别人会直接引用，容错空间小。

提示：质量控制码在文件里通常和要素值配对出现，列的排列是"值、码、值、码"。如果你的解析把码列整体错位了一列，表面上看起来一切正常——因为码本身也是0~9的小整数——但所有质控过滤都作用在了错误的要素上。检查办法是统计码列的取值分布，正常应该高度集中在0附近，如果出现大量8或9混在意料之外的要素上，就该警惕错位了。

3.3 单位换算和"看起来正常其实错"的典型

量纲这件事，说明文档里写着，但实际文件里存在少数年份单位不一致的历史遗留。最常见的是气温和气压：绝大多数年份按0.1为单位存储，需要除以10得到℃和hPa，但个别早期年份可能直接用整数摄氏度。判断办法不是查文档，是看数值分布——一个站一年的平均气温编码如果集中在-300到400之间，那显然是0.1℃；如果集中在-30到40之间，那就是直接摄氏度。

还有一类错误不来自单位，来自统计口径。比如降水量，日值的截断时刻通常是20时，而有些分析者默认它是00时截断，去做"日降水与日平均气温的联合分析"时时间窗就错开了8小时。这8小时在强对流天气分析里足够让结论完全反过来。

地温的深层数据也要特别注意：0cm地表温度、5cm、10cm、15cm、20cm，一直到320cm。深层地温由于热传导的滞后，同一天的记录反映的其实是更早时段的热状况，做土壤热通量或者冻土分析时，这个相位差必须显式建模，不能当作同步观测用。

4. 日值的统计口径：时间序列断裂的根源

前面讲的都是单条记录的清洗，这一章讲的是记录与记录之间的一致性问题。这是做长期序列分析时最隐蔽、破坏力最大的一类坑，因为清洗完的数据看起来每一行都合法，但拼起来就是不连续。

4.1 日平均气温到底平均了几次

日平均值不是"天上掉下来的一个数"，它是从定时观测值算出来的。而这个"几次观测"在历史上变过。

早期阶段，很多站采用3次定时观测（08、14、20时）计算日平均，后来改为4次定时观测（02、08、14、20时）的算术平均。进入自动观测时代之后，部分站改为24次逐时值的平均。三种口径算出来的日平均气温，在日较差大的地区、季节差异明显的地方，能差出零点几度——听起来很小，但放到几十年趋势上，零点几度的口径跳变足以把你的升温速率估计抬高或压低一个可观的比例。

我第一次做区域平均气温年际变化时，曲线在某个年代出现了一个明显的台阶，查了很久以为是数据版本问题，最后定位到就是观测时次改变导致的统计口径跳变。

应对办法有两个方向。一是在元数据里找到口径变更的时间点，对变更前后的序列做分段处理，分别计算趋势再比较，而不是强行拼接。二是用统一的算法重算——如果原始数据同时提供了定时观测值，可以用同一套公式重新计算日平均，把口径统一起来。可惜V3.0的日值数据集本身只给日平均值，不给定时值，所以这条路走不通的话就只能走第一条。

4.2 蒸发皿换代与日照计换代

蒸发量的断点问题比气温严重得多，因为它不是"差一点点"，而是换了测量原理。小型蒸发皿和E601大型蒸发皿的水面面积、安装高度、防鸟设施都不同，测出来的蒸发量存在系统性偏差。业界有一些经验换算系数，但系数随气候区、季节、风速变化很大，直接用一个大系数硬乘，误差能和你要研究的信号一个量级。

我的处理原则是：蒸发量序列不做跨型号拼接。如果要研究长期趋势，只在同一型号的时段内做分析，把型号变更点当作序列的断点；如果必须给出全长序列，那就在图里把断点标出来，明确说明前后不可比。

日照时数的问题类似。传统的暗筒式日照计依靠纸带灼痕记录，换成光电式传感器之后，起测阈值和响应特性都变了。光电式对弱日照更敏感，通常记录到的日照时数会偏多一些。做"日照时数年变化"这类分析时，如果不做均一化订正，会把仪器换代误读成"云量减少、晴天变多"。

4.3 迁站、海拔与城市化

站点搬家是另一个必须处理的问题。观测场从一个位置迁到另一个位置，海拔可能变几十米，周边环境可能从农田变成城区，或者反过来。海拔变化直接影响气压，也影响气温（按干绝热递减率粗略估算，每升高100米气温降约0.6℃）；周边环境变化影响的是观测的代表性，一个原本代表郊区气候的站，随着城市扩张会逐渐变成城市热岛站。

处理迁站问题的标准工具是均一化检验。常用的方法包括基于相邻站参考序列的检验（比如标准正态均一性检验）、多断点检验方法等。这类方法的核心思路都一样：用周边站构造一条"本应如此"的参考序列，然后看目标站和参考序列的差值在哪些年份发生了突变，突变点就是潜在的断点，需要订正。

这件事的工作量不小，但如果你要做的是几十年尺度的趋势分析，它不能省。我见过不少分析直接把全国站的数据粗暴平均，得出"某区域升温速率远高于全球平均"的结论，其中相当一部分来自站点构成变化和未订正的断点，而不是真实的气候信号。

注意：均一化不是万能的。如果断点附近没有足够密的参考站，或者断点前后记录太短，检验的可靠性会大幅下降。这种情况下更诚实的做法是保留断裂标记，在结论里说明不确定性，而不是硬做订正然后当作精确值使用。

5. 一条可复现的清洗流水线

讲了这么多坑，得给一条能直接跑的路。下面这条流水线我用了几个项目，结构稳定，改改参数就能适配。

5.1 目录规范与批量转parquet

原始TXT文件数量多、体积大、每次重读都慢。第一步应该做的事情是一次性转换成分区存储的列式格式，之后所有分析都读parquet，不再碰原始文件。

目录设计上我按"要素/年份"两级分区：

raw/ 原始TXT，只读，永不修改 staging/ 中间产物 curated/ element=temp_avg/ year=1980/part.parquet year=1981/part.parquet element=precip/ year=1980/part.parquet

转换脚本的骨架：

import pandas as pd import pathlib def parse_one(path: pathlib.Path) -> pd.DataFrame: df = pd.read_csv( path, encoding="gbk", sep=",", header=None, names=COLUMNS, dtype=dtypes, keep_default_na=False, ) # 过滤掉非数据行（站点头、年统计行） df = df[df["month"].between(1, 12) & df["day"].between(1, 31)] df["date"] = pd.to_datetime( dict(year=df["year"], month=df["month"], day=df["day"]), errors="coerce", ) return df.dropna(subset=["date"]) for path in pathlib.Path("raw").rglob("*.TXT"): df = parse_one(path) if df.empty: continue year = int(df["year"].mode().iloc[0]) out = pathlib.Path(f"curated/element={element}/year={year}") out.mkdir(parents=True, exist_ok=True) df.to_parquet(out / f"{path.stem}.parquet", index=False)

几个细节值得说。month.between(1,12)这个过滤能干净地干掉文件尾部的年统计行和一些脏行。pd.to_datetime加errors="coerce"是为了兜住2月30日这种明显错误，转换不了的置为NaT然后丢掉。分区写入时保留原文件名，方便日后追溯是哪一份原始文件出的问题。

5.2 缺测率筛查与异常值三级过滤

数据进了parquet之后，第一件事不是分析，是体检。我通常做三张表：

异常值过滤我分三级，逐级加严，每一级都记录被剔除的数量，方便回查：

第一级是物理界限检查，用不可能出现的值做硬性拦截，比如相对湿度超出0–100、气压低于某个下限。这一级的规则要宽，只拦明显不可能的值，宁可漏掉也不误杀。

第二级是气候学界限检查，用同站同月的历史分位数做参照，超出某一分位范围较远的标记为可疑。这一级要分月做，因为气温的年循环太强，用全年分位数会误杀冬季正常低温和夏季正常高温。

第三级是空间一致性检查，看同一区域相邻站在同一天的记录是不是孤立异常。这一级最能抓住仪器故障造成的错误，但对站点稀疏的地区无能为力。

三级过滤的执行顺序不能反。如果先做空间检查，一个区域内如果有多个站同时出现相同故障（比如同一批设备），空间检查会认为"大家都有这个值，说明是对的"。先做单站物理和气候学检查，先清掉单站层面的问题，再做空间检查，效率更高也更准。

5.3 面板化与按需聚合

清洗完之后，把数据组织成(station_id, date)双索引的面板（panel）结构，是最方便后续操作的形态。用pandas的话可以用 MultiIndex，用xarray的话可以组织成(station, time)的二维数组。

panel = ( cleaned .set_index(["station_id", "date"]) [["temp_avg", "precip", "rh"]] .sort_index() )

组织成面板之后，很多操作就变成了一行：

# 各站点年降水总量（要求该年有效日数 >= 330） yearly_precip = ( panel["precip"] .groupby([panel.index.get_level_values("station_id"), panel.index.get_level_values("date").year]) .agg(total="sum", valid_days="count") .query("valid_days >= 330") )

这里valid_days的阈值设成330是我的习惯，留出一个月左右的容错空间。阈值设太高会把大量站点筛掉，设太低会把有明显缺测的年份带进来，把年总量算小。这个参数跟研究区域关系很大，湿润区降水日数多，缺几天影响小；干旱区一旦漏掉一场大雨，年总量能差一小半，阈值反而应该更严。

提示：sum()在pandas里默认跳过NaN，这在小样本上问题不大，但在缺测较多时会把年总量系统性低估。务必同时统计有效日数，用有效日数做质量门槛，而不是盲目相信sum()的结果。

6. 落到具体场景：几个真实用法和对应陷阱

数据清洗到位之后，剩下的就是应用。不同场景对数据的敏感点完全不同，这一章按场景说说各自的要点。

6.1 气候趋势：为什么你的升温速率可能偏高

做区域气温趋势，最常见的流水线是：筛出连续记录超过40年的站 → 计算每站年均温 → 做线性回归得速率 → 空间插值出图。这条流水线里至少有四个地方能让你得出偏高的速率。

第一个是站点筛选偏差。能连续观测40年以上的站，很多位于城市或城市周边，因为那里观测条件维持得好。而真正代表区域的乡村站、山区站，可能中途撤了或者记录不连续。筛完剩下的是一个偏向城市的子集，城市热岛效应被算进了"气候变化"。

第二个是起止年份的选取。如果起点选在一个冷年、终点选在一个暖年，速率会被显著放大。规范的做应做的是对整段序列做回归，而不是拿首尾两点的差去除以年数。

第三个是缺测填补方式。简单的线性插值会抹平极值，让方差变小但均值不受影响；如果缺测集中出现在某种天气类型下（比如大雾天观测困难），插值就会引入偏差。

第四个就是前面讲的均一化。未订正的断点可能系统性地抬高趋势。

我个人的经验是：任何超过30年的趋势结论，都要换两三种站点子集、两三种起止年份重新算一遍，看结果稳不稳定。如果不同设定下速率差了一倍以上，那结论本身就不该发表。

6.2 农业积温与物候：日均值口径的连锁影响

农业上用这套数据算积温、算无霜期、算越冬期低温，这些指标对日平均值口径的敏感度比气候趋势高得多。

以 ≥10℃ 活动积温为例，它是把全年所有日平均气温超过10℃的日子加起来。口径变化造成的零点几度差异，会在临界日期附近决定某一天到底"算进去"还是"不算进去"，而这个决定会被累加到全年的积温里。一个生长季如果边界上有三五天在临界值附近摆动，积温能差出几十度日，足以影响品种布局的推荐结论。

无霜期也有类似问题。日最低气温接近0℃的日子，口径差异可能让"有霜"变"无霜"，直接改变无霜期长度。

我的做法是在做农业指标之前，先用同期有定时观测值的年限做一次交叉验证：用定时值按不同公式算日平均，再和数据集给的日平均值对比，量化出偏差的量级，然后在结论里给出误差区间，而不是给一个精确到小数点的积温数字。农业上的应用方通常更关心"这个品种能不能种"，而不是"积温是3250还是3260"，给出区间反而更有用。

6.3 光伏与建筑负荷：辐射缺失时怎么用气温和日照补

能源评估里，最想要的是太阳辐射数据，但这套日值数据集里没有直接的水平面总辐射。能做的是用日照时数结合气温、湿度去估算。

经典的 Ångström–Prescott 类模型利用"日照时数 / 可照时数"的比值来估算日总辐射，可照时数由纬度和日序算出来。这条路在这套数据上可行，但有几个前提必须满足：日照时数本身得是可靠的、且序列内没有仪器断点；站点必须有日照观测（前面说过，覆盖并不完整）；估算出来的辐射精度通常在10%–20%量级，只适合做区域资源评估，不适合做单个电站的发电量预测。

建筑负荷评估用这套数据相对友好一些，因为它主要靠气温。常用的做法是用日平均气温驱动度日数（HDD/CDD）计算，或者用日最高最低气温构造简化的日温度曲线，再配合湿度算焓值。这里最需要注意的是城市站和郊区站的差别——建筑负荷本身就和城市热岛强相关，用郊区站的数据做城市建筑能耗评估会偏低，用城市站做区域评估会偏高。选站的时候要带着目的选，不能图省事用"离项目最近的站"。

注意：如果项目所在地周边50公里内没有国家级站，而这套数据里又只有国家级站，那么用最近站做空间外推的误差可能非常大，尤其在山区和沿海。这种情况下应该明确说明数据代表性不足，而不是硬套一个插值结果。

最后分享一个我自己踩过好几次、后来养成习惯的小动作：每次拿到一份新下载的数据，先花十分钟做一份"数据日志"——记录下载日期、文件总数、总记录数、站点数、年份范围、各要素缺测率Top10的站、异常值过滤各阶段剔除的行数。这份日志放在项目目录里，半年后回头看，能立刻判断出当时的数据处理有没有出错，也能在发现问题时快速定位是哪一批文件引入的。看起来是额外工作，实际上每次都能省下更多排查时间。数据这东西，你不对它做记录，它就会在你最需要它的时候给你制造麻烦。