張東,宋偉偉, 樓益棟,陳亮
(武漢大學 衛(wèi)星導航定位技術研究中心,湖北 武漢 430079)
作為精密單點定位(PPP)的主要誤差源之一,衛(wèi)星鐘差成為近年來全球衛(wèi)星導航系統(tǒng)(GNSS)領域研究的熱點. 目前,實時、高精度、高頻的衛(wèi)星鐘差產(chǎn)品已經(jīng)廣泛應用于實時精密定位、低軌衛(wèi)星定軌以及精密時頻傳遞等領域[1-2]. 傳統(tǒng)的實時衛(wèi)星鐘差估計側(cè)重于精密定位的應用,鐘差產(chǎn)品中的時間基準偏差只需要優(yōu)于10-6s,便可認為不影響定位精度[3]. 衛(wèi)星鐘差產(chǎn)品的時間基準可以通過引入單個測站鐘或者廣播星歷改正項參考的時間基準進行確定[4]. 但是,單個測站鐘引入的時間基準顯然不夠可靠,而廣播星歷自身精度也會影響引入的時間基準. 隨著PPP授時的發(fā)展,傳統(tǒng)方法引入的實時衛(wèi)星鐘差產(chǎn)品的時間基準已不能滿足亞納秒級授時的需求. 因此,需要對實時衛(wèi)星鐘差中的時間基準進行精化.
目前,國際上具有代表性的時間基準建立方法主要有ALGOS時間尺度算法、AT1時間尺度算法以及基于原子鐘噪聲特性的Kalman濾波時間尺度算法[5-6]. ALGOS法是一種事后的時間尺度算法,不適合建立實時的時間基準. 在實時的時間尺度算法中,Kalman濾波算法可以同時建模并抑制多種原子鐘噪聲,能夠克服AT1法無法抑制多種原子鐘噪聲的問題[6]. 因此,本文基于Kalman濾波時間尺度算法,提出一種利用高精度原子鐘精化時間基準的方法,以改善時間基準偏差對實時精密授時的影響. 文章的內(nèi)容主要包括IGS站原子鐘性能評價與測站選擇以及時間基準的改正量處理兩部分. 最后,本文給出修正后的時間基準與IGS時間尺度(IGST)的比較結(jié)果.
時間基準精化的處理方法主要分為:1)采用阿倫方差對觀測網(wǎng)中備選的各原子鐘進行穩(wěn)定度分析,選取部分穩(wěn)定度高的原子鐘用以修正時間基準,在此基礎上, 根據(jù)阿倫方差確定各原子鐘的權(quán)比.然后,利用阿倫方差與原子鐘噪聲的關系,通過選取不同平滑時間間隔的阿倫方差來確定每臺鐘調(diào)相白噪聲(WPM)、調(diào)頻白噪聲(WFM)以及調(diào)頻隨機游走噪聲(RWFM)的過程噪聲參數(shù).2)時間基準的改正量(相對原始時間基準)處理依據(jù)原子鐘變化特性建立數(shù)學模型,并根據(jù)各原子鐘權(quán)比聯(lián)合原子鐘組計算出時間基準的改正量,處理流程如圖1所示.
由于氫原子鐘的短期穩(wěn)定性優(yōu)于銣原子鐘和銫原子鐘,因此本文優(yōu)先采用氫原子鐘修正實時衛(wèi)星鐘差產(chǎn)品中的時間基準. 基準改正量計算的觀測方程可表示為
(1)
IGS跟蹤站配置的原子鐘存在不定期重啟或者間斷工作的現(xiàn)象,這將會導致獲取的原子鐘鐘差數(shù)據(jù)存在較多誤差[7]. 因此,在分析IGS跟蹤站配置原子鐘性能之前,有必要對原子鐘數(shù)據(jù)進行預處理. 一般而言,原子鐘數(shù)據(jù)異常主要表現(xiàn)為相位跳變和頻率跳變. 本文采用基于中位數(shù)(MAD)的探測方法和頻率偏差之差探測法分別對相位跳變和頻率跳變進行探測與修復[8-10].
原子鐘的鐘差數(shù)學模型可以表達為確定性分量、隨機性分量以及觀測噪聲三部分,如下[11]:
(t-t0)2+φ(t),
(2)
式中:x(t0)為初始時差;y(t0)為初始頻差;d為線性頻漂;φ(t)則對應原子鐘的隨機性分量和觀測噪聲分量,φ(t)可以用五種獨立的能量譜噪聲來描述,分別為RWFM、調(diào)頻閃變噪聲(FFM)、WFM、調(diào)相閃變噪聲(FPM)和WPM. 對于一定的平滑時間而言,通常只有一種或兩種噪聲起主導作用.
由于FPM等有色噪聲的影響,原子鐘頻率的隨機過程呈現(xiàn)非平穩(wěn)的特性,這使得φ(t)的標準差會隨著采樣個數(shù)的增加而發(fā)散. 為了解決這個問題,阿倫提出了使用阿倫方差表征原子鐘時域頻率穩(wěn)定度. 目前已經(jīng)證明,對常見的能量譜噪聲(-2≤α≤2),阿倫方差是收斂的[12]. 基于頻率偏差數(shù)據(jù){yn,n=1,2,…,M}的阿倫方差可表示為
(3)
(4)
文獻[13-14]給出了阿倫方差與各種噪聲參數(shù)的關系:
(5)
式中:q0、q1、q2和q3分別代表WPM、調(diào)相隨機游走噪聲、RWFM和調(diào)頻隨機奔跑噪聲對應的過程噪聲參數(shù). 通過式(5)可知,在雙對數(shù)阿倫方差圖中,q0、q1、q2和q3對應的斜率分別為-1、-1/2、1/2和3/2. 因此,阿倫方差可識別在某段平滑時間內(nèi)起主導作用的過程噪聲. 對于氫鐘和銫鐘而言,一般不考慮甚低頻噪聲q3的影響. 通過式(3)可以獲得不同時間間隔的阿倫方差,進而可以采用最小二乘的方法來獲取過程噪聲系數(shù)q0、q1、q2和q3.
根據(jù)(1)式,認為未出現(xiàn)頻率跳變的原子鐘在同一收斂弧段內(nèi)鐘速均相同,將相鄰歷元的觀測值作差,網(wǎng)解處理的觀測方程可以整理為
(6)
式中:τ為觀測數(shù)據(jù)的時間間隔,T(ti)-T(ti-1)為基準改正量的歷元間變化量. 由于每引入一個測站的觀測方程,必然會引入該站的鐘速和基準改正量的歷元間變化量. 因此觀測方程秩虧,此時需引入一個測站的鐘速作為約束才能求解.
(7)
則聯(lián)立式(6)與式(7)可組成新的觀測方程,求解出鐘速向量,將估計出的鐘速向量記為
(8)
由于氫原子鐘的鐘速容易受到WFM以及隨機游走噪聲的影響,鐘速估計的結(jié)果需要使用Kalman濾波進一步處理. Kalman濾波的狀態(tài)方程和觀測方程可分別表示為
(9)
式中:估計量α(ti)為各原子鐘鐘速;向量e(ti)代表的是估計量的RWFM過程;向量n(ti)為鐘速的觀測白噪聲,即WFM. 相應的過程噪聲矩陣Q和觀測白噪聲可分別表示為
(10)
(11)
式中:pn表示相應的原子鐘對應的權(quán)重;ΔT(0)為一常量,不會影響基準改正后的穩(wěn)定度,因此可令ΔT(0)=0.
本文選取的實驗數(shù)據(jù)為基于無電離層組合觀測值采用歷元間差分方法估計的實時衛(wèi)星鐘差產(chǎn)品(2017年4月27日-2017年5月2日). 在鐘差估計中,觀測數(shù)據(jù)由IGS實時數(shù)據(jù)流提供,測站坐標從IGS單天解的SNX文件中獲得,衛(wèi)星軌道由 IGS提供的超快速星歷確定,鐘差產(chǎn)品的原始時間基準通過廣播星歷改正項引入. 但廣播星歷的鐘差預報精度有限,引入的時間基準將不可避免地存在誤差. 目前,該策略引入的鐘差基準精度為亞納秒級,萬秒穩(wěn)量級為10-14. 本文依據(jù)鐘差基準精化方法,結(jié)合IGS站原子鐘數(shù)據(jù)對鐘差基準進行改正,并將修正后的鐘差基準與IGS精密鐘差產(chǎn)品定義的時間基準進行比對.
目前,IGS跟蹤站中具有約150個外接原子鐘,原子鐘類型包括氫鐘、銫鐘以及銣鐘3種,本文優(yōu)先選取數(shù)據(jù)質(zhì)量完整的64臺氫鐘進行穩(wěn)定度分析. 首先,采用各IGS站2017年1月1日-2017年6月30日觀測數(shù)據(jù),使用IGS最終精密星歷以及鐘差產(chǎn)品進行事后PPP解算,得到相對于IGST基準的各測站接收機鐘差,并對各站鐘差長期時間序列進行阿倫方差分析. 需要指出的是,由于IGS最終精密星歷和鐘差產(chǎn)品以天為單位給出,并且IGST基準在天與天間可能存在跳躍的情況[7]. 因此,本文同樣采用阿倫方差對單天鐘差時間序列進行分析. 依據(jù)各站2 min平滑時間的阿倫方差結(jié)果,本文初步選取45臺短期穩(wěn)定度較高的測站原子鐘,2 min平滑時間的阿倫方差結(jié)果如圖2所示.
圖2 IGS跟蹤站氫原子鐘阿倫方差對比圖
另外,選取原子鐘時同樣需要考慮測站的分布位置與觀測質(zhì)量,上述45臺氫鐘測站分布如圖3所示,可以看到大部分分布在北美、歐洲地區(qū),而亞太及非洲區(qū)域分布較少. 綜合考慮測站分布的均勻性,本文最終選取了全球分布的27個測站原子鐘用以修正時間基準,如圖3中五角星標記所示.
圖3 部分IGS跟蹤站外接氫原子鐘分布圖
選取的27臺氫原子鐘的阿倫方差雙對數(shù)圖如圖4所示.
圖4 IGS跟蹤站原子鐘阿倫方差雙對數(shù)圖
通過阿倫方差雙對數(shù)圖的斜率可以看出,影響這些原子鐘穩(wěn)定度的基本為WPM、調(diào)相隨機游走噪聲、FFM以及RWFM. 但明顯可以看出每臺原子鐘對應的雙對數(shù)曲線走勢并不完全一致,因此需要對各臺原子鐘相應的過程噪聲系數(shù)分別估計. 現(xiàn)分別取不同平滑時間對應的阿倫方差στ,再根據(jù)式(5)使用最小二乘對WPM、調(diào)相隨機游走噪聲、RWFM對應的噪聲系數(shù)q0、q1和q2進行多項式擬合. 需要指出的是,如果同時擬合三種隨機噪聲系數(shù),方程的系數(shù)陣容易病態(tài),不利于各項噪聲系數(shù)的估計. 由于對于一定平滑時間的阿倫方差,通常只有一種噪聲起主導作用. 因此,本文采用 “先識別、再估計”的方法對某一時段內(nèi)起主導作用的噪聲系數(shù)進行擬合[15]. 根據(jù)各站的阿倫方差通過式(4)可以給出修正時間基準的原子鐘權(quán)比. 具體結(jié)果如表1所示.
表1 IGS跟蹤站噪聲參數(shù)及權(quán)比列表
表1(續(xù))
根據(jù)記錄下來的實時衛(wèi)星鐘差產(chǎn)品,再結(jié)合選取的27個IGS站觀測數(shù)據(jù)進行PPP解算,獲取接收機鐘差,即上文對應的Lj(ti),各站的原始鐘差序列如圖5所示.
圖5 IGS跟蹤站原始鐘差序列
如圖5所示:不同測站的原子鐘相對于衛(wèi)星鐘差基準變化具有明顯的線性趨勢項,即不同的鐘速. 當各站原始鐘差序列扣除線性趨勢項后,得到的接收機鐘差時間序列如圖6所示.
圖6 去線性趨勢項后各測站的鐘差時間序列
由圖6可知:去除線性趨勢項后各站鐘差序列存在較為一致的變化趨勢,這種一致性說明了接收機鐘差序列的變化包含了基準差異部分.此外,CEDU、TID1、YAR2以及NIST四個站存在鐘速重收斂的現(xiàn)象. 這是因為CEDU、TID1、YAR2以及NIST四個站原始的接收機鐘差序列存在頻率跳變的現(xiàn)象(如圖5所示). 由于濾波中的狀態(tài)方程認為鐘速短時間內(nèi)保持不變,因此接收機鐘的頻率跳變會導致短時間內(nèi)鐘速序列的異常. 當測站原始鐘差恢復穩(wěn)定時,剔除鐘速后的鐘差序列重新收斂,并與其余鐘差序列的變化保持一致,可重新用來精化時間基準. 考慮到測站鐘速跳變的現(xiàn)象,本文采用中位數(shù)探測的方法對各站得出的基準改正量進行探測,異常站基準改正量將被中位數(shù)代替. 各站得出的基準改正量按照相應的權(quán)比加權(quán)平均,便可得到最終的鐘差基準改正量如圖7所示.
圖7 時間基準改正量
IGS為了提高鐘差產(chǎn)品的穩(wěn)定性,建立了事后時間尺度IGST,其短期穩(wěn)定度明顯優(yōu)于廣播星歷參考的時間基準GPST[16]. 因此,可以將IGST作為參考,用來評估鐘差產(chǎn)品時間基準的精度以及穩(wěn)定度. 本文采用YAO提出的以整體衛(wèi)星星座基準作為參考基準的鐘差評估方法[17],將基準改正前后的實時鐘差分別與IGS(30 s采樣間隔)最終精密鐘差比較,得出改正前后的基準相對于IGST的偏差如圖8所示.
圖8 基準改正前后與IGST比較的結(jié)果
從圖8可以看出,若以IGS最終精密衛(wèi)星鐘差產(chǎn)品定義的時間基準(IGST)作為參考基準,可以觀察到原始鐘差基準中一些異常抖動得以消除. 同時,精化后的鐘差基準明顯更為平滑,這可以說明濾波以及原子鐘的加權(quán)平均等方式能夠有效抑制鐘差基準的噪聲. 另外,鐘差基準改正前后的STD以及萬秒穩(wěn)如表2所示.
表2 基準改正前后精度統(tǒng)計表
經(jīng)過修正后的鐘差基準單天內(nèi)的STD值可以達到0.1 ns以內(nèi),精度相比于基準改正前提高了37%至93%不等. 基準改正后的穩(wěn)定度較之前也有顯著提高,天內(nèi)萬秒穩(wěn)可以達到10-15量級,實現(xiàn)了一個量級的提高.
在PPP的定位精度分析中,衛(wèi)星鐘差可認為由基準偏差、初始衛(wèi)星鐘差偏差以及相對鐘差組成[4].由于衛(wèi)星鐘與接收機鐘的強相關性,基準偏差可以被接收機鐘差吸收從而不影響定位精度.而初始衛(wèi)星鐘差可以被模糊度參數(shù)吸收,同樣不會影響PPP收斂后的定位精度.因此,衛(wèi)星鐘差中的相對鐘差是影響定位精度的關鍵. 在扣除基準偏差和初始衛(wèi)星鐘差偏差后,基準修正前后的實時衛(wèi)星鐘差產(chǎn)品與IGS最終精密產(chǎn)品相對鐘差互差的STD統(tǒng)計結(jié)果如圖9所示.
圖9 衛(wèi)星鐘差相對鐘差精度對比圖
由圖9可知,基準修正前后衛(wèi)星鐘差的相對精度保持一致,這表明鐘差基準修正不影響PPP的定位精度.
本文基于阿倫方差對IGS跟蹤站的外接氫原子鐘進行了穩(wěn)定度分析,并綜合原子鐘分布及數(shù)據(jù)質(zhì)量情況選取了精化時間基準的原子鐘組. 同時,基于阿倫方差與原子鐘噪聲的關系,估算了所選取原子鐘的噪聲系數(shù). 在此基礎上,本文依據(jù)原子鐘變化特性建立時間基準的表達模型,實現(xiàn)了鐘差基準的精化. 通過與IGS最終精密鐘差產(chǎn)品連續(xù)6天的比較,精化后的鐘差基準單天內(nèi)STD可以達到0.1 ns以內(nèi),精度最高提升了93%;基準改正后的萬秒穩(wěn)達到10-15量級,實現(xiàn)了一個量級的提高. 此外,通過相對鐘差精度的分析,表明鐘差基準修正不影響PPP的定位精度.