CSR mascon数据处理实战：从GRACE卫星重力数据到区域水储量时间序列

简介卫星重力技术为地球水循环研究提供了独特视角GRACE及GRACE-FO任务通过双星测距原理捕捉全球重力场变化从而反演陆地水储量迁移。其中等效水高EWH是衡量储水量变化的核心指标。实际应用中传统球谐系数方法存在条带噪声与平滑泄漏问题而mascon质量集中产品直接提供网格化的等效水高时间序列降低了处理门槛。CSR mascon作为代表性产品融合GRACE/GRACE-FO统一解算适合区域水文与气候研究。本文从数据处理视角出发介绍基于Python与xarray的NetCDF文件读取、纬度加权区域平均、趋势与季节项提取等关键步骤并讨论可视化验证与常见坑点助力构建可靠的水储量时间序列。 搞地球物理和水文研究的人大概率绕不开GRACE卫星重力数据。很长一段时间里处理GRACE数据的传统流程相当折腾先下载球谐系数选截断阶数做高斯平滑再想办法去除条带误差一整套流程下来新手基本要磨掉一两个星期。后来mascon产品逐渐成熟CSR mascon数据成了很多人直接上手的首选——打开就是网格化的等效水高时间序列信号干净处理链路大幅缩短。这篇文章就围绕CSR mascon数据加工数据集这件事从原理到实操讲清楚它是什么、怎么拿到、怎么加工成能用的区域水储量时间序列再把这些年踩过的坑也一并填上。内容适合刚入门卫星重力数据处理的研究生也适合想用GRACE/GRACE-FO数据做区域水文分析但不想跟球谐函数死磕的同行。1. 先把概念讲清楚mascon到底是什么CSR这版凭什么值得用1.1 GRACE卫星是怎么“称”地球重量的GRACE任务的基本原理通俗点说就是两颗卫星“互相感应”。GRACE和GRACE-FO都是双星编队前后两颗卫星相距约220公里当飞过某个质量异常区域时前面的卫星先感受到引力变化、速度发生微小改变后面的卫星再跟着变这个时间差导致两星之间的距离产生微米级的变化。卫星上搭载的微波测距系统GRACE或激光干涉测距系统GRACE-FO把这个距离变化记录下来通过长时间连续观测就能反演全球重力场的时空变化。由于水的密度低、在地球表面迁移量最大在月到年这样的时间尺度上地球重力场短期变化的主要贡献者就是水——包括地下水、土壤水、冰盖冰川和海水。所以GRACE数据最核心的应用场景就是水循环研究某片流域的地下水在减少格陵兰的冰在融化亚马逊河流域的储水量在每个雨季上升多少这些变化都能在重力信号里体现出来。1.2 mascon方法与球谐方法的核心差异传统球谐方法是在全球范围内把重力场用球谐函数展开相当于用一组无限长的“正弦波”去拟合全球重力场的空间分布。它的优点是数学上完备、漂亮但问题也很突出高阶项噪声大需要截断球谐展开的截断误差会表现为南北方向的条带噪声业内俗称“条纹”为了压制这些噪声通常要做高斯平滑但平滑在去噪的同时也把真实的空间分辨率削掉了一截还会造成信号从一个区域泄漏到另一个区域。mascon是mass concentration的缩写思路完全换了一个方向把地球表面离散成一个个小区域每个小区域就是一个质量块直接在块内求解质量变化率。这种方法本质上是从“全局拟合”变成“局部拟合约束”通过空间约束和时间正则化等手段能显著压制噪声和条带误差。它的好处非常直观——你要研究亚马逊流域就直接给一片区域区域内的信号就是被约束出来的信号泄漏小边界清晰。1.3 CSR mascon产品的定位与特色目前全球用得比较多的mascon产品有三家CSR德克萨斯大学奥斯汀分校空间研究中心、JPLNASA喷气推进实验室、GSFCNASA戈达德太空飞行中心。三家思路有差异结果自然也有细微不同。CSR mascon在区域水文信号恢复方面表现比较均衡约束相对温和不会把真实的强信号削得太厉害而且产品是直接网格化的NetCDF文件用户拿来就能用不需要再自己滤波。另外要提一下CSR mascon是GRACE和GRACE-FO统一解算的也就是说2002年到2017年的GRACE数据和2018年之后的GRACE-FO数据在同一个产品框架下发布时间序列接得上不需要自己做两个时间段的拼接。这一点在后期处理中非常省心。2. 数据获取与文件结构先搞懂你手里的东西2.1 从哪里下载、下载哪个版本CSR mascon数据可以在德州大学空间研究中心官网下载也能在NASA的PO.DAAC数据中心找到。文件名一般是CSR_GRACE_GRACE-FO_RL06_Mascon_v02.nc这样的格式关键词含义是RL06代表第六版数据释放v02代表第二个正式发布版本。下载时建议优先选RL06版本原因是它基于更新的背景重力场模型和更精确的轨道数据比RL05精度有明显提升。另外要注意看数据的起止时间GRACE和GRACE-FO之间有一个约11个月的断档2017年中到2018年中这是正常的不是文件损坏。2.2 NetCDF文件里的关键变量打开NetCDF文件后核心变量就几个但每个都要搞清楚含义否则后面全都白算lat纬度单位degree_northlon经度单位degree_easttime时间变量CSR的产品通常以“days since 2002-01-01”为参考数值单位是天ewh等效水高Equivalent Water Height单位是cm这是核心数据正值代表该区域质量增加负值代表质量减少uncertainty或error数据不确定度如果文件里有建议一并读取用Python打开文件看一眼结构非常直观在下一节会给出完整代码。2.3 理解产品的时间分辨率和空间分辨率CSR mascon的空间网格一般是0.25度看起来分辨率不低但心里要有数GRACE卫星的物理分辨率大概是300公里左右0.25度网格只是输出格式的细化不代表你能在0.25度尺度上解读信号。很多人刚开始用的时候会被网格精度带偏以为能看到流域内某个小区域的变化这是不对的。时间分辨率是月尺度也就是每个月一个全球网格解。由于GRACE/GRACE-FO卫星的轨道设计和数据处理策略月解之间的噪声水平会有波动比如夏季北半球水量变化大、信号强冬季信号弱噪声占比相对高。处理时间序列时心里绷着这根弦后面做滤波和趋势提取时会清醒很多。3. 核心加工流程从原始文件到可用的区域时间序列3.1 环境准备与依赖库处理CSR mascon数据最顺手的工具组合是Python配xarray再加上numpy、pandas、scipy和matplotlib。xarray处理NetCDF文件的能力很强维度名和坐标信息直接保留做区域切片和加权平均都很方便。如果还没有安装用conda或者pip一行命令搞定conda install -c conda-forge xarray netcdf4 pandas scipy matplotlib建议使用conda创建独立环境避免依赖冲突。我自己的习惯是单独建一个grace_env环境专门跑重力数据相关的脚本不跟日常的Python环境混在一起。3.2 数据读取与基本检查读取CSR mascon文件的代码非常简洁import xarray as xr import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy import stats # 读取数据 ds xr.open_dataset(CSR_GRACE_GRACE-FO_RL06_Mascon_v02.nc) # 查看数据结构 print(ds)输出中能看到各个变量的维度、单位和数值范围。拿到文件后第一件事是检查数据范围lat是否从-90到90lon是0到360还是-180到180time的时间跨度是否覆盖你需要的时段。如果lon是0到360而你习惯用-180到180需要做一个坐标转换# 将经度从0~360转换为-180~180 ds ds.assign_coords(lon(((ds.lon 180) % 360) - 180)).sortby(lon)这个转换看似简单但漏掉会导致后面区域选择时切片范围对不上尤其是研究区域跨越本初子午线或180度线的时候坑很大。3.3 时间变量的处理CSR mascon的time变量通常是数值型单位是“天”参考日期是2002-01-01。需要转换成标准日期格式方便后续按年份和月份分析# 将时间转换为pandas的DatetimeIndex dates pd.to_datetime(ds[time].values, unitD, originpd.Timestamp(2002-01-01)) ds[time] dates转换之后打印ds[time]确认日期范围通常是从2002年4月左右到最新月份。时间序列中会有一个明显的断档GRACE与GRACE-FO之间的空隙这是正常的。3.4 区域平均提取目标流域的信号提取某个区域的水储量变化是mascon数据最常用的加工方式。这里以亚马逊流域为例经纬度范围大致是南纬20度到北纬10度、西经80度到西经50度。注意在等经纬网格上做区域平均时高纬度网格代表的实际面积比低纬度小所以必须做纬度余弦加权否则结果会被高维度区域带偏。# 选择区域 lat_min, lat_max -20, 10 lon_min, lon_max -80, -50 region ds.sel(latslice(lat_min, lat_max), lonslice(lon_min, lon_max)) # 纬度余弦权重 weights np.cos(np.deg2rad(region[lat])) weights weights.fillna(0) # 面积加权平均 # weights的形状是(lat,)需要广播到(lon,lat)用xarray的广播机制自动完成 regional_ewh (region[ewh] * weights).sum(dimlat).sum(dimlon) / (weights.sum() * len(region[lon])) # 转换为DataFrame方便后续分析 ts regional_ewh.to_dataframe(nameewh).drop(columns[lat, lon]) ts ts.dropna()加权平均这一步是区域时间序列提取的核心很多人会用简单的算术平均结果在高纬度地区会偏高这里一定要用余弦权重。3.5 趋势与季节项提取拿到区域时间序列后下一步通常是估计长期趋势和季节项。最直接的方法是线性回归加年周期、半年周期的正弦拟合# 准备自变量时间以年为单位 t_year (dates - dates[0]).days / 365.25 # 设计矩阵趋势 年周期 半年周期 X np.column_stack([ np.ones_like(t_year), t_year, np.sin(2 * np.pi * t_year), np.cos(2 * np.pi * t_year), np.sin(4 * np.pi * t_year), np.cos(4 * np.pi * t_year) ]) # 最小二乘拟合 coef, res, _, _ np.linalg.lstsq(X, ts[ewh].values, rcondNone) # 趋势项cm/year trend_cm_year coef[1] trend_mm_year trend_cm_year * 10 print(f长期趋势: {trend_mm_year:.2f} mm/year)拟合结果解读趋势项的单位是cm/year通常论文里用mm/year所以要乘以10。年周期项反映区域储水量的季节性波动幅度。亚马逊流域的典型特征就是明显的雨季/旱季周期年振幅可能在十几厘米甚至更高。如果只想看扣除季节项后的残差信号就用原始时间序列减去拟合的周期项。残差里往往隐藏着极端水文事件比如2010年亚马逊干旱造成的储水量异常下降。3.6 空间网格数据的进一步处理除了区域平均有时需要保留空间分布信息比如画某一年全球或区域的水储量变化图。这时可以直接对三维数据time, lat, lon做时间维度上的线性拟合# 对每个网格点做线性回归 t_years (ds[time].values - ds[time].values[0]).astype(timedelta64[D]).astype(float) / 365.25 # 用apply_ufunc批量计算每个格点的趋势 def linear_trend(y): mask ~np.isnan(y) if np.sum(mask) 12: return np.nan slope, intercept, r_value, p_value, std_err stats.linregress(t_years[mask], y[mask]) return slope trend_grid xr.apply_ufunc( linear_trend, ds[ewh], input_core_dims[[time]], vectorizeTrue, daskparallelized, output_dtypes[float] )然后直接trend_grid.plot()就能画全球趋势图。注意GRACE数据在某些区域如极地冰盖中心可能由于卫星轨道覆盖问题导致数据质量差异画图时留意异常值。4. 可视化与结果解读让数据“说话”的方式与验证4.1 时间序列图先看全貌再看细节区域平均后的时间序列图信息量非常大。建议把原始月值、3个月滑动平均、拟合趋势线画在同一张图上一眼就能看出趋势、季节项和异常事件的相对大小。fig, ax plt.subplots(figsize(12, 5)) ax.plot(ts.index, ts[ewh], o-, markersize3, labelMonthly, alpha0.7) ax.plot(ts.index, ts[ewh].rolling(3, centerTrue).mean(), g-, linewidth2, label3-month running mean) ax.plot(ts.index, fitted, r-, linewidth2, labelTrend seasonal fit) ax.axhline(0, colorgray, linestyle--, linewidth0.8) ax.set_ylabel(Equivalent Water Height (cm)) ax.set_title(Amazon Basin Water Storage Anomaly from CSR mascon) ax.legend()这里的fitted是上一节用最小二乘拟合出来的重构信号。图上如果趋势线的斜率和原始序列的长期变化方向一致说明信号稳健。如果原始序列前几年和后几年的均值明显不同可能是真实的气候或水文信号也可能是不同时段GRACE数据质量有差异需要谨慎评估。4.2 空间趋势图看出空间分布的不均匀性空间趋势图的价值在于揭示区域内部的不均一性。比如整个亚马逊流域可能平均趋势接近零但南部可能偏干、北部可能偏湿空间视角能避免被区域平均掩盖结构。用trend_grid.plot()时建议设置cbar_kwargs{label: Trend (cm/year)}单位直观。画图时的配色选择也有一点讲究在白色背景上使用红蓝发散色带如RdYlBu、RdBu能清楚区分正负趋势。需要控制颜色对称否则一条幅度的正趋势配一条大负趋势时视觉上会失真vmax np.nanmax(np.abs(trend_grid.values)) trend_grid.plot(cmapRdBu_r, vmin-vmax, vmaxvmax)4.3 结果验证和别的产品比一比自己处理后得到的结果一定要做交叉验证。最简单的办法是同一时间段、同一区域用JPL mascon或者GSFC mascon的公开数据重复一次区域平均流程对比趋势和季节振幅是否一致。不同产品的差异通常在10%~20%以内如果差别太大多半是处理流程有bug。还有一个低频的检验办法用GRACE结果和地面观测如地下水位井数据、水文模型输出做相关分析。以地下水研究为例GRACE反演的水储量变化和地下水位井观测的波动应该有一定相关性。但注意GRACE测的是总水储量变化土壤水地下水地表水与单一地下水位深度的直接对比需要谨慎。5. 常见问题与避坑指南5.1 尺度因子到底用不用JPL mascon官方明确说明需要应用scale factor尺度因子而CSR mascon一般来说不需要额外的全局尺度因子因为它已经通过约束处理把信号恢复到接近真实水平。但不同版本的CSR mascon产品如果文件里包含scale factor变量建议查看官方文档后决定是否使用。不要想当然拿到文件直接看有没有这个变量再看文档有没有相关说明。5.2 为什么我算出来的趋势跟论文里的对不上这是新手最容易遇到的情况。排查思路按顺序走单位是cm还是mm经纬度范围是否一致时间跨度是否一致比如论文是2003到2016你的是2002到2024两端数据会影响趋势估计有没有做GIA冰后回弹校正虽然CSR mascon一般去掉了GIA影响但对比时也要看对方用的产品是否相同。5.3 GRACE和GRACE-FO接缝处的信号不连续GRACE和GRACE-FO之间有一个约11个月的数据缺失接缝前后如果出现信号跳变不要慌先判断是不是真实水文事件。处理上建议不要直接在断口做插值而是保持缺测状态在趋势拟合和滤波时正确跳过NaN。强行插值会引入虚假的低频信号后面怎么解释都会很别扭。5.4 常见问题速查表问题可能原因解决方式趋势符号不对单位或符号约定搞反检查ewh定义和单位空间图有大量条带未做滤波或使用了未处理球谐产品确认使用的是mascon网格产品区域平均结果跳动剧烈区域面积小或信号弱噪声占比大适当扩大区域或用更大平滑时间序列尾部数值异常GRACE-FO初期解算质量不稳定检查数据质量必要时截掉早期几个月折线图断断续续GRACE/GRACE-FO数据缺测使用缺测值不要fillna后直接画图5.5 数据量不大但别忽略数据版本CSR mascon文件本身不算特别大一般几百MB但数据版本处理要记录下来。建议在脚本开头用注释或配置写清楚产品名称、版本号、下载时间和处理参数。做科研的最终交付物里要能回溯到原始数据版本否则后期复现和审稿意见回应的成本会很高。5.6 批量处理的思路如果研究区域不止一个或者需要按不同经纬度窗口反复提取时间序列建议把第3节的区域平均流程封装成函数def extract_regional_series(ds, lat_range, lon_range): 从CSR mascon数据中提取指定区域的等效水高时间序列 region ds.sel(latslice(*lat_range), lonslice(*lon_range)) weights np.cos(np.deg2rad(region[lat])).fillna(0) series (region[ewh] * weights).sum(dim[lat, lon]) series series / (weights.sum() * len(region[lon])) return series.to_dataframe(nameewh).drop(columns[lat, lon])这样后续对任意流域批量处理时只需要循环传入经纬度范围即可。用xarray的接口保留维度坐标比用numpy裸数组挨个索引要清晰很多也不容易犯行列顺序搞错的低级错误。5.7 误差估计的简单思路CSR mascon产品文件中如果带误差变量可以直接用来做加权区域平均时的权重但更多时候我们需要自己估计时间序列的不确定性。一个实用的做法是计算残差原始值减趋势减季节项的标准差作为单月不确定度然后趋势不确定度用回归系数的标准误差。也可以用bootstrap重采样法对时间序列做1000次有放回抽样每次拟合趋势取趋势分布的标准差作为不确定度。这个方法很笨但可靠在论文里写起来也容易被审稿人接受。结尾这套CSR mascon数据处理流程我自己前前后后改过好几版从最早用MATLAB读文本格式到后来换成Python脚本加xarray时间开销省了不止一半。最大的体会是任何重力卫星数据产品拿到手不要急着算结果先把文件在NetCDF层面理解透把单位和经纬度约定确认好后面自然顺。另一个想强调的点是数据加工程序再简单也要保留原始数据版本和处理记录的注释不然半年后回来看脚本连当初用的是哪个版本的产品、坐标怎么转换的都忘干净了。这份流程只是一个起点你完全可以把它改造成适合自己研究区域的批量处理工具接上水位观测数据或者水文模型输出做联合分析。做数据加工这活讲究的就是一个清晰可控——能把每一步都讲清楚你的数据集就成功了一大半。本文还有配套的精品资源点击获取