GRACE无约束解和约束解区域TWSA时间序列RMSE计算——以亚马逊流域为例

2024年最后一天了，本来不想推文了，但想了想好久没发了，今天在这里水上一篇（笑哭）。主要是自己做实验中遇到的一点点小问题，在这里做个分享。例如针对无约束解和约束解（选取DDK3滤波解），求亚马逊流域一段时间内水储量的变化，以官方机构JPL和CSR发布的mascon产品为参考，求解均方根误差（RMSE）。

无约束解（左上）、DDK3滤波解（右上）、JPL-mascon解（左下）和CSR-mascon解（右下）（2004年3月）

 1、绘制各种解在亚马逊流域的水储量变化时间序列如下：

 2、我们以JPL和CSR mascon解为参考，分别计算无约束解和DDK3滤波解时间序列的相关系数和RMSE：

解类型	JPL-mas参考的相关系数	CSR-mas参考的相关系数
无约束解	0.996343847670052	0.996327393952841
DDK3滤波解	0.998829604846770	0.998873910682446

解类型	JPL-mas参考的RMSE	CSR-mas参考的RMSE
无约束解	1.24827346862736	1.22046093875777
DDK3滤波解	0.719650617974905	0.688636641136356

问题出现：明明从全球格网图上很直观的看出，无约束解具有很严重的条带误差，但就时间序列、相关系数和RMSE可以发现，无约束解和滤波解之间的差距并不大，甚至压根体现不出来差异。

可能的原因分析：在计算亚马逊流域整体水储量变化时，存在大量的正负抵消现象，导致总水储量的变化基本不收条带误差的影响。

如果我们想体现出无约束解和滤波解之间的差异，急需从逐个格网点出发，计算每个格网点的相关系数和RMSE。核心代码如下：

grid_refer=grid_mas_region; %%选择一个mascon解作为参考  之后逐个格网点求解相关系数R和RMSE
for ii=1:size(grid_refer,1) %%num_file代表时间序列的长度
    for jj=1:size(grid_refer,2)
        r=corrcoef(reshape(grid_refer(ii,jj,:),[num_file,1]),reshape(grid_unfilter_region(ii,jj,:),[num_file,1]));
        R(ii,jj,1)=r(2,1);
        r=corrcoef(reshape(grid_refer(ii,jj,:),[num_file,1]),reshape(grid_DDK3_region(ii,jj,:),[num_file,1]));
        R(ii,jj,2)=r(2,1);
        RMSE(ii,jj,1)=sqrt(mean((reshape(grid_refer(ii,jj,:),[num_file,1])-...
            reshape(grid_unfilter_region(ii,jj,:),[num_file,1])).^2));
        RMSE(ii,jj,2)=sqrt(mean((reshape(grid_refer(ii,jj,:),[num_file,1])-...
            reshape(grid_DDK3_region(ii,jj,:),[num_file,1])).^2));       
    end
end

 下面展示以JPL-mascon为参考，逐格网点求解的RMSE和相关系数R的空间分布结果：

当我们逐个格网点求解之后，无约束解和滤波解之间的差异还是很明显的。在统计亚马逊流域内所有格网点的RMSE和相关系数。相关系数，无约束解：0.455527711340854，滤波解：0.932328001525175；RMSE，无约束解：29.7580778643129，滤波解：6.07860700210584。

总结：当我们计算整个区域的一些评价指标后发现，滤不滤波的差异不明显的时候，不妨试一试逐个格网点进行求解，说不定会有符合预期的效果。

欢迎点赞、收藏、转发支持，多多交流，我们2025一起变得更好。To be better！！！

新年快乐，诸事顺遂