主從同步晶振頻率誤差的卡爾曼收斂算法

鄧長順，牛蘭杰，黃 輝

(西安機電信息技術研究所,陜西 西安 710065)

0 引言

彈藥組網技術具有廣闊的使用場景，戰場環境下高動態平臺和武器間的關系往往無法預先規劃，系統應具有自組織網絡特性[1]。無線自組網(wireless ad-hoc network)能在陸、海、空基作戰系統上應用，對于機動協同作戰有重要意義[2]。這使得高動態協同攻擊對無線組網時間同步精度控制提出了更高要求，時間一致性對提高信息時效性及協同毀傷效能至關重要，所以時間同步精度是組網的關鍵指標。

無線組網時間同步算法優化趨勢：1) 分析并降低同步消息交換過程中的各類時延，實現同步精度的提高；2) 引入分簇思想對網絡拓撲改進，在多跳網絡中有效克服誤差隨跳距增大而累計的缺陷且平衡網絡能量消耗，可上升至MAC層路由算法層面；3) 融入線性擬合等數學理論，對同步數據進行計算調整[3]。近期學者提出雙向交互同步算法(timing-sync protocol for sensor network，TPSN)，從路由角度入手軟件控制，對網絡拓撲進行層級劃分，逐層進行時間同步實現全局同步，是分簇思想的代表[4]；文獻[5]在增強無線傳感器網絡時間同步中提出了低功耗組網時統方法，單從低頻晶體降低能耗分析；文獻[6]在論文中研究衛星同步高穩晶振的影響，但對采樣噪聲考慮較少。

本文針對當前高動態無線平臺時間控制精度要求較高，現有時間同步算法忽略了無線信道數據包沖突導致的硬件物理時鐘延遲抖動及時間戳的采集等噪聲影響，提出基于主從同步晶振頻率誤差的卡爾曼收斂算法。

1 無線組網時間同步架構

1.1 無線組網時間同步層次劃分

時間同步是指在分布式系統中，物理上彼此獨立的節點為了保持時間一致性，通過無線通信協議(如時分多址TDMA)將本地時鐘調整到統一時間點的過程。無線自組網大多基于報文傳輸同步信息的協議都面臨傳輸過程中的不確定性，有些不確定時間誤差可能高于系統要求的同步精度。無線組網收發分解示意圖見圖1[7]。

圖1 同步數據收發分解示意圖Fig.1 Schematic diagram of synchronous data sending and receiving decomposition

信道接入時間是收發過程中傳輸時延中波動性最大的一部分，消息等待信道空閑所需時間，取決于MAC協議和網絡負載狀況。節點從應用層發出數據包，到接收方應用層收到數據包，每一層之間都存在著一定的噪聲和誤差，越靠近物理層輸出的時間戳標記誤差波動越小，通訊過程中協議設計收發數據包的雙方在MAC層標記時間戳[8]，然后交換各自時間戳后進行統計再根據同步算法調整本地時鐘，兩個節點完成一次“時間握手”。開放式通信系統互聯參考模型(open system interconnection)，是國際標準化組織(ISO)提出的一個使各種計算機在世界范圍內互連為網絡的標準框架，簡稱OSI[9]。主要包括七層，如圖2(a)，本文組網架構層次如圖2(b)所示。

圖2 網絡架構層次Fig.2 Architecture hierarchy of network

1.2 硬件時鐘物理建模

從物理角度來看，由于節點的石英晶體振蕩頻率不是恒定的，溫度、輻射等會造成節點的本地硬件模擬時鐘出現抖動和偏移，導致軟件時鐘發生頻率漂移和時鐘偏移，所以時間同步需將不同節點內的晶振頻率收斂到一個基準點。節點是通過晶體諧振器(如圖3所示)和振蕩器自身組裝的中斷進行計數來實現的，頻率可通過如溫補晶振、改變負載電容CL來調整頻率輸出。每秒產生中斷計數的值理想狀態是不變的，晶體初始定值和頻率誤差不同，長期穩定性不好。圖4為硬件時鐘過慢(或過快)和時鐘漂移與本地標準時鐘的差異。

圖3 晶體諧振器模型Fig.3 Crystal resonator model

圖4 時鐘誤差模型Fig.4 Clock error model

定義T為實際物理時間，T(t)為節點的本地時間。理想時鐘的時間變化率等于常數1。由于晶體振蕩器存在一定的漂移，可定義的時間變化率為f(t)=dC(t)/dt，節點i局部時間Ti(t)從物理時間t0到t的變化可用式(1)表示，ψi(t0)為節點i在物理時刻t0的初始相位值。

(1)

式(1)通過泰勒級數展開為以下方程：

Ti(t)=βi+αit+γit2+…，

(2)

式(2)中，βi和αi分別代表實際物理時間t的初始相位偏移和頻率偏移。若式(2)中二次項以上系數均為0，則可簡化為線性預測時鐘模型，近似為:

Ti(t)=βi+αit。

(3)

基于式(3)描述的線性時鐘模型，時鐘關系取決于時鐘頻率參數αi和初始相位參數βi。節點啟動時，參數βi被抑制短時間不變，而參數αi隨著時間變化而變化，這主要是由晶體振蕩器的漂移效用引起的。根據式(3)，節點A和B之間的時鐘關系可以表示為：

TB(t)=βAB+αAB·TA(t)，

(4)

式(4)中，βAB和αAB分別表示節點A相對節點B之間的相位偏移和頻率偏移。當βAB=0和αAB=1時，節點A和節點B可以實現時鐘同步。而當βAB≠0和αAB≠1時，這兩個參數隨時間的積累變化，會產生如表1的結果。

表1 時鐘參數變化影響Tab.1 Impact of clock parameter changes

1.3 主從同步周期時間步驟

基于主從同步周期的運行步驟如圖5所示。

圖5 子節點時鐘同步模型Fig.5 Child node clock synchronization model

基于上述模型可分析兩個節點之間的時延與頻偏參數：

(5)

(6)

式(5)、式(6)中，w表示時間偏移率，T′表示相位偏差，d為傳輸噪聲誤差，X為外部隨機噪聲時延。經多次同步過程后，若不考慮各類誤差，可以進一步估算k時刻兩個節點之間的時間偏移量T(k)與時間偏移率W(k)：

T(k)=[(Tc(2,k)-Tv(1,k))-(Tv(4,k)-TC(3,k))]/2，

(7)

W(k)=[T(k)-T(k-1)]/Δt。

(8)

2 主從同步晶振頻率誤差卡爾曼收斂 算法

線性濾波預測理論，是在線性狀態空間的基礎上對噪聲和觀測信號進行處理：系統內部擾動和觀測噪聲對它們的統計性質作適當假定，對含噪聲的觀測信號進行狀態觀測反饋。狀態反饋器工作模式：x為系統輸入值，Y為輸出的真實測量值，測量誤差會影響其準確性，y′為輸出變量的估計值。狀態變量的真實值通過狀態觀測器對估計值與系統測量值(含噪音)進行比較輸出后經回路K反饋，對未知狀態估計進行校正后得到一個可靠狀態估計值。

圖6 狀態觀測器Fig.6 Status observer

由于子節點初始化具有隨機性，時鐘初始相位是不同的，將節點相位之間的差值稱之為時鐘偏移(offset)，同時稱晶體振蕩器實際工作頻率與額定頻率之間存在的偏差為時鐘偏移率。同步過程中，節點的時鐘偏移率會隨著時間不斷累積，最終使節點間產生較大的時鐘偏移。因此需周期性地對節點進行校驗，保證時間同步誤差在可控范圍內。基于主從時間同步及對理論時鐘模型變量分析，根據卡爾曼濾波運動學原理，設計了圖7的算法流程圖，通過對之前的時間統計，各種噪音(頻率噪聲、隨機噪聲等)的估計輸入進行濾波，逐漸貼近真實時間戳輸出值，提高系統的時間同步精度。

圖7 主從同步晶振頻率誤差卡爾曼收斂算法流程圖Fig.7 Master-slave synchronous crystal oscillator frequency error Kalman convergence algorithm flowchart

2.1 時間同步狀態轉移

狀態轉移是根據觀測數據對隨機量進行定量推斷，動態行為的狀態估計能實現實時運行狀態的估計和預測。式(9)為通常狀態下的卡爾曼濾波狀態轉移方程，式(10)為參數協方差轉移矩陣。

(9)

P=APAT+Q。

(10)

在時間同步系統中，X為時延參數矩陣，含時間偏移以及時間偏移率，A和B是系統參數，對于多模型系統時為矩陣。H為參數測量轉換矩陣，系統過程Z和測量的噪聲被假設成高斯白噪聲，它們的協方差(covariance)分別是Q，R。對于主從時間同步系統，將式(5)與式(6)作差可得式(11)，推導得到時間同步狀態轉移方程(13)：

(11)

T=HX+S，

(12)

Tk=Tk-1+Δt·Wk-1，

(13)

Wk=Wk-1，

(14)

式(12)是節點時延的向量表達，T為時延參數，S為觀測噪聲如振蕩器內部噪聲誤差等。X選用主從時鐘的偏移量T和時鐘偏移率W為狀態變量，對時鐘進行數學建模，推導可得時間同步的狀態觀測轉移矩陣：

(15)

式(15)中，Δt為時間同步周期。

2.2 基于卡爾曼濾波的主從時間估計同步算法

分析可知，實際描述的系統傳輸狀態轉移延遲抖動較大，存在環境等觀測噪聲干擾，且時間戳獲取存在誤差，故采用卡爾曼濾波，估計這些誤差的影響后在同步過程中消除。該算法具體操作步驟為：預測過程，X(k)是k時刻下的系統狀態，A為系統參數觀測轉移矩陣，B為時間同步控制矩陣，u為同步校正矩陣，L為過程狀態噪聲。

X(k|k-1)=A·x(k-1)+B·u(k-1)+L。

(16)

最小預測協方差矩陣為P(k|k-1)，可以初始化為Q：

P(k|k-1)=AP(k-1)AT+Q

(17)

在卡爾曼濾波中的時間同步卡爾曼增益矩陣為：

K(k)=P(k|k-1)HT[HP(k|k-1)HT+R]-1

(18)

綜合考慮時間戳參數與噪聲影響，得到狀態修正式(19)。可利用估計得到的主從時鐘偏移量估計量對時鐘模型進行修正，從而提高時鐘同步精度。

X(k|k)=X(k|k-1)+K(k)[z(k)-HX(k|k-1)]。

(19)

P(k|k-1)是X(k|k-1)對應的協方差，AT表示A的轉置矩陣，Q是系統過程的協方差，H為觀測量轉換矩陣。得到k狀態下最優的估算值X(k|k)，但是為了要令卡爾曼濾波器不斷運行下去直到系統過程結束，還要更新k狀態下X(k|k)的協方差：

P(k)=(I-K(k)H)P(k|k-1)。

(20)

通過上述系統迭代公式對于時間戳估計誤差的兩個關鍵參數時間偏移和時間偏移率作出了跟蹤估計，使之與真實偏差收斂到最小。子節點根據估計的時間偏移和時間偏移率作出調整，并將新的時間參數輸入到系統中進行迭代，達到時間同步目的。

3 基于算法的仿真實驗驗證

為了驗證本算法對于節點時鐘同步精度的提高，分別對時間偏移和時間偏移率的跟蹤消除，根據第2章介紹的算法使用Matlab軟件進行仿真。本實驗采用卡爾曼濾波，綜合節點的物理時鐘特性，且考慮測量、過程噪聲等因素的影響來驗證本算法。

3.1 仿真實驗子節點參數

本仿真實驗的子節點參數設置：硬件時鐘理論誤差很小，大多數硬件時鐘的頻率都在10～50 MHz，所對應的理論時鐘精度是在納秒級別。但評價時鐘源性能的指標包括頻率準確度(在標稱電壓、標稱負載、基準溫度，晶體振蕩器的頻率相對于其規定標稱的最大允許偏差，公式換算：時鐘偏移率=頻率偏移/標稱頻率=(fmax-fmin)/f0)以及頻率穩定度(在規定溫度范圍內晶體振蕩器輸出頻率的最大變化量相對于溫度范圍內輸出頻率極值之和的允許頻偏值)[13]。溫補晶振及恒溫晶振是通過對溫度特性進行補償或者構建恒溫環境以分別達到±0.01及±0.001級別的頻率穩定度[14]。頻率準確度上，普通的晶體振蕩器一般在±50 ppm以內，即對于標稱頻率為10 MHz的石英晶體而言，其實際頻率與標稱頻率的差別不超過±500 Hz，即兩個標稱頻率相同的晶體振蕩器每秒最大可能時鐘偏移為100 μs。圖8為基于主從同步模型的子節點時鐘模型。

圖8 子節點時鐘同步模型Fig.8 Child node clock synchronization model

3.2 本文濾波估計算法仿真實驗

3.2.1算法有效性驗證

在基于卡爾曼濾波估計時間同步算法中，根據晶振參數與時鐘偏移的模型特性，將初始的時間偏移offset設為100 μs，將初始的時鐘偏移率drift設為50 ppm；采樣次數為100次，進行卡爾曼濾波時間估計同步次數為40次來進行仿真驗證，各類噪聲設定為高斯白噪聲。每進行一次濾波后，對從時鐘進行精度調制，再次測量后輸出。

由于延遲抖動以及獲取時間戳帶來的各類噪聲誤差，系統經過卡爾曼濾波的估計值會在真實值附近震蕩，存在一定的誤差，可利用估計得到的主從時鐘偏移估計量以及估計時鐘偏移率對子節點時鐘模型進行修正，從而提高時鐘同步精度。從仿真數據圖9—圖11看出，經過卡爾曼濾波約20次迭代后，時鐘偏差的差值可穩定控制在±3 μs以內，頻率偏移率經過約20次迭代后穩定控制在差值±3 ppm以內。由圖12統計同步后采樣點的偏差分布，可以看出經過40次迭代同步后偏移分布與偏移率分布大致符合均值為0的一個高斯分布，offset在0.259 3 μs、drift在0.018 8附近震蕩。綜合計算考慮，對于自組網時間同步精度可由百微秒級別提升至十微秒以內量級，且消除了一些噪聲。

圖9 時間偏移與時間偏移率第一次同步后的效果Fig.9 Time offset and time offset rate after the first synchronization

3.2.2算法可靠性分析

仿真基于算法模型增大時間偏移與時間偏移率后分別進行了運算分析，隨著運算次數的增多，增加到1 000 μs誤差級別也逐漸穩定在20～40多次以內迭代同步。濾波各步都在同一個運算級范圍內，說明多次運算不會嚴重影響組網同步信息處理的實時性；且控制迭代次數有限的情況下(20次)，逐步增大輸入偏差，其最終輸出偏差均值可穩定在5 μs以內(見圖12)，符合我們的同步精度提升要求。

圖10 時間偏移與時間偏移率隨同步次數的變化Fig.10 Time offset and time offset rate vary with the number of synchronization times

圖11 同步完成后時間偏移與時間 偏移率采樣點統計分布Fig.11 Statistical distribution of time offset and time offset rate at sampling points after synchronization

圖12 輸入偏差值經過有限次迭代后與輸出均值的關系Fig.12 The relationship of the input deviation value to the mean output deviation after a finite number of iterations

通過調整觀測噪聲與過程狀態噪聲，探尋算法的噪聲適應邊界條件：仿真中對觀測噪聲與過程狀態噪聲標準差設定從1逐步增大到5，同步100次后，最終偏移收斂但方差仍較大，如圖13所示，綜合造成的時間偏移仍有數十微秒，無法滿足所需同步精度提升量級。從數據側面也能看出觀測噪聲相比過程噪聲的影響更大，本算法對抑制噪聲有一定的能力，但當觀測噪聲與過程狀態噪聲影響過大時，未能達到預想的同步效果。

圖13 觀測噪聲與過程噪聲方差對offset 與drift方差影響Fig.13 Observe the variances of noise and process noise on the variances of offer and drift

綜上，將硬件晶振的物理模型轉化為數學模型，經過卡爾曼濾波并對同步過程中的各類噪聲加以分析，子節點時間戳統計量隨著多次迭代時間偏移量逐步降低。通過卡爾曼濾波消除了一些難以避免的由于無線信道數據包碰撞導致的各類噪聲與抖動誤差的影響，將組網的時統精度控制提高一個數量級。

4 結論

本文提出基于主從同步的卡爾曼濾波晶振頻率收斂算法。該算法利用卡爾曼濾波對于主從節點無線組網進行時間戳統計估計，并考慮過程及觀測噪聲影響，對物理硬件時鐘晶振頻率參數處理，分析時間偏移與時間偏移率參數隨多次仿真同步的收斂規律。仿真實驗驗證結果表明，經過卡爾曼濾波時間估計后，在有限的同步次數下，時間同步精度可從百微秒量級提升至十微秒以內，同步精度方差穩定在1以內，有效降低了晶振頻率偏差對時間同步精度的影響。本方法在無線組網硬件設計實現還具有很大的實踐空間，可以進一步調試濾波器對于多節點組網時延精度以及功耗的影響。