一種優(yōu)化的頻率駕馭算法研究

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

①(中國科學(xué)院國家授時(shí)中心 西安710600)

②(中國科學(xué)院時(shí)間頻率基準(zhǔn)重點(diǎn)實(shí)驗(yàn)室 西安 710600)

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

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

1 引言

時(shí)間是一個(gè)國家科技、經(jīng)濟(jì)、軍事和社會生活中至關(guān)重要的參量，其應(yīng)用范圍從基礎(chǔ)研究(相對論驗(yàn)證、基礎(chǔ)物理常數(shù)測量等)，滲透到工程技術(shù)應(yīng)用領(lǐng)域(導(dǎo)航定位、深空探測等)。特別在現(xiàn)代高科技領(lǐng)域和國防建設(shè)中，如全球衛(wèi)星導(dǎo)航系統(tǒng)(Global Navigation Satellite System,GNSS)、現(xiàn)代軍事通信系統(tǒng)、航空航天系統(tǒng)等，高精度、穩(wěn)定可靠的時(shí)間更是其未來技術(shù)迅速發(fā)展的核心要素。

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

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

本文提出了一種優(yōu)化的頻率駕馭算法，主要分為紙面時(shí)間計(jì)算和實(shí)際物理信號實(shí)現(xiàn)兩個(gè)部分。第1步，基于實(shí)時(shí)的原子鐘數(shù)據(jù)和BIPM公布的Circular T公報(bào)的數(shù)據(jù)，采用ALGOS算法，產(chǎn)生一組氫銫聯(lián)合時(shí)間尺度，發(fā)揮了氫鐘和銫鐘的特性，保證駕馭參考的準(zhǔn)確性和實(shí)時(shí)性。第2步，通過調(diào)整最優(yōu)二次型高斯控制算法中的參數(shù)，使得代價(jià)函數(shù)最小化，再利用Kalman算法計(jì)算最優(yōu)的氫鐘頻率駕馭量輸送至頻率調(diào)整設(shè)備中，保證輸出的實(shí)時(shí)物理信號性能，整個(gè)駕馭過程是優(yōu)化、閉環(huán)的。基于我國時(shí)間基準(zhǔn)保持系統(tǒng)，對一臺氫鐘進(jìn)行了為期140 d的頻率駕馭試驗(yàn)，分別在氫銫聯(lián)合時(shí)間尺度算法、LQG算法參數(shù)選擇以及優(yōu)化的頻率駕馭算法等方面開展試驗(yàn)，分析試驗(yàn)結(jié)果，驗(yàn)證了該算法的有效性。

2 基本原理

2.1 ALGOS算法

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

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

其中，Wmax表示最大權(quán)，A為需多考慮參與計(jì)算的原子鐘的性能的參數(shù)。

2.2 實(shí)際物理信號實(shí)現(xiàn)算法

實(shí)際物理信號實(shí)現(xiàn)算法是Kalman算法和LQG算法的綜合，本節(jié)重點(diǎn)介紹Kalman算法、LQG算法原理和實(shí)際物理信號實(shí)現(xiàn)算法流程。

2.2.1 Kalman算法

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

含有控制量的狀態(tài)方程描述為

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

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

2.2.2最優(yōu)二次型高斯控制算法

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

2.2.3物理信號實(shí)現(xiàn)算法流程

物理信號實(shí)現(xiàn)算法包括預(yù)測和修正步驟。在預(yù)測步驟中，利用Kalman狀態(tài)方程確定下一時(shí)刻的估計(jì)。在修正步驟中，使用最新測量值修正模型預(yù)測的偏差[14,15]。

預(yù)測：

3 優(yōu)化的頻率駕馭算法

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

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

4 測試及分析評估

4.1 氫銫聯(lián)合時(shí)間尺度TA產(chǎn)生

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

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

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

圖1 氫鐘的頻率駕馭框圖

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

4.2 實(shí)時(shí)物理信號產(chǎn)生

4.2.1 LQG算法的參數(shù)確定

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

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

WQ和WR參數(shù)選擇可以利用式(20)，當(dāng)函數(shù)值為0時(shí)，選取參數(shù)的最優(yōu)值[12]。

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

圖3 TA相對于UTC的相位偏差以及穩(wěn)定度結(jié)果

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

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

4.2.2最優(yōu)頻率駕馭量計(jì)算

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

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

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

圖7 采用優(yōu)化的頻率駕馭算法和最小二乘估計(jì)方法，駕馭后的輸出信號與TA的結(jié)果分析

圖8 采用優(yōu)化的頻率駕馭算法和最小二乘估計(jì)方法，駕馭后的輸出信號與TA的穩(wěn)定度結(jié)果

圖9 采用優(yōu)化的頻率駕馭算法，駕馭后的輸出信號與UTC的結(jié)果分析

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

圖10 以UTC作為參考，氫鐘H067，紙面時(shí)間TA和駕馭后實(shí)時(shí)物理信號的穩(wěn)定度結(jié)果

5 結(jié)束語

隨著各領(lǐng)域?qū)Ω呔葧r(shí)間的精度不斷提高的要求，更高精度的時(shí)間產(chǎn)生技術(shù)的研究受到廣泛的關(guān)注。頻率駕馭算法作為高精度時(shí)間產(chǎn)生的重點(diǎn)研究問題，其算法的有效性取決于駕馭參考的時(shí)間尺度、最優(yōu)駕馭量計(jì)算等方面。因此本文提出了一種優(yōu)化的頻率駕馭算法，綜合ALGOS算法、Kalman算法和LQG算法，迭代計(jì)算出最優(yōu)頻率駕馭量，利用頻率調(diào)整設(shè)備，最終實(shí)現(xiàn)高精度時(shí)間信號輸出。

基于我國時(shí)間基準(zhǔn)系統(tǒng)搭建了試驗(yàn)平臺，利用優(yōu)化的頻率駕馭算法，對一臺氫原子鐘進(jìn)行140 d的頻率駕馭試驗(yàn)，試驗(yàn)結(jié)果說明，該算法可以提高輸出實(shí)時(shí)物理信號的準(zhǔn)確度和穩(wěn)定度，駕馭輸出的時(shí)間信號與國際標(biāo)準(zhǔn)時(shí)間UTC相比，相位偏差保持在±3 ns以內(nèi)，且30 d穩(wěn)定度優(yōu)于5×10–16，為提高我國時(shí)間基準(zhǔn)的性能提升奠定了基礎(chǔ)。