甘 雨 隋立芬 戚國賓 溫勁杰

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

2 93787部隊，北京市，100076

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

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

1 MW 組合探測性能分析

MW 組合觀測值為［4］：

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

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

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

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

假設(shè)兩個頻點的偽距精度相同（均為σP），載波相位誤差很?。珊雎裕?，則MW 組合的精度為：

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

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

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

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

σW初值取為0.5周。若

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 EMD 閾值消噪

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

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

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

閾值估計：

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

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

圖4 端點發(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 組合探測周跳

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

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

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

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

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

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

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

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

3 計算與分析

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

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

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

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

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

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