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

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

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

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

4.中國(guó)科學(xué)院大學(xué)，北京100049

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

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

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

1 星間鏈路觀測(cè)原理

1.1 星間鏈路觀測(cè)模型

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

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

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

圖1 歷元?dú)w算示意圖Figure1 Epoch-return diagram

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

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

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

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

式中

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

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

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

1.2 星鐘跳變觀測(cè)量

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

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

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

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

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

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

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

2 Kalman 濾波原理

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

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

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

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

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

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

3 仿真校驗(yàn)

3.1 仿真條件

根據(jù)Kalman 濾波器特性可知，R和Q的相對(duì)大小決定了Kalman 濾波器對(duì)于預(yù)測(cè)值和觀測(cè)值的信任程度.本仿真對(duì)兩者進(jìn)行了折中，選擇觀測(cè)誤差協(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ù)工程需求，虛警率及漏檢率不得超過(guò)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，即在本文仿真中可以檢測(cè)的異常跳變大小不低于3.11×10?12×300=9.33×10?10s.且當(dāng)測(cè)距周期T相同時(shí)，相位跳變和頻率跳變的門限值相同.

若星間鏈路中僅有一顆衛(wèi)星的鐘差數(shù)據(jù)存在異常，則經(jīng)過(guò)兩條星間鏈路的雙向測(cè)距，即最少需要3 顆衛(wèi)星可以判斷出異常衛(wèi)星.選取Sate1、Sate2、Sate3 衛(wèi)星，考慮到星間鏈路的工程實(shí)際情況，選取一天的連續(xù)數(shù)據(jù)，檢測(cè)主體為Sate1，假設(shè)Sate2 的星鐘數(shù)據(jù)出現(xiàn)異常，由于北斗導(dǎo)航衛(wèi)星真實(shí)在軌運(yùn)行過(guò)程中，星鐘異常引起跳變會(huì)導(dǎo)致星鐘異常時(shí)刻的數(shù)據(jù)及之后的數(shù)據(jù)均帶有該異常值，因此仿真時(shí)，在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ù)分兩組添加，通過(guò)檢測(cè)Sate1、Sate2 和Sate1、Sate3 之間的鐘差數(shù)據(jù)，進(jìn)行兩次仿真校驗(yàn).

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

3.2.1 相位跳變

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

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

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

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

圖4 Sate1 和State3 相對(duì)鐘差中的頻率偏差和預(yù)測(cè)誤差Figure4 Frequency deviation and prediction error of relative clock difference of Sate 1 and Sate 3

3.2.2 頻率跳變

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

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

圖5 Sate1 和State2 相對(duì)鐘差中頻率跳變的頻率偏差和預(yù)測(cè)誤差Figure5 Frequency deviation and prediction error of relative clock difference of Sate1 and Sate2

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

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

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

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

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

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

4 結(jié) 語(yǔ)

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