GRACE/GRACE-FO泄漏误差数据处理之尺度因子法

尺度因子法的基本原理，直接上图

用公式表示为

通俗来理解：先用GLDAS水文模型的格网产品经过球谐分析得到球谐系数，再经过截断和高斯滤波或者去相关滤波处理得到滤波后的格网，逐格网点求解滤波前后的差异因子，最后将因子乘以GRACE/GRACE-FO经过截断和滤波后的格网数据得到尺度因子改正后的解。

代码实现：

1、先准备GLDAS（以1°×1°的Noah为例，去除04年1月至09年12月的平均值）和GRACE的数据（时间序列：2002年4月至2009年12月）组合滤波：高斯300+SWENSON去相关

具体GLDAS格网数据的生成可参照 matlab读取GLDAS数据，并绘制全球格网图-CSDN博客

clearvars -except gldas_noahr1_02_09 cs_mss_G300Swen_02_09 P;
addpath E:\Data\CSR\result\RL06;
load E:\Code\GLDAS\result\gldas_noahr1_02_09.mat
%%GLDAS_Noah 1°×1°的格网（四层土壤水+冠层含水量+雪水）
load E:\Data\CSR\result\RL06\cs_mss_G300Swen_02_09.mat;
%%02年4月至09年12月球谐系数转为质量（含缺值）
lmax=60;radius_filter=300;destrip_method='SWENSON';
cs_mass(:,:)=cs_mss_G300Swen_02_09(17,:,:); %%时间序列17对应的索引时间为2003年11月
grid_cs=gmt_cs2grid(cs_mass,0,1);

2、格网数据转换为球谐（质量），再经过组合滤波得到GLDAS滤波解

cs_gldas=gmt_grid2cs(gldas_noahr1_02_09.*0.1,lmax);
[grid_filter]=gmt_cs2grid(cs_gldas,radius_filter,1,destrip_method);toc;
temp(:,:)=grid_filter(:,:,17); 

3、 将GLDAS滤波前和滤波后的格网按时间序列排列进行最小二乘求尺度因子

for ii=1:180
    for jj=1:360
    grid_1=reshape(gldas_noahr1_02_09(ii,jj,:).*0.1,[90,1]);
    grid_2=reshape(grid_filter(ii,jj,:),[90,1]);
    k=polyfit(grid_2,grid_1,1);K(ii,jj)=k(1,1);
    end
end

结果展示：

2003年11月GRACE经过组合滤波后未进行尺度因子改正结果

2003年11月GRACE经过组合滤波后再尺度因子改正结果 

发现经过尺度因子改正之后，泄漏至海洋的信号明显减弱，但是以长江流域为例，求EWH时间序列发现，尺度因子改正的结果反而变差了（猜测是选取GLDAS水文模型的原因和尺度因子方法本身的弊端）任何问题欢迎评论或私信（q:2629512206）交流！ 多多交流  点赞支持

参考文献：

[1] 魏征强. GRACE时变重力数据球谐域滤波方法与质量分布反演应用研究[D]. 中国矿业大学, 2023. DOI:10.27623/d.cnki.gzkyu.2023.000516

[2] 姚朝龙. 联合GRACE和水文气象数据研究自然与人为因素对区域水储量变化的影响[D]. 武汉大学, 2017.

[3] Chen, J.L., Wilson, C.R., Famiglietti, J.S. et al. Attenuation effect on seasonal basin-scale water storage changes from GRACE time-variable gravity. J Geod 81, 237–245 (2007). https://doi.org/10.1007/s00190-006-0104-2.

[4] F,W,Landerer.Accuracy of scaled GRACE terrestrial water storage estimates[J].Water Resources Research, 2012.DOI:10.1029/2011WR011453.