2026/9/17 3:42:52
FEATURED · 精选文章 
简介这套代码专注于处理重力恢复与气候实验卫星的观测数据面向地球物理、水文及相关专业的学生和科研人员用于将重力观测转化为全球水储量变化信息支持地下水消耗、冰川融化、干旱洪水等研究。压缩包包含10个文件总大小仅655KB其中有5份Matlab程序脚本、2份过程示意图、1份算法说明文档和1个来源说明文件覆盖数据读取、重力场展开、扰动计算、水储量反演等环节。该工具包已有2799人学习下载适合需要快速完成GRACE数据处理实验的入门者。根据内容预览主控脚本能够串联球谐系数修正、大地水准面解算与总水储量快速估计的完整流程而随包文档和截图对滤波、去噪、坐标转换等关键细节做了可视化解释可帮助读者理解每一步的物理含义并在此基础上适配自己的研究区域或数据集。整体小而完整兼具教学演示与科研参考价值。1. 用 Matlab 解算 GRACE 水储量的第一道坎GRACE 重力卫星已经退役但它的数据仍然是全球陆地水储量变化监测的黄金标准。做水文、冻土、干旱监测的人拿到 Level-2 球谐系数之后第一步往往不是科学问题而是卡在“怎么把几十 MB 的 GSM 文件变成一条时间序列”上。这个过程的门槛在于球谐系数不能直接用必须先做去相关滤波和高斯平滑再把位系数转成等效水高EWh单位、阶数截断、纬度加权、泄漏误差修正任何一步错了结果都会漂移得没法看。常见做法是用 CSR、GFZ 或 JPL 发布的 Level-2 RL06 数据配合 DDK 滤波或高斯平滑在 Matlab 里逐月处理。新手能跟步骤走熟手则会关心滤波半径选多大、截断阶数取多少、泄漏修正有没有必要。这套流程大概是读 GSM 文件、处理 C 和 S 阶系数、做滤波、球谐合成、换算等效水高、再按流域或格网统计。写成代码并不长但参数背后全是物理。2. GRACE 数据产品选型与预处理CSR、GFZ、JPL 怎么选GSM 文件格式怎么读2.1 机构产品差异与 RL06 版本选择GRACE 数据由三家机构官方发布CSR美国德州大学、GFZ德国地学中心、JPL美国喷气推进实验室。三者 Level-2 产品处理细节不同但球谐系数格式一致都遵循 RL06 标准。RL06 相比 RL05 最大的改进是背景模型替换为 AOD1B RL06同时更新了 C20 和 C30 系数。日常解算水储量我一般选 CSR 或 GFZ 的 RL06 数据二者噪声水平接近JPL 的格式多一列 sigma 误差处理时读取方式略有差异初学者先用 CSR 就好。三家机构的产品在文件命名上有统一格式以CSR_2015_01_0001_GSM_0002为例其中的GSM表示重力场模型月度解0002是处理版本号。RL06 的 C20 系数已被卫星激光测距SLR结果替代文件内会带FLAG字段标记处理时 C20 和 C30 通常直接采用文件内数值不用额外修正。2.2 GSM 文件头与块结构的 Matlab 读取方法GSM 文件是纯文本格式读取时需要先解析头部再读取 4 个数据块GRCOEF重力场系数、GRCOFS系数标准差、GRACEF加速度计标记、CALIBF校准参数。GRCOEF块内的每一行格式为l m C_lm S_lm系数值是无量纲的球谐系数单位是标准重力场模型使用的完全归一化形式。下面给出一个通用读取函数function [lm, C, S] read_gsm(filename) fid fopen(filename, r); % 逐行查找 GRCOEF 块起始位置 header textscan(fid, %s, 1, delimiter, \n); while ~strcmp(header{1}{1}(1:min(6,end)), GRCOEF) header textscan(fid, %s, 1, delimiter, \n); if feof(fid), error(No GRCOEF block found); end end % 读取系数直到下一个块标记或文件结束 lm []; C []; S []; while ~feof(fid) line fgetl(fid); if isempty(line), continue; end if ~strcmp(line(1:min(6,end)), GRCOFS) data sscanf(line, %d %d %f %f); if length(data) 4 lm(end1,:) data(1:2); C(end1) data(3); S(end1) data(4); end else break; end end fclose(fid); end这段代码的核心逻辑是先按行扫描定位GRCOEF块再从块起始位置连续读取四列数值行。sscanf的格式字符串写明了l m C S的顺序Matlab 会跳过空行和注释行。注意line(1:min(6,end))的写法是为了兼容行尾带空格的场景避免字符串索引越界。2.3 阶数截断与 C20 替换的预处理约定读入系数后第一件事不是滤波而是做两个决定截断到多少阶、是否替换 C20。GRACE 高阶系数噪声大常规做法是截断到 60 阶对应空间分辨率约 330 公里保守一点用 40 阶。截断操作直接丢弃 C/S 数组中阶数大于 60 的行即可在构造合成矩阵时也能省内存。C20 替换也是一个可选项RL06 数据里的 C20 已经是 SLR 解直接使用是合理的但如果你拿的是旧版 RL05 数据则必须用 SLR 提供的 C20 时间序列替换。预处理还有一个容易被忽略的细节GRACE 反演的 C00 项恒等于 1C10、C11、S11 在地心参考系下被强制归零。这是因为卫星重力无法独立测定地球质心运动这些项不包含地质信息。处理时不需要特殊移除这些项合成时自然被球谐函数积分掉。3. GRACE 水储量反演的核心公式与两类滤波从球谐系数到等效水高3.1 球谐合成公式为什么位系数要先乘负荷勒夫数等效水高Equivalent Water Height, EWh是全球质量重分布直接换算成水层厚度的结果。球谐域中每个阶次 (l,m) 的位系数变化量 ΔC_lm、ΔS_lm 通过如下公式映射为等效水高变化ΔEWh(θ,λ) (a·ρ_avg)/(3·ρ_w) · Σ_{l0}^{L} Σ_{m0}^{l} (2l1)/(1k_l) · [ ΔC_lm·cos(mλ) ΔS_lm·sin(mλ) ] · P̄_lm(cosθ)其中a是地球平均半径6378.1363 kmρ_avg是地球平均密度5517 kg/m³ρ_w是水密度1000 kg/m³k_l是负荷勒夫数load Love numberP̄_lm是完全归一化连带勒让德函数。乘(2l1)/(1k_l)这一步把“重力场变化”转成了“地表质量变化”是物理上最关键的一步。负荷勒夫数序列是固定常量CSR 官方的 RL06 处理文档中给出了 0-200 阶的表用load_love_numbers.m读入即可。3.2 高斯平滑与扇形滤波去相关去不掉的交给空间滤波球谐域内高阶系数的条带误差会让制图结果出现南北向的条纹经典的解决办法有两类。一类是去相关滤波比如 DDK 1/2/3 系列它利用系数的协方差信息做平滑不改变物理分辨率另一类是空间平滑最常见的是高斯滤波其本质是对球谐系数乘一个阶相关的权重因子W_l。两者的选择依据是目标流域面积大流域用 350-500 km 半径的高斯滤波中小流域用 200-300 km过度滤波会把真实信号也削掉。Matlab 里做高斯平滑不需要显式计算权重因子再逐阶相乘可以用gauss_smooth.m这类现成函数输入半径后输出一个与阶数等长的向量W_l。官方 DDK 滤波直接提供每个月的滤波后系数文件以DDK3命名读入后跳过滤波步骤只做高斯平滑即可。两条路线最终殊途同归去相关除掉条带空间平滑再压一次噪声。3.3 等效水高合成的快速实现用网格求值避开 legendre 循环球谐合成最直观的实现是双层循环逐格网点累加P_lm但这种方法在 1°×1° 网格和 60 阶截断下需要计算 64800×3661 次连带勒让德函数Matlab 循环跑起来可能要几分钟。更快的方式是先用legendre_array.m部分开源包提供例如 MATLAB GRACE Toolbox 中的函数一次性生成所有阶次在某一纬度上的 P_lm 值再利用矩阵乘法逐纬度合成。一个折中方案是在 Matlab 中调用legendre函数按纬度循环每纬度对的 P_lm 是向量仍比逐格点快两个数量级。实际项目中我常用以下实现function [lon, lat, ewh] sph_synth(lm, C, S, Lmax, R_filter) % 生成等间距经纬度网格默认 1 度 lat 89.5:-1:-89.5; lon 0.5:1:359.5; [LON, LAT] meshgrid(lon, lat); % 高斯滤波权重向量W_l 由半径 R_filterkm确定 W gauss_weights(Lmax, R_filter); % 预分配输出矩阵 ewh zeros(size(LON)); % 勒让德归一化常数映射表后续循环中重复使用 norm_factor legendre_norm(lm); for i 1:length(lat) theta (90 - lat(i)) * pi/180; P legendre(Lmax, cos(theta), sch); % 连带勒让德4 列 sum_val 0; for l 0:Lmax % 提取第 l 阶所有 m 的 P_lm P_lm P(l1, 1:l1); Clm C(lm(:,1)l, :); % 取该阶系数 % 合成C 项与 cos、S 项与 sin 结合 cos_term Clm(:,3) .* cos((0:l) .* LON(i,:)); sin_term Clm(:,4) .* sin((0:l) .* LON(i,:)); sum_val sum_val W(l1) * (2*l1)/(1love(l1)) ... * (sum(cos_term) sum(sin_term)); end ewh(i,:) sum_val * (6378.1363 * 5517 / 3000); end end上述代码里的legendre函数采用sch输出 Schmidt 归一化形式而 GRACE 标准用完全归一化两者差一个sqrt(2)因子需要注意。legendre_norm和love是预先从文件加载的表避免循环内重复计算。代码的核心思路是纬度循环里先用legendre获得该纬度的全部阶次值再在阶次循环中组合 C/S 系数和三角基。W(l1)是高斯权重(2*l1)/(1love(l1))是负荷 LOVE 数因子最后一个乘法里的5517/3000是ρ_avg/(3×ρ_w)化简结果。3.4 泄漏误差与信号恢复要不要做尺度因子修正反演出来的 EWh 经滤波后信号幅度会被削弱 30%-60%尤其是局部水储量变化大的区域如青藏高原、亚马逊流域。针对流域平均时间序列标准处理是做一个“尺度因子”修正对每个格网或流域先假设一个真实信号的时变模型比如 GLDAS 水文模型用同样的滤波处理得到“滤波后”的信号再用最小二乘拟合出比例系数k最后将观测值除以k。这个系数对强信号区域在 1.5-2.0 之间对弱信号区接近 1。对于单格网时间序列不推荐强行做尺度因子修正因为信噪比太低拟合出的系数不稳定。更准确的做法是按流域聚合后再修正或者直接与 GLDAS、GLDAS-2.1 的地表水储量输出做对比验证。泄漏误差除了幅度衰减还包括相邻区域信号的串扰比如海洋信号泄漏到沿海陆地这个目前没有一劳永逸的修正方案减少方式是选足够大的研究区。4. 用 Matlab 跑通 GRACE 水储量时间序列的完整流程从原始 GSM 到流域平均4.1 全流程代码数据下载组织、逐月反演与结果存储在动手算之前先约定目录结构。我习惯把不同月份的 GSM 文件统一放在raw/目录命名含年份月份方便批量读取。全流程代码可以拆为三个阶段预处理、滤波与合成、流域平均。下面的主脚本展示完整调度逻辑% GRACE 水储量解算主流程 % 目录raw/*.txt 存放 CSR RL06 GSM 文件 fnames dir(raw/GSM_*.txt); R_filter 300; % 高斯滤波半径单位 km Lmax 60; % 截断阶数 % 预加载负荷勒夫数与高斯权重避免循环内重复计算 love load(love_numbers_lmax200.txt); % 两列阶数k_l W gauss_weights(Lmax, R_filter); % 存放所有月份的 EWh 网格 ewh_all []; dates []; for i 1:length(fnames) % 从文件名解析年月 tok regexp(fnames(i).name, (\d{4})_(\d{2}), tokens); year str2double(tok{1}{1}); month str2double(tok{1}{2}); % 读取球谐系数 [lm, C, S] read_gsm([raw/ fnames(i).name]); % 截断到 Lmax idx lm(:,1) Lmax; lm lm(idx,:); C C(idx); S S(idx); % 计算等效水高 [lon, lat, ewh] sph_synth(lm, C, S, Lmax, R_filter, W, love); % 存入数组后续计算时间序列 ewh_all cat(3, ewh_all, ewh); dates [dates, datetime(year, month, 15)]; end % 保存结果 save(grace_ewh_global_1deg_300km.mat, lon, lat, ewh_all, dates);这段主脚本的调度很清楚先读取文件列表再对每个文件完成“读取-截断-合成-存储”四步。注意gauss_weights里的返回值需要显式传给合成函数避免每层循环重复算高斯因子。文件名的正则解析用regexp提取年月datetime(year, month, 15)选了每月 15 日作为该月时间戳方便后续画时间序列图。4.2 流域平均用经纬度边界框或 Shapefile 掩膜聚合拿到全球 EWh 网格后流域平均是水文应用最常见的一步。最简单的方法是用经纬度矩形框选目标区域适合流域形状接近矩形的区域更精确的方法是用 Shapefile 做掩膜再输出掩膜内的平均 EWh。下面给出掩膜聚合的代码% 流域平均以 Shapefile 区域的格网掩膜为例 shp shaperead(basin_boundary.shp); mask inpolygon(lon, lat, shp.X, shp.Y); % 计算掩膜面积权重纬度余弦加权 area_w cosd(lat); area_w_masked area_w .* mask; for i 1:size(ewh_all, 3) ewh_month ewh_all(:,:,i); % 避免 NaN 区域影响平均 valid ~isnan(ewh_month) mask; basin_ewh(i) sum(ewh_month(valid) .* area_w_masked(valid)) / ... sum(area_w_masked(valid)); end % 时间序列绘图 plot(dates, basin_ewh, LineWidth, 1.5); xlabel(Date); ylabel(EWh (cm)); grid on;这个加权平均逻辑很关键直接用mean会把高纬度格网权重放大因为 1° 格网的物理面积随纬度余弦变化真实流域平均必须用cos(lat)作为权重。inpolygon是基于经纬度平面的多边形判断天然适合 GRACE 1° 网格尺度。这一步输出的时间序列单位是厘米水柱与降水量、蒸散量对比时要注意量纲统一1 cm EWh 10 kg/m²。4.3 时间序列后处理季节性信号分离与长期趋势提取流域平均序列里通常包含年周期、半年周期和长期趋势。用最小二乘拟合提取趋势时常见回归模型中包含趋势项与年/半年谐波项方程如下EWh(t) β_0 β_1·t A_a·cos(2πt) B_a·sin(2πt) A_s·cos(4πt) B_s·sin(4πt) ε拟合结果是β_1即长期趋势单位通常转换为 cm/yrsqrt(A_a²B_a²)是年振幅。要在 Matlab 中实现使用设计矩阵和\运算符即可核心代码如下% 构建设计矩阵趋势 年周期 半年周期 t year(dates) (month(dates)-0.5)/12; X [ones(length(t),1), t, cos(2*pi*t), sin(2*pi*t), ... cos(4*pi*t), sin(4*pi*t)]; % 最小二乘求解 beta X \ basin_ewh; trend_mm_per_yr beta(2) * 1000; % 把 cm/yr 转成 mm/yr % 去季节信号后画残差序列 seasonal X(:,3:6) * beta(3:6); detrended basin_ewh - X(:,1:2)*beta(1:2) - seasonal;设计矩阵的构建方式要考虑时间变量的基准t用小数年表示并以年中点采样避免相位偏差。X \ basin_ewh在 Matlab 中是左除等效于最小二乘解不需要显式求逆。如果时间序列较短少于 3 年年周期的拟合不稳定趋势结果要谨慎解读优先看相对变化而不是绝对速率。5. 参数调优与验证技巧滤波半径、截断阶数与结果的可信度检查5.1 滤波半径与截断阶数的组合选择表GRACE 处理中最影响结果质量的两个参数是高斯平滑半径和截断阶数。两者的作用是叠加的截断阶数决定最短波长的下限L60 对应波长约 330 km高斯半径进一步衰减高阶能量。实践中常用组合如下应用场景截断阶数 L高斯半径 (km)适合的流域尺度全球尺度制图60300 20 万 km²区域干旱监测402505-20 万 km²大流域如亚马逊40500 100 万 km²中小流域研究601501-5 万 km²注意表格中“流域尺度”是保守估计实际可反演的最小流域与信号强度高度相关强信号区域如华北平原地下水开采区3 万 km² 也可能勉强分辨弱信号区 10 万 km² 也不一定稳定。选择参数时不要机械照搬先用研究区 3-4 个月的合成 EWh 图做目视检查条纹是否明显、信号是否连续。5.2 三种快速验证方法与官方 Mascon 对比、与 GLDAS 对比、自洽性检验日常验证中我最常用的是与 CSR/GFZ 官方发布的 Mascon 产品做对比。Mascon 产品是块体质量异常解空间分辨率约 1°但它的滤波处理与球谐法不同不会出现条带。对比方法是将球谐法结果重采样到 Mascon 网格做流域平均后看两条时间序列的相关性和均方根误差相关系数大于 0.8 且 RMSD 小于 3 cm 属于正常范围。另一个验证途径是与 GLDAS 水文模型对比重点看季节振幅和相位。GLDAS 反演的是土壤水、雪水当量、植被冠层水的总和与 GRACE 反演的总水储量差一个地表水项因此差值过大不必慌张看年周期振幅一致性更可靠。自洽性检验则是改变滤波半径如 300 km vs 400 km看流域平均时间序列的差异是否落在误差范围内差异太大说明结果对参数过于敏感需要重新考虑滤波策略。5.3 常见错误清单与排查建议以下问题在 GRACE Matlab 处理中反复出现优先级从高到低排列忘记对legendre输出做归一化转换导致振幅放大或缩小检查手段是看全球 EWh 图像是否出现南北极对称的“香蕉形”假信号。高斯权重未正确映射到阶次导致滤波效果不明显直接检查滤波后 40 阶以上系数能量是否显著下降。流域平均权重忘记cos(lat)高纬结果偏小对比官方 Mascon 流域序列即可发现。C20 替换遗漏或重复RL06 数据已含 SLR C20不要再自行替换重复替换会引入 0.5 cm 级别的趋势误差。NaN 处理不当海洋格网通常在掩膜时被设定为 NaN但流域平均时候选区域若与掩膜有缝隙会出现个别月份平均值为 NaN 的断层需要在聚合前判断有效格网数量并输出标记。5.4 输出成果的存档结构建议一次完整的 GRACE 解算建议落盘三个文件原始系数裁剪后的.mat文件、全球 EWh 月度网格、流域平均时间序列 CSV。.mat里保存lm/C/S和参数快照滤波半径、截断阶数、数据源方便复现。CSV 里加一列“有效格网数”用于诊断不良月份。这样整理后后续做趋势分析、季节分解、极端事件识别都只需要加载时间序列文件不必重新反演。本文还有配套的精品资源点击获取
Matlab GRACE水储量反演：从球谐系数到等效水高的完整流程
简介这套代码专注于处理重力恢复与气候实验卫星的观测数据面向地球物理、水文及相关专业的学生和科研人员用于将重力观测转化为全球水储量变化信息支持地下水消耗、冰川融化、干旱洪水等研究。压缩包包含10个文件总大小仅655KB其中有5份Matlab程序脚本、2份过程示意图、1份算法说明文档和1个来源说明文件覆盖数据读取、重力场展开、扰动计算、水储量反演等环节。该工具包已有2799人学习下载适合需要快速完成GRACE数据处理实验的入门者。根据内容预览主控脚本能够串联球谐系数修正、大地水准面解算与总水储量快速估计的完整流程而随包文档和截图对滤波、去噪、坐标转换等关键细节做了可视化解释可帮助读者理解每一步的物理含义并在此基础上适配自己的研究区域或数据集。整体小而完整兼具教学演示与科研参考价值。1. 用 Matlab 解算 GRACE 水储量的第一道坎GRACE 重力卫星已经退役但它的数据仍然是全球陆地水储量变化监测的黄金标准。做水文、冻土、干旱监测的人拿到 Level-2 球谐系数之后第一步往往不是科学问题而是卡在“怎么把几十 MB 的 GSM 文件变成一条时间序列”上。这个过程的门槛在于球谐系数不能直接用必须先做去相关滤波和高斯平滑再把位系数转成等效水高EWh单位、阶数截断、纬度加权、泄漏误差修正任何一步错了结果都会漂移得没法看。常见做法是用 CSR、GFZ 或 JPL 发布的 Level-2 RL06 数据配合 DDK 滤波或高斯平滑在 Matlab 里逐月处理。新手能跟步骤走熟手则会关心滤波半径选多大、截断阶数取多少、泄漏修正有没有必要。这套流程大概是读 GSM 文件、处理 C 和 S 阶系数、做滤波、球谐合成、换算等效水高、再按流域或格网统计。写成代码并不长但参数背后全是物理。2. GRACE 数据产品选型与预处理CSR、GFZ、JPL 怎么选GSM 文件格式怎么读2.1 机构产品差异与 RL06 版本选择GRACE 数据由三家机构官方发布CSR美国德州大学、GFZ德国地学中心、JPL美国喷气推进实验室。三者 Level-2 产品处理细节不同但球谐系数格式一致都遵循 RL06 标准。RL06 相比 RL05 最大的改进是背景模型替换为 AOD1B RL06同时更新了 C20 和 C30 系数。日常解算水储量我一般选 CSR 或 GFZ 的 RL06 数据二者噪声水平接近JPL 的格式多一列 sigma 误差处理时读取方式略有差异初学者先用 CSR 就好。三家机构的产品在文件命名上有统一格式以CSR_2015_01_0001_GSM_0002为例其中的GSM表示重力场模型月度解0002是处理版本号。RL06 的 C20 系数已被卫星激光测距SLR结果替代文件内会带FLAG字段标记处理时 C20 和 C30 通常直接采用文件内数值不用额外修正。2.2 GSM 文件头与块结构的 Matlab 读取方法GSM 文件是纯文本格式读取时需要先解析头部再读取 4 个数据块GRCOEF重力场系数、GRCOFS系数标准差、GRACEF加速度计标记、CALIBF校准参数。GRCOEF块内的每一行格式为l m C_lm S_lm系数值是无量纲的球谐系数单位是标准重力场模型使用的完全归一化形式。下面给出一个通用读取函数function [lm, C, S] read_gsm(filename) fid fopen(filename, r); % 逐行查找 GRCOEF 块起始位置 header textscan(fid, %s, 1, delimiter, \n); while ~strcmp(header{1}{1}(1:min(6,end)), GRCOEF) header textscan(fid, %s, 1, delimiter, \n); if feof(fid), error(No GRCOEF block found); end end % 读取系数直到下一个块标记或文件结束 lm []; C []; S []; while ~feof(fid) line fgetl(fid); if isempty(line), continue; end if ~strcmp(line(1:min(6,end)), GRCOFS) data sscanf(line, %d %d %f %f); if length(data) 4 lm(end1,:) data(1:2); C(end1) data(3); S(end1) data(4); end else break; end end fclose(fid); end这段代码的核心逻辑是先按行扫描定位GRCOEF块再从块起始位置连续读取四列数值行。sscanf的格式字符串写明了l m C S的顺序Matlab 会跳过空行和注释行。注意line(1:min(6,end))的写法是为了兼容行尾带空格的场景避免字符串索引越界。2.3 阶数截断与 C20 替换的预处理约定读入系数后第一件事不是滤波而是做两个决定截断到多少阶、是否替换 C20。GRACE 高阶系数噪声大常规做法是截断到 60 阶对应空间分辨率约 330 公里保守一点用 40 阶。截断操作直接丢弃 C/S 数组中阶数大于 60 的行即可在构造合成矩阵时也能省内存。C20 替换也是一个可选项RL06 数据里的 C20 已经是 SLR 解直接使用是合理的但如果你拿的是旧版 RL05 数据则必须用 SLR 提供的 C20 时间序列替换。预处理还有一个容易被忽略的细节GRACE 反演的 C00 项恒等于 1C10、C11、S11 在地心参考系下被强制归零。这是因为卫星重力无法独立测定地球质心运动这些项不包含地质信息。处理时不需要特殊移除这些项合成时自然被球谐函数积分掉。3. GRACE 水储量反演的核心公式与两类滤波从球谐系数到等效水高3.1 球谐合成公式为什么位系数要先乘负荷勒夫数等效水高Equivalent Water Height, EWh是全球质量重分布直接换算成水层厚度的结果。球谐域中每个阶次 (l,m) 的位系数变化量 ΔC_lm、ΔS_lm 通过如下公式映射为等效水高变化ΔEWh(θ,λ) (a·ρ_avg)/(3·ρ_w) · Σ_{l0}^{L} Σ_{m0}^{l} (2l1)/(1k_l) · [ ΔC_lm·cos(mλ) ΔS_lm·sin(mλ) ] · P̄_lm(cosθ)其中a是地球平均半径6378.1363 kmρ_avg是地球平均密度5517 kg/m³ρ_w是水密度1000 kg/m³k_l是负荷勒夫数load Love numberP̄_lm是完全归一化连带勒让德函数。乘(2l1)/(1k_l)这一步把“重力场变化”转成了“地表质量变化”是物理上最关键的一步。负荷勒夫数序列是固定常量CSR 官方的 RL06 处理文档中给出了 0-200 阶的表用load_love_numbers.m读入即可。3.2 高斯平滑与扇形滤波去相关去不掉的交给空间滤波球谐域内高阶系数的条带误差会让制图结果出现南北向的条纹经典的解决办法有两类。一类是去相关滤波比如 DDK 1/2/3 系列它利用系数的协方差信息做平滑不改变物理分辨率另一类是空间平滑最常见的是高斯滤波其本质是对球谐系数乘一个阶相关的权重因子W_l。两者的选择依据是目标流域面积大流域用 350-500 km 半径的高斯滤波中小流域用 200-300 km过度滤波会把真实信号也削掉。Matlab 里做高斯平滑不需要显式计算权重因子再逐阶相乘可以用gauss_smooth.m这类现成函数输入半径后输出一个与阶数等长的向量W_l。官方 DDK 滤波直接提供每个月的滤波后系数文件以DDK3命名读入后跳过滤波步骤只做高斯平滑即可。两条路线最终殊途同归去相关除掉条带空间平滑再压一次噪声。3.3 等效水高合成的快速实现用网格求值避开 legendre 循环球谐合成最直观的实现是双层循环逐格网点累加P_lm但这种方法在 1°×1° 网格和 60 阶截断下需要计算 64800×3661 次连带勒让德函数Matlab 循环跑起来可能要几分钟。更快的方式是先用legendre_array.m部分开源包提供例如 MATLAB GRACE Toolbox 中的函数一次性生成所有阶次在某一纬度上的 P_lm 值再利用矩阵乘法逐纬度合成。一个折中方案是在 Matlab 中调用legendre函数按纬度循环每纬度对的 P_lm 是向量仍比逐格点快两个数量级。实际项目中我常用以下实现function [lon, lat, ewh] sph_synth(lm, C, S, Lmax, R_filter) % 生成等间距经纬度网格默认 1 度 lat 89.5:-1:-89.5; lon 0.5:1:359.5; [LON, LAT] meshgrid(lon, lat); % 高斯滤波权重向量W_l 由半径 R_filterkm确定 W gauss_weights(Lmax, R_filter); % 预分配输出矩阵 ewh zeros(size(LON)); % 勒让德归一化常数映射表后续循环中重复使用 norm_factor legendre_norm(lm); for i 1:length(lat) theta (90 - lat(i)) * pi/180; P legendre(Lmax, cos(theta), sch); % 连带勒让德4 列 sum_val 0; for l 0:Lmax % 提取第 l 阶所有 m 的 P_lm P_lm P(l1, 1:l1); Clm C(lm(:,1)l, :); % 取该阶系数 % 合成C 项与 cos、S 项与 sin 结合 cos_term Clm(:,3) .* cos((0:l) .* LON(i,:)); sin_term Clm(:,4) .* sin((0:l) .* LON(i,:)); sum_val sum_val W(l1) * (2*l1)/(1love(l1)) ... * (sum(cos_term) sum(sin_term)); end ewh(i,:) sum_val * (6378.1363 * 5517 / 3000); end end上述代码里的legendre函数采用sch输出 Schmidt 归一化形式而 GRACE 标准用完全归一化两者差一个sqrt(2)因子需要注意。legendre_norm和love是预先从文件加载的表避免循环内重复计算。代码的核心思路是纬度循环里先用legendre获得该纬度的全部阶次值再在阶次循环中组合 C/S 系数和三角基。W(l1)是高斯权重(2*l1)/(1love(l1))是负荷 LOVE 数因子最后一个乘法里的5517/3000是ρ_avg/(3×ρ_w)化简结果。3.4 泄漏误差与信号恢复要不要做尺度因子修正反演出来的 EWh 经滤波后信号幅度会被削弱 30%-60%尤其是局部水储量变化大的区域如青藏高原、亚马逊流域。针对流域平均时间序列标准处理是做一个“尺度因子”修正对每个格网或流域先假设一个真实信号的时变模型比如 GLDAS 水文模型用同样的滤波处理得到“滤波后”的信号再用最小二乘拟合出比例系数k最后将观测值除以k。这个系数对强信号区域在 1.5-2.0 之间对弱信号区接近 1。对于单格网时间序列不推荐强行做尺度因子修正因为信噪比太低拟合出的系数不稳定。更准确的做法是按流域聚合后再修正或者直接与 GLDAS、GLDAS-2.1 的地表水储量输出做对比验证。泄漏误差除了幅度衰减还包括相邻区域信号的串扰比如海洋信号泄漏到沿海陆地这个目前没有一劳永逸的修正方案减少方式是选足够大的研究区。4. 用 Matlab 跑通 GRACE 水储量时间序列的完整流程从原始 GSM 到流域平均4.1 全流程代码数据下载组织、逐月反演与结果存储在动手算之前先约定目录结构。我习惯把不同月份的 GSM 文件统一放在raw/目录命名含年份月份方便批量读取。全流程代码可以拆为三个阶段预处理、滤波与合成、流域平均。下面的主脚本展示完整调度逻辑% GRACE 水储量解算主流程 % 目录raw/*.txt 存放 CSR RL06 GSM 文件 fnames dir(raw/GSM_*.txt); R_filter 300; % 高斯滤波半径单位 km Lmax 60; % 截断阶数 % 预加载负荷勒夫数与高斯权重避免循环内重复计算 love load(love_numbers_lmax200.txt); % 两列阶数k_l W gauss_weights(Lmax, R_filter); % 存放所有月份的 EWh 网格 ewh_all []; dates []; for i 1:length(fnames) % 从文件名解析年月 tok regexp(fnames(i).name, (\d{4})_(\d{2}), tokens); year str2double(tok{1}{1}); month str2double(tok{1}{2}); % 读取球谐系数 [lm, C, S] read_gsm([raw/ fnames(i).name]); % 截断到 Lmax idx lm(:,1) Lmax; lm lm(idx,:); C C(idx); S S(idx); % 计算等效水高 [lon, lat, ewh] sph_synth(lm, C, S, Lmax, R_filter, W, love); % 存入数组后续计算时间序列 ewh_all cat(3, ewh_all, ewh); dates [dates, datetime(year, month, 15)]; end % 保存结果 save(grace_ewh_global_1deg_300km.mat, lon, lat, ewh_all, dates);这段主脚本的调度很清楚先读取文件列表再对每个文件完成“读取-截断-合成-存储”四步。注意gauss_weights里的返回值需要显式传给合成函数避免每层循环重复算高斯因子。文件名的正则解析用regexp提取年月datetime(year, month, 15)选了每月 15 日作为该月时间戳方便后续画时间序列图。4.2 流域平均用经纬度边界框或 Shapefile 掩膜聚合拿到全球 EWh 网格后流域平均是水文应用最常见的一步。最简单的方法是用经纬度矩形框选目标区域适合流域形状接近矩形的区域更精确的方法是用 Shapefile 做掩膜再输出掩膜内的平均 EWh。下面给出掩膜聚合的代码% 流域平均以 Shapefile 区域的格网掩膜为例 shp shaperead(basin_boundary.shp); mask inpolygon(lon, lat, shp.X, shp.Y); % 计算掩膜面积权重纬度余弦加权 area_w cosd(lat); area_w_masked area_w .* mask; for i 1:size(ewh_all, 3) ewh_month ewh_all(:,:,i); % 避免 NaN 区域影响平均 valid ~isnan(ewh_month) mask; basin_ewh(i) sum(ewh_month(valid) .* area_w_masked(valid)) / ... sum(area_w_masked(valid)); end % 时间序列绘图 plot(dates, basin_ewh, LineWidth, 1.5); xlabel(Date); ylabel(EWh (cm)); grid on;这个加权平均逻辑很关键直接用mean会把高纬度格网权重放大因为 1° 格网的物理面积随纬度余弦变化真实流域平均必须用cos(lat)作为权重。inpolygon是基于经纬度平面的多边形判断天然适合 GRACE 1° 网格尺度。这一步输出的时间序列单位是厘米水柱与降水量、蒸散量对比时要注意量纲统一1 cm EWh 10 kg/m²。4.3 时间序列后处理季节性信号分离与长期趋势提取流域平均序列里通常包含年周期、半年周期和长期趋势。用最小二乘拟合提取趋势时常见回归模型中包含趋势项与年/半年谐波项方程如下EWh(t) β_0 β_1·t A_a·cos(2πt) B_a·sin(2πt) A_s·cos(4πt) B_s·sin(4πt) ε拟合结果是β_1即长期趋势单位通常转换为 cm/yrsqrt(A_a²B_a²)是年振幅。要在 Matlab 中实现使用设计矩阵和\运算符即可核心代码如下% 构建设计矩阵趋势 年周期 半年周期 t year(dates) (month(dates)-0.5)/12; X [ones(length(t),1), t, cos(2*pi*t), sin(2*pi*t), ... cos(4*pi*t), sin(4*pi*t)]; % 最小二乘求解 beta X \ basin_ewh; trend_mm_per_yr beta(2) * 1000; % 把 cm/yr 转成 mm/yr % 去季节信号后画残差序列 seasonal X(:,3:6) * beta(3:6); detrended basin_ewh - X(:,1:2)*beta(1:2) - seasonal;设计矩阵的构建方式要考虑时间变量的基准t用小数年表示并以年中点采样避免相位偏差。X \ basin_ewh在 Matlab 中是左除等效于最小二乘解不需要显式求逆。如果时间序列较短少于 3 年年周期的拟合不稳定趋势结果要谨慎解读优先看相对变化而不是绝对速率。5. 参数调优与验证技巧滤波半径、截断阶数与结果的可信度检查5.1 滤波半径与截断阶数的组合选择表GRACE 处理中最影响结果质量的两个参数是高斯平滑半径和截断阶数。两者的作用是叠加的截断阶数决定最短波长的下限L60 对应波长约 330 km高斯半径进一步衰减高阶能量。实践中常用组合如下应用场景截断阶数 L高斯半径 (km)适合的流域尺度全球尺度制图60300 20 万 km²区域干旱监测402505-20 万 km²大流域如亚马逊40500 100 万 km²中小流域研究601501-5 万 km²注意表格中“流域尺度”是保守估计实际可反演的最小流域与信号强度高度相关强信号区域如华北平原地下水开采区3 万 km² 也可能勉强分辨弱信号区 10 万 km² 也不一定稳定。选择参数时不要机械照搬先用研究区 3-4 个月的合成 EWh 图做目视检查条纹是否明显、信号是否连续。5.2 三种快速验证方法与官方 Mascon 对比、与 GLDAS 对比、自洽性检验日常验证中我最常用的是与 CSR/GFZ 官方发布的 Mascon 产品做对比。Mascon 产品是块体质量异常解空间分辨率约 1°但它的滤波处理与球谐法不同不会出现条带。对比方法是将球谐法结果重采样到 Mascon 网格做流域平均后看两条时间序列的相关性和均方根误差相关系数大于 0.8 且 RMSD 小于 3 cm 属于正常范围。另一个验证途径是与 GLDAS 水文模型对比重点看季节振幅和相位。GLDAS 反演的是土壤水、雪水当量、植被冠层水的总和与 GRACE 反演的总水储量差一个地表水项因此差值过大不必慌张看年周期振幅一致性更可靠。自洽性检验则是改变滤波半径如 300 km vs 400 km看流域平均时间序列的差异是否落在误差范围内差异太大说明结果对参数过于敏感需要重新考虑滤波策略。5.3 常见错误清单与排查建议以下问题在 GRACE Matlab 处理中反复出现优先级从高到低排列忘记对legendre输出做归一化转换导致振幅放大或缩小检查手段是看全球 EWh 图像是否出现南北极对称的“香蕉形”假信号。高斯权重未正确映射到阶次导致滤波效果不明显直接检查滤波后 40 阶以上系数能量是否显著下降。流域平均权重忘记cos(lat)高纬结果偏小对比官方 Mascon 流域序列即可发现。C20 替换遗漏或重复RL06 数据已含 SLR C20不要再自行替换重复替换会引入 0.5 cm 级别的趋势误差。NaN 处理不当海洋格网通常在掩膜时被设定为 NaN但流域平均时候选区域若与掩膜有缝隙会出现个别月份平均值为 NaN 的断层需要在聚合前判断有效格网数量并输出标记。5.4 输出成果的存档结构建议一次完整的 GRACE 解算建议落盘三个文件原始系数裁剪后的.mat文件、全球 EWh 月度网格、流域平均时间序列 CSV。.mat里保存lm/C/S和参数快照滤波半径、截断阶数、数据源方便复现。CSV 里加一列“有效格网数”用于诊断不良月份。这样整理后后续做趋势分析、季节分解、极端事件识别都只需要加载时间序列文件不必重新反演。本文还有配套的精品资源点击获取 
