第一次把 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 或者反过来再滚动数组让区域连续。时间戳更容易出错而且错了不容易发现。格点数据的时间标签通常标的是时段的起点或者终点二者差一个时间步不同版本、不同产品不一样。这个差别在日尺度上不明显但在做日累计的时候会直接导致错位如果标签是时段终点你把 t00:00 到 t21: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 个有效时次才输出否则返回 NaNxarray 里可以用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 小时时间戳标记的是时段终点整体平移一个时间步高纬区域降水明显偏高用了算术平均而非余弦加权改为面积加权暴雨量级普遍偏小融合过程平滑极值接受或用分位数映射校正冬季降水明显偏低固态降水反演困难结合温度判断谨慎使用最后分享一个我自己一直在用的小技巧拿到任何一份新的格点降水数据别急着做分析先挑三个格点——一个在站点密集的平原、一个在山区、一个在海上是——把它们的年际序列和季节循环画出来跟常识对一对。平原站点应该和本地站点吻合得不错山区应该有明显的地形增幅海上应该平滑且季节信号清晰。这三张图能在半小时内告诉你这份数据大概靠不靠谱比跑一堆指标来得快。等到图看着不对劲了再回头去查单位、时区、聚合方式那才是系统地排查。