GRACE卫星重力数据处理：CSR mascon产品原理与Python实践

简介卫星重力观测技术通过监测地球质量分布的时间变化为研究全球水循环与气候变化提供了独特视角。GRACE及其后继GRACE-FO卫星任务利用星间测距精确反演地球重力场从而捕捉陆地水储量、冰盖消融及地下水变化的时空信息。其中mascon质量集中产品凭借局部约束反演策略显著抑制了传统球谐系数产品中的条带噪声和信号泄漏问题成为陆地水文学与冰川学研究的利器。CSR mascon作为主流产品之一以规则的1度网格输出等效水高异常附带尺度因子、空间掩膜和误差估计便于用户直接用于区域水量变化、极端气候事件响应及海冰质量平衡等分析。本文结合Python与xarray工具系统讲解CSR mascon数据的读取、预处理、区域加权平均、趋势分解等完整流程并针对常见错误提供排查方法助力研究者快速掌握卫星重力数据处理的核心技能。1. 项目背景为什么需要专门的mascon数据产品1.1 GRACE/GRACE-FO卫星与重力场反演简介提到CSR mascon数据必须先聊聊GRACE和GRACE-FO这两颗对其实是两个卫星编队改变了整个地球物理研究范式的卫星任务。GRACEGravity Recovery and Climate Experiment从2002年运行到2017年GRACE-FOFollow-On从2018年接棒至今它们做的事情本质上很简单——通过测量两颗卫星之间微波测距的微小变化反演出地球重力场的时空变化。听起来很抽象但落到实际应用上就是能告诉我们某个流域的地下水到底少了多少、格陵兰冰盖每年融化有多快、大地震之后断层附近的地壳是如何调整的。这里有个关键点重力卫星并不直接测量水它测量的是地球质量的空间分布变化。在陆地水文和冰川学领域我们把这个信号解释为等效水高Equivalent Water HeightEWH单位通常用厘米或毫米水柱表示。之所以能这么做是因为在月到十年的尺度上地球内部质量重新分布的信号远小于地表水、冰、雪等成分的变化。这在物理上是合理的近似也是整个mascon数据应用的基础。1.2 球谐系数产品与mascon产品的核心区别我们这一行的人都知道传统的GRACE数据处理会走球谐系数这条路先把卫星轨道和星间测距数据解算成一组全球球谐系数也就是Stokes系数再通过合成公式计算地表质量变化。这个过程看起来很标准但实际使用中会被两个问题折磨到怀疑人生。第一个问题是条带误差stripe noise。球谐解算在南北方向会产生细长的条带误差尤其是高阶项必须用滤波来压制。经典的方案是destriping比如Swenson和Wahr在2006年提出的P3M6滤波或者用高斯平滑半径300公里甚至更大。可一滤波真实信号也被抹掉了——空间分辨率从理论上几百公里直接退化到上千公里。第二个问题是泄漏误差leakage error。球谐展开本质上是用全球基函数拟合区域性信号结果就是某区域的信号会被泄漏到周边包括从陆地泄漏到海洋、从强信号区泄漏到弱信号区。你画一个流域边界去积分质量变化算出来的结果总是比真实值偏小因为一部分信号被平滑滤波搬走了。mascon方法走的是另一条路它不做全球球谐拟合而是在空间上把地球表面划分成一个个质量集中块mass concentration即mascon每个块内部假设质量均匀分布然后直接通过卫星观测数据解算每个块的质量变化。这样做的好处很直接条带误差几乎不存在了因为基函数是局部的泄漏误差也被大幅抑制因为在反演迭代过程中可以施加空间约束比如陆地和海洋分开处理、限制海洋信号向陆地泄漏等。业界主流的mascon产品有三大系NASA JPL的JPL mascon、NASA GSFC的mascon也就是GSFC mascon以及UTCSR的CSR mascon。三者在反演策略、约束方式和后处理上有细微差别但整体思路一致。这几种数据各有拥趸有人喜欢JPL的平滑特性有人觉得GSFC在极区表现更好我自己则因为CSR产品的开源程度高、附带误差估计比较完善加上格式统一、易上手日常处理基本是用CSR居多。今天我重点拆解的就是CSR mascon这条线。1.3 CSR mascon数据产品的独特设计CSR mascon产品由德州大学奥斯汀分校空间研究中心Center for Space Research简称CSR发布目前主流的版本是RL06基于GRACE和GRACE-FO统一重新处理。它有一个很实用的特点分辨率虽然是规则的1度网格采样但实际反演不是等网格的而是在高纬度和海岸线附近做了加密处理。这带来一个好处就是在极地冰盖和沿海区域你可以拿到比纯粹1度球谐产品更细腻的空间分布。另一个让我觉得踏实的地方是CSR mascon产品附带了一个spatial mask文件明确区分了陆地、海洋、冰盖、湖泊等掩膜信息。因为反演过程中本身就把陆地和海洋分开处理过拿到数据后做区域分析时不需要自己再去手动扣海陆边界错误概率低了不少。要知道GRACE信号在海洋上主要是海水质量变化而在陆地上主要是水量变化两者物理意义不同如果掩膜搞错了后面的分析基本就没法看。此外CSR mascon自带不确定性估计。这个误差并不只是拟合残差的统计误差还包括了轨道误差、仪器误差以及反演约束导致的误差传播。虽然这个不确定性估计偏乐观实际表现往往比它报的数字略大一些但至少给了一个量化的参考框架。对写论文的人来说有误差条总比没有好。2. 数据获取与文件格式完全解读2.1 官方数据源与版本选择先解决去哪儿下载的问题。CSR mascon RL06 v02的数据可以从UTCSR官网直接下载或者通过NASA的Physical Oceanography Distributed Active Archive CenterPO.DAAC访问。如果你是学生或者刚入行建议直接走官网的ftp或https目录因为目录结构很直观不需要额外申请权限。版本方面目前CSR mascon已经更新到了RL06 v02/v03/v04等版本。如果不确定选哪个我的建议是优先选最新版本因为新版本通常会修复早期算法中约束参数、海冰掩膜或者陆地掩膜的若干问题。但有一点要提醒如果后续做结果对比最好统一用同一个版本不要混着用因为不同版本之间的绝对数值可能有几个厘米等效水高的差异这在某些水文研究中足以影响结论。下载时注意区分月度产品和季度产品的文件命名规则。CSR mascon的文件名一般长这样CSR_Mascon_RL06_V02_yyyy_mm.nc其中yyyy_mm对应数据月份。文件是netCDF格式直接用Python的xarray或netCDF4库就能读取。一个典型文件体积不大大概几十MB包含全球网格的月度质量异常场不需要做复杂的预处理就能拿来画图或做时间序列提取。2.2 NetCDF文件内部结构详解拿到.nc文件之后第一件事就是用 xarray 打开它看看结构。我习惯先用一行命令快速浏览文件内容和维度信息import xarray as xr ds xr.open_dataset(CSR_Mascon_RL06_V02_2022_01.nc) print(ds)输出内容大致是这样不同小版本命名会略有差别Dimensions: time: 1 lat: 180 lon: 360 Coordinates: * time: (time) datetime64[ns] 2022-01-01 * lat: (lat) float64 -89.5 ... 89.5 * lon: (lon) float64 -179.5 ... 179.5 Data variables: solution: (time, lat, lon) float64 uncertainty: (time, lat, lon) float64 spatial_mask: (lat, lon) int16 scl: (lat, lon) float64这里面有几个变量要重点说明solution就是每个网格点的等效水高异常值单位是厘米水柱cm w.e.也就是cm of water equivalent。注意是厘米不是毫米很多新手在这里单位换算出错。uncertainty对应网格点的误差估计单位同样是厘米水柱。spatial_mask整数值掩膜不同数字代表不同地表类型比如0表示海洋、1表示陆地、2表示冰盖、3表示湖泊等。scl尺度因子scaling factor用来恢复反演过程中被约束压制掉的信号。如果要做区域总量分析需要把solution乘以scl得到恢复后的等效水高值再参与积分。不乘的话某些信号偏弱的区域会低估。维度lat和lon是规则的1度网格从-89.5到89.5和-179.5到179.5。注意这里不是经纬度从0开始而是以格点中心表示的。所以当你用经纬度范围选取区域时要留出半个格点的余量否则边界上的格点会被切掉。2.3 核心参数与单位换算在开始处理数据之前先统一单位、统一概念这一步能少踩不少坑。首先是solution的单位。CSR RL06 v02的官方文档一般写的是cm w.e.但这并不意味着全球所有网格的物理含义都相同。在陆地上它代表的是陆地水储量异常土壤水、地下水、雪、地表水等总量的变化在海洋上它代表的是海水质量的异常变化在冰盖上则代表冰质量的累积或消融。看起来都是等效水高但物理含义差别很大做解释的时候必须分区域讨论。其次是尺度因子scl。有些教程会提醒你分析流域时只需要把solution乘以scl再做区域平均。为什么需要尺度因子因为mascon反演本身是有约束的约束会把真实信号往零均值方向拉拢导致振幅被低估。尺度因子就是用来补偿这个低估的它可以通过在约束反演时用合成的先验信号来标定。实际使用中scl的取值范围大概在1到3之间在信号弱的地方会更大一些。再次是掩膜spatial_mask。这个字段在整个产品里属于小细节、大作用的类型。比如你想算某个陆地流域的质量变化直接用掩膜值筛掉海洋网格就行。不需要再去额外下载一个海岸线形状文件——省事也降低出错概率。最后说时间坐标。CSR mascon的时间是月中的时间戳直接用datetime格式给出的读取后可以直接用于时间序列的画图。但要注意GRACE数据中间有若干个月缺失比如2011年、2016年、2017年出现过数据断层GRACE和GRACE-FO交接之间也有空白期。处理时间序列时不需要手动插值填补空洞直接画散点图即可后期如果需要算趋势可以考虑用时间回归模型来规避缺失月的影响。3. Python处理流程的完整实现3.1 环境搭建与依赖库处理CSR mascon数据强烈推荐Python生态尤其是xarray、numpy、matplotlib或者cartopy这三个组合。其中xarray是核心它能直接操作带坐标标签的netCDF数据比手动操作数组方便太多。建议在conda环境里一次性装好conda create -n grace python3.10 conda activate grace conda install -c conda-forge xarray netcdf4 numpy scipy matplotlib cartopy装好之后下面这套流程就可以直接跑了。3.2 读取CSR mascon数据的标准姿势一次性读取多个月的CSR mascon文件并合并成一个xarray Dataset是我最常用的方式。这样可以避免循环里反复打开文件处理起来也快得多。import xarray as xr import glob file_list sorted(glob.glob(CSR_Mascon_RL06_V02_*.nc)) ds xr.open_mfdataset(file_list, combineby_coords)这样读取之后ds[solution]的维度就是(time, lat, lon)可以直接做区域切片和统计。如果你的机器内存比较小也可以先用ds.load()检查一下内存占用不行的话分月份循环处理。我实测过CSR mascon全球的月度场一个月是180×360个格点几十年的数据叠在一起体量也就几百MB到1GB左右。现代PC跑起来完全没问题不需要分布式那套东西。接下来是做基本的数据预处理乘尺度因子、选陆地或海洋掩膜。solution ds[solution] * ds[scl] mask ds[spatial_mask] # 只看陆地 solution_land solution.where(mask 1)因为scl本身只有空间维度没有时间维度xarray的广播机制会自动在时间维度上做乘法。这一步处理完之后solution_land里面陆地区域是数值海洋区域是NaN后续做区域统计时会自动跳过NaN。3.3 区域平均时间序列提取流域尺度的水量变化估计是mascon数据最常用的场景之一。操作要点是先根据经纬度范围切出目标区域再做加权平均。这里有一个容易忽略的问题在1度规则的经纬度网格上纬度越高的格点实际覆盖面积越小。如果你算区域平均时不考虑这个效应高纬度的格点会被过度加权结果偏得离谱。解决方案是给每个格点乘以cos(lat)权重。我自己写了一个简单函数可以设置经纬度范围然后计算区域加权平均。顺手把多个区域放在一起遍历处理。import numpy as np def regional_series(lat_range, lon_range, solution_da): lat_sel slice(lat_range[0], lat_range[1]) lon_sel slice(lon_range[0], lon_range[1]) sub solution_da.sel(latlat_sel, lonlon_sel) weights np.cos(np.deg2rad(sub[lat])) # 网格面积权重 weight_sum weights.sum() series (sub * weights).sum(dim[lat, lon]) / weight_sum return series然后把区域范围内的数据代入即可。举个例子如果研究华北平原的地下水变化范围大致取(30, 42)度纬度和(110, 122)度经度。调用函数后会得到一个以时间为索引的一维数组然后就可以顺手画图或者做趋势拟合。3.4 趋势与季节信号提取拿到了时间序列下一步通常要分解出趋势项、季节项和残差。我推荐用一个简单的多元线性回归模型自变量包括时间线性趋势、年周期sin/cos和半年周期因变量就是区域的平均等效水高。import pandas as pd from numpy.polynomial import polynomial as P def decompose_series(time, series): # time: numpy datetime64 array # series: 等效水高序列 t (time - time.min()).astype(timedelta64[M]).astype(float) / 12.0 X np.column_stack([ np.ones_like(t), t, np.sin(2 * np.pi * t), np.cos(2 * np.pi * t), np.sin(4 * np.pi * t), np.cos(4 * np.pi * t), ]) coeff, _, _, _ np.linalg.lstsq(X, series, rcondNone) trend coeff[1] * 12 # 厘米/年 seasonal X[:, 2:4] coeff[2:4] residual series - X coeff return trend, seasonal, residual这个模型虽然简单但在处理月尺度水文信号时非常实用。趋势项trend就是区域等效水高的长期变化速率单位是cm/year。如果你想换算成水量体积变化再乘以区域面积就行注意单位统一。这里有个经验之谈算趋势的时候别上来就对原始数据做线性拟合。水文信号里季节性波动振幅常常比长期趋势大一个数量级如果不把季节项去掉拟合误差会非常大趋势估计的置信区间也很宽。先用多元回归把季节项剥离剩下的信号再做趋势估计就稳很多。3.5 常见可视化绘制处理完数据画图是必须的。用cartopy画全球或区域的等效水高分布图配合pcolormesh是最快的。下面是一个直接能跑的模板可以做某个月份的全球陆地水储量异常图。import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature fig plt.figure(figsize(12, 6)) ax fig.add_subplot(1, 1, 1, projectionccrs.PlateCarree()) data_plot solution.sel(time2022-01-01) data_plot.plot.pcolormesh( axax, transformccrs.PlateCarree(), cmapRdBu, vmin-20, vmax20, cbar_kwargs{label: Equivalent Water Height (cm w.e.), shrink: 0.8}, ) ax.add_feature(cfeature.COASTLINE, linewidth0.5) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.5) plt.title(CSR Mascon RL06 EWH Anomaly (2022-01)) plt.tight_layout() plt.show()如果要画时间序列图直接matplotlib就行趋势线、季节线和原始序列画在一起加个图例即可。这个视觉组合我觉得对读者来说是最友好的信息密度也够。4. 实际应用场景与结果解读4.1 区域水量变化监测CSR mascon最常见的应用是做流域尺度的陆地水储量变化监测。这种监测的独到之处在于地表实测井网往往分布不均只覆盖有限点位而卫星重力数据天然是空间平均的能够给出整个流域的大尺度感观。举个具体的例子假如关心长江源区的水储量变化区域范围大概取东经90°到100°、北纬32°到36°。我用上一节的函数提取序列然后做趋势分解能看到明显的降水和冰川消融信号混合在一起。如果要进一步分离一般会结合降水数据、冰川物质平衡观测和积雪覆盖数据分析——不过这就是另一个故事了。有一点要注意CSR mascon在山区因为反演约束本身更保守数值可能比预期偏小所以在解释绝对值时要谨慎看相对变化和空间格局更可靠。4.2 极端气候事件响应mascon数据另一个被高频使用的场景是分析极端气候事件比如干旱和洪水对区域总水储量的冲击。2011年德州干旱、2019年澳大利亚东南部干旱、2020年中国南方汛期都能在GRACE/GRACE-FO数据里看到明显的水储量异常信号。处理方法很简单先把多年每个月的信号做逐月气候态然后拿目标月份减去气候态得到的距平场就能直观展现出这次事件的强度和空间范围。CSR mascon的优势是在部分区域信号恢复得好距平场更加锐利不会像球谐产品那样糊成一片。当然GRACE的时间分辨率是月级的无法用来监测日尺度的洪水过程它适合的是事后评估和趋势分析。如果有研究人员拿它来做实时洪水预警这个思路本身就不太现实需要和陆面模型、遥感数据配合。4.3 结果验证方法用mascon数据做研究免不了要和别的数据互相验证。我常用的验证思路有两类第一类是内部自洽验证。CSR mascon、JPL mascon、GSFC mascon三者做同一个区域的时间序列如果三者的趋势和季节相位一致那说明信号是真实的如果某一个产品单独异常就要仔细检查约束参数带来的影响。从实用角度讲论文里同时展示三个产品结果是个不错的加分项也能提升结论的可信度。第二类是外部数据验证。对流域来说可以拿GRACE结果和陆面模型GLDAS、ERA5-Land对比对地下水区域可以和井水位观测对比注意单位换算井水位是深度单位质量变化需要用水位变化乘以给水度近似对冰盖区域可以和高程变化如ICESat-2测高对比冰厚度变化乘以冰密度近似后应和mascon总质量变化一致。这些交叉验证过程并不能说明哪一个产品绝对正确但能帮助我们判断某个信号是否可靠、是否需要做后处理修正。做研究方法论上的严谨比一张好看的图重要得多。5. 常见问题与排查技巧实录5.1 数据读取与维度不匹配现象用xr.open_mfdataset读取多个文件时报维度不一致的错误。原因不同月份的CSR文件有时因为格式版本微调会导致某一维度命名或者坐标值出现细微差异比如lat坐标从-89.5变成了-89.0导致全部月份无法合并。排查方法先单独打开一个文件检查坐标ds_one xr.open_dataset(CSR_Mascon_RL06_V02_2018_01.nc) print(ds_one.coords)如果确实存在坐标不一致要么重新下载统一版本的文件要么先用循环读取再把solution字段抽出来拼成numpy数组再用xr.DataArray重新封装。在我的经验里重新下载统一版本是最干净的办法不建议临时去对齐坐标。5.2 区域平均面积加权错误现象算出来的区域平均值在高纬度一直明显偏大或偏小。原因没有考虑cos(lat)面积权重导致等面积平均在高纬度格点权重过大。解决办法在上面给出的regional_series函数里已经包含面积加权。如果自己写代码一定要注意 lat 单位必须从角度转为弧度再取cos。5.3 符号正负与物理意义混乱现象某地区的等效水高趋势为负但实际地面观测显示地下水位在上升。原因分析大多数情况下不是数据错了而是符号约定和参考基准的问题。CSR mascon的solution场默认参考的是2004到2009年的平均场所以所有值都是相对这个基期的异常。如果地面观测的基期不同两者就不能直接比较。此外GRACE观测的是总水柱变化地下水只是其中一个分量土壤水和地表水的反向变化会掩盖地下水信号这也是常见误解。对策处理前先确定自己的研究基期并确认产品的参考期在论文方法部分明确写出来。否则审稿人一眼就能看出基期不匹配的问题。5.4 尺度因子未使用的典型错误现象区域平均序列画出来后振幅特别小和文献差了一个量级。原因没有乘以scl尺度因子。这也是CSR mascon使用中最常见的一个问题——很多用户下载数据之后直接就用solution忽略了尺度因子。对策无论做全球还是区域分析开头先做一步solution * scl把它作为一个默认操作固定下来。如果想严谨地验证是否用了尺度因子可以检查数据量级比如华北平原地下水储量异常的月振幅通常在5到10厘米水柱左右如果算出来只有2厘米以下多半是没乘因子。5.5 冰盖与海岸线区域信号泄漏现象在格陵兰或南极沿岸的网格出现相邻格点一正一负的椒盐图案。原因mascon反演在海岸线附近由于掩膜处理硬性分离了陆地和海洋沿岸格的噪声会增加。这是反演策略的固有特点所有mascon产品都有这个问题。对策在做极地区域时空分析时尽量用spatial_mask过滤掉边缘网格或者对沿岸做一圈缓冲区剔除。同时建议和GSFC mascon或JPL mascon对比如果信号只在某个产品里出现大概率是噪声而非真实信号。6. 后续扩展思路与个人建议我在处理CSR mascon数据的过程中最大的感受是这套产品在易用性和物理合理性之间取得了很好的平衡。对刚接触卫星重力数据的人来说它比球谐系数产品友好得多对资深用户来说它的误差估计和掩膜设计又足够支持深入分析。如果你要做长时间序列的陆地水储量变化、干旱监测或冰川质量平衡CSR mascon是一个很稳妥的起步选择。最后分享一个小技巧如果你准备把mascon结果写进论文建议在方法部分明确写出产品版本、下载时间、参考期、是否乘了尺度因子、是否做了区域面积加权这五要素。审稿人对GRACE数据的处理透明性要求很高这五要素写清楚很多数据质疑根本不会出现。个人实际使用中这五要素全部写清后论文材料审查阶段基本没人再纠结数据问题。这个内容后续还可以从几个方向扩展一是做多产品交叉验证把CSR、JPL、GSFC同时纳入不确定度分析二是结合机器学习或数据同化方法把mascon信号降尺度到流域内的高分辨率网格三是把处理流程封装成自动化的Python包方便团队内部共享。无论哪条路CSR mascon数据集都是一个值得深耕的起点。本文还有配套的精品资源点击获取