利用SSA方法插值GRACE資料的空缺資料

目錄

 1.SSA方法簡介

2.SSA缺失值填補方法

3.基于累積分布函數(CDF)的模型選擇

4.論文公開的程式包下載下傳

5.程式包調試

國際各大機構釋出GRACE産品存在資料空缺的問題，尤其是在GRACE和GRACE-FO之間的11個月的資料GAP，一般可以使用周期函數拟合進行資料插值,下圖展示了GRACE資料缺失值的時間分布(Yi et al., 2021,JGR)。Yi and Sneeuw(2021)利用奇異譜分析的方法研究了如何對GRACE缺失的資料進行插值處理。

 1.SSA方法簡介

假設有一個均勻采樣的時間序列

，是以其軌迹矩陣為

矩陣的列數為

，Y的每個上升斜對角線都是相同的值。在之前的研究中，基于Broomhead和King(1986)提出的方法，使用矩陣Y形成滞後協方差矩陣

，并采用主成分分析(PCA)程式進行處理。一種等效但更簡潔的替代方法是通過奇異值分解(SVD)直接對Y進行分解:

其中的下表代表矩陣大小，

 是對角矩陣，

是正交的，

也是滞後協方差矩陣的特征向量。主成分分析中經驗正交函數(EOFs)和主成分(PCs)的輸出可以表示如下：

其中

和

的第

列，以及

的第

個對角線元素。請注意，如果延遲協方差矩陣定義為

，則EOF和PC的定義應該交換。

矩陣可以通過模态

的和來重建，其中每個模态Z都等于EOF與PC的乘法:

每個模态

與

具有相同的結構。是以，我們可以用斜對角元素的平均值來表示一個新的時間序列，稱為重構分量(RCs，記為

):

其中

中的所有元素滿足

。注意，這裡的上标是位置索引。所有

的總和等于原始輸入時間序列:

已經按照其奇異值(即信号強度)降序排列，可以表示長期和振蕩分量或噪聲。是以，通常不使用包含信号(和噪聲)全部資訊的原始

值，而是降低

到

的值以保留所需的

:

2.SSA缺失值填補方法

SSA缺失值填補方法是利用可用樣本的時間相關性來處理缺失樣本的時間序列。這裡，我們采用Kondrashov和Ghil(2006)提出的疊代政策。下圖(Yi et al., 2021,JGR)給出了SSA缺失值填補方法的流程圖。這個想法是在原來的SSA方法中添加兩個循環，并逐漸疊代更新缺失值

。當

的假設值趨于穩定時，内層循環将停止，外層循環将逐漸增加重構序列的複雜度，複雜度由

階模數控制(式6)。

初始置零，并保留更新後的值用于下一輪内外層循環。

 SSA方法有兩個關鍵參數:視窗寬度

和重建成分數

，K的含義很容易了解，但參數M的影響就不那麼直接了，它的值往往是由實驗經驗确定的。M的最優值是采用 (Khan & Poskitt, 2011)的計算方法為

,其中的

，是以對于GRACE的15年觀測資料，其為月采樣，M的取值為12月到140個月。在這裡，我們從軌迹矩陣Y的角度來讨論這種情況。我們通常認為

，但同樣适用于

，因為矩陣是轉置的。在重建方程(5),每個元素的平均時間相關的滞後時間從1到M(除了端點,最大時間間隔短于M)。增大

具有不同特點的相反的效果:它提高了信噪比水準高度重複的信号,而它往往抑制瞬态信号和變頻振動自相關性主要存在于短延遲。從這個角度來看，較大的

不一定是更好的選擇。另一方面，

的長度為

，長度不能太小，以保證

時頻分析的可靠性。

RCs按方差降序排序，通常保留前幾個

。然而，當序列噪聲很大時，前幾個

可能不能代表有用的資訊。是以，我們提出了一個累積分布函數(CDF)測試來判斷一個模式是否為噪聲，隻有通過測試的

才會被保留。我們将在下一節中詳細介紹CDF測試。

GRACE資料缺失類型可根據難度和不确定度分為SSA-filling-a和SSA-filling-b兩類，GRACE内的空白标記為SSA-filling-a, GRACE- FO内11個月的空白标記為SSA-filling-b。在填充SSA-filling-b的缺失資料時，需要使用較少的

數，是以我們采用了兩種不同的政策。在SSA-filling-a中，我們設定所有系數M = 24, K = 12。根據我們的經驗，對于一兩個月的間隔來說，長達兩年的相關性就足夠了。通過省略不同M值(設定為整年，因為結果對M的微小變化不敏感)的一部分已知觀測值的交叉驗證測試表明，M = 12或24總是比較大M值更好地恢複被省略的觀測值。選擇K = 12總是足以保持信号，因為超過K = 12的RCs要麼微不足道，要麼屬于高頻振蕩(HFOs)的範疇，因為較高的RCs通常具有較高的噪聲水準。

在完成SSA-filling-a步驟後，我們将處理SSA-filling-b步驟。SSA-filling-b的M和K參數分别通過交叉驗證确定每個SHC系列。在交叉驗證過程中，使用了2003年至2016年的解，但我們先後遺漏了2004年至2015年的一年資料集;也就是說，我們為12個實驗留出了一年的間隔。我們通過基于M和K的各種組合的SSA gap-fill方法恢複了差距，M在12個月的間隔内從12到72個月，K在1到12之間變化。用12個實驗的均方根值對每組參數進行評價。當恢複值與觀測值的內插補點達到最小時，選擇M和K的最優集合。此外，交叉驗證步驟得到的差異歸因于SSA-filling-b結果的誤差。

3.基于累積分布函數(CDF)的模型選擇

SSA方法對時間序列進行分解，得到各模态的EOF和PC。對某一模态的PC進行傅裡葉變換，得到其功率譜密度，再将功率譜密度按頻率累加，得到頻譜CDF。CDF圖的形狀顯示了其組成頻率的階躍變化，揭示了時間序列的周期性特征。例如，在任何頻率上都沒有偏好的白噪聲，其CDF曲線類似于對角線，這是白噪聲Kolmogorov-Smirnov檢驗的基礎(Massey Jr, 1951;Wouters & Schrama, 2007)。在這裡，我們不隻是拒絕白噪音:我們也拒絕HFOs(稍後定義)。這是必要的一步，因為高頻信号容易受到噪聲的影響。此外，目前的大多數研究很少對HFOs進行調查，主要集中在年度和半年的變化。為了減少資料噪聲，犧牲高頻頻率是值得的。除了這兩個原因外，由于缺乏周期限制，hfo将在較大的資料間隙中劇烈振蕩。是以，我們提出了一個特别的準則，隻保留至少90%的累積能量在頻率低于每年3個周期(即周期≥4個月)的模态，這要求CDF曲線在每年3個周期的頻率下超過0.9。

4.論文公開的程式包下載下傳

在論文的網頁端的附件處，我們可以找到公開分享的程式包。

跳轉至相應的界面，我們可以找到對應的程式包，下載下傳調試。

5.程式包調試

下載下傳程式包進行讀取，我發現代碼缺失一個函數，

% [tt1,X1] = uniform_time(ser(:,1),ser(:,2), [2002,4,2020,8]);
           

這個函數我沒有實作，但是其基本的功能是：把原來有缺失并且連續的數值給對應的缺失資料處以NaN代替。

運作結果

參考資料：

Yi, S., & Sneeuw, N. (2021). Filling the data gaps within GRACE missions using Singular Spectrum Analysis. Journal of Geophysical Research: Solid Earth, 126, e2020JB021227. https://doi. org/10.1029/2020JB021227