GRACE Level-2數據的RTS狀態平滑

鳳 勇,常國賓,錢妮佳,魏征強,宦越洋,楊一帆

中國礦業大學環境與測繪學院,江蘇 徐州 221116

全球氣候變暖導致的冰蓋融化、海平面上升是制約人類社會發展的重大問題。GRACE衛星自2002年3月發射到2017年7月結束任務,對全球重力場變化進行了十幾年的觀測。目前GRACE衛星重力技術廣泛應用于陸地水存儲量變化、南極和格陵蘭島冰蓋質量變化及全球海平面變化等研究[1-2]。后續GRACE-FO(GRACE follow-on)衛星的發射為繼續研究地表物質遷移提供基礎數據。

球諧系數反演法計算簡單,應用較為廣泛,但受到衛星載荷的儀器測量誤差、混頻誤差、衛星軌道、先驗模型誤差和時變重力場截斷誤差的影響[3-5],不同研究中心發布的GRACE level-2產品解算出來的地表物質遷移反演結果存在顯著的南北方向條帶。除此之外,GRACE時變重力場模型球諧系數存在明顯的高頻噪聲,尤其是40階后的重力場位系數,其噪聲大于信號,主要原因是GRACE衛星軌道設計缺陷、背景模型誤差及星上載荷測量精度導致的觀測誤差等的共同作用[6-7]。本文利用狀態空間模型對GRACE level-2產品處理以削弱條帶誤差。用到GRACE level-2產品包括截斷至60階的球諧系數以及球諧系數對應的形式誤差協方差矩陣。

為削弱誤差對真實信號的影響,往往會對數據施加某種約束,這種約束大體上可劃分為兩類:一是空間約束,利用數據在空間域上存在的聯系,對數據加以約束,例如各向同性或各向異性平滑[8]、經驗去相關[9]、基于反演或正則化的去相關去噪核(decorrelation and denoising kernel,DDK)濾波[10-12]、獨立分量分析(independent component analysis,ICA)[13]等;二是時間約束,對同一組觀測量進行連續觀測構成時間序列,運用合適的方法表征觀測量在時間域上的特性,比如小波分析[14]、時差正則化[15]等。有些方法不僅存在空間約束,還同時兼顧時間約束,比如,主成分分析(principal component analysis,PCA)[16-17]、奇異譜分析(singular spectrum analysis,SSA)[18]等。容易發現,上述時間域約束是一類定性層面的非參數化方法,沒有將形式誤差協方差矩陣考慮進來,因而損失大量統計信息。參數化的方法也不乏學者研究,比如用關于時間的多項式來參數化重力信號(球諧系數),用關于球諧系數次數的多項式函數來參數化條帶誤差,所有參數都用加權最小二乘法同時進行估計。本文所研究的狀態空間模型法也是一種參數化方法,采用狀態空間模型對信號進行參數化表示,進而采用卡爾曼濾波對GRACE Level-2球諧系數進行處理,以達到降噪和去條帶的效果。本文方法可以綜合利用真實地球物理信號的協方差矩陣及條帶誤差的協方差矩陣,不僅反映條帶誤差的空間域特征,而且能反映地球物理信號在時間域上的相關性。

首先,狀態向量中僅包含表示真實地球物理信號的球諧系數,GRACE Level-2產品中的測量誤差協方差矩陣被用作狀態空間模型中觀測誤差協方差矩陣。然后,狀態向量的協方差矩陣用一個僅與階數有關的冪律模型來表示,以便提供月與月之間球諧系數的不確定性度量;引入一個尺度因子(方差分量)表示該不確定性的總體幅度,并根據數據對該因子進行自適應調整。最后,將平滑解而非濾波解作為最終的結果,進一步改善數據處理效果[19-20]。

DDK濾波代價函數由正則化項和殘差加權平方和項兩部分組成,分別涉及信號協方差矩陣和誤差協方差矩陣,恰與狀態空間DDK濾波方法(state space decorrelation and denoising kernel,SS-DDK)中的過程噪聲協方差矩陣和觀測噪聲協方差矩陣相對應,只不過SS-DDK中的過程噪聲協方差矩陣表示球諧系數在相鄰月份變化的統計信息,而DDK中的信號協方差矩陣表示的是信號在單個月份的統計信息,這一相似性是本文方法名稱中“DDK”的來源。

1 數據和方法

1.1 數據

采用美國得克薩斯大學空間研究中心(Center for Space Research,University of Texas at Austin)提供的2002年4月至2014年6月共135個月GRACE Level-2(Release-05)數據。該數據為60階無約束球諧系數,并扣除了非潮汐大氣、高頻海洋信號、各種潮汐、固體潮和固體極潮等影響。數據包括球諧系數yk,球諧系數對應的完整形式誤差協方差矩陣Rk,下標“k”表示月份。丟失數據的月份直接舍棄,若第二個月數據丟失,可以簡單地將下標“2”用來表示第三個月的數據。

1.2 狀態空間模型

表示真實地球物理信號的球諧系數xk是待估計的變量,與yk有相同的維數。狀態空間模型的觀測模型表示為

(1)

式中,觀測矩陣Hk是單位陣。觀測誤差vk一般包含高頻誤差和條帶噪聲等。條帶噪聲的統計規律與時間無關,是一種時間域純隨機誤差,即白噪聲,將條帶噪聲納入vk中加以考慮,而不是將其作為狀態向量的一部分,這更符合實際,也可降低狀態向量的維數,認為條帶噪聲的統計信息全部包括在觀測噪聲協方差矩陣中(事實也是如此),則無須對條帶噪聲統計信息進行其他人為設置,簡化了建模工作量,提高模型可靠性,這是本文狀態空間法與前人最大的不同之處。

狀態轉移方程(過程模型)表示為

(2)

式中,狀態轉移矩陣Fk是單位陣;Q是過程噪聲協方差矩陣,一般為對角陣,考慮到GRACE存在月份缺失問題,Qk=mkQ,mk表示yk與yk-1之間的月數差,Q的對角線元素根據冪律模型得到,冪律模型的冪指數表示為-μ,底數為階數l;α表示方差分量。式(2)表示對任意階次的球諧系數采用的均為隨機游走過程模型。DDK和SS-DDK濾波解分別為

(3)

(4)

不難發現,式(3)和式(4)存在相似之處,均包含觀測誤差所構成的代價函數這一項。式(3)中引入正則化因子λ和信號協方差矩陣S[12],并根據Roelof Rietbroek發布在Github上關于DDK濾波器的軟件(https:∥github.com/strawpants/GRACE-filter),進一步得知從DDK1到DDK8,λ取值介于5×109～1×1014之間,矩陣S的對角線元素為l-v,l表示階數,指數v一般取4。式(4)引入方差分量因子α和過程噪聲協方差矩陣Q,其中α根據公式迭代計算得到,矩陣Q則和矩陣S設計相同,均由冪律模型得到,冪指數μ取4。式(3)和式(4)所體現出的相似性正是將本文方法命名為“SS-DDK”的原因。

1.3 狀態估計

給定建模參數,最優狀態向量估計一般選用卡爾曼濾波或平滑[21]。這里給出卡爾曼濾波的計算公式

(5)

(6)

(7)

(8)

(9)

RTS(Rauch-Tung-Striebel)平滑是一種最為常用的固定間隔平滑方法,公式為[22]

(10)

(11)

(12)

(13)

PN,N-1|N=(E-KkHk)Fk-1Pk-1|k-1

(14)

1.4 參數估計

卡爾曼濾波解和RTS平滑解容易得到,然而建模參數α是未知的。本文試圖同時找到狀態向量和建模參數的最大似然估計。對應的代價函數,即負對數似然函數為

(15)

注意Qk中包含冪指數μ。當使用EM算法進行參數估計時,E-step輸出只含建模參數為未知數的近似代價函數[23]

(16)

M-step中找到使得近似代價函數式(16)取得最小值的參數α。經證明,存在如下封閉形式的解

(17)

其中Dk的計算公式如下

(18)

如圖1所示,詳細闡述SS-DDK實施步驟。步驟1對GRACE Level-2數據進行預處理,根據給定的冪指數和階數構建過程噪聲協方差矩陣Q;步驟2是人為給定初值和初始協方差矩陣,并假定尺度因子(方差分量)α的初值為1;步驟3按式(5)—式(9)逐月計算狀態向量的卡爾曼濾波估計值和對應協方差矩陣;步驟4按式(10)—式(14)計算RTS平滑解以及相應協方差矩陣;步驟5運用卡爾曼濾波解和RTS平滑解求得尺度因子α的更新值;步驟6判斷參數α是否收斂,若收斂,輸出RTS平滑解作為最終處理結果,否則,更新尺度因子α繼續執行步驟3—步驟6,直至迭代停止。

圖1 狀態空間建模流程

按狀態空間流程圖(圖1)進行迭代處理,得到迭代50次的結果,最后參數α大約收斂至1×10-19(圖2)。

圖2 狀態空間模型中方差分量因子α隨迭代次數的變化

2 試驗和結果討論

2.1 濾波質量異常分析

美國得克薩斯大學空間研究中心目前只發布RL05球諧系數解的協方差矩陣,故沒有選用RL06數據進行處理。RL05數據可以從網站(http:∥download.csr.utexas.edu/outgoing/grace/)下載。該網站僅提供2002年4月至2014年6月的數據,后續GRACE在軌運行期間的數據暫未提供。這135個月數據足以驗證SS-DDK方法的性能。試驗前先做如下預處理:從球諧系數中扣除該時段的平均值;用衛星激光測距(satellite laser ranging,SLR)[24]結果替換球諧系數C20項;地心1階項由GRACE衛星和海洋模型聯合估計得到;利用ICE-6G_D模型進行冰川均衡化調整(glacial isostatic adjustment,GIA)。選用DDK1—8濾波與本文SS-DDK對比。CSR RL06 v2.1 mascon解(https:∥www2.csr.utexas.edu/grace/RL06_mascons.html)可以得到局部較好的質量變化估計[25-26],故將其作為參照數據。

大地水準面是一個重力等位面,大地水準面高可由時變重力場球諧系數計算得到,地球表面及其內部的物質運動引起大地水準面高度的變化(ΔN),公式為

(19)

ΔSlmsin(mλ))

(20)

式(20)得到全球范圍內等效水柱高變化[27],其中ΔH表示等效水柱高變化;ρave為地球平均密度,取5517 kg/m3;ρω為水的密度,取1000 kg/m3;kl為l階的負荷勒夫數。在扣除冰川均衡化調整并且沒有地震影響的情況下,EWH變化即表征地表質量變化。

本文以2006年1月、3月、6月、12月為例,圖3為SS-DDK、DDK3和DDK5濾波解及0.25°×0.25° mascon解的EWH變化,圖4為球諧系數無約束解和其余DDK濾波解的EWH變化。從圖3和圖4中不難發現,DDK1—8濾波的約束能力依次減弱,DDK1的約束能力最強,全球范圍內南北條帶去除得最干凈,但局部區域出現失真現象,尤其是高緯度地區,如格陵蘭島;DDK8的約束能力最弱,局部區域信號保留較好,但是存在明顯的條帶噪聲。一般地,DDK3—5濾波的效果最佳,可以在去除條帶噪聲的同時保留更多有用的信號。重點關注南美洲亞馬孫河流域、非洲剛果河流域、我國長江流域、印度恒河流域和格陵蘭島5大區域的信號保留,以及去條帶情況。經SS-DDK濾波之后,能夠很好地削弱條帶誤差的影響,幾乎看不到明顯的條帶存在,并且在上述5大區域,可以較為清晰地看到EWH的變化,例如在亞馬孫河流域,2006年6月前后等效水柱高EWH為正值,并且變化異常明顯,而在2006年12月前后EWH為負值。再例如格陵蘭島區域,SS-DDK與DDK3—8一樣保存了絕大部分的地球物理信號,不像DDK1—2那樣具有明顯的信號失真。比較SS-DDK和DDK1—8與mascon在格陵蘭島區域的EWH變化發現,任何濾波解都存在不同程度的信號泄漏現象。泄漏誤差不是本文的研究重點,故這里不做過多展開。以mascon解為參考,計算2002年4月至2014年6月整個時間段全球720×1440個格網點EWH的均方根差異(root mean square differences,RMSD)見表1,SS-DDK的均方根最小,DDK3的均方根與SS-DDK最為接近,從全球來看,SS-DDK方法要優于DDK系列濾波。

表1 全球格網點EWH的均方根

圖3 2006年1、3、6和12月分別經DDK3、DDK5和SS-DDK處理后的結果及對應月份的mascon解

為更進一步地論證SS-DDK去條帶和保留信號的能力,選取上述5個區域分別繪制區域整體的EWH變化的折線圖(圖5左側一欄)以及與mascon解差異的折線圖(圖5右側一欄)(繪圖時忽略存在數據缺失的月份)[28]。在圖5左側一欄前4個區域,SS-DDK和DDK系列濾波結果與mascon解在整個時間段內符合地較好;在亞馬孫河流域、剛果河流域和恒河流域部分月份,mascon解的振幅大于其余濾波解。從整個時間序列來看,無論哪一種解前4個區域的EWH都存在周期性變化。在格陵蘭島,SS-DDK和DDK1—8都與mascon解存在不同程度的偏離,這可能是由于信號泄漏所造成的,但SS-DDK與DDK系列濾波之間符合地較好。圖5右側一欄展示SS-DDK和DDK1—8與mascon解的差異,SS-DDK濾波解的質量異常差異(以mascon為參考)與DDK系列濾波解的質量異常差異相近,在部分區域的部分時間月份SS-DDK質量異常差異要小于任何一種DDK,但在亞馬孫河流域和剛果河流域存在初始月份偏差較大的問題,可能由于滾動求解時第一年第一個月的約束解僅求解一次,導致SS-DDK濾波解在初始月份與mascon解的質量異常差異相較于DDK1—8濾波更大一些。

2.2 方差和協方差分析與不確定性估計

測量誤差協方差矩陣R通常是不準確的,往往高估了數據的準確性。為解決這一問題,可以簡單引入方差分量因子δ2來對協方差矩陣R進行縮放。因此,本文通常用δ2R來代替R來表征球諧數據誤差的大小[29]。方差分量因子的估計公式如下

(21)

式中,p表示約束解中非零元素的個數;n表示方程數,球諧系數的最大展開階數為60,由式(1)和式(2)可知,n=3717×2;p取3717。式(20)表征球諧系數變化與等效水柱高變化之間的轉換關系,可簡單縮寫為

ΔH=gTx

(22)

(23)

表2 SS-DDK和DDK2—5濾波在全球和區域的等效水柱高平均不確定度大小

圖6 6大區域SS-DDK和DDK2—5濾波解估計的等效水柱高不確定性隨時間的變化

由圖5—圖6分析可知,SS-DDK估計的等效水柱高不確定性要小于DDK3—5濾波,但大于DDK2。綜合考慮,SS-DDK解的不確定性大小介于DDK2和DDK3之間,與DDK3濾波最為接近。通過對比各約束解的不確定性大小,充分說明SS-DDK的去條帶效果是可靠的。

圖7 經SS-DDK處理后2006—2011年等效水柱高全球格網點的平均不確定性

2.3 時間序列擬合分析

本文選用的GRACE觀測時段是2002年4月至2014年6月(含缺失數據月份),反演135個月GRACE球諧系數的數據,并利用不同濾波方法(SS-DDK、DDK1—8)進行約束求解得到EWH。對全球范圍內所有格網點的EWH采用線性擬合模型進行最小二乘擬合,得到線性年變化率、周年振幅及半周年振幅。線性擬合模型具體為[30]

ΔH=a0+a1(t-t0)+b1cos(ωt)+b2sin(ωt)+b3cos(2ωt)+b4sin(2ωt)+ε

(24)

式中,ΔH為等效水柱高變化;t為時變重力場模型時間(以年為單位);t0為參考時間;a0為常數項;a1為年變化率,表征EWH年變化快慢;ε為擬合殘差;b1和b2的平方和開根號表示周年振幅(annual amplitude,AAMP);b3和b4的平方和開根號表示是半周年振幅;ω為頻率。區域EWH變化的RMSD計算公式為

(25)

圖8展示從2002年4月至2014年6月經濾波之后的周年振幅,考慮到無約束解中已經扣除非潮汐大氣、高頻海洋信號、各種潮汐的影響,所以海洋上信號的周年振幅恰好反映濾波方法的去條帶水平。在海洋區域,SS-DDK方法的周年振幅低于DDK5—8濾波,與DDK3、DDK4最為相似;在信號較強的區域,如亞馬孫河流域、剛果河流域、長江流域及恒河流域,SS-DDK的周年振幅與DDK系列濾波相比沒有發現明顯的條帶誤差。圖9表示全球格網點的RMSD,SS-DDK與DDK3在誤差水平上最為接近。從全球格網點的周年振幅和均方根來看,SS-DDK的去條帶能力強于DDK1—2和DDK5—8。

圖8 2002年4月至2014年6月不同濾波方法EWH的周年振幅

圖9 全球格網點2002年4月至2014年6月不同濾波方法EWH的RMSD

2.4 信號和噪聲水平分析

重力場模型的精度可以用各階大地水準面誤差的平方和開根號來表示[31-32]

(26)

圖10 各種約束解在2007年2月和2010年3月的大地水準面階誤差變化

2.5 仿真試驗分析

由于實際地表質量變化的真值無從得知,無法更加客觀地評價狀態空間模型SS-DDK的去條帶和保留信號的性能,故進行仿真試驗,進一步說明SS-DDK性能的優劣。

假定CSR機構發布的RL06 v2.1 mascon解為對地表質量實際變化的真實反映,將mascon格網點質量轉為最大階數為60的球諧系數,人為在球諧系數中加入條帶誤差,分別利用狀態空間SS-DDK和DDK2—5濾波進行處理得到1°×1°的等效水柱高全球格網圖和區域等效水柱高時間序列差異圖。具體操作如下:

(1) 利用CSR發布的mascon數據得到2006—2011年全球1°×1°的等效水柱高格網圖作為真值,圖11展示2009年上半年mascon的全球格網。

圖11 2009年1—6月CSR RL06 v2.1 mascon解1°×1°等效水柱高全球格網

(2) 在mascon解轉換得到的球諧系數中人為加入條帶誤差,誤差的協方差矩陣已知,圖12展示2009年上半年加入條帶噪聲之后的未濾波全球格網。

圖12 2009年1—6月mascon解人為加入條帶誤差后1°×1°等效水柱高全球格網

(3) 用狀態空間SS-DDK濾波對未濾波解進行處理,圖13展示2009年上半年經SS-DDK濾波之后的全球格網。

圖13 2009年1—6月未濾波解經SS-DDK后1°×1°等效水柱高全球格網

(4) 分別用DDK2—5濾波處理未濾波解,將濾波結果和SS-DDK結果與mascon解(真值)進行比較分析,圖14展示在亞馬孫河流域、剛果河流域、長江流域及格陵蘭島區域各濾波解與mascon解的差異。

圖14 SS-DDK和DDK2—5濾波解相對于mascon解(真值)在亞馬孫、剛果、長江和格陵蘭島的等效水柱高差異

對比圖11—圖13不難發現,2009年上半年mascon解和SS-DDK解的等效水柱高全球格網圖上,條帶誤差幾乎被去除,SS-DDK解在亞馬孫河、剛果河、長江和格陵蘭島等區域的信號與mascon解基本一致。對比4大區域2006—2011年未濾波解和各種濾波解與mascon解的差異發現,各濾波解與真值的差異都在0附近,且SS-DDK解的差異比DDK2—5解更小,說明經狀態空間SS-DDK后的處理結果更加貼近真值。

3 結 論

本文提出一種狀態空間模型方法,用于處理GRACE Level-2產品。該方法與傳統狀態空間法的區別在于:狀態向量只包含球諧系數,而將GRACE Level-2中的協方差矩陣用作觀測誤差協方差矩陣;基于數據對過程噪聲協方差矩陣的尺度因子(方差分量)進行自適應調整;采用平滑解而非濾波解作為最終的數據處理結果。

選取CSR發布的RL05 2002年4月至2014年6月的GRACE Level-2產品,得到經約束后的球諧系數和全球地表質量變化,并與DDK1—8濾波和CSR RL06 v2.1 mascon進行對比分析。分析結果表明:

(1) 分別繪制SS-DDK、DDK1—8及CSR RL06 v2.1 mascon全球范圍的EWH圖。對比發現在高緯度的格陵蘭島,SS-DDK與DDK3—8濾波都能保留更多的信號。在其他區域,例如亞馬孫河流域,在2006年6月和12月前后,SS-DDK方法與DDK3、5濾波一樣具有明顯的等效水柱高(EWH)變化。對比不同方法全球所有網格點,SS-DDK的均方根最小為10.36 cm,即SS-DDK方法的去條帶和保留信號的能力與mascon解最為接近。

(2) 從5大區域入手,繪制約束解和mascon解流域質量異常圖以及各約束解與mascon解的質量異常差異圖。在格陵蘭島,所有約束解均較大地偏離mascon解,且DDK系列濾波解的偏離均大于SS-DDK解,但各約束解之間符合地較好。在其他4個區域,各約束得到的時間序列圖像都存在周期性變化;SS-DDK與DDK3—5有不同程度的相似,而且在某些時段SS-DDK方法的質量異常差異要小于任意一種DDK濾波解。

(3) 分析6大區域的不確定性時間序列圖發現,SS-DDK估計的等效水柱高不確定性要小于DDK3—5濾波,但大于DDK2濾波,綜合考慮全球和各個區域的平均不確定度大小可知,SS-DDK估計結果的不確定度介于DDK2和DDK3濾波之間。需要說明的是,試驗中發現整個時段數據連續處理效果不佳,每5年(60個月)為一個單位進行整體處理最為合適,能夠充分顧及球諧系數的誤差特性。

(4) 對不同約束解和mascon解135個月的EWH進行線性模型擬合,求解流域RMSD和年振幅,并繪制全球格網點圖像。進一步說明SS-DDK解在格陵蘭島的均方根最小,存在更少的信號失真。在海洋區域,SS-DDK解的年振幅低于DDK5—8,與DDK3—4較為相近;在強信號區域,SS-DDK解不存在明顯的條帶誤差;SS-DDK與DDK3在誤差水平上最為接近。

(5) 從不同約束解的大地水準面階誤差來看,SS-DDK解中噪聲水平低于DDK4—8,保留信號的能力與DDK2—3解相當。

(6) 增添仿真試驗發現SS-DDK解與mascon解(真值)的差異相較于DDK系列濾波更加接近,更加充分說明SS-DDK方法在去條帶和保留信號方面具有較好的性能。

后續研究中可以嘗試使過程噪聲協方差矩陣的方差分量參數隨月份變化,以反映不同時間處信號變化幅度不一致的現象,這需要對方差分量進行逐月調整。