混合半徑高斯濾波算法在去除GRACE條帶誤差中的應(yīng)用

付 林 趙東明 付林威

1 信息工程大學(xué)地理空間信息學(xué)院,鄭州市科學(xué)大道62號(hào),450001

重力場(chǎng)是地球的基本物理場(chǎng),通過(guò)確定地球重力場(chǎng)及其變化可以反映出地球內(nèi)部質(zhì)量分布以及地表物質(zhì)遷移情況[1]。GRACE(gravity recovery and climate experiment)重力衛(wèi)星為監(jiān)測(cè)大尺度地球重力場(chǎng)變化、恢復(fù)地球重力場(chǎng)提供支持,它的出現(xiàn)解決了人們無(wú)法持續(xù)有效觀測(cè)地球重力場(chǎng)這一難題,為地球科學(xué)研究作出重要貢獻(xiàn)[2]。

由于GRACE采用南、北雙星跟蹤飛行模式,導(dǎo)致其觀測(cè)信息也是呈南北條形分布, GRACE獲得的地球時(shí)變重力場(chǎng)信息會(huì)出現(xiàn)嚴(yán)重的南北條帶誤差,這在很大程度上掩蓋了真實(shí)重力信號(hào)[3]。目前,國(guó)內(nèi)外學(xué)者處理?xiàng)l帶誤差的主要手段分為兩大類:一是通過(guò)降低含有噪聲的高階項(xiàng)權(quán)重來(lái)降低噪聲,比如高斯濾波等,該類方法雖然可以較好地濾除噪聲,但由于是通過(guò)降低高階權(quán)重實(shí)現(xiàn)的,不可避免地會(huì)損失掉一些真實(shí)高頻信號(hào),而且隨著濾波半徑的增大,信號(hào)損失會(huì)更嚴(yán)重;二是通過(guò)濾除重力場(chǎng)球諧系數(shù)奇偶項(xiàng)之間的相關(guān)性而達(dá)到濾除噪聲的目的[4],如去相關(guān)濾波[5]等,此類方法對(duì)濾波造成的信號(hào)損失較小,但對(duì)于低緯度地區(qū)的濾波效果較差,仍會(huì)存在條帶誤差。因此,為實(shí)現(xiàn)最大程度地保留真實(shí)信號(hào)并且擁有較好的濾波效果,人們通常采用去相關(guān)與高斯濾波相結(jié)合的組合濾波方法。同時(shí),在GRACE運(yùn)行期間存在部分月份數(shù)據(jù)缺失,對(duì)此常用三次樣條插值法對(duì)短期缺失數(shù)據(jù)進(jìn)行填補(bǔ)。而對(duì)于GRACE與GRACE-FO兩代衛(wèi)星之間11個(gè)月的數(shù)據(jù)空白,莫紹興等[6]提出利用貝葉斯卷積神經(jīng)網(wǎng)絡(luò)提取時(shí)變重力場(chǎng)信號(hào)有效特征,從而對(duì)缺失數(shù)據(jù)進(jìn)行預(yù)測(cè);也有專家提出利用奇異譜分析法[7]恢復(fù)兩代衛(wèi)星間的缺失值,均取得較好效果。

1 濾波處理方法

利用GRACE時(shí)變重力場(chǎng)反演地球表面質(zhì)量變化,在大尺度范圍內(nèi),我們通常近似認(rèn)為地表質(zhì)量變化所引起的時(shí)變重力場(chǎng)變化實(shí)際就是地面水質(zhì)量的變化。因而通過(guò)地表質(zhì)量變化模型再除以水的密度便可以得到陸地等效水高,進(jìn)而計(jì)算陸地水儲(chǔ)量變化[8]:

(1)

1.1 Gaussian濾波方法

Wahr等[8]提出一種基于平滑函數(shù)的高斯濾波方法,并推導(dǎo)了一組高斯平滑濾波系數(shù),應(yīng)用于GRACE數(shù)據(jù)處理中。高斯濾波削弱條帶誤差的主要原理就是通過(guò)降低球諧系數(shù)高階項(xiàng)的權(quán)重,將高階項(xiàng)中存在的條帶誤差信號(hào)濾除掉。其高斯平滑核函數(shù)定義為:

(2)

式中,α為球面兩點(diǎn)夾角,r為濾波半徑。給出高斯濾波的遞推系數(shù):

(3)

Chambers等[9]研究發(fā)現(xiàn),當(dāng)以上高斯濾波算法在球諧系數(shù)高于50階時(shí)存在不確定性,因此提出高斯濾波改進(jìn)算法。改進(jìn)算法的平滑函數(shù)如下:

(4)

1.2 去相關(guān)濾波

Swenson和Wahr[10]研究發(fā)現(xiàn),GRACE條帶誤差存在的一個(gè)重要原因就是其時(shí)變重力場(chǎng)球諧系數(shù)奇偶項(xiàng)之間存在明顯的相關(guān)性,因此國(guó)內(nèi)外專家學(xué)者提出利用多項(xiàng)式擬合來(lái)去除相關(guān)性的去相關(guān)濾波方法,其基本原理就是對(duì)球諧系數(shù)項(xiàng)的奇、偶數(shù)階分別進(jìn)行滑動(dòng)固定窗口的多項(xiàng)式擬合,擬合值認(rèn)為是存在的相關(guān)誤差,將原始球諧系數(shù)減去擬合值即為去相關(guān)濾波后的球諧系數(shù)。

Chen等[11]根據(jù)上述基本思想進(jìn)行改進(jìn),提出PnMl方法。該方法原理主要為:對(duì)系數(shù)的前l(fā)階位系數(shù)不作處理,對(duì)大于等于l階的位系數(shù)用n階多項(xiàng)式進(jìn)行奇偶項(xiàng)擬合并扣除相關(guān)值。

不同品種的茼蒿葉片中葉綠素a、b含量、總?cè)~綠素含量及SPAD值都不相同，對(duì)SPAD值及葉綠素含量間的相互關(guān)系進(jìn)行統(tǒng)計(jì)分析，發(fā)現(xiàn)SPAD值與葉綠素含量均呈顯著正相關(guān)性，與水稻、小麥、草莓、園林樹木、烤煙、無(wú)花果等的研究結(jié)果相似[12-14]。而“小葉茼蒿”的4種數(shù)學(xué)模型的相關(guān)系數(shù)均達(dá)到0.6以上，葉綠素a含量與SPAD值的最優(yōu)函數(shù)模型是指數(shù)方程，葉綠素b與SPAD值的最優(yōu)函數(shù)模型是對(duì)數(shù)函數(shù)，而總?cè)~綠素含量與SPAD值的最優(yōu)函數(shù)模型則是線性方程?！按笕~茼蒿”的4種數(shù)學(xué)模型的相關(guān)系數(shù)普遍較低，由此可推斷“大葉茼蒿”的SPAD值與葉綠素含量相關(guān)性較差。

同時(shí)研究發(fā)現(xiàn),GRACE時(shí)變重力場(chǎng)模型的球諧系數(shù)精度與階次數(shù)密切相關(guān),低階次系數(shù)間的相關(guān)誤差小,但球諧系數(shù)的標(biāo)準(zhǔn)差隨著階數(shù)和次數(shù)的增大而增大。Duan等[12]在Swenson等[10]提出的滑動(dòng)固定窗口多項(xiàng)式擬合去相關(guān)濾波的基礎(chǔ)上,建立滑動(dòng)可變窗口多項(xiàng)式擬合去相關(guān)濾波方法。

1.3 組合濾波方法

隨著高斯濾波半徑的增大,對(duì)條帶誤差的濾除效果越好,但由于高斯平滑濾波是通過(guò)降低高階項(xiàng)的權(quán)重系數(shù)來(lái)壓制噪聲,因此會(huì)不可避免地濾除掉高階項(xiàng)真實(shí)信號(hào)。去相關(guān)濾波也可在一定程度上去除條帶誤差,但在中緯度地區(qū)仍存在噪聲信息。因而,為在濾除條帶噪聲的同時(shí)盡量保留真實(shí)信號(hào),國(guó)內(nèi)外學(xué)者發(fā)現(xiàn)組合濾波方法,即去相關(guān)濾波與高斯濾波相結(jié)合的方法可以最大程度地壓制噪聲并保留真實(shí)信號(hào)[4]。

1.4 信噪比估計(jì)

目前對(duì)各種濾波效果的評(píng)估國(guó)內(nèi)外尚未有統(tǒng)一標(biāo)準(zhǔn),通常采用的是通過(guò)分析濾波前后誤差的RMS(root-mean-square)以及陸地與海洋信噪比(signal-to-noise ratio,SNR)方法。經(jīng)研究,GRACE測(cè)量誤差在陸地與海洋上大致處于同一水平,而陸地質(zhì)量變化信號(hào)強(qiáng)于海洋,因此Chen等[11]根據(jù)此原理構(gòu)建出一套基于陸地海洋信噪比最大為準(zhǔn)則的濾波器,其信噪比公式為:

(5)

式中,MASSland為陸地質(zhì)量變化,MASSocean為海洋質(zhì)量變化,Err為GRACE測(cè)量誤差。

2 實(shí)驗(yàn)分析

2.1 數(shù)據(jù)預(yù)處理

本文采用得克薩斯空間研究中心(Center for Space Research, CSR)的GRACE RL06月重力場(chǎng)模型,模型截?cái)嚯A數(shù)為60,選取2004-01～2010-12數(shù)據(jù)作為研究對(duì)象,該范圍內(nèi)GRACE衛(wèi)星運(yùn)行正常,數(shù)據(jù)質(zhì)量較好,且不存在數(shù)據(jù)空白月份。

由于GRACE衛(wèi)星對(duì)代表地心運(yùn)動(dòng)的一階項(xiàng)不敏感,導(dǎo)致一階項(xiàng)缺失,這里采用Swenson等[10]計(jì)算的一階項(xiàng)進(jìn)行補(bǔ)充。受GRACE衛(wèi)星運(yùn)行軌道因素影響,其解算的C20項(xiàng)存在較大不確定性,利用衛(wèi)星激光測(cè)距(SLR)獲取的C20項(xiàng)對(duì)其進(jìn)行替換。為體現(xiàn)出重力場(chǎng)時(shí)變信息,這里選取所有月份球諧系數(shù)平均值作為背景重力場(chǎng),將每月的球諧系數(shù)減去該背景重力場(chǎng),得到各月時(shí)變重力場(chǎng)球諧系數(shù)。

2.2 高斯濾波信噪比分析

郭飛宵等[3]指出,時(shí)變重力場(chǎng)模型低于20階的項(xiàng)主要反映地幔以下的重力場(chǎng)信號(hào),20～30階項(xiàng)體現(xiàn)信號(hào)與噪聲混合信息,而高于30階項(xiàng)則包含更多的噪聲信息。由表1可見,高斯濾波半徑為0 km(即對(duì)條帶誤差不進(jìn)行處理)時(shí),隨著階數(shù)的不斷增大,信噪比不斷變小,在球諧系數(shù)大于30階時(shí)信噪比減小較快,至60階時(shí)降低為1.252。這也證實(shí)高階項(xiàng)以噪聲為主,對(duì)真實(shí)信號(hào)掩蓋嚴(yán)重,造成信噪比急劇減小。隨著高斯濾波方法的加入,GRACE球諧系數(shù)中的噪聲誤差被減弱,信噪比有較明顯提升,從不同階下的信噪比情況也可以看出,當(dāng)高斯濾波半徑為400 km時(shí),其信噪比在各階均處于最優(yōu)水平,因此在僅用高斯濾波去除條帶誤差時(shí),常選取400 km半徑的高斯濾波作為最優(yōu)濾波半徑。

表1 各階數(shù)下不同濾波半徑的高斯濾波信噪比

2.3 組合濾波信噪比分析

研究發(fā)現(xiàn),GRACE球諧系數(shù)奇偶項(xiàng)之間具有相關(guān)性,因此認(rèn)為GRACE球諧系數(shù)誤差中存在相關(guān)性誤差。為有效去除GRACE存在的誤差,目前常采用組合濾波的方法,即去相關(guān)濾波+高斯濾波方法,可以最大程度地減弱條帶誤差項(xiàng)。這里選取Duan等[12]提出的去相關(guān)濾波與高斯濾波結(jié)合的組合濾波方法,不同階數(shù)以及高斯濾波半徑下的信噪比如表2所示。可以看出,在不進(jìn)行高斯濾波的情況下,經(jīng)去相關(guān)濾波后的時(shí)變重力場(chǎng)信噪比有明顯提升,這也印證了相關(guān)性誤差確實(shí)存在于時(shí)變重力場(chǎng)球諧系數(shù)中。同時(shí),分析去相關(guān)濾波后不同半徑下的高斯濾波信噪比也可以看出,由于去相關(guān)濾波減弱了相關(guān)性誤差,因此在300 km高斯濾波半徑下其信噪比在各階下均處于最優(yōu)水平。

表2 各階數(shù)下不同濾波半徑的組合濾波信噪比

2.4 混合半徑高斯濾波算法

不同半徑下的高斯濾波權(quán)重系數(shù)如圖1所示,可以看出,隨著濾波半徑的增大,高階項(xiàng)的權(quán)重系數(shù)下降得越來(lái)越快。雖然壓制了高階項(xiàng)噪聲,但由于權(quán)重系數(shù)下降過(guò)快,也會(huì)導(dǎo)致信號(hào)損失過(guò)多,尤其在含誤差較少的低階項(xiàng)部分,由于濾波半徑增大,權(quán)重系數(shù)減小較快會(huì)導(dǎo)致真實(shí)低頻信號(hào)損失。

圖1 不同濾波半徑高斯權(quán)重系數(shù)Fig.1 Gaussian weight coefficients with different filtering radius

分析表1可知,當(dāng)球諧系數(shù)階數(shù)較低(如20階)時(shí),里面所包含的噪聲信息較少,不同高斯濾波半徑與原始球諧系數(shù)相比信噪比變化較小,且隨著濾波半徑的增大,權(quán)重系數(shù)下降快,會(huì)造成真實(shí)信號(hào)較快損失。為應(yīng)對(duì)這種情況,本文提出一種混合半徑的高斯濾波方法,即通過(guò)調(diào)整球諧系數(shù)不同階次項(xiàng)的高斯權(quán)重系數(shù)進(jìn)行分段高斯濾波,實(shí)現(xiàn)在盡可能濾除高頻誤差的同時(shí)盡量多地保留低頻信息。

首先確定階數(shù)L、濾波半徑R和分段步數(shù)S,則可以將階數(shù)和高斯濾波半徑按下式劃分:

(6)

式中,r0為原始濾波半徑按步長(zhǎng)劃分的濾波半徑段長(zhǎng),l0為給定球諧階數(shù)按步長(zhǎng)劃分后每段階數(shù)大小,ri、li分別為第i段對(duì)應(yīng)的高斯濾波半徑以及所取的部分階數(shù)。

根據(jù)式(4)給出的不同階數(shù)半徑下的改進(jìn)高斯濾波權(quán)重系數(shù)Wl,可計(jì)算出每段高斯濾波半徑對(duì)應(yīng)的濾波半徑下高斯權(quán)重系數(shù):

(7)

提取不同階數(shù)范圍的高斯權(quán)重系數(shù)并進(jìn)行重組,便可得到拼接后的高斯濾波權(quán)重系數(shù)。由于不同階數(shù)對(duì)應(yīng)濾波系數(shù)不同,為保證高斯濾波的平穩(wěn)性,還要對(duì)組合拼接的濾波系數(shù)進(jìn)行移動(dòng)平均平滑處理,最終得到改進(jìn)的高斯濾波權(quán)重系數(shù):

(8)

根據(jù)表1及表2可得出結(jié)論,只進(jìn)行高斯濾波情況下,400 km高斯濾波的信噪比最優(yōu),進(jìn)行組合濾波時(shí),Duan等[12]與300 km半徑高斯濾波信噪比最優(yōu)。混合半徑高斯濾波方法的目的就是在相同濾波半徑下盡可能保留信號(hào)、提升信噪比。利用式(8)改進(jìn)的高斯濾波方法進(jìn)行信噪比分析(表3),可以看出,當(dāng)僅用一步進(jìn)行分段高斯濾波時(shí),分段高斯濾波即為經(jīng)典高斯濾波;在僅進(jìn)行混合半徑高斯濾波的情況下,400 km濾波半徑分兩步進(jìn)行濾波信噪比最大;在進(jìn)行組合濾波時(shí),分3步進(jìn)行的300 km混合半徑高斯濾波信噪比最大。

表3 混合半徑高斯濾波信噪比

為進(jìn)一步驗(yàn)證該方法能否更好地保留真實(shí)信號(hào),本文對(duì)不同濾波半徑的傳統(tǒng)高斯濾波、傳統(tǒng)組合濾波以及混合半徑高斯濾波、去相關(guān)濾波[12]與混合半徑高斯濾波的改進(jìn)組合濾波方法的信噪比進(jìn)行對(duì)比,結(jié)果如圖2所示?？梢钥闯?混合半徑高斯濾波方法在不同濾波半徑下信噪比均處于較優(yōu)水平。此外,隨著傳統(tǒng)高斯濾波半徑的增大,因?yàn)V波造成的信號(hào)損失也越嚴(yán)重,而利用本文的混合半徑高斯濾波方法對(duì)信號(hào)保留效果較好。

圖2 不同濾波方法信噪比Fig.2 SNR of different filtering methods

圖3展示了利用GRACE反演全球陸地水儲(chǔ)量變化的效果,可以看出,不進(jìn)行濾波或只進(jìn)行去相關(guān)濾波都無(wú)法完全提取真實(shí)的等效水高信號(hào),而經(jīng)過(guò)經(jīng)典高斯濾波、經(jīng)典組合濾波以及混合半徑高斯濾波、Duan等[12]與混合半徑高斯濾波的改進(jìn)組合濾波均對(duì)條帶噪聲有明顯壓制作用,真實(shí)等效水高信號(hào)已經(jīng)顯現(xiàn)。

圖3 全球陸地水儲(chǔ)量變化Fig.3 Changes in global land water reserves

同時(shí),對(duì)比分析圖3中高斯濾波(c)和混合半徑高斯濾波(d)、組合濾波(e)和Duan等[12]與混合半徑高斯濾波的改進(jìn)組合濾波(f)兩組全球陸地水儲(chǔ)量變化中水儲(chǔ)量信號(hào)變化較劇烈的亞馬遜流域、剛果盆地地區(qū)的等效水高信號(hào)可以看出,經(jīng)典高斯濾波方法對(duì)信號(hào)平滑效果更強(qiáng),尤其在剛果盆地地區(qū)最為明顯,而對(duì)信號(hào)過(guò)度平滑則會(huì)造成區(qū)域信號(hào)泄漏。相比之下,混合半徑高斯濾波在保證去條帶誤差的同時(shí)還能較好地減弱信號(hào)泄漏,提高區(qū)域反演結(jié)果的準(zhǔn)確性,這也再次印證改進(jìn)高斯濾波能更好地保留水儲(chǔ)量變化信號(hào)的判斷,證明改進(jìn)高斯濾波具有可行性。

3 結(jié) 語(yǔ)

本文對(duì)GRACE時(shí)變重力場(chǎng)條帶誤差的處理方法進(jìn)行分析,著重討論時(shí)變重力場(chǎng)球諧系數(shù)不同階數(shù)下不同半徑的高斯濾波以及組合濾波的信噪比情況。結(jié)果表明,在低于20階時(shí),球諧系數(shù)具有較高信噪比,因此不同半徑高斯濾波效果近似;當(dāng)球諧系數(shù)大于30階時(shí)信噪比較低,因此需要用濾波手段提高信噪比。

通過(guò)不同高斯濾波半徑下的經(jīng)典高斯濾波以及組合濾波信噪比的分析可以發(fā)現(xiàn),僅進(jìn)行高斯濾波時(shí),400 km半徑高斯濾波信噪比最優(yōu),在各階下處于最好水平;在進(jìn)行Duan等[12]去相關(guān)濾波和高斯濾波組合濾波時(shí),由于消除了球諧系數(shù)的相關(guān)性誤差,因此在進(jìn)行300 km半徑高斯濾波時(shí),球諧系數(shù)各階下信噪比處于最好水平。

對(duì)得到的最優(yōu)高斯濾波半徑進(jìn)行改進(jìn),結(jié)果表明,僅進(jìn)行經(jīng)典高斯濾波時(shí),分兩步進(jìn)行的400 km混合半徑的高斯濾波信噪比較經(jīng)典高斯濾波更優(yōu);進(jìn)行組合濾波時(shí),Duan等[12]與分3步進(jìn)行的300 km混合半徑高斯濾波信噪比更優(yōu)。對(duì)全球陸地水儲(chǔ)量變化較為劇烈的亞馬遜流域、剛果盆地地區(qū)信號(hào)分析可以看出,混合半徑高斯濾波能更好地保留真實(shí)信號(hào),減弱因高斯平滑造成的區(qū)域信號(hào)泄漏。