付 林 趙東明 付林威

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)性，與水稻、小麥、草莓、園林樹(shù)木、烤煙、無(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可見(jiàn),高斯濾波半徑為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)?；旌习霃礁咚篂V波方法的目的就是在相同濾波半徑下盡可能保留信號(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)泄漏。