顧及衛(wèi)星間相關性的Kalman短期鐘差預報

楊 旭，王潛心，呂偉才，鐘東升，譚福臨

（1.安徽理工大學礦區(qū)環(huán)境與災害協(xié)同監(jiān)測煤炭行業(yè)工程研究中心，安徽 淮南 232001；

2.安徽理工大學 空間信息與測繪工程學院，安徽 淮南 232001；

3.安徽理工大學礦山采動災害空天地協(xié)同監(jiān)測與預警安徽普通高校重點實驗室，安徽 淮南 232001；4.中國礦業(yè)大學 環(huán)境與測繪學院，江蘇 徐州 221116）

0 引言

全球衛(wèi)星導航系統(tǒng)（global navigation satellite system, GNSS）衛(wèi)星鐘差短期預報（一天以內）是GNSS應用中的一個重要問題，對實時精密單點定位技術（precise point positioning，PPP）、GNSS系統(tǒng)間時間同步等方面研究具有重要作用[1]?，F(xiàn)有衛(wèi)星鐘差短期預報研究已對預報模型本身研究較多，對建模策略、方式正在進一步拓展。文獻[2]利用卡爾曼（Kalman）濾波、遞推最小二乘和灰色模型GM(1,1)進行了一天內的衛(wèi)星鐘差短期預報，結果表明 Kalman濾波預報精度最高。文獻[3]在對比線性模型、二次多項式模型、灰色模型以及Kalman濾波模型預報精度后，也得出Kalman濾波模型的短期預報精度最高。文獻[4]利用哈達瑪（Hadamard）總方差，計算了Kalman濾波預報模型觀測與狀態(tài)噪聲協(xié)方差陣，得到的全球定位系統(tǒng)（global positioning system, GPS）Block IIR銣原子鐘（Rb）的短期預報結果，驗證了該模型的正確性。文獻[5]提出將預報殘差周期項應用到北斗衛(wèi)星導航系統(tǒng)（BeiDou navigation satellite system,BDS）衛(wèi)星超短期鐘差預報模型中，利用德國地學研究中心（Deutsches GeoForschungsZentrum, GFZ）的事后多GNSS精密鐘差（GFZ multi-GNSS precise clock, GBM）進行仿真，其結果表明，三類BDS衛(wèi)星在1 h的預報時間內，其精度整體優(yōu)于0.8 ns，2和 3 h的預報精度較傳統(tǒng)線性模型提高了4%~20%。文獻[6]利用星地雙向時間比對（twoway time transfer, TWTT）與軌道確定與時間同步（orbit determination and time synchronization,ODTS）獲得的鐘差數(shù)據，對BDS衛(wèi)星鐘的隨機特征、周期特征進行了分析，設計了顧及該特征的Kalman濾波算法，實現(xiàn)了TWTT和ODTS鐘差數(shù)據的融合，驗證了該算法應用到衛(wèi)星鐘差短期預報的正確性。文獻[7]提出一種多項式同時結合周期項與反向傳播（back propagation, BP）神經網絡的BDS鐘差預報模型，對國際GNSS監(jiān)測評估系統(tǒng)（international GNSS monitoring and assessment system，iGMAS）與 GFZ的超快速鐘差數(shù)據進行了預報，3、6、12及24 h的預報精度較原產品提高了10.34% ~26.14%。文獻[8-9]研究了法國空間研究中心（Centre National D'Etudes Spatiales,CNES）的狀態(tài)域空間表達（state space representation，SSR）改正數(shù)更新頻率及中斷時間間隔對多系統(tǒng)實時 PPP的影響，結果表明在 PPP服務端，SSR更新頻率低于60 s；在PPP用戶端，SSR改正數(shù)中斷超過360 s時，均會明顯降低PPP收斂時間和定位精度。

由以上研究結果可以看出，對于衛(wèi)星鐘差短期預報，在單一預報模型中，Kalman濾波模型較其他模型具有一定優(yōu)勢；組合模型或單一模型的改進一般從顧及衛(wèi)星鐘物理、周期、隨機特征等方面考慮；衛(wèi)星鐘差實時預報對實時 PPP較為關鍵。但這些模型仍存在一定的局限性：部分單一模型不適合實時預報；組合模型存在模型復雜度高且執(zhí)行效率低等問題；僅針對單一衛(wèi)星鐘差序列進行建模預報時，沒有考慮各衛(wèi)星之間的相關性。針對上述不足，本文進一步拓展Kalman濾波算法在衛(wèi)星鐘差中的建模方式與策略，以此來提高鐘差預報性能；擬對GBM產品采用一種顧及衛(wèi)星間相關性的Kalman濾波預報算法，并結合阿倫（Allan）方差和Hadamard方差來進行四大衛(wèi)星導航系統(tǒng)的衛(wèi)星鐘差超短期（3 h以內）或短期預報，并與傳統(tǒng)的Kalman濾波短期預報方法進行對比分析。

1 衛(wèi)星鐘差預報的Kalman算法模型

原子鐘隨機模型與冪律譜噪聲模型符合較好，由5種獨立的噪聲構成[10]，該噪聲的表現(xiàn)與取樣時間有關：在衛(wèi)星鐘中短期取樣時間內（15 d以內），衛(wèi)星鐘一般只具有1~3種噪聲[11]。本文基于三狀態(tài)（調相隨機游走，調頻隨機游走，調頻隨機奔跑噪聲）Kalman濾波模型研究衛(wèi)星鐘差短期預報。

1.1 鐘差模型

對第k歷元、第j顆衛(wèi)星建立的原子鐘鐘差Kalman濾波模型為

式中：τ為預報時間；x(t)為t時刻的鐘差；y(t)為x(t)的時間導數(shù)；z(t)為y(t)的時間導數(shù)；Δx、Δy、Δz分別為獨立于x(t)、y(t)、z(t)的隨機模型誤差，其過程噪聲分別為q1、q2、q3；wk為模型狀態(tài)噪聲，其協(xié)方差矩陣為Qw；vk為模型的觀測噪聲，其噪聲方差為Rk。

1.2 噪聲方差陣推導

通過Allan方差或Hadamard方差可求取式（3）中與衛(wèi)星鐘頻率穩(wěn)定度相關的Qw。

1）Allan/Hadamard方差計算。銫原子（Cs）鐘的頻率穩(wěn)定性常用Allan方差來描述[12-13]，對于相位偏差數(shù)據{xn,n= 1 ,2,… ,N}，N為采樣個數(shù)，Allan方差為

式中：τ=mτ0為平滑時間；τ0為相鄰鐘差相位數(shù)據采樣間隔；m為平滑因子，一般取1 ≤m≤ in t[(N′?1)/3]；N′為按照τ大小間隔的相位數(shù)據個數(shù)。

Rb鐘頻漂項較為明顯，常用 Hadamard方差描述其頻率穩(wěn)定性，利用相位偏差數(shù)據計算的Hadamard 方差[14]為

同時考慮式（3）中的3種鐘差調頻噪聲與調相白噪聲時，式（5）與式（6）中的方差分別為

式中，q0為調相白噪聲的測量噪聲參數(shù)。

2）噪聲矩陣求解。式（5）、式（6）中，τ取不同值，可得不同的Allan、Hadamard方差值，再聯(lián)立式（7）、式（8），q1、q2、q3和q0可利用最小二乘算法進行估計，進而可確定式（3）、式（4）中的Qw、Rk。

1.3 顧及衛(wèi)星間相關性的Kalman濾波衛(wèi)星鐘差短期預報模型

目前，基于Kalman濾波的鐘差建模與預報多數(shù)情況下都是基于一顆衛(wèi)星進行的，鮮有針對多顆衛(wèi)星進行整體建模與預報。同時，衛(wèi)星軌道等誤差會影響多衛(wèi)星導航融合估計的衛(wèi)星鐘差，衛(wèi)星的空間環(huán)境也會干擾星載原子鐘的運行，這些因素均會使估計的衛(wèi)星鐘差存在一定的相關性。本節(jié)將推導顧及這種相關性影響的情況下，Kalman濾波鐘差預報模型。

假設第k歷元共有n顆衛(wèi)星，顧及衛(wèi)星間相關性的Kalman濾波鐘差預報模型可表達為

式中：sk，ηk分別為模型狀態(tài)噪聲與觀測噪聲，對應的協(xié)方差陣和噪聲方差陣分別為w′Q、Rk，由式（1）至式（4）可知

式中：ijσ、r分別為第i顆衛(wèi)星和第j顆衛(wèi)星之間的協(xié)方差及相關系數(shù)。

r可由皮爾遜相關系數(shù)（pearson correlation coefficient）公式計算[15]，其計算方法為

式中：d為鐘差預報采用的基礎數(shù)據個數(shù)；gi、ui分別為該基礎數(shù)據利用二次多項式擬合后的兩衛(wèi)星擬合殘差的時間序列；分別為gi、ui的均值；r為兩時間序列間相關關系密切程度的強弱，r為0時表示無線性相關關系，r越大，相關程度越高。

以上過程為顧及衛(wèi)星間相關性的Kalman濾波模型，考慮相關性的二次多項式、二次多項式加周期項的預報模型均與以上步驟類似，此處未列出。

2 多系統(tǒng)鐘差數(shù)據特征分析

本文利用事后鐘差，設計顧及衛(wèi)星鐘差隨機特征與周期特征的預報方法進行預報實驗，以此來檢驗文中所提新方法的可靠性和有效性。本文選取了2019年年積日第14—40天共27天GBM鐘差產品（30 s采樣間隔）、共84顆衛(wèi)星參與了數(shù)據特征分析，這些導航衛(wèi)星分別為C01、C02、C03、C04、C05、C06、C07、C08、C09、C10、C11、C12、C13、C14等 14顆 BDS衛(wèi)星；G01、G02、G03、G04、G05、G06、G07、G08、G09、G10、G11、G12、G13、G14、G15、G16、G17、G18、G19、G20、G21、G22、G23、G24、G25、G26、G27、G28、G29、G30、G31、G32等 32顆 GPS衛(wèi)星；R01、R02、R03、R04、R05、R07、R08、R09、R10、R11、R13、R14、R15、R16、R17、R18、R19、R20、R21、R22、R23、R24等22顆俄羅斯格洛納斯衛(wèi)星導航系統(tǒng)（global navigation satellite system,GLONASS）衛(wèi)星；E01、E02、E03、E04、E05、E07、E08、E09、E11、E12、E14、E18、E19、E24、E26、E30等 16顆歐盟伽利略衛(wèi)星導航系統(tǒng)（Galileo navigation satellite system, Galileo）衛(wèi)星。該 84顆衛(wèi)星軌道類型及衛(wèi)星鐘類型參見文獻[16]。以C13、G01、R01及E01為參考衛(wèi)星，利用二次差法計算的標準偏差（standard deviation, STD）進行鐘差精度統(tǒng)計[17]。利用文獻[5]衛(wèi)星鐘單天解基準處理方法，消除GBM天與天產品之間存在的系統(tǒng)性鐘差基準偏差，同時利用文獻[18]改進的巴爾達（Baarda）粗差探測法，對該鐘差序列進行數(shù)據預處理。

2.1 頻率穩(wěn)定性分析

利用 Allan方差、Hadamard方差計算的連續(xù)27天四系統(tǒng)按照衛(wèi)星軌道類型分類的衛(wèi)星鐘千秒穩(wěn)，萬秒穩(wěn)，天穩(wěn)定性結果及衛(wèi)星鐘的噪聲水平q0，q1，q2，q3統(tǒng)計結果見圖1。

由圖1可知：

1）Allan方差和Hadamard方差計算的頻率穩(wěn)定性指標基本一致，GLONASS衛(wèi)星的頻率穩(wěn)定性較差，GPS II-F衛(wèi)星鐘、Galileo完全運行能力（full operational capability, FOC）衛(wèi)星鐘、Galileo在軌驗證（design and on-orbit verification, IOV）衛(wèi)星鐘的頻率穩(wěn)定性較高。

2）Allan方差和Hadamard方差計算的噪聲水平基本一致，各衛(wèi)星系統(tǒng)的原子鐘噪聲水平變化規(guī)律和其頻率穩(wěn)定性變化規(guī)律一致。

2.2 周期特性分析

由于衛(wèi)星軌道和衛(wèi)星鐘差高度相關，使得鐘差數(shù)據有一定的周期性。去除上述27天鐘差數(shù)據的二次趨勢項之后，得到相應的鐘差殘差數(shù)據，利用快速傅里葉變換（fast Fourier transform, FFT）方法獲取鐘差的主要頻率來分析其周期特性。圖2是GBM產品數(shù)據完整率較高的C03，C07，C13，G01，R15，E01衛(wèi)星鐘差主頻率分析結果，橫軸為一天中出現(xiàn)某個周期的次數(shù)。表1是連續(xù) 27天 GBM鐘差產品兩個主周期的統(tǒng)計結果。

圖2 GBM產品C03，C07，C13，G01，R15，E01衛(wèi)星振幅頻譜

表1 GBM產品主周期、第二周期統(tǒng)計結果 單位：h

（續(xù)）

由圖2及表1可知：

GBM鐘差有一定的周期性，且隨著周期的縮短，振幅有逐漸降低的趨勢。對于BDS，GBM產品中有較明顯的12、24 h周期；對于GPS，GBM產品中有較明顯的6、12 h周期；對于GLONASS，GBM產品中有較明顯的11 h周期；對于Galileo，GBM產品中有較明顯的7、14 h周期。

2）GBM鐘差的周期性與衛(wèi)星系統(tǒng)有關；對于每個衛(wèi)星系統(tǒng)，衛(wèi)星的周期性與衛(wèi)星軌道類型的相關性不強。衛(wèi)星鐘差周期不一致的原因有：①軌道周期的耦合影響；②晝夜環(huán)境變化；③擬合殘差存在誤差；④鐘差數(shù)據部分缺失。

3）分析出周期項后，就可以利用考慮周期影響的Kalman濾波，附加周期項的二次多項式等模型進行鐘差預報。

3 多系統(tǒng)鐘差短期預報分析

為了驗證本文所提顧及衛(wèi)星間相關性方法的有效性，與傳統(tǒng)Kalman濾波、二次多項式、附周期項的二次多項式預報模型進行預報性能比較，設計了4種預報實驗：1）方案1，Kalman濾波模型預報；2）方案2，顧及相關性的Kalman濾波模型預報；3）方案3，附加周期項二次多項式預報；4）方案 4，顧及相關性的附加周期項二次多項式預報。Cs鐘，Rb鐘、氫原子（Ph）鐘分別采用Allan方差、Hadamard方差計算 Kalman濾波模型中的狀態(tài)噪聲。

利用前一天鐘差數(shù)據（30 s間隔）求取Kalman濾波模型狀態(tài)噪聲、衛(wèi)星間的相關性、二次多項式擬合系數(shù)與2個主要周期之后，預報后一天1、3、6、12及24 h的數(shù)據（30 s間隔）；鐘差精度統(tǒng)計時，采用與自身GBM產品比較。

圖3為4種方案下連續(xù)26天（2019年年積日第15—40天）四系統(tǒng)衛(wèi)星鐘差5種預報時長下的預報精度，表2為不同預報方案精度統(tǒng)計與對比分析結果。

表2 GBM產品不同方案四系統(tǒng)衛(wèi)星鐘差預報精度及方案間精度對比統(tǒng)計結果

圖3 GBM產品方案1（圓點）、2（三角形）、3（菱形）、4（正方形）四系統(tǒng)衛(wèi)星鐘差1、3、6、12及24 h預報精度

（續(xù)）

由圖3與表2統(tǒng)計結果可知：

衛(wèi)星系統(tǒng)不同，預報精度不同：對于6 h以內的預報，精度最低的是 GLONASS，最高的是Galileo，GPS、BDS次之；對于超過12 h的預報，BDS、GLONASS的預報精度相對較差，次于Galileo、GPS；在BDS，GPS中，不同衛(wèi)星軌道、衛(wèi)星鐘類型的衛(wèi)星鐘差預報精度不盡相同，BDS中的部分傾斜地球同步軌道（inclined geosynchronous orbits, IGSO）衛(wèi)星預報精度優(yōu)于地球靜止軌道（geostationary Earth orbit, GEO）衛(wèi)星；多數(shù)GPS II-R衛(wèi)星預報精度比GPS II-F衛(wèi)星差，多數(shù)配備Rb鐘的衛(wèi)星比Cs鐘的衛(wèi)星（G08，G24）精度要高；在GLONASS、Galileo中，衛(wèi)星軌道、

衛(wèi)星鐘類型對衛(wèi)星預報精度影響較小，除R10、R13、R19、E08及E18衛(wèi)星外，系統(tǒng)內其他衛(wèi)星預報精度基本相同。

預報方案不同，預報精度不同：①整體上，方案3預報精度最差，方案2、方案4的預報精度最好，方案1次之；②方案2相對于方案1，在BDS、Galileo中提升率較高，BDS的提升率最大可達36.79%，GLONASS的提升率較低；同時表現(xiàn)出隨著預報時長增加，預報精度提升率逐漸增加的特征；③方案4相對于方案3，四系統(tǒng)均有較高的提升率，平均提升率為35.19%；④方案2相對于方案4，BDS有較高的提升率，在1、3 h預報中，四系統(tǒng)中均有較好的提升率；對于GPS，GLONASS，Galileo，隨著預報時長的增加，方案4具有比方案2更優(yōu)的預報性能；⑤方案3在預報時長為1、3、6及12 h時，多數(shù)衛(wèi)星的預報精度劣于其他方案，這也說明在12 h內的衛(wèi)星鐘差預報中，Kalman濾波模型和考慮相關性的預報模型更具有預報優(yōu)勢。

由上可知，顧及衛(wèi)星間相關性的Kalman濾波、附加周期項二次多項式模型均比相應不考慮相關性的模型預報精度好，主要原因在于多衛(wèi)星系統(tǒng)鐘差整體估計是基于多衛(wèi)星系統(tǒng)聯(lián)合定軌同時解算軌道和鐘差得到的。對多顆衛(wèi)星進行整體建模與預報，不僅能夠很好地顧及衛(wèi)星鐘差序列之間的一致性，還能夠通過高穩(wěn)定度的衛(wèi)星提高低穩(wěn)定度衛(wèi)星鐘差序列的建模精度，從而提高多顆衛(wèi)星整體預報精度。同時從表1、表2可知，對于四系統(tǒng)鐘差產品，存在6、7、11、12及24 h等周期，在這些周期以內的預報（時長為1、3 h），顧及衛(wèi)星間相關性的Kalman濾波模型較顧及衛(wèi)星間相關性的附加周期項二次多項式模型、附加周期項二次多項式、Kalman濾波模型更具優(yōu)勢，分別提升了8.95 %，46.66 %，7.35 %，最大可提升12.77 %（GLONASS）、57.16 %（Galileo）、28.35 %（BDS），而在大于這些周期以外的預報，顧及衛(wèi)星間相關性的Kalman濾波模型的預報精度絕大多數(shù)情況是次優(yōu)的，可能是由于該模型沒有顧及衛(wèi)星鐘差的周期特征引起的。

4 結束語

本文推導了顧及衛(wèi)星間相關性的Kalman鐘差預報模型，分析了連續(xù)27天GBM最終產品的頻率穩(wěn)定性、周期特性，設計了4種實驗方案來比較分析所提方法在鐘差短期預報中的正確性和有效性，主要實驗結果如下：

1）衛(wèi)星鐘類型、衛(wèi)星軌道類型對衛(wèi)星鐘差頻率穩(wěn)定性有一定影響，配備Cs鐘的GLONASS衛(wèi)星頻率穩(wěn)定性較差，配備 Ph鐘的 Galileo衛(wèi)星頻率穩(wěn)定性較高；GPS II-F，Galileo IOV、Galileo FOC的衛(wèi)星鐘頻率穩(wěn)定性較高；同時對衛(wèi)星預報精度也有影響；

2）GBM的 BDS，GPS，GLONASS及 Galileo衛(wèi)星鐘差都有一定的周期性，分別有較明顯的12、24，6、12，11，7、14 h 周期。

3）在GBM鐘差數(shù)據的預報中：①顧及衛(wèi)星間相關性的Kalman濾波預報模型較傳統(tǒng)Kalman濾波模型精度平均提升了 15.11%，提升率最大值為36.56%（Galileo），BDS提升效果最好；②顧及衛(wèi)星間相關性的附加周期二次多項式預報模型較附加周期二次多項式預報模型精度平均提升了28.93%，最大提升率為51.95%（GPS），四系統(tǒng)提升效果相當；③在1、3 h的預報中，顧及衛(wèi)星間相關性的Kalman濾波較顧及衛(wèi)星間相關性附加周期的二次多項式預報模型精度平均提升了 8.95%，提升率最大值為12.77%（GLONASS）。