企业官网建设流程全解析

前阵子接了个活儿，某半干旱地区要做地下水资源可持续利用规划，第一个问题就是：这块地的地下水补给量到底有多少？水文站网稀疏，监测井资料零散，地面实测这条路基本走不通。我给的方案很直接——用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. 数据选型：降水、蒸散、土壤水、地形这四件套怎么配

多源遥感数据的"多"不是越多越好，而是每一类变量选一个最合适的产品，再留一个备选做交叉验证。我这次搭配下来比较顺手的一套，放在下面这个表格里。

2.1 降水数据：CHIRPS打底、GPM补充

降水是补给量计算里最刚性的输入，它的精度直接决定最终结果的可信度。我习惯用CHIRPS V2.0作为主数据源，原因很实际：序列从1981年延续到现在，空间分辨率0.05°（约5公里），时间分辨率逐日，单位就是mm/day，GEE里一个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表格里，逐个月份对着水量平衡方程手算一遍。这个朴素的检查，能帮你省掉至少一周的返工时间。注册好账号之后，先找个小区域把流程打通，再慢慢放大范围，是我能给的最实用的一条建议。