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

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

1 引 言

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

2 原子時(shí)ALGOS算法

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

2.1 ALGOS算法基本理論

國(guó)際計(jì)量局計(jì)算原子時(shí)(TAI)的算法[1～6]主要有以下兩個(gè)步驟：

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

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

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

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

(1)

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

2.1.1 頻率預(yù)測(cè)算法

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

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

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

(2)

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

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

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

(3)

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

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)驗(yàn)常數(shù)(經(jīng)常取2.5)。選擇合適的A值,達(dá)到利用各原子鐘的優(yōu)點(diǎn)和控制擁有滿(mǎn)權(quán)因子鐘的數(shù)量。

2.2 新的ALGOS算法

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

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

(7)

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

(8)

(9)

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

ai(tk)的估計(jì)值：

(10)

Bip,Ik的估計(jì)值：

(11)

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

(1)氫鐘

(12)

(2)銫鐘

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

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

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

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

算法的具體步驟如下：

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

(13)

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

(4)計(jì)算1年以上的εi,Ik來(lái)確保EAL的長(zhǎng)時(shí)間穩(wěn)定性。

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

(14)

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

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

(15)

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

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

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

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

3 Kalman濾波算法

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

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

(16)

(17)

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

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

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

(18)

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

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

測(cè)量方程為：

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

(19)

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

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

(20)

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

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

4 AT1算法

AT1算法與ALGOS算法的核心思想都是加權(quán)算法；AT1不計(jì)算過(guò)去頻率的真實(shí)值,只考慮頻率的變化,是實(shí)時(shí)的計(jì)算?；驹砣缦耓4,12,14]：

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

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

這臺(tái)鐘在t+τ時(shí)刻相對(duì)于組合鐘的時(shí)間差可以估算為：

(21)

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

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

(22)

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

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

(23)

5 各算法時(shí)間尺度分析與對(duì)比

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

(1)核心思想

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

(2)穩(wěn)定性

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

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

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

6 結(jié) 語(yǔ)

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