利用EMD閾值消噪改進MW 組合周跳探測性能

甘 雨 隋立芬 戚國賓 溫勁杰

1 信息工程大學地理空間信息學院，鄭州市科學大道62號，450001

2 93787部隊，北京市，100076

GPS 動態非差周跳探測一般需使用多頻觀測值，如無幾何距離組合［1，2］、MW組合［3］等。Blewitt提出的TurboEdit算法［4］使用這兩種組合觀測值聯合探測周跳，目前已得到廣泛應用［5］。TurboEdit算法中，MW 組合通過遞推平均來估計均值和方差，由于估值需逐漸收斂，收斂前出現的周跳無法探測［6］。此外，MW 組合包含偽距觀測值，其噪聲遠大于載波相位，受多路徑效應和衛星高度角等因素的影響，可能無法探測1～2周的小周跳［7］。因為遞推平均過程雖然能減小期望和方差估值受噪聲的影響，但是MW 觀測值本身的噪聲水平并沒有降低。蔡昌盛［7］對本歷元使用后向平滑，在其后一個歷元使用前向平滑，以前后歷元差作為探測檢驗量以減少噪聲影響，但是要求后續歷元用于平滑的數據一定不含周跳，降低了實際應用中的效率和自動化程度。部分學者將小波分析等用于周跳探測，但均針對原始GPS觀測值［8－9］。在動態導航定位中，接收機的某些運動可能引起原始觀測值發生類似于周跳的突變。

針對MW 組合噪聲水平較高的情況，使用經驗模分解（empirical mode decomposition，EMD）閾值消噪的方法進行處理。對于EMD 的端點效應問題，提出在MW 遞推均值上增加虛擬噪聲構成虛擬觀測值的方法對數據進行延拓。對MW組合進行分解，估計各層分量的閾值，削弱噪聲影響。分析消噪前后含周跳MW 組合所發生的變化，指出周跳具有“擴散”現象，提出基于消噪后MW 組合的探測方案。

1 MW 組合探測性能分析

MW 組合觀測值為［4］：

其中，Φ表示載波相位觀測值，P表示偽距，f表示頻率，λ表示波長。該組合消除了幾何距離和大部分誤差項，僅受偽距噪聲和多路徑效應影響。

寬巷波長λW表達式為：

c為光速，波長約為86.2cm，相比原始L1和L2增加數倍。直接以下式進行分析：

若無特別說明，這里MW 組合觀測均指NW。

假設兩個頻點的偽距精度相同（均為σP），載波相位誤差很?。珊雎裕瑒tMW 組合的精度為：

單位為周。代入具體頻率和寬巷波長，得：

若以4倍中誤差作為周跳探測的閾值，要探測1周的周跳，需要

才能滿足σW＜1／4＝0.25的需求。受到高度角和多路徑效應影響，式（6）的條件有時難以滿足，在動態條件下尤為嚴重。

MW 組合遞推平均計算均值和方差為［4］：

σW初值取為0.5周。若

則認為k歷元發生周跳。

可以看出，在周跳探測過程中，NW本身受到噪聲和多路徑效應的影響。若誤差水平較高，有可能發生漏檢。

圖1顯示的是某動態測站在一段時間內觀測衛星PRN10（高度角大于55°）無周跳時的MW組合觀測值（這里實際上是，去除了均值）及其閾值?？梢钥闯?，MW 組合的噪聲水平較高，分布在［－0.664，0.598］范圍內。顯然，在某些歷元內發生1周的小周跳時無法被探測。圖2 中，黑色虛線以下的部分若發生＋1周的周跳，無法探測出來；而圖3中，黑色虛線以上的部分無法探測可能發生的－1周的周跳。

圖1 無周跳時MW 組合及閾值Fig.1 MW and thresholds（0cycle slip）

圖2 ＋1周周跳時的漏檢區域Fig.2 Miss detection area（＋1cycle slip）

圖3 －1周周跳時的漏檢區域Fig.3 Miss detection area（－1cycle slip）

如果能夠降低MW 組合觀測值的噪聲水平，同時合理地估計降噪后觀測的均值和方差，就能探測出小周跳，減小漏檢率。

2 EMD 閾值消噪MW 組合的周跳探測

2.1 EMD原理及端點效應的抑制

經驗模分解［10］能把復雜的信號分解為一系列固有模態函數（IMF，intrinsic mode function）分量和余項之和的形式。分解出的各分量頻率由高到低，每個分量為一零均值平穩信號。分解的具體過程為：

1）根據信號x（t）的局部極大值和極小值，用三次樣條曲線插值求出信號上下包絡線，計算上下包絡的均值m1，得到h1＝x（t）－m1。檢測h1是否滿足IMF的兩個條件：①零點數目與極值點數目相同或最多相差1；②局部極大值點構成的包絡線和局部極小值點構成的包絡線的均值為零。若滿足，則h1為第一個分量imf1。

2）若不滿足，重復執行第1步，直至重復k次后h1k滿足條件，則imf1＝h1k，求出原始信號與imf1的差值：r1＝x（t）－imf1。

3）將r1作為新原始信號重復上述過程，逐個提取出n個IMF分量imf2，imf3，…，imfn，得到余項rn＝rn－1－imfn。當imfn或rn小于預先設定的值，或rn已經成為單調函數時，分解完成。信號x（t）分解成n個IMF分量和1個余項的和：

在分解過程中，求包絡線是通過極值點進行插值擬合得到的，在樣條插值時，如果數據兩個端點本身不是極值點，就不能確定端點處的極值點，導致構成包絡線的三次樣條曲線在數據序列的兩端出現發散現象，并且這種發散的結果可能會逐漸向內“污染”，這就是EMD 端點效應［11］。

如果在MW 組合消噪過程中發生周跳，且周跳位置發生在數據段的邊緣，端點效應會對含周跳的組合觀測值進行錯誤的分解，導致周跳被湮沒。如果能對數據進行延拓，即在端點以外再增加一段虛擬觀測數據，則原來發生周跳的數據處于數據中部，受端點效應的影響將大大降低，關鍵在于所增加的數據要和已有數據具有類似的特性。

在第k個歷元，我們由之前k－1個歷元的數據經遞推平均能夠得到數據的均值和標準差σW（k－1），在假定MW 組合觀測符合正態分布的情況下，可以產生n個虛擬數據：

其中，randn產生標準正態分布的隨機數組。此為向后延拓的方向，向前延拓方法類似。

2.2 EMD 閾值消噪

文獻［13－14］提出了基于Hurst參數的閾值消噪方法。對前兩個IMF分量的標準差采用具有抗差性的中位數來估計，以免受到周跳成分的干擾：

而其他IMF分量的方差按下式遞推計算：

式中，H即為Hurst參數。MW 組合觀測值誤差可視為白噪聲，即H＝0.5，則：

閾值估計：

可以看出，各層分量的閾值不同，EMD 閾值消噪是一種自適應閾值消噪方法。

MW 組合觀測值如圖4 所示。對衛星PRN10的一段數據最后一個歷元的MW 組合加入＋1周周跳，分別用原始數據和數據延拓后的數據進行EMD 閾值消噪，結果見圖5、6。可以看出，如果不對數據進行延拓，受端點效應的影響，周跳會被湮沒。

圖4 端點發生周跳的MW 組合Fig.4 MW with cycle slip on the end

圖5 未使用數據延拓的消噪MWFig.5 De－noised MW without data extension

圖6 使用數據延拓的消噪MWFig.6 De－noised MW with data extension

2.3 消噪MW 組合探測周跳

MW 組合中的觀測值若發生周跳，相當于數據發生極強的突變，EMD 閾值消噪將這種突變進行平滑，讓整個數據的變化更為平緩。因此，某個歷元上發生的周跳經過閾值消噪后，會“擴散”到附近的歷元，如圖6。圖7、8是另一段MW 組合觀測值及其消噪的結果，周跳未發生在端點處而是發生在中間，周跳向前后兩個方向擴散。

消噪后周跳的“擴散”現象對于周跳的發現是有利的。假設在k歷元發生周跳，則對k－1歷元及之前歷元的數據，采用EMD 閾值消噪后MW組合未發生突變，即k－1、k－2…等歷元均無法探測到周跳；而使用包含k歷元的數據進行消噪后，組合觀測值在k、k－1、k－2…等歷元均可能超過檢驗閾值。這說明，使用消噪后的MW 組合探測周跳時，不能僅以是否超過檢驗閾值作為定位周跳發生歷元的標準，需結合前后兩個歷元的結果進行對比。

消噪MW 組合進行周跳探測包括以下步驟：

1）基于遞推平均過程估計的均值和方差，對MW 組合觀測值序列進行數據延拓，如式（11）所示。對于消噪后的數據，可取標準差初值為0.2周或更小。

2）對信號進行EMD 分解，并進行閾值消噪。

3）對于消噪后MW 組合超過閾值的點i、i＋1、i＋2、…，由k＝i及之前一段歷元的數據進行步驟1、2，檢驗消噪后k歷元結果是否超過閾值。若超出，則說明k歷元發生周跳，將k歷元結果用遞推平均計算的均值代替，以免影響后續探測，k往前推進，即k＝i＋1；若未超出，說明當前歷元未發生周跳，直接令k＝i＋1。最終使k遍歷所有超過閾值的點。上述步驟也可由后向前遞推。

圖7 發生周跳的MW 組合Fig.7 Un－de－noised MW with cycle slip

圖8 消噪后發生周跳的MW 組合Fig.8 De－noised MW with cycle slip

3 計算與分析

選擇256個歷元數據進行實驗，未消噪MW組合的誤差與圖1類似。其中，第256個歷元的MW 組合觀測值加入＋1 周的周跳，進行EMD閾值消噪前均進行了數據延拓。圖9顯示的是使用前255個歷元數據探測的結果，圖10為256個歷元數據探測的結果。使用少于255個歷元的結果和圖9類似，這里不再給出。

由上述結果可以看出：

1）經過消噪后，MW 組合觀測值的噪聲水平下降，最大值不超過0.2周，能夠迅速探測出小至1周的周跳。

2）對數據進行數據延拓后再使用EMD，能探測到位于端點處的周跳。

3）發生周跳時，若使用包含周跳對應歷元的數據進行消噪并探測，其他歷元也受到周跳的影響，結合使用不含周跳歷元的數據方能最終確定周跳位置。

圖9 發生周跳的MW 組合Fig.9 Un－de－noised MW with cycle slip

圖10 消噪后發生周跳的MW 組合Fig.10 De－noised MW with cycle slip

4 結 語

TurboEdit方法無法降低MW 組合觀測值本身的噪聲水平，在接收機質量不穩定、多路徑效應顯著等情況下容易造成漏檢。使用EMD 閾值消噪的方法能夠對噪聲進行降噪處理，根據噪聲特性得到各層分量的閾值，削弱噪聲影響，提高對周跳的辨識度。EMD 存在端點效應問題，提出在MW 遞推均值上增加虛擬噪聲構成虛擬觀測值的方法對數據進行延拓。進行數據延拓后，即使在端點處發生周跳，消噪后周跳特性仍然得到保留。若MW 組合包含周跳，則消噪后觀測值中的周跳具有“擴散”現象。針對這一現象提出的基于消噪后MW 組合的探測方案，能夠定位周跳發生的具體位置，增強消噪后MW 組合探測周跳的準確性。

［1］方榮新，施闖，魏娜，等.GPS數據質量控制中實時周跳探測研究［J］.武漢大學學報：信息科學版，34（9）：1 094－1 098（Fang Rongxin，Shi Chuang，Wei Na，et al.Realtime Cycle－Slip Detection for Quality Control of GPS Measurements［J］.Geomatics and Information Science of Wuhan University，34（9）：1 094－1 098）

［2］陳品馨，章傳銀，黃昆學.用相位減偽距法和電離層殘差法探測和修復周跳［J］.大地測量與地球動力學，2010，30（2）：120－125（Chen Pinxin，Zhang Chuanyin，Huang Kunxue.Cycle Slips Detecting and Repairing by Use of Phase Reduce Pseudorange Law and Ionized Layer Remant Method of Difference［J］.Journal of Geodesy and Geodynamics，2010，30（2）：120－125）

［3］陽仁貴，歐吉坤，袁運斌.一種GPS相位周跳分段平均組合的自動修復方法［J］.大地測量與地球動力學，2009，29（5）：76－77（Yang Rengui，Ou Jikun，Yuan Yunbin.An Approach of Sub－Average Joint Computation for Auto Restoring GPS Phase Cycle Slip［J］.Journal of Geodesy and Geodynamics，2009，29（5）：76－77）

［4］Blewitt G.An Automatic Editing Algorithm for GPS Data［J］.Geophysical Research Letters，1990，17（3）：199－202

［5］吳繼忠，施闖，方榮新.TurboEdit單站GPS數據周跳探測方法的改進［J］.武漢大學學報：信息科學版，2011，36（1）：29－34（Wu Jizhong，Shi Chuang，Fang Rongxin.Improving the Single Station Data Cycle Slip Detection Approach TurboEdit［J］.Geomatics and Information Science of Wuhan University，2011，36（1）：29－34）

［6］鄭作亞，程宗頤，黃城，等.對Blewitt周跳探測與修復方法的改進［J］.天文學報，2005，46（2）：216－225（Zheng Zuoya，Cheng Zongyi，Huang Cheng，et al.Improving of Cycleslip Detection and Correction of Blewitt Method［J］.Acta Astronomica Sinica，2005，46（2）：216－225）

［7］Cai Changsheng，Liu Zhizhao，Xia Pengfei，et al.Cycle Slip Detection and Repair for Undifferenced GPS Observations under High Ionospheric Activity［J］.GPS Solution，2013，17：247－260

［8］蔡昌盛，高井祥.GPS周跳探測及修復的小波變換法［J］.武漢大學學報：信息科學版，2007，32（1）：39－43（Cai Changsheng，Gao Jingxiang.Cycle－slip Detection and Correction of GPS Data by Wavelet Transform［J］.Geomatics and Information Science of Wuhan University，2007，32（1）：39－43）

［9］李明然，田林亞，熱比亞.小波和Kalman濾波組合在GPS周跳探測與修復中的應用［J］.大地測量與地球動力學，2012，32（4）：76－79（Li Mingran，Tian Linya，Rebiya.Application of Wavelet and Kalman Filtering Combination in GPS Cycle Slip Detection and Restoration［J］.Journal of Geodesy and Geodynamics，2012，32（4）：76－79）

［10］Huang N E，Shen Z，Long S R，et al.The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis［J］.Proceedings of the Royal Society：A，1998，1971，454：903－993

［11］曹沖鋒，楊世錫，楊將新.一種抑制EMD 端點效應新方法及其在信號特征提取中的應用［J］.振動工程學報，2008，21（6）：588－594（Cao Chongfeng，Yang Shixi，Yang Jiangxin.A New Method for Restraining the End Effect of Empirical Mode Decomposition and its Applications to Signal Feature Extraction［J］.Journal of Vibration Engineering，2008，21（6）：588－594）

［12］Huang N E，Shen S S P.Hilbert－Huang Transform and Its Applications［M］.Singapore：World Scientific，2005

［13］甘雨，隋立芬，王冰.經驗模態分解閾值消噪方法及其在慣性導航系統數據處理中的應用［J］.測繪學報，2012，41（4）：504－510（Gan Yu，Sui Lifen，Wang Bing.EMD Threshold De－Noising and Its Applications in INS Data Processing［J］.Acta Geodaetica et Cartographica Sinica，2012，41（4）：504－510）

［14］Gan Yu，Sui Lifen.An EMDThreshold De－Noising Method for Inertial Sensors［J］.Measurement，2014，49：34－41