企业官网建设流程全解析

简介：本资源是面向地球物理、冰川学、水文学及气候学等领域科研人员的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_mfdataset：

import xarray as xr files = sorted(glob.glob("CSRM_BA01_*.nc")) ds_all = xr.open_mfdataset(files, concat_dim="time", combine="nested") ds_all = ds_all.sortby("time")

这里有个细节。open_mfdataset默认会尝试用坐标对齐，但有时候各个文件的lat、lon因为浮点存储精度略有差异，导致拼接时报错或产生NaN。稳妥做法是确认所有文件都来自同一版本，必要时可以先用ds.load()把单个文件读入内存，检查坐标是否一致，再批量拼接。

时间轴重构也很重要。NetCDF里的time通常是自2002-01-01起算的天数，读取后建议显式转换为cftime或pandas的DatetimeIndex：

time_idx = xr.cftime_range(start="2002-01", periods=ds_all.sizes["time"], freq="MS") ds_all = ds_all.assign_coords(time=time_idx)

月度数据的时间戳一般取每月1号。需要注意GRACE数据在某些月份缺失，拼接后时间轴会出现空洞。处理时不要直接用resample把空洞填掉，先明确哪些月份缺失，再决定插值还是保留NaN。

3.3 掩膜与区域提取：流域尺度的裁剪

拿到全球网格后，提取特定流域的方式有两种：一种是直接用经纬度范围裁剪矩形框，另一种是用流域边界矢量做掩膜。

矩形裁剪适合初步查看：

lon_min, lon_max = 90, 122 lat_min, lat_max = 21, 36 region = ds_all.sel( lon=slice(lon_min, lon_max), lat=slice(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)), dtype=bool) # 用点是否落在多边形内来构建掩膜 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], rcond=None) 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, f"Unknown 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是否已经扣除？不要再扣一遍

冰川均衡调整（GIA，Glacial 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_email@example.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_dim="time", combine="nested") 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)), transform=from_origin(lon.min(), lat.max(), abs(lon[1]-lon[0]), abs(lat[1]-lat[0])), fill=0, dtype=np.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], weights=weights_2d[valid_t])

再做趋势拟合就很简单了，沿用第3节里的最小二乘模型。长江流域GRACE时间序列的特点是季节变化非常明显，夏季水储量高、冬季低，但常年趋势相对较小。做趋势分析时建议以年为单位输出趋势，并给出95%置信区间。

6.3 可视化和导出

最终图表我一般画成上下两部分：上面是流域平均等效水高的逐月曲线，下面是用颜色表示趋势的空间分布图。空间图可以用cartopy绘制，叠加流域边界和主要水系。

导出数据时，最终数据集用NetCDF格式保存，时间序列再额外导出一份CSV：

df = pd.DataFrame({ "time": ds_all["time"].values, "basin_mean_ewh_cm": basin_series }) df.to_csv("output/yangtze_basin_ewh_monthly.csv", index=False)

NetCDF导出时别忘了把前面提到的元数据写进去。

6.4 验证思路与常见误区

加工完的数据，一定要做一次外部验证，不能画完图就收工。我常用的验证方式有三个：

一是和卫星水文模型对比。比如把CSR mascon的流域平均水储量变化和GLDAS水文模型、PCR-GLOBWB模型对比，看季节相位和年际波动是否一致。如果mascon趋势和模型趋势符号相反，先别着急说模型不准，回头检查自己的时间序列是否扣错了参考时段。

二是和GRACE官方发布的同区域时间序列对比。CSR官网提供了一些流域或全球平均的时间序列，可以直接拿来对比，数量级应该非常接近。

三是和局部地面监测数据对比。比如大型水库的蓄水量变化、地下水井水位变化，对比时要注意空间尺度差异，地面点数据往往不能完全代表整个流域的平均状态。

常见误区主要是三类：把0.25度网格当真实分辨率，误以为能识别几十公里尺度的局部信号；直接用全球文件做流域平均而不做面积加权，导致高纬度网格贡献偏大；以及把GRACE和GRACE-FO期间的数据直接拼起来做趋势，没有处理任务间的系统差，造成2018年前后出现人为跳变。

数据加工这个环节说起来不复杂，绝大部分时间都耗在“确认单位”、“确认参考时段”、“确认GIA处理”这类的检查上。我个人体会是，把这些琐碎的确认做成必须执行的检查函数，写进处理流程的最前面，后面再大的数据量都能跑得安心。希望这篇记录能帮你少走几步弯路，后面有时间我再单独写写CSR、JPL、GSFC三家mascon产品的对比和实际差异。

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