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

孫同川，王振嶺，孫建設，劉鐵強

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

0 引言

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

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

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

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

1 原子鐘噪聲模型的研究

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

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

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

(1)

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

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

(2)

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

(3)

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

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

假設原子鐘的噪聲模型符合：

(4)

(5)

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

由Kalman濾波原理，得到鐘組的動態模型為：

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

(6)

式中，φ為轉移矩陣，W(t)為系統噪聲矩陣。

鐘組的觀測方程為：

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

(7)

(8)

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

(9)

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

(10)

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

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

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

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

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

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

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

(11)

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

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

(12)

即：

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

(13)

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

假設鐘組內有N臺中，那么可以得到N-1個鐘差值，在對其進行加權平均可以得到：

(14) )

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

(15)

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

(16)

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

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

(17)

式中,Di是頻率漂移率。

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

xj(t+τ)=

(18)

原子鐘的權重計算方式通過時間殘差計算：

(19)

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

(20)

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

(21)

mi是鐘i的指數頻率加權常數，mi由鐘i的調頻白噪聲和隨機閃爍調頻白噪聲決定，用來衡量時間預測值的穩定性。

其中：

(22)

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

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

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

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

(23)

測量模型：

(24)

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

殘差估計：

(25)

其中：

(26)

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

方程的穩態形式中：

(27)

此原子鐘的指數頻率加權常數為：

(28)

頻率預測為：

(29)

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

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

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

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

(30)

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

|yavg-y-L|>

(31)

通過這種方式檢測頻率跳變的發生，當檢測到某臺時鐘發生頻率跳變時，此時算法會將發生頻率跳變的時鐘從這個計算周期中剔除，也就是把發生頻率跳變的時鐘的權重變為0，此后的計算周期中對該時鐘的頻率繼續觀測，直到算法自適應學習了此時鐘的新頻率特性，然后將此時鐘重新加入到原子時的計算當中。這種方法的好處在于原子鐘發生頻率跳變時，此臺鐘對綜合原子時的影響會降低，有利于提高綜合原子時的長期穩定度。

3 守時系統硬件設計及原理

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

圖1 硬件系統結構圖

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

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

4 仿真實驗

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

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

圖2 AT1算法頻率偏差

圖3 改進算法的頻率偏差

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

圖4 跳變銫鐘的頻率

圖5 跳變銫鐘的平均頻率

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

圖6 綜合原子時穩定度

5 結束語

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