經(jīng)驗(yàn)?zāi)B(tài)分解濾波法用于時(shí)變重力場(chǎng)去噪研究

艾尚校， 肖 云

（1.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，西安 710054； 2.地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710054；3.西安測(cè)繪研究所，西安 710054）

GRACE（Gravity Recovery and Climate Experiment）重力衛(wèi)星獲取了大量的全球尺度、高精度、高時(shí)空分辨率的地球重力場(chǎng)信息，為認(rèn)識(shí)和探究地球物質(zhì)質(zhì)量的空間分布及變化規(guī)律提供了新方法，特別是在監(jiān)測(cè)和分析陸地水儲(chǔ)量變化的研究領(lǐng)域取得了大量、優(yōu)異的研究成果［1-7］.

受衛(wèi)星觀測(cè)模式、載荷儀器誤差等因素的影響，GRACE時(shí)變重力場(chǎng)的高階球諧位系數(shù)誤差較大，在空域中表現(xiàn)為嚴(yán)重的南北條帶噪聲，不利于陸地水儲(chǔ)量監(jiān)測(cè)分析研究.因此，時(shí)變重力場(chǎng)濾波算法研究成為一個(gè)熱門研究領(lǐng)域，很多相關(guān)研究機(jī)構(gòu)和科學(xué)家投入大量精力研究時(shí)變重力場(chǎng)濾波方法，產(chǎn)出大量研究成果，實(shí)現(xiàn)有效削弱條帶噪聲目的，獲得可靠、高精度重力場(chǎng)時(shí)變信號(hào)，用于支撐陸地水儲(chǔ)量變化研究分析.

目前研究較成熟的時(shí)變重力場(chǎng)濾波方法可以分為兩類［8］：第一類方法是通過(guò)引入平滑核函數(shù)以降低高階次項(xiàng)位系數(shù)權(quán)重的空間平滑濾波方法，Wahr等［9］最早利用高斯平均核函數(shù)降低高階項(xiàng)位系數(shù)的權(quán)重，建立了高斯平滑濾波.隨后研究發(fā)現(xiàn)，高階項(xiàng)位系數(shù)中高次項(xiàng)位系數(shù)的誤差大于低次項(xiàng)位系數(shù)的誤差.于是Han 等［10］提出了與位系數(shù)階項(xiàng)和次項(xiàng)均相關(guān)的各向異性濾波，有效地提高了空間濾波方法的濾波能力.Zhang等［11］提出了Fan濾波，即對(duì)位系數(shù)階項(xiàng)和次項(xiàng)同時(shí)進(jìn)行高斯平滑濾波處理，是一種更為簡(jiǎn)單的各向異性濾波.這類空間平滑濾波方法在去除噪聲的同時(shí)會(huì)損失部分真實(shí)地球物理信號(hào)，選取的平滑半徑參數(shù)越大，損失的真實(shí)信號(hào)越多，空間分辨率也越低.第二類方法是去相關(guān)濾波方法，Swenson和Wahr［12］發(fā)現(xiàn)球諧系數(shù)中同次的奇（偶）項(xiàng)階存在相關(guān)誤差，并設(shè)計(jì)滑動(dòng)固定窗多項(xiàng)式擬合消除該誤差，即去相關(guān)濾波.Chambers［13］和Chen等［14］提出了對(duì)大于m次的位系數(shù)進(jìn)行n次多項(xiàng)式擬合的PnMm方法，Duan等［15］根據(jù)球諧系數(shù)中誤差分布特點(diǎn)提出了滑動(dòng)可變窗去相關(guān)濾波.這類去相關(guān)濾波方法在低緯度地區(qū)的去噪能力較弱.在實(shí)際時(shí)變重力場(chǎng)濾波應(yīng)用中，通常將空間濾波方法與去相關(guān)濾波方法組合使用，可獲得更好的濾除噪聲效果.盡管國(guó)內(nèi)外大量的相關(guān)研究取得了很好的成果，但是在時(shí)變重力場(chǎng)濾波方法的去噪能力、保留真實(shí)信號(hào)能力等方面還有潛力可以挖掘.基于此，本文引入了經(jīng)驗(yàn)?zāi)B(tài)分解法，應(yīng)用于時(shí)變重力場(chǎng)濾波處理，以期獲得更好的濾波效果，提高時(shí)變重力場(chǎng)數(shù)據(jù)的精度.

本文首先引入經(jīng)驗(yàn)?zāi)B(tài)分解方法，構(gòu)建了適用時(shí)變重力場(chǎng)濾波的函數(shù)模型，分析了技術(shù)可行性，進(jìn)一步與其他成熟時(shí)變重力場(chǎng)濾波方法進(jìn)行比較，評(píng)估本文提出的方法去噪能力、保留真實(shí)信號(hào)能力、反演準(zhǔn)確度等能力指標(biāo)，最后通過(guò)與陸面水文模型GLDAS（Global Land Data assimilation System）進(jìn)行對(duì)比，進(jìn)一步驗(yàn)證該方法提取陸地水儲(chǔ)量的有效性和可靠性.

1 原理模型

1.1 經(jīng)驗(yàn)?zāi)B(tài)分解原理

經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）是一種自適應(yīng)信號(hào)分解方法，將輸入信號(hào)分解成有限個(gè)不同特征尺度的固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF）分量和一個(gè)殘余分量［16］.每一個(gè)IMF分量都要滿足兩個(gè)條件：一是IMF中極值點(diǎn)的個(gè)數(shù)和過(guò)零點(diǎn)的個(gè)數(shù)相等或相差最大不超過(guò)1；二是由極大值和極小值構(gòu)造的上包絡(luò)線和下包絡(luò)線的均值為0［16］.給定一個(gè)待分解的信號(hào)為x(t)，其經(jīng)驗(yàn)?zāi)B(tài)分解過(guò)程如下.

步驟1：構(gòu)造包絡(luò)線.尋找待分解信號(hào)x(t)中的極大值點(diǎn)和極小值點(diǎn)，并通過(guò)極大值點(diǎn)利用三次樣條插值函數(shù)構(gòu)造上包絡(luò)線，用相同方法通過(guò)極小值點(diǎn)構(gòu)造下包絡(luò)線.

步驟2：信號(hào)分解.將信號(hào)x(t)減去上下包絡(luò)線的均值xˉ(t)，得到新信號(hào)y(t)，即

步驟3：條件判斷.若y(t)滿足IMF 的兩個(gè)條件，則認(rèn)定y(t)為首個(gè)IMF，并記為s1；若不滿足，將y(t)作為新的待分解信號(hào)，重復(fù)步驟1～3，直至得到第一個(gè)IMF分量，記為s1.

步驟4：迭代篩選.上述過(guò)程可視為IMF分量的篩選過(guò)程，此時(shí)將待分解信號(hào)x(t)減去y(t)，得到新的待分解信號(hào)g1(t)，即

對(duì)g1(t)重復(fù)IMF分量的篩選過(guò)程，得到第二個(gè)IMF分量，記為s2.如此反復(fù)進(jìn)行，直至gn(t)為單調(diào)信號(hào)時(shí)記為r(t)，經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)束.則待分解信號(hào)x(t)的函數(shù)表達(dá)式為：

式中：si表示信號(hào)x(t)分解出的第i個(gè)IMF分量，頻率由高到低排列，r(t)為分解殘余項(xiàng).

1.2 構(gòu)建適用時(shí)變重力場(chǎng)的濾波方法

未經(jīng)濾波處理的時(shí)變重力場(chǎng)數(shù)據(jù)在空域中存在明顯的南北條帶噪聲.其任意緯度帶剖面數(shù)據(jù)，可以看作是隨經(jīng)度λ變化的緯度帶信號(hào)x(λ) .對(duì)于緯度帶信號(hào)，噪聲主要集中在相對(duì)高的頻率波段，而真實(shí)的重力場(chǎng)信號(hào)集中在相對(duì)低頻率波段［12］.由經(jīng)驗(yàn)?zāi)B(tài)分解原理可知，IMF分量按照頻率由高到低排序，則時(shí)變重力場(chǎng)高頻噪聲集中在排序靠前的模態(tài)分量中，而排序靠后的分量是低頻真實(shí)重力場(chǎng)信號(hào)主導(dǎo)的模態(tài)分量.因此，只需確定區(qū)分高頻噪聲和低頻信號(hào)的分界模態(tài)分量sk，然后重構(gòu)信號(hào)分量就可以實(shí)現(xiàn)緯度帶信號(hào)的濾波處理.

本文使用模態(tài)相關(guān)分選準(zhǔn)則［17］作為選取分界模態(tài)分量的依據(jù)，首先計(jì)算緯度帶的原始信號(hào)x(λ)與每個(gè)模態(tài)分量之間的互相關(guān)系數(shù)，計(jì)算式如下：

式中：N為經(jīng)度采樣點(diǎn)的個(gè)數(shù)，對(duì)于1°×1°的時(shí)變重力場(chǎng)格網(wǎng)數(shù)據(jù)而言，N=360；λ表示經(jīng)度；x(λ)表示以經(jīng)度為參量的原始信號(hào)；xˉ為原始信號(hào)的平均值；si表示第i個(gè)模態(tài)函數(shù)分量；sˉi為第i個(gè)模態(tài)函數(shù)分量的平均值；R為互相關(guān)系數(shù).

上述所得的互相關(guān)系數(shù)中，第一個(gè)局部極小值點(diǎn)對(duì)應(yīng)的模態(tài)分量即為分界模態(tài)分量sk，則對(duì)于sk之后的模態(tài)分量即為緯度帶中的低頻信號(hào).重構(gòu)低頻信號(hào)分量，并引入包含著部分細(xì)節(jié)信號(hào)信息sk和sk-1分量，就得到濾波后的緯度帶信號(hào)x?(λ)：

2 時(shí)變重力場(chǎng)經(jīng)驗(yàn)?zāi)B(tài)分解試驗(yàn)

以2004 年11 月的時(shí)變重力場(chǎng)數(shù)據(jù)為例，分析時(shí)變重力場(chǎng)經(jīng)驗(yàn)?zāi)B(tài)分解的可行性和有效性.在時(shí)變重力場(chǎng)信號(hào)中選取0°N緯度帶、30°N緯度帶、40°N 緯度帶等幾個(gè)具有代表性的時(shí)變信號(hào)x(λ)，其信號(hào)特征見(jiàn)圖1.圖1（a）～圖1（c）分別表示0°N 緯度帶、30°N 緯度帶、40°N 緯度帶的信號(hào)，其中左側(cè)圖為緯度帶信號(hào)隨經(jīng)度的變化，右側(cè)為相應(yīng)的歸一化功率譜（下同）.由圖1 可見(jiàn)，不同緯度帶的噪聲呈現(xiàn)出不同特征，緯度越低，噪聲的頻率越高.

圖1 不同緯度帶信號(hào)及其歸一化功率譜Fig.1 The signals and normalized power spectra of different latitudes

以40°N 緯度帶為例分析上述方法的可行性和有效性.將信號(hào)x(λ)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，得到的結(jié)果如圖2所示.

圖2 經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)果Fig.2 Empirical mode decomposition results

從圖2 中可以看出，信號(hào)x(λ)分解后得到了9 個(gè)固有模態(tài)函數(shù)分量和1 個(gè)殘余項(xiàng).根據(jù)式（4）計(jì)算x(λ)與各分量之間的互相關(guān)系數(shù)，如圖3所示.

圖3 互相關(guān)系數(shù)曲線Fig.3 The cross correlation coefficient curve

由圖3 可知，固有模態(tài)函數(shù)分量s6為分界模態(tài)，則選擇s5及其之后的模態(tài)分量與殘余分量進(jìn)行重構(gòu)，就得到濾波后信號(hào).以本文方法對(duì)不同緯度帶信號(hào)進(jìn)行濾波，濾波結(jié)果如圖4所示.圖4（a）～圖4（c）分別表示經(jīng)過(guò)經(jīng)驗(yàn)?zāi)B(tài)分解濾波后的0°N緯度帶、30°N緯度帶、40°N緯度帶的信號(hào).對(duì)比圖1和圖4可以發(fā)現(xiàn)，圖4 左側(cè)圖呈現(xiàn)出低頻率特性，右側(cè)圖的信號(hào)也主要集中在低頻帶，因此，本文提出的經(jīng)驗(yàn)?zāi)B(tài)分解濾波對(duì)于不同緯度帶信號(hào)的高頻噪聲均有明顯的濾除效果.

圖4 不同緯度帶信號(hào)濾波后的特征Fig.4 The characteristics of denoised signals in different latitudes

3 應(yīng)用試驗(yàn)

利用本文方法對(duì)實(shí)際的時(shí)變重力場(chǎng)數(shù)據(jù)進(jìn)行濾波處理，分析其濾波效果.以2007年10月的數(shù)據(jù)為例，濾波前后結(jié)果如圖5所示，其中圖5（a）為未經(jīng)濾波處理的時(shí)變重力場(chǎng)數(shù)據(jù)，圖5（b）為EMD濾波后的結(jié)果.由圖5可見(jiàn)EMD濾波有效濾除了時(shí)變重力場(chǎng)的噪聲，濾波后可以明顯辨識(shí)出陸地上亞馬遜流域、奧里諾科河流域、恒河流域等地區(qū)的重力場(chǎng)時(shí)變信號(hào).

圖5 EMD濾波前后結(jié)果對(duì)比Fig.5 Comparison of results before and after EMD filtering

為進(jìn)一步分析本文提出的濾波方法在時(shí)變重力場(chǎng)濾波的有效性，選取P3M6多項(xiàng)式擬合（P3M6）和滑動(dòng)可變窗多項(xiàng)式擬合（Duan）兩種廣泛應(yīng)用的去相關(guān)濾波作為對(duì)比方法，在相同的高斯平滑半徑下對(duì)比分析各方法的濾波效果.本文采用兩個(gè)評(píng)價(jià)指標(biāo)比較各濾波方法效果，一是選取濾波后海洋上殘余信號(hào)的均方根（Root Mean Square，RMS）作為評(píng)價(jià)濾波方法去噪能力的指標(biāo)［1］，二是選擇濾波后陸地信號(hào)和海洋殘余信號(hào)的RMS比值（信噪比）作為評(píng)價(jià)濾波方法保留真實(shí)信號(hào)能力的指標(biāo)［18］.計(jì)算三種方法濾波結(jié)果的均方根指標(biāo)和信噪比指標(biāo)，結(jié)果如表1所示.

由表1可知，在相同的平滑半徑下，本文方法濾波后的海洋殘余信號(hào)均方根值最小，且信噪比值最高，說(shuō)明該方法的去噪能力和保留真實(shí)信號(hào)的能力均優(yōu)于其他兩種方法.另外，由表可見(jiàn)平滑半徑是影響濾波效果的重要因素，當(dāng)平滑半徑為300 km時(shí)，三種方法均取到了最高信噪比值，分別為2.90（P3M6）、3.01（Duan）和3.14（EMD）.選擇最優(yōu)信噪比條件，即平滑半徑設(shè)置為300 km，利用三種濾波方法分別對(duì)圖5（a）數(shù)據(jù)進(jìn)行處理，結(jié)果如圖6所示.

圖6（a）～圖6（c）分別為當(dāng)平滑半徑為300 km時(shí)P3M6、Duan和EMD濾波的結(jié)果.對(duì)比圖5（a）和圖6，三種濾波方法均較好地濾除了南北條帶噪聲，可以清晰地辨識(shí)時(shí)變重力場(chǎng)信息.其中本文濾波方法結(jié)果中殘余噪聲低于其他兩種方法，說(shuō)明本文濾波方法的去噪能力更強(qiáng)，與表1結(jié)論相吻合.為評(píng)價(jià)濾波結(jié)果之間的一致性，計(jì)算三種方法濾波結(jié)果的兩兩差值，并統(tǒng)計(jì)差值的均方根，結(jié)果分別為1.2 cm（P3M6和Duan作差）、1.3 cm（EMD和P3M6作差）和0.92 cm（EMD和Duan作差），證明本文濾波方法具有較高的可靠性.此外，相比于其他兩種方法，本文濾波方法在格陵蘭島西部和南美洲中部的信號(hào)泄漏誤差最小，反演結(jié)果的準(zhǔn)確度更高.

為了進(jìn)一步驗(yàn)證本文濾波方法的可靠性，選取亞馬遜流域、密西西比河流域和長(zhǎng)江流域作為研究區(qū)域，比較三個(gè)區(qū)域重力衛(wèi)星反演的時(shí)變重力場(chǎng)濾波結(jié)果與GLDAS 水文模型估計(jì)的水資源變化結(jié)果一致性.利用本文提出濾波方法及上述兩種濾波方法反演選定區(qū)域陸地水儲(chǔ)量變化，時(shí)間跨度為2003年1月至2014 年12 月，同時(shí)利用GLDAS 水文模型反演相同時(shí)段內(nèi)研究區(qū)域水儲(chǔ)量的變化，為保證空間分辨率一致性，對(duì)GLDAS 數(shù)據(jù)進(jìn)行300 km 高斯濾波.各流域水儲(chǔ)量變化的時(shí)間序列如圖7 所示.通過(guò)最小二乘法擬合估算各濾波方法反演結(jié)果的周年振幅和周年相位，并計(jì)算各濾波方法與水文模型反演結(jié)果之間的相關(guān)系數(shù)，結(jié)果如表2所示.

由圖7 和表2 可知，本文提出的濾波方法與其他兩個(gè)傳統(tǒng)成熟算法相比較，結(jié)果之間具有較高的一致性，時(shí)間序列曲線幾乎重合.在周年振幅和相位方面，本文濾波方法與Duan濾波結(jié)果基本一致，充分說(shuō)明本文方法的可靠性.進(jìn)一步與水文數(shù)據(jù)GLDAS反演的水儲(chǔ)量信號(hào)比較，發(fā)現(xiàn)本文方法結(jié)果與GLDAS 反演的水儲(chǔ)量相關(guān)性最高，其相關(guān)系數(shù)均大于0.86；但是重力衛(wèi)星反演的水儲(chǔ)量信號(hào)振幅明顯大于水文模型反演的水儲(chǔ)量信號(hào)振幅，究其原因是水文模型反演的水儲(chǔ)量?jī)H包含淺層地表水的變化，而重力衛(wèi)星反演的水儲(chǔ)量囊括了土壤水、地下水、湖泊、濕地等更綜合的水文變化信號(hào).

圖7 研究區(qū)域水儲(chǔ)量變化Fig.7 Water storage changes in study areas

表2 各濾波方法結(jié)果的周年振幅、相位和與GLDAS的相關(guān)系數(shù)Tab.2 Anniversary amplitudes，phases and correlation coefficients with GLDAS for different filter results

4 結(jié)語(yǔ)

本文引入了一種新的時(shí)變重力場(chǎng)空域?yàn)V波方法——經(jīng)驗(yàn)?zāi)B(tài)分解濾波法，介紹了經(jīng)驗(yàn)?zāi)B(tài)分解基本原理，構(gòu)建了匹配時(shí)變重力場(chǎng)濾波的濾波函數(shù)，分析了經(jīng)驗(yàn)?zāi)B(tài)分解濾除時(shí)變重力場(chǎng)條帶噪聲的有效性，得到主要結(jié)論如下：

1）本文提出時(shí)變重力場(chǎng)經(jīng)驗(yàn)?zāi)B(tài)分解濾波法是一種有效濾波方法，與現(xiàn)有成熟的濾波方法比較，有更優(yōu)的去噪能力.試驗(yàn)分析表明，經(jīng)驗(yàn)?zāi)B(tài)分解濾波方法的信噪比達(dá)到3.14，而P3M6濾波和Duan濾波的信噪比分別為2.90和3.01，因此本文濾波方法具有較好的保留真實(shí)信號(hào)能力.

2）本文提出方法與其他兩種成熟方法在反演亞馬遜流域、密西西比河流域和長(zhǎng)江流域等地區(qū)的水儲(chǔ)量方面具有很好一致性，與Duan濾波結(jié)果的周年振幅和相位基本相同.進(jìn)一步與水文數(shù)據(jù)GLDAS反演的水儲(chǔ)量信號(hào)比較，反演水儲(chǔ)量相關(guān)性最高，其相關(guān)系數(shù)均大于0.86.