尺度因子法的基本原理,直接上图
用公式表示为
通俗来理解:先用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.