陳婷婷，林寶軍，龔文斌，馮 磊，常家超，林 夏

1.中國科學院上海微系統(tǒng)與信息技術(shù)研究所，上海200050

2.上海科技大學信息科學與技術(shù)學院，上海201210

3.中國科學院微小衛(wèi)星創(chuàng)新研究院，上海200120

4.中國科學院大學，北京100049

目前，四大全球衛(wèi)星導航系統(tǒng)（global navigation satellite system,GNSS）——美國的GPS,俄羅斯的GLONASS，歐盟的GALILEO 以及中國的北斗衛(wèi)星導航系統(tǒng)（BeiDou navigation satellite system,BDS）已經(jīng)得到了廣泛的應(yīng)用.四大系統(tǒng)均采用基于地面站的完好性監(jiān)測處理機制[1]進行工作，該機制下，若想保障全球范圍內(nèi)的完好性監(jiān)測性能，需要實現(xiàn)地面系統(tǒng)高水平密集布站.而對于北斗系統(tǒng)來說，全球布站存在極大困難，短時間內(nèi)無法得到良好效果.

通過衛(wèi)星自主運行的方式在空間段實現(xiàn)系統(tǒng)完好性監(jiān)測可以極大程度解決全球布站困難的問題.自主運行的自主守時和自主時間同步是指在不依賴地面站的情況下，通過星間鏈路測量使用導航星座衛(wèi)星原子鐘組共同維持北斗時間基準和維持衛(wèi)星之間時間同步的過程.隨著GNSS 系統(tǒng)的不斷進步，科研人員已開始探討如何利用星間鏈路所獲取的雙向測距數(shù)據(jù)來獲取衛(wèi)星星歷參數(shù)并以此監(jiān)測星歷參數(shù)的完好性.

目前利用星間鏈路實現(xiàn)完好性監(jiān)測的方法有：基于星間鏈路的單星完好性監(jiān)測法、利用平穩(wěn)時間序列的χ2檢驗法[2]、基于星間鏈路檢驗統(tǒng)計量和最小監(jiān)測精度分析的完好性監(jiān)測法.但是，這些方法存在以下不足：1）問題的研究尚停留在理論算法上，并沒有真實星上數(shù)據(jù)進行檢驗；2）沒有開展對自主導航時間系統(tǒng)的影響分析；3）沒有充分結(jié)合現(xiàn)有的星間鏈路體制，需要在當前星上系統(tǒng)中增加星載設(shè)備，或者對當前星上系統(tǒng)提出需求和改動卻無法真正實現(xiàn).針對以上不足，本文提出了一種基于星間鏈路的星上時間自主完好性監(jiān)測技術(shù)，利用真實星間鏈路的測距值，解算出鐘差信息，從而利用二狀態(tài)的Kalman 濾波算法，實現(xiàn)異常星鐘的快速檢測，同時可分辨出該星鐘異常是相位跳變還是頻率跳變.

1 星間鏈路觀測原理

1.1 星間鏈路觀測模型

不同于衛(wèi)星與地面的通信方式，衛(wèi)星通過輪詢的方式完成彼此之間的雙向測距，且雙向測距采用時分體制，載波頻率處于Ka波段.假設(shè)衛(wèi)星i在ts時刻發(fā)送，衛(wèi)星j在tr時刻接受；衛(wèi)星j在時刻發(fā)送，衛(wèi)星i在時刻接收.則兩顆衛(wèi)星之間的觀測方程如下[3]：

式中，ρij和ρji為測距偽量，rij和rji為衛(wèi)星之間的空間距離，c為光速，δi(ts)為衛(wèi)星i在信號發(fā)送時刻的鐘差，δj(tr)衛(wèi)星j在信號接收時刻的鐘差，δj()為衛(wèi)星在信號發(fā)送時刻的鐘差，δi()為衛(wèi)星i在信號接收時刻的鐘差，D為衛(wèi)星之間發(fā)射通道和接受通道總時延，I為電離層延遲誤差，T為對流層延遲誤差，ε為各種測距隨機誤差和[4-5].由文獻[5]可知，在星間鏈路的情況下，電離層誤差I(lǐng)由于頻率較高誤差較小，對流層誤差T由于軌道較高誤差較小，這兩種均可不考慮.

歷元歸算是指將當前時刻的星間鏈路測距值通過歷元分析規(guī)劃到某一固定時刻的計算過程.在北斗Ka 星間鏈路中，為實現(xiàn)實時定軌算法，需要將不同時刻的Ka 星間鏈路觀測量統(tǒng)一至同一時間起點[6].歷元歸算原理如圖1所示.

圖1 歷元歸算示意圖Figure1 Epoch-return diagram

根據(jù)文獻[7]，可將式(1)和(2)改寫為

式中，和為改化偽距，ri(tr)和rj(ts)分別為衛(wèi)星i和j在tr和ts時刻的真實位置向量，εij和εji為偽距觀測誤差.記歷元歸算的目標時刻為te，?t為衛(wèi)星信號接收時刻與歷元歸算目標時刻之間的時間間隔.以式(3)為例，在實際歷元歸算過程中，?t的絕對值一般小于60 s，即

由于歸算時間間隔?t為小量，因此可僅考慮衛(wèi)星鐘差的一階項，假設(shè)衛(wèi)星i和衛(wèi)星j的頻漂分別為d1和d2，并用τ來表示信號傳輸時延，一般情況下τ0.25 s，因此有

將式(6)和(7)代入式(3)有

式中

在式(9)和(10)中，ρ(te)為歷元歸算目標時刻對應(yīng)的偽距，?ρ即為歸算偽距修正量.d1?t相對于?t為小量，因此有?t=tr ?δi(te)?te，在考慮了衛(wèi)星本身鐘差之后的歸算偽距修正量為

式中，tr為衛(wèi)星i接收信號時刻鐘面時；δi(te)為歷元歸算目標時刻衛(wèi)星i的鐘差，可以用衛(wèi)星i實時估計的鐘差和鐘漂進行預測.ri(te),ri(tr ?δi(te))為衛(wèi)星i在te時刻的位置，rj(te),rj(tr ?δi(te)?τ)為衛(wèi)星j在te時刻的位置，該值可通過北斗衛(wèi)星廣播星歷計算.τ為信號傳播時延[8-9]，可以根據(jù)文獻[10]求得.經(jīng)過上述歷元歸算后，可將式(3)和(4)統(tǒng)一歸算至te時刻，結(jié)果為

將式(3)和式(4)相減，可得到用于時間同步計算的觀測方程，由于歸算會引入部分誤差，因此觀測誤差替換為v來表示[10].

1.2 星鐘跳變觀測量

假設(shè)衛(wèi)星i在某段時間內(nèi)與L顆衛(wèi)星互相可見，且與其中的m顆可見衛(wèi)星建鏈.在一段時間內(nèi)，衛(wèi)星i將與衛(wèi)星j建鏈且進行雙向測距.假設(shè)第n次觀測值與第n?1次觀測值的時間間隔為T，即測距周期.

假設(shè)衛(wèi)星i與j在第n ?1 次和第n次雙向測距的觀測方程值分別為δtrs,n和δtrs,n?1，則在一個測距周期T內(nèi)的平均頻率偏差如式(15)所示，本文利用該式對星鐘的相位和頻率跳變進行監(jiān)測.

從式(14)中可以看出，星鐘相對鐘差觀測量中的測量噪聲是影響相位和頻率跳變監(jiān)測靈敏度的最重要因素.本文用北斗新一代導航衛(wèi)星獲得的星間鏈路實測數(shù)據(jù)對星鐘相對鐘差的測量噪聲進行分析.為保證數(shù)據(jù)的完整性，實驗使用2018年10月25日到2018年10月31日4 顆新一代北斗導航衛(wèi)星（Sate1、Sate2、Sate3、Sate4）之間的雙向測量數(shù)據(jù)對測量噪聲進行分析.

首先，按照式(3)～(14)對星間雙向測量數(shù)據(jù)進行處理，獲得星間相對鐘差.然后，對相對鐘差進行一階線性擬合，以其擬合殘差評估星間測量精度，4 顆衛(wèi)星之間的相對鐘差擬合殘差如圖2所示，統(tǒng)計結(jié)果如表1所示.由于衛(wèi)星軌道位置關(guān)系約束和工程原因，某些衛(wèi)星之間無星間測量數(shù)據(jù)，同時某些星間測量數(shù)據(jù)不連續(xù)，使得圖2中擬合殘差不連續(xù)[11].

圖2 星鐘相對鐘差擬合殘差Figure2 Fitting residual of relative clock bias

表1 星間鏈路擬合殘差統(tǒng)計Table1 Inter-satellite link fitting residual statistics

由圖2和表1可以看出，北斗系統(tǒng)星間鏈路測量噪聲小于0.15 ns，由于測量殘差均值為0，因此其測量噪聲標準差小于0.15ns，則式(14)中測量噪聲v服從N(0,0.152×10?18s2).由于前后2 個時刻的鐘差測量值誤差不相關(guān)，則式(15)中vn服從N(0,4.5×10?20s2/T2).

2 Kalman 濾波原理

Kalman 濾波為離散時間系統(tǒng)內(nèi)常用的濾波方法，常被用于離散時間系統(tǒng)參數(shù)最優(yōu)估計[10]，其基本思想是：首先基于系統(tǒng)狀態(tài)方程對系統(tǒng)狀態(tài)量進行時間更新，然后根據(jù)系統(tǒng)測量方程使用觀測量對狀態(tài)量進行測量更新，獲得狀態(tài)量的最優(yōu)估計[11].其中，作為Kalman 濾波算法的先驗統(tǒng)計信息，系統(tǒng)狀態(tài)方程的過程噪聲協(xié)方差矩陣和測量方程的觀測噪聲協(xié)方差矩陣決定了對時間更新狀態(tài)量和觀測量的信任程度[12].Kalman 濾波器包含歷史數(shù)據(jù)及未來觀測量在每個歷元的觀測值及估計值，通過對歷史數(shù)據(jù)的觀測量，Kalman 濾波器通過預測方程及狀態(tài)方程對未來時刻進行估計，該過程稱為預測過程，并通過上一時刻得到校驗的先驗估計值對其進行校驗，該步驟為更新過程.

在狀態(tài)預測步驟，假設(shè)xk表示系統(tǒng)狀態(tài)，= [x0x1] 代表第k個歷元時Kalman 濾波器對xk的最優(yōu)估計，此處，x0和x1分別代表鐘差和鐘速，uk?1為輸入量，A(tk,tk?1)代表tk?1到tk歷元的狀態(tài)轉(zhuǎn)移矩陣，x?k為沒有得到測量值yk的驗證值，即先驗估計值.在下文中，用上標“”表示估計值，右上標“–”表示先驗估計值，右下標“k”代表第k個歷元.則狀態(tài)方程為

由文獻[12]可知，P為狀態(tài)先驗估計值的均方誤差陣，在Kalman 濾波過程中，由于過程噪聲w具有未知性，所以在P的計算過程中加入過程噪聲的協(xié)方差陣Q，用來降低狀態(tài)估計值的可靠性且保持均方誤差陣的正確性，則在k歷元時

在狀態(tài)更新步驟中，假設(shè)K代表Kalman 濾波增益，R為觀測噪聲協(xié)方差矩陣，C代表觀測量與系統(tǒng)狀態(tài)之間的關(guān)系矩陣，v代表噪聲向量，則有

假設(shè)y代表觀測矩陣，則觀測方程為[13]

本文選取二狀態(tài)的Kalman 濾波器，且設(shè)置Kalman 濾波器參數(shù)為

3 仿真校驗

3.1 仿真條件

根據(jù)Kalman 濾波器特性可知，R和Q的相對大小決定了Kalman 濾波器對于預測值和觀測值的信任程度.本仿真對兩者進行了折中，選擇觀測誤差協(xié)方差R=(4.5×10?11/T)2，系統(tǒng)誤差協(xié)方差Q=(4.5×10?11/T)2/100.于是，vn ～N(0,4.5×10?20s2/T2)，根據(jù)工程需求，虛警率及漏檢率不得超過10?5，根據(jù)正態(tài)分布查表法可知，門限值不得高于4.40σ，由于本文中噪聲分布服從N(0,4.5×10?20T2s2s?2)，而在仿真中T=5 min=300 s，則因此門限值不得高于3.11×10?12s·s?1，即在本文仿真中可以檢測的異常跳變大小不低于3.11×10?12×300=9.33×10?10s.且當測距周期T相同時，相位跳變和頻率跳變的門限值相同.

若星間鏈路中僅有一顆衛(wèi)星的鐘差數(shù)據(jù)存在異常，則經(jīng)過兩條星間鏈路的雙向測距，即最少需要3 顆衛(wèi)星可以判斷出異常衛(wèi)星.選取Sate1、Sate2、Sate3 衛(wèi)星，考慮到星間鏈路的工程實際情況，選取一天的連續(xù)數(shù)據(jù)，檢測主體為Sate1，假設(shè)Sate2 的星鐘數(shù)據(jù)出現(xiàn)異常，由于北斗導航衛(wèi)星真實在軌運行過程中，星鐘異常引起跳變會導致星鐘異常時刻的數(shù)據(jù)及之后的數(shù)據(jù)均帶有該異常值，因此仿真時，在Sate2 的鐘差數(shù)據(jù)中添加P= 4 ns 作為相位跳變異常，添加到第235 min 及之后所有數(shù)據(jù)中.在Sate2 的鐘差數(shù)據(jù)中添加?f= 8×10?9s·s?1作為頻率異常跳變，添加到第50 min 及之后所有數(shù)據(jù)中.相位跳變和頻率跳變數(shù)據(jù)分兩組添加，通過檢測Sate1、Sate2 和Sate1、Sate3 之間的鐘差數(shù)據(jù)，進行兩次仿真校驗.

3.2 仿真結(jié)果及分析

3.2.1 相位跳變

從測距值中利用式(1)～(14)解算出鐘差信息，按照上文添加異常數(shù)據(jù)，計算Sate1、Sate2間的相對鐘差，與Sate1、Sate3 間的相對鐘差，通過Kalman 濾波器，結(jié)果如圖3和4 所示.圖3為Sate1、Sate2 之間的相對鐘差在發(fā)生相位跳變時的頻率偏差和預測誤差，圖4為Sate1、Sate3 之間的相對鐘差的頻率偏差和預測誤差，通過對比兩圖可以得出在第235 min時，Sate1、Sate2 之間的相對鐘差存在異常，而Sate1、Sate3 并未出現(xiàn)該情況，由于僅存在一個星鐘異常，則Sate1 可以判定是Sate2 在第235 min 出現(xiàn)了星鐘異常.

從圖3可以看出，頻率偏差在發(fā)生跳變后很快恢復到之前的值，根據(jù)式(15)及鐘差數(shù)據(jù)的性質(zhì)可以判定，該異常為相位異常.

如圖3和4 所示，此時預測誤差為真值與預測值之差，可以根據(jù)預測誤差設(shè)置Kalman 濾波的閾值.當觀測誤差越大時，所需閾值越大，相當于噪聲越大，越不容易監(jiān)測出跳變.

圖3 Sate1 和Sate2 相對鐘差中相位跳變的頻率偏差和預測誤差Figure3 Frequency deviation and prediction error of relative clock difference of Sate 1 and State 2

圖4 Sate1 和State3 相對鐘差中的頻率偏差和預測誤差Figure4 Frequency deviation and prediction error of relative clock difference of Sate 1 and Sate 3

3.2.2 頻率跳變

圖5為Sate1 和Sate2 之間的相對鐘差在發(fā)生頻率跳變時的頻率偏差，對比圖4的Sate1和Sate3 之間的相對鐘差的頻率偏差，可以得出在第50 min 時，Sate1 和Sate2 之間的相對鐘差存在異常，而Sate1 和Sate3 并未出現(xiàn)該情況，在一個星鐘異常，則Sate1 可以判定是Sate2在第50 min 出現(xiàn)了星鐘異常.

從圖5可以看出，頻率偏差在發(fā)生跳變后逐漸在新的頻率處達到穩(wěn)定，根據(jù)式(15)及鐘差數(shù)據(jù)的性質(zhì)可以判定，該異常為頻率異常.

圖5 Sate1 和State2 相對鐘差中頻率跳變的頻率偏差和預測誤差Figure5 Frequency deviation and prediction error of relative clock difference of Sate1 and Sate2

如圖4和5 所示，此時預測誤差為真值與預測值之差，可以根據(jù)預測誤差設(shè)置Kalman 濾波的閾值.當觀測誤差越大時所需閾值越大，相當于噪聲越大越不容易監(jiān)測出跳變.

為了考察所設(shè)置門限值，本文針對不同量級的鐘跳進行了仿真，以相位跳變?yōu)槔?，對比之前的仿真，?35 min 時在Sate2 的鐘差數(shù)據(jù)中添加P=0.5 ns 和P=1.6 ns 作為相位跳變異常.結(jié)果如圖6和7 所示.

圖6 Sate1 和Sate2 相位跳變的預測誤差(P =0.5 ns)Figure6 Prediction error of relative clock difference of Sate1 and Sate2 when P =0.5 ns

圖7 Sate1 和Sate2 相位跳變的預測誤差(P =1.6 ns)Figure7 Prediction error of relative clock difference of Sate1 and Sate2 when P =1.6 ns

對比圖3、6、7 可知，針對本文中所設(shè)置的門限值，當跳變大小低于1 ns 時，Kalman 濾波器無法檢測出異常，當跳變大小高于2 ns 時，檢測概率高于99.999%.

對比圖3、4、5 可知，對于相位跳變，相對鐘差的平均頻率偏差在發(fā)生躍變之后會很快回到之前的值，而對于頻率跳變，相對鐘差的平均頻率偏差在發(fā)生躍變之后將在新的頻率處達到平穩(wěn)，因此可以根據(jù)相對鐘差的平均頻率偏差變化情況來判斷發(fā)生的異常為相位跳變還是頻率跳變.

4 結(jié) 語

針對北斗衛(wèi)星系統(tǒng)難以在全球布站的現(xiàn)實情況，提出了一種星上時間的自主完好性監(jiān)測方法，基于北斗衛(wèi)星真實在軌數(shù)據(jù)與北斗衛(wèi)星星間鏈路的測距與數(shù)傳功能，對雙向測距值進行歷元歸算得到兩顆衛(wèi)星的鐘差，作為Kalman 濾波器的觀測量實現(xiàn)了對相位跳變和頻率跳變的監(jiān)測和識別.同時根據(jù)工程實際情況及測距周期，提出工程合適的門限值.仿真結(jié)果表明，當在建鏈衛(wèi)星中僅存在一個異常星鐘時，根據(jù)3顆衛(wèi)星的鐘差數(shù)據(jù)，利用Kalman 濾波器可以監(jiān)測并判斷出異常星鐘，且可以根據(jù)Kalman 濾波器的濾波結(jié)果分辨相位跳變和頻率跳變.當存在超過一個異常星鐘時，3顆衛(wèi)星無法判斷出異常星鐘，此時需要根據(jù)具體情況分析所需衛(wèi)星數(shù).