利用EMD閾值消噪改進(jìn)MW 組合周跳探測(cè)性能

甘 雨 隋立芬 戚國賓 溫勁杰

1 信息工程大學(xué)地理空間信息學(xué)院，鄭州市科學(xué)大道62號(hào)，450001

2 93787部隊(duì)，北京市，100076

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

針對(duì)MW 組合噪聲水平較高的情況，使用經(jīng)驗(yàn)?zāi)７纸猓╡mpirical mode decomposition，EMD）閾值消噪的方法進(jìn)行處理。對(duì)于EMD 的端點(diǎn)效應(yīng)問題，提出在MW 遞推均值上增加虛擬噪聲構(gòu)成虛擬觀測(cè)值的方法對(duì)數(shù)據(jù)進(jìn)行延拓。對(duì)MW組合進(jìn)行分解，估計(jì)各層分量的閾值，削弱噪聲影響。分析消噪前后含周跳MW 組合所發(fā)生的變化，指出周跳具有“擴(kuò)散”現(xiàn)象，提出基于消噪后MW 組合的探測(cè)方案。

1 MW 組合探測(cè)性能分析

MW 組合觀測(cè)值為［4］：

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

寬巷波長λW表達(dá)式為：

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

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

假設(shè)兩個(gè)頻點(diǎn)的偽距精度相同（均為σP），載波相位誤差很小（可忽略），則MW 組合的精度為：

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

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

才能滿足σW＜1／4＝0.25的需求。受到高度角和多路徑效應(yīng)影響，式（6）的條件有時(shí)難以滿足，在動(dòng)態(tài)條件下尤為嚴(yán)重。

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

σW初值取為0.5周。若

則認(rèn)為k歷元發(fā)生周跳。

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

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

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

圖2 ＋1周周跳時(shí)的漏檢區(qū)域Fig.2 Miss detection area（＋1cycle slip）

圖3 －1周周跳時(shí)的漏檢區(qū)域Fig.3 Miss detection area（－1cycle slip）

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

2 EMD 閾值消噪MW 組合的周跳探測(cè)

2.1 EMD原理及端點(diǎn)效應(yīng)的抑制

經(jīng)驗(yàn)?zāi)７纸猓?0］能把復(fù)雜的信號(hào)分解為一系列固有模態(tài)函數(shù)（IMF，intrinsic mode function）分量和余項(xiàng)之和的形式。分解出的各分量頻率由高到低，每個(gè)分量為一零均值平穩(wěn)信號(hào)。分解的具體過程為：

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

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

3）將r1作為新原始信號(hào)重復(fù)上述過程，逐個(gè)提取出n個(gè)IMF分量imf2，imf3，…，imfn，得到余項(xiàng)rn＝rn－1－imfn。當(dāng)imfn或rn小于預(yù)先設(shè)定的值，或rn已經(jīng)成為單調(diào)函數(shù)時(shí)，分解完成。信號(hào)x（t）分解成n個(gè)IMF分量和1個(gè)余項(xiàng)的和：

在分解過程中，求包絡(luò)線是通過極值點(diǎn)進(jìn)行插值擬合得到的，在樣條插值時(shí)，如果數(shù)據(jù)兩個(gè)端點(diǎn)本身不是極值點(diǎn)，就不能確定端點(diǎn)處的極值點(diǎn)，導(dǎo)致構(gòu)成包絡(luò)線的三次樣條曲線在數(shù)據(jù)序列的兩端出現(xiàn)發(fā)散現(xiàn)象，并且這種發(fā)散的結(jié)果可能會(huì)逐漸向內(nèi)“污染”，這就是EMD 端點(diǎn)效應(yīng)［11］。

如果在MW 組合消噪過程中發(fā)生周跳，且周跳位置發(fā)生在數(shù)據(jù)段的邊緣，端點(diǎn)效應(yīng)會(huì)對(duì)含周跳的組合觀測(cè)值進(jìn)行錯(cuò)誤的分解，導(dǎo)致周跳被湮沒。如果能對(duì)數(shù)據(jù)進(jìn)行延拓，即在端點(diǎn)以外再增加一段虛擬觀測(cè)數(shù)據(jù)，則原來發(fā)生周跳的數(shù)據(jù)處于數(shù)據(jù)中部，受端點(diǎn)效應(yīng)的影響將大大降低，關(guān)鍵在于所增加的數(shù)據(jù)要和已有數(shù)據(jù)具有類似的特性。

在第k個(gè)歷元，我們由之前k－1個(gè)歷元的數(shù)據(jù)經(jīng)遞推平均能夠得到數(shù)據(jù)的均值和標(biāo)準(zhǔn)差σW（k－1），在假定MW 組合觀測(cè)符合正態(tài)分布的情況下，可以產(chǎn)生n個(gè)虛擬數(shù)據(jù)：

其中，randn產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)組。此為向后延拓的方向，向前延拓方法類似。

2.2 EMD 閾值消噪

文獻(xiàn)［13－14］提出了基于Hurst參數(shù)的閾值消噪方法。對(duì)前兩個(gè)IMF分量的標(biāo)準(zhǔn)差采用具有抗差性的中位數(shù)來估計(jì)，以免受到周跳成分的干擾：

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

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

閾值估計(jì)：

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

MW 組合觀測(cè)值如圖4 所示。對(duì)衛(wèi)星PRN10的一段數(shù)據(jù)最后一個(gè)歷元的MW 組合加入＋1周周跳，分別用原始數(shù)據(jù)和數(shù)據(jù)延拓后的數(shù)據(jù)進(jìn)行EMD 閾值消噪，結(jié)果見圖5、6。可以看出，如果不對(duì)數(shù)據(jù)進(jìn)行延拓，受端點(diǎn)效應(yīng)的影響，周跳會(huì)被湮沒。

圖4 端點(diǎn)發(fā)生周跳的MW 組合Fig.4 MW with cycle slip on the end

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

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

2.3 消噪MW 組合探測(cè)周跳

MW 組合中的觀測(cè)值若發(fā)生周跳，相當(dāng)于數(shù)據(jù)發(fā)生極強(qiáng)的突變，EMD 閾值消噪將這種突變進(jìn)行平滑，讓整個(gè)數(shù)據(jù)的變化更為平緩。因此，某個(gè)歷元上發(fā)生的周跳經(jīng)過閾值消噪后，會(huì)“擴(kuò)散”到附近的歷元，如圖6。圖7、8是另一段MW 組合觀測(cè)值及其消噪的結(jié)果，周跳未發(fā)生在端點(diǎn)處而是發(fā)生在中間，周跳向前后兩個(gè)方向擴(kuò)散。

消噪后周跳的“擴(kuò)散”現(xiàn)象對(duì)于周跳的發(fā)現(xiàn)是有利的。假設(shè)在k歷元發(fā)生周跳，則對(duì)k－1歷元及之前歷元的數(shù)據(jù)，采用EMD 閾值消噪后MW組合未發(fā)生突變，即k－1、k－2…等歷元均無法探測(cè)到周跳；而使用包含k歷元的數(shù)據(jù)進(jìn)行消噪后，組合觀測(cè)值在k、k－1、k－2…等歷元均可能超過檢驗(yàn)閾值。這說明，使用消噪后的MW 組合探測(cè)周跳時(shí)，不能僅以是否超過檢驗(yàn)閾值作為定位周跳發(fā)生歷元的標(biāo)準(zhǔn)，需結(jié)合前后兩個(gè)歷元的結(jié)果進(jìn)行對(duì)比。

消噪MW 組合進(jìn)行周跳探測(cè)包括以下步驟：

1）基于遞推平均過程估計(jì)的均值和方差，對(duì)MW 組合觀測(cè)值序列進(jìn)行數(shù)據(jù)延拓，如式（11）所示。對(duì)于消噪后的數(shù)據(jù)，可取標(biāo)準(zhǔn)差初值為0.2周或更小。

2）對(duì)信號(hào)進(jìn)行EMD 分解，并進(jìn)行閾值消噪。

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

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

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

3 計(jì)算與分析

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

由上述結(jié)果可以看出：

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

2）對(duì)數(shù)據(jù)進(jìn)行數(shù)據(jù)延拓后再使用EMD，能探測(cè)到位于端點(diǎn)處的周跳。

3）發(fā)生周跳時(shí)，若使用包含周跳對(duì)應(yīng)歷元的數(shù)據(jù)進(jìn)行消噪并探測(cè)，其他歷元也受到周跳的影響，結(jié)合使用不含周跳歷元的數(shù)據(jù)方能最終確定周跳位置。

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

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

4 結(jié) 語

TurboEdit方法無法降低MW 組合觀測(cè)值本身的噪聲水平，在接收機(jī)質(zhì)量不穩(wěn)定、多路徑效應(yīng)顯著等情況下容易造成漏檢。使用EMD 閾值消噪的方法能夠?qū)υ肼曔M(jìn)行降噪處理，根據(jù)噪聲特性得到各層分量的閾值，削弱噪聲影響，提高對(duì)周跳的辨識(shí)度。EMD 存在端點(diǎn)效應(yīng)問題，提出在MW 遞推均值上增加虛擬噪聲構(gòu)成虛擬觀測(cè)值的方法對(duì)數(shù)據(jù)進(jìn)行延拓。進(jìn)行數(shù)據(jù)延拓后，即使在端點(diǎn)處發(fā)生周跳，消噪后周跳特性仍然得到保留。若MW 組合包含周跳，則消噪后觀測(cè)值中的周跳具有“擴(kuò)散”現(xiàn)象。針對(duì)這一現(xiàn)象提出的基于消噪后MW 組合的探測(cè)方案，能夠定位周跳發(fā)生的具體位置，增強(qiáng)消噪后MW 組合探測(cè)周跳的準(zhǔn)確性。

［1］方榮新，施闖，魏娜，等.GPS數(shù)據(jù)質(zhì)量控制中實(shí)時(shí)周跳探測(cè)研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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］陳品馨，章傳銀，黃昆學(xué).用相位減偽距法和電離層殘差法探測(cè)和修復(fù)周跳［J］.大地測(cè)量與地球動(dòng)力學(xué)，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］陽仁貴，歐吉坤，袁運(yùn)斌.一種GPS相位周跳分段平均組合的自動(dòng)修復(fù)方法［J］.大地測(cè)量與地球動(dòng)力學(xué)，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數(shù)據(jù)周跳探測(cè)方法的改進(jìn)［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2011，36（1）：29－34（Wu Jizhong，Shi Chuang，F(xiàn)ang 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］鄭作亞，程宗頤，黃城，等.對(duì)Blewitt周跳探測(cè)與修復(fù)方法的改進(jìn)［J］.天文學(xué)報(bào)，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周跳探測(cè)及修復(fù)的小波變換法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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周跳探測(cè)與修復(fù)中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，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 端點(diǎn)效應(yīng)新方法及其在信號(hào)特征提取中的應(yīng)用［J］.振動(dòng)工程學(xué)報(bào)，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īng)驗(yàn)?zāi)B(tài)分解閾值消噪方法及其在慣性導(dǎo)航系統(tǒng)數(shù)據(jù)處理中的應(yīng)用［J］.測(cè)繪學(xué)報(bào)，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