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

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

0 引 言

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

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

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

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

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

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

SD=α0Acos(φ),

(1)

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

(2)

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

(3)

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

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

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

(4)

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

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

(5)

(6)

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

(7)

多路徑效應(yīng)的大小與衛(wèi)星高度角有關(guān),衛(wèi)星高度越低,影響就越大．多路徑誤差具有周日重復(fù)性,這是由于衛(wèi)星位置、反射源與接收機(jī)天線三者的空間位置密切相關(guān)所決定的,但是多路徑效應(yīng)在各站之間沒有相關(guān)性．當(dāng)部分反射源的距離大于50 m時(shí),干擾信號(hào)基本可以忽略不計(jì)．通過對(duì)多路徑誤差頻譜特性的研究,驗(yàn)證了近距離反射源多路徑誤差呈低頻性,遠(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時(shí)間序列時(shí)域和頻域同時(shí)進(jìn)行分析的方法,小波分析方法在某種程度上通過函數(shù)的伸縮變換改變GNSS時(shí)間序列的頻域改變時(shí)間序列的時(shí)間窗大小,從而達(dá)到分離噪聲的目的．小波降噪最早是由Weaver等人開創(chuàng)性提出,在此基礎(chǔ)上Donoho教授系統(tǒng)性地對(duì)小波降噪進(jìn)行研究,完善小波的降噪體系．小波變換可以通過將GNSS時(shí)間序列分析逐層分解提取有用的信息從而達(dá)到降噪的目的．

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

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

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

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

圖2 小波閾值降噪流程圖

2.2 EMD基本原理

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

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

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

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

(8)

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

EMD的具體流程如下:

1)對(duì)于輸入的待分解信號(hào)x(t),尋找其中所有局部極值點(diǎn);

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

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

(9)

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

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

(10)

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

5)計(jì)算殘余信號(hào)分量ri(t):

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

(11)

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

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

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

(12)

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

(13)

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

2.3 ASF基本原理

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

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

(14)

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

GM=3.986005×1014m3/s2

衛(wèi)星沿軌道運(yùn)行一周的時(shí)間TS為

(15)

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

T0=2TS.

(16)

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

t=86400-T0.

(17)

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

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

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

(18)

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

(19)

式中,n為序列長(zhǎng)度．

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

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

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

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

3.2 數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

圖4 DOY203原始觀測(cè)序列和降噪后序列

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

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

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

由EMD降噪并使用ASF剔除多路徑誤差后的坐標(biāo)殘差序列,相比Wavelet降噪方法并剔除多路徑誤差得到的坐標(biāo)殘差序列更加穩(wěn)定,如圖6所示．同時(shí),EMD多尺度分解對(duì)精度改善程度較為明顯,該方法對(duì)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)序列的精度．小波分解并不能剔除瞬時(shí)強(qiáng)噪聲,而EMD能夠較大程度地削弱隨機(jī)噪聲．

4 結(jié)束語(yǔ)

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