GEE极端降水分析：从数据源选型到NDWI淹没范围

简介：一套基于Google Earth Engine与NASA GPM IMERG V07降水数据的极端降水分析代码包，面向洪涝灾害风险评估、流域降水过程及气候变化研究场景，适合从事水文气象、环境遥感等方向且有一定GEE基础的开发者或研究者。压缩包仅6KB，共3个文件，分别承载分析主代码、可视化展示页面和版本管理配置，结构轻量，可直接导入GEE环境运行。代码覆盖研究区加载、数据预处理、日降水构造、阈值天数统计、三天滑动窗口累计与峰值日期提取等核心环节，运行后能输出平均日降水、最大日降水、极端降水天数等地图图层和图表，并附参数调整与ROI资产准备说明，便于按研究需求修改阈值和区域；整体流程完整，可复现从输入数据到结果导出的全过程。目前已有120人学习浏览，可作为快速上手GEE极端降水分析的参考实现。

1. GEE 极端降水分析：为什么不先下载数据再本地算

2021 年河南“7·20”暴雨之后，我接过一个复盘任务：把那次极端降水过程的强降雨中心、三小时雨强峰值和累计雨量算出来，给防汛复盘会出图。传统做法是先下载 GPM 数据，再本地解压、重投影、拼接、裁剪，最后还要跟站点数据做交叉验证，一套流程下来小半天就没了。后来我把整套分析直接搬到了 GEE（Google Earth Engine）上，从加载降水数据到输出极端降水日数和淹没范围统计，压缩到十分钟级别。这份 GEE 极端降水分析代码包，就是当时这套流程的完整沉淀：覆盖逐日降水合成、95 分位极端阈值划定、极端降水指数计算、NDWI 水体提取和淹没面积统计，适合水文、气象、遥感和防灾减灾方向的从业者直接改参数复现。

2. 数据源怎么选：GPM IMERG 与 ERA5-Land 的取舍

2.1 为什么直接上 GEE，而不是继续用本地下载的老流程

本地下载的处理链路并不算复杂，但每次极端降水复盘，你同时要准备降水、土地利用、高程、行政边界和卫星影像。GEE 最大的优势是这些数据已经在同一套坐标系里对齐，前端 Console 可见即可得，不用手动管理下载失败的断点。注册入口只需要一个 Gmail 账号，进入 Code Editor 后界面分成脚本区、Console 区、Inspector 区和地图区，这些面板就是整个分析的主战场。

很多教程里提到的 GEE Admin，其实指的是左侧 Assets 面板：你把自定义矢量或栅格上传之后，管理和共享权限都在这里。新建脚本后，右上角的搜索框可以直接搜数据集，输入 GPM/IMERG 会同时列出 Early、Late、Final 三个版本，这就是你后面做实时监测和事件复盘的数据中心。新手最容易忽略 Inspector 面板：点地图上的像元，能直接看到该位置的波段值，判断数据有没有加载成功，比反复 print 影像靠谱得多。

2.2 四个常用降水数据源的参数对比

极端降水分析里，数据源选错比参数调错更致命。以下是我在不同任务里反复对比后保留的四个候选集合，参数全部来自 GEE 官方元数据：

四个集合里，IMERG 的 precipitationCal 波段单位是 mm/hr，半小时一个时相；ERA5-Land 的 total_precipitation 单位是米；CHIRPS 的 precipitation 单位是 mm/day。单位不统一是交叉验证时最容易翻车的点，后面脚本里我会专门讲怎么换算。

2.3 选型逻辑：先回答三个问题再动手

我一般不会一上来就选数据源，而是先问三个问题：要复盘单次暴雨还是算气候态？需要近实时结果还是允许滞后？区域里有没有可靠的站点观测？这三个问题的答案基本能锁定数据源。单次暴雨复盘用 GPM IMERG Final，近实时监测用 IMERG Early，长序列年最大日雨量统计用 ERA5-Land Daily，干旱半干旱区跟站点数据互检用 CHIRPS。

还有一个不成文的习惯：任何极端降水过程都不建议只用单一数据源。IMERG 的卫星反演在复杂地形区有系统偏差，CHIRPS 在站点稀疏区又过度依赖插值，两个源算出来的过程总量如果差超过 30%，我会回头找当地气象站的人工雨量做基准。这套代码包默认主源是 GPM IMERG Final，CHIRPS 作为校验源，你可以在参数区直接切换两个集合 ID。

3. 核心脚本：逐日降水合成与 95 分位极端指数

3.1 逐日降水合成：listSequence + map 的写法

GPM IMERG 半小时一个时相，一天 48 个文件。直接把整段事件期的时相丢进 reducer 求和，内存压力大且不好排查异常景。我习惯先把每个日期的所有时相各自合成一景日降水影像，用 ee.List.sequence 加 map 实现。下面这段代码是代码包里最核心的预处理片段：

// 选择研究区和事件期 var roi = ee.Geometry.Polygon( [[[113.5, 34.0], [114.5, 34.0], [114.5, 35.0], [113.5, 35.0]]] ); var start = ee.Date('2021-07-17'); var end = ee.Date('2021-07-25'); // 加载 IMERG Final，过滤到研究区 var imerg = ee.ImageCollection('GPM/IMERG/V06/FINAL') .filterBounds(roi) .filterDate(start, end) .select('precipitationCal'); // 以天为单位生成序列，逐日合成降水 var days = ee.List.sequence(0, 7).map(function(n) { var day = start.advance(n, 'day'); var dayEnd = day.advance(1, 'day'); var daySum = imerg.filterDate(day, dayEnd) .sum() .multiply(0.5) .rename('precip_mm'); return daySum.set('system:time_start', day.millis()); }); var daily = ee.ImageCollection(days); print(daily);

这段代码的逻辑是先把日期偏移量 0 到 7 映射成 8 个日期，再对每个日期做一次 filterDate。.sum()把该日期的所有时相累加，.multiply(0.5)则是关键：IMERG 的 precipitationCal 单位是 mm/hr，半小时步长乘以 0.5 小时才换算成毫米。.rename('precip_mm')统一了波段名，后面跟阈值影像比较时保证同名波段对齐。

参数上你只需要改 roi、start、end 和 ee.List.sequence 里的天数上限。sequence 的第二个参数是偏移数量，写成 7 表示 2021-07-17 到 2021-07-24 共 8 天，如果事件期是 10 天，这里要改成 9。我经常因为只改了日期没改天数上限，导致事件期最后一天缺失，这种低级问题最容易在复盘汇报时被当场指出来。

3.2 95 分位阈值：先算气候态，再换事件期

极端降水指数里，R95P 指的是日降水超过气候态 95 分位阈值的那部分降水总量。算 R95P 必须先有一条历史序列做基准。这里有个数据边界要特别注意：IMERG Final 在 GEE 里的 V06 集合最早只能回溯到 2000 年 6 月，所以气候态窗口不能写 1998。我会用 2001 到 2020 这 20 年做历史基准，避免边界年份数据不完整：

// 历史期逐日降水，先合成再求 95 分位 var histStart = '2001-01-01'; var histEnd = '2020-12-31'; var histDaily = ee.ImageCollection('GPM/IMERG/V06/FINAL') .filterBounds(roi) .filterDate(histStart, histEnd) .select('precipitationCal') .map(function(img) { // 半小时数据转日累计，mm/hr * 0.5h return img.multiply(0.5).rename('precip_mm'); }); // 这里为了控制计算量，直接用时相求百分位做近似 var p95 = histDaily.reduce(ee.Reducer.percentile([95])) .rename('p95_threshold'); print('95th percentile threshold: ', p95);

严格做法是先按 3.1 的逐日合成思路把 20 年历史数据全部合成日序列，再对日序列求百分位。但在 GEE 免费配额下，20 年逐日合成非常耗时，所以这个脚本默认直接用原始时相求百分位，算出来的阈值会比标准 R95P 略微偏低。如果你的区域极端降水事件集中在一两个月份，建议把历史窗口改成对应季节，比如只取 6 到 8 月，这样阈值更贴近汛期实际。

p95 输出是一景栅格影像，不是全局数值。不同地理位置的 95 分位阈值不一样，山区和盆地的阈值可能相差一倍以上，这也是极端降水分析必须做栅格化阈值而不是用站点单点阈值的原因。

3.3 极端事件识别与导出：scale 参数决定成败

有了日降水序列和 p95 阈值影像，识别极端事件就很直接了。逐日影像跟阈值影像做比较，大于阈值记为 1，然后累加得到极端降水天数。这里要保证两个影像的波段名一致，GEE 的gte比较是按同名波段逐像元运算的：

// 逐日与阈值比较，统计极端降水总天数 var extremeFlag = daily.map(function(img) { var flag = img.gte(p95).rename('extreme_flag'); return flag.set('system:time_start', img.get('system:time_start')); }); var extremeDays = extremeFlag.sum().rename('extreme_days'); // 导出到 Drive，注意 scale 不要小于数据源原始分辨率 Export.image.toDrive({ image: extremeDays.float(), description: 'extreme_precip_days_202107', folder: 'GEE_extreme_export', region: roi, scale: 10000, crs: 'EPSG:4326', maxPixels: 1e10 });

导出的 scale 参数我设置成 10000 米，正好对应 IMERG 0.1° 的原始分辨率。把 scale 调到 1000 米并不会带来更高精度，只会让像元数量暴增一百倍，触发 GEE 内存超限。如果你的研究区跨多个经纬度带，crs 建议换成区域投影而不是默认的 EPSG:4326，比如中国中东部用 UTM 50N，这样导出后的面积统计更准确。

4. NDWI 水体联动：暴雨后的淹没范围到底在哪

4.1 为什么极端降水分析要叠加 NDWI

降水数据只能告诉你“下了多少雨”，淹没范围则需要另一套证据链。暴雨过后，淹没区水体在可见光绿色波段反射率升高，近红外波段反射率骤降，NDWI 能把这个差异拉满。NDWI 的计算公式是 (Green - NIR) / (Green + NIR)，水体像元通常大于 0，而土壤和植被像元明显小于 0。代码包里默认用 Sentinel-2 的 B3 绿波段和 B8 近红外波段做 NDWI，空间分辨率 10 米，比降水的 10 公里网格精细得多。

但 NDWI 不是银弹。暴雨过程往往伴随云覆盖，而 Sentinel-2 重访周期约 5 天，最好的观测窗口是雨后 1 到 3 天、云量低于 30% 的景。如果窗口内云太多，NDWI 会漏掉大片水体。另一个坑是山体阴影和高层建筑阴影的 NDWI 值容易与水体重叠，所以阈值不能照抄论文里的固定值。

4.2 云掩膜与 NDWI 计算

代码包里已经写好了 Sentinel-2 云掩膜函数。S2 SR 产品的 QA60 波段用位编码记录云和卷云信息，处理时需要按位判断：

// Sentinel-2 云掩膜函数 function maskS2clouds(image) { var qa = image.select('QA60'); var cloudBitMask = 1 << 10; var cirrusBitMask = 1 << 11; var mask = qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(cirrusBitMask).eq(0)); return image.updateMask(mask) .copyProperties(image, ['system:time_start']); } // 加载雨后的 Sentinel-2 影像 var s2 = ee.ImageCollection('COPERNICUS/S2_SR') .filterBounds(roi) .filterDate('2021-07-22', '2021-07-25') .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 70)) .map(maskS2clouds); // 计算 NDWI 并取最大值合成，降低单景噪声 var ndwi = s2.map(function(img) { return img.normalizedDifference(['B3', 'B8']).rename('NDWI'); }).max().clip(roi); var water = ndwi.gt(0.2); Map.addLayer(water, {palette: ['blue']}, 'Flood Water');

位掩膜代码里的1 << 10和1 << 11分别对应 QA60 波段第 10 和第 11 位。把云和卷云的位值清零后，剩余像元才参与 NDWI 计算。.filterDate('2021-07-22', '2021-07-25')是刻意选的雨后窗口：7 月 20 日强降雨发生后，22 日到 25 日之间影像既包括刚退去的积水，又避开了降雨当天的密云。

NDWI 阈值 0.2 是代码包的默认值，适合平原河网区。山区建议先做一次椒盐试验：把阈值从 0 到 0.4 每次加 0.05，分别统计水体面积，画出面积变化曲线，找到拐点对应的阈值。这个拐点通常就是区分水体和阴影的最佳值。

4.3 淹没面积统计与逐日变化

水体面积统计不能直接用reduceRegion的 sum reducer 数像元个数，因为不同纬度像元实际面积不一样。正确做法是把水体掩膜乘以ee.Image.pixelArea()，让每个像元乘上自己的面积，再求和：

// 统计淹没面积，单位从平方米转为平方公里 var areaKm2 = water.multiply(ee.Image.pixelArea()) .reduceRegion({ reducer: ee.Reducer.sum(), geometry: roi, scale: 10, maxPixels: 1e13 }) .getNumber('NDWI') .divide(1e6); print('Flood area (km2):', areaKm2);

如果你需要逐日淹没面积变化，不要每天单独导出一张影像再在 GIS 里拼。代码包里提供了批量统计模板：把 Sentinel-2 影像集合按日期分组，每组算一个水体面积，输出一张时间序列表。我一般会把 7 月 22 日到 30 日的面积曲线和降水柱状图叠在一起看，面积峰值一般会比降雨峰值滞后 12 到 24 小时，这个滞后量对洪水预警有直接参考价值。

5. 避坑专题：GEE 极端降水分析最常见的五个翻车点

5.1 现象：事件当天在 IMERG Final 里取不到数据

你用filterDate('2021-07-20', '2021-07-21')查 IMERG Final，结果某一天影像数是零。这不是代码 bug，而是 IMERG Final 产品需要融合地面站点观测，在 GEE 上的发布滞后大约 14 天。事件刚发生就去拉 Final，必然空手而归。解决方法是：实时监测用 EARLY 或 LATE 集合，事件结束两到三周后再用 FINAL 做复盘级分析。EARLY 滞后约 4 小时，LATE 滞后约 12 小时，两者精度差距不大，但监测场景足够用。

5.2 现象：ERA5-Land 的 total_precipitation 数值小得离谱

ERA5-Land Daily 的降水单位是米，不是毫米。同一场 50 毫米的暴雨，在这个数据集里是 0.05。如果直接拿去跟 IMERG 的毫米值比较，你会得出“这场雨只有 5 厘米”的荒谬结论。原因就是单位换算没做。解决方法是读取集合前先看 bands 信息，total_precipitation统一乘以 1000 转成毫米再入后续流程。这是整个代码包里我唯一建议写死单位转换的地方。

5.3 现象：reduceRegion 报错内存超限

报错信息通常是Computed image is too large。原因有两类：一是研究区跨度过大，二是 scale 设得比数据源原始分辨率小很多。比如用 IMERG 数据却设 scale 100 米，GEE 就得在一个本来 10 公里分辨率的像元里强行细分出无数个伪像元。解决方法是先把 scale 对齐到数据源原始分辨率，IMERG 用 10000，ERA5 用 11132，Sentinel-2 用 10。如果影像还是太大，再配合maxPixels: 1e10兜底。这条经验我几乎每次跑大区域都要用上。

5.4 现象：历史气候态算出来全是 0

原因多半是历史日期范围写到了 GPM 数据还没开始的时间。GPM IMERG 在 GEE 里最早是 2000 年 6 月，如果你把历史窗口写成 1990 到 2000 年，filter 出来的集合是空的，reduce 结果自然全是 0。解决方法是查清楚每个数据集的存档起止时间再写窗口。TRMM 覆盖 1998 到 2019，IMERG 覆盖 2000 年 6 月至今，两者中间有重叠期，但绝对不能混着当同一套气候基准。

5.5 现象：NDWI 水体提取结果里全是斑点

原因是雨后地表湿润，土壤和低矮植被的 NDWI 也会短暂升高，加上 Sentinel-2 影像受云掩膜影响，水体边界破碎。解决方法是加两步后处理：第一步用focal_min加focal_max做形态学闭运算，填补水体内部空洞；第二步用connectedComponents加面积滤波，删除小于 30 万平方米的碎块。水面面积阈值可以按研究区灵活调，但要记住一个原则：宁可漏掉小水塘，也不要让噪声淹没真正的主河道淹没范围。

6. 可复现模板：把参数抽出来以后批量导出

这套代码包用得越久，我越觉得参数抽离比算法本身更重要。极端降水分析有一个明显特点：每次任务的日期、区域、阈值都不同，但处理链完全一致。我把代码包里的核心处理链封装成一个模板函数，每次新任务只改函数入参：

function runExtremePrecipAnalysis(startDate, endDate, roi, scale, thresholdMode) { var imerg = ee.ImageCollection('GPM/IMERG/V06/FINAL') .filterBounds(roi) .filterDate(startDate, endDate) .select('precipitationCal'); // 逐日合成逻辑与第 3 章一致，这里精简为直接取日累计 var daily = ee.ImageCollection( ee.List.sequence(0, ee.Date(endDate).difference(ee.Date(startDate), 'day')) .map(function(n) { var day = ee.Date(startDate).advance(n, 'day'); return imerg.filterDate(day, day.advance(1, 'day')) .sum().multiply(0.5).rename('precip_mm'); }) ); // 阈值模式：clim 用历史 95 分位，event 用固定阈值 var threshold; if (thresholdMode === 'clim') { threshold = ee.Image('projects/your_asset/p95_2001_2020'); } else { threshold = ee.Image.constant(50); } var extremeDays = daily.map(function(img) { return img.gte(threshold).rename('extreme_flag'); }).sum(); Export.image.toDrive({ image: extremeDays.float(), description: 'extreme_days_' + startDate, region: roi, scale: scale, crs: 'EPSG:4326', maxPixels: 1e10 }); }

模板函数里的thresholdMode是核心设计：复盘一次具体暴雨过程时，固定阈值 50 毫米比气候态 95 分位更直观；做区域极端性评估时又切回气候态模式。我建议你每次接到新任务，都强制走一遍这套流程：先传一个最小 roi 做 dry-run，打印出daily.size()确认影像数量与预期天数一致，再查看p95栅格的最小值和最大值是否在合理范围，最后才全量导出。这一步能拦住 80% 的低级错误。

另外，给未来的自己留个“后悔药”：所有 Task 导出任务的 description 里必须带日期，比如extreme_days_20210720_v2。GEE 的 Task 列表是按时间滚动的，同名任务一旦重复提交，你根本分不清哪个是最新结果。我在一次连续分析了三个台风事件时，因为两个任务的 description 都叫extreme_days，导出后完全对不上号，白跑了一下午。从那以后，我每次都强制在 description 里拼上日期和版本号，并在导出前用Export.table.toDrive导一份像元统计 CSV 做校验。希望帮到你。

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