劉建成，楊睿峰，徐 赟，王茂磊，桑懷勝

（北京環(huán)球信息應(yīng)用開發(fā)中心，北京100094）

0 引 言

星地時(shí)間同步技術(shù)對(duì)衛(wèi)星導(dǎo)航系統(tǒng)的導(dǎo)航、授時(shí)以及定位精度有著直接的影響。星地?zé)o線電雙向時(shí)間比對(duì)的工作原理是地面和衛(wèi)星同時(shí)進(jìn)行星地偽碼測(cè)距，根據(jù)兩組測(cè)距信息解算出衛(wèi)星鐘相對(duì)于地面站基準(zhǔn)鐘的鐘差。這種方法抵消了兩組測(cè)距共同的誤差，所以，鐘差測(cè)量精度很高。大多數(shù)實(shí)時(shí)導(dǎo)航用戶采用的是廣播星歷給出的鐘差，其精度將影響導(dǎo)航定位的精度，因此，研究星地?zé)o線電雙向時(shí)間比對(duì)體制下的衛(wèi)星鐘差預(yù)報(bào)精度有著重要意義。

常用的鐘差預(yù)報(bào)方法包括最小二乘方法［1－2］、FIR 濾 波［1，3］、Kalman 濾 波［1，4－7］。 當(dāng) 鐘 差 狀 態(tài) 模型是線性的，且狀態(tài)噪聲統(tǒng)計(jì)特性和測(cè)量噪聲統(tǒng)計(jì)特性已知時(shí)，Kalman濾波性能最優(yōu)，Kalman濾波在鐘差預(yù)報(bào)方面獲得了越來越多的應(yīng)用。對(duì)于Kalman濾波器，估計(jì)誤差方差隨著時(shí)間的遞推逐漸減小，直至達(dá)到穩(wěn)態(tài)，此時(shí)估計(jì)誤差方差最小，在穩(wěn)態(tài)條件下得到的預(yù)報(bào)誤差也最小，因此，穩(wěn)態(tài)情況下的鐘差預(yù)報(bào)精度可作為分析星地時(shí)間同步性能的重要依據(jù)。

Meditch研究了衛(wèi)星鐘差狀態(tài)噪聲為頻率白噪聲情況下的Kalman濾波器鐘差預(yù)報(bào)精度，利用Kalman濾波器狀態(tài)轉(zhuǎn)移矩陣得到了最小預(yù)報(bào)誤差的解析表達(dá)式［8］。Beisner研究了衛(wèi)星鐘狀態(tài)噪聲為頻率白噪聲、頻率閃爍噪聲和頻率隨機(jī)游走噪聲且不考慮測(cè)量噪聲情況下的Wiener預(yù)測(cè)器的最小鐘差預(yù)報(bào)誤差，利用Wiener方法直接從頻域得到了比Meditch得到的更具有普適性的解析形式。目前，Kalman濾波器采用的鐘差狀態(tài)噪聲逐漸擴(kuò)展到3階［9－11］，不但包括頻率白噪聲、頻率隨機(jī)游走噪聲，還包括頻率漂移隨機(jī)游走噪聲。研究鐘差狀態(tài)噪聲包括3階白噪聲且測(cè)量噪聲為白噪聲情況下的基于Kalman濾波器的衛(wèi)星鐘差穩(wěn)態(tài)精度。

根據(jù)Kalman濾波器達(dá)到穩(wěn)態(tài)時(shí)的條件建立方程組，然后解該方程組，可以得到穩(wěn)態(tài)解。該方法的優(yōu)點(diǎn)是概念直觀，但當(dāng)目標(biāo)狀態(tài)變量增多時(shí)解方程組存在很大困難。通過消除Kalman遞歸形式中的增益項(xiàng)可得到預(yù)測(cè)協(xié)方差矩陣的遞歸公式，即Riccati方程，Vaughan解決了離散Riccati方程的非遞歸代數(shù)解問題［12］，利用這個(gè)結(jié)果可直接得到穩(wěn)態(tài)解，繞開了解方程組的問題。在建立衛(wèi)星鐘差的狀態(tài)方程和測(cè)量方程的基礎(chǔ)上，利用這種方法得到Kalman濾波器的穩(wěn)態(tài)解，進(jìn)一步得到預(yù)報(bào)誤差，建立起穩(wěn)態(tài)情況下的鐘差預(yù)報(bào)誤差與衛(wèi)星鐘狀態(tài)誤差、測(cè)量誤差及采樣時(shí)間間隔之間的解析關(guān)系，最后做了典型參數(shù)情況下的數(shù)值分析。

1 衛(wèi)星鐘差Kalman濾波方程

采用3狀態(tài)衛(wèi)星鐘差模型，包括相位、鐘速和頻率漂移值?？紤]3類鐘差噪聲［10］，包括頻率白噪聲、頻率隨機(jī)游走噪聲和頻率漂移隨機(jī)游走噪聲。

若采樣時(shí)間間隔為T，衛(wèi)星鐘差狀態(tài)向量為xk＝ ［a0a1a2］T，其中a0為鐘差，a1為鐘速，a2為頻率漂移量，那么建立多項(xiàng)式狀態(tài)方程為［10］

式中狀態(tài)轉(zhuǎn)移矩陣誤差矩陣為

式中q1、q2、q3分別為頻率白噪聲、頻率隨機(jī)游走噪聲和頻率漂移隨機(jī)游走噪聲的譜密度。

根據(jù)星地?zé)o線電雙向時(shí)間傳遞的原理，鐘差測(cè)量方程用xk表示為

從狀態(tài)方程和測(cè)量方程可以看出，兩者均為線性方程，則最優(yōu)狀態(tài)估計(jì)就是卡爾曼濾波，其誤差協(xié)方差陣遞推過程為

穩(wěn)態(tài)協(xié)方差矩陣?P定義為

式中:σ211為鐘差預(yù)測(cè)誤差方差；σ212為鐘差、鐘速預(yù)測(cè)誤差相關(guān)系數(shù)；σ213為鐘差、頻率漂移量預(yù)測(cè)誤差相關(guān)系數(shù)；σ222為鐘速預(yù)測(cè)誤差方差；σ223為鐘速、頻率漂移量預(yù)測(cè)誤差相關(guān)系數(shù)；σ233為頻率漂移量預(yù)測(cè)誤差方差。

2 衛(wèi)星鐘差的穩(wěn)態(tài)精度和預(yù)報(bào)精度

推導(dǎo)基于Kalman濾波器的衛(wèi)星鐘差穩(wěn)態(tài)精度，得到衛(wèi)星鐘差預(yù)報(bào)誤差協(xié)方差陣與穩(wěn)態(tài)協(xié)方差矩陣和衛(wèi)星鐘狀態(tài)誤差矩陣的關(guān)系。

2．1 Kalman濾波器的穩(wěn)態(tài)解

根據(jù)文獻(xiàn)［12］提出的離散Riccati方程的非遞歸代數(shù)解算法，構(gòu)造矩陣

式中F、H、R、Q由式（1）和式（3）確定。

由于星鐘的系統(tǒng)模型為三階，所以Kf為6階特征值。若矩陣的特征值為λi（i＝1，2，3），則其對(duì)應(yīng)的特征向量為

穩(wěn)態(tài)解是S1＝λ1＋λ2＋λ3、S2＝λ1λ2＋λ2λ3＋λ3λ1和S3＝λ1λ2λ3的函數(shù)。下面說明如何確定S1、S2、S3．

矩陣Kf的特征方程為

由矩陣Kf的性質(zhì)得到

式中λ1、λ2、λ3均為模大于1的特征值。

把式（11）展開，并與式（10）比較，得到

令S3＝x2，采用文獻(xiàn)［13］的做法，解上式方程組得到

式中的x是式（13）代入式（12）而得到的一元四次方程的解的最大值

從式（14）可以看出，基于Kalman濾波器的衛(wèi)星鐘差穩(wěn)態(tài)精度取決于參數(shù)q1、q2、q3、T、σ2．

2．2 衛(wèi)星鐘差預(yù)報(bào)精度

假設(shè)Kalman濾波器開始預(yù)測(cè)工作前已經(jīng)達(dá)

等式右端第一部分反映了Kalman濾波穩(wěn)態(tài)誤差對(duì)預(yù)報(bào)精度的影響，取決于參數(shù)q1、q2、q3、T、N、σ2．第二部分反映了衛(wèi)星鐘的狀態(tài)誤差對(duì)預(yù)報(bào)精度的影響，取決于參數(shù)q1、q2、q3、T、N．

如果定義間隔為T時(shí)的狀態(tài)噪聲矩陣函數(shù)為Q［T］，由式（2）確定，那么NT 預(yù)報(bào)誤差協(xié)方差陣還可表示為

式中NT衛(wèi)星鐘差的預(yù)報(bào)誤差方差為σ211（N）＝P?N（1，1）．

當(dāng)不考慮頻率隨機(jī)游走噪聲和頻率漂移隨機(jī)游走噪聲時(shí)，即q2＝q3＝0，此時(shí)β＝γ＝0，x＝到穩(wěn)態(tài)，則鐘差參數(shù)預(yù)報(bào)初始誤差協(xié)方差陣為?P1＝?P．根據(jù)Kalman濾波器預(yù)測(cè)方程，預(yù)報(bào)誤差協(xié)方差陣的遞推關(guān)系為

根據(jù)此遞推關(guān)系，得到衛(wèi)星鐘差狀態(tài)向量的N步預(yù)報(bào)誤差協(xié)方差陣，即衛(wèi)星鐘差狀態(tài)向量的NT預(yù)報(bào)誤差協(xié)方差陣為度，再根據(jù)式（17）或（18）得到鐘差預(yù)報(bào)誤差方差，結(jié)果與文獻(xiàn)［8］一致。

3 數(shù)值分析

根據(jù)上面得出的理論結(jié)果，采樣時(shí)間間隔、偽碼距離測(cè)量誤差均方根對(duì)基于Kalman濾波器的鐘差穩(wěn)態(tài)誤差均方根和鐘差預(yù)報(bào)誤差均方根的影響做了數(shù)值分析。數(shù)值分析采用的衛(wèi)星鐘狀態(tài)誤差參數(shù)為［11］:q1＝1．11×10－22s2／s、q2＝2．22×10－32s2／s3、q3＝6．66×10－45s2／s5．

3．1 鐘差穩(wěn)態(tài)精度與偽距測(cè)量精度的關(guān)系

分析了采樣時(shí)間間隔分別為5s、10s、15s時(shí)，基于Kalman濾波器的鐘差穩(wěn)態(tài)誤差均方根隨偽距測(cè)量誤差均方根的變化關(guān)系，如圖1所示。從圖1可以看出，偽距測(cè)量誤差均方根越大，鐘差穩(wěn)態(tài)誤差均方根就越大。采樣時(shí)間間隔越大，鐘差穩(wěn)態(tài)誤差均方根就越大。因此，提高鐘差穩(wěn)態(tài)精度的方法包括提高偽距測(cè)量精度和提高數(shù)據(jù)率。

3．2 鐘差預(yù)報(bào)精度與偽距測(cè)量精度的關(guān)系

圖1 鐘差穩(wěn)態(tài)精度與偽距測(cè)量精度的關(guān)系曲線

分析了采樣時(shí)間間隔分別為5s、10s、15s時(shí)，基于Kalman濾波器的鐘差預(yù)報(bào)誤差均方根隨偽距測(cè)量誤差均方根的變化關(guān)系。1h和2h的鐘差預(yù)報(bào)誤差均方根與偽碼測(cè)量誤差均方根的關(guān)系曲線分別如圖2（a）和圖2（b）所示。

從圖2可以看出，偽距測(cè)量誤差均方根越大，鐘差預(yù)報(bào)誤差均方根就越大。采樣時(shí)間間隔越大，鐘差預(yù)報(bào)誤差均方根就越大。

圖2 鐘差預(yù)報(bào)精度與偽距測(cè)量精度的關(guān)系曲線

4 結(jié) 論

研究了衛(wèi)星鐘差包括3階白噪聲情況下的基于Kalman濾波器的衛(wèi)星鐘差穩(wěn)態(tài)精度。在建立衛(wèi)星鐘差的狀態(tài)方程和測(cè)量方程的基礎(chǔ)上，利用離散Riccati方程的非遞歸代數(shù)解，得到了Kalman濾波器的穩(wěn)態(tài)解，進(jìn)一步得到了衛(wèi)星鐘差預(yù)報(bào)誤差，建立起鐘差預(yù)報(bào)誤差與衛(wèi)星鐘狀態(tài)誤差、測(cè)量誤差及采樣時(shí)間間隔之間的解析關(guān)系。做了典型參數(shù)情況下的數(shù)值分析。下一步的工作是分析衛(wèi)星鐘的閃爍噪聲對(duì)基于Kalman濾波器的衛(wèi)星鐘差預(yù)報(bào)精度的影響程度。

［1］ 劉 利．相對(duì)論時(shí)間比對(duì)理論與高精度時(shí)間同步技術(shù)［D］．鄭州:解放軍信息工程大學(xué)，2004．

［2］ DAVIS J A，HARRIS P M，et al．Least－squares analysis of two－way Satellite time and frequency Transfer Measurements［C］∥Proceedings of the 33thAnnual Precise Time and Time Interval（PTTI）Applications and Planning Meeting，Long Beach，California，USA，2001（11）:121－132．

［3］ SHMAL I Y，IBARRA－MANZHNO O．An optimal FIR filtering algorithm for a time error model of a precise clock［C］∥Proceedings of the 34thAnnual Precise Time and Time Interval（PTTI）Applications and Planning Meeting，Reston，Virginia，USA，2002（12）:53－68．

［4］ DAVIS J A，STACEY P W，et al．Combining time transfer measurements using a Kalman filter［C］∥Proceedings of the 34thAnnual Precise Time and Time Interval（PTTI）Applications and Planning Meeting，Reston，Virginia，USA，2002（12）:53－68．

［5］ DAVIS J A，GREENGHALL C A，STACEY P W．A Kalman filter clock algorithm for use in the presence of flicker frequency modulation noise［C］∥Proceedings of the 35thAnnual Precise Time and Time Interval（PTTI）Applications and Planning Meeting，Reston，Virginia，USA，2003（12）:53－68．

［6］ 陳小敏，李孝輝．基于自適應(yīng)Kalman濾波器的原子鐘信號(hào)預(yù)測(cè)方法［J］．吉林大學(xué)學(xué)報(bào)·理學(xué)版，2009，47（3）:591－593．

［7］ GALLEANI L，TAVELLA P．Time and the Kalman filter－applications of optimal estimation to atomic timing［J］．IEEE control systems magazine，2010，30（4）:44－65．

［8］ BEUSMER H M．Clock error with a Wiener predic－tor and by numerical Calculation［J］．IEEE Transactions on Instrumentation and measurement，1980，29（2）:105－113．

［9］ HUTSELL S T，REID W G，GRUM J D，et al．Operational use of the Hadamard variance in GPS［C］∥Proceedings of the 28thAnnual Precise Time and Time Interval（PTTI）Applications and Planning Meeting，Reston，Virginia，USA，1996（12）:201－214．

［10］ HUTSELL S T．Relating the Hadamard variance to MCS Kalman filter clock estimation［C］∥Proceedings of the 27thAnnual Precise Time and Time Interval（PTTI）Applications and Planning Meeting，San Diego，California，USA，1995（12）:291－302．

［11］ HUTSELL S T．Fine tuning GPS clock estimation in the MCS［C］∥Proceedings of the 26thAnnual Precise Time and Time Interval（PTTI）Applications and Planning Meeting，Reston，Virginia，USA，1994（12）:63－74．

［12］ VAUGHAN D R．A non－recursive algebraic solution for the discrete riccati equation［J］．IEEE Transactions on Automatic Control，1970，15（1）:597－599．

［13］ RAMACHANDRA K V，MOHAN B R，GEETHA B R．A three－state Kalman tracker using position and rate measurements［J］．IEEE Transactions on Aerospace and Electronic Systems，1993，29（1）:215－222．