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

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

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

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

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

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

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

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

1.1 數(shù)據(jù)

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

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

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

(1)

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

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

(2)

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

(3)

(4)

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

1.3 狀態(tài)估計

給定建模參數(shù),最優(yōu)狀態(tài)向量估計一般選用卡爾曼濾波或平滑[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 參數(shù)估計

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

(15)

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

(16)

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

(17)

其中Dk的計算公式如下

(18)

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

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

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

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

2 試驗和結果討論

2.1 濾波質量異常分析

美國得克薩斯大學空間研究中心目前只發(fā)布RL05球諧系數(shù)解的協(xié)方差矩陣,故沒有選用RL06數(shù)據(jù)進行處理。RL05數(shù)據(jù)可以從網(wǎng)站(http:∥download.csr.utexas.edu/outgoing/grace/)下載。該網(wǎng)站僅提供2002年4月至2014年6月的數(shù)據(jù),后續(xù)GRACE在軌運行期間的數(shù)據(jù)暫未提供。這135個月數(shù)據(jù)足以驗證SS-DDK方法的性能。試驗前先做如下預處理:從球諧系數(shù)中扣除該時段的平均值;用衛(wèi)星激光測距(satellite laser ranging,SLR)[24]結果替換球諧系數(shù)C20項;地心1階項由GRACE衛(wèi)星和海洋模型聯(lián)合估計得到;利用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],故將其作為參照數(shù)據(jù)。

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

(19)

ΔSlmsin(mλ))

(20)

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

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

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

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

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

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

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

(21)

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

ΔH=gTx

(22)

(23)

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

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

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

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

2.3 時間序列擬合分析

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

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

(24)

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

(25)

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

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

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

2.4 信號和噪聲水平分析

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

(26)

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

2.5 仿真試驗分析

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

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

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

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

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

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

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

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

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

3 結 論

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

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

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

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

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

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

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

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

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