趙元元，蘇宗躍，劉 元，穆大鵬，3

(1.華北地質(zhì)勘查局五一九大隊，河北保定071051;2.山東科技大學測繪科學與工程學院，山東青島266590;3.中國科學院遙感與數(shù)字地球研究所，北京100094)

一、引 言

GRACE采取衛(wèi)星跟蹤衛(wèi)星的觀測模式，能夠確定地球重力場的時變特征，其時間分辨率為10 d到30 d，空間分辨率為 300 km 到 400 km［1］。GRACE衛(wèi)星極大的提高了人類對地球表層質(zhì)量遷移和再分布的認識，尤其是兩極地區(qū)冰川消融［2］和大尺度流域水文變化［3］。

由于受到衛(wèi)星軌道誤差、觀測誤差、模型誤差以及數(shù)據(jù)處理造成的誤差等影響，使用GRACE數(shù)據(jù)研制的時變地球重力場模型含有較多噪聲，尤其是位系數(shù)的高階部分。為了抑制GRACE含有的噪聲，一般通過設(shè)計平滑核函數(shù)，對位系數(shù)采取某種限制措施，比如經(jīng)典的高斯平滑核函數(shù)［4］，通過降低高階位系數(shù)的權(quán)重，來降低GRACE的隨機噪聲;Fan濾波是在高斯濾波基礎(chǔ)上［5］，不僅降低高階位系數(shù)的權(quán)重，而且也降低高次位系數(shù)的權(quán)重，從而進一步降低GRACE的噪聲;DDK2濾波是通過貝葉斯估計［6］，利用GRACE兩顆衛(wèi)星的軌道誤差來設(shè)計誤差矩陣和模型信息設(shè)計信號矩陣，來提高位系數(shù)解的精度。GRACE濾波的方法還有很多，比如將衛(wèi)星觀測誤差和泄露誤差最小化來設(shè)計適合研究區(qū)域的平滑核函數(shù)［7］，基于高斯濾波的其他非各向同性濾波［8］，基于時間序列的統(tǒng)計方法［9］和經(jīng)驗正交函數(shù)方法［10］，以及被廣泛應(yīng)用的經(jīng)驗去相關(guān)濾波等［11］。

本文將使用高斯濾波、Fan濾波和DDK2濾波這3種方法來反演地表質(zhì)量變化，并分析比較它們的結(jié)果異同。GRACE數(shù)據(jù)使用的是由CSR(Center for Space Research)發(fā)布的Level-2 RL05版本，用戶可以從http:∥isdc.gfz-potsdam.de/grace進行下載，DDK2濾波的結(jié)果由Kushce提供，下載網(wǎng)址是http:∥icgem.gfz-potsdam.de/ICGEM/TimeSeries.html。

二、濾波方法

GRACE確定的時變重力場可以反演地球表層的質(zhì)量變化，由于其位系數(shù)含有較大噪聲，直接反演的結(jié)果很難分辨出需要的信號，需要對位系數(shù)進行某種限制。經(jīng)過平滑之后，地表的質(zhì)量異常可以由下式計算［12］

式中，θ和λ分別為地心緯度和地心經(jīng)度;Δσ表示質(zhì)量異常;a為地球平均半徑;ρave為地球平均密度(5517 kg/m3)，?Plm是規(guī)格化締合勒讓德函數(shù)，kl表示負荷勒夫數(shù)，Wlm為平滑核函數(shù)，ΔClm和ΔSlm為GRACE位系數(shù)與平均值的差值。Δσ除以水密度ρw(=1000 kg/m3)就得到了以等效水高表示的質(zhì)量異常。

對于平滑核函數(shù)Wlm而言，其不同的設(shè)計構(gòu)成了不同的濾波器。Wahr最早引入高斯平滑核函數(shù)來抑制高階位系數(shù)的噪聲，該平滑核函數(shù)僅對位系數(shù)的階部分起作用，即Wlm退化為

由于高斯濾波的實質(zhì)是對位系數(shù)的不同階賦予不同的權(quán)重，同一階下的不同次的權(quán)重一樣，文獻［5］在高斯濾波基礎(chǔ)上，構(gòu)造了Fan濾波

式中，Wl和Wm均為高斯平滑核函數(shù)，這樣Fan濾波對位系數(shù)的階和次同時起作用。

文獻［6］通過貝葉斯估計，使用GRACE兩顆衛(wèi)星的幾何軌道誤差來近似估計誤差矩陣E，使用模型確定的信息來估計GRACE的信號協(xié)方差矩陣S，這樣，平滑核函數(shù)可以由下式得到

式中，a表示正則化因子。

三、結(jié)果分析

本文使用CSR發(fā)布的Level-2 RL05數(shù)據(jù)，時間跨度為2003年1月到2012年12月，共120個月數(shù)據(jù)，缺失月份的數(shù)據(jù)通過線性內(nèi)插得到，利用這120個月的數(shù)據(jù)計算GRACE位系數(shù)的平均值，再計算每個月位系數(shù)的殘差。利用式(1)，使用高斯濾波、Fan濾波和DDK2濾波計算了2010年4月和10月的質(zhì)量異常，并轉(zhuǎn)化成等效水高值，其中高斯濾波和Fan濾波的半徑均為400 km。選擇4月份和10月份的原因是，在赤道南北兩側(cè)附近，這兩個月份的質(zhì)量異常變化較其他月份強烈，并且趨勢相反。

從等效水高圖看，在4月份，3種濾波方法均可以在亞馬遜流域、剛果河流域、澳洲北部以及南極和格陵蘭島等地區(qū)觀察到強烈的質(zhì)量變化信號。對于殘余的南北條帶狀噪聲，高斯濾波要明顯大于Fan濾波和DDK2濾波。Fan濾波不僅抑制高階的位系數(shù)，而且抑制高次的位系數(shù)，相對于高斯濾波，F(xiàn)an濾波能更有效的降低噪聲，但同時也會造成信號的損失;在10月份，3種濾波方法的南北條帶狀噪聲均要比各自4月份的結(jié)果要大，但Fan濾波和DDK2濾波的噪聲殘余仍然比高斯濾波小(見表1)。此外，在亞馬遜地區(qū)，F(xiàn)an濾波結(jié)果的振幅明顯要小于高斯濾波和DDK2濾波，這也說明Fan濾波雖然能有效抑制高階次的噪聲，但也會造成該部分信號的損失。

表1 3種濾波方法結(jié)果統(tǒng)計

3種濾波的統(tǒng)計結(jié)果也有較大差異，如圖1所示，高斯濾波和Fan濾波結(jié)果較為接近，后者的最大值、最小值以及均方根都要略小于前者，這也是前面指出的Fan濾波對于高階次位系數(shù)的抑制。對于DDK2濾波，在4月份，其最大值要高于高斯濾波48 mm，高于Fan濾波71 mm，而最小值差距則超過了260 mm，均方根也較前兩者有9 mm和12 mm的差距。在10月份，DDK2濾波與高斯濾波和Fan濾波的結(jié)果差異進一步拉大，其最小值是Fan濾波的2倍，而均方根的差距則達到了13 mm和24 mm。造成這種較大差異的原因在于，DDK2濾波使用的是貝葉斯估計方法，依賴于使用的先驗信息，也就是GRACE衛(wèi)星的幾何軌道誤差矩陣和模型確定的信號矩陣，它在有效降低南北條帶噪聲的同時，也可能會造成信號的失真。

圖1 2010年等效水高圖

四、結(jié) 論

本文使用高斯濾波、Fan濾波和DDK2濾波反演了2010年4月和10月GRACE時變重力場模型確定的地表質(zhì)量異常，并作了比較分析。高斯濾波只對GRACE位系數(shù)的不同階起降權(quán)作用，也就是同一階里面的不同次的權(quán)重一樣;Fan濾波的構(gòu)造是在高斯濾波基礎(chǔ)上，它不僅對位系數(shù)的階起作用，而且對次也同時起作用，進一步壓縮GRACE高階次的噪聲;DDK2濾波是基于貝葉斯估計，使用先驗信息構(gòu)造平滑核函數(shù)，其濾波的結(jié)果依賴于先驗信息的近似程度，文獻［6］使用GRACE兩顆衛(wèi)星的軌道誤差來設(shè)計誤差矩陣，使用模型來計算信號矩陣，有效地降低了南北條帶噪聲。

對于4月份和10份的反演結(jié)果，從等效水高圖和統(tǒng)計結(jié)果看，不同方法之間有著較大差異:Fan濾波和DDK2濾波之后的南北噪聲殘余要小于高斯濾波，3種方法的最大值、最小值和均方根也有著較大區(qū)別，尤其是DDK2與前兩者之間;此外，對于4月份和10月份這兩個不同的時間，3種濾波方法的對于南北條帶噪聲的抑制作用也很不同，4月份要明顯優(yōu)于10月份，可能原因是GRACE軌道運行以及地面跟蹤、GPS跟蹤在不同時間的誤差水平不同。

［1］TAPLEY B D，BETTADPUR S，RIES J C，et al.GRACE Measurements of Mass Variability in the Earth System［J］.Science，2004，305(5683):503-505.

［2］LUTHCKE SB，ZWALLY H J，ABDALATI W，et al.Recent Greenland Ice Mass Loss by Drainage System from Gravity Observations［J］.Science，2006，314(5803):1286-1289.

［3］馮偉，萊莫尼 JM，鐘敏，等.利用重力衛(wèi)星GRACE監(jiān)測亞馬遜流域2002-2010年的陸地水變化［J］.地球物理學報，2012，55(3):814-821.

［4］詹金剛，王勇，郝曉光.GRACE時變重力位系數(shù)誤差的改進去相關(guān)算法［J］.測繪學報，2011，40(4):442-446.

［5］ZHANG Z Z，CHAO B F，LU Y，et al.An Effective Filtering for GRACE Time-variable Gravity:Fan Filter［J］.Geophysical Research Letters，2009(36):L17311.

［6］KUSCHE J.Approximate De-correlation and Non-isotropic Smoothing of Time-variable GRACE-type Gravity Field Models［J］.Journal of Geodesy，2007(81):733-749.

［7］SWENSON S，WAHR J.Methods for Inferring Regional Surface-mass Anomalies form Gravity Recovery and Climate Experiment(GRACE)Measurements of Time-variable Gravity［J］.Journal of Geophysical Research，2002，107(B9):109-117.

［8］HAN SC，SHUM CK，JEKELI C，et al.Non-isotropic Filtering of GRACE Temporal Gravity for Geophysical Signal Enhancement［J］.Geophysical Journal International，2005，163(1):18-25.

［9］DAVISJM，TAMISIEA E，ELOSEGUI P，et al.A Statistical Filtering Approach for Gravity Recovery and Climate Experiment(GRACE)Gravity Data［J］.Journal of Geophysical Research，2008，113(01):117-202.

［10］WOUTERS B，SCHRAMA E O.Improved Accuracy of GRACE Gravity Solutions through Empirical Orthogonal Function Filter of Spherical Harmonics［J］.Geophysical Research Letters，2007(34):711-715.

［11］SWENSON S，WAHR J.Post-processing Removal of Correlated Errors in GRACE Data［J］.Geophysical Research Letters，2006(33):553-561.

［12］WAHR J，MOLENAAR M，BRYAN F.Time Variability of the Earth’s Gravity Filed:Hydrological and Oceanic Effects and Their Possible Detection Using GRACE［J］.Journal Geophysical Research，1998，103(B12):30205-30229.