王 銳, 袁 靜, 班 亞, 楊 帆,羅 浩, 江 力, 羅亞雄

(1.重慶市計量質(zhì)量檢測研究院,重慶 401123; 2.瑞士洛桑聯(lián)邦理工大學(xué), Lausanne, CH 1015，Switzerland)

1 引 言

計算國際原子時的算法是為了保證得到一個可靠,長期穩(wěn)定和準(zhǔn)確的時間尺度。近些年隨著全球各地的高性能原子鐘的引入,使得地方原子時時間尺度以及國際原子時(international atomic time，TAI)的準(zhǔn)確程度越來越高。但是隨著硬件設(shè)備的性能瓶頸，原子時的算法改進則成為提高原子時的主要方向之一。ALGOS是最經(jīng)典的加權(quán)算法,Kalman濾波算法解決了噪聲強弱的影響,AT1(NIST)算法具有實時性的預(yù)測。本文從理論上分析對比了這3種算法。

2 原子時ALGOS算法

原子時的經(jīng)典算法ALGOS是由國際計量局(BIPM)設(shè)計與開發(fā)的一種時間尺度算法，其測量了全球各地實驗室的原子鐘時間尺度,然后進行加權(quán)平均,從而保證優(yōu)化時間尺度的可靠性、頻率的長時間穩(wěn)定性以及頻率的準(zhǔn)確性。

2.1 ALGOS算法基本理論

國際計量局計算原子時(TAI)的算法[1～6]主要有以下兩個步驟：

(1)計算自由原子時的時間尺度(EAL)。利用ALGOS算法,對世界各地超過400臺原子鐘的數(shù)據(jù)進行計算,得出權(quán)重平均數(shù)。

(2)利用基準(zhǔn)頻標(biāo)PFS/SFS校準(zhǔn)EAL的頻率,從而得到原子時(TAI)。

ALGOS算法是經(jīng)典的加權(quán)平均法,定義原子自由時間尺度EAL見式(1)。

xj(t)=[EAL(t)-hj(t)]

(1)

式中:N代表在周期1個月中參與測試的鐘的數(shù)量;wi是鐘Hi對應(yīng)的權(quán)重值,權(quán)重的設(shè)定是為了反映長時間的穩(wěn)定性;hi是鐘Hi在時刻t的測量值,為了保障時間尺度的連續(xù)性,h′i是對鐘Hi的測量值的預(yù)測修正量;xi,j是鐘j和鐘i的相位差。

2.1.1 頻率預(yù)測算法

頻率的預(yù)測是為了在有新的鐘加入或移出時,避免頻率的跳躍和不穩(wěn)定。

具體的方法為: 在兩個連續(xù)的時間段Ik-1=[tk-1,tk]、Ik=[tk,tk+1],為保證頻率和相位的連續(xù)性,需要預(yù)測鐘Hi在tk時刻的修正量h′i。BIMP對遠(yuǎn)程的時間數(shù)據(jù)處理采用的是后置處理模式,獲得各個實驗室的原子鐘的相位差數(shù)據(jù)。在1個月的測試周期內(nèi),每6天產(chǎn)生1個測試數(shù)據(jù)點, 每個數(shù)據(jù)點對應(yīng)的時刻是BIPM對應(yīng)的UTC 00：00:00,故可以得到原子自由時間尺度EAL與每臺鐘測量值的差值xi(tk+nT/6),n=0,…,6;然后對這些數(shù)據(jù)進行最小二乘優(yōu)化擬合得到鐘Hi對于EAL的頻率Bi,Ik(tk+T);從而在1個月的周期Ik=[tk,tk+T], 鐘Hi的修正量h′i為：

h′(t)=ai,Ik(tk)+Bip,Ik(t)(t-tk)

(2)

式中:i是鐘的編號;p是上一段時間預(yù)測的參數(shù);ai,Ik(tk)是在tk時刻,鐘Hi對于EAL的相位偏差。

ai,Ik(tk)的估計值為：

ai,Ik(tk)=EAL(tk)-hi(tk)=xi(tk)

(3)

Bip,Ik(t)的估計值為：

Bip,Ik(t)=Bi(tk),t=tk+nT/6

(4)

2.1.2 權(quán)重算法

(5)

為了避免獲得較大權(quán)重的原子鐘出現(xiàn)不穩(wěn)定情況,最大權(quán)限制理論[7]提出了一種限制：

wi(t)

(6)

式中:wmax為最大權(quán)重比例;A是經(jīng)驗常數(shù)(經(jīng)常取2.5)。選擇合適的A值,達到利用各原子鐘的優(yōu)點和控制擁有滿權(quán)因子鐘的數(shù)量。

2.2 新的ALGOS算法

2.2.1 新的頻率預(yù)測算法

在原子時時間尺度的線性頻率預(yù)測中,假設(shè)原子的頻率不發(fā)生改變，但由于氫鐘和銫鐘都存在頻率漂移,所以線性預(yù)測是不夠精確的。2011年9月,BIPM增加了對頻率漂移量的預(yù)測,提出了新的頻率預(yù)測[8]：

(7)

為了估計式(7)的參數(shù),假定在tk時刻,修正量h′i(t)可以寫成以下形式：

(8)

(9)

從式(9)可以看出：頻率的漂移量是恒定的,而頻率卻不再是個定值。

ai(tk)的估計值：

(10)

Bip,Ik的估計值：

(11)

EAL(tk)-hi(tk)被用來估計相位和頻率,鐘Hi關(guān)于地球時TT的頻率值yTT-hi被用來估計頻率漂移。需要假定在1個月的測試期間內(nèi),鐘的頻率漂移保持不變。氫鐘和銫鐘的漂移量的計算如下：

(1)氫鐘

(12)

(2)銫鐘

當(dāng)采用線性預(yù)測模型時,每個月都需要微調(diào)修正EAL的頻率;而采用新的模型,加入二次項頻率漂移量的模型,EAL幾乎不需要修正。圖1展示了在不同數(shù)據(jù)長度,不同預(yù)報模型條件下,計算得到EAL的穩(wěn)定度;數(shù)據(jù)獲取的時間是2006年和2008年;從圖1可以看出新的二次模型對短期和長期的穩(wěn)定性都有提高。

圖1 不同數(shù)據(jù)長度、不同模型的自由時間尺度穩(wěn)定度Fig.1 The frequency stability of the data with different time scales and models

2.2.2 新的權(quán)重算法

考慮到在鐘的時間尺度系統(tǒng)中,只有鐘的編號和它們的預(yù)測值不同,BIMP在2014年3月對權(quán)重算法提出了新的改進[9],認(rèn)為可預(yù)測的鐘才是一個好鐘。研究發(fā)現(xiàn)如果對頻率預(yù)測準(zhǔn)確,頻率漂移和老化等因素不會對EAL的穩(wěn)定度造成影響。因此,新的權(quán)重方法將對頻率進行預(yù)測,根據(jù)各實驗室的原子鐘的每個月的預(yù)測頻率和真實頻率的差值來預(yù)測權(quán)重因子,從而保證EAL和UTC的長期穩(wěn)定性。

算法的具體步驟如下：

(1)通過相關(guān)的權(quán)重集可以得到EAL-hi的值。在第一次迭代中,采用的是上個時間區(qū)間所歸一化后的權(quán)重；后面的迭代采用前一次迭代的值。

(13)

(3)計算得到步驟(2)中差值的平方。

(4)計算1年以上的εi,Ik來確保EAL的長時間穩(wěn)定性。

(5)在計算權(quán)重因子時,采用了更具有統(tǒng)計意義的濾波方法進行分配；該方法會分配給新加入的原子鐘數(shù)值更大的權(quán)重。

(14)

式中：j表示計算時間間隔；M表示月份,取值范圍是5～12，因UTC的計算至少要5個月的數(shù)據(jù),1年是標(biāo)準(zhǔn)的測試周期。

(6)鐘的權(quán)重因子表達式為：

(15)

BIMP用新的權(quán)重算法對2006～2013年的8年數(shù)據(jù)進行計算，結(jié)果見表1。各原子鐘的最大權(quán)的個數(shù)和占總權(quán)重的比例都發(fā)生了改變,全世界氫鐘最大權(quán)重從10%增加到35%，解決了氫鐘比例過小,銫鐘的權(quán)重比例過大的問題，平衡了原子鐘權(quán)的比例分布。

表1 新舊權(quán)重算法權(quán)重分布對比Tab.1 Statistics relative to the clock distribution in UTC computation according to the chosen algorithm

計算各時間段的Allan方差,對比得到如表2所示的結(jié)果；新的權(quán)重方法Allan方差變小,提高了EAL的短時間和長時間的頻率穩(wěn)定性。

表2 新舊權(quán)重算法Allan方差Tab.2 The values of the overlapping Allan deviation σy(τ )

3 Kalman濾波算法

針對ALGOS算法的權(quán)重大小與某一種噪聲過程的強弱有關(guān)的問題,1982年Barnes 提出了Kalman算法。其核心思想是估值理論而不是權(quán)重概念,主要是對參考鐘和理想時間之差做最佳優(yōu)化,以估計值為修正量,從而計算時間尺度。Kalman算法的模型是將鐘組內(nèi)的所有單鐘模型疊加起來,構(gòu)成鐘組的Kalman濾波器;輸入量是鐘組各鐘與參考鐘的鐘差測量值,輸出量是鐘組內(nèi)所有鐘的相位和頻率的估計值。本節(jié)將詳細(xì)描述Kalman算法的基本原理。[10～12]

設(shè)有N臺鐘,均符合下述模型：

(16)

(17)

式中：xi(t)是第t次測量時,鐘i的測量值與理想?yún)⒖肩姷牟钪?；yi(t)是它們的頻率差值；T是時間間隔；Qi是系統(tǒng)噪聲矩陣；ξi是白色調(diào)頻噪聲,方差是Ei;ηi是隨機游走調(diào)頻噪聲,方差是Hi。

根據(jù)Kalman的濾波原理,鐘組的動態(tài)模型為：

X(t)=φX(t-1)+W(t)

(18)

式(18)中各向量具體表示為：

式中:φ為轉(zhuǎn)移矩陣;W(t)為系統(tǒng)噪聲矩陣。

測量方程為：

Z(t)=HX(t)+V(t)

(19)

式中:Z(t)為測量向量;H為測量矩陣;V(t)是測量噪聲矩陣。

Kalman的濾波方程組的迭代步驟見式(20)：

(20)

式中:P為X的誤差方差矩陣;Q為系統(tǒng)噪聲矩陣;K為Kalman增益;R為測量噪聲矩陣。

科研人員采用18臺HP5071A銫原子鐘的比對數(shù)據(jù),TA(Kalman)和TA(ALGOS)的結(jié)果相似,長期穩(wěn)定度相當(dāng);但是Kalman算法的不完全性會導(dǎo)致估值結(jié)果的誤差隨時間無窮增大。

4 AT1算法

AT1算法與ALGOS算法的核心思想都是加權(quán)算法；AT1不計算過去頻率的真實值,只考慮頻率的變化,是實時的計算。基本原理如下[4,12,14]：

在原子鐘組中選出穩(wěn)定性最好的一臺鐘作為參考鐘,將參考鐘的頻率連接到相位微躍計,通過原子時算法控制微躍器的速率,這樣的一種組合等效為一臺組合鐘,其時間是在每次物理測量完成后計算出來的。

令xi(t),yi(t)分別代表第i臺鐘相對于組合鐘在時刻t的時間差和速率,設(shè)τ為測量中的時間間隔,xi(t)、yi(t)都需要預(yù)先設(shè)定初始值。

這臺鐘在t+τ時刻相對于組合鐘的時間差可以估算為：

(21)

式中:yi(t)·τ表示從t時刻到t+τ時刻相對于組合鐘的時間變化量。

參考鐘相對于組合鐘的時間和速率也遵循式(21)，區(qū)別只是參考鐘用角標(biāo)r。硬件在時刻t+τ測量得到的第i臺鐘和參考鐘的時間差,稱之為ti(t+τ)。 參考鐘相對于組合鐘的鐘差就可以通過第i臺鐘估算得到。

(22)

式(22)表示用第i臺鐘在t時刻的時間差和速率計算出組合鐘相對于參考鐘時間差的估算值。

通過之前得到的時間和速率數(shù)據(jù),再和現(xiàn)在的測量結(jié)果結(jié)合,每臺鐘都可以提供參考鐘相對組合鐘的時差估算。如果鐘有N臺,可重復(fù)N-1次獨立估算,組合鐘相對于參考鐘的時間可由加權(quán)平均得出。

(23)

5 各算法時間尺度分析與對比

理想的原子鐘時間尺度性質(zhì)是具有長期穩(wěn)定性、實時性和準(zhǔn)確性，因此,很多算法都使用加權(quán)平均組合鐘的時間尺度,這就等于增強了時間尺度的穩(wěn)定性。原子鐘時間尺度的噪聲也是由各原子鐘噪聲的線性疊加,這些噪聲可以用一些現(xiàn)有的數(shù)學(xué)模型去可靠地估計。根據(jù)前面章節(jié)對ALGOS算法、Kalman濾波算法和AT1算法的思路分析,可以得到以下結(jié)論：

(1)核心思想

ALGOS算法和AT1算法的核心思想是加權(quán)平均法,都分配給每臺原子鐘一個權(quán)重，從而增加時間尺度的穩(wěn)定性；但是沒有消除或者抑制原子鐘的噪聲,不能進一步提高原子時間尺度的精度和可靠性。Kalman算法經(jīng)過一系列的濾波算法,對于多原子鐘上各種噪聲進行數(shù)學(xué)建模從而達到部分的抑制。

(2)穩(wěn)定性

Kalman濾波算法測得的時間尺度具有長期和短期的穩(wěn)定性, ALGOS算法主要表明原子鐘具有長期穩(wěn)定性；但Kalman濾波算法的不完全性會導(dǎo)致估值結(jié)果的誤差隨時間無窮增大,在應(yīng)用Kalman算法時,必須考慮Kalman算法的發(fā)散性問題,合理選擇估值區(qū)間

(3)實時性與應(yīng)用范圍

AT1算法和Kalman濾波算法都具有實時性,而ALGOS是一種滯后算法,因為其測量時間為1個月；帶有滯后的算法不適用于很多場景,如衛(wèi)星導(dǎo)航系統(tǒng)需要實時地預(yù)測星載原子鐘的狀態(tài),并對其做出修正。若滯后1個月來進行修正,衛(wèi)星鐘所積累的誤差對于衛(wèi)星高精度的測距和估時是毀滅性的。因此，AT1算法和Kalman濾波算法適應(yīng)性相對更廣。

6 結(jié) 語

本文詳細(xì)地分析了目前常用的3種基本的原子鐘時間尺度的計算方法,從理論上分析了ALGOS算法、Kalman濾波算法和AT1算法數(shù)學(xué)模型；對比了這3種算法的優(yōu)缺點,以及其使用范圍。結(jié)果表明：ALGOS,AT1是加權(quán)平均算法,Kalman是濾波算法；AT1和Kalman具有實時性,并且應(yīng)用范圍更廣;Kalman濾波算法能更好地處理噪聲的影響,適用于衛(wèi)星導(dǎo)航系統(tǒng)。