崔立魯，張 誠，何明睿

(成都大學 建筑與土木工程學院，四川 成都 610106)

0 引 言

隨著全球經(jīng)濟快速發(fā)展和人口的快速增加，導致全球氣候變暖，極端天氣頻發(fā)，水資源在時空分布上呈現(xiàn)出嚴重的不均勻現(xiàn)象，極大地影響了人類生存和發(fā)展[1].隨著2002年4月重力恢復與氣候實驗(Gravity Recovery and Climate Explorer,GRACE)衛(wèi)星正式提供服務，一種全新的對地觀測手段出現(xiàn)在各國科學家面前.相對于傳統(tǒng)監(jiān)測技術手段，該手段可以提供低成本、持續(xù)性、大面積的觀測數(shù)據(jù)，且這些觀測數(shù)據(jù)能夠反映出陸地水儲量變化(Terrestrial Water Storage Change,TWSC)總的情況[2].因此，科學家們利用GRACE對衛(wèi)星觀測數(shù)據(jù)進行了大量的水文應用研究[3-5].

由于衛(wèi)星軌道誤差、衛(wèi)星載荷誤差、儀器測量誤差以及重力場模型本身缺陷等，由GRACE時變重力場模型所得到的TWSC結果中出現(xiàn)了顯著的南北條帶誤差，所以必須對GRACE時變重力場模型進行濾波處理以削弱上述誤差.常用的濾波算法有高斯濾波[6]、Fan濾波[7]、去相關濾波[8-9]、Han濾波[10]等，不同濾波算法的處理效果不盡相同，并直接影響到反演結果.稀慧等[11]采用多種濾波方法推算全球平均海水質量變化，與去相關滑動濾波的比較結果顯示，由不同濾波算法處理得到全球平均海水質量變化結果具有顯著的區(qū)別；張青全等[12]比較了4種不同濾波算法在GRACE反演西南地區(qū)巖溶區(qū)陸地水儲量變化的影響，指出不同濾波反演結果在空間分布上差異較大.但是上述比較是以其他濾波方法結果和降水數(shù)據(jù)為基準，并不能反映出真實的陸地水儲量變化.

本研究采用2002年4月至2017年6月連續(xù)15年的GRACE時變重力場模型數(shù)據(jù)反演全球陸地水儲量變化，并利用高斯濾波、Fan濾波、去相關濾波和Han濾波分別對相關誤差進行處理.為了對比4種濾波算法的處理結果，本研究從時空分布，長期趨勢變化和季節(jié)性變化3個方面分別進行闡述，并將4種濾波算法結果與NASA提供的結果進行比較.

1 數(shù)據(jù)與算法

1.1 GRACE數(shù)據(jù)處理

采用由美國德克薩斯大學空間研究中心提供的2002年4月至2017年6月的GRACE RL06時變重力場球諧系數(shù)，其截斷階數(shù)為60.在進行反演前需要對模型的球諧系數(shù)進行一系列的預處理，具體步驟如下：1)采用衛(wèi)星激光測距(Satellite Laser Ranging,SLR)獲取的高精度C20項數(shù)據(jù)對重力場模型相應項的系數(shù)進行替換[13]；2)利用文獻[14]的成果對模型一階項進行地心變化改正；3)采用濾波算法對相關誤差進行處理.

利用GRACE時變重力場反演TWSC的計算公式如下，用等效水高(Equivalent Water High,EWH)表示TWSC[7]，有，

ΔSlmsin(mλ))

(1)

1.2 濾波算法原理

1.2.1 高斯濾波

該濾波是將每個點的密度變化用所有點密度變化的加權平均值替代，以達到平滑效果.當式(1)中Wlm=Wl時，即空間平滑函數(shù)只與階數(shù)相關，可根據(jù)遞歸關系計算高斯濾波系數(shù)Wi(i=1,…,l)，具體如下：

(2)

式中，b=ln2/[1-cos(r/a)]；r為濾波半徑/km;e為自然常數(shù).

1.2.2 Fan濾波

該算法本質上是各向異性濾波，即對球諧系數(shù)的階數(shù)和次數(shù)均采用與高斯濾波相同的處理方式，即式(1)中Wlm=Wl·Wm，具體表達式如下:

ΔSlmsin(mλ))

(3)

式中，Wj(j=1,…,m)的求解方法與式(2)完全一致.

1.2.3 各向異性高斯濾波

考慮到GRACE時變重力場模型誤差是各向異性[10]，因此，高斯濾波僅對階相關誤差進行處理的方法存在著缺陷，各向異性高斯濾波系數(shù)的具體表達式如下：

(4)

式中，r0和r1為濾波半徑,/km，且r0

1.2.4 去相關濾波

去相關濾波的基本原理是保持重力場模型前l(fā)×l階的位系數(shù)不變，將n階多項式擬合得到的大于或等于m階次位系數(shù)從原模型中扣除，以消除誤差[15-16].

1.3 線性擬合模型

為了詳細分析TWSC在時域中的變化規(guī)律，一般采用線性自回歸方程從TWSC時間序列中提取長期趨勢項、長期趨勢加速度項、周年項和半周年項，其具體表達式如下[17]：

ΔEWH(t)=a+bt+ccos(2πt)+dsin(2πt)+ecos(4πt)+fsin(4πt)+ε

(5)

式中，a為常數(shù)，b為長期趨勢，c和d為周年，e和f為半周年，ε為殘差，t為時間.采用最小二乘配置法求解上述方程.同時，TWSC時間序列中的周年項和半周年項的振幅(Aann和Asemi-ann)、相位(Φann和Φsemi-ann)計算公式如下[18]：

(6)

2 實驗結果與分析

本研究計算了2014年10月全球陸地水儲量變化，結果如圖1所示.由圖1(a)可明顯看出陸地水儲量變化結果呈現(xiàn)出顯著的南北分布條帶誤差，由于誤差的干擾造成了正常信號提取的困難，因此對條帶誤差的處理是必要的.分別采用350 km 高斯濾波，300 km Fan濾波，各向異性高斯濾波和P3M6多項式濾波對條帶誤差進行了處理，結果如圖1(b)～(e)所示.對比可知，F(xiàn)an濾波和各向異性高斯濾波的處理效果最為顯著，幾乎看不到條帶誤差的痕跡，但是在削弱條帶誤差影響的同時，真實信號也受到了影響.比較圖(c)和(d)可知，在相同誤差處理效果的前提下，各向異性高斯濾波比Fan濾波能更多地保留真實信號.結合圖(b)和(e)可知，P3M6多項式濾波處理中低緯度地區(qū)的條帶誤差處理效果較好，而高斯濾波則對高緯度地區(qū)的誤差處理效果較好.綜上所述，各向異性高斯濾波既能很好地處理條帶誤差的影響，又可以最大限度的保留真實信號.

圖1 全球陸地水儲量變化濾波結果(2014年10月)

為了進一步比較上述4種濾波算法的處理效果，本研究計算了2002年4月至2017年6月的全球陸地水儲量變化的時間序列如圖2(a)～(d)所示.由于4種算法的結果非常接近，由圖可知4組時間序列的變化趨勢基本一致，同時在數(shù)量級上，除了Fan濾波結果以外，其他濾波也完全相同.比較圖2(b)和圖2(a)、(c)、(d)發(fā)現(xiàn)，F(xiàn)an濾波結果要比其他3種濾波算法結果小一個數(shù)據(jù)量級，這再次印證了Fan濾波除了削弱條帶誤差之外，對真實信號的削弱效果也較其他3種濾波算法更為顯著.本研究計算了4組時間序列的長期趨勢變化、周年振幅、周年相位、半周年振幅和半周年相位，結果如表1所示.

圖2 全球陸地水儲量變化時間序列(2002年4月至2017年6月)

由表1可知，在長期趨勢變化、周年振幅、周年相位、半周年振幅和半周年相位上，350 km 高斯濾波、300 km Fan濾波和P3M6多項式濾波3者結果較為接近.而各向異性高斯濾波則在長期趨勢變化、周年振幅和半周年振幅三方面存在著一定的差別，但是考慮到3種指標的數(shù)量級都很小，這種差別是可以忽略不計的.在周年相位和半周年相位兩種指標方面，350 km 高斯濾波、300 km Fan濾波和各向異性高斯濾波的結果一致，但是300 km Fan濾波與上述3種濾波算法的結果則存在著一定的差異，這與圖2所得的結果相互印證.

表1 全球陸地水儲量變化長期趨勢變化、周年項和半周年項

3 結 論

針對GRACE時變重力場模型反演陸地水儲量變化中出現(xiàn)的條帶誤差問題，本研究采用4種常見的濾波算法(高斯濾波、Fan濾波、各向異性高斯濾波和P3M6多項式濾波)分別對條帶誤差進行處理.為了比較4種濾波算法的誤差處理效果，分別比較了4種濾波算法計算得到的全球陸地水儲量變化及其時間序列，以及長期變化趨勢、周年變化項和半周年變化項.結果表明，各向異性高斯濾波和Fan濾波在進行誤差處理方面要優(yōu)于其他兩種算法，在保留真實信號方面，各向異性高斯濾波要優(yōu)于Fan濾波;在長期變化趨勢、周年變化和半周年變化方面，除了Fan濾波以外，其他3種濾波算法的結果基本上完全一致.綜上所述，各向異性高斯濾波更適合用于對GRACE模型條帶誤差的處理.