宋會杰,董紹武,2,王翔，王燕平,張繼海,屈俐俐，趙書紅,張首剛,2

(1.中國科學院 國家授時中心，西安 710600；2.中國科學院大學 天文與空間科學學院，北京 100049)

0 引言

原子鐘的噪聲及守時性能分析是守時工作中非常重要的一項研究工作。原子鐘噪聲是影響時間穩(wěn)定性的重要因素，分析原子鐘的噪聲特性，研究降低原子鐘噪聲的算法可以提升守時性能[1-4]。守時性能主要表現(xiàn)為產(chǎn)生標準時間的準確性、穩(wěn)定性和實時性。利用原子鐘的比對測量數(shù)據(jù)分析原子鐘的性能，研究守時算法提升標準時間長期性能[5-9]。當前國內(nèi)外研究人員對守時算法的研究不斷取得新的進展，比如2014年后，國際權度局(BIPM)對原子鐘進行二次模型預報并根據(jù)可預測性方法取權，提升了標準時間的整體性能[10-11]。目前我國采用的守時算法也是基于可預測性方法取權。但是上述研究和算法都是基于國外原子鐘，比如美國的銫原子鐘、氫原子鐘和俄羅斯的氫原子鐘等，由于國產(chǎn)原子鐘的工作原理、材料、制造工藝與國外原子鐘不同，導致噪聲特性存在一定差異，因此研究國產(chǎn)鐘的守時性能，提出符合國產(chǎn)鐘特性的守時算法，解決我國標準時間自主產(chǎn)生和時間基準保持的核心技術是當前面臨的迫切問題。文中基于小波分解研究了國產(chǎn)原子鐘不同尺度下的噪聲情況，并與國外原子鐘進行了比較，計算出不同尺度下的噪聲強度。通過經(jīng)典加權算法與Kalman濾波算法分別計算基于國產(chǎn)鐘的時間尺度，并比較了時間尺度的準確度、穩(wěn)定度和實時性，整體上分析了基于國產(chǎn)鐘的守時性能。

1 國產(chǎn)原子鐘的噪聲分析

小波變換作為一種時頻局部化方法，其窗口是可變的，即在低頻部分具有較低的時間分辨率和較高的頻率分辨率，具有對信號的自適應性。應用于國產(chǎn)原子鐘的測量數(shù)據(jù)，分析不同尺度的噪聲強度。

1.1 小波變換的原理

(1)

‖Ψjk(t)‖=‖Ψ(t)‖j,k∈Z，

(2)

通過對國產(chǎn)鐘測量數(shù)據(jù)進行離散小波變換，可以對數(shù)據(jù)進行時頻分析，即測量數(shù)據(jù)在某一時間上，對應某一尺度a0的成分。相應的離散小波變換可表示為

(3)

對于離散小波變換，一個3尺度的分解如圖1所示。從圖1可以看出小波變換分離信號的不同成分，其中不同尺度的低頻系數(shù)表示去噪后的數(shù)據(jù)，不同尺度的高頻系數(shù)表示噪聲，對應國產(chǎn)原子鐘不同頻率段的噪聲。

圖1 3尺度分解的組織形式

1.2 基于小波變換的國產(chǎn)鐘噪聲分析

基于小波變換理論對國產(chǎn)鐘的測量數(shù)據(jù)進行噪聲分析，并與國外原子鐘的測量噪聲進行比較。國產(chǎn)原子鐘編號選為Cs2025、Cs3050和Cs3059，國外原子鐘編號選為Cs3089。測量數(shù)據(jù)時間為2020年8月至2021年6月，測量參考為UTC(NTSC)。利用db5小波對國產(chǎn)鐘測量數(shù)據(jù)進行6尺度分解，得到不同尺度的原子鐘噪聲，其結果如圖2，圖3，圖4和圖5所示。為了清晰，圖中只表示了3尺度分解。通過小波多尺度分解得出不同尺度噪聲的標準差，其結果如表1和表2所示。結合圖與表可以看出尺度越低，噪聲越小，隨著尺度的增大，噪聲變大。小波分解中尺度越小，對應的頻率越高，尺度越大，對應的頻率越低。隨著頻率的降低，噪聲變大，說明國產(chǎn)鐘在低頻段噪聲明顯，低頻段的噪聲對國產(chǎn)鐘性能的影響起主要作用。守時算法需重點考慮降低國產(chǎn)鐘低頻噪聲的影響。比較國產(chǎn)鐘Cs2025，Cs3050和Cs3059不同尺度的噪聲大小，在尺度1至尺度3，Cs2025和Cs3050的噪聲明顯低于Cs3059的噪聲，說明國產(chǎn)鐘Cs3059的高頻噪聲相對偏大。分析尺度4至6，Cs2025的噪聲明顯高于Cs3050和Cs3059，說明國產(chǎn)鐘Cs2025的低頻噪聲相對偏大。3臺原子鐘，Cs3050的低頻噪聲與高頻噪聲都相對偏小，表示Cs3050的長期穩(wěn)定度和短期穩(wěn)定度最高。通過上述分析，得出不同原子鐘的噪聲情況不同，部分國產(chǎn)鐘低頻噪聲相對偏大，部分國產(chǎn)鐘高頻噪聲相對偏大。守時過程中，根據(jù)國產(chǎn)原子鐘不同尺度的噪聲情況進行加權處理，能夠同時降低不同尺度噪聲的影響，提升守時的整體性能。

圖2 Cs2025不同尺度的噪聲

圖3 Cs3050不同尺度的噪聲

圖4 Cs3059不同尺度的噪聲

圖5 Cs3089不同尺度的噪聲

表1 原子鐘尺度(1～3)的噪聲標準差 單位：ns

表2 原子鐘尺度(4～6)的噪聲標準差 單位：ns

2 國產(chǎn)原子鐘的守時性能研究

守時性能通常表現(xiàn)為多臺原子鐘產(chǎn)生時間尺度的準確性、穩(wěn)定性和實時性。在時間保持工作中，時間尺度相當于一臺虛擬原子鐘，這臺虛擬原子鐘是通過測量鐘差實現(xiàn)的物理鐘的加權平均。虛擬鐘通過與某一臺鐘的偏差來表示。一般認為虛擬鐘的性能優(yōu)于單臺原子鐘的性能[5]。下面主要通過時間尺度算法研究國產(chǎn)原子鐘的守時性能，重點研究一種Kalman濾波時間尺度算法用于國產(chǎn)鐘的守時性能，并與一種經(jīng)典的加權平均算法的守時性能進行比較。

2.1 鐘的模型

Kalman濾波時間尺度算法是一種實時的原子鐘狀態(tài)估計方法，估計參考鐘和理想鐘之間的差，并將此值作為修正量，計算時間尺度。算法原理如下,原子鐘的狀態(tài)方程可表示為：

(4)

y(tk+1)=y(tk)+δω(tk)+η(tk),

(5)

ω(tk+1)=ω(tk)+α(tk),

(6)

式(4)至(6)中,x(tk)，y(tk)和ω(tk)分別表示鐘在tk時刻的相位偏差，頻率偏差和頻率漂移偏差。tk和tk+1是相鄰取樣時間，具有如下關系，tk+1=tk+δ。ε(tk)，η(tk)和α(tk)分別表示tk時刻原子鐘的調(diào)頻白噪聲，調(diào)頻隨機游走噪聲和頻漂隨機游走噪聲。原子在tk時刻的測量模型表示為

zij(tk)=xi(tk)-xj(tk)+v(tk),

(7)

式(7)中，zij表示第ith臺鐘和第jth臺鐘之間的相位差測量，v表示測量噪聲。

Kalman濾波時間尺度問題是估計每臺原子鐘的狀態(tài)，但是實際測量數(shù)據(jù)是鐘差數(shù)據(jù)，因此建立鐘差的狀態(tài)方程估計原子鐘的狀態(tài)。通過方程(4)，(5)和(6)對鐘i和鐘j做差的鐘差方程表示為：

(8)

yij(tk+1)=yij(tk)+δωij(tk)+ηij(tk),

(9)

ωij(tk+1)=ωij(tk)+αij(tk)。

(10)

根據(jù)N-1個鐘差測量估計N臺鐘的狀態(tài)，這種情況的解是不穩(wěn)定的。根據(jù)中心極限定理，當原子鐘數(shù)量很多時，噪聲的加權和為0。即有下列輔助方程：

(11)

(12)

(13)

式(11)至(13)中,ai(tk)，bi(tk)和ci(tk)分別表示tk時刻鐘i噪聲的權重。根據(jù)鐘差方程和輔助方程估計單臺原子鐘的狀態(tài)。

2.2 Kalman濾波時間尺度算法原理

首先定義tk鐘差狀態(tài)變量和噪聲向量為：

(14)

鐘差方程的向量形式表示為

(15)

(16)

式(16)中，H表示測量矩陣，R是測量標準偏差。Kalman濾波的遞推方程表示為：

(17)

(18)

(19)

(20)

Kij表示Kalman增益，表示為一個3-元素列向量：

(21)

式(20)結合式(14)的估計噪聲向量為

(22)

2.3 國產(chǎn)原子鐘的守時性能研究

研究國產(chǎn)原子鐘的守時性能，基于不同時間尺度算法分析國產(chǎn)原子鐘的準確性、穩(wěn)定性和實時性。選取國產(chǎn)原子鐘，編號為Cs3050，Cs2025和Cs3059計算時間尺度，計算數(shù)據(jù)時間段為：2020年8月至2021年6月，數(shù)據(jù)采樣周期為1 h，國產(chǎn)鐘參考UTC(NTSC)的相對頻率偏差數(shù)據(jù)如圖6所示，為了表示清楚，Cs3059的相對頻率偏差都減去1×10-12。從圖6可以看出Cs2025的相對頻率偏差相比Cs3050與Cs3059明顯偏大，Cs3050與Cs3059的相對頻率偏差相當。文中采用上述給出的Kalman濾波算法和經(jīng)典加權平均(類ALGOS)算法計算國產(chǎn)原子鐘的時間尺度，并分析時間尺度的準確度、穩(wěn)定度和實時性?；趦煞N算法計算的時間尺度如圖7所示，得出本文采用的類ALGOS算法計算的時間尺度長期波動較大，基于上述Kalman濾波算法的時間尺度長期波動相比本文采用的類ALGOS算法明顯減小。分析原因，本文采用的類ALGOS算法采用前一個月的速率作為本月的預報值，如果本月速率相比前一個月的速率有一個較大的變化，會引起時間尺度較大的波動。Kalman濾波算法是一種實時算法，給定初始狀態(tài)下根據(jù)測量鐘差數(shù)據(jù)實時修正原子鐘的狀態(tài)估計，提高了時間尺度估計的準確性。表3給出了基于兩種算法的時間尺度準確度指標，本文采用的類ALGOS算法的時間尺度與參考UTC(NTSC)的最大偏差為216.42 ns,明顯大于基于Kalman算法的最大偏差128.42 ns；另外，類ALGOS算法的時間尺度與參考UTC(NTSC)偏差的平均值為33.19 ns，基于Kalman算法的時間尺度與參考UTC(NTSC)偏差的平均值為39.40 ns，兩者相差不大；最后，基于本文類ALGOS算法的時間尺度的RMS為89.17 ns，明顯大于基于Kalman濾波算法的RMS值52.63 ns，說明基于本文的類ALGOS算法的時間尺度長期波動大?；谝陨戏治觯疚牟捎玫腒alman濾波算法優(yōu)于類ALGOS算法計算的時間尺度性能。

圖6 國產(chǎn)原子鐘的相對頻率偏差數(shù)據(jù)

圖7 基于兩種算法的國產(chǎn)鐘時間尺度

表3 基于兩種算法的時間尺度的準確度 單位：ns

分析基于兩種算法的時間尺度的穩(wěn)定度，穩(wěn)定度通過重疊Allan偏差表示，重疊Allan偏差越小穩(wěn)定度越高,其結果如圖8所示。從圖8可以得出，不同國產(chǎn)原子鐘穩(wěn)定度差別明顯，Cs2025和Cs3050穩(wěn)定度性能10天內(nèi)相差不大。10～83 d天內(nèi)，Cs3050的穩(wěn)定度明顯優(yōu)于Cs2025,Cs3059穩(wěn)定度在相應取樣時間差于Cs2025和Cs3050。說明相比較于Cs3050，國產(chǎn)鐘Cs2025具有較差的長期穩(wěn)定度，Cs3059短期穩(wěn)定度和長期穩(wěn)定度都明顯低于其余兩臺國產(chǎn)原子鐘。基于文中Kalman濾波算法計算國產(chǎn)鐘的時間尺度，算法建立了鐘差的相位偏差，頻率偏差和頻率漂移偏差的狀態(tài)方程，并且狀態(tài)方程中包括鐘差的調(diào)頻白噪聲，調(diào)頻隨機游走噪聲和頻漂隨機游走噪聲。Kalman濾波算法通過使線性方差取極小值，同時降低了鐘差狀態(tài)方程中的多種噪聲，提升了時間尺度的穩(wěn)定度。圖8給出了本文采用的類ALGOS算法計算的時間尺度的重疊Allan偏差和文中Kalman濾波算法計算的時間尺度的重疊Allan偏差。本文采用的類ALGOS算法計算時間尺度的重疊Allan偏差在相應取樣時間均小于單臺原子鐘的重疊Allan偏差，30 d的重疊Allan偏差為2.9×10-14，83 d的重疊Allan偏差為1.1×10-14?；贙alman濾波算法計算時間尺度在取樣時間超過5 d后，重疊Allan偏差明顯減小，穩(wěn)定度顯著提高，30 d的重疊Allan偏差為1.0×10-14，83 d的重疊Allan偏差為5.5×10-15。對于取樣時間小于5 d，基于Kalman濾波算法計算時間尺度的重疊Allan偏差高于類ALGOS算法計算的重疊Allan偏差，甚至高于單臺原子鐘的重疊Allan偏差。分析原因，采用文中Kalman濾波時間尺度算法，首先需要給定鐘差數(shù)據(jù)的初始狀態(tài)，基于Kalman濾波算法估計預測狀態(tài)，根據(jù)量測數(shù)據(jù)實時修正預測狀態(tài)，給出鐘差的估計，這是一個遞推的過程。當原子鐘可預測性變差時，Kalman濾波中量測數(shù)據(jù)不斷的修正預測值，用于提高估計的準確性，這樣就降低了時間尺度的短期穩(wěn)定度。本文采用的類ALGOS算法計算時間尺度，根據(jù)前一個月的速率作為當前月的速率預報值，速率是恒定的，然后根據(jù)月速率的變化計算權重，這樣降低了對時間尺度短期穩(wěn)定度的影響。

圖8 國產(chǎn)原子鐘與時間尺度的重疊Allan偏差

對于實時性方面，本文采用的類ALGOS算法需要利用前一個月的速率作為當月速率的預報值，是一種滯后的時間尺度算法。Kalman濾波算法根據(jù)預測方程和量測方程遞推的估計原子鐘的狀態(tài)，是一種實時的時間尺度算法。

3 結語

本文基于小波分解得到國產(chǎn)原子鐘和國外原子鐘的相應分解尺度的噪聲，比較不同分解尺度的噪聲大小，分析國產(chǎn)鐘的性能。守時過程中，通過對多臺國產(chǎn)鐘在不同分解尺度加權處理，能夠降低不同分解尺度的噪聲，提升守時的性能。時間尺度的準確度、穩(wěn)定度和實時性是守時性能的一個重要表現(xiàn)，時間尺度是由多臺原子鐘通過算法產(chǎn)生的一個紙面時間，國產(chǎn)鐘守時的性能與時間尺度的性能有直接的關系。本文研究了基于國產(chǎn)鐘的Kalman濾波時間尺度算法和本文采用的類ALGOS時間尺度算法，通過分析得出，Kalman濾波算法計算得到的時間尺度的準確度優(yōu)于類ALGOS算法的準確度。5 d內(nèi)，類ALGOS算法計算的時間尺度的穩(wěn)定度優(yōu)于文中Kalman濾波算法計算的穩(wěn)定度；5～80 d，Kalman濾波算法計算的時間尺度的穩(wěn)定度優(yōu)于類ALGOS算法計算的穩(wěn)定度。在具體實施方面，Kalman濾波算法具有較強的實時性。