雷 雨，趙丹寧

（1.中國科學(xué)院國家授時(shí)中心，陜西 西安 710600；2.中國科學(xué)院時(shí)間頻率基準(zhǔn)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710600；3.中國科學(xué)院精密導(dǎo)航定位與定時(shí)技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710600；4.中國科學(xué)院研究生院，北京 100039）

GPS精密單點(diǎn)定位（Precise Point Positioning，PPP）時(shí)間傳遞就是利用國際GNSS服務(wù)IGS（International GNSS Service）提供的精密星歷與鐘差產(chǎn)品，對單臺測地型時(shí)間傳遞接收機(jī)接收的雙頻載波相位、碼觀測量進(jìn)行解算，獲得本地鐘與IGS時(shí)間尺度（IGST）的差值，兩個(gè)測站再通過一個(gè)簡單的差分運(yùn)算便可實(shí)現(xiàn)納秒至亞納秒級精度的時(shí)間比對［1］。與傳統(tǒng)的GPS CV（Common View）方法相比，PPP的優(yōu)點(diǎn)是測站不受地理位置的限制，它允許任何裝備測地型時(shí)間傳遞接收機(jī)的時(shí)間實(shí)驗(yàn)室參加國際原子時(shí)TAI（International Atomic Time）的計(jì)算。鑒于該方法廣泛的適用性以及高傳遞性能，國際計(jì)量局BIPM（Bureau International des Poids et Mesures）于2009年9月將該技術(shù)用于國際原子時(shí)TAI（International Atomic Time）計(jì)算［2］，并定期發(fā)布于時(shí)間比對鏈路的FTP數(shù)據(jù)庫ftp：//tai.bipm.org/TimeLink/TAIPPP/。

盡管PPP采用了各種精確的誤差改正模型［3］，但解算的鐘差值仍受到未模型化誤差的影響，另外原子頻標(biāo)本身也受到各種噪聲的影響［4］，兩者均表現(xiàn)出隨機(jī)噪聲的特性，因此，對鐘差數(shù)據(jù)進(jìn)行平滑濾波以消除它們的影響是時(shí)間傳遞中一項(xiàng)重要的工作［5］。因?yàn)閮膳_鐘比對結(jié)果的擬合函數(shù)形式難以確定，所以使用曲線擬合或滑動(dòng)平均等方法無法取得良好效果。由于Vondrak平滑方法可以在觀測資料擬合函數(shù)形式未知的情形下，也能有效地進(jìn)行平滑［6］，因此，在GPS CV和TWSTFT（Two Way Satellite Time and Frequency Transfer）中經(jīng)常使用該方法平滑比對結(jié)果［7－8］。Vondrak平滑法的關(guān)鍵在于平滑因子的選?。?］。本文將Vondrak平滑法應(yīng)用于PPP時(shí)間傳遞數(shù)據(jù)的平滑中，利用觀測誤差法選取合理的平滑因子，并對實(shí)測的鐘差比對數(shù)據(jù)進(jìn)行平滑處理，從精確度和穩(wěn)定度兩方面衡量平滑效果。

1 Vondrak平滑方法

1.1 基本原理

設(shè)有觀測資料x（ti）（i＝1～N），Vondrak平滑的基本準(zhǔn)則為

其中：

式中：x′（ti）為待求的平滑值，pi為觀測值x（ti）的權(quán)重。F為Vondrak平滑法的擬合度，S表示對平滑值三階差分平方和的要求，它反映待求平滑曲線總體平滑程度，故稱為平滑度。λ2為待定的正系數(shù)，它在平滑過程中調(diào)整著擬合度與平滑度之間的關(guān)系，在0和∞兩個(gè)極端邊界值之間取值。當(dāng)λ2→0，要使式（1）取最小，必須使F→0，此時(shí)得到的是一條逼近測量數(shù)據(jù)的光滑曲線；相反，當(dāng)λ2→∞時(shí)，若使式（1）達(dá)最小值，必須讓S→0，這時(shí)對應(yīng)的是一條十分光滑的拋物線。由此可知，Vondrak平滑法就是尋求介于對觀測值絕對擬合與絕對平滑之間的一條折中的曲線，折中的程度由平滑因子ε＝1/λ2控制。

1.2 平滑因子選取

Vondrak平滑法的關(guān)鍵在于如何選取合適的平滑因子，它的大小決定了平滑程度的強(qiáng)弱，ε值越小，平滑度越強(qiáng)，相反地，ε值越大，平滑度越弱。目前，常用的選取ε的方法有交叉檢驗(yàn)法［10］、頻率響應(yīng)法和觀測誤差法［11］，本文使用觀測誤差法來選取ε。

首先，用不同的平滑因子對觀測資料x（ti）進(jìn)行Vondrak平滑，得到對應(yīng)的平滑值x′（ti），由式（4）計(jì)算出平滑值的均方誤差

對一系列的σ（ε），選取使σ（ε）≈σm的平滑因子作為最后確定的平滑因子，其中σm代表觀測資料誤差。

然而，該方法在觀測資料精度未知的情況下無法使用。一種辦法是根據(jù)均方誤差σ（ε）隨平滑因子ε的變化規(guī)律來選取對應(yīng)的ε值，即按式（4）計(jì)算出對應(yīng)不同ε的平滑值的σ（ε），以ε為橫坐標(biāo)，σ（ε）為縱坐標(biāo)繪出σ（ε）的曲線圖，觀測σ（ε）曲線的變化規(guī)律，選取σ（ε）變化最緩慢時(shí)所對應(yīng)的ε作為最終選取的平滑因子。

2 實(shí)例分析

2.1 平滑因子選取實(shí)例

采用NTSC（中國科學(xué)院國家授時(shí)中心）與PTB（德國技術(shù)物理研究所）、KRISS（韓國標(biāo)準(zhǔn)科學(xué)研究院）兩條鏈路之間的PPP時(shí)間傳遞結(jié)果，數(shù)據(jù)從BIPM ftp下載，時(shí)間跨度為2011－10－28（MJD 55862）至2011－11－11（MJD 55876）。采用觀測誤差法，通過分析比較選出較為合理的平滑因子。圖1、圖2分別為根據(jù)觀測誤差法對PTB－NTSC與KRISS－NTSC結(jié)果計(jì)算得到的σ（ε）值與ε的對應(yīng)情況。

圖1 根據(jù)PTB－NTSC鏈路得到σ（ε）的變化曲線

圖2 根據(jù)KRISS－NTSC鏈路得到σ（ε）的變化曲線

從圖中可以看出，當(dāng)ε取500～600時(shí)，σ（ε）的變化較為緩慢，因此取ε＝550。

2.2 平滑結(jié)果

為了檢驗(yàn)Vondrak平滑效果，仍選取2.1節(jié)的數(shù)據(jù)做實(shí)驗(yàn)，從精確度和穩(wěn)定度兩個(gè)方面衡量平滑效果。

2.2.1 精確度分析

取平滑因子ε＝550，利用Vondrak平滑法對上述兩條鏈路PPP時(shí)間傳遞系列平滑，原始值和平滑值（為了清楚表示消噪前后的區(qū)別，平滑結(jié)果下移了5ns）的比較如圖3、圖4所示，圖5、圖6給出了局部放大的平滑比較結(jié)果。RMS統(tǒng)計(jì)結(jié)果見表1。

圖3 PTB－NTSC鏈路比對結(jié)果平滑比較

圖4 KRISS－NTSC鏈路比對結(jié)果平滑比較

圖5 PTB－NTSC鏈路比對結(jié)果平滑局部放大圖

圖6 KRISS－NTSC鏈路比對結(jié)果平滑局部放大圖

表1 兩條鏈路平滑前后精確度對照 ns

結(jié)合圖3～6和表1可以看出，對于變化復(fù)雜無規(guī)律的PTB－NTSC、KRISS－NTSC比對序列，平滑濾除了比對結(jié)果中的隨機(jī)噪聲，比對精度有所提高。

2.2.2 穩(wěn)定度分析

采用式（5）表達(dá)的重疊阿倫方差（Overlapping Allan Variance）來估計(jì)PPP時(shí)間傳遞結(jié)果所體現(xiàn)的頻率穩(wěn)定度

式中：xi為第i個(gè)測量結(jié)果（樣本），N為樣本總數(shù)，τ0為樣本采樣間隔，τ為計(jì)算阿倫方差的采樣時(shí)間，τ＝nτ0。

根據(jù)式（5），分別計(jì)算上述兩條鏈路比對結(jié)果平滑前后的穩(wěn)定度，結(jié)果如表2所示。

表2 兩條鏈路平滑前后穩(wěn)定度對照

由表2可見，Vondrak平滑對PTB－NTSC、KRISS－NTSC鏈路1～24h的穩(wěn)定度有較明顯的改善。

3 結(jié) 論

Vondrak平滑法可以在擬合函數(shù)形式未知的情況下對觀測資料進(jìn)行平滑，其平滑程度取決于平滑因子的選取，因此，平滑因子的選取是Vondrak平滑法的關(guān)鍵。利用觀測誤差法，選取了較為合理的平滑因子。

通過Vondrak平滑法對不同時(shí)間實(shí)驗(yàn)室之間的精密單點(diǎn)定位（PPP）時(shí)間傳遞結(jié)果進(jìn)行平滑處理，驗(yàn)證了Vondrak平滑法的有效性，主要得到以下結(jié)論：

1）Vondrak平滑法可有效消除PPP傳遞結(jié)果中的隨機(jī)噪聲，不僅可以提高時(shí)間傳遞的精確度，也能明顯改善PPP傳遞所體現(xiàn)的頻率穩(wěn)定度；

2）利用觀測誤差法選取平滑因子簡單直觀，即使在觀測資料精度未知的情況下也是適用的，但平滑因子的確定還要結(jié)合數(shù)據(jù)處理經(jīng)驗(yàn)和實(shí)際情況綜合予以考慮。

［1］GéRARD PETIT，ZHIHENG JIANG.Precise Point Positioning for TAI Computation［J］.International Journal of Navigation and Observation，2008：1－8.

［2］E F ARIAS，G PANFILO，G PETIT.Timescales at the BIPM［J］.Metrologia，2011，48（4）：145－153.

［3］JAN KOUBA，PIERRE HéROUX.Precise Point Positioning Using IGS Orbit and Clock Products［J］.GPS Solutions，2001，5（2）：12－28.

［4］李變，屈俐俐，高玉平.Kalman算法在原子時(shí)計(jì)算中的應(yīng)用研究［J］.時(shí)間頻率學(xué)報(bào)，2010，33（1）：11－15.

［5］管健安，董立剛.基于Vondark模型平滑精密單點(diǎn)定位（PPP）授時(shí)鐘差解［J］.測繪工程，2012，21（2）：25－28.

［6］VONDRáK J.A contribution to the problem of smoothing observational data［J］.Bulletin of the Astronomical Institute of Czechoslovakia，1969，20（6）：349－355.

［7］李變.利用Vondrak方法處理GPS CV觀測數(shù)據(jù)的隨機(jī)噪聲［J］.計(jì)算機(jī)測量與控制，2006，14（7）：953－954.

［8］Z.JIANG.Interpolation of TW time transfer from measured points onto standard MJD for UTC generation［Z］.France：BIPM TM167，2009

［9］聶士忠.Vondrak數(shù)據(jù)平滑方法及其在微機(jī)上的實(shí)現(xiàn)［J］.石油大學(xué)學(xué)報(bào)：自然科學(xué)版，1994，18（4）：111－114.

［10］朱書博.外測試驗(yàn)中的高采樣數(shù)據(jù)處理方法研究［D］.長沙：國防科學(xué)技術(shù)大學(xué)，2009.

［11］丁月蓉.天文數(shù)據(jù)處理方法［M］.南京：南京大學(xué)出版社，1998：256－257.