简介：全球陆面数据同化系统GLDAS是NASA推出的陆地水循环再分析数据集，融合卫星遥感、地面观测与陆面过程模型，为土壤湿度、雪水当量、冠层蓄水等关键变量提供连续时空覆盖。该数据免费公开，支持多分辨率与多模型输出，已成为区域水文研究、干旱监测及地下水替代分析的重要基础设施。在实际工程中，处理GLDAS往往涉及NetCDF格式解析、坐标方向检查、单位换算与区域裁剪等环节，尤其是土壤水分单位kg/m²与等效水深mm之间的关系，直接影响水储量计算结果。通过Python批量读取与栅格重采样，可高效提取流域尺度水储量变化序列，并进一步与GRACE重力卫星的陆地水储量异常进行交叉验证，为理解水文过程提供数据支撑。本文围绕GLDAS数据获取、预处理及水储量计算展开，完整梳理从原始数据到可分析时间序列的技术路线。 GLDAS这套数据，说白了就是NASA搞的一套全球陆地数据同化产品，把卫星遥感、地面观测和陆面过程模型揉在一起，输出一套连续、时空一致的陆地水循环变量。做水储量研究的人，尤其是搞区域尺度水文、干旱监测和地下水替代研究的，基本都绕不开它。几个月前我帮同事处理一批GLDAS.zip数据，踩了不少坑，也把从下载、解压、单位换算到计算水储量的整套流程重新捋了一遍，今天干脆把这些经验完整写下来。
这套数据好在哪？一是免费公开，分辨率覆盖1°和0.25°，时间能回溯到1948年；二是有多个陆面模型输出（NOAH、CLM、VIC、MOSAIC），你想对比模型不确定性也有素材。最关键的是，变量里有土壤湿度、雪水当量、冠层水储量，这些直接是水储量变化的组成部分，能用来和GRACE重力卫星的陆地水储量异常做交叉验证。如果你是刚接触GLDAS的研究生，或者是从GRACE数据处理转过来想补充陆地变量的人，这篇文章应该能帮你少走很多弯路。
1. GLDAS数据整体认知：为什么它成了水储量研究的标配
1.1 GLDAS到底是什么，能用来做什么
GLDAS的全称是Global Land Data Assimilation System，全球陆面数据同化系统。它做的事情可以用一句话概括：把所有能获取的陆地水循环观测信息都塞进陆面过程模型，让模型输出的状态变量尽可能接近真实。你可以把GLDAS理解成一个"重新分析"的数据集，类似气象里的ERA5，只不过它重点模拟的是土壤、积雪、蒸散发这些陆地过程，而不是大气环流。
它输出的变量很多，常见的包括地表和土壤各层体积含水量（Soil Moisture）、雪水当量（Snow Water Equivalent）、冠层蓄水量（Canopy Water Storage）、蒸散发（Evapotranspiration）、径流（Runoff）、地表和土壤温度（Temperature）等。这些变量说起来复杂，实际落地就是一个个NetCDF文件里面带单位的多维数组。
对做水储量的人来说，GLDAS最大的价值在于它能给出"陆地水储量"中除了地下水以外的几个关键分量。GRACE卫星测得的是陆地总水储量异常，包含土壤水、雪水、地表水和地下水变化。GLDAS虽然不含地下水，但它能提供土壤水和雪水的部分，所以在GRACE数据处理里，大家常用模型输出做尺度校正或者扣除某些分量。反过来，如果没有GRACE数据但只有GLDAS，你也能用GLDAS单独看某个区域土壤水和雪水储量的相对变化，这在干旱监测里非常实用。
1.2 版本与产品怎么选：GLDAS-2.1与GLDAS-2.2的取舍
GLDAS现在主流的分两代。第一代GLDAS-1有1°分辨率，时间范围从1979年到当前近实时，产品包括NOAH、CLM、VIC和MOSAIC四个模型。第二代GLDAS-2又分了两套：GLDAS-2.0是纯水文模拟，从1948年到2014年，没有同化气象观测；GLDAS-2.1从2000年到现在，同化了大量观测数据，表现更接近真实，也是目前用得最多的。
我个人的经验是，做现代水文分析和水储量变化研究，首选GLDAS-2.1的NOAH模型0.25°产品。原因很简单：2.1版本同化信息更丰富，0.25°空间细节比1°强很多，尤其在流域尺度上能看出空间差异。如果研究历史长序列、要看几十年的变化趋势，那GLDAS-2.0的1948年起点就很有价值，但由于它没同化观测，单年数值精度不如2.1。
在看数据的时候需要注意，GLDAS-2.1里面0.25°和1°两个分辨率的产品都存在。经常有人下载的时候没注意，混用了两个分辨率还浑然不觉。0.25°文件在变量名和维度上通常更细，存储量也大不少。我的建议是，如果电脑配置允许而且研究的不是全球尺度，就直接用0.25°，后面做流域平均时空间代表性更好。
1.3 空间分辨率与时间尺度的选择思路
GLDAS数据的时间分辨率有3小时、6小时、日、月等几种。做水储量变化研究，月尺度是最常用的，因为水储量本身季节性强，月均值可以滤掉不少日内高频噪声。3小时产品一般用来驱动模型或分析极端降水事件，日常水储量分析用不上。
空间上做全球尺度的分析选高分辨率没太大必要，1°全球范围能省不少内存和处理时间。但如果做某个流域，比如长江上游、华北平原、美国中部大平原，0.25°和1°的差别就很明显——0.25°能看到更细致的土壤水空间分布，1°则会把河谷和山坡平均成一个值，细节丢失严重。
我自己的经验是先用0.25°数据做分析，如果后续要跟粗分辨率GRACE做对比，再把GLDAS重采样到GRACE的1°或0.5°网格。重采样的过程其实也是一种空间滤波，能减少尺度不匹配带来的误差。这个思路在下面实操部分还会细讲。
2. GLDAS数据格式与单位深度解析：别让单位坑了你
2.1 NetCDF文件结构与变量命名规则
GLDAS数据在21世纪以后基本都是NetCDF格式，少数旧版本是二进制或GRIB。NetCDF格式的好处是自描述性强，你打开文件就能看到每个变量的名字、单位、维度、坐标轴信息，不需要额外查说明书。
每个GLDAS的NetCDF文件通常包含lat、lon、time（或是time向量）、然后就是一堆2D或3D变量。维度顺序一般是time在最前，然后是lat和lon。注意这里有个非常容易踩的坑：GLDAS的纬度方向是从北到南还是从南到北？不同时期产品并不完全一致，有的lon从-180到180，有的从0到360。你第一次画图如果发现图像上下颠倒或者左右错位，别急着改图，先去查一下文件的坐标范围。
变量命名看着复杂，其实比较规律。土壤湿度通常叫SoilMoi0_10cm_inst、SoilMoi10_40cm_inst、SoilMoi40_100cm_inst、SoilMoi100_200cm_inst。如果你用的是NOAH模型，这些名字基本就是固定写法。雪水当量叫SWE_inst，冠层水叫CanopInt_inst。后面带不带_inst取决于数据是瞬时值还是平均通量，"inst"是instantaneous，表示该时刻的瞬时状态。
2.2 核心变量单位与物理含义
GLDAS变量的单位在NetCDF文件里通常有标准属性，但中文资料里经常把它搞混，这里一定要搞清楚。
土壤湿度SoilMoi的各层单位是kg/m²，这个单位在数值上等于等效水深毫米（mm）。为什么？因为水的密度是1000 kg/m³，如果1平方米面积上、某一层土壤含水相当于水柱高度 h 毫米，那么水分质量为 0.001 h m × 1 m² × 1000 kg/m³ = h kg，正好1 kg/m²对应1 mm。所以GLDAS输出的土壤湿度从单位看是kg/m²，你完全可以当成"等效水深mm"来读数。
雪水当量SWE_inst单位也是kg/m²，同样等效于mm水深。冠层水CanopInt_inst单位是kg/m²，也按同样方式理解。这个单位统一性对后面计算水储量太重要了——你想叠加土壤水、雪水和冠层水，直接做加法就行，不需要额外换算。
温度变量比如土壤温度SoilTemp_inst单位是K，蒸散发变量单位则是W/m²或kg/(m²·s)，这都是通量类的，算累计量时需要乘以时间步长。我们做水储量变化只用状态变量，不涉及蒸散发单位转换，但如果你之后想算水平衡闭合度，就要仔细处理这些通量单位。
2.3 单位换算实操：从kg/m²到mm水深的秘密
前面说了kg/m²和mm水深数值相同，但很多同学看到NetCDF里soil moisture的单位写着kg/m²还是忍不住怀疑，担心是不是要除以土壤密度或者乘以某层厚度。这里我直接从定义再推一遍，保证你踏实：
层内含水量对应等效水深 h（mm）。某层土壤含水量体积分数为 θ（无量纲），层厚度为 D（mm），那等效水深就是 θ × D。GLDAS输出的SoilMoi就是 θ × D 的积分结果，单位按kg/m²标。因为水的密度是1000 kg/m³，把等效水深mm换算到质量除以面积时，1 mm水深 = 1 kg/m²。所以两者刚好数值相等。
实际操作里我从来不做任何单位除法，直接读出来的SoilMoi值就按mm算，用累加来做土壤剖面总含水量。比如NOAH模型常见的四层土壤，深度范围是0-10、10-40、40-100、100-200 cm，四层含水量相加，就得到0-200 cm土壤剖面总含水量，单位是mm。这个数据再和SWE相加，就是考虑了土壤和积雪的地表水储量。
有一个小提醒：有些旧版数据或转换工具可能输出的是体积含水量（m³/m³），这种就没法和雪水当量直接相加了。判断方法很简单，看值的范围：如果是体积含水量，数值一般0到0.6之间；如果是kg/m²或mm，常见数值在10到500之间。看到数值范围不对，要先追查Data变量单位，别硬算。
2.4 数据缺失值与有效范围的识别
GLDAS的NetCDF里一般没有NaN，但个别区域、个别时间可能出现填充值。填充值通常是一个特别大的负数，比如-9999，或者明显超出物理范围的数值。你不处理它，后面求区域平均或画图时会得到一片诡异的黑色或极大值。
我的习惯是打开文件后先看一眼变量的min和max。如果min是-9999或者比-100还小的数值，说明有填充值，需要先用掩膜处理。netCDF4库里读取时可以用variable[:].data，把填充值替换成numpy.nan，然后再算统计量。判断完填充值后还要看有效范围，比如土壤水分不可能为负，如果出现负值大概率是数据噪声，区域统计时建议把它mask掉。
3. GLDAS数据获取与预处理实操：从下载到能用数据
3.1 下载前的准备：账号、数据订购与批量下载
GLDAS数据在NASA的Earthdata平台下载，地址是earthdata.nasa.gov。首次使用必须注册账号，整个过程跟着邮箱确认走就行，不需要审批。注册登录后会跳到GES DISC，这是戈达德地球科学数据与信息服务中心，GLDAS数据主要在这托管。
下载路径上有两种选择：网页手动点选，或者用脚本批量下载。手动下载适合临时拿几个时间片试试，做研究要一整段时间序列，建议直接用脚本。GES DISC支持HTTPS链接，获得下载链接列表后用wget或curl加认证信息就行。我用的是GES DISC提供的批量下载txt文件，里面每一行是一个文件的完整URL。然后写一个for循环，逐行用wget下载，记住加--user和--password参数。
一个重要的经验：加个限速或断点续传参数。NASA服务器在高峰期经常断流，wget的-c参数能断点续传，--tries=5能自动重试。我下载两年逐月数据大概40多个文件，高峰期遇到过3次中断，有了-c参数一次跑完，没再回头补数据。
3.2 解压GLDAS.zip后：文件命名与时间序列的匹配
下载下来的文件通常长这个样子：
GLDAS_NOAH025_M.A200001.021.CSL07.20220301.tar.gz GLDAS_NOAH025_M.A200002.021.CSL07.20220301.tar.gz其中A200001代表2000年1月，"021"代表GLDAS-2.1版本，后面是模型和发布时间戳。这些是月度文件，每个文件内部解压出来一般是一个NetCDF。如果你下载的是3小时产品，文件名里会出现A2000010100之类的字样，也就是2000年1月1日00时。
解压命令很简单：
tar -xzf GLDAS_NOAH025_M.A200001.021.CSL07.20220301.tar.gz解压完记得检查文件大小，有些NetCDF压缩后只有几MB，但解压出来有几十MB。把这个过程重复40次以后，你会意识到批量脚本的必要性——手动解压40个文件，虽然不累，但浪费时间，还容易漏。
我建议把解压和重命名写成一个bash循环，逐个解压后按照年份月份规律放到一个统一目录。这样后面Python读取时会非常方便，只需要用glob.glob挨个处理。文件名里其实自带时间信息，尽量不要自己重命名成a、b、c这类无意义名称，否则处理几十个文件的时候，对应关系会乱成一团。
3.3 区域裁剪与时间聚合：一整套Python处理逻辑
处理GLDAS最频繁的需求是"提取某个经纬度范围、某段时间序列的区域平均值"。以下是我常用的一个示例脚本，完整展示了打开文件、裁剪、取区域平均、叠加变量、输出CSV的流程：
import netCDF4 as nc import numpy as np import pandas as pd # 设定研究区域：比如中国华北平原附近 lon_min, lon_max = 110, 120 lat_min, lat_max = 32, 42 file_path = "GLDAS_NOAH025_M.A200001.021.CSL07.20220301.nc4" ds = nc.Dataset(file_path) lons = ds.variables['lon'][:] lats = ds.variables['lat'][:] # 找最近格点索引 lon_idx = np.where((lons >= lon_min) & (lons <= lon_max))[0] lat_idx = np.where((lats >= lat_min) & (lats <= lat_max))[0] # 读取各层土壤水和雪水当量 soil_10 = ds.variables['SoilMoi0_10cm_inst'][:, lat_idx, lon_idx] soil_40 = ds.variables['SoilMoi10_40cm_inst'][:, lat_idx, lon_idx] soil_100 = ds.variables['SoilMoi40_100cm_inst'][:, lat_idx, lon_idx] soil_200 = ds.variables['SoilMoi100_200cm_inst'][:, lat_idx, lon_idx] swe = ds.variables['SWE_inst'][:, lat_idx, lon_idx] canop = ds.variables['CanopInt_inst'][:, lat_idx, lon_idx] # 区域平均，注意先求时间轴上的平均值再求空间平均 soil_total = (soil_10 + soil_40 + soil_100 + soil_200 + swe + canop).mean(axis=(1, 2)) print(f"{ds.variables['time'].units} -> 时间长度 {len(soil_total)}") ds.close()这个脚本读出来的是该文件内所有时间步的剖面水储量。如果是月数据，每个文件通常只有一个时间点，区域平均后就得到一个数值。把所有月份文件输出拼起来，就是时间序列。如果文件是3小时的，你可以按日、按月再聚合，用pandas的resample就行。
3.4 重采样与插值的取舍建议
GLDAS自带的空间分辨率是0.25°或1°，但有时候做区域分析需要统一到5km、10km网格，或者要和GRACE的1°网格对齐。这时就要做重采样，方式有两种：栅格数据重投影插值，或者网格面积权重平均。
我自己更推荐对气象水文数据用面积权重平均，而不是双线性插值。原因是双线性插值在数值场里会制造人为的平滑和振荡，而面积权重平均能保留原始网格的物理量守恒性质。比如把0.25°的GLDAS处理成0.5°，用xESMF库是最省事的：
import xarray as xr import xesmf as xe ds = xr.open_dataset("GLDAS_xxx.nc4") ds_out = xr.Dataset({ "lat": (["lat"], np.arange(32, 42, 0.5)), "lon": (["lon"], np.arange(110, 120, 0.5)), }) regridder = xe.Regridder(ds, ds_out, "bilinear") ds_regrid = regridder(ds)xESMF底层调用ESMF，支持很多插值方案，其中"conservative"方案能严格守恒，适合通量变量。但对土壤湿度这类状态变量，bilinear和conservative差别不大，用bilinear就够快够稳。需要注意插值后会生成一些非物理的边界噪声，建议在重采样后对覆盖范围再检查一遍max和min，防止出现离谱的极值。
4. 利用GLDAS计算区域水储量变化：一条完整的技术路线
4.1 水储量变化的基本公式与变量叠加
从GLDAS提取水储量，本质就是把陆地水储量中模型能描述的部分都加起来。GLDAS语境下，可写的总水储量（Terrestrial Water Storage, TWS）通常定义为：
TWS = 冠层水 + 雪水当量 + 各层土壤水
注意这个公式不含地下水，也不含河流湖泊的地表水。GLDAS模型有的版本包含简单的河道汇流蓄水量，但标准输出里没有单独给出，你算的时候也用不到。
单位上全部是等效水深mm，因此直接加和即可。如果要看水储量变化，而不是绝对储量，就做去均值化处理，即从每个格点的月值减去多年月平均或全序列平均。这样得到的就是Water Storage Anomaly，可以和GRACE的结果直接对比。
做华北平原这种地下水超采区时会发现，GLDAS得到的TWS趋势和GRACE相比明显偏弱，因为地下水变化在GLDAS里体现不出来。这个不是数据问题，而是物理过程缺失，写文章的时候一定要注明。
4.2 从GLDAS提取逐月水储量变化序列：完整示例
把前面裁剪逻辑放到循环里，逐月读取所有文件，输出一个区域平均水储量时间序列，用pandas组织成DataFrame：
import glob import pandas as pd import numpy as np file_list = sorted(glob.glob("GLDAS_NOAH025_M.A*.nc4")) ts = [] for fp in file_list: with nc.Dataset(fp) as ds: # 提取年月 date_str = fp.split(".")[1] # A200001 -> 200001 date = pd.to_datetime(date_str[1:], format="%Y%m") # 读取所有分量 soil_10 = ds.variables["SoilMoi0_10cm_inst"][0, lat_idx, lon_idx] soil_40 = ds.variables["SoilMoi10_40cm_inst"][0, lat_idx, lon_idx] soil_100 = ds.variables["SoilMoi40_100cm_inst"][0, lat_idx, lon_idx] soil_200 = ds.variables["SoilMoi100_200cm_inst"][0, lat_idx, lon_idx] swe = ds.variables["SWE_inst"][0, lat_idx, lon_idx] canop = ds.variables["CanopInt_inst"][0, lat_idx, lon_idx] tws = (soil_10 + soil_40 + soil_100 + soil_200 + swe + canop).mean() ts.append({"time": date, "tws_mm": tws}) df = pd.DataFrame(ts).set_index("time") df["tws_anom"] = df["tws_mm"] - df["tws_mm"].mean() df.to_csv("tws_anom.csv")整个流程下来数据结构非常干净，后面画图、做趋势分析都很顺手。需要注意文件读取时用with自动关闭，避免打开过多句柄导致内存泄漏。如果你处理的是几十年的月数据，一次读全部文件也没问题，NetCDF文件本身还在磁盘，Python只是读取需要的切片。
4.3 与GRACE/GRACE-FO数据的对比验证
得到GLDAS的TWS anomaly后，最典型的操作就是和GRACE/GRACE-FO的数据对比。GRACE的分辨率很粗，而且输出的是全球0.5°或1°网格的陆地水储量异常。对比前至少要做两步：
第一步，统一空间范围和时间分辨率。GRACE月度数据最好和GLDAS在同一个月平均，并且都重采样到1°网格。第二步，去掉长期趋势和季节性周期的影响。GRACE数据里已经做了球谐系数解算，自带低阶信号，GLDAS的异常值在对比时也会去均值化，两边都统一口径。
实际操作中，很多论文直接用相关系数和RMSD来衡量两个数据之间的一致性。在华北、亚马逊、撒哈拉这些地下水变化显著的地区，两者差距会比较大。这时候可以看一下GLDAS缺了哪个分量——比如地下水变化如果是主项，GRACE会明显大于GLDAS；如果是土壤湿度主导（比如半干旱草原），两个数据的一致性就会很高。这种对比本身就是很好的科学分析。
4.4 结果可视化与信号去噪
画GLDAS水储量时间序列时，建议直接画月值加一个12个月的滑动平均，既能看季节循环，也能看出年际变化。用matplotlib画两条线就行，原始月值用浅色细线，滑动平均用深色粗线。这样写论文很直观，审稿人一眼能看到趋势。
空间图上要选好colormap，水储量变化类的变量建议使用BlueWhiteRed类的发散色标，中间值设为0。要注意GLDAS数据在极地和高寒区域可能出现较大的雪水当量极端值，如果不将色标范围设置合理，很多区域会变成一片红色或蓝色，信息全丢。我一般用np.percentile把色标上下限设为数据的5%和95%分位数，能自动适应大多数区域。
5. 常见问题与排查技巧实录
5.1 GLDAS读取常见报错与修复方法
用netCDF4读GLDAS最常见的报错是File format not recognized或Variable not found。前者大多数是因为文件没解压完，或者下载过程中出错了只剩个不完整镜像。我建议先看文件大小，如果明显小于正常值，直接重下，不要尝试修复。
Variable not found则多半是读错了变量名。不同GLDAS版本变量名有差异，GLDAS-2.0和2.1的土壤分层命名几乎一致，但如果你下载的是VIC模型输出，变量名会变成SoilMoist_tavg之类，并不完全相同。读取前先print一下ds.variables.keys()，看到实际名称再写代码，这能省下很多排查时间。
5.2 单位换算中的经典陷阱：kg/m²与mm到底有没有区别
这个陷阱我在2.3节已经详细讲过了，但实际操作中还会遇到一个衍生版本：有人把kg/m²直接乘以1000当成mm，结果数值大了1000倍，画出的图极不自然。这里再强调一次，GLDAS里的kg/m²就是等效水深mm的数值，不要做额外变换。
如果真要换算，唯一的场景是当你需要把土壤水体积分数（m³/m³）转成mm时，才需要用体积分数乘以该层厚度（mm）。GLDAS输出已经帮你做了这个积分，你直接拿去叠加就行。这个"已经帮你做了"是关键，很多人总想自己再算一遍，反而算错。
5.3 坐标系与网格方向怎么检查
GLDAS数据的经纬度坐标并非所有版本都一致，有的版本lat从-59.875到89.875递增，有的则是从90到-60递减。如果你读文件时不检查坐标方向，直接按索引取区域，会选到完全相反的区域。
我的检查方法是，在读取后立刻输出lons[0], lons[-1], lats[0], lats[-1]。如果是lats从大到小，脚本里要么用np.flip翻转lat轴，要么用切片的时候调整索引位置。还有一点，lon范围有的是-180到180，有的是0到360，如果你研究区域跨越180°经线（比如太平洋），就需要做换经度范围的处理。用xarray时直接ds.assign_coords(lon=(ds.lon % 360))可以快速把-180~180转成0~360，但要注意转完之后可能需要按lon排序。
5.4 数据缺失与时间戳错位处理
GLDAS月度产品偶尔会有两三个月数据缺失的情况，尤其在一些新版本发布前的数据空洞。处理方式看你的分析目的：做趋势分析建议直接忽略缺失月份，不要插值填充，因为插值会人为拉平趋势；做季节分析则可以按多年平均填充，把缺失值替换成气候态值。
时间戳错位是一个更隐蔽的问题。有些GLDAS文件的时间变量是自1850年1月1日以来的天数，有些则是模拟开始时间起的小时数。读取时我习惯打印ds.variables['time'].units，然后根据units做转换。用netCDF4自带的num2date很方便，能生成标准的datetime对象。如果你用pandas解析，直接pd.to_datetime(file_str里的时间信息)反而更麻烦，不如用读出来的真实时间。
5.5 下载速度慢与断流问题
NASA服务器在国内访问速度时好时坏，这个局面短期没啥好办法。我的经验是把下载拆成小批次，一次不要超过20个文件，同时在命令行里加上--wait和--random-wait来降低服务器压力。如果中途断掉，用wget的-c续传即可。
批量下载的脚本建议保存为bash脚本，设好用户名密码变量，别把密码写死在博客或仓库里，避免漏出。下载完成后立刻检查文件个数是否和预期一致，缺哪个补哪个，别等处理数据时才发现前一年数据少了一个月。
6. 进一步扩展：从GLDAS延伸到更完整的水储量分析
6.1 GLDAS与SMAP、SMOS的融合思路
GLDAS虽然好，但土壤水文模拟仍有不确定性。如果你想提高浅层土壤湿度的精度，可以把SMAP或SMOS的卫星土壤湿度产品同化进来，或者做简单的数据融合。最粗糙但有效的办法是：用SMAP观测的0-5cm土壤湿度来校正GLDAS的0-10cm土壤湿度，按比例调整深层。这有点tricky，但对干旱监测来说效果立竿见影。
我个人在区域干旱监测项目里试过，GLDAS+SMAP融合后，表层土壤湿度的空间分布和站点观测的一致性明显提升。如果要做近实时监测，这套方案比单纯用GLDAS更可靠，因为GLDAS模型的近实时输出有一定滞后和误差积累。
6.2 结合GRACE做地下水储量变化的估算思路
既然GLDAS缺少地下水，而GRACE能看到总水储量，那自然可以做个减法：地下水储量变化约等于GRACE总水储量异常减去GLDAS模型中陆地水储量异常。这个思路被很多论文称作"GRACE minus GLDAS"方法。
实际做的时候有几个坑：GRACE的分辨率很粗，直接减去0.25°的GLDAS会引入尺度和泄漏误差，建议先把GLDAS重采样到GRACE网格再减。减去后结果中还存在大量噪声，需要用高斯滤波或者去相关滤波处理。至少我在华北平原做出来的地下水枯竭信号，数值上和官方公报中的埋深变化趋势能对上，说明方法在区域尺度上是可用的。
6.3 未来GLDAS数据更新与替代方案
GLDAS不会一直是唯一选择。现在欧洲的ERA5-Land、国内的CLDAS，都提供类似的陆地水循环变量，分辨率和质量各不相同。如果你做全球尺度，ERA5-Land的0.1°分辨率比GLDAS的0.25°更有优势；如果做中国区域，CLDAS在站点同化上有自己独特的优势。
但GLDAS仍然有它的不可替代性：产品稳定、版本序列长、模型多样、变量体系完善。在论文里做基准数据或交叉验证时，GLDAS依然是最稳妥的起点。如果你把多套数据的对比结果写进论文，审稿人的接受度也会更高。
个人实操感受
踩了这么多坑之后，我对GLDAS的判断很简单：它不完美，但你一定要把它用到顺手。单位、坐标、命名规则这些坑，踩过一次之后记得写在笔记里，下次就是闭眼操作。数据下载和预处理这一步最枯燥，也最值得花时间写脚本自动化，一劳永逸。
最后分享一个小技巧：处理完GLDAS后，一定要顺手把时间序列画出来，先求一个区域平均值的曲线，看看峰值和谷值出现的月份是否符合流域气候特征。我见过太多次因为坐标范围写反导致整条曲线反相、却还继续往下算的例子。如果曲线形状不对，先别急着写结论，回头检查数据读取和坐标选择，往往能省下几天返工时间。
本文还有配套的精品资源，点击获取
