WH_G_Y的博客

可以发现在主要的模态中的2004年底存在明显的阶越信号（阶跃信号在地震监测中大概率由地震引起，主要因其直接反映地壳突发性永久位移）。

到目前为止，大多数GRACE研究都集中在制作质量变化图上，很少讨论这些图中的误差。然而，如果GRACE要被用作评估和改进水文和海洋模型的工具，误差估计是必要的。此外，只有误差估计才能确定数据的某些特征是否真实，以及GRACE确定该特征的准确性。上来以示尊重，先给出经典参考文献：Wahr J, Swenson S, Velicogna I. Accuracy of GRACE mass estimates[J]. Geophysical Research Letters, 2006, 33(6).

主成分分析是一种技术，广泛用于数据分析和机器学习中。它的主要目标是通过找到数据的主成分，将高维数据投影到一个低维空间，同时尽可能保留原数据的。PCA 的核心思想是将数据从原始坐标系转换到一个新的坐标系，这个新坐标系由数据的构成。

例如针对无约束解和约束解（选取DDK3滤波解），求亚马逊流域一段时间内水储量的变化，以官方机构JPL和CSR发布的mascon产品为参考，求解均方根误差（RMSE）。：明明从全球格网图上很直观的看出，无约束解具有很严重的条带误差，但就时间序列、相关系数和RMSE可以发现，无约束解和滤波解之间的差距并不大，甚至压根体现不出来差异。：当我们计算整个区域的一些评价指标后发现，滤不滤波的差异不明显的时候，不妨试一试逐个格网点进行求解，说不定会有符合预期的效果。

本博文旨在针对仿真生成的棋盘格网，分别经过球谐分析得到截断至60、120、200阶数的结果，来清晰地展示截断对最后信号大小的影响。选择美国西部一片矩形区域（经度从西经106°-西经126°，北纬38°-北纬52°） 具体的边界经纬度坐标值存入文后的网盘链接中。是一种简单而有效的图像质量评估方法，它可以帮助用户了解图像处理和显示设备对图像质量的影响，从而选择最佳的图像处理方法和显示设备。最后，将格网转为球谐分别进行200、120和60阶截断，之后再转为区域格网绘制。

主成分分析-维基百科，自由的百科全书 --- Robust principal component analysis - Wikipedia。[CC字幕]流体动力学机器学习：11.鲁棒主成分分析 (RPCA)_哔哩哔哩_bilibili。变量L为低秩矩阵，S为稀疏矩阵，count为迭代次数。

Institute of Geodesy at Graz University of Technolog（格拉茨理工大学大地测量研究所） 2018版本，最高60阶球谐系数。代码中将所有球谐系数均人为截断至60阶，只对球谐系数做最基本的替换低阶项和扣除平均重力场处理，其余包括GIA改正、滤波处理和泄漏改正等均未涉及；：时变重力场解决方案组合服务解RL01版本，最高96阶球谐系数。：武汉大学 基于卫星间位势差的约束GRACE月重力场模型。：同济大学2022版本，最高96阶球谐系数。

以CSR RL05球谐系数为例，生成无约束格网解（grid_CSR）、高斯300km滤波解（grid300_CSR）、去相关滤波解（gridP4M6_CSR）和组合滤波解（grid300P4M6_CSR），并以CSR mascon解（C_mas）作为参考。

从结果可以看出40阶的Swarm球谐系数反演的噪声太大，即使经过强约束的高斯滤波后效果仍然不理想，不大清楚是Swarm数据本身所致还是代码实现有误？还请有懂的友友帮忙解答下，不胜感激！ICGEM官网下载 COST-G发布的40×40的球谐系数。

Kusche 等（2007）从时变重力场模型球谐系数求解出发，对球谐系数的最小二乘解的法方程组应用贝叶斯估计，构造了 DDK 滤波。

需下载计算海底压力、海水密度以及其他有关海水信息的“seawater”代码包，还需下载Yan-Ning Kuo编写的用于计算比容海平面变化的代码包（“steric_height_calculation.m”）按照上述盐容和热容海平面变化的计算公式，计算热容时，保持盐度为平均盐度不变；3、读取IPRC温盐数据，并绘制比容海平面变化全球趋势图与盐容、热容以及比容在太平洋某处的时间序列图。

JPL：https://podaac.jpl.nasa.gov/dataset/TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.1_V3 (可能需要科学上网)由于下载的JPL mascon数据为“.nc4”文件，于是将其修改为“.nc”文件；并将文件名开头修改为“JPL_”，以与其它机构的mascon文件名对应起来。以2002年4月的全球EWH格网点为例，

先用GLDAS水文模型的格网产品经过球谐分析得到球谐系数，再经过截断和高斯滤波或者去相关滤波处理得到滤波后的格网，逐格网点求解滤波前后的差异因子。发现经过尺度因子改正之后，泄漏至海洋的信号明显减弱，但是以长江流域为例，求EWH时间序列发现，乘以GRACE/GRACE-FO经过截断和滤波后的格网数据得到尺度因子改正后的解。3、 将GLDAS滤波前和滤波后的格网按时间序列排列进行最小二乘求尺度因子。2、格网数据转换为球谐（质量），再经过组合滤波得到GLDAS滤波解。尺度因子法的基本原理，直接上图。

上述公式出自 Chen J L, Wilson C R, Seo K W. Optimized smoothing of Gravity Recovery and Climate Experiment (GRACE) time‐variable gravity observations[J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B6). 中的公式（6）在JGR: Solid Earth本篇文章网站显示所有引用这篇文章的学者论文。

提前准备好：（1）CSR机构发布的Mascon产品，这里以CSR-Mascon RL06数据为例，准备了从2002年4月至2021年7月 共197个月的数据（实际试验的时候可以并不需要这么多数据）（2）全球陆地格网点的纬度加权值，即代码中的grid_weight。另一方面，南极区域外非冰盖质量变化信号泄漏至南极区域内，形成内泄漏，两者统称为泄漏误差。①M_obs为Mascon解经过球谐分析和500km高斯滤波后的结果，第一次循环开始前，假设M_tru=M_obs，每次循环的过程中，更新M_tru即可；

首先以等效水柱高形式绘制测量信号的阶方差图像，再运用上述前两个公式进行最小二乘拟合，求得拟合参数a、b、c、d，最后利用第三公式计算滤波系数h（注意：对于最大阶数60的球谐产品，上述第一个公式适用于0-20阶，第二个公式适用于31-60阶数）示例代码中滤波半斤r0=200km、r1=300km、选定的次m1=15，绘制滤波系数图如下。

①直接从机构数据读取，经过一阶项、二阶项和取平均等预处理操作后的61*121的球谐系数sc矩阵；②例如经过高斯300km之后的等效水柱高，再利用gmt_grid2cs函数转化为球谐系数。2、利用pcolor函数绘制球谐系数格网图。

1、首先下载区域的边界经纬度坐标（如下图展示的为全球主要流域的边界经纬度），解压RBvect文件可以找到index用记事本打开，以编号为0011的Changjiang Yangtze流域边界为例进行说明，2、获取到研究区域边界后，即可参照下面程序进行区域格网图的绘制。出图效果为这样，需要人为手动修改标注条至合理的位置。

1、首先整理好需要计算不确定度的GRACE数据，这里以CSR level-2截断至60阶的球谐系数、DDK3和4滤波解、CSR mascon以及JPL mascon数据为例（以JPL mascon数据为参照）[2]姚朝龙, 李琼, 罗志才,等. 利用广义三角帽方法评估GRACE反演中国大陆地区水储量变化的不确定性[J]. 地球物理学报, 2019, 62(3):15.2、根据TCH中的相关公式，利用迭代初值计算噪声协方差矩阵R(n*n)中的相关元素。函数求解在满足约束条件的情况下使得目标函数最小的参数。