GRACE泄露改正之全球正向建模恢复法（Forward-Modeling）

球谐系数的截断和滤波均会造成信号的泄漏（以南极区域为例）：一方面，南极陆地的冰盖质量变化信号泄漏至南极区域外，形成外泄漏； 另一方面，南极区域外非冰盖质量变化信号泄漏至南极区域内，形成内泄漏，两者统称为泄漏误差。

全球正向建模恢复法由Chen et al.(2013)提出，用于改正全球区域的泄漏误差，其改正效果与初始信号的准确性密切相关。

以CSR-Mascon为模拟数据验证全球正向建模恢复法改正南极外泄漏误差的有效性，

提前准备好：（1）CSR机构发布的Mascon产品，这里以CSR-Mascon RL06数据为例，准备了从2002年4月至2021年7月 共197个月的数据（实际试验的时候可以并不需要这么多数据）（2）全球陆地格网点的纬度加权值，即代码中的grid_weight

①M_obs为Mascon解经过球谐分析和500km高斯滤波后的结果，第一次循环开始前，假设M_tru=M_obs，每次循环的过程中，更新M_tru即可；

②循环中对M_tru进行球谐展开至60阶，再同样经过500km的高斯滤波得到M_pre，计算差值ΔM=M_obs-M_pre；

③更新M_tru=M_tru+ΔM*k (k取1.2)，本代码中以ΔM陆地加权值的大小作为迭代停止的依据；

④迭代停止，输出M_tru作为最后的泄露改正值。

clearvars -except C_mas grid_weight weight_sum
% load E:\Data\result\grid\C_mas.mat;
% load E:\Application\Code\weight.mat;
tic
ii=57;ite=30;%ii为月份数据 ite设置迭代次数
tmp=C_mas(:,:,ii);[M,N,~]=size(tmp);k=1.2;
cs_obs=gmt_grid2cs(tmp,60);
M_obs=gmt_cs2grid(cs_obs,500,0.25);%假设进行500km的高斯滤波
M_tru=M_obs;  %时变重力场的信号作为真实信号 未迭代之前M_tru=M_obs
for j=1:ite
cs_pre=gmt_grid2cs(M_tru,60);
M_pre=gmt_cs2grid(cs_pre,500,0.25);
M_delta=M_obs-M_pre;
M_sum(j,1)=sum(sum(M_delta.*grid_weight))./weight_sum;  %ΔM陆地加权值的大小作为迭代停止的依据
M_tru=M_tru+k.*M_delta;
end
cs_pre=gmt_grid2cs(M_tru,60);
M_pre=gmt_cs2grid(cs_pre,500,0.25);%假设进行500km的高斯滤波
toc

data(:,:,1)=tmp(:,:);
data(:,:,2)=M_tru(:,:);
data(:,:,3)=M_obs(:,:);
data(:,:,4)=M_pre;
data(:,:,5)=data(:,:,3)-data(:,:,4);

lon=0.125:0.25:360.125;
lat=-59.125:-0.25:-89.125;
[LON,LAT]=meshgrid(lon,lat);
figure(1)
for ii=1:5
subplot(2,3,ii)
m_proj('stereographic','lat',-90,'long',0,'radius',30);
data1(:,:)=data_anta(:,:,ii);
m_pcolor(LON,LAT,data1); %绘制假彩色图像preseucolor
shading flat;
TH1 = eval('load(''E:\Application\Code\boundary\b1.txt'')');
for ii = 1:size(TH1,1)
    if TH1(ii,1)<0
        TH1(ii,1) = 360+TH1(ii,1);
    end
end

m_grid('xtick',12,'tickdir','in','ytick', 1,'yticklabels',[],'linest','none','xaxisloc','top' );
% m_grid('xtick',12,'tickdir','out','ytick',[-90 -65.5],'linest','none','xaxisloc','top' );
%南纬65度-90度
caxis([-10,10]);
colormap('jet')
end

第一行：第一幅图为CSR-Mascon解，第二幅图为经过Forward-Modeling泄露改正后的解，第三幅图为CSR-Mascon解经球谐展开截断至60阶，并进行500km高斯滤波的结果；

第二行：第一幅图为倒数第二次迭代解M_tru经球谐展开截断至60阶，并进行500km高斯滤波的结果M_pre；第二幅图为ΔM，即第一行第三幅图与第二行第一幅图的差异。

通过比较泄露改正后的约束解和mascon经截断+高斯滤波解的差异可以发现，说明以CSR-Mascon为模拟数据足以验证全球正向建模恢复法改正南极外泄漏误差的有效性。

References

[1] 邱春洪.基于GRACE时变重力场模型反演南极冰盖质量变化[D].西南交通大学.2016.

[2] 郭飞霄. 地表物质迁移的卫星大地测量反演理论与方法研究[D].战略支援部队信息工程大学.2019.

[3] Chen J L，Wilson C R，Li J，et al. Reducing Leakage Error in GRACE-Observed Long-Term Ice Mass Change: A Case Study in West Antarctica［J］. Journal of Geodesy, 2015, 89（9）：925-940