一種BDS星鐘異常實時探測與處理算法

趙　丹，周佩元，杜　蘭

信息工程大學導航與空天目標工程學院，河南 鄭州，450001

?

一種BDS星鐘異常實時探測與處理算法

趙丹，周佩元，杜蘭

信息工程大學導航與空天目標工程學院，河南 鄭州，450001

通過連續監測鐘差時序數據，可以快速探測星載原子鐘異常，從而確保衛星導航系統的完好性。本文結合廣義似然比檢驗(GLRT)方法和多項式鐘差模型，提出了一種新的實時原子鐘異常探測與處理方法：首先，基于GLRT方法進行實時異常探測；然后采用多項式鐘差模型進行短期鐘差外推，在此基礎上剔除異常值并以估值代替；最后，利用武漢大學分析中心提供的精密BDS鐘差數據進行仿真。驗證結果表明，該方法能夠快速有效地探測和處理原子鐘異常(尤其是頻率跳變異常)，較好地保證了BDS原子鐘數據的完好性。

BDS星鐘；原子鐘異常探測；GLRT探測器；鐘差模型

1　引　言

原子鐘是許多復雜系統的關鍵部件之一，如電信網絡系統、電力系統以及全球衛星導航系統等[1]。原子鐘信號優異的計量學特性是實現這些應用的基礎。但是在原子頻標長期運行測量過程中，常會由于斷電、設備故障和外部干擾等因素引起數據異常，如信號缺失、頻率漂移、相位跳變和頻率跳變等[10]。因此，必須對原子鐘進行實時監測和評估，對其異常進行實時探測和處理，并在短時間內給終端用戶發出告警信號。

近年來，國內外學者在原子鐘異常數據的分析和處理上做了許多有意義的工作，提出很多有效的原子鐘異常探測算法，如最小二乘法、阿倫方差法、基于中位數(MAD)的抗差估計法、卡爾曼濾波器法以及最大似然比檢驗法等[2-7]。但是這些方法在實時性上都存在一些不足，其中，最小二乘法、阿倫方差法和基于中位數的抗差估計法主要用于對事后鐘差數據進行處理，而卡爾曼濾波器法和最大似然比檢驗法的參數選取比較復雜。目前，我國對BDS星載原子鐘完好性的研究相對較少，值得進一步研究。為了保證實時導航定位用戶的需求，必須對原子鐘進行實時的異常識別、解釋和處理。

已有研究表明，原子鐘的頻率數據服從高斯白噪聲分布，其均值和方差與時間有關[4]。因此，本文基于BDS星載原子鐘的時頻特性，結合廣義最大似然比估計和鐘差多項式模型，提出了一種便捷可靠的實時BDS原子鐘異常探測與處理算法。仿真結果表明，該算法能夠有效地對頻率跳變異常進行探測和處理。

2　基于廣義似然比檢驗的原子鐘異常實時探測與處理算法

為了實現對原子鐘異常的實時探測與處理，本文對廣義似然比檢驗進行了改進，提出了一種快速的閾值選取方法；然后利用該方法進行原子鐘異常的探測，并利用實時的鐘差預報模型計算出異常點的估計值，從而實現對原子鐘異常的修復。

2.1廣義似然比檢驗(GLRT)

廣義似然比檢驗是一種快速有效地探測均值和方差變化的方法，文獻[4]的研究表明，正常工作模式下原子鐘的頻率數據服從高斯白噪聲分布。因此，廣義似然比檢驗能滿足實時原子鐘異常探測的要求。

2.1.1GLRT基本原理

基于原子鐘頻率數據服從高斯白噪聲分布的前提，可以用兩個假設H0和H1來描述其模型參數的變化。H0(零假設)表明樣本數據y[n],n=0,…,N-1服從于均值為μ0、方差為σ0的高斯概率密度函數；H1(備擇假設)表明樣本數據的均值和方差發生了改變，均值變為μ1或者方差變為σ1。上述假設可以用下面的等式來描述：

(1)

式中，n0為未知的跳變點，需要對其進行估計。

因此，利用最大似然估計法判斷H1假設成立的準則[4,5]為：

T(y)=log(LG(y))=log(p(y;[μ0,σ0,μ1,σ1,n0],H1))-log(p(y;[μ0,σ0],H0))>γ

(2)

式中，LG(y)為似然比，p(y;[μ0,σ0],H0)和p(y;[μ0,σ0,μ1,σ1,n0],H1)分別為假設H0和H1成立下序列的概率分布，γ為判斷零假設成立與否的閾值。為了便于計算，T(y)取似然比的對數形式。

將式(2)中的未知參數用其最大似然估計值代替，則式(2)可化為：

(3)

需要指出，為了保證算法的實時性，僅對序列的最后一個值進行異常探測，同時，需要將n0值固定。因此，當取定一個γ值后，即可進行異常的判斷。

2.1.2GLRT閾值γ的選取

在進行GLRT方法計算的過程中，關鍵之處在于γ值的選取。常用的方法為通過計算ROC曲線來選取[4,8]，這需要計算正確檢測的概率(PD)以及發出錯誤警告的概率(PFA，即H0假設成立時選擇H1假設)。這種常用的判斷γ的方法雖然能夠確定一個可靠而且精確的閾值，但是選取閾值時需要進行蒙特卡洛實驗，計算量較大，不適合實時異常探測。因此，本文提出了一種新的快速確定γ取值的方法，其步驟如下。

(1)已知數據分段。將已知無異常頻率數據Y分為三部分Y1、Y2和Y3，并假設Y3為“待探測”數據，則數據結構如圖1所示。

圖1　數據結構示意圖

從上圖可以看出，Y1和Y2均服從于同一正態分布N(μ0,σ0)。在已知序列Y2中加入待探測數據y[N]后，若y[N]為無異常數據，Y2仍服從正態分布N(μ0,σ0)，否則Y2將不再服從N(μ0,σ0)。

γ=max(T(y))

(4)

2.2異常處理算法

常見的原子鐘模型有很多，如多項式模型、灰色系統模型、卡爾曼濾波器模型等[9]?？紤]到算法的復雜性與計算量，本文采用最簡單且短期預報精度高的多項式模型對已知數據進行建模。不同星載原子鐘的頻率漂移特性有所不同，故針對星載銫鐘和銣鐘的多項式模型可表示為：

(5)

式中，a、b、c為多項式模型的未知參數，ε為白噪聲。

通過最小二乘擬合求解出多項式模型的參數，再利用該模型進行預報。需要指出，對原子鐘異常的探測常常是基于頻率數據，而對原子鐘的預報采用鐘差(即相位)數據，故需要通過一次差分將鐘差數據轉換為頻率數據。

2.3算法描述

基于前文所述，本文設計了一種基于GLRT方法的實時原子鐘異常探測和基于鐘差多項式模型預報的異常處理算法。其主要步驟如下：

(1)數據截取。從無異常的鐘差歷史序列中截取序列X及其相對應的頻率序列。

(2)確定閾值。利用頻率序列Y計算得到T(y)序列，取其最大值為GLRT探測器的判斷閾值。

(5)異常警告。重復上述步驟，探測下一個待探測頻率數據點；當連續出現異常警告次數超過N/2時，說明原子鐘發生了持久性的頻率跳變，需要重新選擇鐘差序列X和頻率序列Y。

在上述流程中，需要注意幾點：(1)鐘差序列X僅僅在出現異常值時作為已知數據建立鐘差模型，GLRT探測時使用的是頻率數據(一階差分)；(2)在新的鐘差數據輸入時，只需要對其進行一次差y[i]=x[i+1]-x[i]即可。

3　算例驗證與分析

在星載原子鐘常見異常現象中，信號缺失是由于原子鐘發生故障或徹底喪失信號輸出功能，而頻率漂移在多數情況下可以人工建模預測并修正。星載原子鐘單純的相位跳變很少發生，通常相位跳變是頻率跳變的間接反映，因此最常見的原子鐘異常為頻率跳變[10]。本節基于武漢大學分析中心提供的BDS精密鐘差數據，驗證了原子鐘頻率跳變探測與處理的有效性。

3.1數據設計與仿真

頻率跳變通常分為暫時性跳變和持久性跳變。在實時原子鐘異常探測的情況下，持久性跳變可以看作暫時性跳變的累積。因此，本文設計了BDS原子鐘兩種情況的暫時性頻率跳變：(1)單點頻率跳變；(2)連續多點頻率跳變。然后利用本文算法進行了異常探測與處理，以此驗證算法的有效性。

本文使用的樣本數據是武漢大學分析中心提供的2013年1月8日采樣周期為5分鐘的精密鐘差數據，選取BDSC01號衛星作為研究對象，其鐘差和頻率數據如圖2所示。

圖2　BDS C01星載原子鐘的鐘差與一階差分(頻率)序列

從圖2可以看出，其頻率數據基本上滿足高斯白噪聲分布，證實了本文算法的前提條件，即星載原子鐘的頻率數據服從高斯分布。由于原子鐘出現異常是小概率事件，因此，難以得到帶有異常事件的鐘差數據，本文采用文獻[11]的方法，即在正常鐘差數據的疊加異常數據表征異常。在該頻率序列n=150、n=210、n=211、n=212處分別加入人為的“暫時性跳變”，該異常大小為均值的2、3、3、3倍，加入異常后的頻率序列如圖3所示。

圖3　在原始頻率數據中加入人為異常后的頻率序列

3.2算例結果

首先，從歷史無異常鐘差序列中截取長度為140的序列進行閾值確定。令滑動窗口長度為116，其中子窗口Y1、Y2長度分別為100和15(固定n0=100)，“待檢測”序列長度為25。利用Y1和Y2對Y3進行檢驗計算得到T(y)，將窗口滑動25次，取T(y)中最大值作為閾值，則有γ=9.475。

其次，繼續滑動窗口對未知數據逐一進行GLRT探測。首先探測到n=150的異常點，T150(y)=262.05，T150(y)>γ，因此判斷該點為頻率跳變點，發出警告，如圖4所示。為了便于對比，圖中給出了對全部待探測值進行計算的T(y)序列。

圖4　第一次異常警告

為了進行比較，計算異常點附近幾個點的T(y)值，如表1所示。

表1異常點及其附近點的T(y)值(γ=9.475)

點號148149150151152T(y)0.501.05262.05211.20211.54

從上表可以看出，T150(y)遠遠大于閾值γ，并且在異常未經處理的情況下，T150(y)之后的T(y)的值也大于閾值。由此可見，一個異常數據影響了整個序列的分布。

發現異常數據后，實時對異常進行修復處理。 利用已知正常鐘差序列建立多項式模型，并外推得到該點的估計值，從而對異常點進行修復。為了驗證是否對該異常進行了有效的修復，再一次進行GLRT探測。從圖5可以看出，這時在y[150]處并無異常，不發出警告信號。

圖5　修復第1個異常點后的T(y)序列

在修復完y[150]處的異常后，繼續滑動窗口，進行GLRT探測。此段的T(y)值如圖6所示。從圖6可以看出，當檢測到y[210]時，此時T(y)=343.4，T(y)>γ，因此判斷該點為頻率跳變點，發出警告。

圖6　第二次異常警告

對其進行修復后，再滑動窗口，可以檢測出連續的三個頻率跳變點。對這些異常點進行修復，平均不符值大小為3.27×10-10s。表2列出了異常處理前后的T(y)值。從表中可以看出，GLRT探測出了連續的三個頻率跳變點，并通過多項式模型對這三個異常點進行了及時和有效地修復，并將誤差控制在10-10s量級。

表2　異常點及其附近點的T(y)值

圖7　修復所有異常點后的T(y)序列

4　結　論

本文給出了一種基于GLRT探測和鐘差預報修復的實時原子鐘異常探測和處理算法，通過兩種頻率跳變(單點跳躍和連續多點跳躍)的計算可以看出，該算法不僅能夠有效地探測常見的頻率跳變(暫時性和持久性頻率跳變)，而且能及時對異常點進行警告和修復處理。該算法簡單有效，自動化程度較高，對實際工程應用有一定參考意義。需要指出的是，本文算法的前提條件為存在一段已知的無異常數據，因此在不滿足這一前提的條件下，如何啟動該數據處理流程是下一步的研究重點。此外，本文算法中的參數可能需要根據不同數據的特性做出相應的調整，以更好地對異常進行探測和修復。

[1]MORIKAWA T. Spaceborne Atomic Clock [J]. IEIC Technical Report (Institute of Electronics, Information and Communication Engineers), 2001, 101(157): 75-81.

[2]郭海榮. 導航衛星原子鐘時頻特性分析理論與方法研究[D]. 鄭州：信息工程大學，2006.

[3]黃新明，龔航，朱祥維.基于Kalman濾波器的衛星鐘實時異常監測算法[J]. 宇航計測技術，2011，31(5): 6-11.

[4]Nunzi E, Galleani L, Tavella P, et al. Detection of anomalies in the behavior of atomic clocks[J]. Instrumentation and Measurement, 2007, 56(2): 523-528.

[5]Nunzi E, Carbone P. Monitoring signal integrity of atomic clocks by means of the GLRT [J]. Metrologia, 2008, 45(6): 103.

[6]Galleani L, Tavella P. Tracking nonstationarities in clock noises using the dynamic Allan variance[C].Frequency Control Symposium and Exposition, 2005.

[7]Weiss M, Shome P, Beard R. On-board GPS clock monitoring for signal integrity[C]. Proceedings of the 42nd Annual Precise Time and Time Interval (PTTI) Applications and Planning Meeting, 2010.

[8]Nunzi E, D'Ippolito D. A novel theoretical analysis of fault detection for atomic clock[C].Advanced Methods for Uncertainty Estimation in Measurement,2009.

[9]崔先強，焦文海.灰色系統模型在衛星鐘差預報中的應用[J]. 武漢大學學報·信息科學版, 2005, 30(5): 447-450.

[10]唐升，劉婭，李孝輝.星載原子鐘自主完好性監測方法研究[J].宇航學報，2013, 34(1): 39-45.

[11]牛飛. GNSS完好性增強理論與方法研究[D]. 鄭州: 信息工程大學, 2008.

A Real-time BDS Satellite Clock Anomaly Detection and Processing Algorithm

Zhao Dan，Zhou Peiyuan，Du Lan

Institute of Navigation and Aerospace Engineering, Information Engineering University, Zhengzhou 450001, China

Spaceborne atomic clock anomaly can be detected quickly by continuously monitoring of clock bias, thus it ensures the integrity of satellite navigation system. This paper puts forward a new method of real-time atomic clock anomaly detection and processing combined with principles of Generalized Likelihood Ratio Test (GLRT) and polynomial clock model. Firstly, the real time anomaly detection is done based on the GLRT. Then the short-term clock bias is forecasted and extrapolated by using the polynomial clock model, and on this base, the anomaly points are eliminated and replaced by extrapolation values. This paper simulates the examples based on BDS satellite clock bias provided by Wuhan University Analysis Center. The results show that the method can detect and handle atomic clock anomalies effectively (especially abnormal frequency jumps), and it can guarantee the integrity of BDS atomic clock data.

BDS satellite clock; atomic clock anomaly detection; GLRT detector; clock bias model

2015-07-02。

國家自然科學基金資助項目(41174025)。

趙丹(1986—)，女，碩士研究生，主要從事導航數據處理算法研究。

P236

A