利用CEEMDAN進(jìn)行GNSS-MR雪深反演

薛張芳 劉立龍 吳昊艦 張 志 劉睿國

1 桂林理工大學(xué)測(cè)繪地理信息學(xué)院，桂林市雁山街319號(hào)，5410062 廣西空間信息與測(cè)繪重點(diǎn)實(shí)驗(yàn)室，桂林市雁山街319號(hào)，541006

基于多路徑效應(yīng)的全球?qū)Ш叫l(wèi)星系統(tǒng)多路徑反射測(cè)量(GNSS multipath reflectometry，GNSS-MR)技術(shù)具有時(shí)間分辨率高、全天候、實(shí)時(shí)自動(dòng)化、覆蓋面積廣等優(yōu)點(diǎn)，在積雪深度探測(cè)方面具有較大潛力。張雙成等[1]利用實(shí)測(cè)GPS數(shù)據(jù)對(duì)基于信噪比(signal-to-noise, SNR)的GNSS-MR探測(cè)雪深的算法進(jìn)行了初步驗(yàn)證。Ozeki等[2]提出無幾何距離的線性組L4觀測(cè)值反演積雪深度，與SNR觀測(cè)量探測(cè)的結(jié)果較為一致。但上述研究均采用二次多項(xiàng)式擬合方法獲取殘差序列，無法準(zhǔn)確擬合SNR信號(hào)趨勢(shì)，同時(shí)由于衛(wèi)星高度角的增加，使天線接收的反射信號(hào)來自不同地表和位置，干擾有效反射信號(hào)，導(dǎo)致反演結(jié)果極易出現(xiàn)跳變現(xiàn)象。本文利用自適應(yīng)噪聲的完全集合經(jīng)驗(yàn)?zāi)B(tài)分解[3](complete ensemble empirical mode decomposition with adaptive noise, CEEMDAN)代替二次多項(xiàng)式獲取殘差序列，并以美國科羅拉多州NWOT測(cè)站數(shù)據(jù)為例進(jìn)行雪深反演，驗(yàn)證本文方法的有效性。

1 GNSS-MR反演雪深原理

GNSS天線接收機(jī)可以同時(shí)接收衛(wèi)星發(fā)射的直射信號(hào)Ad和經(jīng)地表面反射的反射信號(hào)Ar。圖1是基于信噪比觀測(cè)值的GNSS-MR雪深探測(cè)示意圖。

圖1 GNSS-MR雪深探測(cè)示意圖Fig.1 Diagram of GNSS-MR for snow depth

直射信號(hào)與反射信號(hào)的振幅存在如下關(guān)系：

Ad?Ar

(1)

記Ac為合成信號(hào)的振幅,φ為直射信號(hào)與經(jīng)地表積雪面反射的信號(hào)的夾角，其關(guān)系為[4]：

(2)

由于Ad?Ar，通常采用低階多項(xiàng)式擬合方法來去除趨勢(shì)項(xiàng)，獲得低衛(wèi)星高度角下SNR殘差序列[1]。多路徑效應(yīng)下的反射信號(hào)振幅可表示為：

(3)

式中，λ為載波波長(zhǎng)，θ為衛(wèi)星高度角，h為垂直反射距離。取f=2h/λ，對(duì)SNR殘差序列進(jìn)行Lomb-Scargle頻譜分析，可得到GPS 多路徑反射信號(hào)的頻率f，從而得到h；再結(jié)合圖1，可獲取積雪厚度。

2 CEEMDAN

CEEMDAN是對(duì)經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)改進(jìn)后得到的一種變換形式[5],其每個(gè)本征模態(tài)函數(shù)(intrinsic mode function, IMF)對(duì)應(yīng)一個(gè)頻率，具有各自的物理意義。因此，可通過CEEMDAN將原始SNR信號(hào)分解成各階IMF并重構(gòu)，從而獲取高質(zhì)量的殘差序列。

設(shè)原始SNR信號(hào)為x(t)，Ej(·)是對(duì)信號(hào)進(jìn)行EMD分解后的第j階IMF分量，ωi(t)為第i次添加的高斯白噪聲(i=1，2，…，I)，εk為幅值(其中k=0，1，…，K，K是模態(tài)總體數(shù)量)。CEEMDAN分解步驟如下：

1)與EEMD[6]算法相同，在x(t)中加入白噪聲εkωi(t)，對(duì)其進(jìn)行EMD分解得到多個(gè)IMF分量并求集合平均值。第1階最終分量IMF1(t)為：

(4)

2)當(dāng)k=1時(shí)，計(jì)算1階殘差r1(t):

r1(t)=x(t)-IMF1(t)

(5)

3)對(duì)1階殘差r1(t)添加經(jīng)EMD分解后的噪聲分量，得到一個(gè)新信號(hào)，即r1(t)+ε1E1[ωi(t)]，再對(duì)其分解，得到的集合平均值為該階的IMF最終分量：

IMF2(t)=

(6)

4)當(dāng)k=2，…，K時(shí),計(jì)算第k階殘差rk(t)：

首先，要多主體參與教學(xué)評(píng)價(jià)環(huán)節(jié)，學(xué)習(xí)者、管理者、企業(yè)、教師多角度、多方面進(jìn)行教學(xué)評(píng)價(jià)，對(duì)評(píng)價(jià)信息及時(shí)反饋，實(shí)行教學(xué)改進(jìn)，引導(dǎo)教學(xué)貼近大眾多樣化需求；

rk(t)=rk-1(t)-IMFk(t)

(7)

5)對(duì)k階殘差rk(t)添加噪聲分量εkEk[ωi(t)]，獲得一個(gè)新信號(hào)，計(jì)算集合平均值，得第k+1階的最終IMF為：

IMFk+1(t)=

(8)

6)重復(fù)第4步，直到獲得的殘差不可再被分解。殘差最后滿足：

(9)

x(t)最終被分解為K個(gè)IMF分量和1個(gè)殘余分量：

(10)

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

本文觀測(cè)數(shù)據(jù)來自美國科羅拉多州NWOT觀測(cè)站(https:∥www.unavco.org)，地理坐標(biāo)為(40°03′19.4″N,105°35′25.8″W) ，接收機(jī)類型為大地測(cè)量型Trimble NetRS，天線類型為TRM41249，測(cè)站高3 m，數(shù)據(jù)采集間隔為5 s。實(shí)測(cè)數(shù)據(jù)來自U-C氣象站(40°01′48″N,105°34′48″W)，其距NWOT測(cè)站2.2 km，由美國國家氣候中心提供(https:∥gis.ncdc.noaa.gov/maps/ncei/summaries/daily )。

本實(shí)驗(yàn)收集了NWOT測(cè)站2016-11～2017-03連續(xù)125 d的GPS觀測(cè)原始數(shù)據(jù)。為對(duì)比CEEMDAN算法和二次多項(xiàng)式擬合算法反演雪深的精度差異，以平均誤差(ME)、均方根誤差(RMSE)和相關(guān)系數(shù)(R)作為精度評(píng)價(jià)指標(biāo)：

(11)

限于篇幅，僅以2016-12-16 GPS 32號(hào)衛(wèi)星原始信噪比序列為例說明，截取高度角為5°～25°，波段為L(zhǎng)2，利用CEEMDAN對(duì)其分解，結(jié)果如圖2所示。分解獲得的殘余分量就是SNR信號(hào)的趨勢(shì)項(xiàng)，作為自適應(yīng)結(jié)果，克服了二次多項(xiàng)式擬合帶來的振幅異常情況。

圖2 CEEMDAN分解SNR序列Fig.2 Decomposition of SNR sequence by CEEMDAN

考慮到IMF分量與原信噪比序列的相關(guān)性[7]和振幅變化趨勢(shì)，表1列出了部分IMF分量與原始信號(hào)的相關(guān)系數(shù)?？梢钥闯?，對(duì)于不同年積日和不同衛(wèi)星，對(duì)應(yīng)的IMF個(gè)數(shù)及與原始信號(hào)的相關(guān)性也不同。經(jīng)過對(duì)比，在2016-12-16將IMF6和IMF7組合作為有效殘差序列，隨后對(duì)其組合進(jìn)行Lomb-Scargle譜分析得到頻率f。圖3是在2016-12-16用2種方法獲取的殘差序列，可以看出，二次多項(xiàng)式擬合方法在高度角正弦值為0.17～0.19之間出現(xiàn)高頻振蕩時(shí)沒能準(zhǔn)確地?cái)M合出SNR變化趨勢(shì)；而CEEMDAN方法提取的殘差序列更為光滑，無高頻振蕩部分。另外，該方法避免了模態(tài)混疊現(xiàn)象，能有效防止其他信號(hào)的干擾。圖4表示2種方法獲取信噪比殘差序列的Lomb-Scargle譜分析圖，橫軸表示天線相位中心到積雪表面的距離，縱軸表示SNR反射信號(hào)頻譜振幅，位于振幅最高峰值對(duì)應(yīng)的距離就是垂直反射距離。由圖可見，采用二次多項(xiàng)式方法時(shí)，最高峰值低于準(zhǔn)確峰值，使得反演結(jié)果出現(xiàn)負(fù)值；雪深實(shí)測(cè)值為0.61 m，本文方法反演值為0.597 m，更接近實(shí)測(cè)雪深值。

表1 不同IMF分量與原始信號(hào)的相關(guān)系數(shù)

圖3 2種方法獲得的殘差序列Fig.3 Residual sequences obtained by two methods

圖4 SNR殘差序列Lomb-Scargle譜分析圖Fig.4 Lomb-Scargle spectrum analysis diagramof SNR residual sequence

統(tǒng)計(jì)3顆GPS衛(wèi)星在2016-11-16～2017-03-21的觀測(cè)數(shù)據(jù)，衛(wèi)星運(yùn)行周期為11 h 58 min，1 d內(nèi)至少有2次經(jīng)過同一測(cè)站上空，因此會(huì)出現(xiàn)多個(gè)SNR數(shù)據(jù)，以均值作為當(dāng)天有效雪深值，2種方法的雪深反演結(jié)果如圖5所示。由圖5(a)看出，PRN28衛(wèi)星在2016-12-06、2016-12-13～14和2016-12-24～27雪深值均為負(fù)值，尤其在12-06雪深探測(cè)值為-0.18 m，與其他反演值相差太大，且在其他衛(wèi)星也出現(xiàn)該異?，F(xiàn)象，應(yīng)剔除該異常值，將其余雪深值作為有效雪深探測(cè)值。相比之下，圖5(b)反演結(jié)果更佳，有效避免了異常跳變的情況。圖6表示2種方法在多星融合下反演的雪深變化，可以看出，利用CEEMDAN進(jìn)行GNSS-MR雪深反演的雪深值與實(shí)測(cè)雪深有較高的吻合度。結(jié)合表2，相比二次多項(xiàng)式擬合方法，CEEMDAN均方根誤差降低了30.7%，與實(shí)測(cè)值相關(guān)系數(shù)為0.965，該結(jié)果進(jìn)一步說明本文方法用于雪深探測(cè)具有良好的精度。

圖5 2種方法多衛(wèi)星反演雪深與實(shí)測(cè)雪深的對(duì)比Fig.5 Comparison between actual snow depth and multi-satellite inversion of snow depthbased on two methods

圖6 2種方法多衛(wèi)星融合反演雪深與實(shí)測(cè)雪深Fig.6 Comparison between actual snow depth and multi-satellite fusion inversion results of snow depth based on two methods

表2 2種方法反演雪深與實(shí)際雪深的平均誤差與均方根誤差

4 結(jié) 語

1) CEEMDAN分解得到的殘余分量可準(zhǔn)確體現(xiàn)原始SNR信號(hào)的趨勢(shì)，由實(shí)驗(yàn)數(shù)據(jù)分析可知，該方法可有效改善二次多項(xiàng)式擬合方法出現(xiàn)的異常跳變現(xiàn)象。

2) 本文方法可精確地提取SNR反射信號(hào)與直射信號(hào)的干涉信號(hào)，反演結(jié)果與實(shí)測(cè)數(shù)據(jù)的相關(guān)系數(shù)大于0.96，RMSE比使用傳統(tǒng)二次多項(xiàng)式法降低30.7%，ME優(yōu)于3 cm，驗(yàn)證了該方法用于GNSS-MR雪深反演的優(yōu)越性。