詹金剛，王 勇，郝曉光

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

1 引 言

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

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

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

2 滑動窗多項(xiàng)式擬合去相關(guān)誤差方法的改進(jìn)

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

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

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

3 效果對比與檢驗(yàn)

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

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

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

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

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

4 結(jié) 論

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

（2）將改進(jìn)方法應(yīng)用于全球土壤濕度變化數(shù)據(jù)，計(jì)算土壤濕度變化對全球重力場變化趨勢的影響。計(jì)算結(jié)果表明，本文提出的改進(jìn)方法僅對高緯度區(qū)域內(nèi)信號的強(qiáng)度即振幅變化有稍微消弱，而對信號的位置和形狀并沒有明顯的改變，也沒有產(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］.中國科學(xué)：D輯地球科學(xué)，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］.科學(xué)通報(bào)，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時(shí)變重力場的三峽水庫補(bǔ)給水系水儲量變化［J］.地球物理學(xué)報(bào)，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］.地球物理學(xué)報(bào)，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星間速度恢復(fù)地球重力場［J］.測繪學(xué)報(bào)，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)跟蹤觀測技術(shù)反演地球重力場的空間分辨率［J］.地球物理學(xué)報(bào)，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］.地球物理學(xué)報(bào)，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］.測繪學(xué)報(bào)，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］.測繪學(xué)報(bào)，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星間速度恢復(fù)地球重力場［J］.測繪學(xué)報(bào)，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.