詹金剛，王 勇，郝曉光

1.中國科學院測量與地球物理研究所，湖北武漢430077；2.中國科學院研究生院，北京100049

1 引 言

重力場及其變化反映了地球表層和內(nèi)部的物質(zhì)密度分布及運動狀態(tài)。因此，根據(jù)重力場的時空變化，可以推演和監(jiān)測地球系統(tǒng)物質(zhì)運動和交換過程。過去，人們主要通過地面重復重力測量獲得重力場隨時間的變化信息。由于地面連續(xù)重力觀測站點相對較少，高精度的重復重力觀測又較耗時，難以獲得連續(xù)的時空重力場信息。2002年GRACE的成功運行使得研究地球重力場隨時間的變化信息成為現(xiàn)實。起初，由于GRACE資料解算得到的時變重力場的精度較差以及資料積累的周期較短，主要利用GRACE資料驗證和研究大尺度上的重力季節(jié)變化［1－9］，國內(nèi)也有一些學者利用GRACE資料對重力場的恢復方法、仿真及在海平面變化中的應用等方面進行研究并取得一定成果［10－12］。隨著重力場恢復方法的改進和濾波技術的進一步提高，目前由一個月的GRACE重力衛(wèi)星資料反演大地水準面的精度在400km尺度上已達到1～2cm［13－14］。

針對GRACE衛(wèi)星資料，目前展開的工作主要包含兩方面：一是利用GRACE資料解算時變重力場模型，二是時變重力場的應用。由于目前解算的GRACE時變重力場模型的高階系數(shù)存在較大的誤差，由模型系數(shù)恢復地球時變重力場結果中表現(xiàn)為嚴重的條帶誤差，給地球物理信號的識別帶來一定的困難，因此必須進行濾波提高信噪比。針對消除條帶誤差常用的處理方法，按照其濾波思想主要分為兩類：一類是通過引入濾波因子，降低模型高階系數(shù)的權重，從而達到平滑的效果，這類方法主要包括各向同性和非各向同性的高斯濾波、維納濾波以及扇形濾波［15－18］，其實質(zhì)是以犧牲模型的空間分辨率為代價來換取空間平滑的效果，濾波半徑越大，模型空間分辨率的犧牲程度也越大；另一類方法是通過消除模型系數(shù)之間的相關誤差，從而達到去條帶誤差的效果，這類方法主要包括文獻［19］的基于階l的變滑動窗多項式擬合去相關誤差方法，文獻［20—21］提出的多項式擬合去相關誤差方法以及文獻［22］提出的基于階l次m的變滑動窗多項式擬合去相關誤差方法。其中尤以文獻［19］提出的滑動窗多項式擬合去相關誤差方法效果最為顯著，文獻［22］的濾波方法計算過于復雜，而濾波效果和文獻［19］的滑動窗多項式擬合去相關誤差方法并無顯著差異。這類方法的優(yōu)點是直接消除系數(shù)間的相關誤差，而不會犧牲模型的空間分辨率，其缺點是在赤道附近區(qū)域去條帶誤差的效果并不顯著?；谏鲜鲈颍ǔ２捎脤深悶V波技術相結合的組合濾波方法。

第二類去相關誤差方法效果的優(yōu)劣，直接決定著高斯濾波半徑選取的大小。去相關誤差效果越顯著，則高斯平滑半徑的選取就越小，對模型空間分辨率的損失也就越?。环粗?，則對模型空間分辨率的損失越大。文獻［19］提出的變滑動窗多項式擬合去相關誤差方法在赤道兩側區(qū)域取得了很好的去條帶誤差效果，但因其在文獻中并沒有公開具體的實現(xiàn)步驟，因而許多學者按照其思想濾波時并未實現(xiàn)其文獻中的理想效果，目前國際上根據(jù)其思想常用的做法是采用多項式擬合去相關技術，而沒有采用滑動窗［20－21］。根據(jù)文獻［19］的思想，針對滑動窗的特點，詳細分析了滑動窗在時變重力場數(shù)據(jù)處理中的不足，提出滑動窗多項式擬合去相關誤差的改進方法。將改進方法應用于GRACE時變重力場數(shù)據(jù)處理時，在赤道兩側區(qū)域取得了較為顯著的去條帶濾波效果。

2 滑動窗多項式擬合去相關誤差方法的改進

在抑制條帶噪聲方面，文獻［19］提出的滑動窗多項式擬合去相關方法的思想是：固定時變重力場模型的次m，對模型系數(shù)的偶數(shù)階系數(shù)序列（l＝0，2，4，…）以及奇數(shù)階系數(shù)序列（l＝1，3，5，…）分別采用滑動窗多項式擬合，并以該擬合值作為相應系數(shù)的相關誤差改正值。這里以窗口寬度為7點，m＝11和45時為例說明（假定模型最大展開為60階）：當m＝11時，系數(shù)序列為C11，11、C13，11、C15，11、C17，11、C19，11、…、C51，11、C53，11、C55，11、C57，11、C59，11。由于窗口的寬度為7點，按照滑動窗的基本原理，得到去相關誤差改正的系數(shù)依次是C17，11、C19，11、…、C51，11、C53，11，這樣數(shù)據(jù)序列的兩端各有三個系數(shù)沒有得到去相關誤差改正。當m＝45時，系數(shù)序列C45，45、C47，45、C49，45、C51，45、C53，45、C55，45、C57，45、C59，45中僅有C51，45、C53，45兩個系數(shù)得到去相關誤差改正，而數(shù)據(jù)序列兩端的其他6個系數(shù)并沒有得到去相關誤差改正。這種情況下，滑動窗多項式擬合去相關誤差改正方法實際上已經(jīng)沒有多大意義，因為沒有得到誤差改正或者說舍棄的系數(shù)數(shù)量已遠大于改正系數(shù)的數(shù)量。顯然滑動窗去條帶誤差方法在m＞45時，已經(jīng)失去了去相關誤差改正的作用和意義。

為減少兩端數(shù)據(jù)的舍棄，同時使得滑動窗多項式擬合去相關誤差方法在模型系數(shù)的更高次數(shù)m上得到應用，本文根據(jù)相鄰高階系數(shù)間具有相關誤差這一特點，在數(shù)據(jù)處理時，先對數(shù)據(jù)序列的兩端做反向邊界延拓，向邊界兩端各延伸1／2窗寬的數(shù)據(jù)。仍然以窗口寬度為7點，m＝11為例說明（假定模型最大展開為60階），經(jīng)反向邊界延拓后，系數(shù)序列變?yōu)椋瑿19，11、－C15，11、－C13，11、C11，11、C13，11、C15，11、C17，11、C19，11、…、C51，11、C53，11、C55，11、C57，11、C59，11、－C57，11、－C55，11、－C53，11。這樣，在計算C13，11系數(shù)的誤差改正時，用到了－C15，11、－C13，11、C11，11、C13，11、C15，11、C17，11、C19，117個系數(shù)間的相關性。這樣處理后，不但使得原數(shù)據(jù)序列兩端的數(shù)據(jù)得到滑動窗多項式擬合去相關誤差改正，同時也使得滑動窗去相關誤差方法可以在更高階次系數(shù)中得到應用。圖1是滑動窗多項式擬合去相關誤差方法改進前后的系數(shù)變化曲線。圖1（a）與圖1（b）為普通滑動窗結果，（a）為每階系數(shù)變化曲線，（b）為奇偶階系數(shù)變化曲線；圖1（c）與圖1（d）為改進后的滑動窗結果，（c）為每階系數(shù)變化曲線，（d）為奇偶階系數(shù)變化曲線。從圖1中可以明顯看出，模型系數(shù)經(jīng)原滑動窗多項式擬合去相關誤差方法改正后，高階次的系數(shù)仍然存在較大的誤差，表現(xiàn)為高階次的系數(shù)變化不收斂，例如m＝15、18時尤為明顯；圖1（c）、（d）顯示，模型系數(shù)經(jīng)改進后的滑動窗多項式擬合去相關誤差方法處理后，模型系數(shù)隨著階l的增大趨于穩(wěn)定和收斂，表明系數(shù)受誤差的影響較小。例如當m＝15時，經(jīng)改進后的滑動窗多項式擬合去相關誤差方法處理后，高階次系數(shù)誤差得到顯著改善，曲線變化的幅度減小了一個量級。

圖1 去相關誤差方法改進前后系數(shù)變化曲線對比Fig.1 Coefficients change before and after improvement of de－correlation method

3 效果對比與檢驗

為檢驗改進后數(shù)據(jù)處理方法的效果，以CSR RL04數(shù)據(jù)為例進行驗證。由于在RL04數(shù)據(jù)的解算中，引入了更高精度的地球重力場背景模型、海潮模型以及極潮模型等，使得模型的低階次系數(shù)的精度得到進一步提高，系數(shù)變化曲線顯示系數(shù)間相關誤差在m≥15時開始越來越明顯；其次，CSR RL04數(shù)據(jù)的階能譜顯示模型系數(shù)約在l≥15時，系數(shù)誤差逐漸震蕩放大。綜合以上因素，本文選取在m≥15時，開始做去相關誤差改正，而對于低階系數(shù)盡量不做改動。文獻［19］和文獻［22］中，對低階次系數(shù)去相關誤差改正采用了較大的窗口寬度。這主要是因為模型低階次系數(shù)間的相關誤差小或者說其精度相對較高，利用大的窗寬進行弱相關的多項式擬合，以減小對低階次系數(shù)的改正量。由于本文采用精度較高的RL04數(shù)據(jù)且起算階數(shù)取為15，所以采用固定的窗口寬度。大量數(shù)據(jù)試算后發(fā)現(xiàn)，窗口寬度越大，去條帶誤差效果越差，表現(xiàn)為結果中的條帶噪聲越明顯，而當窗口寬度為5時，計算結果過于平滑，有效信號的幅度削弱過多。因此，本文計算中選擇窗口寬度為7點，多項式擬合次數(shù)為3。

圖2為隨機抽取的2008年1～4月GRACE時變重力場模型計算結果。其中，圖2（a）為沒有進行任何濾波的月重力異常結果，從圖中可以看出月重力異常結果主要表現(xiàn)為條帶噪聲。圖2（b）、（c）分別為滑動窗去相關誤差方法改進前后的重力異常圖。從圖2（b）、（c）中可以明顯看出，改進后的方法在赤道兩側區(qū)域抑制條帶噪聲的效果得到顯著改善，說明了新方法在赤道兩側區(qū)域抑制條帶誤差的有效性。此外，對2003年1月至2008年12月共72個月的GRACE數(shù)據(jù)進行了驗證對比和統(tǒng)計。統(tǒng)計結果表明：新方法在赤道兩側區(qū)域抑制條帶噪聲的效果較原方法均具有顯著的改善。這一結果有效證明新方法在GRACE數(shù)據(jù)處理中抑制條帶誤差的有效性。

圖2結果說明了滑動窗去相關誤差的改進方法在抑制條帶誤差方面的效果及有效性。在證明新方法去條帶誤差的有效性之后，通常還需要驗證該方法的正確性，即新方法是否會對真實信號產(chǎn)生扭曲以及是否會產(chǎn)生虛假信號。為此，采用全球陸地資料同化系統(tǒng)GLDAS土壤濕度模型數(shù)據(jù)作為驗證。GLDAS模型同化了四個不同的全球水文模型，并采用NASA新一代的空間對地觀測技術和地面數(shù)據(jù)來約束地球表面的狀態(tài)，進而獲得地球表面的近實時信息，是目前最好的全球水文模型之一［23］。土壤濕度變化對全球重力場變化趨勢的影響如圖3所示，其中圖3（a）為沒有采用任何濾波器的結果，圖3（b）為只采用滑動窗去相關誤差改進方法的濾波結果。從圖3中可以看出，滑動窗去相關誤差改進方法除在高緯度地區(qū)對信號的強度或振幅稍微有削弱外，并沒有改變信號的位置和形狀，也沒有產(chǎn)生虛假信號。圖3結果證明本文提出的滑動窗去相關誤差改進方法的可靠性和正確性。

圖2 GRACE月重力異常圖Fig.2 Maps of monthly anomaly of GRACE gravity field（1Gal＝1cm／s2）

圖3 去相關誤差改進方法對重力場變化趨勢的影響Fig.3 Effects of improved de－correlation method on gravity trend

4 結 論

（1）在滑動窗多項式擬合數(shù)據(jù)處理技術上做了相應改進，改進后的數(shù)據(jù)處理方法在不影響原滑動窗去相關誤差方法的前提下，不僅使得數(shù)據(jù)序列兩端的系數(shù)得到去相關誤差改正，而且使得滑動窗去條帶誤差技術能夠應用到模型的更高階次。將改進后的滑動窗多項式擬合去條帶誤差方法應用于CSR解算的GRACE數(shù)據(jù)時，在抑制條帶噪聲的效果和誤差的改善范圍上較改進前方法結果均有顯著提高，表明新方法對于消除模型系數(shù)間的相關誤差具有顯著效果。

（2）將改進方法應用于全球土壤濕度變化數(shù)據(jù)，計算土壤濕度變化對全球重力場變化趨勢的影響。計算結果表明，本文提出的改進方法僅對高緯度區(qū)域內(nèi)信號的強度即振幅變化有稍微消弱，而對信號的位置和形狀并沒有明顯的改變，也沒有產(chǎn)生虛假信號。

［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］ VELICOGNA I，WAHR J.Measurements of Time－variable Gravity Show Mass Loss in Antarctica［J］.Science，2006，311（5768）：1745－1756.

［3］ HU Xiaogong，CHEN Jianli，ZHOU Yonghong，et al.Seasonal Water Storage Change of the Yangtze River Basin Detected by GRACE［J］.Science in China：Series D Earth Sciences，2006，49（5）：483－491.（胡小工，陳劍利，周永宏，等.利用GRACE空間重力測量監(jiān)測長江流域水儲量的季節(jié)性變化［J］.中國科學：D輯地球科學，2006，36（3）：225－232.）

［4］ ZHONG Min，DUAN Jianbin，XU Houze，et al.Trend of China Land Water Storage Redistribution at Medi－and Large－spatial Scales in Recent Five Years by Satellite Gravity Observations［J］.Chinese Science Bulletin，2009，54（5）：816－821.（鐘敏，段建賓，許厚澤，等.利用衛(wèi)星重力觀測研究近5年中國陸地水量中長空間尺度的變化趨勢［J］.科學通報，2009，54（9）：1290－1294.）

［5］ WANG Hansheng，WANG Zhiyong，YUAN Xudong，et al.Water Storage Changes in Three Gorges Water Systems Area Inferred from GRACE Time－variable Gravity Data［J］.Chinese Journal of Geophysics，2007，50（3）：730－736.（汪漢勝，王志勇，袁旭東，等.基于GRACE時變重力場的三峽水庫補給水系水儲量變化［J］.地球物理學報，2007，50（3）：730－736.）

［6］ E Dongchen，Yang Yuangde，Chao Dingbo.The Sea Level Change from the Antarctic Ice Sheet Based on GRACE［J］.Chinese Journal of Geophysics，2009，52（9）：2222－2228.（鄂棟臣，楊元德，晁定波.基于GRACE資料研究南極冰蓋消減對海平面的影響［J］.地球物理學報，2009，52（9）：2222－2228.）

［7］ XIAO Yun，XIA Zheren，WANG Xingtao.Recovering the Earth Gravity from Inter－satellite Range－rate of GRACE［J］.Acta Geodaetica et Cartographica Sinica，2007，36（1）：19－25.（肖云，夏哲仁，王興濤.用GRACE星間速度恢復地球重力場［J］.測繪學報，2007，36（1）：19－25.）

［8］ ZHOU Xuhua，WU Bin，XU Houze，et al.Resolution Estimation of Earth Gravity Field Recovery through the Low－low Satellite to Satellite Technology by Numerical Simulation［J］.Chinese Journal of Geophysics，2005，48（2）：282－287.（周旭華，吳斌，許厚澤，等.數(shù)值模擬估算低低衛(wèi)－衛(wèi)跟蹤觀測技術反演地球重力場的空間分辨率［J］.地球物理學報，2005，48（2）：282－287.）

［9］ ZHOU Xuhua，XU Houze，WU Bin，et al.Earth’s Gravity Field Derived from GRACE Satellite Tracking Data［J］.Chinese Journal of Geophysics，2006，49（3）：718－723.（周旭華，許厚澤，吳斌，等.用GRACE衛(wèi)星跟蹤數(shù)據(jù)反演地球重力場［J］.地球物理學報，2006，49（3）：718－723.）

［10］ ZOU Xiancai，LI Jiancheng，JIANG Weiping，et al.Research on the Simultaneous Solution Method for Satellite Gravity Data Analysis and Its Simulation［J］.Acta Geodaetica et Cartographica Sinica，2010，39（4）：344－348.（鄒賢才，李建成，姜衛(wèi)平，等.衛(wèi)星重力資料分析的同解法研究及其仿真［J］.測繪學報，2010，39（4）：344－348.）

［11］ JIANG Tao，LI Jiancheng，WANG Zhengtao，et al.Global Sea Level Variations from Combined Jason－1and GRACE Data［J］.Acta Geodaetica et Cartographica Sinica，2010，39（2）：135－140.（蔣濤，李建成，王正濤，等.聯(lián)合Jason－1與GRACE衛(wèi)星數(shù)據(jù)研究全球海平面變化［J］.測繪學報，2010，39（2）：135－140.）

［12］ XIAO Yun，XIA Zheren，WANG Xingtao.Recovering the Earth Gravity Field from Inter－satellite Range－rate of GRACE［J］.Acta Geodaetica et Cartographica Sinica，2007，36（1）：19－25.（肖云，夏哲仁，王興濤.用GRACE星間速度恢復地球重力場［J］.測繪學報，2007，36（1）：19－25.）

［13］ BETTADPUR S.Level－2Gravity Field Product User Handbook［EB／OL］.Austin，Texas：University of Texas.2003－12－01［2010－03－10］.http：∥www.csr.utexas.edu／grace／publications／handbook／L2－User－Handbook＿v1.0.pdf.

［14］ TAPLEY B D，BETTADPUR S，WATKINS M M，et al.The Gravity Recovery and Climate Experiment：Mission Overview and Early Results［J／OL］.Geophy Research Letters，2004，31：1－4［2010－03－10］.http：∥www.spaceweather.ac.cn／publication／jgrs／2004／Geophysical% 20Research%20Letters／may／2004GL019920.pdf.

［15］ WAHR J，MOLENAAR M，BRYAN F.Time Variability of the Earth’s Gravity Field：Hydrological and Oceanic Effects and Their Possible Detection Using GRACE［J］.Journal of Geophysical Research，1998，103（B12）：30205－30229.

［16］ HAN S C，SHUM C K，JEKELI C，et al.（2005），Nonisotropic Filtering of GRACE Temporal Gravity for Geophysical—Signal Enhancement［J］.Geophysical Journal International，2005，163（9）：18－25.

［17］ SASGEN I，MARTINEC Z，F(xiàn)LEMING K.Wiener Optimal Filtering of GRACE Data［J］.Studia Geophysica et Geodaetica，2006，50（4）：499－508.

［18］ ZHANG Z H，CHAO B F，LU Y，et al.An Effective Filtering for GRACE Time－variable Gravity：Fan Filter［J］.Geophysical Research Letters，2009，36：1－6.

［19］ SWENSON S，WAHR J.Post－processing Removal of Correlated Errors in GRACE Data［J］.Geophysical Research Letters，2006，33：1－6.

［20］ CHAMBERS D P.Evaluation of New GRACE Timevariable Gravity Data over the Ocean［J］.Geophysical Research Letters，2006，33：12－16.

［21］ CHEN J L，WILSON C R，TAPLEY B D，et al.GRACE Detects Coseismic and Postseismic Deformation from the Sumatra－Andaman Earthquake［J］.Geophysical Research Letters，2007，34：33－37.

［22］ DUAN X J，GUO J Y，SHUM C K.et al.On the Post－processing Removal of Correlated Errors in GRACE Temporal Gravity Field Solutions［J］.Journal of Geodesy，2009，83（11）：1095－1106.

［23］ RODELL M，HOUSER P R，JAMBOR U.et al.The Global Land Data Assimilation System［J］.Bulletin of American Meteorological Society，2004，85：381－394.