第一次把 3小时尺度降水数据集-MSWEPV2 解压出来的时候,我盯着那个 0.1° 的全球网格看了好一会儿——一个格点大约 11 公里见方,每 3 小时一个值,从 1979 年一路铺到今天。对做流域水文、城市内涝或者新能源功率预测的人来说,这几乎等于把"雨到底什么时候下、下在哪儿"这个问题摊开在桌面上。但真把它接进模型、跑完第一轮验证之后,我发现多数人对这类数据的期待是错位的:它不是测雨雷达,更不是雨量站,它是一份经过多源融合的背景降水场,用得好是骨架,用不好就是一堆看起来很精致的假精度。
这篇东西不讲虚的,就讲我从下载、读取、聚合、验证到批量工程化这条链路上踩过的坑和总结出来的口径。适合三类人看:手上有一份格点降水数据但不确定该怎么用的人、需要给水文模型或作物模型配驱动的人、以及被"3 小时"这个时间分辨率吸引但还没想清楚它意味着什么的人。基础知识我会顺带补,不会假定你熟悉 NetCDF 或者 xarray。
1. 3 小时这个步长,到底解决了日尺度的什么问题
1.1 日降水把洪峰"摊平"了,这是换数据的直接动机
早几年我用日尺度的格点降水驱动过一个山区小流域的产汇流模型,结果非常别扭:模拟出来的洪水过程线永远是"胖"的,峰值偏低、峰现时间偏晚,而且不管怎么调参数都调不回来。后来想明白了,问题不在模型,在输入。日降水是 24 小时的累计量,把它平均分配到 24 个小时里,等于人为假设这场雨从早到晚匀速下——可现实里一场对流雨可能两个小时就下完了 40 毫米,剩下 22 小时是毛毛雨甚至不下。
这两种情形在日总量上完全一样,但产流机制天差地别。短历时高强度降水更容易超过下渗能力、形成地表径流,洪峰又尖又陡;而均匀小雨大部分被土壤吸收,产流系数低。你拿一个被抹平的降水序列去驱动非线性模型,模型只能通过调整参数来"补偿",最后参数被调得面目全非,换一场雨就崩。
3 小时步长把这个问题缓解了一大半。它依然不够细——真正的城市排水设计用的是 5 分钟到 1 小时的雨强——但对中小流域的洪水预报、日内的灌溉调度、风电光伏的短时出力波动评估来说,3 小时是一个性价比很高的折中点:既能捕捉到降水过程的形态和峰现时间,又不至于因为时间步太细而让数据噪声和缺测问题放大。
1.2 多源加权集合,加的是三类源的"互补性"
理解这份数据之前,得先接受一个前提:全球任何单一观测手段都测不准降水。雨量站最准,但它只代表一个点,海上和高原基本是空白;卫星覆盖广,但红外通道反演的是云顶温度,跟地面雨强的物理关系很间接,微波通道好一些但 revisit 周期有限;再分析产品时空连续、变量齐全,可它的降水本质上是模式算出来的,强对流经常被平滑掉。
这套数据集的做法,是把雨量计、卫星、再分析这三类源放在一起,用一个时空变化的权重去合成。权重的逻辑大致是:在站点密集的地区和时段,给雨量计更大的话语权;在站点稀疏的海洋和山区,转向卫星和再分析;同时根据每个源在历史时段内与"真值"(站点观测)的误差特征,动态调整它在不同月份、不同区域的分量。公开的技术文档里能看到这种"分源、分时空权重"的思路,具体算法在不同版本间有调整,用之前建议把对应版本的说明文档翻一遍。
这里有个非常重要的推论,也是很多人忽略的:这套产品的日总量比它的 3 小时分布更可靠。因为在融合过程中,日尺度的约束更强,而 3 小时尺度上的时间分配往往要借助卫星或再分析的日内形态来做降尺度。所以如果你用它算年降水量、月水资源量,放心用;如果你用它去论证"某地某天下午 3 点出现了 50 年一遇的短历时强降水",那就要非常谨慎,最好有本地站点或雷达佐证。
1.3 这套组合的适用边界,比参数本身更值得记
我一般会把它能干的活分三档。
第一档是"放着放心"的场景:流域尺度的日径流模拟、月度水资源评价、干旱指数(SPI/SPEI)计算、作物生长期的降水盈亏分析、大区域的水量平衡核算。这些场景对单点误差不敏感,对空间连续性和长时间序列的完整性要求高,正是融合产品的强项。
第二档是"能用但要验证"的场景:中小流域的洪水预警、城市内涝的风险区划初筛、水电和风电的来水/出力中长期评估、农业灌溉的旬尺度调度。这些场景开始关心降水的时序形态和极值,必须拿本地站点做一轮对比,看偏差有多大、能不能接受。
第三档是"别指望它"的场景:城市排水管网设计(需要分钟级雨强和设计暴雨公式)、单沟单坡的滑坡临界雨量阈值(需要站点级别的实时雨量)、短时临近预报(需要雷达外推)。这几个场景对时间和空间精度的要求远超 3 小时 / 0.1°,硬用会出事。
记住一句话:这份数据是"面上看得清、点上不一定准"。它的价值在空间格局和长时间序列,不在单点极值的精确刻画。
2. 文件到手之后,读取环节就有三个硬坑
2.1 先摸清文件组织和变量维度,再写代码
拿到数据包,别急着写循环。先随便挑一个文件用命令行工具看一眼结构,把下面这几件事确认清楚:文件是按年组织还是按月组织、变量名是什么、维度顺序是 (time, lat, lon) 还是别的、经纬度是递增还是递减、时间轴的单位和起点是什么。
这些看起来琐碎,但任何一个搞错,后面全是错的。我见过最典型的两种情况:一是纬度从北往南递减(90 到 -90),直接按数组下标做区域切片,切出来的正好是南北颠倒的一块;二是经度用 0–360 表示,而你的研究区边界文件是 -180–180,两者相减之后出现负值或者跳变,裁剪结果只剩半个区域。
比较稳妥的做法是先写一小段探查代码,把 shape、坐标范围、第一个和最后一个时间戳、变量属性全部打印出来,人工确认一遍。这一步花十分钟,能省后面十个小时。对 NetCDF 来说,xarray.open_dataset打开之后直接看 repr 就能看到维度、坐标和属性,非常省事。
2.2 short 型存储背后的换算,换算错了整套数据都废
这类全球长序列数据为了压体积,几乎都用 2 字节整型存储,而不是浮点。也就是说,文件里存的不是直接的降水值,而是一个整数编码,真实物理量要经过类似value * scale_factor + add_offset的换算,缺测位置则用一个约定的填充值(常见的是 -32768 这类极小值)标记。
这意味着如果你直接把原始数组当降水量用,会得到一堆莫名其妙的负数或者巨大的数。正确姿势是让读取库自动处理:xarray 在打开时会识别scale_factor、add_offset和_FillValue属性,并自动做解码(除非你手动关掉了 mask_and_scale)。但自动解码不是万能的,我遇到过填充值属性缺失、或者属性名不规范导致解码失败的情况,这时候就得手动设。
还有一个细节值得单独说:这份数据里会出现少量负值。这不是 bug,是插值和融合算法在极干燥地区留下的数值副产品,量级通常很小。处理方式很简单,做聚合之前先做一个下界截断,把负值置为 0。如果不做这一步,你在做月累计的时候会看到某个格点的月降水比周边低好几毫米,查半天查不出原因。
2.3 经度范围和时间戳,是两个一上手就会错的地方
经度的事刚才提了,这里补一个更隐蔽的:当一个区域的经度跨越 180° 时(比如太平洋上的岛国或者跨日期变更线的区域),你按 -180–180 去切,会切出两块互不相连的区域,中间空洞。解决办法是先统一到一个表示体系,必要时做一次经度回绕(把 0–360 减去 360 或者反过来),再滚动数组让区域连续。
时间戳更容易出错,而且错了不容易发现。格点数据的时间标签通常标的是时段的起点或者终点,二者差一个时间步,不同版本、不同产品不一样。这个差别在日尺度上不明显,但在做日累计的时候会直接导致错位:如果标签是时段终点,你把 t=00:00 到 t=21:00 的八个时次加起来,得到的其实是前一天 21:00 到当天 21:00 的雨量,整体偏移 3 小时。
我的处理办法是找一次已知的强降水过程——最好本地有逐小时站点记录——把格点序列和站点序列并排画出来,看峰值落在哪个时次,偏移多少一目了然。这个校验花不了多少时间,但能让后面所有分析都站得住。
3. 从原始格点到"能进模型的结果",中间有四道工序
3.1 单位是毫米每 3 小时,不是毫米每小时
这是最容易犯、后果最严重的一个错误。这份数据的每个值代表的是该 3 小时时段内的累计降水量,单位是毫米,不是降水强度。如果你把它当强度去做乘法(比如乘以 3 算日总量),结果会偏大三倍;反过来如果你把它当强度直接送进模型而模型期望的是时段累计量,结果会偏小。
判断单位有个简单粗暴的办法:看数值范围。全球范围内单格点单时次的降水,量级一般在 0 到几十毫米之间,极端情况可能上百。如果是强度(mm/h),典型值应该在 0 到十几之间,很少超过 50。拿这个量级去对照一下,基本能判断出单位。
在代码里我习惯在读入之后立刻把语义写进变量名或属性,比如把变量重命名成precip_mm_3h,这样后面任何人接手都不会误解。这种命名习惯看着啰嗦,但在多人协作的项目里能救命。
3.2 聚合顺序:先截断,再求和,最后判缺测
把 3 小时数据聚合成日、旬、月,看似一个resample就完事,实际有讲究。我总结的顺序是三步。
第一步,负值截断,前面说过。第二步,求和。这里的关键是求和函数对缺测的处理方式。默认情况下,很多库在求和时会跳过 NaN,把剩下的有效值加起来——这在降水里是灾难性的:一天 8 个时次里如果缺了 2 个,你得到的是一个偏低的日总量,而且你根本看不出来它偏低。正确的做法是显式要求"必须凑齐 8 个有效时次才输出,否则返回 NaN",xarray 里可以用min_count参数控制。第三步,标记缺测,把 NaN 数量超过阈值的时段单独打上质量标记,后面统计时排除。
日变月、月变年同理。年尺度上还要注意闰年,2 月 29 日这一天在按"每年 365 天"硬编码的处理流程里会被悄悄丢掉,一年少一天的降水,在湿润地区可能有几十毫米的偏差。
3.3 时区转换:做日变化分析绕不过去的一道坎
格点数据的时间一律是 UTC,这点没有例外。而降水有明显的日变化节律——很多地区对流雨集中在午后到傍晚,山区可能有夜雨现象。如果你直接用 UTC 去做"逐时次平均降水量"的图,得到的峰值时间是 UTC 的峰值,换算成地方时可能要加减好几个时区,结论完全反过来。
具体操作上,先算出研究区中心经度对应的时区偏移(经度除以 15,四舍五入),然后在聚合之前把时间轴整体平移。注意平移之后会出现"日界不整齐"的问题:平移 8 小时之后,原本 00:00–21:00 的八个时次不再对应一个自然日。处理办法是要么接受这种偏移(做日变化分析时没关系),要么在平移到小时级之后再重新按当地日切分(做日总量时必须这样)。
我做日变化分析的习惯是:先把 3 小时序列插值或展开到小时,再按地方时分组统计,最后画 24 小时的平均降水强度曲线。这样得到的峰值时间可以直接和本地气象站的观测对照,说服力强很多。
3.4 面雨量的 cos(lat) 加权,以及流域掩膜的正确做法
把格点降水平均到一个流域,最直觉的做法是算术平均,但这是错的。等经纬度网格上,每个格点代表的实际面积随纬度变化,高纬度的格点面积远小于赤道附近。全球平均用算术平均会严重高估高纬降水的权重。正确做法是按纬度的余弦加权,或者更严格地按格点实际面积加权。
流域掩膜也有细节。常见做法是拿流域边界去裁剪,落在边界内的格点赋 1,边界外赋 0。但 0.1° 网格和流域边界往往对不齐,边缘格点可能只有一部分在流域内。对小流域(面积小于几个格点)来说,这种"全有全无"的处理会造成很大的面积误差。更稳妥的是按覆盖比例加权,也就是算每个格点被流域覆盖的面积占比,用它做权重。
还有一个经常被忽略的点:流域面积和降水量相乘得到的是体积,单位要统一。1 毫米降水落在 1 平方公里上,是 1000 立方米。这个换算在水文模型率定和水量平衡校核里天天用,但每年都有人在这里错一个数量级。
4. 验证这一步,决定了这份数据在你的区域能不能用
4.1 连续型指标为什么会"骗人"
拿到数据和站点做对比,很多人第一反应是算相关系数和 RMSE。这两个指标在降水上都有明显的局限。
降水是高度偏态分布,大量时次是 0,少数时次是几十毫米。这种情况下相关系数会被那些"共同为 0"的时次撑得很高——你哪怕预测全是 0,相关系数也可能有 0.6 以上,但实际毫无价值。RMSE 同理,它对大值极度敏感,一场暴雨的偏差就能主导整个指标,让你看不清日常降水的表现。
我的做法是组合使用。除了相关系数和 RMSE,一定加上偏差百分比(PBIAS),看总量是系统性偏高还是偏低;再加上 KGE 及其三个分解项(相关性、变率比、偏差比),能区分"形态像但幅度不对"和"幅度对但形态不像"这两种完全不同的错误。变率比小于 1 说明数据把降水的波动压平了,这在融合产品里很常见,尤其是在强对流频发的地区。
4.2 阈值化分类检验:POD、FAR、CSI 怎么读
对预警类应用来说,"有没有下雨"比"下了多少"更重要。这时候要做阈值化检验:设一个阈值(比如 3 小时 5 毫米、10 毫米、25 毫米),把格点和站点都二值化,然后算命中率、空报率和临界成功指数。
| 指标 | 含义 | 偏高说明什么 | 偏低说明什么 |
|---|---|---|---|
| POD(命中率) | 实际有雨且被报出的比例 | 数据对降水事件敏感 | 漏报多,小雨被平滑掉 |
| FAR(空报率) | 报有雨但实际没雨的比例 | 数据"虚胖",毛毛雨过多 | 数据偏保守 |
| CSI(临界成功指数) | 综合考虑命中与空报 | 整体判别能力好 | 判别能力弱 |
| PBIAS(偏差) | 总量偏差百分比 | 系统性高估 | 系统性低估 |
实测经验是:这类融合产品在中雨量级(3 小时 5–15 毫米)表现最好,CSI 通常最高;小雨量级容易空报,因为格点代表的是 11 公里范围的平均,实际落区可能只有其中一小块,站点没测到但格点有值,严格说这不算错,只是尺度不匹配;暴雨量级容易漏报,因为融合过程会把极值平滑掉。所以做检验时一定要分阈值段报告,只报一个总体指标是没有意义的。
4.3 一次山区暴雨的复盘:偏差到底出在哪个环节
说个具体的。某次山区强降水过程,站点记录 6 小时下了 80 多毫米,格点数据在同一时段只给出 40 毫米左右,差了整整一半。一开始我怀疑是单位或者时间对齐的问题,逐项排查完发现都对,那就是数据本身的问题。
往下追,发现三个因素叠加。第一是地形,这个流域高差大,站点在山谷,而 0.1° 格点平均了山谷和山脊,山脊的强降水被稀释;第二是融合权重,这一带站点稀疏,权重更多给了卫星和再分析,而卫星对地形云的降水反演本来就偏弱;第三是时间平滑,日尺度融合结果向下分配到 3 小时时,峰值被削平了。
这个复盘给我的教训是:遇到格点与站点不符,先别急着做偏差校正,先搞清楚偏差的来源。如果是尺度不匹配(格点平均 vs 点观测),校正也校正不好,只能接受;如果是系统性偏差(全区都偏低 20%),那用简单的比例因子就能改善;如果是随机误差,校正反而可能引入新的问题。
4.4 偏差校正:线性缩放之外还有别的选择
如果确认是系统性偏差,最省事的办法是分位数映射或者简单的月尺度比例因子。做法是:按月份统计历史期格点与站点的累计量比值,得到 12 个系数,然后对未来的格点值做缩放。这个方法简单、可解释、不容易过拟合,对水量平衡类应用足够。
更精细一点的是分位数映射:把格点分布的各个分位点对齐到站点分布的对应分位点,这样不仅修正总量,还能修正分布形态,尤其是极值端。但要注意,分位数映射在样本量小的站点上会不稳定,尾部容易出问题,做的时候要对极值段做平滑处理。
我的建议是:先做最简单的比例校正,看效果;如果关键指标(比如洪峰流量)还是不达标,再上分位数映射。别一上来就上复杂的机器学习校正,那玩意儿在样本不足的时候,大概率是把噪声学到了模型里。
5. 跑通单点之后,工程化还有一堆细节
5.1 全球数据本地化:先裁剪再计算,能省掉 90% 的算力
全球 0.1° 的 3 小时数据,一年大概是几个 GB 的量级,四十多年下来是好几百 GB。绝大多数人的研究区只是其中一小块,直接全量读进来再裁剪,既慢又吃内存。
正确的顺序是:先按经纬度范围做切片,再读入内存。如果格式支持,最好在读取阶段就带上范围参数(很多库支持 SEL 或者索引切片,能在解压层面就只读需要的块)。切完之后,一个中等流域的数据量可能只有几十 MB,处理起来飞快。
时间维度上也一样。如果只研究汛期,就只读汛期那几个月的文件;如果只做某十年的分析,就别把四十年的都打开。
5.2 批量处理的并行与内存:按年切块最省心
要处理长序列,我的固定套路是按年循环,每年一个任务,任务内部再按需要的空间范围切片。这样做的好处是内存占用可控、失败可重试、进度可追踪。用多进程并行的时候,进程数不要开太多,一般等于物理核数或者略少,因为瓶颈往往在磁盘 IO 而不是 CPU。
如果是用 dask 这类延迟计算框架,一定要显式设置分块大小(chunk),并且沿着时间维分块,因为大多数聚合操作(求和、平均)是沿时间做的。分块设置不当会导致每个任务都要读一遍整个文件,反而更慢。我踩过这个坑:默认分块下跑一个月的聚合花了两小时,手动指定 chunk 之后降到八分钟。
还有一个容易忽略的点是输出的存储格式。如果中间结果要反复读写,别用 CSV 存格点数据,用 NetCDF 或者压缩的二进制格式,体积能小一个数量级,读写也快。
5.3 版本、元数据与可复现:三个月后你会感谢自己
这类数据集会持续更新和修订,不同版本之间可能有系统性差异。所以第一件事是把用的版本号、下载时间、文件校验信息记下来,写进项目文档。第二件事是在输出的结果文件里保留完整的处理链条元数据:原始版本、裁剪范围、聚合方式、缺测处理规则、是否做过偏差校正、校正系数是多少。
我吃过亏:半年前跑的一套结果,回头想复现,发现当时用的是哪个版本记不清了,聚合时用的min_count是多少也想不起来,只能重跑一遍,还得重新验证。现在我固定写一个处理脚本,把所有参数放在文件顶部的配置区,每次运行自动把配置和版本信息写进输出文件的属性里。三个月后回头看,一目了然。
5.4 常见问题速查
| 现象 | 最可能的原因 | 处理方式 |
|---|---|---|
| 数据里出现负数 | 融合算法的数值副产品 | 聚合前截断为 0 |
| 区域裁剪结果只有一半 | 经度范围 0–360 与 -180–180 不匹配 | 统一经度体系,必要时回绕 |
| 日总量系统性偏低 | 求和时跳过 NaN 或时间戳偏移 | 设置 min_count,校验时间标签 |
| 峰值时间比实际晚 3 小时 | 时间戳标记的是时段终点 | 整体平移一个时间步 |
| 高纬区域降水明显偏高 | 用了算术平均而非余弦加权 | 改为面积加权 |
| 暴雨量级普遍偏小 | 融合过程平滑极值 | 接受,或用分位数映射校正 |
| 冬季降水明显偏低 | 固态降水反演困难 | 结合温度判断,谨慎使用 |
最后分享一个我自己一直在用的小技巧:拿到任何一份新的格点降水数据,别急着做分析,先挑三个格点——一个在站点密集的平原、一个在山区、一个在海上是——把它们的年际序列和季节循环画出来,跟常识对一对。平原站点应该和本地站点吻合得不错,山区应该有明显的地形增幅,海上应该平滑且季节信号清晰。这三张图能在半小时内告诉你这份数据"大概靠不靠谱",比跑一堆指标来得快。等到图看着不对劲了,再回头去查单位、时区、聚合方式,那才是系统地排查。