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（德克萨斯大学奥斯汀分校空间研究中心）、JPL（NASA喷气推进实验室）、GSFC（NASA戈达德太空飞行中心）。三家思路有差异，结果自然也有细微不同。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_north
• lon：经度，单位degree_east
• time：时间变量，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到90，lon是0到360还是-180到180，time的时间跨度是否覆盖你需要的时段。如果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, unit='D', origin=pd.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(lat=slice(lat_min, lat_max), lon=slice(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(dim='lat').sum(dim='lon') / (weights.sum() * len(region['lon'])) # 转换为DataFrame，方便后续分析 ts = regional_ewh.to_dataframe(name='ewh').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, rcond=None) # 趋势项（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']], vectorize=True, dask='parallelized', 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-', markersize=3, label='Monthly', alpha=0.7) ax.plot(ts.index, ts['ewh'].rolling(3, center=True).mean(), 'g-', linewidth=2, label='3-month running mean') ax.plot(ts.index, fitted, 'r-', linewidth=2, label='Trend + seasonal fit') ax.axhline(0, color='gray', linestyle='--', linewidth=0.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(cmap='RdBu_r', vmin=-vmax, vmax=vmax)

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 常见问题速查表

5.5 数据量不大，但别忽略数据版本

CSR mascon文件本身不算特别大，一般几百MB，但数据版本处理要记录下来。建议在脚本开头用注释或配置写清楚产品名称、版本号、下载时间和处理参数。做科研的最终交付物里要能回溯到原始数据版本，否则后期复现和审稿意见回应的成本会很高。

5.6 批量处理的思路

如果研究区域不止一个，或者需要按不同经纬度窗口反复提取时间序列，建议把第3节的区域平均流程封装成函数：

def extract_regional_series(ds, lat_range, lon_range): """从CSR mascon数据中提取指定区域的等效水高时间序列""" region = ds.sel(lat=slice(*lat_range), lon=slice(*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(name='ewh').drop(columns=['lat', 'lon'])

这样后续对任意流域批量处理时，只需要循环传入经纬度范围即可。用xarray的接口保留维度坐标，比用numpy裸数组挨个索引要清晰很多，也不容易犯行列顺序搞错的低级错误。

5.7 误差估计的简单思路

CSR mascon产品文件中如果带误差变量，可以直接用来做加权区域平均时的权重，但更多时候我们需要自己估计时间序列的不确定性。一个实用的做法是计算残差（原始值减趋势减季节项）的标准差作为单月不确定度，然后趋势不确定度用回归系数的标准误差。也可以用bootstrap重采样法，对时间序列做1000次有放回抽样，每次拟合趋势，取趋势分布的标准差作为不确定度。这个方法很笨但可靠，在论文里写起来也容易被审稿人接受。

结尾

这套CSR mascon数据处理流程，我自己前前后后改过好几版，从最早用MATLAB读文本格式，到后来换成Python脚本加xarray，时间开销省了不止一半。最大的体会是，任何重力卫星数据产品，拿到手不要急着算结果，先把文件在NetCDF层面理解透，把单位和经纬度约定确认好，后面自然顺。另一个想强调的点是，数据加工程序再简单，也要保留原始数据版本和处理记录的注释，不然半年后回来看脚本，连当初用的是哪个版本的产品、坐标怎么转换的都忘干净了。这份流程只是一个起点，你完全可以把它改造成适合自己研究区域的批量处理工具，接上水位观测数据或者水文模型输出做联合分析。做数据加工这活，讲究的就是一个清晰可控——能把每一步都讲清楚，你的数据集就成功了一大半。

本文还有配套的精品资源，点击获取