简介这是一份面向GRACE/GLDAS水文数据处理学习者的Matlab程序包旨在解决GLDAS数据从读取、网格化到球谐展开的处理流程便于与GRACE反演结果结合计算地下水储量变化。包内依据水量平衡方程组织陆地水储量、地表水储量、冰后回弹改正及地下水储量等计算环节其中地表水储量部分依托GLDAS水文模型数据实现程序可直接运行并输出结果文件需要成图时可用GMT绘制脚本函数划分清晰便于二次修改。压缩包共57个文件包含48个GLDAS NOAH月值nc4数据文件、7个Matlab脚本有主程序、去相关滤波、高斯滤波、球谐展开等函数以及2个txt说明文档整体约84.79MB自带测试数据便于对照验证。目前已有1031人学习下载适合正在开展GRACE与GLDAS结合研究的科研人员、研究生及工程师通过模块化脚本可快速掌握从原始nc数据处理、滤波平滑到球谐展开的完整链路并可与GRACE反演结果对比估算地下水储量变化。 记得我第一次拿到GRACE反演的陆地水储量变化结果时盯着那张全球网格图想了一个下午这颗卫星测的是重力场变化反演出来的是整个陆水柱的总储水量变化可它根本没法告诉你这些水储量的增加或减少究竟是发生在地下水层还是只是浅层土壤湿度变了。如果我想单独拿出地下水储量变化就必须借助一套外部数据把非地下水的部分扣掉而GLDAS就是那个最常用的数据源。所谓结合GRACE和GLDAS计算地下水储量变化核心思路就是残差法从GRACE得到的陆地水储量变化中减去GLDAS模拟的土壤水、雪水、冠层水等地表储水分量剩下的残差就解释为地下水储量变化。这套流程里GLDAS数据的读取、单位转换、空间重采样和时间对齐每一步都有坑。这篇文章就按我实际跑通这条链路的顺序把Matlab处理GLDAS数据的完整思路和可复用代码记录下来。适合正在做GRACE水文学应用、水资源变化研究的研究生和科研人员参考也适合刚接触遥感水文数据处理的初学者。1. GLDAS在残差法里的位置为什么算地下水要先搞清楚它1.1 一个公式把GRACE和GLDAS串起来GRACE反演得到的陆地水储量变化TWSC本质上是某一区域内所有水储量的总和随时间的变化。简化表达就是TWSC 土壤水变化 雪水当量变化 冠层截留水变化 地表水变化 地下水变化这个公式里GRACE只给出等式左边的总量。右边的前几项正好可以由GLDAS提供。GLDAS是同化卫星和地面观测资料驱动的陆面过程模型输出NOAH模型等可以提供多层土壤湿度、雪水当量SWE、冠层截留水CanopInt、径流、蒸散发等变量。用GLDAS作为辅助数据把GRACE总量中属于浅层土壤和地表的水量扣除剩下那个残差就认为是地下水储量变化。写成Matlab里最直观的计算式就是GWS GRACE_TWSC - (SoilMoi_total SWE CanopInt RiverStorage)之所以能这么写是因为GLDAS的陆面模型本身不显式模拟深层地下水土壤水分层通常到2米或更浅雪水和冠层水都是地表过程。它恰恰缺的那部分就是你关心的地下水。1.2 GLDAS有哪些分量能扣哪些不能扣使用残差法之前必须先确认GLDAS输出变量里到底打包了哪些水分量。不同陆面模型的输出变量命名不一样但物理含义基本一致。以最常用的NOAH v2.1产品为例SoilMoi0_10cm_inst、SoilMoi10_40cm_inst、SoilMoi40_100cm_inst、SoilMoi100_200cm_inst四层土壤湿度单位kg/m²合起来是0到200厘米土柱的水量。SWE_inst雪水当量单位kg/m²代表积雪融化后对应的水柱高度。CanopInt_inst冠层截留水单位kg/m²植物叶片表面附着的水分。Qs_acc、Qsb_acc地表径流和地下径流累计量一般用来做水量平衡闭合残差法里不一定要扣。Evap_tavg蒸散发属于通量不是储量不能直接拿来扣。这里有个容易混淆的点有的版本GLDAS产品会直接给出总的雪水有的会给出不同积雪深度。另外不同模型对径流的处理也不一样NOAH把地表和地下径流分开但如果你用CLSM模拟水分变量就可能变成另外一套命名。所以拿到数据第一件事一定是看变量列表和单位而不是直接套公式。还需要记住的是GLDAS不含水库蓄水、湖泊蓄水变化、冰川变化和深层地下水。这些分量如果存在会全部混进残差法的地下水结果里。也就是说你用这个方法得到的GWS其实是人类活动和深层储水量的综合信号尤其在地表水库密集或存在明显冰川变化的区域解读结果时必须加个小心。2. 动手前先定好三件事数据产品、网格、时间范围2.1 GLDAS产品怎么选GLDAS有多个版本和陆面模型。做地下水储量变化研究我建议优先用NOAH模型因为它历史最长、变量最完整、使用的人最多出问题好查文献。NOAH产品常见版本包括GLDAS-2.1 NOAH空间分辨率0.25度时间范围2000年至今3小时或月平均产品都有。GLDAS-2.0 NOAH空间分辨率1度覆盖1948到2014年适合做长时间序列趋势分析。GLDAS-2.1 CLSM、VIC同样0.25度但土壤分层和输出变量不同。具体选哪个主要看你要研究和GRACE哪个时间段匹配。GRACE有效数据期是2002到2017年GRACE后续有GRACE-FO延续所以GLDAS-2.1 NOAH 0.25度是大多数人的选择。如果GRACE用的球谐系数产品经过平滑处理后有效分辨率在3度左右那GLDAS用1度还是0.25度对最终结果的影响其实不大但推荐先用0.25度做处理再重采样到GRACE网格。产品分辨率时间跨度适用场景GLDAS-2.1 NOAH0.25°2000年至今与GRACE主任务期匹配最常用GLDAS-2.0 NOAH1.0°1948-2014长趋势重建覆盖更早GLDAS-2.1 CLSM0.25°2000年至今对地下水文过程模拟更好但变量说明更复杂2.2 GRACE产品和GLDAS的分辨率怎么对齐GRACE数据产品主要分成两类球谐系数SH产品和mascon产品。SH产品需要自己截断、去相关滤波、高斯平滑处理流程长结果对平滑半径很敏感。mascon产品是官方处理好的网格产品例如JPL mascon和CSR mascon已经做了泄漏校正和尺度因子处理空间分辨率约0.5度或1度。如果你目标是快速算地下水储量变化我个人强烈建议直接用JPL mascon省去大量滤波调参的功夫把精力放在GLDAS处理和残差法本身。分辨率对齐上两种常见做法把GLDAS 0.25度数据通过面积加权或双线性插值重采样到GRACE mascon的0.5度或1度网格上。直接下载GLDAS 1度产品天然和很多SH产品或1度mascon网格一致。更严谨的做法是先对GLDAS做空间粗化再和GRACE相减而不是先相减再插值。因为GRACE信号本身代表大尺度平均如果先在高分辨率上算GWS再粗化容易把局部土壤水信号当作区域地下水信号带入结果。先把GLDAS粗化到和GRACE同一个空间支撑物理学上更自洽。2.3 研究区掩膜和时间范围的准备处理这类数据之前我习惯先把研究区边界经纬度范围、陆海掩膜和GRACE数据覆盖时段列出来。掩膜这点很关键。GRACE mascon在海洋区域会输出零值或负值如果你直接拿去和GLDAS做差海岸线附近的像元会瞬间出现一个很大的伪异常。最好准备一份标准land mask对所有包含海洋像元的结果做掩膜处理。掩膜可以用GLDAS提取、用GRACE的系数文件自带掩膜也可以用常见的陆海掩膜数据。时间范围上GRACE和GLDAS不需要完全一致的起止时间但重叠期内不能有太大缺测否则趋势拟合时会出问题。3. Matlab批量读GLDAS框架与可复用代码3.1 先学会看数据再写读取代码很多同学的第一个错误是拿到NetCDF文件就直接ncread也不管变量叫什么、单位是什么、时间坐标怎么编码。GLDAS文件变量命名在不同版本之间差异很大必须先看元数据。% 查看单个GLDAS文件的基本信息 info ncinfo(GLDAS_NOAH025_3H.A20020101.0000.021.nc4); % 列出全部变量名 for i 1:length(info.Variables) fprintf(%s: %s\n, info.Variables(i).Name, info.Variables(i).Attributes); end运行之后重点确认几件事维度顺序是不是(lon, lat, time)变量单位是不是kg/m²时间坐标的参考日期是什么很多NetCDF的time是自1900-01-01起的天数或小时数。这些信息看起来琐碎但单位弄错一个后面的GWS就等于白算。3.2 批量读取的时间循环与结构化存储GLDAS数据按小时或3小时输出一个月有几十甚至上百个文件。最稳妥的做法是逐文件读取在内部完成空间范围裁剪后再聚合到月尺度。我习惯用Matlab的dir函数先拿到所有文件列表再按文件名解析日期最后用结构体或三维数组存储结果。rootDir D:/gldas_data/; files dir(fullfile(rootDir, GLDAS_NOAH025_3H.A*.nc4)); % 按文件名提取时间 for i 1:length(files) tokens regexp(files(i).name, A(\d{4})(\d{2})(\d{2}), tokens); if ~isempty(tokens) yyyy str2double(tokens{1}{1}); mm str2double(tokens{1}{2}); dd str2double(tokens{1}{3}); files(i).dateNum datenum(yyyy, mm, dd); end end读取变量时只需要读取研究区范围内的数据即可。比如你研究区经度70到100度纬度25到45度那就用ncread的start和count参数做裁剪避免把全球数据都读进内存尤其是0.25度数据的时间序列累计下来内存压力很大。lonRng [70, 100]; latRng [25, 45]; lon ncread(fname, lon); lat ncread(fname, lat); ilon find(lon lonRng(1) lon lonRng(2)); ilat find(lat latRng(1) lat latRng(2)); sm0 ncread(fname, SoilMoi0_10cm_inst, [ilon(1), ilat(1), 1], [length(ilon), length(ilat), 1]);3.3 从瞬时值到月平均的处理逻辑GLDAS的原始产品很多是3小时分辨率要做月平均最简单的逻辑是每天该变量所有时次取平均得到日值再对一个月内所有日值求平均。但实际处理时如果文件量很大我建议直接用GLDAS官方提供的月平均产品文件少很多精度也能满足大多数研究。如果必须自己聚合注意区分变量类型带_inst后缀的是瞬时值带_tavg后缀的是时间段平均通量这两类不能混在一起取平均。% 聚合土壤水瞬时值到月平均 % 假设filesInMonth存放该月所有文件lon和lat已裁剪 monthData zeros(length(ilon), length(ilat), length(filesInMonth)); for k 1:length(filesInMonth) fname fullfile(rootDir, filesInMonth(k).name); monthData(:,:,k) ncread(fname, SoilMoi0_10cm_inst, ... [ilon(1), ilat(1), 1], [length(ilon), length(ilat), 1]); end sm0_month mean(monthData, 3);这一步的输出建议保存成Matlab的.mat文件或者标准NetCDF格式方便后续多个实验复用。不要每次都从原始文件重新读时间成本太高。4. 单位对齐和空间对齐残差法最容易翻车的地方4.1 kg/m²与mm等效水高别小看scale_factorGLDAS的土壤湿度和雪水当量单位是kg/m²。因为水的密度是1000 kg/m³1 kg/m²的等效水高恰好就是1 mm。也就是说数值上kg/m²可以直接读成mm等效水高不需要额外转换。而GRACE mascon产品的单位通常是cm等效水高有些是mm。这就出现了一个非常常见的低级错误GLDAS是mmGRACE是cm两者直接相减地下水储量变化结果差了10倍。我处理JPL mascon时习惯先统一成mmGRACE_mm GRACE_cm * 10;如果是球谐系数反演出来的TWSC单位取决于你自己的换算方式通常也是cm等效水高或mm等效水高务必先确认清楚再相减。另一个坑是部分NetCDF变量带scale_factor和add_offset属性。GLDAS官方产品一般已经换算好但有些衍生产品或重发布数据会有缩放因子。读取后建议立刻检查数据范围是否合理土壤水0到几百kg/m²雪水当量在干燥区接近0在高山区可以到几百。如果数据范围不合理大概率是scale_factor没处理。4.2 把GLDAS网格重采样到GRACE网格空间对齐推荐的做法是把较高分辨率的GLDAS数据重采样到GRACE mascon的低分辨率网格。常见的插值函数是interp2但要注意GLDAS经度通常是0到360而GRACE可能是-180到180需要先统一坐标。% 把GLDAS 0.25度数据插值到GRACE 0.5度网格 [glon_grid, glat_grid] meshgrid(gldas_lon, gldas_lat); [glon_target, glat_target] meshgrid(grace_lon, grace_lat); gldas_regrid interp2(glon_grid, glat_grid, sm0_month, ... glon_target, glat_target, linear);这段代码里sm0_month的维度是(lon, lat)转置是因为meshgrid在Matlab里第一个输出是沿x方向扩展和数据的行列顺序要对清楚。插值完成后一定要加一步陆海掩膜把海洋区域重新设成NaN不然海岸带会出现插值导致的阶梯状伪值。4.3 月份对不上怎么办缺测月与重采样GRACE在2011年初、2016年底等时段有大量缺失月份GLDAS则相对连续。残差法要求两个数据集在时间上一一对应缺了的GRACE月份GLDAS那几个月也必须去掉否则趋势拟合会产生系统性偏移。处理方式是先读GRACE时间列表再根据该列表去筛选GLDAS月度数据。% 假设graceTime是GRACE月份序列gldasTime是GLDAS月份序列 validIdx ismember(gldasTime, graceTime); gldas_filt gldas_all(:,:,validIdx);这里有一个更隐蔽的问题GLDAS月平均文件的时间戳一般取该月最后一天或下月第一天GRACE mascon时间戳通常是月中或月首。如果直接用datenum匹配很可能全部错开。我建议只比较年和月不考虑具体日用MATLAB的year和month函数提取后再匹配。gldasMonthNum year(gldasTime) * 12 month(gldasTime); graceMonthNum year(graceTime) * 12 month(graceTime); [validIdx] ismember(gldasMonthNum, graceMonthNum);5. 残差法核心计算从逐网格GWS到区域时间序列5.1 逐像元计算地下水储量变化当GLDAS和GRACE在时间、空间、单位三个维度上都对齐之后残差法本身就是一个很简单的减法。我把GLDAS地表储水总量定义为土壤水四层之和、雪水当量、冠层截留水三者相加记住是三部分不是把径流蒸散发也混进来。% GLDAS地表储水总量 gldas_storage_mm sm0_month_mm sm10_month_mm sm40_month_mm ... sm100_month_mm swe_month_mm canop_month_mm; % 残差法计算地下水储量变化 GWS_anomaly_mm GRACE_TWSC_mm - gldas_storage_mm;这里需要明确一个概念GRACE TWSC和GLDAS存储量本身都有一个长期平均基线两个基线的绝对量没有可比性所以直接相减得到的GWS是一个异常值序列通常需要扣除该序列的时间均值得到相对于平均态的偏离。GWS_anomaly_detrended GWS_anomaly_mm - mean(GWS_anomaly_mm, 3, omitnan);5.2 区域平均与异常值计算逐像元GWS结果通常空间噪声还比较大实际研究中更多是分析区域平均时间序列。区域平均就是对研究区内所有有效像元求空间平均注意每个像元的纬度权重。因为等经纬度网格中高纬度像元代表的实际面积更小如果研究区跨越的纬度范围超过几度建议按纬度的余弦加权。% 纬度权重 latWeight cosd(grace_lat); latWeight2D repmat(latWeight, 1, length(grace_lon)); % 区域平均 gws_region squeeze(mean(GWS_anomaly_detrended .* latWeight2D, [1 2], omitnan) ./ ... mean(latWeight2D, [1 2]));5.3 趋势和季节项的拟合与结果输出区域GWS时间序列里同时包含长期趋势、年周期、半年周期和残余噪声。要单独提取地下水储量变化的长期趋势一般是做最小二乘拟合模型包含线性趋势加上年周期和半年周期的正弦余弦项TWSC(t) a b·t c₁·cos(2πt) s₁·sin(2πt) c₂·cos(4πt) s₂·sin(4πt)这个模型写进Matlab可以用简单的回归矩阵实现t (1:length(gws_region)); X [ones(length(t),1), t, ... cos(2*pi*t/12), sin(2*pi*t/12), ... cos(4*pi*t/12), sin(4*pi*t/12)]; b X \ gws_region; trend_mm_per_month b(2); trend_mm_per_year trend_mm_per_month * 12;b的第二个元素就是每月变化速率乘以12得到年变化速率。很多文献里报告的地下水储量变化速率就是单位mm/year就是这个值。拟合之前记得剔除NaN月份否则回归矩阵无法求解。输出结果时我习惯把逐像元GWS保存为NetCDF或GeoTIFF区域平均序列保存为CSV方便后续画图和进一步分析。6. 跑数据过程中我踩过的几个典型坑6.1 重复累加土壤水分的四层问题最初我在写GLDAS脚本时看到NOAH的土壤湿度有四层就习惯性把四层相加作为总土壤水。后来核对数据时发现某些版本的产品除了分层土壤湿度之外还会直接给出一个总土壤水变量有的还带SoilMoi0_100cm_inst这种从0到1米的整合量。如果手滑把分层相加的结果又加上总土壤水那地下水储量变化就被系统性高估了。建议每次处理前先用ncinfo确认变量列表不要凭经验想当然。6.2 插值后海洋像元出现NaN用interp2把GLDAS插值到GRACE网格时如果GLDAS的原始陆地掩膜和GRACE的掩膜不一致海岸线附近的像元可能在插值后变成NaN也有可能被插出一个跨越海洋的假信号。尤其是和GRACE相减后那个NaN会被带入GWS结果区域平均时如果不加omitnan整个时间序列就全是NaN。我后来一律在插值完成后立刻用isnan检查覆盖范围再用GRACE自身的海陆掩膜做一次过滤。6.3 GLDAS和GRACE的版本差异对结论的影响同一个研究区如果用GLDAS-2.0和GLDAS-2.1分别处理得到的土壤水长期趋势可能是不同的尤其在植被覆盖和灌溉影响显著的地区。GRACE这边JPL mascon和CSR mascon的结果在某些区域也会有几个毫米每年的差异。这不是程序写错了而是不同陆面模型参数化方案和重力场反演处理策略带来的真实不确定性。建议在论文或报告里至少用两组数据做交叉验证不要只跑一条链路就下结论。还有一个小细节GRACE mascon产品中某些像元会带有质量因子或uncertainty字段在提取区域平均时可以把不确定性作为权重。不过如果只是做趋势分析通常不强制要求但至少应该看一眼误差量级避免结果被个别异常像元主导。最后再分享一个处理技巧如果你和我一样需要反复调试参数建议把整个流程拆成几个独立的函数readGLDASMonth、regridToGrace、computeGWS、fitTrend。每个函数输入输出都写在注释里路径和数据版本用变量统一管理。这样每次换研究区或换GRACE产品时只需要改输入参数不需要重写核心逻辑。我初期把所有步骤写在同一个脚本里改一次参数从头跑一遍浪费了很多时间。拆成函数之后效率提升非常明显。另外处理完的数据一定要存一份中间结果尤其是重采样后的GLDAS月平均数据和月份对齐后的GRACE数据。很多问题是在后处理阶段才发现的如果没有中间结果就得重新下载和读取原始文件那时间成本就很可观了。这个流程跑通之后你会发现真正花时间的不是减法那一步而是数据下载、格式检查和版本匹配这些看起来不起眼的前置工作。把这些记录下来后面做GRACE-FO数据的延续分析时就轻松多了。本文还有配套的精品资源点击获取