简介本资源是一套面向地球物理与水文遥感研究者的MATLAB工具集聚焦GRACE重力卫星数据与GLDAS陆面模型的协同分析专为水储量变化反演任务设计适用于具备基础MATLAB编程能力及地球重力场球谐展开知识的科研人员与研究生。压缩包含16个文件以9个核心MATLAB脚本如main.m、gravityDisturbance_fast.m、totalWaterStorage_fast.m为主体辅以2个GRACE球谐系数gfc文件、2份球谐函数理论PDF讲义、1个海岸线dat数据及配套txt参数说明整体仅1.52MB轻量易部署。已有796人学习下载提供从GLDAS数据读取、勒让德函数计算、重力扰动快速求解、大地水准面校正到总水储量反演的完整处理链代码模块清晰、注释充分并包含filterCoefficientsGaussian.m等实用预处理工具显著降低GRACE水文信号提取的技术门槛。 直接说项目背景。做水文遥感或水资源方向的人大概率都绕不开GRACE——那个从太空称重地球水资源的卫星任务。而GRACE水储量解算这个活儿最有意思的地方在于它从来不是一套数据能单打独斗的。尤其是当你需要把重力信号中的地表水、土壤水、雪水这部分区分出来时GLDAS陆地同化数据几乎是绕不过去的搭档。这个项目的标题一上来就很直白read_gldasGLDAS读取后面还跟了个IWant!IWant_use那种迫切感就是我在第一次动手处理GLDAS数据时的心情——网上下载的原始文件是净CDF格式变量名和单位都跟论文里对不上时间步长还要自己拼真想立刻有一套能直接跑通的读取代码。这篇文章就把这套流程掰开讲清楚GRACE水储量解算是什么原理为什么一定要读GLDASGLDAS数据从哪里下载、怎么读、怎么转换单位以及GRACE和GLDAS怎么配合做水储量组分解算。内容面向正在被数据处理折磨的师弟师妹也面向初次接触卫星水文学的从业者保证是能直接照着操作的经验总结。1. 项目整体解读GRACE水储量解算到底在算什么1.1 GRACE卫星的核心原理GRACE全称是Gravity Recovery and Climate Experiment中文通常叫重力恢复与气候实验。2002年发射2017年退役后来由GRACE-FO接棒。它的原理并不复杂但非常巧妙两颗卫星在同一轨道上前后飞行相距大概220公里。当地球重力场发生局部变化时前面的卫星会先产生微小加速或减速导致两颗卫星之间的距离发生变化。通过连续测量这个距离变化科学家就可以反演出地球重力场的时空变化。而重力场的这些变化在地表上最主要的体现之一就是陆地水储量的变化。水是有质量的哪个地方地下水位上升了、或者某个流域里蓄了大量地表水那个地方上方的重力信号就会增强。所以GRACE做的水储量解算本质上是把重力异常转换成等效水高Equivalent Water Height, EWH单位是毫米或厘米物理意义就是某个区域在某个时间段内平均每平方米面积上的水量相对于长期基准值高了多少毫米。需要强调一点GRACE结果不是点上的数据它的空间分辨率很粗一般在300公里左右。所以它给出的水储量变化是一个大区域的平均情况比如整个华北平原、整个长江流域这种尺度。它跟站点观测是两回事但正好适合用来分析区域尺度的水资源变化。1.2 为什么这一步要引入GLDASGRACE测的是总水储量变化也就是从地表到地下深处所有水的总和变化。但我们在水资源分析里经常需要知道具体是哪一部分水在变是土壤水增加了还是地下水在超采还是高纬度地区积雪在减少GRACE自己区分不了这些。这时候就该GLDAS出场了。GLDAS是Global Land Data Assimilation System的缩写由NASA戈达德空间飞行中心主导它做的事情是把陆面过程模型如Noah、CLM和观测数据降水、辐射、气温等同化在一起输出全球陆地表面的土壤湿度、雪水当量、冠层水、蒸发、径流等变量。所以一个很经典的分离方案是用GRACE得到的总水储量变化减去GLDAS给出的地表以上部分和根区土壤水的贡献剩下的残差主要就是地下水储量变化。这在地下水补给、枯水期水储量评估等研究里非常常用。用数学表达就是(\Delta GWSA \Delta TWSA_{GRACE} - (\Delta SM_{GLDAS} \Delta SWE_{GLDAS} \Delta Canop_{GLDAS}))其中TWSA是总水储量异常SM是土壤水SWE是雪水当量Canop是冠层水。1.3 read_gldas在整个流程中的角色看完上面的公式你应该明白了GRACE水储量解算里GLDAS不是可有可无的陪衬而是分离地下水信号的核心输入。所以read_gldas这个模块是整个流程的第一步也是决定后面所有分析质量的关键一环。在实际项目中我见过不少人在这一步翻车。GLDAS官网提供的原始数据是NetCDF格式文件名里带了一长串时间戳和版本号变量名是SoilMoi0_10cm_inst这种看起来很直白但实际上需要仔细对照文档才能确认单位。如果照着别人的博客代码读进来就直接用很可能踩到单位换算的坑——比如GLDAS的土壤水单位其实是kg/m²但如果你当成mm直接用量级会差出一倍以上因为1 kg/m²的水深实际上是1 mm但如果涉及土壤体积含水量转换还要考虑土壤密度和分层厚度一不留神就错了。所以这个项目的read_gldas模块解决的是最基础也最关键的三个问题一是解析NetCDF文件结构和变量二是把不同分层的土壤水数据合并成总土壤水储量三是统一单位并重采样到和GRACE一致的时空格网。后面第三部分我会给出完整的代码框架和操作流程。2. GLDAS数据读取从下载到预处理的关键细节2.1 GLDAS产品版本与格式选择GLDAS目前常用的产品有两个版本GLDAS-2.1和GLDAS-2.2。如果做水和旱分析2.1系列的用户最多它是基于NOAH陆面模型输出的时间范围从2000年到现在时间分辨率有3小时和月的两种。做水储量解算建议直接用月产品因为GRACE本身就是月分辨率数据两者天然匹配。数据下载方式主要有两种第一种是NASA GES DISC的官网按条件点选下载适合小数据量、偶尔用的人。但如果你是做长时间序列比如2002年到2020年一个月的全球文件大约几十MB全部下载下来有十几个GB一个个点选简直要人命。这时候就得上第二种方式用wget脚本或者Python的requests库基于GES DISC的OpenDAP接口或者HTTPS目录批量下载。注意GES DISC现在要求用户登录HTTP下载链接里需要带上认证信息。建议在GES DISC官网注册一个账号然后生成一个.netrc文件放在用户目录下内容格式是machine urs.earthdata.nasa.gov login 你的用户名 password 你的密码之后用wget加上--auth-no-challenge参数就能实现无人值守批量下载了。2.2 核心变量及其物理含义GLDAS-2.1 NOAH月产品里和GRACE组分解算直接相关的变量主要有以下几项变量名含义单位层次SoilMoi0_10cm_inst0~10cm土壤含水量kg/m²第1层SoilMoi10_40cm_inst10~40cm土壤含水量kg/m²第2层SoilMoi40_100cm_inst40~100cm土壤含水量kg/m²第3层SoilMoi100_200cm_inst100~200cm土壤含水量kg/m²第4层SWE_inst雪水当量kg/m²地表CanopInt_inst冠层拦截水kg/m²植被冠层这里有个很关键的细节GLDAS中土壤水是按4层给的深度范围分别是0~10cm、10~40cm、40~100cm、100~200cm全部加起来才是2米深度的总土壤水储量。在计算水储量变化时如果只用前两层来代表整个土壤水那系统的误差会很大尤其在干旱半干旱地区深层土壤水的变化不可忽略。单位问题也要单独说一遍。GLDAS里这些变量的单位都是kg/m²但很多人在做图的时候会把它直接当成mm。实际上在标准重力条件下1 kg/m²的水就等于1 mm的水深所以数值上两者确实相等。但如果数据经过了缩放或插值处理这个关系可能会被打破所以建议在代码里先除以1000转成标准的等效水高存储也就是米再参与后续计算避免后期集成时单位混乱。2.3 读取代码框架与单位换算下面这套代码是我在项目中反复用过的可以直接作为read_gldas模块的骨架。它负责读取一个月度的GLDAS NetCDF文件提取四个土壤层、雪水和冠层水的变量输出总水储量单位统一为毫米。import netCDF4 as nc import numpy as np import glob import os def read_gldas_file(filepath): 读取单个月度GLDAS-2.1 NOAH NetCDF文件 返回: 总水储量包含土壤水雪水冠层水单位mm ds nc.Dataset(filepath) # GLDAS-2.1月度产品的纬度是从北到南排列的GRACE是南到北 lat ds.variables[lat][:] lon ds.variables[lon][:] # 读取土壤水分层 sm1 ds.variables[SoilMoi0_10cm_inst][0, :, :] # kg/m2 sm2 ds.variables[SoilMoi10_40cm_inst][0, :, :] sm3 ds.variables[SoilMoi40_100cm_inst][0, :, :] sm4 ds.variables[SoilMoi100_200cm_inst][0, :, :] # 雪水和冠层水 swe ds.variables[SWE_inst][0, :, :] canopy ds.variables[CanopInt_inst][0, :, :] # 合并为总水储量单位转为mm total_storage (sm1 sm2 sm3 sm4 swe canopy) # 等效mm # 如果纬度是降序需要翻转让数组是南到北 if lat[0] lat[-1]: total_storage np.flipud(total_storage) lat lat[::-1] ds.close() return lon, lat, total_storage说明几个比较值得注意的细节第一NetCDF文件读取之后建议立刻关闭ds.close()不然文件句柄会在一会儿大批量处理时耗尽。第二GLDAS-2.1的月产品时间维度只有一层所以索引写成[0, :, :]但如果是3小时产品时间维度就是8层写法不一样别直接套。第三经纬度数组最好在读取后立刻确认方向和范围GLDAS纬度默认是从北纬90到南纬-60GRACE数据通常是-90到90如果不统一方向后面做空间重采样时地图方向会出问题这个坑我在初学阶段踩过不止一次。3. GRACE水储量反演完整流程3.1 数据准备与时间配准GRACE官方数据产品主要来自三家机构CSR德克萨斯大学奥斯汀分校、JPL喷气推进实验室、GFZ波茨坦地学研究中心。产品形式有两种第一种是球谐系数Spherical Harmonics格式每个月的文件是一个文本或者NetCDF文件包含几十阶的Stokes系数。这属于最原始的产品使用前要做不少处理去C20项、替换C20、加回C20不同机构得到的结果略有差异替换C20已经成为默认做法、去C21/S21、去地心项等等处理起来异常繁琐。第二种是Mascon产品其中JPL的Mascon RL06是目前最常用的。它以0.5度格网直接给出等效水高异常单位是厘米不需要做Gaussian平滑和去相关滤波因为Mascon算法本身已经通过空间约束解决了GRACE数据的条带误差问题。做组分解算的话首选Mascon产品因为它的格式和GLDAS差距小处理起来直接。JPL Mascon数据的NetCDF文件里有一个变量叫lwe_thickness它已经是等效水高异常了单位是cm直接除以100就能转成m几乎可以说是开箱即用。然后做时间配准。GRACE月产品的时间标签是每个月的中间时刻比如2005年3月的GRACE数据标示的是2005-03-15而GLDAS月产品时间是2005-03-01或者2005-03-31时间上略有偏移。但从月尺度来看这种差异对水储量变化的计算影响很小因为半个月的偏差在水文信号里基本可以忽略。真正要注意的是GRACE在部分月份有数据缺失比如2011年、2012年、2016年某些时段因电池问题无数据这些月份的处理策略要么直接剔除要么用插值填补建议直接用线性插值补但注意补出的数据不能用于趋势反转分析只能作为展示用。3.2 水储量异常的计算与滤波如果是用JPL Mascon数据水储量异常直接在数据里有了不需要额外计算。但如果是用球谐系数产品需要经过下面几个步骤第一步从低阶Stokes系数中提取每个月相对长期平均的变化量。第二步对球谐系数做Gaussian平滑。Gaussian平滑半径一般取300km左右它相当于一个低通滤波器把GRACE信号中的高频噪声主要是南北向条带滤掉。第三步做去相关滤波最常用的是PCA方法或者Swenson-Wahr提出的经验正交函数滤波。做完这两步再通过球谐综合求逆得到每个格网上的等效水高异常。滤波半径的选取是个需要反复权衡的事。选大一点信号更平滑噪声更小但真实信号也会被压扁尤其是一些小流域的储水信号可能直接被抹掉。选小一点空间分辨率更高但条带噪声会更明显。我的经验是如果目标是全国或大区域尺度300km就够如果关注某个具体的小流域可以试试250km同时对比不同滤波半径的结果看关键结论是否稳定。Mascon产品不需要这些操作但有它自己的问题南北纬60度以上地区的质量守恒约束做得比较紧冰川均衡调整GIA信号剔除也比较依赖模型所以在高纬度应用时要额外小心建议去读一下JPL Mascon的文档了解它适用的地理范围。3.3 GRACE和GLDAS结果的对比验证拿到GRACE和GLDAS两个水储量序列之后第一件事不是急着相减而是先做对比验证。做法是把GLDAS总水储量序列和GRACE总水储量异常序列画在同一张图里看它们的变化趋势是否一致。正常情况下在湿润地区两者应该有很强的相关性因为GRACE测到的总水储量变化大部分由土壤水和地表水主导而在干旱半干旱地区由于地下水占比大GLDAS里没有地下水模块所以两者差异会比较明显甚至趋势都可能相反。这种差异本身就是一个重要信号它意味着该区域地下水正在发生显著变化。以我处理过的一个华北地区项目为例GRACE显示该地区2003到2015年总水储量平均每年下降约37毫米等效水高而GLDAS给出的土壤水雪水变化几乎是一条水平线两者相减之后的地下水储量变化就成了主导每年下降约35毫米左右。这个结果和当地地下水监测井数据的趋势基本吻合验证之后才敢放心把结果用于水资源变化报告。对比的时候还有一个参数要注意相关系数和标准化均方根误差NRMSE都要算。其中相关系数反映的是两个时间序列的变化同步性NRMSE反映的是量级上的差异很多时候两个序列相关系数很高都超过0.8但量级差得很远这说明GRACE和GLDAS的绝对水量基准不一样但变化趋势一致。这种趋势一致、基准不同的情况在做残差分析时完全OK因为相减的时候偏差就自动消掉了。4. 常见问题与排错实录4.1 单位与量级问题我见过最多的错误就是单位问题。一个很经典的案例是某同学从GRACE官网下载了等效水高数据单位是cm但他误认为mm直接填进模型里计算速率结果趋势值大了10倍导致水资源评价结论完全失真。GLDAS这边也一样土壤水分层变量的单位是kg/m²雪水当量也是kg/m²这两者直接相加没有问题。但如果你从GES DISC下载的GLDAS-2.2产品单位可能是kg m-2 s-1通量型变量那含义就完全不同了不能直接用来算储量得先做时间积分。还有一个比较隐蔽的坑JPL Mascon的lwe_thickness变量的scale factor。NetCDF文件里通常还带一个scale_factor属性某些版本需要乘上这个因子才是真实的等效水高。如果你的结果量级看起来不对——比如一个盆地级的月异常居然有几百毫米——先检查是不是漏了这个因子。4.2 格网投影与空间分辨率匹配GRACE Mascon格网是0.5度经纬度规则格网GLDAS-2.1也是0.25度规则格网所以两者需要在统一格网上对齐。最常见的方法是把GLDAS从0.25度降采样到0.5度。降采样要注意方法的选择。我建议用面积加权平均area-weighted average或者简单的双线性插值而不是直接取邻近点。因为取邻近点会丢失大量空间信息尤其是GLDAS的高分辨率土壤水细节在降采样后会被抹掉这会导致最终的残差序列方差偏大。实际处理时我通常用xarray的interp方法或者scipy的griddata来做操作很简单代码就几行import xarray as xr import numpy as np # gldas_ds 是xarray DataArray0.25度 gldas_05 gldas_ds.interp(latnew_lat_05, lonnew_lon_05, methodlinear)需要注意xarray的interp默认用线性插值对规则格网完全够用。插值完成后最好跟原始0.25度数据做一下对比确认总水量守恒——也就是插值前后区域均值变化不超过1%。4.3 数据缺失与时间序列拼接GRACE在17年的寿命中有大约一百个月的正常数据但有十几个月只有单星观测数据质量较差通常不建议使用。GLDAS的时间直接从2000年到现在基本无缝所以时间拼接主要问题出在GRACE不是最大的问题出在GRACE和GLDAS时间基准不一致上。GRACE水储量异常是基于2004年1月到2009年12月的平均基线算出来的而你自己基于GLDAS计算的水储量异常是基于你自己选的时间段平均。如果把两个序列直接放在一起比较整体会有一个系统性偏移这个偏移会让两者绝对值对不上。解决办法是在计算GLDAS水储量异常时选和GRACE一致的基线时段2004-2009并把两者都做去均值处理。这样两者比较的可信度才高。做完去均值后再做相关性分析就不会出现明明趋势一样但两条线隔着十万八千里的尴尬局面。4.4 常见问题速查表症状可能原因解决办法GRACE结果量级偏大10倍单位cm误认为mm先确认产品文档统一换算成mmGLDAS总土壤水与GRACE序列趋势相反GLDAS没有地下水模块区域地下水主导变化做差分离地下水单独分析不要直接拿GLDAS与GRACE比较空间插值后区域均值变化超过5%插值方法不合适或使用了取邻近点改用面积加权平均或双线性插值并做水量守恒检查GRACE有效月份数量太少单星观测期数据质量差剔除单星数据或用插值补齐但补齐数据不用于趋势分析截取子区域时方向颠倒纬度数组降序统一翻转确保所有数据集是南到北排列两条序列趋势相同但绝对值差很多时间基准不一致选取相同基线时段各自去均值后再比较4.5 经验总结实战中的几条体会再多说几句实操体会。第一永远把原始数据处理放在最前面。我曾经因为偷懒直接在Excel里手动改了一列GRACE数据后来发现有个数值的小数点看错了导致整个区域的趋势分析报废。数据处理流程要写脚本、要可复现不要手动改。第二GRACE和GLDAS是互补关系不是替代关系。GRACE的优势在总水储量变化GLDAS的优势在组分分离两者的分辨率、精度、误差来源都不同。不要指望GLDAS能替代GRACE做地下水监测但也别指望GRACE能告诉你地下水在上层还是下层。第三如果做时间跨度非常长的分析比如2002-2020建议把GRACE、GRACE-FO两代卫星的数据都处理进来同时在拼接处做重叠时段的交叉验证确保两代卫星数据趋势一致。GRACE和GRACE-FO之间有大约一年的数据空缺2017年下半年到2018年上半年这个空缺只能插值没别的办法。第四建议把数据版本说明写进代码注释和论文的方法部分。GRACE数据版本RL05还是RL06、GLDAS版本2.1还是2.2、滤波半径、插值方法这些细节直接决定了最终结论的数值。读者在复现你的工作时第一个要确认的就是这些参数。最后分享一个小技巧。在写read_gldas的时候我习惯把读取逻辑和计算逻辑拆成两个模块读取模块只管读取和维度统一计算模块只负责水量计算和单位转换。这样当你需要换GLDAS版本、或者换GRACE产品时只需要改读取模块计算模块完全不用动。这个设计让我后来在对比不同GLDAS版本时少改了大量代码实测下来非常省心。本文还有配套的精品资源点击获取 GRACE水储量解算中的GLDAS数据读取与处理实战指南
简介本资源是一套面向地球物理与水文遥感研究者的MATLAB工具集聚焦GRACE重力卫星数据与GLDAS陆面模型的协同分析专为水储量变化反演任务设计适用于具备基础MATLAB编程能力及地球重力场球谐展开知识的科研人员与研究生。压缩包含16个文件以9个核心MATLAB脚本如main.m、gravityDisturbance_fast.m、totalWaterStorage_fast.m为主体辅以2个GRACE球谐系数gfc文件、2份球谐函数理论PDF讲义、1个海岸线dat数据及配套txt参数说明整体仅1.52MB轻量易部署。已有796人学习下载提供从GLDAS数据读取、勒让德函数计算、重力扰动快速求解、大地水准面校正到总水储量反演的完整处理链代码模块清晰、注释充分并包含filterCoefficientsGaussian.m等实用预处理工具显著降低GRACE水文信号提取的技术门槛。 直接说项目背景。做水文遥感或水资源方向的人大概率都绕不开GRACE——那个从太空称重地球水资源的卫星任务。而GRACE水储量解算这个活儿最有意思的地方在于它从来不是一套数据能单打独斗的。尤其是当你需要把重力信号中的地表水、土壤水、雪水这部分区分出来时GLDAS陆地同化数据几乎是绕不过去的搭档。这个项目的标题一上来就很直白read_gldasGLDAS读取后面还跟了个IWant!IWant_use那种迫切感就是我在第一次动手处理GLDAS数据时的心情——网上下载的原始文件是净CDF格式变量名和单位都跟论文里对不上时间步长还要自己拼真想立刻有一套能直接跑通的读取代码。这篇文章就把这套流程掰开讲清楚GRACE水储量解算是什么原理为什么一定要读GLDASGLDAS数据从哪里下载、怎么读、怎么转换单位以及GRACE和GLDAS怎么配合做水储量组分解算。内容面向正在被数据处理折磨的师弟师妹也面向初次接触卫星水文学的从业者保证是能直接照着操作的经验总结。1. 项目整体解读GRACE水储量解算到底在算什么1.1 GRACE卫星的核心原理GRACE全称是Gravity Recovery and Climate Experiment中文通常叫重力恢复与气候实验。2002年发射2017年退役后来由GRACE-FO接棒。它的原理并不复杂但非常巧妙两颗卫星在同一轨道上前后飞行相距大概220公里。当地球重力场发生局部变化时前面的卫星会先产生微小加速或减速导致两颗卫星之间的距离发生变化。通过连续测量这个距离变化科学家就可以反演出地球重力场的时空变化。而重力场的这些变化在地表上最主要的体现之一就是陆地水储量的变化。水是有质量的哪个地方地下水位上升了、或者某个流域里蓄了大量地表水那个地方上方的重力信号就会增强。所以GRACE做的水储量解算本质上是把重力异常转换成等效水高Equivalent Water Height, EWH单位是毫米或厘米物理意义就是某个区域在某个时间段内平均每平方米面积上的水量相对于长期基准值高了多少毫米。需要强调一点GRACE结果不是点上的数据它的空间分辨率很粗一般在300公里左右。所以它给出的水储量变化是一个大区域的平均情况比如整个华北平原、整个长江流域这种尺度。它跟站点观测是两回事但正好适合用来分析区域尺度的水资源变化。1.2 为什么这一步要引入GLDASGRACE测的是总水储量变化也就是从地表到地下深处所有水的总和变化。但我们在水资源分析里经常需要知道具体是哪一部分水在变是土壤水增加了还是地下水在超采还是高纬度地区积雪在减少GRACE自己区分不了这些。这时候就该GLDAS出场了。GLDAS是Global Land Data Assimilation System的缩写由NASA戈达德空间飞行中心主导它做的事情是把陆面过程模型如Noah、CLM和观测数据降水、辐射、气温等同化在一起输出全球陆地表面的土壤湿度、雪水当量、冠层水、蒸发、径流等变量。所以一个很经典的分离方案是用GRACE得到的总水储量变化减去GLDAS给出的地表以上部分和根区土壤水的贡献剩下的残差主要就是地下水储量变化。这在地下水补给、枯水期水储量评估等研究里非常常用。用数学表达就是(\Delta GWSA \Delta TWSA_{GRACE} - (\Delta SM_{GLDAS} \Delta SWE_{GLDAS} \Delta Canop_{GLDAS}))其中TWSA是总水储量异常SM是土壤水SWE是雪水当量Canop是冠层水。1.3 read_gldas在整个流程中的角色看完上面的公式你应该明白了GRACE水储量解算里GLDAS不是可有可无的陪衬而是分离地下水信号的核心输入。所以read_gldas这个模块是整个流程的第一步也是决定后面所有分析质量的关键一环。在实际项目中我见过不少人在这一步翻车。GLDAS官网提供的原始数据是NetCDF格式文件名里带了一长串时间戳和版本号变量名是SoilMoi0_10cm_inst这种看起来很直白但实际上需要仔细对照文档才能确认单位。如果照着别人的博客代码读进来就直接用很可能踩到单位换算的坑——比如GLDAS的土壤水单位其实是kg/m²但如果你当成mm直接用量级会差出一倍以上因为1 kg/m²的水深实际上是1 mm但如果涉及土壤体积含水量转换还要考虑土壤密度和分层厚度一不留神就错了。所以这个项目的read_gldas模块解决的是最基础也最关键的三个问题一是解析NetCDF文件结构和变量二是把不同分层的土壤水数据合并成总土壤水储量三是统一单位并重采样到和GRACE一致的时空格网。后面第三部分我会给出完整的代码框架和操作流程。2. GLDAS数据读取从下载到预处理的关键细节2.1 GLDAS产品版本与格式选择GLDAS目前常用的产品有两个版本GLDAS-2.1和GLDAS-2.2。如果做水和旱分析2.1系列的用户最多它是基于NOAH陆面模型输出的时间范围从2000年到现在时间分辨率有3小时和月的两种。做水储量解算建议直接用月产品因为GRACE本身就是月分辨率数据两者天然匹配。数据下载方式主要有两种第一种是NASA GES DISC的官网按条件点选下载适合小数据量、偶尔用的人。但如果你是做长时间序列比如2002年到2020年一个月的全球文件大约几十MB全部下载下来有十几个GB一个个点选简直要人命。这时候就得上第二种方式用wget脚本或者Python的requests库基于GES DISC的OpenDAP接口或者HTTPS目录批量下载。注意GES DISC现在要求用户登录HTTP下载链接里需要带上认证信息。建议在GES DISC官网注册一个账号然后生成一个.netrc文件放在用户目录下内容格式是machine urs.earthdata.nasa.gov login 你的用户名 password 你的密码之后用wget加上--auth-no-challenge参数就能实现无人值守批量下载了。2.2 核心变量及其物理含义GLDAS-2.1 NOAH月产品里和GRACE组分解算直接相关的变量主要有以下几项变量名含义单位层次SoilMoi0_10cm_inst0~10cm土壤含水量kg/m²第1层SoilMoi10_40cm_inst10~40cm土壤含水量kg/m²第2层SoilMoi40_100cm_inst40~100cm土壤含水量kg/m²第3层SoilMoi100_200cm_inst100~200cm土壤含水量kg/m²第4层SWE_inst雪水当量kg/m²地表CanopInt_inst冠层拦截水kg/m²植被冠层这里有个很关键的细节GLDAS中土壤水是按4层给的深度范围分别是0~10cm、10~40cm、40~100cm、100~200cm全部加起来才是2米深度的总土壤水储量。在计算水储量变化时如果只用前两层来代表整个土壤水那系统的误差会很大尤其在干旱半干旱地区深层土壤水的变化不可忽略。单位问题也要单独说一遍。GLDAS里这些变量的单位都是kg/m²但很多人在做图的时候会把它直接当成mm。实际上在标准重力条件下1 kg/m²的水就等于1 mm的水深所以数值上两者确实相等。但如果数据经过了缩放或插值处理这个关系可能会被打破所以建议在代码里先除以1000转成标准的等效水高存储也就是米再参与后续计算避免后期集成时单位混乱。2.3 读取代码框架与单位换算下面这套代码是我在项目中反复用过的可以直接作为read_gldas模块的骨架。它负责读取一个月度的GLDAS NetCDF文件提取四个土壤层、雪水和冠层水的变量输出总水储量单位统一为毫米。import netCDF4 as nc import numpy as np import glob import os def read_gldas_file(filepath): 读取单个月度GLDAS-2.1 NOAH NetCDF文件 返回: 总水储量包含土壤水雪水冠层水单位mm ds nc.Dataset(filepath) # GLDAS-2.1月度产品的纬度是从北到南排列的GRACE是南到北 lat ds.variables[lat][:] lon ds.variables[lon][:] # 读取土壤水分层 sm1 ds.variables[SoilMoi0_10cm_inst][0, :, :] # kg/m2 sm2 ds.variables[SoilMoi10_40cm_inst][0, :, :] sm3 ds.variables[SoilMoi40_100cm_inst][0, :, :] sm4 ds.variables[SoilMoi100_200cm_inst][0, :, :] # 雪水和冠层水 swe ds.variables[SWE_inst][0, :, :] canopy ds.variables[CanopInt_inst][0, :, :] # 合并为总水储量单位转为mm total_storage (sm1 sm2 sm3 sm4 swe canopy) # 等效mm # 如果纬度是降序需要翻转让数组是南到北 if lat[0] lat[-1]: total_storage np.flipud(total_storage) lat lat[::-1] ds.close() return lon, lat, total_storage说明几个比较值得注意的细节第一NetCDF文件读取之后建议立刻关闭ds.close()不然文件句柄会在一会儿大批量处理时耗尽。第二GLDAS-2.1的月产品时间维度只有一层所以索引写成[0, :, :]但如果是3小时产品时间维度就是8层写法不一样别直接套。第三经纬度数组最好在读取后立刻确认方向和范围GLDAS纬度默认是从北纬90到南纬-60GRACE数据通常是-90到90如果不统一方向后面做空间重采样时地图方向会出问题这个坑我在初学阶段踩过不止一次。3. GRACE水储量反演完整流程3.1 数据准备与时间配准GRACE官方数据产品主要来自三家机构CSR德克萨斯大学奥斯汀分校、JPL喷气推进实验室、GFZ波茨坦地学研究中心。产品形式有两种第一种是球谐系数Spherical Harmonics格式每个月的文件是一个文本或者NetCDF文件包含几十阶的Stokes系数。这属于最原始的产品使用前要做不少处理去C20项、替换C20、加回C20不同机构得到的结果略有差异替换C20已经成为默认做法、去C21/S21、去地心项等等处理起来异常繁琐。第二种是Mascon产品其中JPL的Mascon RL06是目前最常用的。它以0.5度格网直接给出等效水高异常单位是厘米不需要做Gaussian平滑和去相关滤波因为Mascon算法本身已经通过空间约束解决了GRACE数据的条带误差问题。做组分解算的话首选Mascon产品因为它的格式和GLDAS差距小处理起来直接。JPL Mascon数据的NetCDF文件里有一个变量叫lwe_thickness它已经是等效水高异常了单位是cm直接除以100就能转成m几乎可以说是开箱即用。然后做时间配准。GRACE月产品的时间标签是每个月的中间时刻比如2005年3月的GRACE数据标示的是2005-03-15而GLDAS月产品时间是2005-03-01或者2005-03-31时间上略有偏移。但从月尺度来看这种差异对水储量变化的计算影响很小因为半个月的偏差在水文信号里基本可以忽略。真正要注意的是GRACE在部分月份有数据缺失比如2011年、2012年、2016年某些时段因电池问题无数据这些月份的处理策略要么直接剔除要么用插值填补建议直接用线性插值补但注意补出的数据不能用于趋势反转分析只能作为展示用。3.2 水储量异常的计算与滤波如果是用JPL Mascon数据水储量异常直接在数据里有了不需要额外计算。但如果是用球谐系数产品需要经过下面几个步骤第一步从低阶Stokes系数中提取每个月相对长期平均的变化量。第二步对球谐系数做Gaussian平滑。Gaussian平滑半径一般取300km左右它相当于一个低通滤波器把GRACE信号中的高频噪声主要是南北向条带滤掉。第三步做去相关滤波最常用的是PCA方法或者Swenson-Wahr提出的经验正交函数滤波。做完这两步再通过球谐综合求逆得到每个格网上的等效水高异常。滤波半径的选取是个需要反复权衡的事。选大一点信号更平滑噪声更小但真实信号也会被压扁尤其是一些小流域的储水信号可能直接被抹掉。选小一点空间分辨率更高但条带噪声会更明显。我的经验是如果目标是全国或大区域尺度300km就够如果关注某个具体的小流域可以试试250km同时对比不同滤波半径的结果看关键结论是否稳定。Mascon产品不需要这些操作但有它自己的问题南北纬60度以上地区的质量守恒约束做得比较紧冰川均衡调整GIA信号剔除也比较依赖模型所以在高纬度应用时要额外小心建议去读一下JPL Mascon的文档了解它适用的地理范围。3.3 GRACE和GLDAS结果的对比验证拿到GRACE和GLDAS两个水储量序列之后第一件事不是急着相减而是先做对比验证。做法是把GLDAS总水储量序列和GRACE总水储量异常序列画在同一张图里看它们的变化趋势是否一致。正常情况下在湿润地区两者应该有很强的相关性因为GRACE测到的总水储量变化大部分由土壤水和地表水主导而在干旱半干旱地区由于地下水占比大GLDAS里没有地下水模块所以两者差异会比较明显甚至趋势都可能相反。这种差异本身就是一个重要信号它意味着该区域地下水正在发生显著变化。以我处理过的一个华北地区项目为例GRACE显示该地区2003到2015年总水储量平均每年下降约37毫米等效水高而GLDAS给出的土壤水雪水变化几乎是一条水平线两者相减之后的地下水储量变化就成了主导每年下降约35毫米左右。这个结果和当地地下水监测井数据的趋势基本吻合验证之后才敢放心把结果用于水资源变化报告。对比的时候还有一个参数要注意相关系数和标准化均方根误差NRMSE都要算。其中相关系数反映的是两个时间序列的变化同步性NRMSE反映的是量级上的差异很多时候两个序列相关系数很高都超过0.8但量级差得很远这说明GRACE和GLDAS的绝对水量基准不一样但变化趋势一致。这种趋势一致、基准不同的情况在做残差分析时完全OK因为相减的时候偏差就自动消掉了。4. 常见问题与排错实录4.1 单位与量级问题我见过最多的错误就是单位问题。一个很经典的案例是某同学从GRACE官网下载了等效水高数据单位是cm但他误认为mm直接填进模型里计算速率结果趋势值大了10倍导致水资源评价结论完全失真。GLDAS这边也一样土壤水分层变量的单位是kg/m²雪水当量也是kg/m²这两者直接相加没有问题。但如果你从GES DISC下载的GLDAS-2.2产品单位可能是kg m-2 s-1通量型变量那含义就完全不同了不能直接用来算储量得先做时间积分。还有一个比较隐蔽的坑JPL Mascon的lwe_thickness变量的scale factor。NetCDF文件里通常还带一个scale_factor属性某些版本需要乘上这个因子才是真实的等效水高。如果你的结果量级看起来不对——比如一个盆地级的月异常居然有几百毫米——先检查是不是漏了这个因子。4.2 格网投影与空间分辨率匹配GRACE Mascon格网是0.5度经纬度规则格网GLDAS-2.1也是0.25度规则格网所以两者需要在统一格网上对齐。最常见的方法是把GLDAS从0.25度降采样到0.5度。降采样要注意方法的选择。我建议用面积加权平均area-weighted average或者简单的双线性插值而不是直接取邻近点。因为取邻近点会丢失大量空间信息尤其是GLDAS的高分辨率土壤水细节在降采样后会被抹掉这会导致最终的残差序列方差偏大。实际处理时我通常用xarray的interp方法或者scipy的griddata来做操作很简单代码就几行import xarray as xr import numpy as np # gldas_ds 是xarray DataArray0.25度 gldas_05 gldas_ds.interp(latnew_lat_05, lonnew_lon_05, methodlinear)需要注意xarray的interp默认用线性插值对规则格网完全够用。插值完成后最好跟原始0.25度数据做一下对比确认总水量守恒——也就是插值前后区域均值变化不超过1%。4.3 数据缺失与时间序列拼接GRACE在17年的寿命中有大约一百个月的正常数据但有十几个月只有单星观测数据质量较差通常不建议使用。GLDAS的时间直接从2000年到现在基本无缝所以时间拼接主要问题出在GRACE不是最大的问题出在GRACE和GLDAS时间基准不一致上。GRACE水储量异常是基于2004年1月到2009年12月的平均基线算出来的而你自己基于GLDAS计算的水储量异常是基于你自己选的时间段平均。如果把两个序列直接放在一起比较整体会有一个系统性偏移这个偏移会让两者绝对值对不上。解决办法是在计算GLDAS水储量异常时选和GRACE一致的基线时段2004-2009并把两者都做去均值处理。这样两者比较的可信度才高。做完去均值后再做相关性分析就不会出现明明趋势一样但两条线隔着十万八千里的尴尬局面。4.4 常见问题速查表症状可能原因解决办法GRACE结果量级偏大10倍单位cm误认为mm先确认产品文档统一换算成mmGLDAS总土壤水与GRACE序列趋势相反GLDAS没有地下水模块区域地下水主导变化做差分离地下水单独分析不要直接拿GLDAS与GRACE比较空间插值后区域均值变化超过5%插值方法不合适或使用了取邻近点改用面积加权平均或双线性插值并做水量守恒检查GRACE有效月份数量太少单星观测期数据质量差剔除单星数据或用插值补齐但补齐数据不用于趋势分析截取子区域时方向颠倒纬度数组降序统一翻转确保所有数据集是南到北排列两条序列趋势相同但绝对值差很多时间基准不一致选取相同基线时段各自去均值后再比较4.5 经验总结实战中的几条体会再多说几句实操体会。第一永远把原始数据处理放在最前面。我曾经因为偷懒直接在Excel里手动改了一列GRACE数据后来发现有个数值的小数点看错了导致整个区域的趋势分析报废。数据处理流程要写脚本、要可复现不要手动改。第二GRACE和GLDAS是互补关系不是替代关系。GRACE的优势在总水储量变化GLDAS的优势在组分分离两者的分辨率、精度、误差来源都不同。不要指望GLDAS能替代GRACE做地下水监测但也别指望GRACE能告诉你地下水在上层还是下层。第三如果做时间跨度非常长的分析比如2002-2020建议把GRACE、GRACE-FO两代卫星的数据都处理进来同时在拼接处做重叠时段的交叉验证确保两代卫星数据趋势一致。GRACE和GRACE-FO之间有大约一年的数据空缺2017年下半年到2018年上半年这个空缺只能插值没别的办法。第四建议把数据版本说明写进代码注释和论文的方法部分。GRACE数据版本RL05还是RL06、GLDAS版本2.1还是2.2、滤波半径、插值方法这些细节直接决定了最终结论的数值。读者在复现你的工作时第一个要确认的就是这些参数。最后分享一个小技巧。在写read_gldas的时候我习惯把读取逻辑和计算逻辑拆成两个模块读取模块只管读取和维度统一计算模块只负责水量计算和单位转换。这样当你需要换GLDAS版本、或者换GRACE产品时只需要改读取模块计算模块完全不用动。这个设计让我后来在对比不同GLDAS版本时少改了大量代码实测下来非常省心。本文还有配套的精品资源点击获取