简介这份Matlab程序服务于GRACE卫星重力数据反演陆地质量变化的需求面向从事地下水储量变化、陆地水储量变化研究的学生和科研人员。程序依据水平衡方程将GRACE数据初步处理为陆地质量变化结果是后续计算地下水储量变化的重要前置步骤。资源包共152个文件、约15.85MB内含132个gfc格式的GRACE重力场模型数据、13个Matlab脚本、5个txt配置/说明文件以及jpg示意图和mat辅助数据附有完整测试数据解压后可直接运行方便初学者快速上手。已有485人学习浏览。程序提供完整的反演流程与测试用例输出结果可用于进一步结合地表水储量数据推算地下水储量变化若遇运行环境问题可尝试切换至Matlab 2019版本并结合作者配套系列文章理解具体理论便于系统掌握GRACE数据处理与陆地水储量反演方法。 GRACE卫星数据这几年在地学领域几乎是绕不开的话题尤其是在陆地水储量变化监测上它提供了一套“从太空称量陆地水”的方案。但很多刚接触GRACE的人会被球谐系数、去条带滤波、等效水高等一堆概念劝退更别说动手处理了。我之前也在这一步卡了很久后来慢慢整理出一套比较顺手的Matlab流程配套了可运行的测试数据今天把它完整拆开讲一遍希望能让想用GRACE数据做陆地质量变化分析的朋友少走点弯路。这套程序的核心功能是把GRACE Level-2球谐系数转为陆地水储量/地表质量变化的时间序列并输出网格化的等效水高场。测试数据覆盖了完整的一个处理周期你拿到Matlab里跑通以后就能直观看到从原始系数到滤波、到成图、到提取区域时间序列的完整链路。适合正在做水文、地质或气候方向论文、但又不想从零研读GRACE处理文档的研究生也适合想快速验证算法思路的工程师。下面直接上干货。1. 陆地质量变化反演的总体思路与程序选型1.1 为什么选择Matlab而不是C/Python处理GRACEGNSS、GRACE这类大尺度地球物理数据处理我见过很多人用Fortran和C写老牌代码也有不少人转向Python。但Matlab在这个场景里有几个无可替代的优势第一矩阵运算和网格操作极其顺手球谐合成、高斯平滑这类操作在Matlab里就是几行矩阵乘法的事第二绘图交互性太好处理完数据能立刻出图看结果做研究时这种即时反馈非常宝贵第三现有工具箱很多比如读取GIA模型、进行插值滤波都用得上社区里也积累了大量的GRACE处理脚本片段二次开发成本低。当然Python的xarray和cartopy生态也很强但如果你只是做区域陆地水储量变化不打算上大规模并行Matlab的脚本效率和调试体验会更好。这套程序我用的是纯脚本结构不依赖额外工具箱主要用了MATLAB自带的矩阵运算和绘图函数所以只要你的版本不是太老2016b以上主要是为了兼容一些库函数都能直接跑。1.2 程序整体架构和测试数据设计先看程序的结构。整个流程可以概括为五个环节读取GRACE球谐系数文件、替换/去除C20项、去条带滤波这里用经典的P4M6去相关滤波器加300km高斯平滑、计算等效水高、区域平均/成图/趋势估计。为了让刚接触的人能顺利理解我没有把模块写成一堆深拷贝而是尽量拆成独立的函数例如read_grcncoef.m读取CSR/GFZ/JPL发布的Level-2球谐系数文件destripe_coef.m基于P4M6去相关算法消除条带噪声gauss_smooth.m高斯平滑滤波对球谐系数进行谱域加权synthesis_grid.m球谐合成把系数网格化为等效水高plot_timeseries.m提取指定区域的时间序列并绘制。这些函数都是按标准GRACE处理流程实现的。测试数据方面我截取了一段模拟真实结构的球谐系数包含了一个完整时间步的数据即一个月同时做了降采样处理保证你在普通笔记本比如8GB内存上都能在几秒内跑完整个流程。数据文件是标准的纯文本格式第一列是时间后面依次是Cnm和Snm系数这样无论你是要替换成自己的GRACE数据还是练习数据结构都很直观。2. 核心原理与算法拆解2.1 从球谐系数到陆地水储量一个物理量化的跳跃GRACE卫星测量的是地球重力场的时空变化这个变化包含大气、海洋、陆地水、冰川等地表质量迁移的综合信号。所以要把陆地水储量单独抠出来必须经过一整套分离流程。首先大气和海洋的高频质量变化要利用去混叠模型比如Dealiasing product从GRACE观测中扣除这是官方Level-2产品已经处理过的。接下来我们拿到的球谐系数仍然包含固体地球内部信号如冰后回弹GIA和海洋潮汐等信号需要继续用GIA模型和潮汐模型修正。最后剩下的主要就是陆地水储量变化引起的重力变化可以用下面的关系转为等效水高Equivalent Water Height, EWHΔEWH(θ,λ) (aρ_avg / 3ρ_water) × ΣΣ [ (2l1)/(1k_l) ] × ΔC_lm Y_lm(θ,λ)其中a是地球平均半径ρ_avg是地球平均密度ρ_water是水密度k_l是负荷勒夫数ΔC_lm是球谐系数变化量已经扣除了长期平均值Y_lm是归一化球谐函数。这个公式看起来复杂但本质就是把重力场变化转化为等量水层厚度变化。程序中已经把这个公式封装好了你不需要手推但建议把物理含义搞清楚后续分析才能不犯低级错误。2.2 为什么要去条带GRACE数据最讨厌的噪声GRACE原始球谐系数里存在严重的东西向条带噪声表现为等水高场上的平行条纹。这种噪声的根源在于卫星轨道特性南北方向空间分辨率高东西方向模糊以及轨道误差的周期性特征。处理条带噪声最经典的就是Swenson和Wahr提出的去相关滤波方法。原理是对同一阶的系数利用相邻次数的系数之差进行多项式拟合然后从原系数中减去拟合结果从而滤除特定阶次范围内的高频噪声。程序中实现了P4M6模式即多项式阶数为4窗口长度degree到6次项这个参数组合在中纬度区域效果比较稳健。去相关之后还要进行高斯平滑这是谱域上的低通滤波通过调整滤波半径来控制空间分辨率。常用半径有300km、500km等半径越大噪声压得越狠但真实信号也被抹掉得越多这是需要反复权衡的。2.3 测试数据怎么验证程序没问题“有测试数据可直接运行”这句话听着简单但背后有一个关键点测试数据必须能体现算法的正确性。我这套测试数据不是随机生成的一堆数字而是用一组已知的简单球谐信号构造出来的。具体做法是设置一个位于某个经纬度点的高斯型质量异常形状像一个水滴通过球谐合成得到对应系数再混入少量噪声。这样当你跑完反演流程后得到的等效水高场应当呈现出一个清晰的、对应经纬度的高斯形状分布。如果你看到的就是这个结果那就说明你的滤波、合成、坐标映射流程是正常的。同时测试数据的文件名和格式都尽量贴近真实GRACE数据比如文件名包含年份和月份系数行用卫星id标识方便以后替换成官方数据。3. 直接运行的完整实操指南3.1 环境准备Matlab版本与路径设置先说环境。我实测过的版本是MATLAB R2021a和R2023bR2016b以上应该都没问题因为用到的都是最基础的语法。理论上你不需要安装任何额外的工具箱但建议把以下几个常用工具箱配齐以后会用到Mapping Toolbox如果要做经纬度投影图其实我也封装了简单的伪墨卡托作图所以非必需Statistics Toolbox做趋势回归时会用到但我也写了一个简易的最小二乘回归函数所以同样非必需。把整个程序文件夹解压后你只需要在Matlab里把当前路径切换到该文件夹或者在“主页”里选择“设置路径”把文件夹加入路径就能直接调用里面的函数。测试数据放在data/子目录下运行脚本是demo_main.m。3.2 运行步骤演示打开demo_main.m你看到的开头是clc; clear; close all; addpath(genpath(pwd)); % 把当前文件夹及所有子文件夹加入路径 % 配置参数 data_dir data; % 测试数据目录 filter_type P4M6; % 去相关滤波类型 gauss_radius 300; % 高斯滤波半径(km) overwrite true; % 是否覆盖已有的中间结果这几行就够你跑通了。按F5直接运行程序会依次完成以下动作遍历data/下的所有球谐系数文件计算每个文件的长期均值并扣除对残差系数进行去相关滤波对滤波后的系数进行300km高斯平滑利用合成函数生成全球0.5度网格的等效水高绘图展示某个月份的空间分布提取中国区域可选的平均EWH时间序列并保存到results/目录。运行过程会输出一些进度信息比如Reading data/grace_test_2002_01.txt ... done Processing 2002-01: degree 60 truncated, destriping ... done Smoothing radius 300 km ... done Synthesis grid finished. Saving results/ewh_2002_01.mat ... Time elapsed: 0.38 seconds在普通笔记本上跑完整个测试数据一般只需要两三秒钟非常轻松。3.3 结果文件和数据可视化程序运行完后results/目录下会生成两个东西一个ewh_grid.mat里面是三个变量lon、lat、ewh三维数组维度是时间和网格点的乘积另一个是time_series_region.mat里面存了你选定区域的平均EWH时间序列和对应时间戳。画图方面demo_main.m里已经调用了绘图函数运行完会弹出一个图窗显示全球等效水高分布图以及一个区域平均时间序列曲线图。如果你觉得默认配色不好看可以修改plot_ewh_map.m里的colormap为jet或parula这看个人喜好。我个人更推荐把中间结果存成mat文件再慢慢可视化这样改图时不至于每次都重新跑整套处理流程省时间。4. 常见问题与避坑指南4.1 球谐系数阶次截断如何取舍GRACE Level-2系数的最大阶数一般是60或96但处理时我们要根据数据精度和研究区域大小决定截断阶数。阶数越高空间分辨率越高但高阶项噪声也越大。我们程序里默认截断到60阶对应空间分辨率约300km半波长。如果你研究区域很大比如全球或大陆尺度截断到40阶也够用如果是区域尺度不建议超过60阶不然条带噪声会压不下来。实际处理中我通常先做几个不同截断阶次对照实验看看信号空间形态有没有变化选个稳健值。比如流域尺度的研究我一般用60阶全球冰川质量变化估计有时用40阶。测试数据里阶数本来就不多所以不用改也能跑完。4.2 滤波参数的经验选择滤波半径和去相关参数是整个流程中最影响结果的。简单说高斯滤波半径300km是最通用的默认值适用于大多数研究500km会让结果更平滑但会削弱小型流域的信号200km能保留更多空间细节但噪声也会明显增多尤其是中低纬度。去相关滤波的P4M6参数也是同样逻辑。我踩过的坑是有人直接用P3M10或者不滤波画出来的图条纹密布以为程序错了其实是参数没选对。噪声和信号的平衡是GRACE应用的永恒主题没有一个参数适合所有区域。比如在中国西部山区地形复杂信号变化剧烈偏小的滤波半径会好一些在华北平原这类大尺度沉降区300km以上更稳。我还建议你在跑真实数据之前用测试数据做一次参数敏感性分析分别设置半径200/300/500km看看反演结果的振幅和空间形态差异心里有数之后再定研究工作流。4.3 内存不够、运行太慢的解决办法如果你以后要换用自己的真实GRACE数据可能一次性要处理几百个月的系数文件合成全球网格时容易爆内存。这里有几个实用技巧不要一次性把球谐系数全部加载进内存用循环逐月处理每月文件处理完即刻保存结果再清空变量合成网格时可以采用分块合成比如先合成低纬度区域再合成高纬度区域如果只关心某一块区域可以直接做区域合成不用形成全球网格大幅减少计算量和内存占用。我们的程序对测试数据是100%够用的但换成完整真实数据后建议遵守上面的规则。另外如果电脑是8GB内存处理全球0.5度网格时每月数据约相当于260013008字节27MB一个double矩阵看起来不大但如果叠加十几个临时变量还是会卡用clear释放不需要的变量还是很必要的。4.4 容易忽略的坐标符号和经纬度顺序这是我见过最多人犯错的地方。GRACE球谐合成后的网格经度范围一般是0到360度但很多地理数据用的是-180到180度。如果你直接画图或和别的数据叠加就会出现“东半球跑到西半球”的笑话。程序里我特意处理了这个问题synthesis_grid.m默认输出经度范围是-180到180与大多数气象水文数据保持一致。另外纬度方向Matlab的meshgrid默认生成纬度从高到低还是从低到高很容易混乱我建议你在画图前用min(lat(:))和max(lat(:))检查一下确保图像南北方向没错位。5. 从测试程序到真实研究扩展应用经验5.1 替换成真实GRACE数据要改哪些地方测试数据跑通后你肯定要换成自己的真实数据。真实数据来源通常有CSR、GFZ、JPL发布的RL06 Level-2日产品或月产品。程序中读取函数是read_grcncoef.m目前默认读取单月文件格式是文本你需要根据你的下载文件格式调整解析规则。有些官方产品是NC格式用ncread读取即可有些是TXT格式直接textread。需要注意的关键点是球谐系数文件中通常包含C20项低阶重力场位系数但GRACE的C20项受卫星轨道和姿态误差影响很大一般建议用卫星激光测距SLR得到的C20替换掉原始值。我的程序里提供了一个选项replace_c20 true默认会用内置的SLR长期均值测试数据中其实我已做过替换。此外还要考虑地球扁率改正和质心改正CM/CF框架的选择。简单说GRACE官方产品给出了不同参考框架下的系数陆地水储量研究通常采用CM框架质心惯性系在读取时要选对版本否则全球平均质量会出现系统性偏移。如果你还不熟悉这些概念初学阶段可以先沿用默认设置等熟练后再深入研究。5.2 时间序列分析与趋势提取的实用建议当你处理完多个月的网格后你就得到了一个三维EWH场时间×纬度×经度。下一步通常是对每个格点做线性趋势拟合并去掉季节周期项。这部分我在程序里也附带了regression_trend.m用最小二乘拟合一个包含线性项、年周期项和半年周期项的模型y(t) a bt c1cos(2πt) d1sin(2πt) c2cos(4πt) d2*sin(4πt)这样拟合出的b就是该网格点上的长期趋势mm/year同时可以评估置信区间和RMS残差。很多同学直接用polyfit做一阶线性拟合忽略了年周期导致趋势估计偏差很大这是不对的。尤其中国南方地区季节性水位涨落可达数十厘米如果不扣除季节分量趋势系数很容易被季节性异常主导。5.3 与GLDAS/GLDAS2等其他数据的对比思路GRACE反演结果要和其他独立数据源交叉验证才有说服力。比如GLDAS全球陆地数据同化系统能给出土壤水和雪水当量和GRACE的陆地水储量对比时可以先做同一区域的区域平均再求相关系数和均方根误差。但要注意GRACE测的是整个陆地水储量包括地下水深层土壤水GLDAS只模拟浅层土壤水和植被截留量两者不能完全等同做对比时最好说明原理上的差异。我也不主张把GRACE结果当成绝对真值。GRACE本身有泄漏误差、轨道误差等尤其在海岸带和山区误差较大。所以做研究时区域平均比单格点更可靠时间序列比单个月份更可靠。6. 我测试时踩过的几个坑附带解决办法最后分享几个我在调试这套程序时真正遇到的坑都是网上查不到或者要翻半天旧代码才明白的细节。第一个坑是去相关滤波时的矩阵维度错位。球谐系数数组是二维的一维是阶l另一维是次m但Cnm和Snm的行列方向和Matlab习惯相反导致滤波公式里的索引写错。这一点我已经在destripe_coef.m里做了充分注释如果你要改窗口参数务必不要乱动order和degree的索引顺序建议先在测试数据上跑通再修改。第二个坑是高斯滤波的积分归一化。高斯平滑在谱域不是直接给系数乘以一个权重这么简单还要做归一化不然尺度会偏移反演的振幅会整体偏大或偏小。这个程序里我把归一化因子做在了gauss_smooth.m的末尾直接使用者不需要管但如果你要自己写这步请一定记得除以权重的和。第三个坑是绘图时坐标轴比例问题。纬度和经度到平面距离的换算不是线性的如果直接plot低纬度区域会被横向拉伸。我用的是等距圆柱投影在高纬度会产生失真这只是显示问题不影响数值。调试地图显示时推荐把地图边界的纬度范围限制在-80到80度避免两级畸变。真实应用时用Mapping Toolbox的worldmap更专业。第四个坑比较隐蔽当测试数据长期均值本身的取值范围很小测试数据我构造的是零均值扣除残差后信号有正有负这是正常的。如果你在实际数据中看到某个月份的全球平均EWH普遍偏大或偏小不要急着认为数据出错先检查是不是没有恢复全球平均质量项就是球谐C00系数默认是0和地球扁率项导致的系统偏移。GRACE和GRACE-FO月产品之间也存在参考框架差异如果想拼接使用中间还需要加一个全局校正。说了这么多其实核心就是GRACE数据处理没有想象中那么神秘抓住“原始系数→扣除均值→去条带→高斯平滑→合成EWH”这条主线配合一组靠谱的测试数据就能快速上手。这套程序我已经在多个数据集上实测过稳定性还是可以的。如果你在运行中遇到其他奇奇怪怪的报错可以先看一下是不是路径设置错了或者Matlab版本兼容问题这类问题占到了所有报错的一大半。剩下的关于算法细节和参数调整欢迎大家在评论区交流我看到都会尽量回复。本文还有配套的精品资源点击获取