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

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

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

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

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

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

GRACE；濾波方法；去相關(guān)濾波；高斯濾波；重力場變化

一、引 言

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

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

本文用3種不同的濾波方法獲得了2007年相對于2005年的衛(wèi)星重力變化圖像，并與同期地面重力測量的結(jié)果[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為所選重力場模型相應(yīng)階次的位系數(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時變重力場的常用濾波方法有高斯濾波和去相關(guān)平滑濾波。

1.高斯濾波

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

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

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

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

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

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

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

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

2.去相關(guān)平滑濾波

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

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

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

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

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

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

四、時變重力場處理方案

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

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

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

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

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

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

五、結(jié) 論

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

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

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

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

[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].大地測量與地球動力學(xué)，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

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

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

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

10.13474/j.cnki.11-2246.2014. 0350