GRACE CSR mascon数据加工：从原始产品到流域水储量时间序列

简介本资源是面向地球物理、冰川学、水文学及气候学等领域科研人员的CSR mascon加工数据集聚焦于GRACE/GRACE-FO卫星反演的地球质量变化分析显著降低重力场数据处理门槛。压缩包共6个文件225.42MB含1个核心NetCDF格式mascon数据文件含时间序列与空间网格化质量变化信息、2个MATLAB脚本分别用于区域格网计算与测试调用、1个MATLAB数据文件预存处理中间结果、2个文本说明文件含长江流域案例参数与使用指南。已有2308人学习下载配套代码完整覆盖nc文件读取、区域提取、时序分析与可视化流程无需额外开发即可开展冰川消融、陆地水储量变化或海平面响应等典型研究特别适合具备基础MATLAB编程能力的研究生与青年科研工作者快速启动课题分析。 开篇先交代一下背景。前两年接了个和GRACE重力卫星数据相关的项目核心任务是做陆地水储量变化的分析。项目里面最绕不开的一环就是怎么把GRACE Level-2球谐系数产品加工成真正能用的区域时间序列。当时在CSR mascon、JPL mascon、GSFC mascon这三个产品之间反复横跳最后选了CSR mascon作为主力数据源顺手把整个加工流程做成了一套标准化数据集。这篇就把我自己的处理路线包括下载、读取、裁剪、掩膜、趋势提取、数据集规范化这些步骤完整梳理出来给后面啃重力卫星数据的朋友做个参考。CSR mascon数据加工数据集从GRACE原始产品到可直接分析的流域水储量时间序列1. 为什么非要自己加工一手CSR mascon数据1.1 GRACE时代成品数据其实分了两条路线很多人第一次接触GRACE、GRACE-FO数据的时候会直接被“球谐系数”“Stokes系数”“去相关滤波”“高斯平滑”这一串名词劝退。确实传统做法是先从官网下载RL06的Level-2 GSM文件拿到一长串球谐系数然后自己去做维度滤波、条带噪声消除、高斯平滑再换算成等效水高。这一套流程跑下来中间任何一个参数选得不对最后的信号都会变形。好在后来各家机构推出了mascon产品也就是“质量集中”反演结果。它的思路不再用球谐系数描述全球重力场而是把地表划分成一个个球冠或网格单元直接求解每个单元的质量变化。对应用户来说拿到手的文件里已经是空间上的质量分布不需要自己再处理球谐展开和滤波省掉了一大截工作量。但产品归产品距离“能用”还是差了一步。CSR发布的mascon原始文件里虽然有全球逐月的质量变化但文件命名、网格坐标、时间轴、单位、参考时段这些细节并不会自动适配你的研究区域。比如我只关心某个流域就得从全球网格里裁出区域再把逐月时间序列整理成DataFrame或者NetCDF还要把GIA、尺度因子、参考时段这些物理含义搞清楚否则求出来的趋势可能差出好几厘米。1.2 CSR mascon相对球谐系数的三大优势先说结论。我选择CSR mascon而不是自己处理球谐系数主要看中三点。第一省掉滤波环节。传统球谐系数产品需要做条带滤波比如Swenson-Wahr滤波然后还要做高斯平滑典型半径150到300公里。这两个操作看起来是降噪实际上也会把真实信号磨平一部分。CSR mascon用球冠基函数加空间约束反演条带噪声在反演过程里就被压住了不用再做高斯平滑空间细节保得更好。第二误差特性更清晰。CSR mascon产品提供了每个网格的误差估计这个在做区域平均和分析显著性时非常有用。传统球谐系数自己做滤波误差传播很难算得准。第三产品自带的物理校正比较完整。CSR的RL06 mascon产品已经处理了地心运动、C20/C30替换、GIA校正等一堆容易踩坑的项用户只要确认自己使用的版本对应的参考时段就能直接进入应用分析省心很多。当然这不代表mascon是万能药。每家mascon产品因为约束和基函数的设定不同结果会有差异。下一篇我会专门对比CSR、JPL、GSFC这里先重点讲讲CSR mascon的数据加工流程。2. 搞到原始产品下载入口、文件类型与目录结构2.1 CSR官网能拿到的文件清单CSR mascon数据目前是公开下载的入口在德克萨斯大学空间研究中心的官网。进入页面后会看到RL06版本的mascon数据集合文件名通常带着CSRM_BA01或者CSR_Mascon_global_v02之类的标记。我习惯按下面的类型去理解这些文件逐月NetCDF文件每个月份一个文件文件名像CSRM_BA01_200204_200230_0004_UTC_sl.nc。这类文件包含全球网格的质量变化。组合时间序列文件把全时段数据打包在一个文件里可能是NetCDF也可能是MAT文件。做长时间趋势分析时比较方便。辅助文件比如球冠定义文件、网格权重文件、GRACE和GRACE-FO衔接说明等这些不直接参与绘图但做数据处理时能用来核对坐标系和掩膜。我第一次下载的时候一度搞不清_sl.nc和_gsm.nc的区别。_sl.nc里的sl是surface load也就是表面质量负载通常已经是等效水高_gsm.nc则更接近重力场模型解里面是Stokes系数或者网格化的重力场变化。实际做水储量分析直接使用_sl.nc即可。2.2 netCDF内部到底存了什么打开一个逐月NetCDF文件里面变量不多但每一个都要确认清楚。常见的变量包括lat、lon全球网格的纬度、经度数组。time时间标记单位一般是从某个参考日期起算的天数。slev等效水高单位可能是厘米或毫米。这个变量是关键后面所有分析和出图都基于它。误差相关变量有的版本会附带误差估计。除了变量全局属性里通常还会写参考时段、GIA模型、单位、数据源等信息。不要小看这些属性后面排查趋势异常的时候全靠它们。有一点要特别提醒CSR的全球mascon网格通常以0.25度或0.5度间隔输出但这不是球谐系数产品那种“真实分辨率”球冠本身的平滑尺度比网格间距大得多。做流域平均时网格大小带来的误差并不等于空间分辨率别用0.25度网格宽度去解释高频细节。2.3 下载前的版本选择RL06、CRI、filtered别搞混CSR mascon页面上会出现好几个版本号或者后缀比如RL06、RL05还可能出现像CRI这样的标记。我第一次就踩过坑下载了RL05的旧文件后面换到RL06后发现趋势差异不小。选版本时我建议直接定RL06因为这是目前广泛验证的版本和GRACE-FO数据的衔接也做得好。如果你要对比GRACE时期2002到2017和GRACE-FO时期2018至今的长期趋势务必使用同一处理版本的数据不要混用RL05和RL06。另外CSR页面可能会区分“filtered”和“unfiltered”数据。mascon产品虽然不需要用户自己滤波但某些后处理文件会先做平滑或者时间域滤波。对大多数分析目的选择常规的、官方推荐给用户做水储量分析的版本即可。如果下载页面没说明建议直接看数据的README或者引用说明那里会写清楚哪个文件是面向最终应用的。3. 核心加工流程从原始文件到可用的水文数据集3.1 环境准备与读取NetCDF我的处理环境是Python 3.10核心库是xarray、numpy、pandas、matplotlib区域绘图还会用到cartopy和geopandas。其中xarray一定要会因为它对带时间维、经纬度维的NetCDF数据支持非常顺滑。读取一个逐月NetCDF文件的代码非常简单import xarray as xr # 以CSR mascon的一个逐月文件为例 ds xr.open_dataset(CSRM_BA01_200204_200230_0004_UTC_sl.nc) print(ds)打印出来的数据结构里能看到维度、坐标和变量。slev变量通常形如(time, lat, lon)如果只有一个时间点time维也可能是空的但坐标信息还在。如果某个月份的文件读取后坐标不是单调递增可以用ds.sortby(lat)处理一下。有些处理工具生成的文件坐标顺序是反的特别是纬度从北极到南极排列后续绘图时会让子区域提取变得混乱。3.2 时间轴与空间轴的重构逐月文件单独看都正常但要做长期分析就必须把几十个甚至一百多个文件拼成一个完整的时间序列。最稳妥的方式是用xarray的open_mfdatasetimport xarray as xr files sorted(glob.glob(CSRM_BA01_*.nc)) ds_all xr.open_mfdataset(files, concat_dimtime, combinenested) ds_all ds_all.sortby(time)这里有个细节。open_mfdataset默认会尝试用坐标对齐但有时候各个文件的lat、lon因为浮点存储精度略有差异导致拼接时报错或产生NaN。稳妥做法是确认所有文件都来自同一版本必要时可以先用ds.load()把单个文件读入内存检查坐标是否一致再批量拼接。时间轴重构也很重要。NetCDF里的time通常是自2002-01-01起算的天数读取后建议显式转换为cftime或pandas的DatetimeIndextime_idx xr.cftime_range(start2002-01, periodsds_all.sizes[time], freqMS) ds_all ds_all.assign_coords(timetime_idx)月度数据的时间戳一般取每月1号。需要注意GRACE数据在某些月份缺失拼接后时间轴会出现空洞。处理时不要直接用resample把空洞填掉先明确哪些月份缺失再决定插值还是保留NaN。3.3 掩膜与区域提取流域尺度的裁剪拿到全球网格后提取特定流域的方式有两种一种是直接用经纬度范围裁剪矩形框另一种是用流域边界矢量做掩膜。矩形裁剪适合初步查看lon_min, lon_max 90, 122 lat_min, lat_max 21, 36 region ds_all.sel( lonslice(lon_min, lon_max), latslice(lat_min, lat_max) )但矩形框会包含大量非目标区域。研究长江流域、黄河流域这种不规则形状时需要把流域边界矢量转成网格掩膜。我是这样做的先读取流域矢量shp文件用geopandas处理再生成一个和全球网格同样形状的布尔掩膜数组import geopandas as gpd import numpy as np basin gpd.read_file(yangtze_basin.shp).to_crs(EPSG:4326) lons ds_all[lon].values lats ds_all[lat].values mask np.zeros((len(lats), len(lons)), dtypebool) # 用点是否落在多边形内来构建掩膜 from shapely.geometry import Point polygon basin.geometry.unary_union for i, lat in enumerate(lats): for j, lon in enumerate(lons): if polygon.contains(Point(lon, lat)): mask[i, j] True这个双重循环在网格很密时会很慢。更高效的办法是用rasterio.features.rasterize直接对矢量栅格化速度能快几个量级。拿到掩膜后就能用where提取区域数据region_ds ds_all.where(mask)然后对空间维度做加权平均得到流域平均的时间序列。3.4 从时间序列去趋势到信号量级分析区域平均得到的是一个随月份变化的等效水高序列。我一般会先画原始曲线看季节波动和年际变化然后用最小二乘拟合一个包含趋势、年周期和半年周期的模型import numpy as np import pandas as pd y region_series.values # 等效水高 t np.arange(len(y)) mask_valid ~np.isnan(y) # 设计矩阵趋势 年周期cos/sin 半年周期 A np.column_stack([ t[mask_valid], np.cos(2*np.pi*t[mask_valid]/12), np.sin(2*np.pi*t[mask_valid]/12), np.cos(4*np.pi*t[mask_valid]/12), np.sin(4*np.pi*t[mask_valid]/12), np.ones(mask_valid.sum()) ]) coef, _, _, _ np.linalg.lstsq(A, y[mask_valid], rcondNone) trend_mm_per_month coef[0] * 30.44 # 由每月趋势换算成mm/月趋势换算成毫米每年时记得乘以月份数。CSR mascon常用单位是cm等效水高算趋势时把cm转成mm乘上10再乘年月份数。这个模型虽然简单但对GRACE时间序列很实用。它能把长期趋势和季节性信号分开方便判断区域整体是增水还是减水。4. 加工过程中的四个高频坑含排查链路4.1 单位与符号cm、mm、slev到底是不是等效水高先说一个最容易踩的坑单位。CSR mascon不同文件版本、不同变量的单位并不统一。有的变量叫slev注释写的是“equivalent water thickness”单位是cm有的文件里的变量是lwe_thickness单位却是mm。如果你不打印ds直接开算趋势很容易差10倍。我排查过的一个案例同事下载了一批数据算出来的流域趋势是每年几百毫米明显不合理。我让他先打印变量属性结果发现他把原始单位mm误当成了cm处理导致趋势偏大10倍。建议是拿到数据后第一步永远先执行print(ds)并写一小段断言来检查变量单位unit ds[slev].attrs.get(units, ) assert cm in unit or mm in unit, fUnknown unit: {unit}符号问题也要注意。slev通常代表的是地表质量负载对应的等效水高变化符号约定一般是正值表示水储量增加。但某些辅助文件可能是“重力位变化”或“等效水高变化”的反号。用前先画一张全球图看看青藏高原、亚马逊等已知水储量变化区域的正负是否符合常识。4.2 参考时段的剩余信号不要对长期平均再求一次平均CSR mascon产品的数值是相对某个参考时段的变化量。RL06版本的参考时段一般是2004年1月到2010年12月也就是这一段的长期平均被设为0。这是官方明确说明的。但这带来一个常见错误有人下载数据后为了得到“相对多年平均的异常”会再对2004到2010或者全时段再求一次平均然后扣除。如果参考时段恰好和你的研究时段重叠这个二次平均会把真实的长期趋势部分抹掉。正确的做法是直接使用产品已有的数值不再额外扣除全时段平均。如果你确实需要相对某一特定基准的异常比如相对2010年以后的平均值那可以自己在指定时段上扣平均但要明白这会把趋势的“原点”移动不是消除趋势。我之前在做一个流域分析时全时段的平均值并不为0检查后发现是GRACE和GRACE-FO的衔接期数据导致了若干月份的跳变。这种跳变不是参考时段的问题而是两个任务之间的仪器差异做长序列分析时需要留意2017年底到2018年初的衔接。4.3 GIA是否已经扣除不要再扣一遍冰川均衡调整GIAGlacial Isostatic Adjustment是固体地球对末次冰期以来地表负载变化的缓慢响应。在GRACE重力信号里GIA会造成很大的长波趋势比如在格陵兰和南极地区GIA趋势占信号的比例非常高。不同数据处理机构的产品处理方式不同。CSR mascon产品在官方说明里写明了推荐使用的GIA模型并且发布的版本通常已经包含了GIA校正。也就是说用户拿到的数据里已经减掉了GIA效应不需要再额外扣除一次。我见过从事极地研究的同学拿到CSR mascon后又用ICE-6G模型扣了一次GIA导致最终趋势严重偏低。排查链路是这样的先看数据产品README确认是否已扣除GIA再看自己研究区域是否在强GIA区最后对比同一站点附近的GPS垂直速度看看趋势数量级是否合理。如果你要对比不同机构的mascon产品务必确认各自的GIA处理是否一致。CSR、JPL、GSFC对GIA模型的选择不同直接对比会导致几毫米每年的趋势差异。4.4 尺度因子的误用CSR、JPL、GSFC三家的差异凡是做过GRACE数据的人一定听说过“尺度因子”scale factor。这个概念源自JPL mascon产品由于数据处理中的约束和滤波真实信号会被衰减需要用尺度因子乘回原信号。JPL会根据每个网格的时间变化幅度计算一个最优尺度因子并随数据一起发布。CSR mascon产品的情况不太一样。CSR采用的球冠基函数和约束方式在设计时尽量减少了空间平滑造成的信号衰减因此官网通常不要求用户额外乘以尺度因子。但这不等于CSR没有尺度信息。在处理某些特定区域时球冠平滑仍然会带来一定程度的信号泄漏。这里的坑在于把JPL的尺度因子处理习惯直接套用到CSR产品上会给信号乘上一个不合理的增益。我第一次处理时顺手把JPL的scale因子乘到CSR数据上结果青藏高原好几个网格的趋势大得离谱后来仔细排查才发现产品之间的处理逻辑根本不一样。正确的做法是使用CSR产品时直接按官方提供的等效水高数值使用如果做跨产品对比先把各家的尺度因子处理方式统一比如都先把数据还原到“未加尺度因子”状态再一起比较。5. 把加工结果做成“数据集”命名、元数据与质量检查5.1 数据集的目录规范与文件命名项目做大以后加工出来的数据如果不做规范化管理一个月后再看就会乱成一锅粥。我后来把加工产物整理成标准数据集目录结构如下csr_mascon_basin_dataset/ ├── data/ │ ├── raw/ # 原始CSR mascon文件 │ ├── processed/ # 裁剪、掩膜后的区域数据 │ └── final/ # 最终时间序列和趋势结果 ├── shapefiles/ # 流域边界矢量 ├── scripts/ # 处理脚本 ├── metadata/ │ ├── dataset_metadata.json │ └── qc_report.html └── README.md文件命名我统一采用{机构}_{版本}_{区域}_{变量}_{时间范围}.nc的风格。例如CSR_RL06_Yangtze_ewh_200204_202308.nc。这样从文件名就能看出数据来源、版本、空间范围、变量和时间跨度比final_v3_final_2.nc不知道高明到哪里去了。如果是逐月数据我建议用basin_ewh_200204.nc这种按月存放的方式再配合目录下的时间索引文件后期追加新数据也方便。5.2 元数据字段设计数据文件里如果只存数值没有元数据后续协作或者一年后自己回头看时很多物理含义都会丢失。我一般在NetCDF文件里顺手写入以下属性ds_out.attrs[title] CSR RL06 mascon equivalent water height over Yangtze Basin ds_out.attrs[institution] Your Lab ds_out.attrs[data_source] CSR Mascon RL06 v02 ds_out.attrs[reference_period] 2004-01 to 2010-12 ds_out.attrs[gia_correction] ICE-6G-D applied by data provider ds_out.attrs[scale_factor_applied] None, per CSR official guidance ds_out.attrs[processing_date] 2025-01-15 ds_out.attrs[contact] your_emailexample.com如果要做成JSON元数据文件也是同样的字段结构。这些信息对论文方法部分非常有用写的时候直接照抄就行不用再回忆当时是怎么处理的。5.3 质量检查清单异常值、缺失值、变量边界数据集发布前我跑一套简单的质量检查脚本包括这几项缺失月份统计GRACE和GRACE-FO都有若干缺失月份比如2017年下半年的GRACE数据并不完整。把这些缺失月份明确列出来。空间范围检查掩膜后区域内不应有大量NaN除非网格本身在海洋或无数据区。数值合理性全球等效水高的年振幅通常在几十厘米以内如果出现绝对值超过3米的网格多半是掩膜没做好或单位搞错了。趋势显著性对每个网格做最小二乘趋势估计如果趋势的置信区间跨越0在出图时可以降低透明度免得读者误读。时间连续性检查时间轴是否等间隔是否存在重复时间。写一个简单的QA报告把每个网格的趋势、季节振幅、缺失比例都画成图能快速发现加工中的bug。6. 一个完整案例长江流域水储量变化分析6.1 读取与裁剪为了演示效果我拿长江流域作为案例。脚本流程是读取所有逐月NetCDF文件按经纬度裁剪出一个大的矩形范围比如经度90到122纬度21到36然后叠加长江流域矢量掩膜。裁剪代码示例import glob import xarray as xr import geopandas as gpd import numpy as np from rasterio.features import rasterize # 1. 读取全部逐月文件并拼接 files sorted(glob.glob(data/raw/CSRM_BA01_*.nc)) ds_all xr.open_mfdataset(files, concat_dimtime, combinenested) ds_all ds_all.sortby(time) # 2. 构建流域掩膜 lat ds_all[lat].values lon ds_all[lon].values basin gpd.read_file(shapefiles/yangtze_basin.shp).to_crs(EPSG:4326) # 利用经纬度网格中心点生成0/1数组 shapes [(geom, 1) for geom in basin.geometry] mask rasterize( shapes, out_shape(len(lat), len(lon)), transformfrom_origin(lon.min(), lat.max(), abs(lon[1]-lon[0]), abs(lat[1]-lat[0])), fill0, dtypenp.uint8 ) # 3. 应用掩膜 ds_region ds_all.where(mask 1)这个处理里有一个关键点rasterize的transform参数必须严格对应卫星数据的经纬度网格起点和分辨率否则掩膜会整体偏移把流域外的网格也选进来。我实际工作中曾因为from_origin的经度起点不精确导致半个掩膜漂移到流域外画出来的时间序列趋势直接失真。6.2 时间序列计算与趋势提取掩膜裁剪后按空间维度做面积加权平均得到流域逐月平均的等效水高序列# 如果网格面积随纬度变化需要按cos(lat)加权 weights np.cos(np.deg2rad(lat)) weights_2d np.broadcast_to(weights[:, None], (len(lat), len(lon))) weights_2d np.where(mask 1, weights_2d, 0) # 对每个时间步做加权平均 values ds_region[slev].values valid ~np.isnan(values) basin_series np.full(values.shape[0], np.nan) for t in range(values.shape[0]): valid_t valid[t] (mask 1) if valid_t.sum() 10: basin_series[t] np.average(values[t][valid_t], weightsweights_2d[valid_t])再做趋势拟合就很简单了沿用第3节里的最小二乘模型。长江流域GRACE时间序列的特点是季节变化非常明显夏季水储量高、冬季低但常年趋势相对较小。做趋势分析时建议以年为单位输出趋势并给出95%置信区间。6.3 可视化和导出最终图表我一般画成上下两部分上面是流域平均等效水高的逐月曲线下面是用颜色表示趋势的空间分布图。空间图可以用cartopy绘制叠加流域边界和主要水系。导出数据时最终数据集用NetCDF格式保存时间序列再额外导出一份CSVdf pd.DataFrame({ time: ds_all[time].values, basin_mean_ewh_cm: basin_series }) df.to_csv(output/yangtze_basin_ewh_monthly.csv, indexFalse)NetCDF导出时别忘了把前面提到的元数据写进去。6.4 验证思路与常见误区加工完的数据一定要做一次外部验证不能画完图就收工。我常用的验证方式有三个一是和卫星水文模型对比。比如把CSR mascon的流域平均水储量变化和GLDAS水文模型、PCR-GLOBWB模型对比看季节相位和年际波动是否一致。如果mascon趋势和模型趋势符号相反先别着急说模型不准回头检查自己的时间序列是否扣错了参考时段。二是和GRACE官方发布的同区域时间序列对比。CSR官网提供了一些流域或全球平均的时间序列可以直接拿来对比数量级应该非常接近。三是和局部地面监测数据对比。比如大型水库的蓄水量变化、地下水井水位变化对比时要注意空间尺度差异地面点数据往往不能完全代表整个流域的平均状态。常见误区主要是三类把0.25度网格当真实分辨率误以为能识别几十公里尺度的局部信号直接用全球文件做流域平均而不做面积加权导致高纬度网格贡献偏大以及把GRACE和GRACE-FO期间的数据直接拼起来做趋势没有处理任务间的系统差造成2018年前后出现人为跳变。数据加工这个环节说起来不复杂绝大部分时间都耗在“确认单位”、“确认参考时段”、“确认GIA处理”这类的检查上。我个人体会是把这些琐碎的确认做成必须执行的检查函数写进处理流程的最前面后面再大的数据量都能跑得安心。希望这篇记录能帮你少走几步弯路后面有时间我再单独写写CSR、JPL、GSFC三家mascon产品的对比和实际差异。本文还有配套的精品资源点击获取