劉 迪，武巖波，朱 敏，李 棟

(1. 中國科學院聲學研究所海洋聲學技術(shù)中心，北京100190；2.中國科學院大學，北京100190；3. 北京市海洋聲學裝備工程技術(shù)研究中心，北京100190；4. 中國科學院聲學研究所聲場聲信息國家重點實驗室，北京100190)

0 引 言

基于水聲通信技術(shù)的無線傳感器網(wǎng)絡(luò)作為海底觀測網(wǎng)的重要組成部分，在維護國家海洋權(quán)益、開發(fā)和保護海洋資源方面有重要意義，在民用、軍用海洋信息開發(fā)領(lǐng)域具有巨大的發(fā)展?jié)摿?，吸引了大量學者進行研究[1]。時鐘同步是水聲無線網(wǎng)各個節(jié)點協(xié)同觀測的關(guān)鍵技術(shù)，對于某些聲學或地震信號傳感器，在長期工作過程中需要保持較高的同步精度才能檢測和分析地震、海嘯等突發(fā)事件。另外，為節(jié)省傳感器節(jié)點能量，延長水下無線傳感網(wǎng)絡(luò)的生命周期，大多數(shù)工作節(jié)點采用睡眠-喚醒機制[2]，因此水下授時顯得更為重要。

一般來說，水下獲取精確時間基準的主要方式有：(1) 利用原子鐘等高穩(wěn)定度時鐘授時，但該方法主要有兩個缺陷。通常授時精度越高的系統(tǒng)，功耗、體積也越大。2016年，Woods Hole實驗室[3]實驗中使用的原子鐘能耗在幾十毫瓦量級，與普通石英晶振微瓦量級能耗間存在較大差異，而現(xiàn)階段水下固定平臺主要依靠自身攜帶的電池供電，能量非常有限，能耗是系統(tǒng)設(shè)計要考慮的首要問題。另外這種方式的時間累積誤差也是不斷積累的，長期工作也必須通過其他方式校準。(2) 水下聲波協(xié)議授時，這種方式通常通過已知精確時間的節(jié)點和需要進行時間校準的傳感器節(jié)點進行水聲信息交互，進而實現(xiàn)授時功能。水下聲波協(xié)議授時是本文研究的主要內(nèi)容。

傳感器節(jié)點使用晶振產(chǎn)生本地時鐘。在實際使用中，出于對降低成本的考慮，普遍使用頻率準確度低和穩(wěn)定性差的廉價晶體振蕩器[4]。因制造工藝、溫度等因素的影響，不同傳感器節(jié)點晶振的輸出頻率會不同。節(jié)點時鐘輸出頻率的差異稱為時鐘漂移率。若兩時鐘之間的時鐘漂移率為α，意味著記錄相同一段時間T，第一個時鐘相對于第二個時鐘記錄的時間的偏差為αT。時鐘漂移率的存在會導致本地時鐘之間產(chǎn)生累積誤差，稱為時鐘偏差[5]。這就要求授時協(xié)議要同時獲得時鐘漂移率和時鐘偏差并進行補償。

現(xiàn)已有多種針對陸地無線傳感網(wǎng)絡(luò)的授時算法[6-9]。其中參考廣播時鐘同步協(xié)議(Reference Broadcast Synchronization, RBS)算法[6]、洪泛時間同步協(xié)議(Flooding Time Synchronization Protocol,FTSP)算法[7]、無線精確時間同步協(xié)議(Wireless Precision Time Protocol, WPTP)算法[8]認為信號的發(fā)送和接收可以瞬時完成，(Timing-sync Protocol for Sensor Networks, TPSN)[9]沒有考慮同步消息交換過程中時鐘漂移率的影響，而水聲通信帶寬窄、通信比特速率低[10]、傳播時延大，這些特點使許多用于地上傳感器網(wǎng)絡(luò)的授時協(xié)議不能直接應(yīng)用于水下。

針對水下傳感網(wǎng)絡(luò)授時問題，研究者們做了一系列的研究，并取得了一些成果。2006年，Syed等[11]提出高時延時間同步(Time Synchronization for High Latency, TSHL)水聲網(wǎng)絡(luò)授時協(xié)議，第一個考慮了時鐘漂移率在高時延授時過程中對授時誤差的影響。2009年，Tian等[12]提出的Tri-message授時協(xié)議在TSHL算法開銷方面進行了改進，通過3次時間戳交換過程修正時鐘漂移率和時鐘偏差。2012年，美國海軍研究生院在TSHL授時協(xié)議的基礎(chǔ)上考慮了信道多徑對授時精度造成的影響[13]，并對授時過程進行了湖試實驗，實驗中時鐘漂移率最小估計誤差為 9×10-6s。以上研究中，水聲授時協(xié)議的設(shè)計都基于授時過程中時鐘漂移率為常數(shù)的假設(shè)。2016年，德國Evologics公司研制具有授時功能的S2CR Series Modem，并用其在地中海進行了授時實驗[14]。實驗中接收節(jié)點與發(fā)送節(jié)點之間的時間戳交換過程基于雙程信息交換模型，接收節(jié)點計算并記錄本地時鐘偏差隨時間的變化情況。該水聲授時實驗進行的時間很短(5～10 min)，并沒有對時鐘漂移率進行修正，也未進行時鐘漂移率的長期變化特性的觀察和研究。然而，在實際應(yīng)用中，溫度是影響時鐘輸出頻率最明顯的因素。傳感器節(jié)點若長時間運行在水下，時鐘漂移率受周圍環(huán)境的影響會產(chǎn)生變化。Sathiya等[15]為觀察時鐘漂移率變化特性與時間和溫度的關(guān)系進行了實驗并得出結(jié)論：各個節(jié)點的時鐘受到外部條件的影響，時鐘漂移率隨溫度和時間的變化呈現(xiàn)非線性特征。實時估計和補償時鐘漂移率，能提高授時精度并減少授時次數(shù)。然而，時鐘漂移率的測量會受到噪聲的影響，同時消息包的傳輸延時也具有不確定性，這些因素都會導致待同步節(jié)點無法對標準時鐘節(jié)點進行準確測量和估計。

卡爾曼(Kalman)濾波算法可以較好地解決實時狀態(tài)估計問題，但常規(guī)Kalman算法需要根據(jù)經(jīng)驗或先驗知識設(shè)定濾波參數(shù)。Sage-Husa算法在Kalman濾波算法的基礎(chǔ)上，可實時更新濾波參數(shù)來提高估計精度[16]。文獻[17]中采用一階線性模型對時鐘漂移率進行建模，并利用Kalman濾波對時鐘狀態(tài)進行估計和糾正，但該文獻未對變時鐘漂移率進行分析和研究。Yang等[18]提出使用交互式多模型Kalman濾波對時鐘漂移率進行估計，但處理的是高斯白噪聲數(shù)據(jù)，濾波參數(shù)固定，且時鐘漂移率的變化模型的提出也沒有實驗數(shù)據(jù)支撐。據(jù)我們所知，目前在水聲通信授時領(lǐng)域，沒有相關(guān)研究考慮長時時鐘漂移率變化對授時誤差的影響。

針對現(xiàn)有水聲授時協(xié)議未考慮授時后時鐘漂移率變化的問題，本文在對時鐘模型深入分析的基礎(chǔ)上，提出一種適用于水聲傳感網(wǎng)絡(luò)授時的時鐘漂移率估計方案。為了找到合適的時鐘漂移率變化模型，本文首先根據(jù)海試實驗數(shù)據(jù)確定時鐘漂移率的漂移范圍和趨勢，使用多個模型描述其變化行為；其次，在使用現(xiàn)有授時協(xié)議仿真水聲通信過程獲得時鐘漂移率的基礎(chǔ)上，使用交互式多模型 Kalman濾波方法對變時鐘狀態(tài)向量進行估計，并在濾波過程中采用 Sage-Husa自適應(yīng)方法動態(tài)確定濾波參數(shù)，提高估計精度；最后，通過仿真實驗驗證了該授時方法的有效性。

1 授時過程基本模型

本文假設(shè)場景中有一個標準的時鐘節(jié)點(即已同步節(jié)點)A和一個為待同步的本地時鐘節(jié)點B。假設(shè)t代表標準時間，T代表節(jié)點B的本地時間，α和b分別代表本地時鐘的時鐘漂移率和時鐘偏差。節(jié)點的本地時間和標準時間之間的關(guān)系可以表示為

時鐘同步的基本思想是節(jié)點間通過數(shù)據(jù)交換的方式，估計時鐘漂移率和偏差，修正本地時鐘，使得本地時鐘與標準時鐘保持同步[19]。

兩節(jié)點間數(shù)據(jù)交換示意圖如圖1所示。B發(fā)送同步消息，該消息包含其發(fā)送時間T1，經(jīng)過時間Δt后，信息到達A，A記錄下該信息到達時間t2；兩節(jié)點交換收發(fā)角色，A向B發(fā)送同步信息，信息中包含其發(fā)送時間t3和上次接收時間t2，經(jīng)過時間Δt后，信息達到B，B記錄下該信息到達時間T4。這樣，B共獲得4個時間戳，時鐘偏差計算可得：

由式(3)可知，授時誤差和時鐘漂移率及數(shù)據(jù)包交換時間成比例增長，數(shù)據(jù)包交換時間越長，由時鐘漂移率α導致的授時誤差e越大。若能較好地對時鐘漂移率進行補償，就能很大程度提高授時精度。

圖1 兩節(jié)點間數(shù)據(jù)交換示意圖Fig.1 Schematic diagram of data exchange between two nodes

2 基于現(xiàn)有水聲授時協(xié)議的時鐘漂移率計算方法

現(xiàn)階段，水聲通信領(lǐng)域較為經(jīng)典的授時協(xié)議主要包括 TSHL[11]和 Tri-message[12]等。其中，TSHL采用 25個時間戳消息線性擬合的方式估計時鐘漂移率，估計精度高，但耗能高。Tri-message采用3次消息交換估計時鐘漂移率，該方式估計精度較TSHL低，但耗能低。本文提出了組合通信方案，在保證一定精度的同時，多數(shù)情況下采用兩種授時方式組合計算時鐘漂移率。本文中采用兩種方式組合使用的方式，示意圖如圖2所示。

圖2 時鐘漂移率組合方式計算方法示意圖Fig.2 Schematic diagram of combined clock skew calculation method

時鐘漂移率具體計算方式如下：

(1) TSHL方法：標準節(jié)點A向待同步節(jié)點B連續(xù)發(fā)送M個標準時間信息(k表示第k次TSHL授時過程，k=1 , 2,??)，節(jié)點 B 根據(jù)接收到的M個消息時的本地時間戳和相應(yīng)的標準時間戳，在以下數(shù)據(jù)點進行線性擬合：，得到本地時鐘B的時鐘漂移率。

(2) Tri-message方法：采用3次信息交換，節(jié)點B獲得6個時間戳后，用式(4)估計本地的時鐘漂移率：

式中：l表示第l次Tri-message授時過程。本文假設(shè)使用 TSHL方式估計時鐘漂移率的時間間隔為I1，使用Tri-message的時間間隔為I2，I1、I2的單位為h。TSHL與Tri-message方式交替使用，具體通信過程圖示如圖 2所示。參數(shù)滿足M=25，2,??)。

對3種授時方法的精度和耗能方面進行分析：授時方法的耗能主要來自于對時鐘漂移率α的估計，α由接收節(jié)點對時間戳消息線性擬合得到。利用文獻[12]中的方法對 TSHL、Tri-messgae和本文采用的組合授時方法進行線性擬合算法復雜度分析，結(jié)果如表1所示。

表1 不同授時方法線性擬合計算復雜度Table 1 The complexities of linear fitting computation of different time synchronization methods

授時精度方面：由于TSHL授時方法發(fā)送的時間戳數(shù)量多，接收節(jié)點可以利用到的時間信息較為豐富，利用該授時方法得到的授時誤差較小，而Tri-message授時方法利用的時間戳數(shù)量少，該授時方法誤差較大，本文采用的組合應(yīng)用方法將兩種授時方案結(jié)合使用，在保證一定授時精度的同時減少耗能。對24個大小不同的時鐘漂移率α，3種不同授時方法的發(fā)送端分別向接收端發(fā)送不同數(shù)量時間戳，接收端通過對接收到的時間戳消息進行線性擬合對α進行估計，每次估計出的α值是100次實驗的平均值。表2是通過上述3種方法計算的授時方法精度誤差均值對比表。

表2 不同授時方法授時誤差均值比較Table 2 The average errors of different time synchronization methods

由表1和表2可知，本文采用的組合授時方法在保證一定同步精度的同時減少了耗能。

3 變時鐘漂移率估計方法

通過使用第2節(jié)描述的方法獲得時鐘漂移率，本文采用文獻[18]中提出的估計方法，將其作為交互式多模型Kalman濾波的輸入。另外，為了更準確地追蹤時鐘漂移率的變化，提高授時精度，本文在濾波過程中采用Sage-Husa算法實時更新濾波參數(shù)。

時鐘漂移率隨時間和溫度變化的情況是復雜的，不會保持在同一種狀態(tài)。本文使用多個模型描述時鐘漂移率的變化。交互式多模型Kalman濾波算法中的多個濾波器對應(yīng)于不同的模型，各模型描述不同的時鐘漂移率的變化特性。設(shè)時鐘離散狀態(tài)向量c(k) =[θ(k)α(k)ρ(k)]T，其中，θ(k)表示時鐘偏差，α(k)表示時鐘漂移率，ρ(k)表示時鐘漂移率變化率；使用恒定加速模型表征時鐘漂移率的變化情況：

式中：τ(k)表示采樣間隔，ωα(k)表示噪聲量。將時鐘模型改寫為矩陣向量表達式：

式中：A為狀態(tài)轉(zhuǎn)移矩陣，ω為系統(tǒng)噪聲向量。定義z(k)=α(k)為觀測值，則觀測模型可寫為

式中：H為觀測矩陣，ν(k)為觀測噪聲。ω和ν分別表示系統(tǒng)噪聲和測量噪聲，其協(xié)方差矩陣分別為Q和R。并設(shè)模型間轉(zhuǎn)移矩陣為馬爾科夫矩陣。對具有n個模型的系統(tǒng)，系統(tǒng)的狀態(tài)方程和測量方程如式(8)和(9)所示：

使用交互式多模型估計時鐘漂移率的過程如圖3所示。

圖3 交互式多模型Kalman濾波結(jié)構(gòu)圖Fig.3 Structure diagram of interactive multi-model Kalman filter

交互式多模型估計時鐘漂移率主要分為3個步驟：輸入交互、Kalman濾波和輸出組合。

(1) 輸入交互。計算匹配模型mj(k)的混合初始條件。mj(k)表示模型j在k時刻的概率，pij表示模型i到模型j的狀態(tài)轉(zhuǎn)移矩陣。模型預測概率和模型轉(zhuǎn)移概率μi/j分別表示為

輸入交互后的混合輸入表示為

(2) Kalman濾波。對模型j，將及z(k)作為輸入進行Kalman濾波。在濾波過程中，z(k)是通過仿真現(xiàn)有水聲授時協(xié)議的時間戳交換過程計算得到，并不符合經(jīng)典 Kalman濾波假設(shè)的高斯的統(tǒng)計特性，且實際應(yīng)用系統(tǒng)中系統(tǒng)模型一般是復雜的，噪聲統(tǒng)計特性往往是未知而且時變的，先驗數(shù)據(jù)會因此失去意義，想獲得噪聲統(tǒng)計特性十分困難。于是，本文采用在工程上廣泛應(yīng)用的Sage-Husa方法[16]，根據(jù)新獲得的觀測數(shù)據(jù)，在線估計出系統(tǒng)噪聲特性，實時更新濾波參數(shù)Q和R，保證濾波穩(wěn)定性，提高濾波收斂速率和效果。更新參數(shù)過程如下：

式(13)～(15)中：b為遺忘因子；ε(k)為信息矩陣。

(3) 輸出組合?；谒迫缓瘮?shù)計算k時刻模型mj的動態(tài)權(quán)重Sj(k)，并由此獲得多模型濾波器的輸出。對m個模型，首先分別計算模型mj的更新概率：

則最終輸出的狀態(tài)估計可表示為

采用該濾波算法可以得到變時鐘漂移率模型的估計結(jié)果，用估計結(jié)果對時鐘漂移率進行修正和補償，之后可進行雙程信息交換，修正時鐘偏差。

4 仿真結(jié)果和分析

2013年7月，在海南某淺海海域采集了一天內(nèi)的溫度變化數(shù)據(jù)。由這些數(shù)據(jù)可得水下節(jié)點晶振輸出頻率隨時間的變化曲線，由溫度變化導致時鐘漂移率產(chǎn)生的變化最大可達10-5s左右，如圖4中虛線所示。圖4中階梯狀的變化曲線表示現(xiàn)有水聲授時協(xié)議對變時鐘漂移率的處理：每一次授時后不考慮時鐘漂移率的變化(假設(shè)一天內(nèi)進行了 6次授時過程)。這會導致在時鐘漂移率變化程度較為明顯時，授時誤差大幅增加。極端情況下，一天時間內(nèi)由時鐘漂移率變化造成的授時誤差可達20 s，整個時鐘無法使用。

圖4 時鐘漂移率的海上試驗數(shù)據(jù)推算結(jié)果(虛線)及現(xiàn)有授時模型的處理結(jié)果(實線)Fig.4 The clock skew calculated from sea test data (dashed line) and the clock skew processed by existing time synchronization model (real line)

將根據(jù)實時溫度推算的時鐘漂移率變化曲線作為參考曲線，本節(jié)對提出的變時鐘漂移率估計方法進行了模擬仿真分析。估計方法使用 Sage-Husa自適應(yīng)Kalman濾波對數(shù)據(jù)進行處理，可根據(jù)輸入的數(shù)據(jù)實時估計系統(tǒng)的過程噪聲和測量噪聲協(xié)方差矩陣，追蹤輸入數(shù)據(jù)的變化并對濾波器參數(shù)進行修正，以此提高濾波效果。另外，用方法1表示基于現(xiàn)有水聲授時協(xié)議的時鐘漂移率估計方法，具體計算方案與參數(shù)設(shè)置詳見第2節(jié)。方法2表示在通過現(xiàn)有水聲授時協(xié)議估計時鐘漂移率的基礎(chǔ)上，使用不調(diào)整濾波參數(shù)的多模型Kalman濾波時鐘漂移率估計方法。本節(jié)對提出的估計方法與方法1和方法2進行了對比實驗。

圖 5和圖 6分別表示使用參數(shù)不變的標準Kalman濾波方法(簡稱：標準Kalman濾波)和Sage-Husa自適應(yīng)濾波方法的濾波結(jié)果對比圖和累計誤差對比圖。設(shè)時鐘漂移率真實變化模型如圖5中紅色實線所示，對該模型的觀察時間為200 s，前100 s模型噪聲為高斯白噪聲，后100 s用混合高斯噪聲仿真濾波的Q、R參數(shù)的改變。之后分別用參數(shù)不變的標準 Kalman濾波和 Sage-Husa 濾波方法對該模型進行濾波處理。

圖5 兩種濾波方法估計結(jié)果對比圖Fig.5 Comparison of estimation results of two filtering methods

圖6 兩種濾波方法累積誤差對比圖Fig.6 Comparison of accumulation results of two filtering methods

由圖5和圖6可以看出，使用參數(shù)固定不變的標準Kalman方法在噪聲發(fā)生變化的情況下，濾波效果不如濾波參數(shù)可以自適應(yīng)調(diào)節(jié)的濾波方法。

之后，對實驗數(shù)據(jù)進行分析、處理及仿真。仿真實驗中，設(shè)待同步節(jié)點初始時鐘偏差10 μs，時間戳接收抖動分布由均值為 0、標準差為 15 μs的高斯分布引入，時鐘的仿真粒度為1 μs，標準節(jié)點發(fā)送時間戳的時間間隔為1 s。

圖7顯示的是一天內(nèi)，采用本文授時算法、方法1和方法2對變時鐘漂移率的估計情況。由圖7可以看出，在0～4 h、8～12 h、16～20 h等時鐘漂移率出現(xiàn)較大程度變化的時間段，可以對其進行較好的估計和跟蹤，本文提出的時鐘漂移率計算方法結(jié)果表現(xiàn)優(yōu)于方法1和方法2。方法1因能耗因素，在大多數(shù)情況下采用 Tri-message授時方法計算時鐘漂移率，精度不高，且該方法單次授時結(jié)束后，并未考慮授時后時鐘漂移率的變化情況；方法2假設(shè)進行Kalman濾波的濾波參數(shù)是固定的，這在已知噪聲分布的情況下是可以的，但本文授時方法的時鐘漂移率的值是通過仿真現(xiàn)有水聲授時協(xié)議過程得到的，噪聲分布并不是已知的，故方法2只能跟蹤時鐘漂移率變化的大致趨勢，對時鐘漂移率模型判斷不準確，授時誤差較大。圖8給出了3種時鐘漂移率估計方法累積授時誤差變化圖。這3種方法的授時誤差隨著時間的推移都會有不同程度的增加，采用本文提出的時鐘漂移率估計方案對其進行估計和修正，授時精度最高，一天內(nèi)現(xiàn)有授時協(xié)議累積的授時誤差大約為本文提出的估計方法的1.5倍。圖9給出了3種時鐘漂移率估計方法結(jié)果的均方誤差變化圖。本文使用的估計方法與其他兩種授時方法相比，誤差增加速度最小，一天內(nèi)的授時均方誤差可由3.0×10-9s2降低到5×10-10s2。

圖7 3種授時方法估計時鐘漂移率的對比圖Fig.7 Comparison between the clock skews estimated by three time synchronization methods

圖8 3種授時方法累積授時誤差對比圖Fig.8 Comparison between the accumulated timing errors of three time synchronization methods

圖9 3種授時方法的均方誤差對比圖Fig.9 Comparison between the mean square errors of three time synchronization methods

綜上所述，該算法能實現(xiàn)對變時鐘漂移率進行預測和追蹤，有效地提高水聲通信節(jié)點的授時精度，具有較好的性能。

5 結(jié) 論

本文針對長時時鐘漂移率隨周圍環(huán)境產(chǎn)生變化的問題，提出了一種適用于水下傳感網(wǎng)絡(luò)的時鐘漂移率跟蹤算法。本文在水聲通信授時領(lǐng)域提出，在進行水聲授時的同時考慮長時變時鐘漂移率對授時精度產(chǎn)生的影響。使用多個模型描述時鐘漂移率的變化，在利用現(xiàn)有水聲授時協(xié)議計算獲得時鐘漂移率的基礎(chǔ)上，采用交互式多模型Kalman濾波對變時鐘漂移率進行估計，并在濾波過程中實時更新濾波參數(shù)，提高時鐘漂移率估計的準確性，進而提高水下節(jié)點的授時精度。仿真結(jié)果表明，該方法與現(xiàn)有的授時方法相比，累積授時誤差小，授時均方誤差可從 3.0×10-9s2降低到 5×10-10s2。