使用GRACE時(shí)變重力場(chǎng)中的維納濾波

郭飛霄，苗岳旺，肖 云

1.西安測(cè)繪研究所，陜西 西安，710054；2.地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安，710054；3.測(cè)繪信息技術(shù)總站，陜西 西安，710054

?

使用GRACE時(shí)變重力場(chǎng)中的維納濾波

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

1.西安測(cè)繪研究所，陜西 西安，710054；2.地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安，710054；3.測(cè)繪信息技術(shù)總站，陜西 西安，710054

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

GRACE衛(wèi)星，時(shí)變重力場(chǎng)，維納濾波，水儲(chǔ)量變化，衛(wèi)星重力

1 引 言

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

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

2 維納濾波基本原理

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

(1)

(2)

(3)

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

(4)

(5)

維納濾波滿足以下條件：

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

(6)

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

(7)

(2)測(cè)量信號(hào)x(Ω)由真實(shí)信號(hào)s(Ω)和噪聲n(Ω)組成：

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

(8)

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

(9)

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

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

(10)

(11)

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

(12)

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

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

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

圖1 GRACE時(shí)變重力場(chǎng)模型平均階方差

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

(13)

(14)

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

(15)

(16)

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

(17)

4 實(shí)驗(yàn)分析

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

(18)

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

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

(19)

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

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

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

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

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

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

5 結(jié) 論

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

[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時(shí)變重力場(chǎng)探測(cè)2010年中國(guó)西南干旱陸地水儲(chǔ)量變化[J],地球物理學(xué)報(bào)， 2013,56(6):1843-1849.

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

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

[5]朱廣彬，李建成，文漢江等．利用GRACE時(shí)變位模型研究南極冰蓋質(zhì)量變化[J].武漢大學(xué)學(xué)報(bào)信息科學(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。

國(guó)家自然科學(xué)基金資助項(xiàng)目(41374083)，國(guó)家973計(jì)劃資助項(xiàng)目(2013CB733303)，大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目(SKLGED2013-3-3-E)。

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

P228

A