小波分析和經(jīng)驗?zāi)B(tài)分解對BDS多路徑誤差削弱對比研究

陶遠(yuǎn),鄧永春,胡豪杰,高翔

(1.安徽理工大學(xué) 測繪學(xué)院,安徽 淮南 232001; 2.礦山采動災(zāi)害空天地協(xié)同監(jiān)測與預(yù)警安徽省教育廳重點實驗室,安徽 淮南 232001)

0 引 言

全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS),作為高精度定位技術(shù)被廣泛應(yīng)用在大地測量學(xué)、地球動力學(xué)、地質(zhì)學(xué)、氣象學(xué)和定位導(dǎo)航授時等技術(shù)[1-4]．隨著工程領(lǐng)域?qū)NSS測量精度要求不斷的提升,眾多學(xué)者對觀測數(shù)據(jù)中各項誤差處理的方式也有不斷的改進(jìn)[5]．GNSS測量中存在較多誤差,如軌道誤差、接收機鐘差、衛(wèi)星鐘差、對流層延遲、電離層延遲、濕延遲、多路徑誤差以及噪聲等,這些誤差在短基線的雙差觀測方程運算中能夠大部分消除[6-8]．但是,外部環(huán)境的改變以及衛(wèi)星的運行會導(dǎo)致多路徑誤差的變化,難以建立一個固有模型來解決不同環(huán)境下的多路徑誤差[9]．北斗衛(wèi)星導(dǎo)航系統(tǒng)(BDS)也存在同樣的難題,眾多學(xué)者紛紛投入到解決BDS的多路徑誤差研究中[10]．

目前,在GNSS實時監(jiān)測應(yīng)用中,由于受到多路徑誤差等效應(yīng)影響,相對定位精度一般為cm級,而要實現(xiàn)mm級乃至亞mm級的高精度變形監(jiān)測,必須限制多路徑效應(yīng)影響[11]．多路徑效應(yīng)會對觀測信號的傳播距離及路徑、相位和幅度產(chǎn)生影響,并且降低信噪比．本文使用改進(jìn)恒星日濾波(ASF)和經(jīng)驗?zāi)B(tài)分解(EMD)方法對原始觀測序列進(jìn)行數(shù)據(jù)處理,對兩種方法的優(yōu)劣勢進(jìn)行分析對比．

1 多路徑產(chǎn)生機理

在GNSS測量中,被測站附近的反射源所反射或衍射的干擾信號被接收機天線接收,將會和直接信號產(chǎn)生干涉,從而使觀測值產(chǎn)生偏差或失真,產(chǎn)生多路徑誤差．這種干擾信號可能是來自于墻面、斜坡、水體和植被等因素的影響,這種由多路徑的信號傳播所引起的干涉時延效應(yīng)被稱作多路徑效應(yīng)．多路徑效應(yīng)將嚴(yán)重?fù)p害GNSS測量的精度,甚至將產(chǎn)生信號的失鎖．多路徑誤差是GNSS測量中一個重要的誤差源．

從產(chǎn)生多路徑誤差的機理可以看出,多路徑效應(yīng)相對于接收機天線的空間關(guān)系以及周圍地物的環(huán)境有關(guān)．

數(shù)學(xué)形式表達(dá)為

SD=α0Acos(φ),

(1)

SM=α0Acos(φ+Δ0).

(2)

式中:SM為多路徑等效距離誤差,即為圖1中AB與BC之和,且SM=2Dcosβ,其中：D為接收機天線至反射面垂直距離；β為衛(wèi)星信號入射角;α0為反射面的反射系數(shù);A為衛(wèi)星信號振幅;φ為衛(wèi)星信號相位;Δ0為反射信號引起的相位偏移量,且如式(3),其中,λ為衛(wèi)星信號的波長．

(3)

圖1 多路徑效應(yīng)示意圖

可以得到直射信號與反射信號共同構(gòu)成的信號為:

ST=SD+SM=ATcos(φ+ψ),

(4)

式中,合成信號的振幅AT與多路徑誤差造成的相位延遲量ψ可分別表示為

AT=A(1+2αcos Δ0+α2)1/2,

(5)

(6)

從公式(6)可得,當(dāng)反射面不變時,反射系數(shù)α為固定值,將公式(6)對Δ0求導(dǎo)

(7)

多路徑效應(yīng)的大小與衛(wèi)星高度角有關(guān),衛(wèi)星高度越低,影響就越大．多路徑誤差具有周日重復(fù)性,這是由于衛(wèi)星位置、反射源與接收機天線三者的空間位置密切相關(guān)所決定的,但是多路徑效應(yīng)在各站之間沒有相關(guān)性．當(dāng)部分反射源的距離大于50 m時,干擾信號基本可以忽略不計．通過對多路徑誤差頻譜特性的研究,驗證了近距離反射源多路徑誤差呈低頻性,遠(yuǎn)距離多路徑誤差呈高頻性[12-13]．這一特性導(dǎo)致了多路徑誤差難以建立固定確切的數(shù)學(xué)模型進(jìn)行改正或消除．

2 Wavelet、EMD和ASF原理

2.1 Wavelet方法基本原理

20世紀(jì)80年代隨著函數(shù)的伸縮和平移的出現(xiàn),小波分析及其應(yīng)用逐漸出現(xiàn)．小波分析方法(Wavelet)隨著Fourier變換的發(fā)展逐步廣泛為人們熟知．小波變換是一種能夠?qū)NSS時間序列時域和頻域同時進(jìn)行分析的方法,小波分析方法在某種程度上通過函數(shù)的伸縮變換改變GNSS時間序列的頻域改變時間序列的時間窗大小,從而達(dá)到分離噪聲的目的．小波降噪最早是由Weaver等人開創(chuàng)性提出,在此基礎(chǔ)上Donoho教授系統(tǒng)性地對小波降噪進(jìn)行研究,完善小波的降噪體系．小波變換可以通過將GNSS時間序列分析逐層分解提取有用的信息從而達(dá)到降噪的目的．

GNSS時間序列小波降噪的具體步驟如圖2所示,介紹如下:

1)小波分解:選擇合適的小波基,對含噪的GNSS時間序列小波分解,得到各個尺度上的小波系數(shù), 目前常用的小波函數(shù)有db系列小波、bior系列小波、morlet小波、sym系列小波、harr小波函數(shù)等;

2)作用閾值過程:對各個尺度上的信號進(jìn)行非線性閾值降噪處理,濾去小波系數(shù)為0的小波項,常用的閾值選取規(guī)則有:硬閾值規(guī)則、軟閾值規(guī)則、固定閾值規(guī)則(sqtwolog)、啟發(fā)式閾值規(guī)則(heursure)、極大極小原理(minimaxi)和使用Stein無偏似然估計原理的自適應(yīng)規(guī)則(rigrsure);

3)小波重構(gòu):降噪處理后的GNSS時間序列小波系數(shù),得到降噪后的GNSS時間序列．

圖2 小波閾值降噪流程圖

2.2 EMD基本原理

EMD的基本思想是,將復(fù)雜的信號看成一些互異的、簡單的、非正弦函數(shù)的多尺度信號分量．根據(jù)這一特性可將復(fù)雜的信號分解成高頻到低頻的若干信號分量,即固有模態(tài)函數(shù)(IMF)．因此每個IMF具有以下兩個特性:

1)極大值點與極小值點數(shù)目之和與零點數(shù)目之差相等或至多相差一個(簡稱零點條件)；

2)由局部極大值構(gòu)成的上包絡(luò)和由局部極小值構(gòu)成的下包絡(luò)的平均值為零(簡稱均值條件).

EMD算法的核心是篩選分解得到IMF的過程,對于信號x(t)的分解可以表示為

(8)

式中:imfi(t)為第i個IMF分量;rn(t)為單調(diào)殘差序列．

EMD的具體流程如下:

1)對于輸入的待分解信號x(t),尋找其中所有局部極值點;

2)利用樣條插值對信號x(t)中的極大值點(極小值點)建立上包絡(luò)線xu(t)(下包絡(luò)線xd(t));

3)計算上下包絡(luò)線的均值:

(9)

4)提取細(xì)節(jié):

h(t)=x(t)-xa(t).

(10)

式中：若h(t)滿足均值條件,均值為零,則此時h(t)即是滿足條件的一個IMF分量,否則對h(t)重復(fù)步驟1)～4),直至h(t)滿足均值條件;

5)計算殘余信號分量ri(t):

ri(t)=x(t)-imfi(t).

(11)

如果ri(t)的極值點數(shù)大于2,則跳轉(zhuǎn)至步驟1),按照上述步驟繼續(xù)進(jìn)行分解,否則分解結(jié)束．

各個IMF分量的能量密度與平均周期的乘積為固定常數(shù),從該性質(zhì)出發(fā),設(shè)計一種自適應(yīng)選擇IMF分量的具體方法,即在分解出自高頻到低頻的i個IMF分量和1個殘余項,計算RP系數(shù),以此判斷信號重構(gòu)尺度[14]．具體步驟如下:

1)計算每個IMF分量的能力密度及其平均周期之積:

(12)

2)計算RPj系數(shù):

(13)

當(dāng)RPj≥1時,第j個IMF分量的Pj相當(dāng)于前j-1個IMF的Pj的平均值成倍增大,即認(rèn)為前j-1個IMF分量的能量密度與平均周期之積為固定常量,第j個IMF分量即為重構(gòu)信號的尺度．

2.3 ASF基本原理

多路徑誤差主要由衛(wèi)星、接收機天線和反射源三者的幾何關(guān)系決定,而在動態(tài)變形監(jiān)測中,認(rèn)為地表發(fā)生較為微小的形變且周圍環(huán)境較為穩(wěn)定,因此,多路徑誤差主要與衛(wèi)星運行軌道有關(guān)．雖然BDS衛(wèi)星有三類不同的運行軌道衛(wèi)星,中軌道衛(wèi)星(MEO)、傾斜地球同步軌道衛(wèi)星(IGSO)和地球靜止軌道衛(wèi)星(GEO),相對于接收機在空間中不斷地運行變化,存在微小的軌道誤差,但運行軌道基本穩(wěn)定不變．

根據(jù)開普勒第三定律,衛(wèi)星運行軌道的周期的平方與其橢圓軌道長半軸的立方成正比,可以得到衛(wèi)星運行周期的平方與其橢圓軌道長半軸之比為一常數(shù)．根據(jù)各系統(tǒng)衛(wèi)星廣播星歷所提供的軌道參數(shù),計算得到各衛(wèi)星的運行周期,即幾何重復(fù)期,稱為廣播星歷法．廣播星歷中提供軌道長半軸平方根aS與平均角速度改正參數(shù)dn則由開普勒第三定律,根據(jù)aS、dn求得衛(wèi)星運行平均角速度:

(14)

式中:G為萬有引力常數(shù);M為地球質(zhì)量,參考值為

GM=3.986005×1014m3/s2

衛(wèi)星沿軌道運行一周的時間TS為

(15)

BDS衛(wèi)星的幾何重復(fù)周期為衛(wèi)星運行兩周的時間T0:

T0=2TS.

(16)

由此可以得到BDS衛(wèi)星第二天相對于第一天提前時間t為

t=86400-T0.

(17)

這種由于相鄰兩天衛(wèi)星運行提前,造成固定時刻多路徑誤差周期以一個固定周期平移的特性,就稱為多路徑誤差的周日重復(fù)性．利用周日重復(fù)性建立模型達(dá)到消除多路徑誤差的方法叫做ASF,因此恒星日濾波至少需要兩天的觀測數(shù)據(jù),通過提取第一天多路徑誤差建立模型,利用周日重復(fù)性將模型應(yīng)用于之后的觀測數(shù)據(jù)中達(dá)到消除多路徑誤差的目的．

通過最大相關(guān)系數(shù)法求得兩天對應(yīng)的時間延遲或者利用IGS站提供的廣播星歷能夠計算出衛(wèi)星的回歸周期根據(jù)衛(wèi)星準(zhǔn)確的回歸周期和每天準(zhǔn)確的時間延遲,建立恒星日濾波,稱為ASF[15]．ASF能夠充分利用多路徑誤差的周日重復(fù)性進(jìn)行多路徑誤差建模．但是,頻繁地利用廣播星歷計算衛(wèi)星的回歸周期較為繁瑣且效率低下．可使用最大相關(guān)系數(shù)法,計算出準(zhǔn)確的相鄰年積日(DOY)的時間延遲,較為高效地進(jìn)行多路徑誤差剔除．

相關(guān)系數(shù)計算公式為

(18)

式中:rxx(0)與ryy(0)分別為x序列與y序列的方差;rxy(l)為x序列與延遲l歷元的y序列的協(xié)方差:

(19)

式中,n為序列長度．

3 實驗分析

3.1 數(shù)據(jù)采集和分析

實驗數(shù)據(jù)采集于安徽理工大學(xué)某學(xué)院樓頂,于2018年DOY203-204兩天采用兩臺中海達(dá)V9接收機并構(gòu)建長約44 m的短基線進(jìn)行觀測．天氣狀況穩(wěn)定,微風(fēng),采樣間隔1 Hz,衛(wèi)星截止高度角10°．對觀測數(shù)據(jù)進(jìn)行動態(tài)基線解算,得到BDS系統(tǒng)的三維動態(tài)基線坐標(biāo)E、N、U方向發(fā)坐標(biāo)序列,使用長期靜態(tài)觀測的坐標(biāo)解為真值并進(jìn)行數(shù)據(jù)驗證．

由于該基線為短基線且處于同一樓頂上,顧及在高程方向上無明顯差異,可認(rèn)為無對流層延遲和電離層延遲,且使用同一款接收機,可認(rèn)為無天線相位偏差．綜上,認(rèn)為該原始觀測數(shù)據(jù)僅存在噪聲和多路徑誤差．

3.2 數(shù)據(jù)實驗

本文分別使用Wavelet和EMD方法對兩天原始BDS數(shù)據(jù)的E、N、U方向進(jìn)行降噪,基于多路徑誤差的周日重復(fù)性對降噪后的坐標(biāo)序列求取最大相關(guān)系數(shù),最大相關(guān)系數(shù)對應(yīng)的時間延遲即為ASF的時間延遲,使用公式(14)計算DOY204對DOY203的E、N、U對應(yīng)的相關(guān)系數(shù)曲線和最大相關(guān)系數(shù)對應(yīng)的時間延遲值如表1所示,圖3示出了兩種方法降噪后DOY204與DOY203的E方向的相關(guān)系數(shù)曲線．對兩天的降噪后序列使用ASF方法,得到剔除多路徑誤差的坐標(biāo)序列．

表1 坐標(biāo)序列的影響系數(shù)

圖3 DOY203與DOY204降噪序列的相關(guān)系數(shù)曲線

本文小波方法使用heursure小波coif4小波基對兩天的原始觀測序列降噪與EMD降噪,其中DOY203降噪結(jié)果如圖4所示．

圖4 DOY203原始觀測序列和降噪后序列

使用EMD對原始序列的E方向坐標(biāo)序列進(jìn)行多尺度分解,建立多路徑誤差模型,由公式(11)可知，其中第14個IMF分量為殘余信號分量,多尺度分解得模態(tài)分量結(jié)果并基于公式(13)進(jìn)行RP系數(shù)運算,計算的RP系數(shù)在第5個IMF分量大于1,從第5個IMF分量到第13個IMF分量進(jìn)行模態(tài)重構(gòu),即為EMD降噪后的多路徑誤差模型,分量結(jié)果如圖5所示．EMD降噪后的序列如圖6所示．

圖5 DOY204的E方向坐標(biāo)序列的EMD分量結(jié)果

第一天的降噪序列作為多路徑誤差模型進(jìn)行第二天的多路徑誤差剔除,提取對應(yīng)降噪后的多路徑誤差模型進(jìn)行坐標(biāo)序列改正,Wavelet對E、N、U方向的精度改善程度分別為38.6%、59.1%和57.8%,改正后坐標(biāo)序列如圖7所示．

由EMD降噪并使用ASF剔除多路徑誤差后的坐標(biāo)殘差序列,相比Wavelet降噪方法并剔除多路徑誤差得到的坐標(biāo)殘差序列更加穩(wěn)定,如圖6所示．同時,EMD多尺度分解對精度改善程度較為明顯,該方法對E、N、U改善程度分別為40.8%、61.0%和57.9%,如表1所示．

圖6 EMD剔除多路徑后第2天坐標(biāo)序列

圖7 Wavelet剔除多路徑后第2天坐標(biāo)序列

EMD與小波相比,前者能得到的坐標(biāo)序列較為平滑穩(wěn)定,能夠較好地改善E、N、U坐標(biāo)序列的精度．小波分解并不能剔除瞬時強噪聲,而EMD能夠較大程度地削弱隨機噪聲．

4 結(jié)束語

通過計算兩種方法的RMSE可知,在使用EMD剔除多路徑誤差后的E、N、U方向坐標(biāo)序列穩(wěn)定性和Wavelet方法均有提高,較優(yōu)且有效地削弱多路徑誤差,坐標(biāo)精度得到提升,可提高BDS的定位精度;EMD方法能夠較大程度地削弱觀測噪聲,留下“干凈”的坐標(biāo)序列,而且在序列的平滑程度上EMD降噪方法更有優(yōu)勢．由于BDS星座異構(gòu),GEO和IGSO的運行回歸周期與MEO相差較大,但是MEO在接收機視野范圍內(nèi)只有2～3 h,不會對整體的周日重復(fù)性有明顯影響．