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

郭飛霄，苗岳旺，肖 云

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

?

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

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

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

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

GRACE衛星，時變重力場，維納濾波，水儲量變化，衛星重力

1 引 言

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

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

2 維納濾波基本原理

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

(1)

(2)

(3)

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

(4)

(5)

維納濾波滿足以下條件：

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

(6)

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

(7)

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

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

(8)

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

(9)

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

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

(10)

(11)

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

(12)

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

3 利用GRACE時變重力場數據確定維納濾波系數

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

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

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

(13)

(14)

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

(15)

(16)

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

(17)

4 實驗分析

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

(18)

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

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

(19)

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

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

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

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

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

圖4 不同半徑高斯濾波和維納濾波系數

5 結 論

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

[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],地球物理學報， 2013,56(6):1843-1849.

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

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

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

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

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

P228

A