一種優化的頻率駕馭算法研究

趙書紅 董紹武*③④ 白杉杉 高 喆

①(中國科學院國家授時中心 西安710600)

②(中國科學院時間頻率基準重點實驗室 西安 710600)

③(中國科學院大學 北京100049)

④(中國科學院大學天文與空間科學學院 北京100049)

1 引言

時間是一個國家科技、經濟、軍事和社會生活中至關重要的參量，其應用范圍從基礎研究(相對論驗證、基礎物理常數測量等)，滲透到工程技術應用領域(導航定位、深空探測等)。特別在現代高科技領域和國防建設中，如全球衛星導航系統(Global Navigation Satellite System,GNSS)、現代軍事通信系統、航空航天系統等，高精度、穩定可靠的時間更是其未來技術迅速發展的核心要素。

目前國際標準時間是協調世界時(Coordinated Universal Time,UTC)，由國際權度局(International Bureau of Weights and Measures,BIPM)以Circular T公報形式發布。但UTC是“紙面時間”，滯后40～45 d，不滿足用戶對高精度時間信號的實時性要求，因此各守時實驗室需要獨立產生和保持一個穩定的標準時間系統UTC(k)，即UTC在本地的物理實現[1]。目前，全球多數時間實驗室利用一定數量的原子鐘(氫鐘、銫鐘等)聯合計算得到地方原子時尺度TA(k)，即“紙面時間”。以TA(k)作為駕馭參考，采用最小二乘估計方法預測待駕馭原子鐘的頻率，使其盡可能接近或符合紙面時間尺度的頻率，依據該值對原子鐘進行頻率駕馭，最終實現實時物理信號UTC(k)的輸出[2]。

實時物理信號UTC(k)產生的過程中，往往會有各種各樣的隨機因素存在，比如觀測噪聲、設備噪聲以及外部干擾等，同時這些隨機因素也是隨時間而變化的，但最小二乘估計方法沒有很好消除噪聲影響，且易受異常數據和異常信號的影響，因此現有的頻率駕馭算法還有改進的空間。自適應的控制方法是基于運行期間的測量值來改變增益矩陣，而不再是依據經驗來選擇，可以更好地控制系統達到期望的目標，因此本文采用基于自適應的控制方法，即最優二次型高斯控制算法(Linear Quadratic Gaussian control,LQG)，利用最小化二次代價函數不斷地逼近最優控制。再通過Kalman濾波方法，構造增益矩陣的遞推公式，將控制值和濾波器不斷迭代直至收斂，最終得到最優濾波器[3]。

本文提出了一種優化的頻率駕馭算法，主要分為紙面時間計算和實際物理信號實現兩個部分。第1步，基于實時的原子鐘數據和BIPM公布的Circular T公報的數據，采用ALGOS算法，產生一組氫銫聯合時間尺度，發揮了氫鐘和銫鐘的特性，保證駕馭參考的準確性和實時性。第2步，通過調整最優二次型高斯控制算法中的參數，使得代價函數最小化，再利用Kalman算法計算最優的氫鐘頻率駕馭量輸送至頻率調整設備中，保證輸出的實時物理信號性能，整個駕馭過程是優化、閉環的?；谖覈鴷r間基準保持系統，對一臺氫鐘進行了為期140 d的頻率駕馭試驗，分別在氫銫聯合時間尺度算法、LQG算法參數選擇以及優化的頻率駕馭算法等方面開展試驗，分析試驗結果，驗證了該算法的有效性。

2 基本原理

2.1 ALGOS算法

BIPM提出了ALGOS算法，目的是通過收集全世界時間實驗室的原子鐘數據，產生一個穩定的自由原子時(échelle Atomique Libre,EAL)，并利用多臺基準頻標的頻率(進行廣義相對論和黑體輻射改正后)的加權平均，對EAL的頻率進行駕馭得到既穩定又準確的時間尺度——國際原子時(International Atomic Time,TAI)。

若有N臺原子鐘，其讀數為hi(t)，i=1,2,···,N，利用加權平均算法，建立一個時間尺度TA(t)。ALGOS算法的基本公式如式(1)所示[4,5]

其中，Wmax表示最大權，A為需多考慮參與計算的原子鐘的性能的參數。

2.2 實際物理信號實現算法

實際物理信號實現算法是Kalman算法和LQG算法的綜合，本節重點介紹Kalman算法、LQG算法原理和實際物理信號實現算法流程。

2.2.1 Kalman算法

Kalman算法的狀態方程，如果假設為線性微分方程，僅描述了自由運轉的原子鐘模型，與頻率駕馭量的實際應用情況不符，因此狀態模型不能簡單用線性微分方程表示，必須增加控制項b(τ)u(t k)[6–8]。

含有控制量的狀態方程描述為

其中，X為原子鐘相位和頻率的狀態向量，t k為離散時間，τ=t k+1?t k為時間間隔，?(τ)為狀態轉移矩陣，b(τ)為 傳播項，u(t k)為tk時刻的頻率駕馭值，w(t k)是 原子鐘噪聲，在時域上互不相關，w(t k)～N(0,Q(τ))。

原子鐘的狀態方程僅考慮原子鐘本身的噪聲，不考慮測量過程中引入的噪聲。測量矩陣h描述系統測量過程，即描述待駕馭原子鐘與TA的相位偏差過程，并且測量噪聲v(t k)～N(0,R(t k))。Kalman算法的測量方程為

2.2.2最優二次型高斯控制算法

最優二次型高斯控制算法(LQG算法)是一種較為靈活的算法，通過調整WQ和WR參數值，優化頻率駕馭值。為保證LQG算法計算獲得最優的主鐘頻率駕馭值，必須保證代價函數J最小[9–11]。

2.2.3物理信號實現算法流程

物理信號實現算法包括預測和修正步驟。在預測步驟中，利用Kalman狀態方程確定下一時刻的估計。在修正步驟中，使用最新測量值修正模型預測的偏差[14,15]。

預測：

3 優化的頻率駕馭算法

優化的駕馭方法通過尋找駕馭函數u(t)，在滿足狀態方程式(7)和測量方程式(10)的條件下，使得代價函數J達到最小。針對實時物理信號UTC(k)產生的隨機控制問題，其全局最優頻率駕馭函數的確定是比較困難的。利用LQG算法使得最小化二次代價函數不斷地逼近最優控制，再通過Kalman濾波方法，構造增益矩陣的遞推公式，將控制值和濾波器不斷迭代直至收斂，最終得到最優濾波器。該算法可以很好地降低各種類型噪聲對系統的影響，并在最短時間可以消除異常值對頻率駕馭的影響，最終提高頻率駕馭后信號的穩定性和系統的可靠性。

優化的頻率駕馭算法分為紙面時間和實際物理信號實現兩個部分。紙面時間實現算法，基于氫原子鐘和銫原子鐘的原子鐘數據，采用ALGOS算法產生一個準確、穩定的紙面時間TA(k)，該紙面時間尺度的穩定度優于任何一臺原子鐘的穩定度。物理信號實現算法，以氫鐘與紙面時間TA(k)的相位偏差數據作為一組輸入數據，考慮到氫鐘的外部頻率駕馭操作，綜合Kalman算法和LQG算法計算主鐘的頻率駕馭量，將該值輸送至相位微調器中，最終實現物理信號輸出。該物理信號作為本地原子鐘數據采集的輸入信號，參與紙面時間計算，循環往復，整個過程是閉環的。氫鐘的頻率駕馭過程如圖1所示，圖中MC為氫鐘經頻率駕馭后的輸出信號。

4 測試及分析評估

4.1 氫銫聯合時間尺度TA產生

利用我國時間基準系統的原子鐘測量比對系統，選擇2019年1月1日(修正的儒略日期(Modified Julian Date,MJD)=58484)至2019年5月23日(MJD=58623)的原子鐘比對數據，采用在線連續運行的20臺原子鐘(13臺銫鐘和7臺氫鐘)與標準時間UTC(NTSC)的相位比對數據，采樣間隔為1 h。經鐘差預測、頻率預報以及權重估計等，利用ALGOS算法，計算得到紙面時間尺度TA。

由于環境異?；蛟O備故障等因素，可能造成原子鐘比對數據中出現缺失、異常數據等現象。為保障紙面時間尺度的準確性和連續性，必須對原始的原子鐘數據進行數據預處理，試驗中選用滑動的3σ粗差剔除法和三次樣條插值方法解決異常數據和缺失數據。

數據預處理后，基于ALGOS方法，利用7臺氫原子鐘，產生一個氫原子鐘組成的時間尺度TAHM。同時又利用13臺銫原子鐘組以及經過TAI速率校正后綜合產生的原子時尺度TACs。這兩個時間尺度互為參考，利用TAHM來減小銫原子鐘的短期波動，而利用TACs來修正TAHM的長期漂移，綜合兩者得到穩定度和準確度最優的時間尺度TA。由于參與計算的原子鐘數據為實時的，且計算前后的時間尺度經過鏈接，因此最終計算獲得的時間尺度TA是連續且實時的[16,17]。

圖1 氫鐘的頻率駕馭框圖

從圖2可以看出，氫銫聯合時間尺度TA與標準時間UTC(NTSC)的相位偏差保持在±5 ns以內，對比國際標準時間UTC與標準時間UTC(NTSC)的結果(數據來源于BIPM的Circular T公報)，兩組結果基本吻合。從圖3可以看出，氫銫聯合時間尺度TA與UTC的相位偏差均值為–0.38 ns，均方根誤差(Root Mean Square,RMS)為0.78 ns，并且30 d穩定度達到4.7×10–16[18,19]，因此氫銫聯合時間尺度TA是一個準確可靠的駕馭參考。

4.2 實時物理信號產生

4.2.1 LQG算法的參數確定

LQG算法通過調整WQ和WR參數，可以計算獲得最優的主鐘頻率駕馭量，其中WQ是具有兩個非零對角線的2×2矩陣，有兩個可變參數，WR是1維可變參數。給定WQ和WR參數，利用2.2.2節中式(12)、式(13)，可以計算出LQG算法的增益G?0,G?0=(g x,g y)，其中gx為 相位增益，gy為頻率增益。

圖4(a)和圖4(b)分別顯示了WR=109,WQ的兩個參數選擇范圍從(10–14,1014)取值時，相位增益g x和頻率增益gy的3維曲面圖。從圖中可以看出，當gx和gy取值接近零值時，隨著WR的增大，WQ(1,1)和WQ(2,2)的取值范圍變大。

WQ和WR參數選擇可以利用式(20)，當函數值為0時，選取參數的最優值[12]。

圖2 TA、UTC與UTC(NTSC)相位偏差結果

圖3 TA相對于UTC的相位偏差以及穩定度結果

圖4 WR值固定，WQ不同取值的增益

圖5 WQ值固定，WR不同取值的函數f變化曲線

4.2.2最優頻率駕馭量計算

圖7中黑色虛線代表氫鐘H067相較于紙面時間TA的相位偏差(扣除了二次趨勢)，紅色實線代表應用優化的頻率駕馭算法后的輸出信號與紙面時間TA的相位偏差。從圖中可以看出，采用優化的頻率駕馭算法，駕馭后信號一直穩定在零處。相比于優化的頻率駕馭算法，采用最小二乘估計方法駕馭，扣除了二次趨勢后，駕馭后的信號波動明顯增大，由此表明優化的頻率駕馭算法可以很好地實現原子鐘的頻率駕馭。

圖6 氫原子鐘的頻率駕馭量

當時間間隔在100～10000 s時，氫鐘噪聲以頻率白噪聲為主；當時間間隔在10000 s～7 d時，氫鐘噪聲以頻率閃爍噪聲為主；從圖8可以看出，利用優化的頻率駕馭算法，有效地抑制了頻率白噪聲和頻率閃爍噪聲，同時很好地消除了氫鐘的2次項趨勢，使得駕馭后的信號與氫銫聯合時間尺度TA的相位偏差控制在±1 ns以內，并提高了駕馭后信號的頻率穩定度[20,21]。

圖7 采用優化的頻率駕馭算法和最小二乘估計方法，駕馭后的輸出信號與TA的結果分析

圖8 采用優化的頻率駕馭算法和最小二乘估計方法，駕馭后的輸出信號與TA的穩定度結果

圖9 采用優化的頻率駕馭算法，駕馭后的輸出信號與UTC的結果分析

利用國家授時中心的溯源鏈路，將駕馭后的時間信號溯源至國際標準時間UTC。從圖9可以看出，駕馭后的實時物理信號與國際標準時間UTC相比，保持在±3 ns以內，且實時物理信號30 d的穩定度優于5×10–16[15]。圖10中帶圈虛線代表氫鐘H067的穩定度，實線代表了紙面時間TA的穩定度，以及帶框實線代表可駕馭后輸出信號的穩定度。從圖中可以看出，氫鐘受頻漂和隨機游走噪聲等的影響，頻率穩定度隨時間推移變差。多臺原子鐘綜合計算獲得的紙面時間尺度TA，很大程度上消除了頻漂和隨機游走噪聲的影響，短期穩定度和長期穩定度相比于氫鐘的穩定度均得到了提高。而應用優化的頻率駕馭算法，駕馭輸出后的物理信號穩定度與紙面時間TA的穩定度基本保持一致，駕馭后實時物理信號的穩定度優于氫鐘的自身穩定度。

圖10 以UTC作為參考，氫鐘H067，紙面時間TA和駕馭后實時物理信號的穩定度結果

5 結束語

隨著各領域對高精度時間的精度不斷提高的要求，更高精度的時間產生技術的研究受到廣泛的關注。頻率駕馭算法作為高精度時間產生的重點研究問題，其算法的有效性取決于駕馭參考的時間尺度、最優駕馭量計算等方面。因此本文提出了一種優化的頻率駕馭算法，綜合ALGOS算法、Kalman算法和LQG算法，迭代計算出最優頻率駕馭量，利用頻率調整設備，最終實現高精度時間信號輸出。

基于我國時間基準系統搭建了試驗平臺，利用優化的頻率駕馭算法，對一臺氫原子鐘進行140 d的頻率駕馭試驗，試驗結果說明，該算法可以提高輸出實時物理信號的準確度和穩定度，駕馭輸出的時間信號與國際標準時間UTC相比，相位偏差保持在±3 ns以內，且30 d穩定度優于5×10–16，為提高我國時間基準的性能提升奠定了基礎。