简介一套基于Google Earth Engine与NASA GPM IMERG V07降水数据的极端降水分析代码包面向洪涝灾害风险评估、流域降水过程及气候变化研究场景适合从事水文气象、环境遥感等方向且有一定GEE基础的开发者或研究者。压缩包仅6KB共3个文件分别承载分析主代码、可视化展示页面和版本管理配置结构轻量可直接导入GEE环境运行。代码覆盖研究区加载、数据预处理、日降水构造、阈值天数统计、三天滑动窗口累计与峰值日期提取等核心环节运行后能输出平均日降水、最大日降水、极端降水天数等地图图层和图表并附参数调整与ROI资产准备说明便于按研究需求修改阈值和区域整体流程完整可复现从输入数据到结果导出的全过程。目前已有120人学习浏览可作为快速上手GEE极端降水分析的参考实现。1. GEE 极端降水分析为什么不先下载数据再本地算2021 年河南“7·20”暴雨之后我接过一个复盘任务把那次极端降水过程的强降雨中心、三小时雨强峰值和累计雨量算出来给防汛复盘会出图。传统做法是先下载 GPM 数据再本地解压、重投影、拼接、裁剪最后还要跟站点数据做交叉验证一套流程下来小半天就没了。后来我把整套分析直接搬到了 GEEGoogle 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 官方元数据数据源GEE 集合 ID空间分辨率时间范围更新滞后适合场景TRMM 3B42TRMM/3B420.25°1998–2019已停止2019 年以前旧事件回算GPM IMERG FinalGPM/IMERG/V06/FINAL0.1°2000-06 至今约 14 天暴雨事件精细复盘ERA5-Land DailyECMWF/ERA5_LAND/DAILY_AGGR0.1°1950 至今约 2–3 个月长序列气候态统计CHIRPS DailyUCSB-CHG/CHIRPS/DAILY0.05°1981 至今约 3 周站点稀疏区交叉验证四个集合里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 FinalCHIRPS 作为校验源你可以在参数区直接切换两个集合 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 用 10000ERA5 用 11132Sentinel-2 用 10。如果影像还是太大再配合maxPixels: 1e10兜底。这条经验我几乎每次跑大区域都要用上。5.4 现象历史气候态算出来全是 0原因多半是历史日期范围写到了 GPM 数据还没开始的时间。GPM IMERG 在 GEE 里最早是 2000 年 6 月如果你把历史窗口写成 1990 到 2000 年filter 出来的集合是空的reduce 结果自然全是 0。解决方法是查清楚每个数据集的存档起止时间再写窗口。TRMM 覆盖 1998 到 2019IMERG 覆盖 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 做校验。希望帮到你。本文还有配套的精品资源点击获取