何文彥HE Wen-yan；雷正偉LEI Zheng-wei；甄紅濤ZHEN Hong-tao；李凱麗LI Kai-li；劉雅迪LIU Ya-di

（①西安航光衛(wèi)星測控技術(shù)有限公司，西安 710100；②西安電子科技大學(xué)，西安 710100）

0 引言

在我國軍事和航天領(lǐng)域中，導(dǎo)航定位和授時(shí)系統(tǒng)在北斗信息服務(wù)中占據(jù)著越來越重要的地位。而時(shí)鐘同步是精密授時(shí)的根本保證[1]。銣鐘短期穩(wěn)定度較好，是常用的同步原子時(shí)鐘[2-4]。時(shí)間統(tǒng)一系統(tǒng)包含時(shí)統(tǒng)中心和若干時(shí)統(tǒng)分中心，是為測控系統(tǒng)提供統(tǒng)一標(biāo)準(zhǔn)時(shí)間信號和標(biāo)準(zhǔn)頻率信號的系統(tǒng)，其設(shè)備由無線電接收機(jī)、原子頻率標(biāo)準(zhǔn)、標(biāo)準(zhǔn)信號發(fā)生器和放大分配設(shè)備組成，無線電接收機(jī)通過接收衛(wèi)星信號校準(zhǔn)本地銣鐘頻率[5]。世界各國以性能優(yōu)異的原子鐘建立并保持著本國的時(shí)間體系。高精度時(shí)間統(tǒng)一系統(tǒng)設(shè)備作為授時(shí)設(shè)備，能夠通過不同的鏈路接口對不同的設(shè)備進(jìn)行授時(shí)。本文分析了一種北斗時(shí)統(tǒng)設(shè)備銣鐘馴服技術(shù)，該技術(shù)基于FPGA 獲取銣鐘鐘頻差值，并利用vondrak 濾波和kalman 濾波算法對鐘差值進(jìn)行處理，然后利用PID 控制算法對銣鐘電壓進(jìn)行控制，調(diào)節(jié)銣鐘頻率，實(shí)現(xiàn)時(shí)間同步、提高守時(shí)設(shè)備精度。

1 銣鐘馴服技術(shù)實(shí)現(xiàn)原理

銣鐘馴服技術(shù)由FPGA 和MCU 實(shí)現(xiàn)。FPGA 根據(jù)時(shí)間間隔測量原理，計(jì)算北斗接收機(jī)的秒脈沖與銣鐘分頻的秒脈沖差值，并采用kalman 算法和vondrak 算法聯(lián)合對此鐘差值濾波，再通過uart 發(fā)給MCU；MCU 通過PID 控制算法，對DAC 轉(zhuǎn)換芯片進(jìn)行精準(zhǔn)的電壓控制，從而準(zhǔn)確地改變銣鐘頻率，達(dá)到馴服銣鐘的目的。

1.1 時(shí)間間隔測量原理 時(shí)間間隔測量由FPGA 實(shí)現(xiàn)，原理如圖1 所示，TX為測量時(shí)間間隔值，T0為采用時(shí)鐘計(jì)數(shù)法計(jì)算出的時(shí)間間隔值，計(jì)數(shù)法存在量化誤差：T1（開始信號與計(jì)數(shù)時(shí)鐘信號不同步部分）和T2（結(jié)束信號與計(jì)數(shù)時(shí)鐘信號不同步部分），關(guān)系如下：

圖1 時(shí)間間隔測量原理圖

T1和T2兩值可由時(shí)間-幅度轉(zhuǎn)換法計(jì)算。計(jì)數(shù)法參考時(shí)鐘為銣鐘倍頻的100MHz 時(shí)鐘，若其在由開始信號和結(jié)束信號構(gòu)成的時(shí)間閘門內(nèi)計(jì)得時(shí)鐘個(gè)數(shù)為N0，則T0=N0×T，這里T 為10ns。若時(shí)間閘門內(nèi)時(shí)鐘個(gè)數(shù)理論值為N，則鐘差值為：±（N-N0）/TX。

1.2 kalman 濾波原理 kalman 濾波算法以最小均方誤差為準(zhǔn)則，根據(jù)前一估計(jì)值和最近一觀察數(shù)據(jù)進(jìn)行估計(jì)信號當(dāng)前值。所以該算法無需知道過去全部觀察數(shù)據(jù)，采用狀態(tài)方程和迭代方法，以估計(jì)值的形式給出濾波值。

圖2 遞歸描繪最小均方估計(jì)問題

kalman 濾波算法有時(shí)間更新和測量更新兩部分，時(shí)間更新有兩個(gè)預(yù)測方程：

測試更新有三個(gè)狀態(tài)方程：

本項(xiàng)目中，需要根據(jù)鐘差實(shí)測值和濾波效果，調(diào)節(jié)過程噪聲協(xié)方差Q 和測量噪聲協(xié)方差R。時(shí)間更新的兩個(gè)公式通過在上一次的鐘差數(shù)據(jù)計(jì)算最優(yōu)化估計(jì)和其協(xié)方差的基礎(chǔ)上，進(jìn)行預(yù)測計(jì)算，得出此次的鐘差數(shù)據(jù)最優(yōu)估計(jì)和其協(xié)方差的預(yù)測值。在式（2）中，為鐘差數(shù)據(jù)的預(yù)測值，A 為狀態(tài)矩陣，為上一次鐘差數(shù)據(jù)得到的最優(yōu)估計(jì)，B 為控制矩陣，uK-1為控制輸入。在式（3）中為鐘差數(shù)據(jù)的預(yù)測值的先驗(yàn)估計(jì)協(xié)方差，A 為狀態(tài)轉(zhuǎn)移矩陣，為了簡便計(jì)算，A 和AT取1。Q 為過程噪聲的方差，是需要根據(jù)實(shí)驗(yàn)和仿真濾波效果，對此參數(shù)進(jìn)行確定。

狀態(tài)更新的三個(gè)公式，分別計(jì)算了本次鐘差數(shù)據(jù)的卡爾曼增益、最優(yōu)估計(jì)以及最優(yōu)估計(jì)的協(xié)方差矩陣。在式（4）中，因?yàn)椋柭鲆鍷K可以擴(kuò)展為：

為了簡便運(yùn)算，取維數(shù)矩陣H 和狀態(tài)矩陣A 為1，因此HT和AT也為1，可得簡化的卡爾曼增益為：

可根據(jù)式（8）對Q 和R 參數(shù)進(jìn)行調(diào)節(jié)，從而可以調(diào)整卡爾曼增益KK的大小。在式（5）中為此次鐘差數(shù)據(jù)的最優(yōu)估計(jì)，調(diào)整卡爾曼增益KK可以調(diào)整預(yù)測值和觀測值yK的占比，從而調(diào)整出最優(yōu)估計(jì)值。在式（6）中，最優(yōu)估計(jì)協(xié)方差矩陣PK可以通過卡爾曼增益KK和最優(yōu)估計(jì)協(xié)方差矩陣估計(jì)值得到[6]。

此濾波算法在FPGA 中實(shí)現(xiàn)，INT_FLOAT 模塊將輸入鐘差數(shù)據(jù)轉(zhuǎn)換為符合IEEE 754 標(biāo)準(zhǔn)的浮點(diǎn)數(shù)，Kalman_forecastOne 模塊和Kalman_forecastTwo 模塊分別完成了對鐘差數(shù)據(jù)的預(yù)測值以及鐘差數(shù)據(jù)的預(yù)測值的先驗(yàn)估計(jì)協(xié)方差和卡爾曼增益的計(jì)算，Kalman_update 模塊完成鐘差數(shù)據(jù)的最優(yōu)估計(jì)以及其協(xié)方差矩陣的計(jì)算，最后通過模塊FLOAT_TO_INT 模塊將計(jì)算得出的最優(yōu)估計(jì)轉(zhuǎn)換為整數(shù)輸出。

1.3 vondrak 濾波原理 vondrak 濾波在不知道觀測數(shù)據(jù)變化規(guī)律及其擬合函數(shù)時(shí)，就能對觀測數(shù)據(jù)平滑濾波。基本思想：通過不同的平滑因子控制數(shù)據(jù)的平滑程度，在絕對擬合與絕對平滑間選擇一條折衷的曲線，保留有用信號，濾除噪聲信號。

修勻數(shù)學(xué)中指出一個(gè)修勻序列滿足光滑性要求，所謂光滑性要求，就是在滿足S 很小的情況下，所得出的修勻值位于一條光滑的曲線上[7]。

實(shí)際處理時(shí)，過度追求光滑性會(huì)使估計(jì)值偏離原始觀測值，則處理的結(jié)果與真值的偏離程度增大?？深A(yù)先采用一定的處理方法，去除粗大誤差，使得觀測值與真值的偏離程度很小。

考慮到整體加權(quán)及正負(fù)抵消的情況下，令

上式中F 為擬合度，yi為觀測值，為修勻值。當(dāng)F 為0 時(shí)，說明與yi完全重合，則沒有起到修勻效果。

Whittaker 于1923 年將式（9）、式（10）做了線性組合，發(fā)展為Whittaker 修勻：

其中：pi為觀測值的權(quán)序列，yi為觀測值序列，為平滑值序列，h 是正常數(shù)。上述準(zhǔn)則下得到的修勻值既滿足位于一條光滑的曲線上，又與原觀測序列具有一定的擬合度。

捷克天文學(xué)家J.Vondrak 對Whittaker 修勻進(jìn)行改進(jìn)為Vondrak 濾波：

1976 年再對上式進(jìn)行了改進(jìn)，將擬合度和平滑度分別用平均值來代替，得：

其中，pi為觀測值的權(quán)序列，yi為觀測值序列，為平滑序列，λ 是正常數(shù)為平滑因子（一個(gè)給定的無量綱的正數(shù)），在觀測數(shù)據(jù)的絕對擬合和絕對平滑之間起著平衡的作用，ε 越小，曲線的平滑程度越強(qiáng)[8]。

1.4 PID 控制原理 PID 控制是一種線性控制器，將系統(tǒng)的輸出偏差進(jìn)行線性組合構(gòu)成控制量，對系統(tǒng)進(jìn)行反饋調(diào)節(jié)，原理框圖如圖3 所示。常用的PID 算法模型包括增量式算法、位置式算法、微分先行法三種。本次銣鐘馴服采用增量式算法對銣鐘控制電壓進(jìn)行調(diào)節(jié)。

圖3 PID 控制原理框圖

PID 算法模型：

u（k）為控制器的輸出值，e（k）為控制器輸入值與設(shè)定值之間的誤差。KP為比例系數(shù)，KI為積分時(shí)間常數(shù)，KD為微分時(shí)間常數(shù)。用式（14）減去式（15），可得：

式（16）稱為增量式PID 算法，Δu（k）為控制器輸出控制量的變化。

對于銣鐘電路上DAC 模塊控制電壓的輸入數(shù)據(jù)，首先，設(shè)定KP、KI、KD的初值以及壓控電壓初始值。對于馴服銣鐘來說，e（k）為當(dāng)前測量得到的相差，即需獲取的測量輸入值。再計(jì)算增量誤差：

積分誤差：

微分誤差：

得到本次PID 的輸出值：

得出更新誤差值：

PID 參數(shù)互相影響，調(diào)試過程中可根據(jù)情況對積分系數(shù)、比例系數(shù)、微分系統(tǒng)進(jìn)行細(xì)調(diào)。經(jīng)實(shí)際試驗(yàn)與測試，設(shè)備PID 參數(shù)設(shè)置：比例系數(shù)P=0.1，積分系數(shù)I=0.15，微分系數(shù)D=0.1，即P=0.1，I=0.15，D=0.1。

1.5 銣鐘內(nèi)部原理 銣原子鐘內(nèi)部原理如圖4 所示，由量子部分和壓控晶體振蕩器組成。壓控晶振的頻率經(jīng)過倍頻和頻率合成器送到量子系統(tǒng)與銣原子躍遷頻率進(jìn)行比較，并將誤差信號送回壓控晶振，使其鎖定在銣原子特定的能級躍遷至對應(yīng)的頻率上[9]。

圖4 銣鐘內(nèi)部原理圖

2 銣鐘馴服模塊硬件組成

北斗接收機(jī)對接收的衛(wèi)星信號進(jìn)行解碼、濾波，輸出標(biāo)準(zhǔn)1PPS 信號至FPGA。FPGA 測量出標(biāo)準(zhǔn)1PPS 與銣原子鐘分頻產(chǎn)生的1PPS 之間的差值，并將濾波處理的真實(shí)差值發(fā)送給MCU，MCU 通過控制算法得到銣原子鐘控制電壓的控制量，并發(fā)送給高精度D/A 轉(zhuǎn)換器，轉(zhuǎn)換得電壓值發(fā)送至銣原子鐘，繼而調(diào)節(jié)銣鐘輸出頻率，使其產(chǎn)生的1PPS 信號與北斗衛(wèi)星信號的1PPS 信號誤差不超過一個(gè)時(shí)間周期，實(shí)現(xiàn)對銣鐘頻率的馴服。

圖5 銣鐘馴服原理圖

3 測試數(shù)據(jù)分析

本守時(shí)設(shè)備以北斗1PPS 信號作為頻標(biāo)信號進(jìn)行比較和分析，對銣原子鐘分頻產(chǎn)生的1PPS 信號進(jìn)行準(zhǔn)確度分析和對比。圖6 為kalman+vondrak 濾波后鐘差數(shù)據(jù)、單獨(dú)kalman 濾波后鐘差數(shù)據(jù)與原始鐘差數(shù)據(jù)對比。對比結(jié)果表明，kalman+vondrak 濾波對原始鐘差數(shù)據(jù)進(jìn)行了較好的平滑和預(yù)測，使其在衛(wèi)星信號失效情況下仍能較為準(zhǔn)確地對銣鐘進(jìn)行電壓控制。室溫環(huán)境下，測試本時(shí)統(tǒng)設(shè)備一天的守時(shí)精度。設(shè)備測試起始測量值208ns，設(shè)備測試24 小時(shí)之后的測量值為4.40us。根據(jù)測量出的數(shù)據(jù)可得，室溫環(huán)境下設(shè)備的守時(shí)精度為4.2us。

圖6 kalman+vondrak 濾波后鐘差數(shù)據(jù)、kalman濾波后鐘差數(shù)據(jù)、原始數(shù)據(jù)數(shù)據(jù)對比

4 結(jié)論

通過研究銣原子鐘馴服原理，提出了基于FPGA 的北斗時(shí)統(tǒng)設(shè)備銣鐘馴服實(shí)現(xiàn)方法。以北斗接收機(jī)輸出時(shí)間頻率為參照，測量銣鐘鐘差值，對鐘差值濾波前后數(shù)據(jù)進(jìn)行多次比對測試。測試數(shù)據(jù)表明，結(jié)合kalman 濾波和Vondrak 濾波算法對鐘差值數(shù)據(jù)的濾波效果，比單獨(dú)采用kalman 濾波算法更好，為銣原子鐘馴服提供了技術(shù)依據(jù)，可以較好地實(shí)現(xiàn)時(shí)間同步和提高守時(shí)設(shè)備精度。