基于Kalman濾波的原子時算法研究

孫同川，王振嶺，孫建設(shè)，劉鐵強

(1.中國電子科技集團 第54研究所，石家莊 050051; 2.中國人民解放軍61711部隊，新疆 喀什 844000)

0 引言

時間對一個國家的經(jīng)濟、社會生活、軍事等方面有著十分重要的作用，時間的穩(wěn)定關(guān)系著國家和社會的安全穩(wěn)定，時間頻率體系建設(shè)作為國家的重大基礎(chǔ)設(shè)施，具有非常重要的意義。衛(wèi)星導(dǎo)航、通信互聯(lián)、聯(lián)合作戰(zhàn)、反導(dǎo)攔截、金融、電力等對時間頻率系統(tǒng)的準確度和可靠性的需求越來越高。只有一個高精度的時頻體系，才能提供可靠的時間頻率服務(wù)。隨著當前國產(chǎn)原子鐘技術(shù)水平的突破[1]，為了提高自主可控能力，采用國產(chǎn)原子鐘組進行原子時算法的研究具有重要意義。通常情況下，構(gòu)建時間頻率體系的關(guān)鍵是建立和保持系統(tǒng)時間，時間尺度的性能一方面受時間頻率系統(tǒng)內(nèi)硬件(比如原子鐘)的性能水平的影響，一方面受所采用的時間尺度算法的影響[2]。

較為著名的時候尺度方法主要包括:目前國外普遍使用的ALGOS算法[3]、AT1算法[4]等，20世紀80年代，針對傳統(tǒng)加權(quán)平均算法采用Allan方差計算權(quán)重會忽視主要噪聲過程之外的噪聲過程的問題，美國學(xué)者Barnes將Kalman濾波用于鐘組噪聲的處理提出了Kalman算法[5-6],近年來，國家授時中心對原子鐘數(shù)據(jù)異常情況下如何計算時間尺度進行了研究，并采用改進的Kalman濾波計算時間尺度取得了比較好的結(jié)果[14]，國家授時中心還對Vondrak-Cepek濾波進行研究，并研究分析這種方法在時間尺度計算中的應(yīng)用，這種方法計算出的綜合時間尺度穩(wěn)定度進一步提高[18]。ALGOS算法和AT1算法的核心思想都是加權(quán)平均，具體做法是根據(jù)某個計算周期內(nèi)鐘組內(nèi)各個原子鐘的穩(wěn)定度來確定每臺原子鐘在此計算周期內(nèi)的權(quán)重，通關(guān)比較合理的權(quán)重分配，使鐘組內(nèi)的噪聲最小，以此來提高時間尺度的穩(wěn)定度。但是加權(quán)平均算法只考慮起主要影響的噪聲過程，而不論其他噪聲過程如何，在計算時間尺度時并沒有消除或抑制鐘組內(nèi)的噪聲，相反Kalman算法通過Kalman濾波器對原子鐘噪聲建模消除或抑制鐘組內(nèi)的噪聲。從文獻[7]可以知悉Kalman、AT1算法的短期穩(wěn)定度較好，而ALGOS的長期穩(wěn)定度較好，AT1算法和Kalman算法是根據(jù)當前計算周期的鐘差數(shù)據(jù)來預(yù)測下個計算周期的鐘差，具有實時性的特點，而ALGOS算法則是采用過去一段時間的鐘差數(shù)據(jù)來計算時間尺度[8]，具有滯后性，造成滯后的原因主要是因為ALOGS算法的權(quán)重計算和速率預(yù)報需要采用數(shù)據(jù)計算前一個月的鐘差預(yù)測值。Kalman算法的穩(wěn)定度較好，但是存在發(fā)散性的問題[8]。

守時實驗室的目標是產(chǎn)生一個穩(wěn)定、準確、可靠的時間尺度，在二級守時節(jié)點中，守時鐘組類型單一(大多為銫鐘組)，針對此前我國大多數(shù)實驗室采用的原子鐘大多為進口原子鐘的問題，為了提高自主可控水平，經(jīng)典的ALGOS算法主要考慮原子時的長期穩(wěn)定度，在二級守時節(jié)點中綜合原子時應(yīng)同時注重綜合原子時的短期穩(wěn)定度和靈活性，以便守時系統(tǒng)隨時進行調(diào)整。

因此對加權(quán)平均算法進行了研究和改進，改進后的算法將針對傳統(tǒng)加權(quán)平均算法中的噪聲問題，加入Kalman濾波過程用于改善頻率預(yù)測值并加入了頻率跳變檢測，可以進一步抑制原子鐘的噪聲，提高綜合原子時的穩(wěn)定度。首先研究分析加權(quán)平均算法的基本原理，然后對算法進行改進，最后通過搭建的試驗平臺對改進的算法性能進行驗證分析。

1 原子鐘噪聲模型的研究

時間尺度算法目的就是研究原子鐘間的噪聲問題，通過設(shè)定符合實際鐘組情況的參數(shù)使得鐘組內(nèi)的噪聲最小，以此提高時間尺度的穩(wěn)定度。而原子鐘作為一種表征頻率的設(shè)備，其紙面讀數(shù)與絕對時間之間必然存在一定的偏差，同時原子鐘的紙面讀數(shù)還受到觀測噪聲的影響，因此為了分析原子鐘的噪聲，首先要對其建立合適的噪聲模型。

原子鐘的觀測方程可以表示為[20]：

Z(tk)=X1(tk)+σ·ξ(k)

(1)

式中，Z(tk)為觀測量；σ·ξ(k)為觀測噪聲，ξ是標準高斯白噪聲，屬于WPM，σ表示噪聲強度。

原子鐘模型可用隨機微分方程[8-9]描述：

(2)

式中,X1(t)，X2(t)表示原子鐘的兩個狀態(tài)變量，X1表示時差狀態(tài)變量，X2表示頻率偏差的隨機游走噪聲，W1(t)，W2(t)是兩個獨立的維納過程。為了方便分析仿真實驗時將此微分方程的初值設(shè)置為0。J.A.Barnes等人提出原子鐘噪聲的冪律譜模型即原子鐘噪聲是五種噪聲分量線性疊加的結(jié)果，文獻[11]不僅描述了各個噪聲特性還給出了仿真方法，它的譜密度函數(shù)：

(3)

不同的α取值對應(yīng)不同類型的噪聲，具體為[17]：α=-2是隨機游走調(diào)頻噪聲(RWFM)，α=-1是閃變調(diào)頻噪聲(FFM)，α=0是調(diào)頻白噪聲(WFM)，α=1是閃變調(diào)相噪聲(FPM)，α=2是調(diào)相白噪聲(WPM)，不同原子時之間的算法受到各種噪聲不同程度的干擾，但這些噪聲有時不能避免，有時可以采用加入濾波器的方式進行消除。AT1算法中一般都采用了指數(shù)濾波器，指數(shù)濾波器可以比較有效的控制原子鐘的調(diào)頻白噪聲和頻率隨機游走噪聲。本文算法采用的濾波器是卡爾曼形式的，卡爾曼濾波器是最小均方意義下的最優(yōu)估計，可以有效避免觀測量中頻差隨機游走噪聲的影響。

加權(quán)平均算法在計算過程中用Allan方差表征一臺原子鐘的穩(wěn)定度，Allan方差主要與原子鐘的噪聲有關(guān)，每臺鐘的權(quán)重大小主要與某個起主要作用的噪聲過程有關(guān)，針對加權(quán)平均算法中的權(quán)重和噪聲過程問題，Kalman算法從估值的角度出發(fā)，避免了權(quán)重的分配問題，對所選取的主鐘和標準時間的差做最小均方差意義下的最優(yōu)估計。這種方法通過對每臺原子鐘建立一個狀態(tài)方程評估的噪聲模型，來構(gòu)成Kalman濾波器。

假設(shè)原子鐘的噪聲模型符合：

(4)

(5)

式中，xi(t)、yi(t)為第t次測量時的時間和頻率與理想時間尺度的偏差；Qi是系統(tǒng)噪聲矩陣；ηi為隨機游走調(diào)頻噪聲，其方差是Hi；ξi為白色調(diào)頻噪聲，其方差是Ei。

由Kalman濾波原理，得到鐘組的動態(tài)模型為：

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

(6)

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

鐘組的觀測方程為：

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

(7)

(8)

式中,H是測量矩陣，V(t)為測量噪聲

(9)

根據(jù)卡爾曼濾波模型和噪聲矩陣，通過卡爾曼濾波迭代計算得到每臺鐘的狀態(tài)，具體迭代步驟如下：

(10)

式中,P為X的誤差方差矩陣；Q為系統(tǒng)驅(qū)動噪聲矩陣；K為卡爾曼增益，R為測量噪聲矩陣。

2 綜合原子時算法研究與改進

目前守時實驗室常用的綜合原子時算法主要有兩大類：加權(quán)平均算法和各種濾波類算法[15]。時間尺度產(chǎn)生最典型的例子就是世界協(xié)調(diào)時(UTC)的產(chǎn)生，國際計量局(BIPM)采用ALGOS算法對全球的守時實驗室數(shù)據(jù)加權(quán)平均，以此提高時間尺度的可靠性和穩(wěn)定性。上文提到加權(quán)平均算法的基本原理都是通過給每臺原子鐘分配一個權(quán)重來調(diào)整鐘組內(nèi)的噪聲關(guān)系[16]，以此來提高時間尺度的穩(wěn)定性。但ALGOS算法采用連續(xù)一個月的數(shù)據(jù)來預(yù)測前一個月的鐘差值，具有滯后性參數(shù)調(diào)整設(shè)置極其不方便，需要的鐘差數(shù)據(jù)時間太長， 不適合在規(guī)模較小的守時實驗室使用。

理想狀態(tài)下，時間尺度算法的最主要步驟是根據(jù)鐘組內(nèi)N臺原子鐘，利用N-1組鐘差對鐘組內(nèi)的原子鐘進行權(quán)重調(diào)整，以及鐘差預(yù)報等。而時間尺度算法就是通過計算每臺鐘的權(quán)重來調(diào)整鐘組內(nèi)的噪聲關(guān)系使其對鐘組產(chǎn)生的時間尺度影響最小化，由此生成的時間尺度比鐘組中任何一個單獨的原子鐘具有更高的可靠性、穩(wěn)定性和頻率精度。鐘組產(chǎn)生的時間尺度是采用數(shù)學(xué)方法計算得到的，并不是某臺鐘的時間也不是他們簡單的集合，這就是綜合原子時，它的物理實現(xiàn)是通過對鐘組中任意一個原子鐘進行適當?shù)男Ｕ齺韺崿F(xiàn)的，如果沒有測量噪聲的影響，該值與所采用的原子鐘無關(guān)。時間尺度算法的輸入是每個原子鐘與選取的主鐘之間的鐘差，原子鐘作為一個物理系統(tǒng)產(chǎn)生一個頻率，原子鐘的時間是由頻率經(jīng)過人工推導(dǎo)得來的，頻率才是真正的物理量，因此用來衡量原子鐘性能的參數(shù)都包含了對頻率的描述，原子時算法中還需要估計每個原子鐘頻率偏移參數(shù)。

類加權(quán)平均算法在計算時間尺度的過程上大致相同，ALOGOS采用之前的鐘差數(shù)據(jù)進行計算，AT1算法則是根據(jù)當前的鐘差值估計下一個時間的鐘差值。在鐘組內(nèi)選定一臺鐘作為主鐘，主鐘的物理信號輸入相位微躍器，微躍器可以調(diào)整主鐘的物理信號相位，微躍器和主鐘的組合稱為組合鐘，每臺原子鐘和組合鐘的鐘差值都可以通過與主鐘的鐘差間接計算得出。物理主鐘信號與相位微躍器調(diào)整后的主鐘信號的鐘差值可通過鐘i與主鐘的鐘差間接計算得出，每臺鐘的鐘差估值為xri(t+τ)，對鐘差估值xri(t+τ)加權(quán)平均計算出主鐘與組合鐘的鐘差xr(t+τ)，對相位微躍器則參照這個值對主鐘物理信號進行調(diào)整。

假設(shè)鐘i在時刻t的鐘差和速率分別為xi(t)，yi(t)，測量時間間隔為τ。則鐘組內(nèi)某臺鐘t+τ時刻相對于組合鐘的鐘差可由下式計算：

xi(t+τ)=xi(t)+yi(t+τ)·τ

(11)

將計數(shù)器在時刻t+τ測量得到的鐘i和主鐘的鐘差記為ti(t+τ)，則主鐘與組合鐘的鐘差可由下式計算：

xri(t+τ)=xi(t+τ)-ti(t+τ)

(12)

即：

xri(t+τ)=xi(t)+yi(t)·τ-ti(t+τ)

(13)

這樣就通過鐘i的鐘差和速率計算了組合鐘相對主鐘的鐘差值。

假設(shè)鐘組內(nèi)有N臺中，那么可以得到N-1個鐘差值，在對其進行加權(quán)平均可以得到：

(14) )

式中,xri(t+τ)為第i臺鐘間接計算出的鐘差估算值，ωi為它的權(quán)重，xi(t+τ)是鐘組內(nèi)所有原子鐘計算的得出的鐘差值的加權(quán)平均，xi(t+τ)是相位微躍器調(diào)整主鐘輸出信號的相位的參照，由此便生成了綜合原子時的物理信號。為了削弱噪聲對時間尺度穩(wěn)定性的影響，AT1算法需要引入一個指數(shù)濾波器，指數(shù)濾波器的計算方式為：

(15)

式中，k是一個重要參數(shù)，這里取k=30，第i臺鐘的頻率估計值為fi，通過時間差分計算得到：

(16)

由上述計算過程可以看出：AT1算法鐘每個原子鐘有兩個狀態(tài)，原子鐘相對主鐘的時間偏移和頻率偏移，頻率偏移可以用來計算，但算法本身不能對頻率偏移進行估計。改進的算法結(jié)合傳統(tǒng)的AT1算法和Kalman濾波來估計表示頻率隨機游走噪聲和原子鐘的頻率偏移狀態(tài)和方差。對于一個給定的原子鐘，這種狀態(tài)不是一種物理狀態(tài)，而是在白噪聲調(diào)制頻率的情況下對頻率偏移的數(shù)學(xué)估計，假設(shè)這種狀態(tài)稱之為“Y”，那么這個狀態(tài)的方差就表示對偏移的數(shù)學(xué)估計的置信度。

將AT1算法中的鐘i的鐘差寫為：

(17)

式中,Di是頻率漂移率。

如果鐘j在t+τ時刻的刻度鐘差測量值是原子鐘xij(t+τ)，則鐘j在t+τ時刻的鐘差預(yù)報偏差為：

xj(t+τ)=

(18)

原子鐘的權(quán)重計算方式通過時間殘差計算：

(19)

頻率偏移量可通過下式計算：

(20)

把頻率偏移量合并到鐘i當前平均頻率偏移的指數(shù)濾波估計中：

(21)

mi是鐘i的指數(shù)頻率加權(quán)常數(shù)，mi由鐘i的調(diào)頻白噪聲和隨機閃爍調(diào)頻白噪聲決定，用來衡量時間預(yù)測值的穩(wěn)定性。

其中：

(22)

τmin是σy(τ)到達最小值的時間，τ0是用于計算σy(τ)的最小時間。

接下來需要對狀態(tài)Y進行頻率偏移估計，文獻[12]和[13]對原子鐘的頻率穩(wěn)定度和鐘差預(yù)測不確定度進行了估計，采用一個簡單的卡爾曼形式對頻率偏移量進行過濾[21],假設(shè)Q(T)是測量噪聲，對狀態(tài)Y建模，狀態(tài)Y包含隨機游走噪聲和固定漂移，那么Y不是一個由原子鐘直接產(chǎn)生的物理層面的頻率，而是過濾了白色調(diào)頻噪聲的之后的頻率隨機游走分量加上頻率漂移。

將式(4)和(5)代入卡爾曼方程[4]，系統(tǒng)模型為：

yi(t+τ)=yi(t)+Diτ+η(τ)

(23)

測量模型：

(24)

式中,yi是頻率直接測量值，εi是白噪聲。

殘差估計：

(25)

其中：

(26)

這相當于一個自適應(yīng)卡爾曼濾波器，Y的殘差可以從初始值進行變化，Y的指數(shù)濾波器參數(shù)也可隨時間變化。

方程的穩(wěn)態(tài)形式中：

(27)

此原子鐘的指數(shù)頻率加權(quán)常數(shù)為：

(28)

頻率預(yù)測為：

(29)

由此可以看出卡爾曼形式也得到了類似AT1中的指數(shù)濾波器形式[19]，此方法繼承了AT1算法對閃變頻率的建模能力。

頻率跳變的檢測需要觀測檢查一段時間內(nèi)的數(shù)據(jù)，考慮到實時性的操作，可以將某一區(qū)間內(nèi)的平均頻率偏移估計值與該區(qū)間開始時的濾波估計值進行比較。將頻率偏差的估計值作為離群值進行測試，就可以確定頻率跳變是否在最近發(fā)生。

通過測量得到的是時鐘之間的相位差，由于受到頻率閃爍噪聲和隨機游走噪聲的影響，頻率跳變會在某一時刻發(fā)生，有些情況下時鐘本身小的頻率波動會被當成隨機噪聲，為了區(qū)分這種現(xiàn)象，需要對頻率跳變的大小進行限定。

本文的頻率跳變檢測方法是通過迭代每個時鐘的測量范圍來檢測頻率跳變，從當前測量時間到之后的某個時間，此時間間隔為τmin，Lmax是測量次數(shù)的最大值，Lmax是由τmin確定的整數(shù)。在每次測量過程中，對每臺時鐘的頻率跳變進行檢測，定義一個時鐘在測量時間間隔內(nèi)的平均頻率yavg：

(30)

x-L是在t-L時刻的估計值，x-1是在t-1時刻的估計值。通過將這個平均頻率和當前測量時刻的估計值進行比較，實際就是把固定時間區(qū)間上的實時頻率估計值與該區(qū)間的平均頻率估計值進行比較，如果兩者頻率差大于閾值，就可以認為時鐘發(fā)生了跳變,具體判斷方法為：

|yavg-y-L|>

(31)

通過這種方式檢測頻率跳變的發(fā)生，當檢測到某臺時鐘發(fā)生頻率跳變時，此時算法會將發(fā)生頻率跳變的時鐘從這個計算周期中剔除，也就是把發(fā)生頻率跳變的時鐘的權(quán)重變?yōu)?，此后的計算周期中對該時鐘的頻率繼續(xù)觀測，直到算法自適應(yīng)學(xué)習(xí)了此時鐘的新頻率特性，然后將此時鐘重新加入到原子時的計算當中。這種方法的好處在于原子鐘發(fā)生頻率跳變時，此臺鐘對綜合原子時的影響會降低，有利于提高綜合原子時的長期穩(wěn)定度。

3 守時系統(tǒng)硬件設(shè)計及原理

二級守時節(jié)點可分為硬件和軟件兩個部分，守時系統(tǒng)的硬件部分主要包括：原子鐘組、時頻信號產(chǎn)生與分配單元、時差測量單元、溯源比對單元、數(shù)據(jù)存儲分析處理和監(jiān)控單元。硬件部分主要考慮物理信號產(chǎn)生問題、與上級節(jié)點或同級節(jié)點間的數(shù)據(jù)交互問題和可拓展性等；軟件部分主要包括：溯源比對軟件、綜合原子時計算軟件、數(shù)據(jù)采集分析軟件和監(jiān)控軟件。穩(wěn)定可靠的守時系統(tǒng)由硬件和軟件共同支撐，缺一不可。守時系統(tǒng)模塊圖如圖1所示。

圖1 硬件系統(tǒng)結(jié)構(gòu)圖

本次試驗的守時系統(tǒng)采用4臺國產(chǎn)高性能銫原子鐘組成的鐘組，銫原子鐘獨立工作，不對原子鐘的數(shù)據(jù)在頻率和相位上進行改變。4臺銫原子鐘的原始信號直接輸出到頻率信號切換器和多通道計數(shù)器。4臺銫鐘輸出的1PPS信號直接進入多通道計數(shù)器，多通道計數(shù)器的作用就是對鐘差進行測量，測量結(jié)果通過UDP協(xié)議端口傳輸至數(shù)據(jù)庫。綜合原子時軟件部分每隔一定的周期會從數(shù)據(jù)庫中提取需要的鐘差數(shù)據(jù)進行綜合原子時的計算，并將計算結(jié)果和每臺鐘相對于綜合原子時的時差存入數(shù)據(jù)庫。

在參與守時的鐘組4臺銫原子鐘里面選用一個穩(wěn)定度較好的原子鐘作為主鐘也就是參考鐘。4臺銫原子鐘產(chǎn)生的10 M信號輸入頻率切換器，頻率切換器輸出主鐘物理信號，隨后進入相位微躍器，綜合原子時軟件對綜合原子時進行計算，同時相位微躍器依據(jù)綜合原子時與主鐘的鐘差對主鐘物理信號的相位進行調(diào)整，同時對主鐘產(chǎn)生的10 MHz信號進行駕馭，完成主鐘向綜合原子時的溯源。

4 仿真實驗

采用本單位4臺銫鐘的數(shù)據(jù)對本文算法進行了仿真驗證，銫鐘數(shù)據(jù)得采樣周期為60 s，采樣時間為120天數(shù)據(jù)點172 800個,計算周期為30天，軟件算法使用C語言完成，運行環(huán)境是國產(chǎn)麒麟操作系統(tǒng)。鐘差測量設(shè)備一般使用時間間隔計數(shù)器(SR620),SR620的觀測噪聲方差為σ2=1×10-20。分別采用AT1算法和本文改進算法對綜合原子時進行了計算，并對兩種方法的頻率偏差值進行了分析；對頻率跳變檢測進行仿真；對本文算法計算得出的原子時穩(wěn)定度驗證分析。

首先對正常狀態(tài)的銫鐘頻率預(yù)測值進行仿真分析，改進的方法對頻率預(yù)測值有明顯改善，頻率偏差曲線更加平滑,波動更小如圖2～3所示。

圖2 AT1算法頻率偏差

圖3 改進算法的頻率偏差

其中一臺銫鐘在第76天銫鐘頻率發(fā)生跳變，圖4是銫鐘的頻率跳變曲線，圖5是由本文檢測方法計算的平均頻率估計值，通過平均頻率計算的跳變值與閾值進行比較，可以判斷出銫鐘頻率發(fā)生了跳變，因此將發(fā)生跳變的銫鐘從當前原子時計算周期中剔除也即權(quán)重置0。

圖4 跳變銫鐘的頻率

圖5 跳變銫鐘的平均頻率

最后根據(jù)4臺銫鐘組成的鐘組系統(tǒng)，分別采用AT1算法和改進后的算法對時間尺度進行了計算，原子時穩(wěn)定度曲線如圖5所示。因為本文算法改善了頻率預(yù)測值且在一定程度上抑制了頻率跳變帶來的影響，所以理論上來講綜合原子時的中長期穩(wěn)定度比AT1算法更好，從圖6的曲線可以看出理論與仿真基本符合。

圖6 綜合原子時穩(wěn)定度

5 結(jié)束語

本文提出了一種Kalman調(diào)整頻率預(yù)測值和檢測頻率跳變的方法，和一般加權(quán)平均算法不同的是，一般的加權(quán)平均算法更加強調(diào)權(quán)重分配，以此提高時間尺度的穩(wěn)定度。本文則是直接改善頻率的預(yù)測值，本文算法核心思想是在每次頻率預(yù)報的過程中，采用Kalman濾波器對預(yù)測值進行改善，然后再計算綜合原子時，可以看到Kalman濾波器可以很好地抑制WFM和RWFM，隨著鐘組的老化，頻率發(fā)生變化，本文算法還可以對其檢測并重新適應(yīng)原子鐘的新頻率特性，所以算法的中長期穩(wěn)定度更高，但時間尺度的短期穩(wěn)定度提升不大，后續(xù)還需要進一步研究。除此之外，與Kalman算法直接計算出的綜合原子時不同，本文算法計算的時間尺度是連續(xù)、實時、可預(yù)測的，可以滿足規(guī)模較小的守時實驗室需求。