崔立魯，張 誠(chéng)，何明睿

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

0 引 言

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

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

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

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

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

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

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

ΔSlmsin(mλ))

(1)

1.2 濾波算法原理

1.2.1 高斯濾波

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

(2)

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

1.2.2 Fan濾波

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

ΔSlmsin(mλ))

(3)

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

1.2.3 各向異性高斯濾波

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

(4)

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

1.2.4 去相關(guān)濾波

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

1.3 線(xiàn)性擬合模型

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

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

(5)

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

(6)

2 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

表1 全球陸地水儲(chǔ)量變化長(zhǎng)期趨勢(shì)變化、周年項(xiàng)和半周年項(xiàng)

3 結(jié) 論

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