前阵子接了个活儿某半干旱地区要做地下水资源可持续利用规划第一个问题就是这块地的地下水补给量到底有多少水文站网稀疏监测井资料零散地面实测这条路基本走不通。我给的方案很直接——用GEE拉取多源遥感数据在一个计算框架里跑通水量平衡把补给量估算落到每个像元上。这在几年前想都不敢想但现在GEE让这些事变成了一堆能跑的代码。这篇文章是这次案例的完整复盘从数据选型、代码实现到验证思路都会摊开讲适合刚接触遥感水文、或者已经会点GEE但想往水循环方向延伸的朋友。基于多源遥感数据的地下水补给量估算核心是一个听起来很朴素的等式降水进来一部分蒸散发掉一部分形成径流走掉一部分存在土壤里剩下的才轮到地下水。难点从来不在数学公式上而在每一个分量的数据从哪来、误差有多大、怎么在GEE里把这些不同来源的栅格数据对齐到同一个时间尺度和空间范围上。下面我把整个处理链路完整拆一遍。1. 补给量估算的思路先把自己装进水量平衡这口锅里这里说的地下水补给量指的是降水通过包气带到达潜水面、真正进入地下水体的那部分水量单位一般用毫米或者立方米。估算方法其实不少但我这次选的是水量平衡法原因后面会讲。先把这个方法的核心逻辑说透。1.1 水量平衡方程中的每一项从哪里来把目标含水层之上的包气带想象成一个水桶降水P是进水口蒸散发ET是桶口不断蒸发掉的水汽地表径流R是桶沿溢出的部分桶内水位变化对应土壤水储量变化ΔS而真正漏到下一层桶里的水才是地下水补给量ΔG。写出来就是P ET R ΔS ΔG在多年平均尺度上ΔS可以近似为零方程简化为ΔG ≈ P – ET – R。但你要是做逐月或者逐年的估算ΔS必须保留不然半干旱地区土壤的季节性蓄水变化会全被算成补给量误差大得没法看。遥感数据能提供的恰好就是这些变量降水用CHIRPS或GPM蒸散发用MOD16A2或SSEBop土壤水储量变化用GLDAS径流没有直接的遥感产品通常用径流系数或者区域水文模型近似。地形数据则用来界定汇水范围和研究区边界。也就是说补给量本身没有一个卫星能直接测到它是通过总进账减去所有其他出路反推出来的剩余项所以每一个输入变量的误差都会累积到最终结果上。这也是后面为什么要做敏感性分析的根本原因。再说说为什么不选其他方法。基流分割法需要长序列日径流数据很多中小流域根本没有水文站地下水储量变化法比如用GRACE卫星数据看总水储量趋势空间分辨率太粗只能用于几十万平方公里以上的大区域。相比之下水量平衡法在GEE里数据齐备、流程可复用、空间分辨率可控是最适合案例区遥感估算这个场景的选择。1.2 为什么是GEE而不是传统下载处理流程以前的做法是去各个数据中心注册下载降水、蒸散、土壤水各下载几十个G的文件再在本地用Python做重投影、裁剪、时间聚合。一套流程下来光数据管理就耗掉大半时间而且换一个研究区、换一个时间窗所有步骤全部重跑。GEE把多源遥感数据集合成了统一读取、统一计算的云端影像集合代码写一遍换区域和时间范围就是改几个参数的事。另一个容易被忽视的点是投影与分辨率。多源数据之间空间参考不同CHIRPS是WGS84经纬度网格MOD16A2是MODIS的正弦投影GLDAS又是自己的规则格网。本地处理时这一步要极其小心否则裁剪出来的结果在边界上全是锯齿状伪差。GEE在reduce和导出时会基于你指定的投影自动重采样虽然它也有自己的坑后面专门讲但至少把最常见的对齐错误挡掉了一大半这对于非遥感专业出身的水文工作者来说特别友好。2. 数据选型降水、蒸散、土壤水、地形这四件套怎么配多源遥感数据的多不是越多越好而是每一类变量选一个最合适的产品再留一个备选做交叉验证。我这次搭配下来比较顺手的一套放在下面这个表格里。数据产品提供的变量时间分辨率空间分辨率GEE里需要处理的坑CHIRPS V2.0降水mm/day逐日0.05°基本没有直接用GPM IMERG降水逐月底层0.1°时间序列短2000年后才有MOD16A2蒸散发8天合成500m比例因子0.1注意异常填充值SSEBop蒸散发8天/月1km序列较短尺度较粗GLDAS-2.1 Noah土壤水kg/m²3小时/月0.25°分辨率太粗注意土层深度SRTM DEM地形静态30m用于提取流域边界2.1 降水数据CHIRPS打底、GPM补充降水是补给量计算里最刚性的输入它的精度直接决定最终结果的可信度。我习惯用CHIRPS V2.0作为主数据源原因很实际序列从1981年延续到现在空间分辨率0.05°约5公里时间分辨率逐日单位就是mm/dayGEE里一个ImageCollection直接读不需要任何换算。CHIRPS的底层是站点观测和卫星红外反演的融合在站点稀疏的半干旱地区表现比纯粹的卫星降水产品稳定不少。如果你的研究区更关注短时强降水过程或者时间窗口集中在近十年可以考虑GPM IMERG。GPM空间分辨率约0.1°时间分辨率能到半小时对极端降水事件的捕捉能力更强但序列从2000年才开始做长期年际变化分析会显得短。我的实操建议是把两套都加载出来做相互校验如果月尺度上两套降水相差超过15%先别急着做下一步优先排查是不是研究区边界、异常值或者数据版本出了问题。花半天时间做这个校验能省下后面一周的返工。2.2 蒸散发MOD16A2与SSEBop的选择蒸散发通常是水量平衡方程里最大的漏项。MOD16A2是MODIS蒸散产品空间分辨率500米8天合成GEE里对应的影像集ID是MODIS/006/MOD16A2。它的算法基于Penman-Monteith方程需要植被叶面积指数、反照率、气象再分析数据等输入。这里有个非常关键的坑它输出的ET变量单位是kg/m²/8day比例因子是0.1意味着要先乘0.1才能得到毫米水量然后还要除以8天才能得到日均蒸散发。很多教程里把这个环节漏了算出来的补给量全年都是负的。MOD16A2另一个问题是在植被稀疏区往往高估蒸散在重度云覆盖区域会填充异常值。处理时必须把异常值像元先掩膜掉再参与计算。SSEBop是另一个常用选择基于地表能量平衡的简化方法空间分辨率1公里受植被参数误差干扰更小对小流域尺度的实用性也不错。我的习惯是用MOD16A2做主流程把SSEBop作为敏感性分析的参照物看最终结果对ET的选择是否稳健。ET如果一换产品结果就翻倍那这个区域的补给量只能给区间不能给单点值。2.3 土壤水和径流数据的近似处理土壤水储量变化ΔS用GLDAS-2.1 Noah产品。注意两点第一GLDAS空间分辨率只有0.25°约25公里在小流域尺度上根本提供不了逐像元空间细节强行参与逐像元计算只会引入一块巨大的伪信号。我的做法是把GLDAS在研究区上做区域平均取其逐月背景变化值作为一个均匀项去修正总量。第二GLDAS的土壤含水量单位是kg/m²按定义恰好等于毫米等效水深数值上可以直接当mm用。但有些文献习惯用体积含水量去理解它数值上会觉得怎么这么小容易误判数据质量。地表径流没有靠谱的遥感直接产品。在干旱和半干旱地区产流以超渗产流为主径流系数通常在0.05到0.15之间。我建议根据研究区下垫面性质取一个经验系数把径流从降水里扣除并在报告里明确标注这是模型简化。如果研究区有水文站实测径流拿实测值来替换会好很多但GEE这个环节里只能先做一个带径流系数选项的可调参数后续拿到实测数据再校准。这个简化不丢人关键是你要知道它简化在哪并且把不确定性传达到结论里去。3. GEE里跑通完整流程代码拆开揉碎下面这段是整个流程的核心代码我会拆开讲每一段在干什么、为什么这么写。研究区以某个中型流域为例时间范围取2015到2020年按月输出补给量。3.1 研究区与时间窗的设定GEE的第一步总是研究区。如果你有流域边界矢量直接上传到GEE资产Assets然后通过ee.FeatureCollection读取var roi ee.FeatureCollection(projects/yourname/assets/watershed).geometry(); var startYear 2015; var endYear 2020; var startDate ee.Date.fromYMD(startYear, 1, 1); var endDate ee.Date.fromYMD(endYear, 12, 31);如果没有现成边界可以用SRTM DEM派生。在GEE里用ee.Terrain.hydroflow对DEM做填洼和水流方向计算再按出水口坐标快速提取汇水区。这个方法在地形起伏明显的区域非常快比手动勾绘更符合水文意义。边界确定后导出一次即可后面所有数据都按这个ROI裁剪。3.2 逐月水量平衡的计算逻辑降水部分加载CHIRPS逐日集合var chirps ee.ImageCollection(UCSB-CHG/CHIRPS/DAILY) .filterBounds(roi) .filterDate(startDate, endDate);要得到逐月降水可以用下面的月序列映射。这段代码有一个小技巧先把年份和月份交叉成一组元组然后对每个年-月窗口做一次过滤与求和最后用flatten()把嵌套列表摊平避免ee.ImageCollection里套着ImageCollection的问题。var years ee.List.sequence(startYear, endYear); var months ee.List.sequence(1, 12); var monthlyP ee.ImageCollection( ee.List(years).map(function(y) { return ee.List(months).map(function(m) { var img chirps .filter(ee.Filter.calendarRange(y, y, year)) .filter(ee.Filter.calendarRange(m, m, month)) .sum() .rename(P); return img.set({ year: y, month: m, system:time_start: ee.Date.fromYMD(y, m, 1) }); }); }).flatten() );ET部分类似但多了单位和合成周期的处理。MOD16A2是8天合成先把比例因子乘上再除以8得到日均蒸散var etDaily ee.ImageCollection(MODIS/006/MOD16A2) .filterBounds(roi) .filterDate(startDate, endDate) .map(function(img) { // ET原始值乘0.1得到mm/8day再除以8得到mm/day return img.select(ET) .multiply(0.1) .divide(8) .rename(ET_daily) .clip(roi) .set(system:time_start, img.get(system:time_start)); });月度ET总量在理想情况下等于日均ET乘以当月天数。由于8天合成窗口在月边界上是错开的严格做法需要做时间插值。实操里很多项目直接用当月ET影像的均值乘以天数或者按月份过滤后sum再乘系数。我的建议是不要在这里过度追求严谨先把量级跑对后面敏感性分析再检验ET误差的影响。搞一个看起来很精确但输入误差很大的模型没有实际意义。3.3 补给量异常的筛洗与导出将P、ET、R、ΔS都统一成月度影像后补给量计算就简化成了影像间的四则运算var runoffCoeff 0.1; // 区域经验值按研究区下垫面调整 var R monthlyP.multiply(runoffCoeff); // dSM 来自GLDAS区域平均后的逐月序列 var rechar monthlyP .subtract(monthlyET) .subtract(R) .subtract(dSM) .rename(Recharge);计算完成后第一件事不是出图而是先做两级检查。第一级观察逐月补给量时间序列如果大面积出现负值先排查是不是ET大于降水或者GLDAS土壤水变化的符号方向搞反了。第二级把研究区平均补给量与当地文献中的补给率做量级对比差出三倍以上就说明某个核心输入有问题不要硬着头皮往下做。导出可以用Export.table.toDrive把研究区平均逐月补给量导出成CSV方便在Excel里进一步分析和绘图。空间分布上把多年平均补给量导出为GeoTIFF再进GIS出图var meanRechar rechar.mean().clip(roi); Export.image.toDrive({ image: meanRechar, description: mean_recharge_2015_2020, region: roi, scale: 100, maxPixels: 1e13, crs: EPSG:4326 });如果真要逐月导出72张影像别在Client端写for循环直接在Export里把时间维度压缩进波段或者用Export.image.toDrive配合ImageCollection.toBands()否则会白白消耗配额后面细说。4. 我踩过的几个坑数据版本、尺度因子和时间对齐这一节是这次案例里真正花时间的地方。网上教程看着都顺自己一跑全是问题。以下三个坑是高频事故基本每次做类似计算都会撞上至少一个。4.1 MOD16A2的ET千万别忘了缩放系数第一次跑MOD16A2的时候我把ET直接当成毫米数用了结果算出来的补给量全年几乎全是负的——蒸散发普遍二三十毫米每八天降水一降下来根本扛不住。查了一圈才发现MOD16A2的ET需要乘0.1而且乘完之后是每8天毫米还要再除以8才是日均值。这个坑的隐蔽之处在于MOD16A2的LST波段不需要缩放如果你习惯只用其中一个波段再切到ET时特别容易忘记处理。我后来在每个.map()里都加上波段重命名和单位转换代码里明确写成ET_mm_per_day避免时间一长自己都忘了这个转换的存在。4.2 GLDAS土壤水单位与深度的误区GLDAS的土壤水变量有SoilMoi0_10cm、SoilMoi10_40cm、SoilMoi40_100cm单位是kg/m²。第一次拿到数据看到SoilMoi0_10cm这个变量名我误以为它是体积含水量或者10厘米土层的等效水深。实际上它的数值确实可以当mm用但那是整层10厘米的等效水深。更关键的是深度选择建议把0到100厘米三层加起来做ΔS这才大致覆盖根系层和包气带上部的蓄水范围。只取表层10厘米季节变化太剧烈会把降水信号全吸收掉补给量反而变得不可解释。这是一个物理含义问题不是代码问题但对结果的影响比任何语法错误都大。4.3 月度聚合与重投影的坑另一个高频问题是月度聚合。缺测日的存在让sum()的语义发生变化——如果一个月里只有25天有有效影像sum出来的值就不是全月总量而是部分天数的总量。稳妥做法是算有效天数用部分总量 ÷ 有效天数 × 月份天数推估全月或者在月度集合里对逐日影像求均值后再乘以天数。我在代码里选择的是均值乘天数因为均值会稀释单日极端值的影响更适合做区域趋势。两种方法都有误差但至少你要知道自己选择了哪一种误差而不是完全没意识到这个语义问题。重投影方面GEE在reduceRegions或导出时会按你指定的crs和scale重采样但默认采用的插值方式并不总是符合预期。CHIRPS是经纬度网格MOD16A2是MODIS正弦投影两者在流域边界上的像元覆盖会有半像元级别的错位。写代码时最好统一给Export指定crs: EPSG:4326并设置和输入产品匹配的scale避免导出一张看起来细腻、实际上是在多个粗分辨率产品之间插值出来的虚构细节的图。这事不细究看不出来但对空间分析和出图影响很大。顺带说一嘴GEE的配额问题。我这几年的体感是常规市级或中小流域分析免费配额每个月几千次请求基本够用。但如果ROI大、时间跨度长、导出任务多还是可能撞上限制。遇到这种情况别把所有年份一次性塞进一个任务把时间切成两到三年一段分批处理配额利用率会高很多。真到非大规模计算不可的份上再考虑配额升级或者本地重算也不迟。5. 结果靠什么验证监测井与文献互检遥感水量平衡算出来的补给量最怕的是自洽但不真实。一套代码跑通、图也画出来并不代表结果可信。我习惯用三个外部手段做交叉验证都不复杂但能把结论的底气撑起来。5.1 用区域地下水位波动粗校准如果研究区有地下水监测井这是一笔宝贵的外部数据。年尺度上如果区域地下水开采量相对稳定水位年际降幅与累积补给量的对比能反映整体水量均衡关系。月尺度上更直接雨季地下水位的抬升量和储水系数Sy的乘积理论上可以与当月补给量互相印证。公式非常简单ΔG ≈ Sy × ΔH储水系数Sy在砂质含水层一般0.05到0.25如果手头有抽水试验或者前人文献的数值用中值试算。这个验证方法不需要GEE参与但它能帮你在结论里写清楚流域月补给量约XX毫米与监测井水位抬升估算的YY毫米量级一致比单给一组模型输出可信得多。当然要记得扣除人为开采对水位的影响不然雨季水位没升多少你会误判补给量偏低。5.2 长时间均值与已发表研究对比还有一个省力有效的验证把多年平均补给量除以多年平均降水得到降水入渗补给系数。这个系数在不同气候区有很明显的经验范围湿润区0.15到0.30半干旱区0.05到0.15干旱区低于0.05。如果你的研究区是半干旱区算出来补给系数0.35基本可以断定要么ET被低估要么径流系数设得太小。把这个系数和相邻流域已发表的研究值做对比是最快的合理性检查。我这次案例算出来多年平均补给系数0.09勉强落在半干旱区经验范围内。然后去查了邻区两篇文献一篇是0.08一篇是0.11量级对得上。有了这一步至少敢把结果写进报告而不是只当练习。5.3 敏感性分析ET一变结果就变怎么应对我的习惯是固定降水不变把ET分别乘0.8和1.2跑两组敏感性测试看补给量的变化幅度。如果ET在±20%扰动下补给量直接从正值变成负值说明这个区域本身就是补给的边缘区结论就不要写死改成在降水正常年份补给量介于X到Y毫米之间结果对蒸散发算法较为敏感。在这次案例里ET乘0.8时补给量约为降水的11%ET乘1.2时补给量降至5%左右变化幅度可控。这个区间就是我给规划部门的最终口径。GRACE卫星数据也可以在趋势层面做对照但只适用于大区域中小流域别碰分辨率完全对不上。验证这事不用多高大上关键是让结论经得起追问。最后说点实操感受。整个流程从头跑一遍下来最浪费时间的不在GEE代码里而在数据解释上。我遇到过几轮结果异常的排查最后都回到同一个问题某个输入数据的物理含义被自己理解错了。所以第一次跑GEE水文分析的读者在做任何漂亮的专题图之前先把研究区的逐月降水、蒸散、径流、土壤水变化打印到一张Excel表格里逐个月份对着水量平衡方程手算一遍。这个朴素的检查能帮你省掉至少一周的返工时间。注册好账号之后先找个小区域把流程打通再慢慢放大范围是我能给的最实用的一条建议。