小波分析和經驗模態分解對BDS多路徑誤差削弱對比研究

陶遠,鄧永春,胡豪杰,高翔

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

0 引 言

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

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

1 多路徑產生機理

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

從產生多路徑誤差的機理可以看出,多路徑效應相對于接收機天線的空間關系以及周圍地物的環境有關．

數學形式表達為

SD=α0Acos(φ),

(1)

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

(2)

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

(3)

圖1 多路徑效應示意圖

可以得到直射信號與反射信號共同構成的信號為:

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

(4)

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

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

(5)

(6)

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

(7)

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

2 Wavelet、EMD和ASF原理

2.1 Wavelet方法基本原理

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

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

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

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

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

圖2 小波閾值降噪流程圖

2.2 EMD基本原理

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

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

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

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

(8)

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

EMD的具體流程如下:

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

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

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

(9)

4)提取細節:

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

(10)

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

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

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

(11)

如果ri(t)的極值點數大于2,則跳轉至步驟1),按照上述步驟繼續進行分解,否則分解結束．

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

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

(12)

2)計算RPj系數:

(13)

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

2.3 ASF基本原理

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

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

(14)

式中:G為萬有引力常數;M為地球質量,參考值為

GM=3.986005×1014m3/s2

衛星沿軌道運行一周的時間TS為

(15)

BDS衛星的幾何重復周期為衛星運行兩周的時間T0:

T0=2TS.

(16)

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

t=86400-T0.

(17)

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

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

相關系數計算公式為

(18)

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

(19)

式中,n為序列長度．

3 實驗分析

3.1 數據采集和分析

實驗數據采集于安徽理工大學某學院樓頂,于2018年DOY203-204兩天采用兩臺中海達V9接收機并構建長約44 m的短基線進行觀測．天氣狀況穩定,微風,采樣間隔1 Hz,衛星截止高度角10°．對觀測數據進行動態基線解算,得到BDS系統的三維動態基線坐標E、N、U方向發坐標序列,使用長期靜態觀測的坐標解為真值并進行數據驗證．

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

3.2 數據實驗

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

表1 坐標序列的影響系數

圖3 DOY203與DOY204降噪序列的相關系數曲線

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

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

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

圖5 DOY204的E方向坐標序列的EMD分量結果

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

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

圖6 EMD剔除多路徑后第2天坐標序列

圖7 Wavelet剔除多路徑后第2天坐標序列

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

4 結束語

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