形態(tài)濾波在單頻周跳探測(cè)中的應(yīng)用*

池傳國(guó)， 姜 迪， 黃國(guó)勇， 吳建德， 馬 軍

(1.昆明理工大學(xué) 信息工程與自動(dòng)化學(xué)院，云南 昆明 650500；2.云南省礦物管道輸送工程技術(shù)研究中心，云南 昆明 650500；3.昆明理工大學(xué) 機(jī)電工程學(xué)院，云南 昆明 650500)

0 引 言

周跳探測(cè)是北斗衛(wèi)星導(dǎo)航系統(tǒng)(Beidou system,BDS)載波相位高精度定位必須要解決的問題之一[1]。國(guó)內(nèi)外相關(guān)研究領(lǐng)域主要有：高次殘差法、偽距載波相位差法、Blewitt法、無(wú)幾何距離法和多項(xiàng)式擬合法等[2]。但是在上述已有方法中，較有效的探測(cè)方法如無(wú)幾何距離法和Blewitt法等均針對(duì)雙頻觀測(cè)數(shù)據(jù)，對(duì)于單頻數(shù)據(jù)因無(wú)法構(gòu)成無(wú)幾何組合以及MW(Melbourne-Wubbena)組合[3]，基于此，為了快速實(shí)現(xiàn)對(duì)北斗單頻周跳的探測(cè)，采用偽距相位差法與奇異值分解(singular value decomposition,SVD)形態(tài)濾波相結(jié)合的周跳探測(cè)方法，提取周跳的特征信號(hào)，進(jìn)行周跳檢測(cè)。實(shí)驗(yàn)結(jié)果表明：該方法對(duì)可以對(duì)小周跳實(shí)現(xiàn)精確探測(cè)，有效提高了探測(cè)周跳發(fā)生位置的能力[4]。

1 周跳檢測(cè)量與檢測(cè)原理

周跳的產(chǎn)生主要是因?yàn)槎嗥绽沼?jì)數(shù)的短暫中止，一般可導(dǎo)致這種中止的情況主要有：信號(hào)在傳遞過(guò)程中被遮擋;接收機(jī)和衛(wèi)星間的其他影響因素導(dǎo)致的低信噪比;接收機(jī)信號(hào)處理出現(xiàn)問題;衛(wèi)星采集信息時(shí)出現(xiàn)故障。信號(hào)一旦中斷，載波鎖相環(huán)便會(huì)產(chǎn)生短暫失鎖，導(dǎo)致多普勒計(jì)數(shù)中止。當(dāng)跟蹤信號(hào)恢復(fù)后，多普勒計(jì)數(shù)恢復(fù)，期間出現(xiàn)的間斷導(dǎo)致的信號(hào)不連續(xù)即為周跳。

本文首先采用相位偽距二次差法處理原始數(shù)據(jù)[5]。偽距組合觀測(cè)量可定義為

Rabc=aR1+bR2+cR3=ρ+T+δr+ηabcI1+εabc

(1)

載波相位觀測(cè)方程為

λijkφijk=ρ+T+δr-λijkNijk-ηijkI1+λijkεijk

(2)

用偽距減相位法可得周跳估值，即將式(1)減式(2)

Nijk,abc=

(3)

將式(3)在歷元間作差，得到周跳檢驗(yàn)量

(4)

式中ηijk,abc為偽距相位電離層影響系數(shù)；Δφijk為載波相位觀測(cè)量的變化量；Δεijk為載波相位噪聲變化量；ΔI1為電離層延遲變化量[5]。當(dāng)采樣率較高時(shí)，ΔI1可以忽略不計(jì)；選取較小的ηijk,abc可減弱電離層I1的影響，則周跳檢驗(yàn)量為

(5)

周跳估值標(biāo)準(zhǔn)差為

(6)

式中σε為三頻載波相位觀測(cè)噪聲標(biāo)準(zhǔn)差，周；σe為三頻偽距觀測(cè)噪聲標(biāo)準(zhǔn)差， m。

2 SVD形態(tài)濾波算法

2.1 構(gòu)造Hankel矩陣

首先對(duì)由測(cè)量數(shù)據(jù)構(gòu)造的Hankle矩陣進(jìn)行SVD，對(duì)Hankel矩陣取不同維數(shù)時(shí)對(duì)空間劃分和濾波效果的影響進(jìn)行分析和研究[6]。

假定X=[x1,x2,…,xn]為加入周跳的檢測(cè)信號(hào)，使用此信號(hào)構(gòu)造Hankel矩陣如下

(7)

式中 1

選擇的維數(shù)過(guò)高，可能使噪聲信號(hào)無(wú)法較大程度地濾除；過(guò)低時(shí)，可能使得有效信號(hào)被濾除，甚至可能導(dǎo)致信號(hào)波形畸變。因此，本文采用奇異值曲線分析綜合衡量分量信號(hào)中所包含的信息量，從而有效確定Hankel矩陣維數(shù)大小[7]。

2.2 峭度準(zhǔn)則

峭度作為無(wú)量綱參數(shù)可以描述波形尖峰度同時(shí)反映信號(hào)的分布特性，其數(shù)學(xué)描述為

(8)

式中μ為信號(hào)X的均值；σ為信號(hào)X的標(biāo)準(zhǔn)方差[8]。當(dāng)信號(hào)中出現(xiàn)周跳時(shí)，非正常信號(hào)的概率密度會(huì)增加，當(dāng)信號(hào)幅值分布明顯偏離正態(tài)分布時(shí)，峭度值會(huì)隨之增大。 峭度值越大，說(shuō)明信號(hào)中非正常成分所占比重越多，而周跳信息往往包含在這些非正常成分以及由此引起的幅值調(diào)制信號(hào)中。

2.3 形態(tài)濾波信號(hào)選擇

定義在離散域F={0，1，…，N-1}和G={0，1，…，M-1} 的實(shí)函數(shù)f(n)和g(n)，其中f(n) 為原始輸入信號(hào)，g(n)為結(jié)構(gòu)元素。f(n)和g(n)之間存在4種基本運(yùn)算分別為膨脹運(yùn)算、腐蝕運(yùn)算、開運(yùn)算,閉運(yùn)算。

假設(shè)采集到的一維多值信號(hào)f(n)定義域?yàn)镈[f]={0,1,2,3,…,N} ，選擇結(jié)構(gòu)元素序列為g(x)，且定義域?yàn)镈[g]={0,1,2,3,…,P}，P和N為整數(shù)。分別定義腐蝕與膨脹運(yùn)算如下

(fΘg)(n)=min{f(n+x)}-g(x);

x∈D[g],n=(1,2,…,N)}

(9)

(f⊕g)(n)=min{f(n-x)}+g(x);

x∈D[g],n=(1,2,…,N)}

(10)

由膨脹腐蝕運(yùn)算引出的形態(tài)學(xué)開、閉運(yùn)算表達(dá)式為

f°g=f?g⊕g

(11)

f·g=f⊕g?g

(12)

式中 “° ”為開運(yùn)算；“·”為閉運(yùn)算。

(13)

式中N為總采樣點(diǎn)個(gè)數(shù)；yi0為未疊加噪聲的輸入信號(hào)離散采樣點(diǎn)；yi為經(jīng)過(guò)形態(tài)學(xué)濾波器濾波后的輸出信號(hào)離散采樣點(diǎn)。E越小，說(shuō)明濾波后的信號(hào)與未加噪聲的信號(hào)越接近[9]。

濾波效果不僅和結(jié)構(gòu)元素的形狀有關(guān)還和結(jié)構(gòu)元素的尺寸有關(guān)。選取合適的尺寸，對(duì)于抑制信號(hào)內(nèi)部的細(xì)節(jié)差異會(huì)起到較好的效果，而且不會(huì)使邊界弱化。

3 實(shí)驗(yàn)仿真

實(shí)驗(yàn)采用某公司的雙星五頻北斗數(shù)據(jù)信號(hào)，選取B3頻段37 118 357.406～37 125 029.047之間連續(xù)的1 000組歷元。設(shè)置數(shù)據(jù)采樣頻率為1 000 Hz，采樣間隔為1 s，以MATLAB為支持實(shí)現(xiàn)實(shí)驗(yàn)仿真。所選數(shù)據(jù)為原始不含周跳的數(shù)據(jù) ，通過(guò)在所選數(shù)據(jù)的不同歷元加入不同周跳，模擬不同類型的周跳探測(cè)，檢測(cè)本文方法對(duì)周跳的探測(cè)范圍。通過(guò)對(duì)單周跳的高精度探測(cè)，進(jìn)行形態(tài)濾波實(shí)驗(yàn)仿真與傳統(tǒng)SVD實(shí)驗(yàn)仿真數(shù)據(jù)對(duì)比，以檢驗(yàn)本文方法的優(yōu)越性。

通過(guò)偽距與載波相位觀測(cè)值的二次差分處理，然后再在歷元間作差，得到無(wú)鐘差和噪聲干擾的實(shí)驗(yàn)用原始信號(hào)，如圖1所示。

圖1 未加周跳的原始信號(hào)波形

3.1 模擬不同類型的周跳探測(cè)

實(shí)驗(yàn)在不同歷元間加入了不同周跳進(jìn)行探測(cè)，如表1。

表1 多周跳在不同歷元間的添加情況

根據(jù)奇異值曲線，如圖2，可以看出m=6時(shí)奇異點(diǎn)趨于零，選取m=6構(gòu)建分析矩陣。

圖2 奇異值曲線

根據(jù)表1所加周跳對(duì)數(shù)據(jù)進(jìn)行處理，經(jīng)SVD得到如圖3(a)所示的6個(gè)分量信號(hào)，可以看出：濾波前的分量信號(hào)雖然存在周跳信號(hào)，但并不十分突出且存在噪聲信號(hào)干擾。采用形態(tài)濾波算法對(duì)分量信號(hào)進(jìn)行重新生成以達(dá)到最大程度的濾除噪聲信號(hào)突出周跳信號(hào)，如圖3(b)。

圖3 濾波前后分量信號(hào)

分量信號(hào)的峭度值提取結(jié)果分別為4.07，5.33，4.31，3.29，3.24，2.94。選取峭度值較大的分量信號(hào)2和分量信號(hào)3進(jìn)行信號(hào)重構(gòu),結(jié)果如圖4。

圖4 重構(gòu)信號(hào)波形

由圖可以看出，在所加5周和3周周跳處可以看出明顯的沖擊信號(hào)，所加1周周跳處也存在沖擊信號(hào)，但相對(duì)不明顯。實(shí)驗(yàn)清晰檢測(cè)出了5周和3周周跳，對(duì)于1周的周跳信號(hào)，雖然也能檢測(cè)，但幅值較小，效果不明顯。

3.2 單周跳的探測(cè)與對(duì)比

通過(guò)前述實(shí)驗(yàn)分析，在多個(gè)歷元加入周跳時(shí)，對(duì)1周的小周跳的探測(cè)效果并不明顯。為進(jìn)一步驗(yàn)證方法的可行性，在400歷元處加入2周周跳進(jìn)行探測(cè)，并采用傳統(tǒng)SVD進(jìn)行相同探測(cè)，通過(guò)對(duì)探測(cè)效果的比較檢驗(yàn)本方法的優(yōu)越性。如圖5為形態(tài)濾波后的信號(hào)波形。

圖5 形態(tài)濾波后分量信號(hào)波形

根據(jù)濾波后的分量信號(hào)，可提取峭度值信號(hào)，分別為4.01，5.32，5.65，4.19，3.47，10.5。圖6為本文方法對(duì)周跳信號(hào)的探測(cè)結(jié)果和傳統(tǒng)SVD方法對(duì)周跳的探測(cè)結(jié)果。

圖6 形態(tài)濾波和傳統(tǒng)SVD對(duì)周跳信號(hào)探測(cè)的波形

4 結(jié) 論

通過(guò)對(duì)多個(gè)歷元加入不同周跳和在單歷元加入單一周跳與傳統(tǒng)SVD探測(cè)方法進(jìn)行對(duì)比實(shí)驗(yàn)，驗(yàn)證形態(tài)濾波結(jié)合SVD算法較傳統(tǒng)SVD算法在周跳探測(cè)中的優(yōu)越性，實(shí)驗(yàn)結(jié)果表明：將形態(tài)濾波用于SVD分量篩選，能夠獲得更清晰的周跳信號(hào)，提高了周跳探測(cè)的精確性。

[1] 王貴文,王澤民,殷海濤.基于三差觀測(cè)量的實(shí)時(shí)動(dòng)態(tài)GPS周跳修復(fù)方法研究[J].武漢大學(xué)學(xué)報(bào)：信息科學(xué)版,2007,32(8):711-714.

[2] 雷 雨,高玉平.單頻非差相位的周跳檢測(cè)與修復(fù)方法研究[J].儀器儀表學(xué)報(bào),2013,34(11):2484-2490.

[3] 張 亮,岳東杰.相位減偽距法與電離層殘差法探測(cè)和修復(fù)周跳[J].測(cè)繪工程,2014,23(2):36-38.

[4] 滕云龍,師奕兵,鄭 植,等.單頻載波相位的周跳探測(cè)與修復(fù)算法研究[J].儀器儀表學(xué)報(bào),2010,31(8):1700-1705.

[5] 張 波,李健君.基于Hankel矩陣與奇異值分解(SVD)的濾波方法以及在飛機(jī)顫振試驗(yàn)數(shù)據(jù)預(yù)處理中的應(yīng)用[J].振動(dòng)與沖擊,2009,28(2):162-166.

[6] Feng S,Ochieng W,Moore T,et al.Carrier phase-based integrity monitoring for high-accuracy positioning[J].GPS Solutions,2009,13(1):13-22.

[7] 李天云,趙 妍,李 楠.基于EMD的Hilbert變換應(yīng)用于暫態(tài)信號(hào)分析[J].電力系統(tǒng)自動(dòng)化,2005,29(4):49-52.

[8] 胡愛軍,馬萬(wàn)里,唐貴基.基于集成經(jīng)驗(yàn)?zāi)B(tài)分解和峭度準(zhǔn)則的滾動(dòng)軸承故障特征提取方法[J].中國(guó)電機(jī)工程學(xué)報(bào),2012,32(11):106-111.

[9] 王家良，程春玲.一種多層自適應(yīng)形態(tài)濾波算法[J].計(jì)算機(jī)科學(xué)，2015，42(5)：72-77.