利用重力衛星監測全球陸地水儲量變化時不同濾波算法的影響

崔立魯，張 誠，何明睿

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

0 引 言

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

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

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

1 數據與算法

1.1 GRACE數據處理

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

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

ΔSlmsin(mλ))

(1)

1.2 濾波算法原理

1.2.1 高斯濾波

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

(2)

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

1.2.2 Fan濾波

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

ΔSlmsin(mλ))

(3)

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

1.2.3 各向異性高斯濾波

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

(4)

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

1.2.4 去相關濾波

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

1.3 線性擬合模型

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

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

(5)

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

(6)

2 實驗結果與分析

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

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

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

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

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

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

3 結 論

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