简介面向GRACE卫星重力数据应用研究的一套Matlab代码聚焦陆地水储量变化反演适合地球物理、水文学及遥感方向的师生和工程师用于解决从重力场位系数出发解算区域水储量变化的关键问题。压缩包共10个文件、约655KB包含5个m脚本覆盖主流程、重力扰动、大地水准面与总水储量等核心函数并配有快速算法版本另有球谐函数实践PDF讲义、来源说明txt和2张示意图可帮助理解原理与数据来源整个处理链路简单而完整。目前已有2799人学习/下载。代码结构清晰可直接运行或二次开发从原始重力观测文件到区域水储量时间序列的中间环节均有体现配合讲义中的球谐分析示例适合课程设计、论文复现和科研入门实践整体具备较高的参考价值。1. GRACE水储量反演一套Matlab代码如何把重力场变成水文信号GRACEGravity 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: lmax1 阶球谐系数矩阵已扣除均值无量纲 % 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 * (l1) 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²如果画图要转成 mGal1 mGal 1e-5 m/s²项目里截图显示的数值基本都是 mGal。2.3 geoid_fast.m大地水准面高与重力扰动是姊妹关系大地水准面高 N 是重力位扰动 T 除以该点正常重力 γ在球谐域的展开形式比重力扰动少一个(l1)因子。压缩包里的 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×(lmax1) 的矩阵而不是在每个经纬度点上反复调用。第二固定某个 m 时cos(mλ) 和 sin(mλ) 在经度方向上也是定长的数组这样内层循环可以改成矩阵外积。第三把滤波权重 W_l、Love 数因子和 (2l1)/(1k_l) 这些只依赖 l 的系数提前乘到一起避免重复乘法。下表整理了三个核心 m 文件的功能差异方便后面读代码时定位文件输入输出特征gravityDisturbance.m球谐系数、GM、ae、r、lat/lon重力扰动m/s²逐点循环适合验证gravityDisturbance_fast.m同样输入但 lat/lon 可为网格重力扰动网格预计算 Legendre 矩阵geoid_fast.m球谐系数、GM、γ、lat/lon大地水准面高m没有 (l1) 因子totalWaterStorage_fast.m滤波后球谐系数、Love 数、密度等效水高m反演水储量主函数这里还隐藏了一个容易忽略的问题fast 版本的内存占用和 lmax 成正比增长。如果直接把 lmax 设成 120legendre(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 sigmaSGFZ 和 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; % 组装成 lmax1 方阵便于索引 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 * (2l1)/(1k_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 (lmax1) % Clm, Slm : 滤波后球谐系数尺寸 (lmax1) x (lmax1) % 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 % l0,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(:, m1)是一个 nlat×1 的列向量coeff是 1×nlon 的行向量两者外积得到 nlat×nlon 的二维场再累加。scale在前面的循环外算好避免每次内层循环都乘一遍。l1项与整体平移相关l0项与总质量守恒相关这两项在常规 GRACE 水文研究中都取零但如果你处理的是极地冰盖可能需要保留质量守恒的一阶项这就是另一个课题了。下面这张表总结了 main.m 六个步骤对应的常见函数或处理操作流程节点输入常见操作输出读取 L2文本或 netCDFtextscan / readmatrixl, m, C, S单位换算C, S乘 1e-10无量纲系数扣均值所有月系数减去 2004-2009 平均值ΔC_lm, ΔS_lm替换 C20对应项SLR 结果覆盖修正后的 ΔC20滤波全部系数Swenson-Wahr 高斯平滑后的系数反演滤波后系数totalWaterStorage_fastEWH 网格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 parity1 : 2 : n; if length(mIdx) window_len continue; end for kind 1:2 % 分别处理 C 和 S coefVec squeeze(Clm(n1, mIdx1)); p polyfit(mIdx, coefVec, poly_order); fitVal polyval(p, mIdx); Clm(n1, mIdx1) coefVec - fitVal; end end end参数里poly_order一般不超过 3window_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 半径为例l60 的权重已经很小所以再把 lmax 提高到 96 并不会带来分辨率优势只会增加计算量。不同半径的选择可以参考下面这张经验表滤波半径适合场景副作用100-200 km大型河流流域、湖泊条带噪声明显泄漏严重300 km大陆尺度水储量综合平衡500 km跨区域干旱/极端气候研究信号被过度平滑750 km 以上海平面与全球质量迁移短波信号几乎丢失4.3 泄漏误差比滤波更隐蔽的坑滤波会在海岸线附近把陆地信号“涂抹”到海洋上反之亦然。这就叫泄漏误差。一个常见误判是在沿海地区反演出的海洋等效水高异常以为发现了什么新的海平面信号其实是陆地地下水变化的泄漏。处理泄漏的通用方法叫尺度因子法用水文模型如 GLDAS、CPC模拟一组真实水储量场正向合成球谐系数再走一遍完全相同的滤波流程得到恢复后的水储量场最后求真实场与恢复场的比值。这个比值就是尺度因子。totalWaterStorage_fast.m 输出的是未做泄漏校正的 EWH。如果你只做全球陆地区域平均泄漏影响不大但出站点时间序列或盆地平均时必须补这一步。常见做法是保存一个月平均的尺度因子然后把所有月份反演结果都乘上这个因子而不是逐月重新计算。4.4 fast 版本的资源边界fast 版本并不是无代价的。预计算 Legendre 矩阵和三角函数的代价是内存。以 lmax60、经度 1° 分辨率、纬度 1° 分辨率为例Legendre 矩阵约 180×61cos/sin 表约 360×61这个规模在 MATLAB 里很轻松。但如果你把 lmax 提高到 180同时保持 0.5° 网格矩阵尺寸会扩大到 361×181再叠加每个 l 的 Legendre 矩阵列表MATLAB 会开始交换内存。所以跑这份代码时我一般建议 lmax 保持 60与 300 km 高斯滤波半径匹配。GRACE 原始数据虽然在 C20 之外还提供了更高阶系数但时变信号的信噪比在 l60 已经很低强行保留只会把噪声一起放进水储量场。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(:))这样的显式转置可以省掉很多排查时间。本文还有配套的精品资源点击获取