分享一个MSSA插值的GRACE level数据集

1. 背景介绍

我们通常使用的GRACE数据包含球谐数据和mascon数据。而不管是球谐产品还是mascon产品，都存在月份数据的缺失，如下图所示（Yi and Sneeuw, 2021）。本专栏分享了一个利用多通道奇异谱分析（MSSA）方法插值好的GRACE level03数据集。以下是文章的摘要（ChatGPT提供翻译支持！）

从2002年3月到2017年10月，重力恢复和气候实验（GRACE）测量了地球重力场的变化（Tapley等，2004）。GRACE任务包括两颗卫星，它们在低轨道、近圆形、近极地轨道上相互跟随，相距约220公里。当领先的卫星飞越大质量物体时，它会比后继卫星稍微受到更多的引力影响，从而使轨道受到不同的扰动。通过精确测量卫星间距的变化，可以通过对两颗卫星的差分重力吸引来量化地球的质量变化。GRACE证明了其在监测地球大气层、海洋、陆地和冰盖内外质量运动方面的重要性，并迅速成为一个必不可少的工具。事实上，在过去的几十年里，GRACE为各个领域提供了深刻的见解，从地球物理学到水文学。例如，从GRACE派生的质量变化观测已经被用来监测全球和区域的陆地水储量（J. Chen等，2016；Long等，2015；Longuevergne等，2013；Syed等，2008），全球海洋质量变化（Gardner等，2013；Morison等，2007；Wouters等，2011），海底压力（Johnson＆Chambers，2013），或最近的冰融化（Luthcke等，2013；Wouters等，2019）。此外，GRACE还揭示了有关地球固体内部发生的过程的宝贵信息，包括地震周期（Bouih等，2022；J. L. Chen等，2007；Panet等，2007）或冰川同衡调整（GIA; Steffen等，2008；Velicogna＆Wahr，2013）。GRACE任务的成功总体上激发了后续任务，即GRACE-Follow On（GRACE-FO; Flechtner等，2016；Landerer等，2020），于2018年5月发射。然而，两个任务之间存在重要的时间差，除了GRACE任务末期缺失观测的增加。然而，拥有足够长、一致和连续的测量时间序列对于研究固体地球过程中发生的长期重力变化，尤其是监测与气候相关的质量变化，如冰盖和冰川的持续演化以及陆地水储存，是至关重要的。

然而，由于两个任务的轨道几何特性，观测在东西方向上的敏感度较低。因此，仪器误差、不足和大气重力场去偏模型（Seo等，2006，2007）中的不准确性等都会导致具有明显南北条纹图案的特殊噪音，限制了GRACE测量的质量和潜在用途，包括更多的地球物理应用（Han等，2004；Swenson＆Wahr，2006；Thompson等，2004）。为了减少这种特征噪音，已经开发了几种信号处理方法，使用各种数学工具。首先，可以使用不同的滤波方法去除原始GRACE观测得到的重力场中的南北条纹。后处理方法的示例包括：高斯滤波器（Seo等，2007；Wahr等，2004），它们的组合（Guo等，2010），或广泛使用的DDK去相关滤波器（Kusche，2007；Kusche等，2009）。DDK滤波器旨在减小重力场球谐分解的Stokes系数之间的相关性，并作为卷积矩阵提供。由于所有滤波方法都需要在平滑（因此空间分辨率和信号衰减）和减小噪音之间寻求折衷，DDK滤波器提供了一组滤波器（DDK1到DDK8），对应于不同的滤波水平。为了进一步减少GRACE和GRACE-FO派生的重力场中的噪音，部分原因在于处理策略的局限，各处理中心提供的解可以在观测水平上组合（COST-G; Jäggi等，2020），或在后处理期间进行平均（Sakumura等，2014）。与解决方案处理和后处理策略相关的不确定性也可以使用合奏方法来考虑（Blazquez等，2018）。另外，GRACE质量集中（mascons）解决方案已经开发出来，以提供抑制泄漏并且可以立即使用的解决方案（Luthcke等，2013；Save等，2016；Watkins等，2015）。然而，要实现这些解决方案需要引入关于信号或噪声结构的时空分布的潜在偏见先验信息，或在最小二乘重力反演中进行正则化（Loomis等，2019）。

同时，统计信号处理技术，特别是统计分解方法，已经用于识别GRACE时间序列中的变化模式。大多数这些方法旨在保留表示大部分地球物理信号变化的一组模式，以过滤掉由南北条纹主导的不太相关的信号部分。特别是，常被应用于隔离GRACE派生重力场时间序列中地球物理信号的是特征空间技术。首先，主成分分析（PCA; Lorenz，1956），也称为经验正交函数分析，已经被用于从GRACE数据中提取主要的正交模式，无论是用于过滤噪音（Chambers，2006；Chambers＆Willis，2008；Schrama等，2007；Wouters＆Schrama，2007），还是提取感兴趣的信号（De Viron等，2006；Rangelova等，2007；Rangelova＆Sideris，2008；Rieser等，2010）。然而，使用PCA提取的模式的物理解释可能受到时间序列中独立信号源的叠加的影响。因此，独立成分分析（ICA）旨在基于信号源的统计独立性来分离主要模式，已经优于PCA（Forootan＆Kusche，2012；Frappart等，2010）。然而，无论是PCA还是ICA，都仅使用现有时间序列之间的信息，忽略了时间序列之间的滞后相关性，因此受限于平稳过程。虽然它们在分离具有不同时间行为的信号方面效率很高，但捕获包括在GRACE数据中的地球物理信号的时空演化性质仍然具有挑战性（Forootan等，2014）。幸运的是，可以使用奇异谱分析（SSA; Vianna等，2007；X. Wang等，2011）在单一时间序列上包含任何滞后信息。此外，多通道（或多变量）-SSA（例如，Ghil等，2002）是PCA和SSA的综合，它使用具有时间滞后观测和多个时间序列，特别适用于捕获GRACE数据的复杂时空变化模式，并过滤处理特定错误和噪音（Prevost等，2019；Rangelova等，2012；F. Wang等，2020；Zotov＆Shum，2010）。

GRACE任务末期缺失观测数量众多，以及任务之间的11个月的观测间隙限制了GRACE和GRACE-FO数据的充分利用。因此，已经进行了一些工作来填补GRACE重力场的时间观测间隙。首先，独立的观测数据已被用来填补GRACE的数据缺失。特别是，可以利用来自卫星激光测距（SLR）或Swarm卫星上的全球定位系统接收机的直接观测来重建地球低次重力场（Jäggi等，2016；Lück等，2018；Richter等，2021）。通过全球导航卫星系统全球网络测量的变形场反演，也可以通过加载理论获得地球低次重力场估计（Chanard等，2018；Rietbroek等，2014；Wu等，2020）。然而，独立数据可能包含特定的技术相关错误或其他可能影响GRACE重力场填充的物理过程（Dong等，2002；Mémin等，2020）。GRACE的时间间隙可以使用数据自适应的统计技术来重建，如SSA和M-SSA，以将时间序列分解为一组时间或时空组分的子集，然后用于重建缺失的观测（Kondrashov＆Ghil，2006）。SSA已经被用于通过GRACE重力场SH分解系数的时间序列执行迭代式间隙填充（Li等，2019；Prevost等，2019；Yi＆Sneeuw，2021）。M-SSA也已经证明了它能够重建缺失的观测，至少对于地球重力场的低次SH系数，使用Swarm观测（F. Wang等，2021），或一部分的重力变化，即气候驱动的水储量变化，使用降水和温度模型（Humphrey＆Gudmundsson，2019；Yang等，2021）。最近，更先进的机器学习技术也已经用于在GRACE和GRACE-FO观测时段内进行间隙填充。示例包括使用基于水文气候数据驱动的贝叶斯卷积神经网络来重建地球重力场的陆地水分量（Mo等，2022），或结合M-SSA和人工神经网络的算法（Lai等，2022）。然而，这些方法比经典的统计方法更复杂，计算上更具挑战性，通常仅适用于地球水储存应用，丢弃了与固体地球过程相关的质量变化。

2. 论文内容

数据集涉及的文章为Gauer et al.(2023)年发布的。文章的摘要的主要内容为：

重力恢复和气候实验及其继续（GRACE/-FO）全球月度地球重力场测量已经在量化质量转移方面取得了显著进展。然而，任务之间存在重要的时间差，阻碍了对长期质量变化的评估。此外，仪器和处理误差会导致大规模的非物理性南北条纹，污染了地球物理信号。我们使用多通道奇异谱分析（M-SSA）来克服这两个问题，利用GRACE/-FO Level-2解的时空信息，经过DDK7去相关和一种新的互补滤波器的滤波，该滤波器是根据完全处理的数据和观测参数化拟合之间的残余噪声构建的。通过对由格拉茨理工大学空间研究中心、GeoForschungsZentrum、测量学研究所和喷气推进实验室处理的等效水高（EWH）时间序列进行迭代M-SSA，我们替代了缺失的数据和异常值，以获得一个组合的均匀采样解。然后，我们应用M-SSA来提取每个EWH时间序列与其纬度相邻的时间序列之间的共同信号，以进一步减少残余的空间不相关噪声。将GRACE/-FO M-SSA解与卫星激光测距和Swarm低次地球重力场以及水文模型进行比较，证明了它能够令人满意地填补缺失的观测。我们的解决方案在海洋上实现了与质量集中（mascon）解决方案相当的噪声水平（3.0毫米EWH），无需先验信息或正则化。虽然使用高度滤波的球谐函数或mascon解决方案难以捕捉短波信号，但我们展示了我们的技术能够有效地恢复与水库蓄水相关的局部质量变化，这些质量变化是有充分记录的。

关键词：

1.利用多通道奇异谱分析对Gravity Recovery And Climate Experiment/Gravity Recovery And Climate Experiment Follow On的重力场进行了缺口填补和时空滤波。

2.Lobe-Edge谱滤波器，作为广泛使用的DDK去相关滤波器的补充，有助于减少条纹噪音。

3.最终解决方案显示出最小的噪音内容，相对于其他解决方案，具有检测更小尺度信号的潜力。

MSSA数据处理前的数据预处理：在使用GRACE和GRACE-FO数据集进行地球物理应用之前，需要进行几个预处理步骤。首先，使用GeoForschungsZentrum（GFZ）、Jet Propulsion Laboratory（JPL）和Center of Space Research（CSR）解决方案的各自提供的对应系数，根据Technical Note 13（TN-13; Swenson et al., 2008; Sun et al., 2016）以及Institute of Geodezy at Graz University of Technology（GRAZ）解决方案的系数平均值，来考虑不可观测的一阶SH地球重心引力系数。然后，由于接近极地轨道和加速度计误差（Klinger＆Mayer-Gürr，2016）难以观测，将替换C2,0地球扁率和C3,0重力系数，根据Technical Note 14（TN-14; J. Chen et al., 2005; Loomis et al., 2020）与SLR观测数据。

此外，为了研究地球重力场的变化，我们分别去除了每个Level-2解决方案在2003年至2022年期间估计的均值。因此，由于已知现象的重力场修正模型中的仪器误差或缺陷引起的特征非物理性的南北拉长的条纹模式占据了GRACE和GRACE-FO解决方案的主导地位。图1a和1b显示了分别表示为等效水高度（EWH）的2008年7月和2019年7月的GRACE和GRACE-FO重力场的示例。南北条纹伪影的大振幅强调了在进行任何地球物理应用之前需要过滤GRACE和GRACE-FO重力场的必要性（Sakumura等，2014）。在这里，我们首先使用非各向同性的去相关滤波器，称为DDK（Kusche，2007; Kusche等，2009）。DDK基于使用误差和信号协方差信息的正则化。该滤波器产生一个滤波矩阵，该矩阵源自于2003年8月GRACE解决方案的先验误差协方差，我们将其应用于所有GRACE和GRACE-FO月度重力场。该滤波器提供了8个级别，从最强的DDK1到最弱的DDK8级别，主要影响SH分解的高阶系数，其中包含大部分的条纹噪音。增加DDK滤波的级别会导致较大的信号衰减和泄漏，导致地球物理信号扩散到更大的区域。因此，必须在解决方案过滤和降噪之间做出权衡。虽然DDK滤波器的选择取决于应用，但DDK5似乎是最常用于地球物理应用的，而在较早的Level-2 SH解决方案中，平均解决方案已被证明可以减少处理工件（Sakumura等，2014）。因此，我们以四个解决方案（CSR、JPL、GFZ、GRAZ）的平均值作为参考解决方案，然后将它们组合以创建经DDK5滤波的GRACE-MSSA解决方案。DDK5有效地去除了大部分的南北条纹，同时保留了地球物理信号，其空间分辨率可以与等效的180千米高斯滤波半径等效进行比较，通过比较非各向同性DDK滤波器的各向同性部分与高斯滤波器获得（Cambiotti，2020；Kusche等，2009；图1c和1d）。在这里，我们选择应用DDK7滤波器，其空间分辨率相当于145公里各向同性高斯滤波半径等效。下图展示了数据处理的流程（Gauer et al., 2023）：

3.数据读取与存储

从以下的路径可以下载得到数据集，并采用matlab进行数据的读取与存储，并成图展示：

https://dataverse.ipgp.fr/dataset.xhtml?persistentId=doi:10.18715/IPGP.2023.lgquie56

% read raw data
% DATA DOWNLOAD FROM:https://dataverse.ipgp.fr/dataset.xhtml?persistentId=doi:10.18715/IPGP.2023.lgquie56
% 我是水怪的哥-Monster Bros
% 2023-10-22
address = 'C:\Users\UCASer\Desktop\IPGP_GRACE_LEVEL\';
GFA = dir(fullfile(address,'*.xyz'));
k=length(GFA);
for i = 1:237
    year = GFA(i).name;
    yy = str2num(year(12:15));
    mm = str2num(year(17:18));
    GRACE_MSSA.tt(i,1) = time_transfer([yy,mm,15],1);
    Data = importdata([address,GFA(i).name]);
    if(i==1)
        % run only once!
        GRACE_MSSA.lon = reshape(Data.data(:,1),181,361);
        GRACE_MSSA.lat = reshape(Data.data(:,2),181,361);
    end
    GRACE_MSSA.rg(:,:,i)  = reshape(Data.data(:,3),181,361);
    disp(i)
end

并进行成图显示，下图分别展示了2010年一年的水储量变化以及长期的变化趋势。（单位：cm|cm/yr）

同时我们可以将csr mascon产品的数据得到的长期变化趋势成图显示

对比长期趋势，我们知道Gauer et al.(2023)年发布的包含GIA的信号，在数据应用时需要考虑。

参考资料

Gauer, L. M., Chanard, K., & Fleitout, L. (2023). Data‐driven gap filling and spatio‐temporal filtering of the GRACE and GRACE‐FO records. Journal of Geophysical Research: Solid Earth, e2022JB025561.

感谢chatGPT的翻译支持!

欢迎点赞收藏!