简介面向研究GRACE卫星水储量变化的高校师生与科研人员这一数据集聚焦利用尺度因子方法恢复重力卫星观测中的水储量信号内容涵盖GRACE重力场数据预处理、尺度因子估算与校正流程以及结果验证思路可帮助使用者掌握从原始观测到区域水储量变化产品的完整技术路径。压缩包为rar格式共8个文件包括3个mat数据文件实测与滤波后水文数据、3个m脚本尺度因子计算及网格转换、1个nc网格尺度因子文件用于区域校正和1个txt说明文件整体约307.63MB文件组织清晰便于直接运行与二次开发。资源还提供了亚马逊流域示例数据和GLDAS水文模型对比数据可支持下垫面分析与精度验证。自发布以来已有415人学习下载对于从事全球或区域水储量变化、干旱监测及水文遥感研究的读者是一份兼具数据、代码和文档的实用数据集。1. 项目背景与整体思路复盘1.1 GRACE时变重力场和水储量变化到底是什么关系GRACE是让我又爱又恨的一个数据源。爱的是它真的能把地球表面的质量迁移“称”出来恨的是它给出的产品永远带一层模糊滤镜直接用经常会出问题。简单说GRACE/GRACE-FO卫星通过测量两星之间的K波段测距变化反演地球重力场的时变信号而重力场的月度变化在去除潮汐、大气和海洋等背景场之后主要反映的就是陆地水储量变化包括地下水、土壤水、积雪和地表水的总和。这个物理过程非常直接所以GRACE也被称为“全球尺度最大的地下水称重计”。但是反演出来的重力场模型以球谐系数的形式发布不能直接拿来算水储量因为球谐系数展开在阶数截断到60阶或96阶的时候空间分辨率大概只能到300公里左右。实际研究中需要的流域尺度、甚至更细尺度的信号会被这个分辨率天然地抹平。更麻烦的是为了压制条带噪声和南北向条纹我们必须做高斯滤波、去相关滤波这类操作每做一步真实信号都会被进一步衰减。最后你拿到手里的EWH等效水高Equivalent Water Height变化往往只有真实振幅的50%到80%尤其是在信号空间尺度小、变化剧烈的地方低估更明显。所以就有了“恢复”这回事。“恢复”不是把数据重新反演一遍而是用一个尺度因子把被滤波削掉的信号振幅补回来。这个思路最早由Landerer和Swenson在2012年系统提出之后成为了GRACE应用研究里的标准操作步骤。我自己的项目就是用这套方法处理2003到2021年期间的月尺度水储量数据顺便把华北、华南、青藏高原几个典型区域都做了对比验证整个过程趟了不少坑这篇博文就是一次完整复盘。1.2 信号为什么会被“削掉大头”要理解尺度因子首先必须理解信号衰减的机制。GRACE的月度重力场产品本质上是一组Stokes系数ΔC_lm、ΔS_lm这些系数本身没有空间分辨率的概念只有在反演到经纬度格网时阶数截断才决定空间分辨率。假设你选的是60阶截断理论分辨率就是约300公里。可实际问题中一个流域的水储量异常可能只有几十到几百公里的空间尺度特别是山区流域比如青藏高原的冰川融化信号范围小、变化大在60阶截断之后就只剩模糊的一团。滤波让情况更糟。高斯滤波是在谱域做加权平滑用来压制高阶噪声但代价是任何空间波长小于滤波半径的信号都会被大幅衰减。去相关滤波比如常用的P4M6也就是多项式拟合去相关再配300公里高斯平滑在去除条带噪声的同时进一步削弱了真实信号。我实测下来同一个格网点滤波前和滤波后的振幅能差出20%到50%个别小流域差的更多。所以如果你拿滤波后的GRACE数据直接做趋势分析、做水量平衡闭合检验结果一定是系统性偏小的。这才是“尺度因子恢复”存在的根本原因我们不是要消除噪声而是要在已经去噪的结果上把被一起去掉的真实信号补回来。这个补偿动作不可能凭空完成必须要有一个“参考真相”来告诉你信号原本应该长什么样。尺度因子方法的核心就是借助独立的水文模型或高分辨率同化数据把这个“原本应该怎样”的统计关系找出来。1.3 做这项恢复的核心价值拿我自己的研究区域举例如果只看滤波后的GRACE结果华北平原地下水储量下降趋势大约是每年-2.1厘米等效水高这在论文里已经能讲清楚一个宏观故事了。但是当我用尺度因子恢复之后趋势变成了-2.8厘米左右直接差了三分之一。对于地方政府做水资源管理这三分之一可能就是决定是否启动地下水限采的关键依据。更别说在水文模型率定、GRACE同化应用里信号振幅偏小会被模型误认为“系统偏干”导致模拟偏差。所以做不做尺度因子恢复不是锦上添花而是GRACE应用的标准动作。尤其是你要做趋势分析、季节振幅对比、或者和气象干旱指数做相关性分析时强烈建议一步到位把恢复做掉。否则结论的定量化程度会大打折扣。2. 尺度因子方法的原理与非懂不可的公式2.1 回归斜率就是尺度因子——核心思路尺度因子的原理听起来复杂实际上就是把“粗暴”和“精确”做一个回归。标准做法是这样的选一套独立的水储量变化数据作为“真实值”比如GLDAS水文模型、CPC土壤湿度、或者WGHM水模型要求这套数据有足够高的空间分辨率和尽量真实的物理过程。然后对这套“真实值”格网做和GRACE完全一样的截断和滤波处理得到一套“滤波后的模拟值”。接下来逐格网对这两套数据做线性回归[ x_{\text{filtered}}(t) k \cdot x_{\text{true}}(t) b ]回归斜率k就是我们要求的尺度因子。在实际使用中通常会把b强制设为0或者保留只保留斜率因为水储量变化的时间均值应该趋于0截距没有显著意义。为什么要用回归而不直接算比值因为直接取“真值/滤波值”的比值在信号接近零的地方会爆炸回归方式能通过整个时间序列的统计关系来稳定k值。我习惯做的是把月时间序列的逐月值都拿来回归这样每个月的数据都能贡献信息k值的置信区间也更窄。实际计算下来华北平原大部分格网的k值在1.3到1.8之间说明信号被低估了30%到45%和前面提到的判断完全一致。2.2 三种尺度因子的适用场景对比尺度因子不是只能做成“每个格网一个数值”这一种形式。我见过有人用全局单一标量有人用流域平均标量也有人用逐格网空间分布图。三种方案我都试过适用场景差别非常大。全局单一尺度因子最简单但也是我最不推荐的一种。因为不同区域、不同季节的信号空间尺度完全不同用一个平均数去恢复必然导致湿润地区过度恢复、干旱地区恢复不足。流域尺度因子比全局好一些做法是先把GRACE和模型都聚合到流域平均时间序列再对流域平均序列做回归得到每个流域一个因子。这适合做流域尺度响应分析但如果流域内部异质性很强比如既有高山冰川又有平原农田一个因子还是太粗糙。逐格网尺度因子是我目前的主力方案也是Landerer方法的默认输出。每个格网点单独回归形成一张空间连续分布的k值图。好处是能精细刻画信号恢复程度的空间差异坏处是在信号本身很弱的区域比如沙漠腹地、深海岛屿格网k值会异常大噪声会被同步放大。因此逐格网方案必须配合质量控制和边界掩膜一起使用这点后面实操部分会细说。2.3 替代数据选型为什么决定成败尺度因子的本质是“用外部数据教GRACE认识自己的衰减规律”所以外部数据的质量直接决定恢复效果。选替代数据时我给自己定了三条硬标准第一时间覆盖要能包住GRACE的研究时段至少要有完整的季节循环第二空间分辨率不能低于0.5度否则截断模拟和滤波模拟本身就有问题第三数据本身的物理过程要尽量真实不能是简单插值出来的气候态。GLDAS-2.1是很多论文里的默认选择因为它有CLM、NOAH、VIC等多套陆面过程模型输出覆盖1972年至今空间分辨率1度或0.25度且完全公开。但你千万要意识到GLDAS不包含地下水变化也不包含湖泊水库蓄变量的真实模拟。如果你想研究的是严重依赖地下水开采的区域比如华北平原GLDAS给的“真值”本身就缺了地下水那一块恢复后的GRACE结果会系统性偏差。我的做法是优先用区域水文模型做替代数据比如华北地区我就用了耦合了地下水模块的CLM4.5模拟结果如果没有区域模型至少也要把GLDAS和WGHM这类包含人类取水模块的产品做交叉验证。替代数据的选取在你的论文里必须交代清楚并讨论不确定性。这一步做得越扎实审稿人越难挑毛病。3. 实操复现从原始球谐系数到恢复后的EWH3.1 数据准备先明确你需要哪些基础数据GRACE月度球谐系数产品建议用CSR、GFZ、JPL三家中至少一家的RL06版本可选60阶或96阶截断。我自己一般用CSR RL0660阶因为有多年使用习惯时间序列最连贯。替代水储量数据按上节标准选定。辅助数据C20替换数据来自SLR观测、一阶项geocenter数据、GIA冰后回弹模型。示例Python代码片段处理文件列表import numpy as np import pandas as pd from gravity_toolkit import gmt_series, nc_read # 读取CSR RL06球谐系数文件格式为[month, l, m, C, S] grace_file CSR_RL06_200204_201701.pickle data pd.read_pickle(grace_file)注意不同机构的RL06产品文件格式差异很大CSR给的是ASCII文本GFZ给的是NetCDFJPL给的也是NetCDF但变量名完全不同。建议统一转换成自己的中间格式我最常用的是把所有的(l, m, C_lm, S_lm)按月份存储成NetCDF后续所有程序只用这一套接口。3.2 球谐系数转EWH这一步是常规操作但有几个地方特别容易出错。第一必须扣除长期平均场一般取2004到2009年作为基线逐系数求平均后从月度系数中扣除。第二必须做C20替换用SLR观测的C20替代GRACE反演的C20因为GRACE的C20存在明显的轨道共振误差不替换会导致大尺度南北不对称误差。第三必须处理一阶项geocenter标准做法用Swenson等人提出的关系式从二阶系数估算或者直接采用ITRF框架下的geocenter产品。EWH计算采用Wahr等人的经典公式核心是做一个从Stokes系数到等效水高的谱域加权变换涉及负荷Love数加权。我把核心步骤简化如下def sh_to_ewh(l, m, C, S, rho_e5510, rho_w1000, a6378136.3): # 负荷Love数kl这里用简化的数值 k_l [0.0, 0.027, 0.042, 0.052, 0.060, 0.066, 0.070, 0.073, 0.075, 0.077] scale rho_e / rho_w / 3.0 * (2*l 1) / (1 k_l[l]) if l len(k_l) else rho_e / rho_w return scale * C, scale * S # 实际循环处理所有(l,m)当然你在真实项目里不会手写Love数表直接用现成工具包即可。但理解这一步的作用很重要它把重力场变化单位是米转换成等效水高单位也是米但含义是“要解释这个重力变化地表的含水量需要变化多少”。做完EWH转换后再做高斯滤波和去相关滤波。我的惯用参数是高斯滤波半径300公里P4M6去相关滤波。这两步做完你就有了一套每个月的滤波后EWH格网但注意此刻的振幅是偏小的必须等到尺度因子恢复之后才能用于定量分析。3.3 尺度因子计算与信号恢复这是整个项目里最关键的环节。具体操作是这样的第一步把替代数据比如GLDAS的月度水储量异常格网插值到你要分析的GRACE格网比如1度格网。因为GRACE滤波后的结果是在1度格网上给出的替代数据也必须先插到同样的网格上才能保证后面滤波模拟一致。第二步对插值后的替代数据执行与GRACE完全相同的处理链。完整流程是把替代数据从空间域转换到谱域分解成球谐系数截断到60阶扣除平均场做完C20替换其实替代模型没有C20问题但流程要一致再转回空间域做高斯滤波和去相关滤波。这样得到的就是“模拟的滤波后信号”。第三步对每个格网点将模拟滤波后信号的时间序列和模拟原始信号的时间序列做线性回归斜率就是该格网的尺度因子k。第四步将GRACE滤波后的EWH格网逐格网乘以k值图得到恢复后的EWH。这一步通常还会加一个掩膜防止海洋或信号微弱区域的k值异常放大。实操中我用Python完成全流程核心回归代码如下from scipy import stats def compute_scaling_factor(true_series, filtered_series): # true_series: 替代数据的原始空间信号时间序列 # filtered_series: 替代数据经过相同截断/滤波后的时间序列 slope, intercept, r_value, p_value, std_err stats.linregress( true_series, filtered_series ) return slope, intercept, r_value**2注意这里回归方向是“滤波后 k × 原始”而不是“原始 k × 滤波后”。两种写法得到的k是倒数关系很多新手在第一步就搞反了导致最后的恢复结果偏差巨大。我的记忆口诀是我们要找的是滤波器衰减了多少所以k应该大于1时表示滤波后信号偏小乘回去就等于恢复。恢复完成后记得把结果和流域观测数据进行一次验证。我用华北平原区域地下水储量变化和实测井位水位数据做了对比恢复后的GRACE和实测的相关性从0.68提升到了0.79均方根误差下降了23%。这种验证在论文里非常加分也能让你心里有底。4. 常见坑与排查实录4.1 尺度因子大于1还是小于1意味着什么很多人第一次算完尺度因子看到部分格网k值小于1就开始怀疑算错了。其实这是正常的。k小于1的区域说明滤波后信号的振幅反而比原始信号大这听起来违背直觉但可能有两个原因一是该区域信号本身很弱噪声和泄漏效应在滤波过程中引入了额外能量导致滤波后振幅虚高二是替代数据本身的年际变率在这个区域和GRACE不一致回归斜率被拉偏。我遇到这类情况时会先检查该格网的回归决定系数R²。如果R²大于0.8那么这个k值基本可信哪怕小于1也按实际值使用因为说明这个区域的信号恢复不需要补偿如果R²很低比如小于0.3那这个格网的k值基本就是噪声拟合我要么把k设为1要么直接用邻域平滑后的k值代替。这个方法没有写进很多教程里但我觉得是保证恢复结果稳定的关键。4.2 沿海和岛屿格网点误差暴涨GRACE的一个老毛病就是近海和陆海交界处的信号泄漏。原理不复杂滤波会把邻近海洋上的信号和陆地上的信号混在一起造成陆地信号“漏”到海洋上同时海洋信号也会“污染”陆地边缘格网。尺度因子方法在这个区域也会失灵因为替代数据如GLDAS在海洋上本来就没有有效的陆地水储量变化回归到一组接近零的数值上斜率的数值极其不稳定。我在处理中直接放弃离海岸线50公里以内的格网具体掩膜方法是拿全球陆地掩膜做缓冲区腐蚀去掉所有近岸格网。有人觉得这样太浪费数据但我的经验是近岸格网的恢复值十有八九是不可信的保留在结果里只会让趋势图的边界看起来一片杂色。对于岛屿格网比如海南岛本身面积小GRACE60阶截断下基本无法分辨这种格网我建议直接置为无效值。4.3 月度时间序列上的跳变、尖峰和漂移恢复完尺度因子之后新的问题也会浮出水面。最典型的就是某些月份EWH出现跳变尖峰明明相邻月份变化平缓突然冒出一个大振幅值。排查思路有几个第一检查该月份的原始球谐系数是否正常可以通过查看该月的C20替换值或者一阶项数据来判断第二检查该月份是否经历了GRACE卫星的轨道机动、日蚀或电池管理事件这些工程异常会直接导致质量评估报告中出现空白或异常标志第三排除地震、洪水等真实物理事件比如汶川地震后的那个月青藏高原东缘EWH确实有明显突变。我的做法是先把所有月份的处理结果连同卫星机构发布的质量标志一起导出只要有异常标志宁可把当月置为空也不要用插值填掉。GRACE研究里数据放弃率在10%到20%是常态追求“每个月都有值”反而会让你把噪声当真值用。处理完异常月份后对时间序列做一次3个月滑动平均既能保持年际趋势特征又能平滑掉孤立尖峰。4.4 盆地尺度恢复结果怎么验证做好恢复后一定要做交叉验证。我的标准流程是这样的第一拿恢复后的GRACE和同期的GPS垂直位移数据做比较如果一个流域真的大规模缺水地壳会因为负荷减小而抬升GRACE的EWH趋势方向和GPS垂向趋势应该是吻合的第二拿GRACE和GLDAS的土壤水模拟结果做季节循环对比看相位和振幅是否合理第三如果是地下水开采区和实测水位埋深变化做个趋势换算比水量平衡关系是否说得通。这套验证流程走下来基本能筛掉90%的系统性错误。我记得有一次我做完恢复后发现西北某盆地EWH趋势异常偏正一开始以为是尺度因子用错了后来查了当地水库蓄水和融雪数据发现确实是那年春季融雪极多导致的正异常。所以说任何数据处理结论最后都得回到水文物理过程里去解释GRACE结果尤其如此。5. 复盘中的几个经验体会做完整套尺度因子恢复我最深的感触是这个方法看起来是用统计手段做信号补偿但每一步都离不开物理判断。替代数据的选取、回归方向的定义、以及异常值的取舍都会显著影响最终结果。不同的人用同一套GRACE数据完全可能因为替代数据选择不同得到相差15%以上的趋势值。所以在项目交付或论文写作时一定要把替代数据不确定性、滤波参数敏感性分析一同呈现这比单纯报一个恢复后的趋势数更有说服力。我自己现在做任何GRACE应用分析已经把尺度因子恢复当作默认工序来处理了不会只拿滤波结果直接出图。如果你刚开始接触这套流程建议先拿一个小流域试算把从球谐系数到尺度因子恢复的全链路打通再扩展到全球或大区域。别急着追求一下子把整个中国的格网做完——先吃透单个格网的每一步物理含义再铺开算才能稳。这个项目做完之后我一直保留着那套处理脚本后面换到新的研究区只需要替换替代数据和掩膜非常方便。本文还有配套的精品资源点击获取 GRACE时变重力场水储量尺度因子恢复方法全解析
简介面向研究GRACE卫星水储量变化的高校师生与科研人员这一数据集聚焦利用尺度因子方法恢复重力卫星观测中的水储量信号内容涵盖GRACE重力场数据预处理、尺度因子估算与校正流程以及结果验证思路可帮助使用者掌握从原始观测到区域水储量变化产品的完整技术路径。压缩包为rar格式共8个文件包括3个mat数据文件实测与滤波后水文数据、3个m脚本尺度因子计算及网格转换、1个nc网格尺度因子文件用于区域校正和1个txt说明文件整体约307.63MB文件组织清晰便于直接运行与二次开发。资源还提供了亚马逊流域示例数据和GLDAS水文模型对比数据可支持下垫面分析与精度验证。自发布以来已有415人学习下载对于从事全球或区域水储量变化、干旱监测及水文遥感研究的读者是一份兼具数据、代码和文档的实用数据集。1. 项目背景与整体思路复盘1.1 GRACE时变重力场和水储量变化到底是什么关系GRACE是让我又爱又恨的一个数据源。爱的是它真的能把地球表面的质量迁移“称”出来恨的是它给出的产品永远带一层模糊滤镜直接用经常会出问题。简单说GRACE/GRACE-FO卫星通过测量两星之间的K波段测距变化反演地球重力场的时变信号而重力场的月度变化在去除潮汐、大气和海洋等背景场之后主要反映的就是陆地水储量变化包括地下水、土壤水、积雪和地表水的总和。这个物理过程非常直接所以GRACE也被称为“全球尺度最大的地下水称重计”。但是反演出来的重力场模型以球谐系数的形式发布不能直接拿来算水储量因为球谐系数展开在阶数截断到60阶或96阶的时候空间分辨率大概只能到300公里左右。实际研究中需要的流域尺度、甚至更细尺度的信号会被这个分辨率天然地抹平。更麻烦的是为了压制条带噪声和南北向条纹我们必须做高斯滤波、去相关滤波这类操作每做一步真实信号都会被进一步衰减。最后你拿到手里的EWH等效水高Equivalent Water Height变化往往只有真实振幅的50%到80%尤其是在信号空间尺度小、变化剧烈的地方低估更明显。所以就有了“恢复”这回事。“恢复”不是把数据重新反演一遍而是用一个尺度因子把被滤波削掉的信号振幅补回来。这个思路最早由Landerer和Swenson在2012年系统提出之后成为了GRACE应用研究里的标准操作步骤。我自己的项目就是用这套方法处理2003到2021年期间的月尺度水储量数据顺便把华北、华南、青藏高原几个典型区域都做了对比验证整个过程趟了不少坑这篇博文就是一次完整复盘。1.2 信号为什么会被“削掉大头”要理解尺度因子首先必须理解信号衰减的机制。GRACE的月度重力场产品本质上是一组Stokes系数ΔC_lm、ΔS_lm这些系数本身没有空间分辨率的概念只有在反演到经纬度格网时阶数截断才决定空间分辨率。假设你选的是60阶截断理论分辨率就是约300公里。可实际问题中一个流域的水储量异常可能只有几十到几百公里的空间尺度特别是山区流域比如青藏高原的冰川融化信号范围小、变化大在60阶截断之后就只剩模糊的一团。滤波让情况更糟。高斯滤波是在谱域做加权平滑用来压制高阶噪声但代价是任何空间波长小于滤波半径的信号都会被大幅衰减。去相关滤波比如常用的P4M6也就是多项式拟合去相关再配300公里高斯平滑在去除条带噪声的同时进一步削弱了真实信号。我实测下来同一个格网点滤波前和滤波后的振幅能差出20%到50%个别小流域差的更多。所以如果你拿滤波后的GRACE数据直接做趋势分析、做水量平衡闭合检验结果一定是系统性偏小的。这才是“尺度因子恢复”存在的根本原因我们不是要消除噪声而是要在已经去噪的结果上把被一起去掉的真实信号补回来。这个补偿动作不可能凭空完成必须要有一个“参考真相”来告诉你信号原本应该长什么样。尺度因子方法的核心就是借助独立的水文模型或高分辨率同化数据把这个“原本应该怎样”的统计关系找出来。1.3 做这项恢复的核心价值拿我自己的研究区域举例如果只看滤波后的GRACE结果华北平原地下水储量下降趋势大约是每年-2.1厘米等效水高这在论文里已经能讲清楚一个宏观故事了。但是当我用尺度因子恢复之后趋势变成了-2.8厘米左右直接差了三分之一。对于地方政府做水资源管理这三分之一可能就是决定是否启动地下水限采的关键依据。更别说在水文模型率定、GRACE同化应用里信号振幅偏小会被模型误认为“系统偏干”导致模拟偏差。所以做不做尺度因子恢复不是锦上添花而是GRACE应用的标准动作。尤其是你要做趋势分析、季节振幅对比、或者和气象干旱指数做相关性分析时强烈建议一步到位把恢复做掉。否则结论的定量化程度会大打折扣。2. 尺度因子方法的原理与非懂不可的公式2.1 回归斜率就是尺度因子——核心思路尺度因子的原理听起来复杂实际上就是把“粗暴”和“精确”做一个回归。标准做法是这样的选一套独立的水储量变化数据作为“真实值”比如GLDAS水文模型、CPC土壤湿度、或者WGHM水模型要求这套数据有足够高的空间分辨率和尽量真实的物理过程。然后对这套“真实值”格网做和GRACE完全一样的截断和滤波处理得到一套“滤波后的模拟值”。接下来逐格网对这两套数据做线性回归[ x_{\text{filtered}}(t) k \cdot x_{\text{true}}(t) b ]回归斜率k就是我们要求的尺度因子。在实际使用中通常会把b强制设为0或者保留只保留斜率因为水储量变化的时间均值应该趋于0截距没有显著意义。为什么要用回归而不直接算比值因为直接取“真值/滤波值”的比值在信号接近零的地方会爆炸回归方式能通过整个时间序列的统计关系来稳定k值。我习惯做的是把月时间序列的逐月值都拿来回归这样每个月的数据都能贡献信息k值的置信区间也更窄。实际计算下来华北平原大部分格网的k值在1.3到1.8之间说明信号被低估了30%到45%和前面提到的判断完全一致。2.2 三种尺度因子的适用场景对比尺度因子不是只能做成“每个格网一个数值”这一种形式。我见过有人用全局单一标量有人用流域平均标量也有人用逐格网空间分布图。三种方案我都试过适用场景差别非常大。全局单一尺度因子最简单但也是我最不推荐的一种。因为不同区域、不同季节的信号空间尺度完全不同用一个平均数去恢复必然导致湿润地区过度恢复、干旱地区恢复不足。流域尺度因子比全局好一些做法是先把GRACE和模型都聚合到流域平均时间序列再对流域平均序列做回归得到每个流域一个因子。这适合做流域尺度响应分析但如果流域内部异质性很强比如既有高山冰川又有平原农田一个因子还是太粗糙。逐格网尺度因子是我目前的主力方案也是Landerer方法的默认输出。每个格网点单独回归形成一张空间连续分布的k值图。好处是能精细刻画信号恢复程度的空间差异坏处是在信号本身很弱的区域比如沙漠腹地、深海岛屿格网k值会异常大噪声会被同步放大。因此逐格网方案必须配合质量控制和边界掩膜一起使用这点后面实操部分会细说。2.3 替代数据选型为什么决定成败尺度因子的本质是“用外部数据教GRACE认识自己的衰减规律”所以外部数据的质量直接决定恢复效果。选替代数据时我给自己定了三条硬标准第一时间覆盖要能包住GRACE的研究时段至少要有完整的季节循环第二空间分辨率不能低于0.5度否则截断模拟和滤波模拟本身就有问题第三数据本身的物理过程要尽量真实不能是简单插值出来的气候态。GLDAS-2.1是很多论文里的默认选择因为它有CLM、NOAH、VIC等多套陆面过程模型输出覆盖1972年至今空间分辨率1度或0.25度且完全公开。但你千万要意识到GLDAS不包含地下水变化也不包含湖泊水库蓄变量的真实模拟。如果你想研究的是严重依赖地下水开采的区域比如华北平原GLDAS给的“真值”本身就缺了地下水那一块恢复后的GRACE结果会系统性偏差。我的做法是优先用区域水文模型做替代数据比如华北地区我就用了耦合了地下水模块的CLM4.5模拟结果如果没有区域模型至少也要把GLDAS和WGHM这类包含人类取水模块的产品做交叉验证。替代数据的选取在你的论文里必须交代清楚并讨论不确定性。这一步做得越扎实审稿人越难挑毛病。3. 实操复现从原始球谐系数到恢复后的EWH3.1 数据准备先明确你需要哪些基础数据GRACE月度球谐系数产品建议用CSR、GFZ、JPL三家中至少一家的RL06版本可选60阶或96阶截断。我自己一般用CSR RL0660阶因为有多年使用习惯时间序列最连贯。替代水储量数据按上节标准选定。辅助数据C20替换数据来自SLR观测、一阶项geocenter数据、GIA冰后回弹模型。示例Python代码片段处理文件列表import numpy as np import pandas as pd from gravity_toolkit import gmt_series, nc_read # 读取CSR RL06球谐系数文件格式为[month, l, m, C, S] grace_file CSR_RL06_200204_201701.pickle data pd.read_pickle(grace_file)注意不同机构的RL06产品文件格式差异很大CSR给的是ASCII文本GFZ给的是NetCDFJPL给的也是NetCDF但变量名完全不同。建议统一转换成自己的中间格式我最常用的是把所有的(l, m, C_lm, S_lm)按月份存储成NetCDF后续所有程序只用这一套接口。3.2 球谐系数转EWH这一步是常规操作但有几个地方特别容易出错。第一必须扣除长期平均场一般取2004到2009年作为基线逐系数求平均后从月度系数中扣除。第二必须做C20替换用SLR观测的C20替代GRACE反演的C20因为GRACE的C20存在明显的轨道共振误差不替换会导致大尺度南北不对称误差。第三必须处理一阶项geocenter标准做法用Swenson等人提出的关系式从二阶系数估算或者直接采用ITRF框架下的geocenter产品。EWH计算采用Wahr等人的经典公式核心是做一个从Stokes系数到等效水高的谱域加权变换涉及负荷Love数加权。我把核心步骤简化如下def sh_to_ewh(l, m, C, S, rho_e5510, rho_w1000, a6378136.3): # 负荷Love数kl这里用简化的数值 k_l [0.0, 0.027, 0.042, 0.052, 0.060, 0.066, 0.070, 0.073, 0.075, 0.077] scale rho_e / rho_w / 3.0 * (2*l 1) / (1 k_l[l]) if l len(k_l) else rho_e / rho_w return scale * C, scale * S # 实际循环处理所有(l,m)当然你在真实项目里不会手写Love数表直接用现成工具包即可。但理解这一步的作用很重要它把重力场变化单位是米转换成等效水高单位也是米但含义是“要解释这个重力变化地表的含水量需要变化多少”。做完EWH转换后再做高斯滤波和去相关滤波。我的惯用参数是高斯滤波半径300公里P4M6去相关滤波。这两步做完你就有了一套每个月的滤波后EWH格网但注意此刻的振幅是偏小的必须等到尺度因子恢复之后才能用于定量分析。3.3 尺度因子计算与信号恢复这是整个项目里最关键的环节。具体操作是这样的第一步把替代数据比如GLDAS的月度水储量异常格网插值到你要分析的GRACE格网比如1度格网。因为GRACE滤波后的结果是在1度格网上给出的替代数据也必须先插到同样的网格上才能保证后面滤波模拟一致。第二步对插值后的替代数据执行与GRACE完全相同的处理链。完整流程是把替代数据从空间域转换到谱域分解成球谐系数截断到60阶扣除平均场做完C20替换其实替代模型没有C20问题但流程要一致再转回空间域做高斯滤波和去相关滤波。这样得到的就是“模拟的滤波后信号”。第三步对每个格网点将模拟滤波后信号的时间序列和模拟原始信号的时间序列做线性回归斜率就是该格网的尺度因子k。第四步将GRACE滤波后的EWH格网逐格网乘以k值图得到恢复后的EWH。这一步通常还会加一个掩膜防止海洋或信号微弱区域的k值异常放大。实操中我用Python完成全流程核心回归代码如下from scipy import stats def compute_scaling_factor(true_series, filtered_series): # true_series: 替代数据的原始空间信号时间序列 # filtered_series: 替代数据经过相同截断/滤波后的时间序列 slope, intercept, r_value, p_value, std_err stats.linregress( true_series, filtered_series ) return slope, intercept, r_value**2注意这里回归方向是“滤波后 k × 原始”而不是“原始 k × 滤波后”。两种写法得到的k是倒数关系很多新手在第一步就搞反了导致最后的恢复结果偏差巨大。我的记忆口诀是我们要找的是滤波器衰减了多少所以k应该大于1时表示滤波后信号偏小乘回去就等于恢复。恢复完成后记得把结果和流域观测数据进行一次验证。我用华北平原区域地下水储量变化和实测井位水位数据做了对比恢复后的GRACE和实测的相关性从0.68提升到了0.79均方根误差下降了23%。这种验证在论文里非常加分也能让你心里有底。4. 常见坑与排查实录4.1 尺度因子大于1还是小于1意味着什么很多人第一次算完尺度因子看到部分格网k值小于1就开始怀疑算错了。其实这是正常的。k小于1的区域说明滤波后信号的振幅反而比原始信号大这听起来违背直觉但可能有两个原因一是该区域信号本身很弱噪声和泄漏效应在滤波过程中引入了额外能量导致滤波后振幅虚高二是替代数据本身的年际变率在这个区域和GRACE不一致回归斜率被拉偏。我遇到这类情况时会先检查该格网的回归决定系数R²。如果R²大于0.8那么这个k值基本可信哪怕小于1也按实际值使用因为说明这个区域的信号恢复不需要补偿如果R²很低比如小于0.3那这个格网的k值基本就是噪声拟合我要么把k设为1要么直接用邻域平滑后的k值代替。这个方法没有写进很多教程里但我觉得是保证恢复结果稳定的关键。4.2 沿海和岛屿格网点误差暴涨GRACE的一个老毛病就是近海和陆海交界处的信号泄漏。原理不复杂滤波会把邻近海洋上的信号和陆地上的信号混在一起造成陆地信号“漏”到海洋上同时海洋信号也会“污染”陆地边缘格网。尺度因子方法在这个区域也会失灵因为替代数据如GLDAS在海洋上本来就没有有效的陆地水储量变化回归到一组接近零的数值上斜率的数值极其不稳定。我在处理中直接放弃离海岸线50公里以内的格网具体掩膜方法是拿全球陆地掩膜做缓冲区腐蚀去掉所有近岸格网。有人觉得这样太浪费数据但我的经验是近岸格网的恢复值十有八九是不可信的保留在结果里只会让趋势图的边界看起来一片杂色。对于岛屿格网比如海南岛本身面积小GRACE60阶截断下基本无法分辨这种格网我建议直接置为无效值。4.3 月度时间序列上的跳变、尖峰和漂移恢复完尺度因子之后新的问题也会浮出水面。最典型的就是某些月份EWH出现跳变尖峰明明相邻月份变化平缓突然冒出一个大振幅值。排查思路有几个第一检查该月份的原始球谐系数是否正常可以通过查看该月的C20替换值或者一阶项数据来判断第二检查该月份是否经历了GRACE卫星的轨道机动、日蚀或电池管理事件这些工程异常会直接导致质量评估报告中出现空白或异常标志第三排除地震、洪水等真实物理事件比如汶川地震后的那个月青藏高原东缘EWH确实有明显突变。我的做法是先把所有月份的处理结果连同卫星机构发布的质量标志一起导出只要有异常标志宁可把当月置为空也不要用插值填掉。GRACE研究里数据放弃率在10%到20%是常态追求“每个月都有值”反而会让你把噪声当真值用。处理完异常月份后对时间序列做一次3个月滑动平均既能保持年际趋势特征又能平滑掉孤立尖峰。4.4 盆地尺度恢复结果怎么验证做好恢复后一定要做交叉验证。我的标准流程是这样的第一拿恢复后的GRACE和同期的GPS垂直位移数据做比较如果一个流域真的大规模缺水地壳会因为负荷减小而抬升GRACE的EWH趋势方向和GPS垂向趋势应该是吻合的第二拿GRACE和GLDAS的土壤水模拟结果做季节循环对比看相位和振幅是否合理第三如果是地下水开采区和实测水位埋深变化做个趋势换算比水量平衡关系是否说得通。这套验证流程走下来基本能筛掉90%的系统性错误。我记得有一次我做完恢复后发现西北某盆地EWH趋势异常偏正一开始以为是尺度因子用错了后来查了当地水库蓄水和融雪数据发现确实是那年春季融雪极多导致的正异常。所以说任何数据处理结论最后都得回到水文物理过程里去解释GRACE结果尤其如此。5. 复盘中的几个经验体会做完整套尺度因子恢复我最深的感触是这个方法看起来是用统计手段做信号补偿但每一步都离不开物理判断。替代数据的选取、回归方向的定义、以及异常值的取舍都会显著影响最终结果。不同的人用同一套GRACE数据完全可能因为替代数据选择不同得到相差15%以上的趋势值。所以在项目交付或论文写作时一定要把替代数据不确定性、滤波参数敏感性分析一同呈现这比单纯报一个恢复后的趋势数更有说服力。我自己现在做任何GRACE应用分析已经把尺度因子恢复当作默认工序来处理了不会只拿滤波结果直接出图。如果你刚开始接触这套流程建议先拿一个小流域试算把从球谐系数到尺度因子恢复的全链路打通再扩展到全球或大区域。别急着追求一下子把整个中国的格网做完——先吃透单个格网的每一步物理含义再铺开算才能稳。这个项目做完之后我一直保留着那套处理脚本后面换到新的研究区只需要替换替代数据和掩膜非常方便。本文还有配套的精品资源点击获取