形態濾波在單頻周跳探測中的應用*

池傳國， 姜 迪， 黃國勇， 吳建德， 馬 軍

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

0 引 言

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

1 周跳檢測量與檢測原理

周跳的產生主要是因為多普勒計數的短暫中止，一般可導致這種中止的情況主要有：信號在傳遞過程中被遮擋;接收機和衛星間的其他影響因素導致的低信噪比;接收機信號處理出現問題;衛星采集信息時出現故障。信號一旦中斷，載波鎖相環便會產生短暫失鎖，導致多普勒計數中止。當跟蹤信號恢復后，多普勒計數恢復，期間出現的間斷導致的信號不連續即為周跳。

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

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

(1)

載波相位觀測方程為

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

(2)

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

Nijk,abc=

(3)

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

(4)

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

(5)

周跳估值標準差為

(6)

式中σε為三頻載波相位觀測噪聲標準差，周；σe為三頻偽距觀測噪聲標準差， m。

2 SVD形態濾波算法

2.1 構造Hankel矩陣

首先對由測量數據構造的Hankle矩陣進行SVD，對Hankel矩陣取不同維數時對空間劃分和濾波效果的影響進行分析和研究[6]。

假定X=[x1,x2,…,xn]為加入周跳的檢測信號，使用此信號構造Hankel矩陣如下

(7)

式中 1

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

2.2 峭度準則

峭度作為無量綱參數可以描述波形尖峰度同時反映信號的分布特性，其數學描述為

(8)

式中μ為信號X的均值；σ為信號X的標準方差[8]。當信號中出現周跳時，非正常信號的概率密度會增加，當信號幅值分布明顯偏離正態分布時，峭度值會隨之增大。 峭度值越大，說明信號中非正常成分所占比重越多，而周跳信息往往包含在這些非正常成分以及由此引起的幅值調制信號中。

2.3 形態濾波信號選擇

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

假設采集到的一維多值信號f(n)定義域為D[f]={0,1,2,3,…,N} ，選擇結構元素序列為g(x)，且定義域為D[g]={0,1,2,3,…,P}，P和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)

由膨脹腐蝕運算引出的形態學開、閉運算表達式為

f°g=f?g⊕g

(11)

f·g=f⊕g?g

(12)

式中 “° ”為開運算；“·”為閉運算。

(13)

式中N為總采樣點個數；yi0為未疊加噪聲的輸入信號離散采樣點；yi為經過形態學濾波器濾波后的輸出信號離散采樣點。E越小，說明濾波后的信號與未加噪聲的信號越接近[9]。

濾波效果不僅和結構元素的形狀有關還和結構元素的尺寸有關。選取合適的尺寸，對于抑制信號內部的細節差異會起到較好的效果，而且不會使邊界弱化。

3 實驗仿真

實驗采用某公司的雙星五頻北斗數據信號，選取B3頻段37 118 357.406～37 125 029.047之間連續的1 000組歷元。設置數據采樣頻率為1 000 Hz，采樣間隔為1 s，以MATLAB為支持實現實驗仿真。所選數據為原始不含周跳的數據 ，通過在所選數據的不同歷元加入不同周跳，模擬不同類型的周跳探測，檢測本文方法對周跳的探測范圍。通過對單周跳的高精度探測，進行形態濾波實驗仿真與傳統SVD實驗仿真數據對比，以檢驗本文方法的優越性。

通過偽距與載波相位觀測值的二次差分處理，然后再在歷元間作差，得到無鐘差和噪聲干擾的實驗用原始信號，如圖1所示。

圖1 未加周跳的原始信號波形

3.1 模擬不同類型的周跳探測

實驗在不同歷元間加入了不同周跳進行探測，如表1。

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

根據奇異值曲線，如圖2，可以看出m=6時奇異點趨于零，選取m=6構建分析矩陣。

圖2 奇異值曲線

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

圖3 濾波前后分量信號

分量信號的峭度值提取結果分別為4.07，5.33，4.31，3.29，3.24，2.94。選取峭度值較大的分量信號2和分量信號3進行信號重構,結果如圖4。

圖4 重構信號波形

由圖可以看出，在所加5周和3周周跳處可以看出明顯的沖擊信號，所加1周周跳處也存在沖擊信號，但相對不明顯。實驗清晰檢測出了5周和3周周跳，對于1周的周跳信號，雖然也能檢測，但幅值較小，效果不明顯。

3.2 單周跳的探測與對比

通過前述實驗分析，在多個歷元加入周跳時，對1周的小周跳的探測效果并不明顯。為進一步驗證方法的可行性，在400歷元處加入2周周跳進行探測，并采用傳統SVD進行相同探測，通過對探測效果的比較檢驗本方法的優越性。如圖5為形態濾波后的信號波形。

圖5 形態濾波后分量信號波形

根據濾波后的分量信號，可提取峭度值信號，分別為4.01，5.32，5.65，4.19，3.47，10.5。圖6為本文方法對周跳信號的探測結果和傳統SVD方法對周跳的探測結果。

圖6 形態濾波和傳統SVD對周跳信號探測的波形

4 結 論

通過對多個歷元加入不同周跳和在單歷元加入單一周跳與傳統SVD探測方法進行對比實驗，驗證形態濾波結合SVD算法較傳統SVD算法在周跳探測中的優越性，實驗結果表明：將形態濾波用于SVD分量篩選，能夠獲得更清晰的周跳信號，提高了周跳探測的精確性。

[1] 王貴文,王澤民,殷海濤.基于三差觀測量的實時動態GPS周跳修復方法研究[J].武漢大學學報：信息科學版,2007,32(8):711-714.

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

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

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

[5] 張 波,李健君.基于Hankel矩陣與奇異值分解(SVD)的濾波方法以及在飛機顫振試驗數據預處理中的應用[J].振動與沖擊,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變換應用于暫態信號分析[J].電力系統自動化,2005,29(4):49-52.

[8] 胡愛軍,馬萬里,唐貴基.基于集成經驗模態分解和峭度準則的滾動軸承故障特征提取方法[J].中國電機工程學報,2012,32(11):106-111.

[9] 王家良，程春玲.一種多層自適應形態濾波算法[J].計算機科學，2015，42(5)：72-77.