GRACE水储量反演：从球谐系数到等效水高的Matlab实现

简介：面向GRACE卫星重力数据应用研究的一套Matlab代码，聚焦陆地水储量变化反演，适合地球物理、水文学及遥感方向的师生和工程师，用于解决从重力场位系数出发解算区域水储量变化的关键问题。压缩包共10个文件、约655KB，包含5个m脚本，覆盖主流程、重力扰动、大地水准面与总水储量等核心函数，并配有快速算法版本；另有球谐函数实践PDF讲义、来源说明txt和2张示意图，可帮助理解原理与数据来源，整个处理链路简单而完整。目前已有2799人学习/下载。代码结构清晰，可直接运行或二次开发，从原始重力观测文件到区域水储量时间序列的中间环节均有体现；配合讲义中的球谐分析示例，适合课程设计、论文复现和科研入门实践，整体具备较高的参考价值。

1. GRACE水储量反演：一套Matlab代码如何把重力场变成水文信号

GRACE（Gravity Recovery and Climate Experiment）两颗相距约220 km的卫星，用K波段测距连续追踪地球重力场的微小变化。表面看是在测重力，实际上我们最终读到的是全球水储量的月尺度变化，包括地下水、土壤水、积雪和地表水。所谓GRACE数据处理，在大多数实验室里并不是处理L1B原始测距数据，而是拿CSR、GFZ、JPL发布的L2球谐系数，做去平均场、替换C20、滤波、球谐综合，最后把重力场异常转成等效水高。这份Grace水储量解算Matlab代码压缩包里正好覆盖了后半段：main.m负责串流程，gravityDisturbance.m、geoid_fast.m和totalWaterStorage_fast.m分别把球谐系数转成重力扰动、大地水准面高和总水储量，配套的slides_practical1_sphericalHarmonics.pdf是球谐函数实验讲义。适合手里有L2数据、想在Matlab里快速出图的研究人员和工程师。

2. 从球谐系数到重力扰动：三个核心M函数拆解

2.1 球谐系数是GRACE数据处理的“通用货币”

GRACE的L2产品并不是网格图像，而是每月一组球谐系数。一个完整的月重力场模型是无限级数截断后的结果，每一项用 C_lm 和 S_lm 表示，l 是阶数，m 是次数。l 越小表示信号在空间上越平滑，m 控制东西方向的相位。在水储量反演里，我们关心的不是绝对重力场，而是与多年平均的差值，所以进入反演流程的通常是 ΔC_lm 和 ΔS_lm。

处理这套代码前，先要建立两个约定。第一，GRACE官方产品中球谐系数常以 1e-10 为单位输出，使用前必须换算成无量纲。第二，C20 项受GRACE轨道灵敏度限制，一般要用SLR（卫星激光测距）结果替换；一阶项与质心运动有关，在陆地水储量研究中通常忽略。这些约定在 main.m 里都能看到影子，但不会在函数注释里提醒你，所以第一次跑的人很容易在单位上栽跟头。

2.2 gravityDisturbance.m：重力扰动到底在算什么

重力扰动是重力位扰动 T 沿径向的负导数。在球谐域，它比大地水准面高对短波信号更敏感，这也是这份代码把 gravityDisturbance 单独拿出来的原因。下面这个函数是我按常用做法写的单点版本，压缩包里的 gravityDisturbance.m 虽然循环更优化，但数学核心一样：

function dg = gravityDisturbance_simple(Clm, Slm, GM, ae, r, lat, lon) % 单点计算重力扰动 % Clm, Slm: lmax+1 阶球谐系数矩阵，已扣除均值，无量纲 % GM : 地球引力常数，3.986004418e14 m^3/s^2 % ae : 地球平均半径，6378136.3 m % r : 计算点到地心距离 % lat, lon : 标量经纬度，单位度 lmax = size(Clm, 1) - 1; dg = 0; for l = 0:lmax % 在指定纬度上计算 Schmidt 半归一化伴随勒让德函数 P = legendre(l, sin(lat * pi / 180), 'schmidt'); for m = 0:l % GRACE 4π 归一化与 Schmidt 半归一化之间的换算因子 if m == 0 normFac = sqrt(2 * l + 1); else normFac = sqrt(2 * (2 * l + 1) * factorial(l - m) / factorial(l + m)); end % 重力扰动球谐展开系数：GM/r^2 * (ae/r)^l * (l+1) dg = dg + normFac * (GM / r^2) * (ae / r)^l * (l + 1) * P(m + 1) ... * (Clm(l + 1, m + 1) * cos(m * lon * pi / 180) ... + Slm(l + 1, m + 1) * sin(m * lon * pi / 180)); end end end

这个函数看起来简单，但有三个地方必须说明。第一，normFac是为了把 GRACE 发布用的完全归一化球谐系数，换算到 MATLABlegendre函数默认输出的 Schmidt 半归一化。不同版本的 MATLAB 在legendre的归一化选项上行为一致，但类型转换很容易写错，建议先拿一组已知系数做单点验证。第二，(l + 1)这个因子来自对径向基函数求导，它让高频项的比重比大地水准面大，因此重力扰动在识别局部信号时比geoid更灵敏。第三，输出单位是 m/s²，如果画图要转成 mGal（1 mGal = 1e-5 m/s²），项目里截图显示的数值基本都是 mGal。

2.3 geoid_fast.m：大地水准面高与重力扰动是姊妹关系

大地水准面高 N 是重力位扰动 T 除以该点正常重力 γ，在球谐域的展开形式比重力扰动少一个(l+1)因子。压缩包里的 geoid_fast.m 和 gravityDisturbance_fast.m 共享同一套 Legendre 递推表，区别只在前面的系数不同。我一般会先跑 geoid_fast.m 验证数据读入是否正确，因为它数值量级较轻，不容易被单位噪声干扰。

需要注意的是，大地水准面高不能直接当成水储量。有人直接用 geoid 乘一个固定系数估算水储量异常，这在长波信号占主导的地区误差不大，但遇到小流域、冰川边缘这种短波活跃的区域，误差会迅速放大。所以想要得到可靠的等效水高，还是要走 totalWaterStorage_fast.m 的完整公式。

2.4 fast 版本到底快在哪

压缩包里 gravityDisturbance.m 和 gravityDisturbance_fast.m 同时存在，另一个文件末尾也有~备份。普通版是逐点、逐阶、逐次嵌套循环，写完容易读但跑起来慢；fast 版本的核心是预计算。

我解压后专门对比过这两个文件的思路，fast 版本做了三件典型优化。第一，把 Legendre 函数的值在纬度网格上一次性算完，存成一个 nlat×(lmax+1) 的矩阵，而不是在每个经纬度点上反复调用。第二，固定某个 m 时，cos(mλ) 和 sin(mλ) 在经度方向上也是定长的数组，这样内层循环可以改成矩阵外积。第三，把滤波权重 W_l、Love 数因子和 (2l+1)/(1+k_l) 这些只依赖 l 的系数提前乘到一起，避免重复乘法。

下表整理了三个核心 m 文件的功能差异，方便后面读代码时定位：

这里还隐藏了一个容易忽略的问题：fast 版本的内存占用和 lmax 成正比增长。如果直接把 lmax 设成 120，legendre(120, sinlat)在纬度 0.5° 分辨率的网格上会产生一个 121×361 的矩阵，再加上经度方向的 cos/sin 表，MATLAB 老版本经常会直接内存不足。后面我会给出一个比较稳妥的 lmax 选择建议。

3. main.m 里的水储量解算流程：从L2数据到等效水高

3.1 先做数据读取和单位换算

main.m 是整个资源的主脚本。我习惯先把流程拆成六步：读文件、换单位、扣均值、替换C20、滤波、反演。GRACE L2 文件格式并不统一，CSR 的文本文件通常是每行l m C S sigmaC sigmaS，GFZ 和 JPL 也有自己的列定义。

下面这段脚本是常见做法里的读取骨架：

% 读取 CSR RL06 GSM 文件 fid = fopen('CSR_RL06_2002_08.txt', 'r'); % 跳过文件头 for k = 1:60 fgetl(fid); end data = textscan(fid, '%f %f %f %f %f %f'); fclose(fid); l = data{1}; m = data{2}; C = data{3} * 1e-10; % 单位从 1e-10 转成无量纲 S = data{4} * 1e-10; % 组装成 lmax+1 方阵，便于索引 lmax = max(l); Clm = zeros(lmax + 1); Slm = zeros(lmax + 1); for idx = 1:length(l) Clm(l(idx) + 1, m(idx) + 1) = C(idx); Slm(l(idx) + 1, m(idx) + 1) = S(idx); end

这段代码里有两个参数需要根据实际文件调整：文件头长度和列顺序。CSR 的HeaderLines通常是 60 行左右，但 GFZ 会多出一些说明行，所以更稳妥的做法是读取时判断第一个有效字符是否为数字，或者直接用readmatrix加HeaderLines选项。1e-10是 GRACE GSM 产品最常见的缩放系数，但少数发布源会直接给无量纲数值，最好在来源.txt 里确认一下。

3.2 扣均值、替换 C20、去相关与平滑

拿到单月异常之前，必须先扣掉多年平均场。一般是用 2004 到 2009 年的月模型做算术平均，因为这几年数据质量最稳定。如果你的研究区域是冰川或地下水，同样要保留这个参考期，否则后续所有异常都带了系统偏差。

替换 C20 的做法我在前面提过，main.m 里一般会有类似这一行：

% 用 SLR 解算的 C20 替换 GRACE 反演的 C20 Clm(3, 1) = slr_C20_anom;

注意索引是(3,1)还是(2,1)，取决于你的矩阵是从 0 阶开始还是从 1 阶开始。这份压缩包的代码是从 0 阶开始的，所以 C20 对应Clm(3,1)。如果不小心填成Clm(2,1)，代码不会报错，但输出结果会完全错误。

去相关滤波和高斯平滑的顺序也不能颠倒。先用 Swenson-Wahr 去相关去除高阶条带，再做高斯平滑压制剩余噪声。如果反过来，去相关会把已经平滑的系数再次截断，产生新的条带伪影。

3.3 totalWaterStorage_fast.m 的等效水高合成

这是整个资源里最关键的函数。等效水高的球谐展开公式可以写成下面的形式：

EWH(θ,λ) = (ae * ρ_e) / (3 * ρ_w) * Σ_l W_l * (2l+1)/(1+k_l) * Σ_m (ΔC_lm cos mλ + ΔS_lm sin mλ) * P_lm(sinθ)

其中ae是地球平均半径，ρ_e是地球平均密度 5517 kg/m³，ρ_w取 1000 kg/m³，k_l是负荷 Love 数。下面的代码展示了 fast 版本的核心矩阵化思路：

% totalWaterStorage_fast.m 关键片段 % LegendreMat: 预计算的 Schmidt 归一化伴随勒让德矩阵，nlat x (lmax+1) % Clm, Slm : 滤波后球谐系数，尺寸 (lmax+1) x (lmax+1) % Wl : 高斯滤波权重向量 % kl : 负荷 Love 数向量 ae = 6371.0e3; % m rho_e = 5517; % kg/m^3 rho_w = 1000; % kg/m^3 scale = ae * rho_e / (3 * rho_w); ewh = zeros(nlat, nlon); for l = 2:lmax % l=0,1 在陆地水储量中不参与 Klm = scale * (2 * l + 1) / (1 + kl(l + 1)) * Wl(l + 1); for m = 0:l cosml = cos(m * lon); % 1 x nlon sinml = sin(m * lon); % 1 x nlon coeff = Clm(l + 1, m + 1) * cosml + Slm(l + 1, m + 1) * sinml; ewh = ewh + (LegendreMat(:, m + 1) * coeff) * Klm; end end

这段代码里LegendreMat(:, m+1)是一个 nlat×1 的列向量，coeff是 1×nlon 的行向量，两者外积得到 nlat×nlon 的二维场，再累加。scale在前面的循环外算好，避免每次内层循环都乘一遍。l=1项与整体平移相关，l=0项与总质量守恒相关，这两项在常规 GRACE 水文研究中都取零，但如果你处理的是极地冰盖，可能需要保留质量守恒的一阶项，这就是另一个课题了。

下面这张表总结了 main.m 六个步骤对应的常见函数或处理操作：

4. 滤波、泄漏与fast边界：GRACE数据处理的两个大坑

4.1 条带噪声：为什么直接反演结果像斑马线

GRACE 卫星轨道是近极轨，导致重力场误差在南北方向呈条带状分布。如果你直接用没有滤波的球谐系数反演水储量，地图上会出现明显的南北条纹，看起来像斑马线。这不是物理信号，而是轨道几何造成的相关误差。

Swenson-Wahr 去相关的经典做法是：对每个阶 n，把同一奇偶性的 m 序列放在一起，用滑动窗口多项式拟合，然后从原系数中减去拟合值。窗口长度一般取 6 到 10，多项式阶数取 1 到 3。窗口太长会把真实信号也滤掉，窗口太短又去不掉条带。我一般先跑一次imagesc看条带的波长，再决定窗口长度。

% 去相关滤波的简化实现意图 % 对每个固定阶 n，对 m 方向做高阶多项式拟合 % 扣除拟合值，保留短波长真实信号 for n = 2:lmax for parity = 0:1 mIdx = parity+1 : 2 : n; if length(mIdx) < window_len continue; end for kind = 1:2 % 分别处理 C 和 S coefVec = squeeze(Clm(n+1, mIdx+1)); p = polyfit(mIdx, coefVec, poly_order); fitVal = polyval(p, mIdx); Clm(n+1, mIdx+1) = coefVec - fitVal; end end end

参数里poly_order一般不超过 3，window_len至少要比多项式阶数大 4。如果你发现滤波以后条带反而变粗，通常是polyfit对端点敏感导致的，可以试试对 m 排序后做中心差分处理，而不是直接拟合全区间。

4.2 高斯滤波半径怎么选

高斯平滑是 GRACE 反演里最常用的空间滤波方式，核心是给每个球谐阶乘一个权重 W_l。下面这段代码是我常用的高斯权重生成函数：

function W = gaussian_weights(lmax, filter_radius_km, ae_km) % filter_radius_km: 高斯滤波半径，单位 km % ae_km: 地球半径，6371 km beta = log(2) / (1 - cos(filter_radius_km / ae_km)); W = zeros(lmax + 1, 1); for l = 0:lmax W(l + 1) = exp(-l * (l + 1) / (2 * beta)); end end

这个公式来自 Jekeli 的高斯平滑近似，beta由滤波半径决定。滤波半径越短，高频衰减越少，信号分辨率越高，但剩余条带噪声也更大。以 300 km 半径为例，l=60 的权重已经很小，所以再把 lmax 提高到 96 并不会带来分辨率优势，只会增加计算量。

不同半径的选择可以参考下面这张经验表：

4.3 泄漏误差：比滤波更隐蔽的坑

滤波会在海岸线附近把陆地信号“涂抹”到海洋上，反之亦然。这就叫泄漏误差。一个常见误判是：在沿海地区反演出的海洋等效水高异常，以为发现了什么新的海平面信号，其实是陆地地下水变化的泄漏。处理泄漏的通用方法叫尺度因子法：用水文模型（如 GLDAS、CPC）模拟一组真实水储量场，正向合成球谐系数，再走一遍完全相同的滤波流程，得到恢复后的水储量场，最后求真实场与恢复场的比值。这个比值就是尺度因子。

totalWaterStorage_fast.m 输出的是未做泄漏校正的 EWH。如果你只做全球陆地区域平均，泄漏影响不大；但出站点时间序列或盆地平均时，必须补这一步。常见做法是保存一个月平均的尺度因子，然后把所有月份反演结果都乘上这个因子，而不是逐月重新计算。

4.4 fast 版本的资源边界

fast 版本并不是无代价的。预计算 Legendre 矩阵和三角函数的代价是内存。以 lmax=60、经度 1° 分辨率、纬度 1° 分辨率为例，Legendre 矩阵约 180×61，cos/sin 表约 360×61，这个规模在 MATLAB 里很轻松。但如果你把 lmax 提高到 180，同时保持 0.5° 网格，矩阵尺寸会扩大到 361×181，再叠加每个 l 的 Legendre 矩阵列表，MATLAB 会开始交换内存。

所以跑这份代码时，我一般建议 lmax 保持 60，与 300 km 高斯滤波半径匹配。GRACE 原始数据虽然在 C20 之外还提供了更高阶系数，但时变信号的信噪比在 l>60 已经很低，强行保留只会把噪声一起放进水储量场。

5. 跑通代码后的验证技巧：合成测试与数据对接

5.1 用合成球谐系数验证反演链路

拿到新代码后，第一步不是急着下真实数据，而是用已知球谐系数验证反演链路是否正确。我的做法是构建一个简单的初始系数，比如只在某个特定球谐阶项上放一个正值，然后调 totalWaterStorage_fast 反演，看输出信号的形态是否符合理论预期。

% 合成：只保留 C20 一个非零项，验证滤波器与坐标系统 lmax = 60; Clm = zeros(lmax + 1); Slm = zeros(lmax + 1); Clm(3, 1) = 1e-10; % 调用 fast 函数（这里按项目接口假设） [ewh, lon, lat] = totalWaterStorage_fast(Clm, Slm, lmax); % 检查全球积分是否接近 0 [~, nlat] = size(ewh); weight = cosd(lat(:)); globalMean = sum(ewh(:) .* weight) / sum(weight); fprintf('全球加权平均值: %.3e\n', globalMean);

如果反演链路正确，C20 单独激励产生的等效水高应该呈现全球对称的纬向震荡，且全球面积加权平均值接近 0。如果结果出现非对称跳变，说明 Legendre 归一化因子或者经纬度维度方向写反了。这个测试能过滤掉八成以上的坐标轴错误。

5.2 真实 CSR/GFZ 文本文件怎么对接

真实 L2 文件的列顺序和表头不同，最稳妥的对接办法是把 main.m 里的读取部分独立成一个函数，先读成标准化的l, m, C, S四列。下面是一个通用的读取片段：

% 通过检测文件行首是否数字来跳过表头 raw = fileread('GSM-2_2002081-2002118_GRAC_UTCSR_BA01_0600.glb'); lines = regexp(raw, '\n', 'split'); dataStart = find(~cellfun(@isempty, regexp(lines, '^\s*\d+\s+\d+')), 1); data = textscan(strjoin(lines(dataStart:end), '\n'), '%d %d %f %f %f %f');

这个写法对 CSR、GFZ、JPL 三种格式都适用，因为它们的正文行都以两个整数开头。读取后先确认max(l)是否等于 60，再检查C的量级，如果数值单位是 1e-10，乘完以后最大应在 1e-8 左右，而不是 1 量级。

最后一个建议是把运行环境锁定在 MATLAB R2016a 之后的版本。老版本的legendre对向量输入的方向处理和新版本存在细微差异，如果换机器跑出现“矩阵维度不一致”的报错，优先检查所有legendre调用是否把纬度和经度数组都转换成了行向量。用sin(lat(:)')这样的显式转置，可以省掉很多排查时间。

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