GRACE時變重力場濾波方法研究

姜永濤，張永志，王 帥，姚曉偉

(1.長安大學，陜西 西安 710000；2.河南工程學院，河南 鄭州 451191)

GRACE時變重力場濾波方法研究

姜永濤1，2，張永志1，王 帥1，姚曉偉1

(1.長安大學，陜西 西安 710000；2.河南工程學院，河南 鄭州 451191)

用3種不同的濾波方法獲得了2007年相對于2005年的衛(wèi)星重力變化圖像，并與同期地面重力測量的結果進行比較，對比分析GRACE月重力場濾波方法的優(yōu)缺點。結果表明，去相關平滑濾波算法優(yōu)于高斯濾波和直接截斷法，且去相關平滑濾波DDK5處理得到的衛(wèi)星重力動態(tài)變化圖像與地面觀測結果符合最好，表明GRACE衛(wèi)星時變重力場可以用來分析大區(qū)域重力動態(tài)變化。

GRACE；濾波方法；去相關濾波；高斯濾波；重力場變化

一、引 言

GRACE(gravity recovery and climate experiment)雙星系統(tǒng)利用衛(wèi)星軌道(SST-hl，GPS定軌數(shù)據(jù))和星間距離測量數(shù)據(jù)(SST-ll，KBR星間測距數(shù)據(jù))來解算時變重力場模型，利用這些時變重力場模型，國內(nèi)外學者第一次觀測到了世界范圍的海水重新分布[1]、區(qū)域水儲量變化[2]和大震的同震及震后重力變化[3]。由于軌道傾角約為89°，使得月重力場模型位系數(shù)之間具有相關性，且誤差隨階次的增高而增大，在空間域表現(xiàn)為明顯的南—北向條帶狀波紋。因此利用GRACE時變重力場數(shù)據(jù)時，首先要進行重力場位系數(shù)去相關和平滑處理。

將重力場模型截斷到低階，可以消除中高階位系數(shù)誤差對結果的影響，但同時也丟掉了中高階位系數(shù)的重力場信息，而且截斷法獲得的重力場空間分辨率較低，很難用于研究較小尺度(100～300 km左右)的地球動力學現(xiàn)象[3]；高斯濾波法[4-5]是早期研究GRACE時變重力場常用的方法，它通過對模型高階次位系數(shù)的降權處理來降低誤差對重力場變化結果的影響；現(xiàn)在青睞于采用去相關并結合平滑后處理技術處理GRACE時變重力場數(shù)據(jù)，去相關處理是確定并消除模型位系數(shù)之間的相關性誤差，經(jīng)驗分析[6-7]和先驗模型[8-9]方法的去相關核都不是軸對稱的形式，而是呈現(xiàn)出南—北向較短的旁瓣狀的形狀，且其大小與點坐標有關。文獻[8]利用GRACE時變重力場反演過程中的信號和誤差的協(xié)方差陣，基于懲罰加權的方法得到時變重力場位系數(shù)的去相關系數(shù)矩陣，但該矩陣包含的系數(shù)太多，不利于計算。GFZ[10]借鑒文獻[6]簡化了去相關系數(shù)矩陣，最終得到了去相關濾波方法(DDK1—DDK5)，相應的GRACE月重力場模型處理結果可在ICGEM網(wǎng)站上獲得。

本文用3種不同的濾波方法獲得了2007年相對于2005年的衛(wèi)星重力變化圖像，并與同期地面重力測量的結果[11]進行了比較，對比分析了GRACE月重力場濾波方法的優(yōu)缺點。

二、衛(wèi)星年時變重力場變化

由GRACE重力場模型得到的橢球面上的重力變化可表達為

式中，θ、λ分別為地心余緯和經(jīng)度；r≈a(1－fcos2θ)，為橢球面上一點的矢徑，a、f分別為參考橢球的長半軸和扁率；GM為地心引力常數(shù)；R為地球平均半徑；n、m分別為模型位系數(shù)的階和次；nm(cos θ)為完全正則化勒讓德函數(shù)[12]；nm、nm為所選重力場模型相應階次的位系數(shù)差值。

本文將2007年1—3月(冬季，可忽略區(qū)域水儲量變化對重力的影響)的GRACE時變平均重力場(簡稱2007S)和2005年1—3月的GRACE時變平均重力場(簡稱2005S)進行差分計算，得到由GRACE衛(wèi)星重力獲得的2007年相對于2005年的重力變化圖像。鑒于GRACE雙星系統(tǒng)對位系數(shù)C20項不敏感，因此計算中視C20項為常數(shù)，并通過差分消除其影響。式(1)中取，可以得到由GRACE衛(wèi)星重力場模型得到的2007年相對于2005年的差分重力場。

圖1以點(34.5，105)處的重力變化為例計算了差分重力場的各階異常大小及其標準差，可以看出差分重力場階標準差隨著階數(shù)n的增高而增大，到90階時，標準差可達10 μGal(1 μGal＝1×10－8m/s2)。圖2為未經(jīng)濾波處理的2007年相對于2005年的衛(wèi)星重力場變化圖像，圖中可以看到明顯的南—北向的條帶狀誤差，且條紋波動可達上百μGal，因此在利用衛(wèi)星重力場數(shù)據(jù)前必須進行濾波處理。

圖1 重力階變化及其標準差(單位：10－8m/s2)

圖2 未經(jīng)濾波處理的衛(wèi)星重力變化圖像(2007S—2005S，單位：10－8m/s2)

三、時變重力場濾波方法

GRACE時變重力場的常用濾波方法有高斯濾波和去相關平滑濾波。

1.高斯濾波

高斯濾波通過對重力場位系數(shù)(頻率域)的加權運算來實現(xiàn)空間域的平滑處理，對式(1)的時變重力模型，其高斯濾波[5]可寫為

式中，Wnm為高斯權函數(shù)。各向同性高斯平滑權函數(shù)Wnm只與模型階數(shù)有關，即Wnm＝Wn。最初的Wn是根據(jù)勒讓德多項式基于分解公式給出的[13]，實用的權函數(shù)公式是將Wn作為平滑半徑rfilter的函數(shù)，其遞推公式為[4，14]

圖3為各向同性高斯平滑權函數(shù)隨濾波半徑的變化圖，可以看出各向同性高斯平滑是通過壓抑高階面諧函數(shù)的貢獻(同時壓制高階項誤差)來達到濾波的目的。由圖3分析可知，400 km濾波半徑的高斯平滑60階以后的位系數(shù)被嚴重壓制，且30～60階位系數(shù)被降權處理，因此可以得到較為平滑的重力場變化結果，如圖4所示。相對于未經(jīng)濾波的衛(wèi)星重力場變化圖像(如圖2所示)，圖4中已不含南—北向條帶誤差，且量值上大為減少，其原因是400 km的各向同性高斯濾波降低了高于50階位系數(shù)誤差(如圖1所示)對重力場變化的影響。

圖3 不同平滑半徑的各向同性高斯權函數(shù)變化

圖4 各向同性高斯平滑(400 km)的重力變化圖像(2007S—2005S，單位：10－8m/s2)

由于時變重力場位系數(shù)誤差與其階數(shù)n和次數(shù)m均有關，文獻[5]提出各向異性高斯濾波，其遞推式為

式中，r0、r1分別對應于次數(shù)為0和m1的平均濾波半徑；m1為濾波截斷次數(shù)。Wnm的確定與r0、r1、m1的選取有關，其中，r0決定了南—北方向上的平滑半徑；r1、m1決定了在東—西方向上的平滑半徑。因此通過選取合適的各向異性高斯濾波參數(shù)，可以提高時變重力場在南—北向的分辨率。

本文試算了不同r0、r1、m1參數(shù)得到的時變重力場圖像，并從條紋消除和變化細節(jié)保留兩方面分析得出r1＝400 km，r0＝200 km，m1＝60的結果最佳，如圖5所示。圖6給出了上述r1、r0、m1取值下的各向異性高斯濾波權系數(shù)，可以看出各向異性高斯濾波同階的權函數(shù)不再相同，且考慮了模型高階低次(m＞m1)的位系數(shù)對時變重力異常的貢獻。由于r0＝200 km，圖5相比于各向同性高斯平滑(如圖4所示)，在南—北方向上更凸顯了細節(jié)，且量值上略大于400 km平滑半徑的各向同性高斯濾波結果。

2.去相關平滑濾波

去相關平滑濾波(DDK1—DDK5)是GFZ根據(jù)Kusche[8]和Swenson等人[6]的研究得出的一個卷積矩陣較為簡單的去相關平滑方法。不同于高斯濾波，去相關平滑濾波(DDK1—DDK5)考慮了由GRACE軌道數(shù)據(jù)和星間距離測量數(shù)據(jù)反演月重力場模型過程中產(chǎn)生的信號協(xié)方差和誤差協(xié)方差矩陣，通過懲罰加權法構建去相關系數(shù)矩陣，又針對去相關濾波矩陣系數(shù)的特點[6]，發(fā)展出的一個“次卷積”去相關濾波核。位系數(shù)的濾波結果可表示為

式中，W( n，n′，m，a)為去相關濾波核；a決定平滑程度；n′∈parity(n)表示n′與n的奇偶性相同；n′m、n′m為未經(jīng)濾波的時變重力場位系數(shù)。由式(5)可以看出，濾波后的某一階次的位系數(shù)可由與該位系數(shù)同次不同階的位系數(shù)的加權求和，即“次卷積運算”獲得。

圖5 各向異性高斯平滑(400 km，200 km，60)的重力變化圖像(2007S—2005S，單位：10－8m/s2)

圖6 各向異性高斯平滑(400 km，200 km，60)的各階次權系數(shù)

GFZ去相關核W( n，n′，m，a)是通過矩陣的方式提供的，此外經(jīng)過DDK濾波處理的時變重力場數(shù)據(jù)可從ICGEM網(wǎng)站下載。圖7為2007年相對于2005年的衛(wèi)星重力場變化DDK3濾波結果，可以看出，去相關濾波不僅消除了條帶誤差，且量值上較高斯濾波結果大，說明DDK3濾波較有效地保留了時變重力場的高階次信息。

圖7 去相關濾波DDK3處理的重力變化圖像(2007S—2005S，單位：10－8m/s2)

四、時變重力場處理方案

文獻[11]給出了由中國大陸地殼運動網(wǎng)絡、中國數(shù)字地震觀測網(wǎng)絡的重力網(wǎng)流動觀測成果得出的2007年相對于2005年中國大陸動態(tài)重力場變化圖像。本節(jié)參照同期地面重力測量的結果(文獻[11])來對比分析GRACE月重力場濾波方法的優(yōu)缺點。

對于各向同性高斯濾波，通過試算，濾波半徑400 km的各向同性高斯平滑得到的結果最優(yōu)，如圖4所示。由圖3中不同濾波半徑的高斯權系數(shù)可知，400 km的高斯平滑已經(jīng)完全壓制了60階以后的位系數(shù)，且使得30～60階位系數(shù)的貢獻值低于30%，因此雖然在量級上圖4小于直接截斷結果(如圖8所示)，但在時變重力場細節(jié)上，400 km各向同性高斯平滑優(yōu)于直接截斷法，且與地面重力測量得到的動態(tài)變化圖像(文獻[11])高值區(qū)和低值區(qū)的對應性上更好些。

圖8 截斷36階得到的重力變化圖像(2007S—2005S，單位：10－8m/s2)

對于各向異性高斯濾波，試算發(fā)現(xiàn)r1＝400 km，r0＝200 km，m1＝60時結果最優(yōu)，如圖6所示。綜合比較圖6、圖4和圖8的重力場變化高、低值區(qū)可以發(fā)現(xiàn)，相對同性高斯濾波，各向異性高斯濾波確實提高了南—北方向上的重力場變化分辨率，但實際意義兩者沒有大的差異。

利用ICGEM提供的DDK1—DDK5月重力場模型數(shù)據(jù)計算了2007S相對于2005S的年重力場變化，結果表明，DDK3—DDK5的結果與文獻[11]的空間對應較好，原因是DDK1—DDK2處理的重力場等效平滑半徑較大，平滑了重力變化細節(jié)。通過與文獻[11]比較可知，DDK5處理衛(wèi)星重力場變化結果(如圖9所示)最優(yōu)。

圖9 去相關濾波DDK5處理的重力變化圖像(2007S—2005S，單位：10－8m/s2)

五、結 論

本文利用GRACE衛(wèi)星時變重力數(shù)據(jù)，采用不同的濾波方法計算了2007年相對于2005年的衛(wèi)星重力動態(tài)變化圖像，通過與地面重力測量結果的對比分析，得出如下結論：

1)去相關平滑濾波算法優(yōu)于高斯濾波和直接截斷法。直接截斷法刪除了中高階項重力場信息，降低了空間分辨率；高斯濾波法隨階數(shù)增高而降低的權函數(shù)在一定程度上減少了最終結果的量值；去相關濾波由于其濾波核的獲取具有一定的物理意義，因此得到了相對較好的處理結果。

2)各種濾波方法得到的衛(wèi)星重力場變化的高、低值區(qū)空間分布與地面重力觀測結果具有較好的一致性，但量值上均相差約一個量級。高、低值區(qū)客觀上的良好對應關系表明GRACE時變重力數(shù)據(jù)可以較好地應用在地學研究中。

本文得到的結論對水文學和地學相關研究中選定時變重力場濾波方法有一定的參考作用。

[1] CHAMBERS D，WAHR J，NEREM R S.Preliminary Observations of Global Ocean Mass Variations with GRACE[J].Geophysical Research Letters，2004，33 (13).DOI：10.1029/2004GL020461.

[2] SCHMIDT R，SCHWINTZER P，F(xiàn)LECHTER F，et al. GRACE Observations of Changes in Continental Water Storage[J].Global and Planetary Change，2006，50(1-2)：112-126.

[3] HAN S C，SIMONS F J.Spatiospectral Localization of Global Geopotential Fields from the Gravity Recovery and Climate Experiment(GRACE)Reveals the Coseismic Gravity Change Owing to the 2004 Sumatra-Andaman Earthquake[J].Journal of Geophysical Research：Solid Earth，2008，113(B1)：B01405.DOI：10.1029/ 2007JB004927.

[4] WAHR J，MOLENAAR F，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：Solid Earth，1998， 103(B12)：30205-30229.

[5] HAN Shin-Chan，SHUM C K，JEKELI C，et al.Non-isotropic Filtering of GRACE Temporal Gravity for Geophysical Signal Enhancement[J].Geophysical Journal International，2005，163(1)：18-25.

[6] SWENSON S，WAHR J.Post-processing Removal of Correlated Errors in GRACE Data[J].Geophysical Research Letters，2006，33(8)：L08402.DOI：10.1029/ 2005GL025285.

[7] WOUTERS B，SCHRAMA E J O.Improved Accuracy of GRACE Gravity Solutions through Empirical Orthogonal Function Filtering of Spherical Harmonics[J].Geophysical Research Letters，2007，34(23)：L23711.DOI：10. 1029/2007GL032098.

[8] KUSCHE J.Approximate Decorrelation and Non-isotropic Smoothing of Time-variable GRACE-type Gravity Field Models[J].Journal of Geodesy，2007，81(11)：733-749.DOI：10.1007/s00190-007-01433.

[9] KLEES R，REVTOVA E A，GUNTER B，et al.The Design of an Optimal Filter for Monthly GRACE Gravity Models[J].Geophysical Journal International，2008，175 (2)：417-432.DOI：10.1011/j.1365-246X.2008.03922. x.

[10] KUSCHE J，SCHMIDT R，PETROVIC S，et al.Decorrelated GRACE Time-variable Gravity Solutions by GFZ，and Their Validation Using a Hydrological Model[J]. Journal of Geodesy，2009，83(10)：903-913.DOI：10. 1007/s00190-009-0308-3.

[11] 李輝，申重陽，孫少安，等.中國大陸近期重力場動態(tài)變化圖像[J].大地測量與地球動力學，2009，29(3)：1-10.

[12] HOLMES S A，F(xiàn)EATHERSTONE W E.A Unified Approach to the Clenshaw Summation and the Recursive Computation of Very High Degree and Order Normalised Associated Legendre Functions[J].Journal of Geodesy，2002，76(5)：279-299.

[13] HOBSON E W.The Theory of Spherical and Ellipsoidal Harmonics[M].New York：Chelsea Publishing Company，2012.

[14] JEKELI C.Alternative Methods to Smooth the Earth’s Gravity Field[R].Ohio：Ohio State University，1981.

Study on Filtering Methods of GRACE Temporal Gravity Variation

JIANG Yongtao，ZHANG Yongzhi，WANG Shuai，YAO Xiaowei

P223

B

0494-0911(2014)11-0001-05

2013-09-12；

2014-08-29

國家自然科學基金(41374028；41274083；41304013)；國土資源大調(diào)查項目(1212010914015)

姜永濤(1985—)，男，山東菏澤人，博士生，主要研究方向為衛(wèi)星重力。

姜永濤，張永志，王帥，等.GRACE時變重力場濾波方法研究[J].測繪通報，2014(11)：1-5.

10.13474/j.cnki.11-2246.2014. 0350