鳳 勇,常國賓,錢妮佳,魏征強(qiáng),宦越洋,楊一帆

中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院,江蘇 徐州 221116

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

球諧系數(shù)反演法計(jì)算簡單,應(yīng)用較為廣泛,但受到衛(wèi)星載荷的儀器測量誤差、混頻誤差、衛(wèi)星軌道、先驗(yàn)?zāi)Ｐ驼`差和時(shí)變重力場截?cái)嗾`差的影響[3-5],不同研究中心發(fā)布的GRACE level-2產(chǎn)品解算出來的地表物質(zhì)遷移反演結(jié)果存在顯著的南北方向條帶。除此之外,GRACE時(shí)變重力場模型球諧系數(shù)存在明顯的高頻噪聲,尤其是40階后的重力場位系數(shù),其噪聲大于信號(hào),主要原因是GRACE衛(wèi)星軌道設(shè)計(jì)缺陷、背景模型誤差及星上載荷測量精度導(dǎo)致的觀測誤差等的共同作用[6-7]。本文利用狀態(tài)空間模型對(duì)GRACE level-2產(chǎn)品處理以削弱條帶誤差。用到GRACE level-2產(chǎn)品包括截?cái)嘀?0階的球諧系數(shù)以及球諧系數(shù)對(duì)應(yīng)的形式誤差協(xié)方差矩陣。

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

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

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

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)

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

1.2 狀態(tài)空間模型

表示真實(shí)地球物理信號(hào)的球諧系數(shù)xk是待估計(jì)的變量,與yk有相同的維數(shù)。狀態(tài)空間模型的觀測模型表示為

(1)

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

狀態(tài)轉(zhuǎn)移方程(過程模型)表示為

(2)

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

(3)

(4)

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

1.3 狀態(tài)估計(jì)

給定建模參數(shù),最優(yōu)狀態(tài)向量估計(jì)一般選用卡爾曼濾波或平滑[21]。這里給出卡爾曼濾波的計(jì)算公式

(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 參數(shù)估計(jì)

卡爾曼濾波解和RTS平滑解容易得到,然而建模參數(shù)α是未知的。本文試圖同時(shí)找到狀態(tài)向量和建模參數(shù)的最大似然估計(jì)。對(duì)應(yīng)的代價(jià)函數(shù),即負(fù)對(duì)數(shù)似然函數(shù)為

(15)

注意Qk中包含冪指數(shù)μ。當(dāng)使用EM算法進(jìn)行參數(shù)估計(jì)時(shí),E-step輸出只含建模參數(shù)為未知數(shù)的近似代價(jià)函數(shù)[23]

(16)

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

(17)

其中Dk的計(jì)算公式如下

(18)

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

圖1 狀態(tài)空間建模流程

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

圖2 狀態(tài)空間模型中方差分量因子α隨迭代次數(shù)的變化

2 試驗(yàn)和結(jié)果討論

2.1 濾波質(zhì)量異常分析

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

大地水準(zhǔn)面是一個(gè)重力等位面,大地水準(zhǔn)面高可由時(shí)變重力場球諧系數(shù)計(jì)算得到,地球表面及其內(nèi)部的物質(zhì)運(yùn)動(dòng)引起大地水準(zhǔn)面高度的變化(ΔN),公式為

(19)

ΔSlmsin(mλ))

(20)

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

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

表1 全球格網(wǎng)點(diǎn)EWH的均方根

圖3 2006年1、3、6和12月分別經(jīng)DDK3、DDK5和SS-DDK處理后的結(jié)果及對(duì)應(yīng)月份的mascon解

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

2.2 方差和協(xié)方差分析與不確定性估計(jì)

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

(21)

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

ΔH=gTx

(22)

(23)

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

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

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

圖7 經(jīng)SS-DDK處理后2006—2011年等效水柱高全球格網(wǎng)點(diǎn)的平均不確定性

2.3 時(shí)間序列擬合分析

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

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

(24)

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

(25)

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

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

圖9 全球格網(wǎng)點(diǎn)2002年4月至2014年6月不同濾波方法EWH的RMSD

2.4 信號(hào)和噪聲水平分析

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

(26)

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

2.5 仿真試驗(yàn)分析

由于實(shí)際地表質(zhì)量變化的真值無從得知,無法更加客觀地評(píng)價(jià)狀態(tài)空間模型SS-DDK的去條帶和保留信號(hào)的性能,故進(jìn)行仿真試驗(yàn),進(jìn)一步說明SS-DDK性能的優(yōu)劣。

假定CSR機(jī)構(gòu)發(fā)布的RL06 v2.1 mascon解為對(duì)地表質(zhì)量實(shí)際變化的真實(shí)反映,將mascon格網(wǎng)點(diǎn)質(zhì)量轉(zhuǎn)為最大階數(shù)為60的球諧系數(shù),人為在球諧系數(shù)中加入條帶誤差,分別利用狀態(tài)空間SS-DDK和DDK2—5濾波進(jìn)行處理得到1°×1°的等效水柱高全球格網(wǎng)圖和區(qū)域等效水柱高時(shí)間序列差異圖。具體操作如下:

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

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

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

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

(3) 用狀態(tài)空間SS-DDK濾波對(duì)未濾波解進(jìn)行處理,圖13展示2009年上半年經(jīng)SS-DDK濾波之后的全球格網(wǎng)。

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

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

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

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

3 結(jié) 論

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

選取CSR發(fā)布的RL05 2002年4月至2014年6月的GRACE Level-2產(chǎn)品,得到經(jīng)約束后的球諧系數(shù)和全球地表質(zhì)量變化,并與DDK1—8濾波和CSR RL06 v2.1 mascon進(jìn)行對(duì)比分析。分析結(jié)果表明:

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

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

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

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

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

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

后續(xù)研究中可以嘗試使過程噪聲協(xié)方差矩陣的方差分量參數(shù)隨月份變化,以反映不同時(shí)間處信號(hào)變化幅度不一致的現(xiàn)象,這需要對(duì)方差分量進(jìn)行逐月調(diào)整。