郭飛霄，苗岳旺，肖 云

1.西安測繪研究所，陜西 西安，710054；2.地理信息工程國家重點實驗室，陜西 西安，710054；3.測繪信息技術(shù)總站，陜西 西安，710054

?

使用GRACE時變重力場中的維納濾波

郭飛霄1,2，苗岳旺3，肖 云1,2

1.西安測繪研究所，陜西 西安，710054；2.地理信息工程國家重點實驗室，陜西 西安，710054；3.測繪信息技術(shù)總站，陜西 西安，710054

GRACE時變重力場模型誤差隨著階數(shù)增大而迅速增大，必須采用濾波進(jìn)行空間平滑。本文利用GRACE時變重力場模型數(shù)據(jù)計算全球陸地水儲量變化，分析了維納濾波對反演結(jié)果的影響，并與不同半徑高斯濾波反演結(jié)果進(jìn)行了對比。實驗結(jié)果表明：維納濾波是一種低通濾波，濾波系數(shù)僅與時變重力場模型噪聲特性相關(guān)；維納濾波能夠有效消除全球陸地水儲量變化反演結(jié)果中的南北條帶誤差，反演結(jié)果與半徑400km高斯濾波反演結(jié)果相當(dāng)。

GRACE衛(wèi)星，時變重力場，維納濾波，水儲量變化，衛(wèi)星重力

1 引 言

地球重力場及其時變效應(yīng)反映了地球表層及內(nèi)部的物質(zhì)分布和變化。當(dāng)前，全球性環(huán)境問題如海平面上升、極地冰川融化等與地球表層的物質(zhì)遷移緊密相關(guān)。因此，研究地球系統(tǒng)的質(zhì)量遷移和重新分布對監(jiān)測全球環(huán)境和氣候變化具有重要意義。2002年3月，GRACE衛(wèi)星的成功發(fā)射開創(chuàng)了高精度全球重力場觀測與氣候變化試驗新紀(jì)元，該計劃旨在恢復(fù)高精度和高分辨率全球靜態(tài)重力場，并估計地球重力場的時變特征。GRACE衛(wèi)星Level-2數(shù)據(jù)能夠提供空間分辨率約為400km、時間分辨率約為30天的地球重力場時變模型。Wahr等建立了利用時變重力場模型估計地球表層質(zhì)量變化的球諧系數(shù)算法[1]，目前已廣泛應(yīng)用于[2-5]陸地水儲量變化、地震同震、海洋環(huán)流和冰川融化等領(lǐng)域的研究。

GRACE時變重力場模型受衛(wèi)星軌道誤差、星間測距誤差、加速度計誤差影響以及高階項誤差等影響，球諧系數(shù)法恢復(fù)地球時變重力場結(jié)果表現(xiàn)為嚴(yán)重南北方向條帶誤差，因此必須進(jìn)行濾波處理。Wahr等提出了各向同性的高斯濾波方法，對時變重力場模型階項進(jìn)行降權(quán)濾波[1]；張子占等提出了扇形濾波方法，該方法對時變重力場模型的階項和次項同時應(yīng)用一次高斯濾波，其本質(zhì)是一種非各向同性的高斯濾波[6]。但不論是高斯濾波還是扇形濾波，計算時都需要選取合適的濾波半徑。如果濾波半徑過大，對模型空間分辨率的犧牲程度也越大；濾波半徑過小，則無法有效消除條帶誤差，給地球物理信號識別帶來困難。因此，對于GRACE時變重力場空間濾波，要在空間分辨率和噪聲減小之間進(jìn)行優(yōu)化處理。Sasgen等提出了針對時變重力場的維納濾波算法，與高斯濾波不同的是，維納濾波僅與期望信號和噪聲的階方差有關(guān)，能夠提高信噪比[7]。本文對GRACE時變重力場維納濾波算法進(jìn)行了推導(dǎo)分析，將維納濾波應(yīng)用球諧系數(shù)法，反演了全球陸地水儲量變化，并與采用高斯濾波方法反演結(jié)果進(jìn)行了對比分析，總結(jié)分析了維納濾波算法的特點。

2 維納濾波基本原理

維納濾波是線性卷積濾波，其輸出信號y(Ω)由實際測量信號x(Ω)和濾波函數(shù)h(Ω)進(jìn)行空間卷積得到：

(1)

(2)

(3)

Ylm(Ω)=Plm(cosθ)eimλ

(4)

(5)

維納濾波滿足以下條件：

(1)實際輸出信號y(Ω)與期望輸出信號y′(Ω)滿足最小二乘準(zhǔn)則，即：

(6)

由Parseval等式，式(6)可寫為：

(7)

(2)測量信號x(Ω)由真實信號s(Ω)和噪聲n(Ω)組成：

x(Ω)=s(Ω)+n(Ω)

(8)

(3)信號s(Ω)和n(Ω)不相關(guān)，即：

(9)

(4)期望輸出信號y′(Ω)與未污染的原始信號s(Ω)等價。

在以上條件下，最小二乘準(zhǔn)則E2可寫成如下形式：

(10)

(11)

式中，{Cs,lm,Ss,lm}和{Cn,lm,Sn,lm}分別為信號s(Ω)和噪聲n(Ω)的球諧模型展開系數(shù)。因此，式(10)就是求解hl使得E2最小。令E2對hl求偏導(dǎo)數(shù)，并令偏導(dǎo)數(shù)等于0，可得：

(12)

hl即為第l階維納濾波系數(shù)。從式(12)可以看出，當(dāng)噪聲可以忽略時，hl為單位值1；當(dāng)噪聲占主要成分時，hl為0。因此，式(12)給出了在兩個極端情況下的最佳過渡。

3 利用GRACE時變重力場數(shù)據(jù)確定維納濾波系數(shù)

維納濾波是各向同性濾波，濾波系數(shù)僅與信號和噪聲的階方差相關(guān)。利用式(12)求解維納濾波系數(shù)，需要分別估計重力信號和噪聲的階方差。但是，如果沒有額外信息或者假設(shè)條件，便無法從測量信號中估計這兩個值。不過，從GRACE數(shù)據(jù)可以獲取額外信息，來幫助估計重力信號和噪聲的階方差。

圖1 GRACE時變重力場模型平均階方差

圖1所示為GRACE時變重力場模型的平均階方差，以大地水準(zhǔn)面高形式表示，所示結(jié)果為2004年1月至2010年12月時變重力場模型統(tǒng)計結(jié)果。平均階方差計算如下：

(13)

(14)

Sasgen等認(rèn)為時變重力場模型的低階項(l< 20)主要是地幔以下地球重力場時變信號的體現(xiàn)，而高階項(l> 30)則反映了GRACE數(shù)據(jù)的噪聲信息；位于這兩者之間(20≤l≤30 )則是信號和噪聲的混合，體現(xiàn)了地球表層的陸地水儲量和海洋質(zhì)量的變化[7-8]。地球重力場時變信號可由式(15)所示的函數(shù)近似表達(dá)[7]：

(15)

(16)

對式(15)兩邊同時取以10為底數(shù)的對數(shù)，如式(16)所示。由于低階項部分(l<20)主要是地球重力場時變信號的體現(xiàn)，因此，參數(shù)a與參數(shù)b可利用平均階方差低階項的數(shù)據(jù)進(jìn)行最小二乘估計求得。對高階項部分，噪聲信息占主要部分。噪聲部分在以10為底數(shù)的對數(shù)形式下，可用線性函數(shù)[7]近似表達(dá)：

(17)

4 實驗分析

實驗中，GRACE時變重力場模型數(shù)據(jù)采用德國地學(xué)中心(GFZ)發(fā)布的Level-2Realease05數(shù)據(jù)，時間跨度為2004年1月至2010年12月。計算時對GRACE時變重力場模型截斷至60階，低階帶諧項C20項采用SLR觀測得到的C20項進(jìn)行替代[9]。先將所有GRACE時變重力場模型求平均值作為穩(wěn)態(tài)背景場，再計算每個月時變重力場模型相對于穩(wěn)態(tài)背景場的改變量，分別采用維納濾波、高斯濾波，按照式(18)計算全球陸地水儲量變化。

(18)

其中，Δh為陸地水儲量變化的等效水柱高度，{ΔClm,ΔSlm}為球諧系數(shù)的改變量，kl為l階負(fù)荷勒夫數(shù)，a為地球平均半徑，ρa(bǔ)為地球平均密度，ρw為水密度，Wl為濾波系數(shù)。

高斯濾波的濾波半徑分別取300km、400km和500km，濾波系數(shù)Wn按式(19)計算。其中，b=ln2/[1-cos(r/a)]，a為地球平均半徑，r為濾波半徑。

(19)

由GRACE時變重力場模型數(shù)據(jù)對時變重力場信號和噪聲部分按照第3節(jié)所述方法進(jìn)行擬合，所得結(jié)果如圖2所示。圖2中，紅色線為平均階方差，藍(lán)色線為噪聲擬合值，綠色線為時變重力場信號擬合值。所得擬合參數(shù)分別為：a=1.291，b=2.382，c=-3.347，d=0.044。

圖2 信號和噪聲階方差譜(紅色為時變重力場模型平均階方差，藍(lán)色為噪聲部分?jǐn)M合值，綠色為時變重力場信號擬合值)

下面以2010年4月全球陸地水儲量變化反演結(jié)果進(jìn)行分析，如圖3所示。其中，圖3(1)～(4)為采用維納濾波、半徑分別為300km、400km和500km高斯濾波的全球陸地水儲量變化反演結(jié)果。從圖3可以看出，采用半徑300km高斯濾波，雖然大尺度的全球陸地水儲量變化信號明顯，但反演結(jié)果中南北方向條帶誤差也非常明顯；而采用維納濾波和半徑400km高斯濾波，噪聲誤差都得到了有效控制，反演結(jié)果中對南北方向條帶誤差抑制效果顯著，大尺度的水文信號更加清晰；半徑500km高斯濾波反演結(jié)果中，南北方向條帶狀誤差基本消除，但反演結(jié)果分辨率有所下降，局部區(qū)域陸地水儲量變化信號減弱。進(jìn)一步對比圖3(1)和圖3(3)可以看出，維納濾波和半徑400km高斯濾波都基本消除了條帶誤差，反演結(jié)果吻合程度較高，但維納濾波反演陸地水儲量變化信號更強(qiáng)。

圖3 2010年4月全球陸地水儲量變化(單位：cm)

圖4所示為不同濾波半徑高斯濾波和維納濾波系數(shù)的比較結(jié)果。從圖4可以看出：(1)隨著濾波半徑增加，對于同階項高斯濾波系數(shù)更小，并且高斯濾波系數(shù)收斂速度更快；(2)維納濾波和高斯濾波都是低通濾波，低階項所占比重大，隨著階數(shù)的增加，濾波系數(shù)迅速減??；(3)與高斯濾波相比，維納濾波在低階項部分(l<20)的濾波系數(shù)權(quán)值較大，在高階項(l>30)的濾波系數(shù)迅速減小收斂。而GRACE時變重力場模型低階項主要是重力信號，高階項部分則主要反映了噪聲誤差。因此，與高斯濾波相比，維納濾波能夠更好地反映時變重力信號，這與圖3所示全球陸地水儲量反演結(jié)果相符。

圖4 不同半徑高斯濾波和維納濾波系數(shù)

5 結(jié) 論

由于GRACE時變重力場模型高階項部分存在較大誤差，利用該時變重力場模型反演陸地水儲量變化時表現(xiàn)為嚴(yán)重的條帶誤差，因此必須采用濾波進(jìn)行空間平滑。維納濾波是基于信號和噪聲的平均階方差譜建立信號和噪聲函數(shù)，其提高了低階項權(quán)重，同時降低了高階項權(quán)重。由實驗結(jié)果可得出以下結(jié)論：(1)維納濾波是一種低通濾波，與高斯濾波不同，維納濾波不需要確定濾波半徑，濾波系數(shù)僅與時變重力場模型的噪聲特性有關(guān)；(2)維納濾波能夠有效抑制條狀誤差，全球陸地水儲量反演結(jié)果與半徑400km的高斯濾波反演結(jié)果相當(dāng)；(3)與高斯濾波相比，維納濾波系數(shù)在低階項權(quán)重更大、高階項權(quán)重更小，因此能夠更好地反映時變重力信號。

[1]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]. J.Geophys.Res, 1998, 103(B12):30205-30229.

[2]李瓊，羅志才，鐘波．利用GRACE時變重力場探測2010年中國西南干旱陸地水儲量變化[J],地球物理學(xué)報， 2013,56(6):1843-1849.

[3]王武星，石耀霖，顧國華等．GRACE衛(wèi)星觀測到的與汶川Ms8.0地震有關(guān)的重力變化[J]．地球物理學(xué)報，2010, 53(8):1767-1776.

[4]蔣濤，李建成，王正濤等．聯(lián)合Jason-1與GRACE衛(wèi)星數(shù)據(jù)研究全球海平面變化[J].測繪學(xué)報，2010，39(2)：134-149.

[5]朱廣彬，李建成，文漢江等．利用GRACE時變位模型研究南極冰蓋質(zhì)量變化[J].武漢大學(xué)學(xué)報信息科學(xué)版，2009，34(10)：1185-1190.

[6]Zi-zhan Zhang, B.F.Chao, Yang Lu, et al.An effective filtering for GRACE time-variable gravity: Fan filter [J]. Geophysical Research Letters, 2009, 36(17):1397-1413.

[7]Sasgen I, Martince Z, Fleming K.Winener optimal filtering of GRACE data [J]. Stud.Geophys.Geod.2006,50:499-508.

[8]Wahr J, Swenson S, Zlotnicki V, et al.Time-variable gravity from GRACE: First results[J].Geophys.Res.Lett,2004,31(11):293-317.

[9]J L Chen，M Rodell，C R Wilson，et al．Low degree spherical harmonic influences on Gravity Recovery and Climate Experiment(GRACE) water storage estimates [J].Geophysical Research Letters, 2005,32(14):57-76.

Wiener Filtering in Use of GRACE Time-variable Gravity Field

Guo Feixiao1, 2,Miao Yuewang3,Xiao Yun1, 2

1. Xi’an Research Institute of Surveying and Mapping, Xi’an 710054, China 2. State Key Laboratory of Geo-information Engineering, Xi’an 710054, China 3.Technical Division of Surveying and Mapping，Xi’an 710054, China

The error of the GRACE time-variable gravity model increases with the increasing of degree. This problem must be solved by spatial smoothing. Global water storage variation is recovered by GRACE time-variable gravity data. The influence of Wiener filtering on the inversion result is analyzed and compared with Gauss filtering of different radius. The results show that Wiener filtering is a low-pass filtering, and weight of Wiener filtering only relys on the noise character of GRACE time-variable gravity models. Stripe error in global water storage variation recovery is removed effectively by Wiener filtering. The inversion result is consistent with radius of 400km Gauss filtering.

GRACE satellites, time-variable gravity field, Wiener Filtering, water storage variation, satellite gravimetry

2015-01-27。

國家自然科學(xué)基金資助項目(41374083)，國家973計劃資助項目(2013CB733303)，大地測量與地球動力學(xué)國家重點實驗室開放基金資助項目(SKLGED2013-3-3-E)。

郭飛霄(1988—)，男，碩士研究生，主要從事重力測量數(shù)據(jù)處理方面的研究。

P228

A