伍貽威

北京衛(wèi)星導(dǎo)航中心，北京 100094

守時試驗室和全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)都需要利用多臺原子鐘組成鐘組，建立和保持一個時間基準[1-4]。守時本質(zhì)上是個預(yù)測問題，即需要通過原子鐘的歷史表現(xiàn)來預(yù)測其未來表現(xiàn)[5-8]。國際上已經(jīng)達成共識[8]：一個好的鐘應(yīng)該是一個可預(yù)測的鐘(a good clock is a predictable clock)。可預(yù)測性的一個重要體現(xiàn)是原子鐘信號沒有異常[8]。建立一個實時的、可預(yù)測的紙面時間，需要在異常發(fā)生后盡可能短的時間內(nèi)檢測出異常；根據(jù)異常的嚴重等級，合理降低其權(quán)重甚至剔除出鐘組，盡可能降低異常對于紙面時間的影響[8]。

頻率異常檢測的性能是與原子鐘的噪聲水平高度相關(guān)的。典型氫鐘的閃爍底(flicker floor，即Allan偏差最小值)一般在10-15～10-16量級，而典型銫鐘的閃爍底只能達到10-14～10-15量級。這就導(dǎo)致相對于氫鐘，銫鐘10-14量級的頻率異?？赡苎蜎]在噪聲之中，從而更不容易檢測，因而其對紙面時間的影響也更大。

國內(nèi)外檢測頻率異常的研究方法包括：動態(tài)Allan方差法[9-11]、時-頻譜分析法[9,12]、假設(shè)檢驗法[11,13-16]。動態(tài)Allan方差法對微小頻率異常的反應(yīng)時間較長；時-頻譜分析法比較適用于檢測時變信號或是周期性波動等在頻譜上出現(xiàn)明顯譜線的異常；假設(shè)檢驗法較適用于頻率異常的檢測。

假設(shè)檢驗法的核心是設(shè)計合理的檢驗統(tǒng)計量和相應(yīng)的檢測門限值。檢測門限值是由檢驗統(tǒng)計量的統(tǒng)計特性以及虛警概率(PFA)決定的。進一步，可以推導(dǎo)出檢測概率(PD)和PFA、平均頻率跳變幅度(|Ya|)之間的函數(shù)關(guān)系。一般認為檢驗統(tǒng)計量的統(tǒng)計特性是先驗已知的，這在大部分情況下都是可行的。因為對于某一臺原子鐘，其噪聲的統(tǒng)計特性一般不會發(fā)生很大的改變。當檢驗統(tǒng)計量的統(tǒng)計特性未知，可以采用廣義似然比(generalized likelihood ratio test，GLRT)的方法，在檢驗統(tǒng)計量的分布函數(shù)中用最大似然估計結(jié)果代替未知參數(shù)。因此，GLRT法[11]可以看成是假設(shè)檢驗法的特例。

目前的研究中，文獻[13]只分析了不同PFA的試驗效果，也就是只進行了第一步，沒有從理論上進一步推導(dǎo)不同|Ya|對應(yīng)的PD；包括采用GLRT法在內(nèi)的大部分假設(shè)檢驗法的檢驗統(tǒng)計量都是頻差[11,14-16]，例如文獻[14—15]采用的檢驗統(tǒng)計量分別為Kalman濾波器輸出的頻差的一步和多步預(yù)測誤差，沒有充分利用頻率異常在時差上的積累效應(yīng)。

本文基于文獻[17—18]等描述的典型銫鐘的鐘差預(yù)測不確定度的解析表達式，將時差的預(yù)測誤差作為檢驗統(tǒng)計量，這一點與文獻[11，14—16]有所不同。本文認為原子鐘的噪聲特性可以先驗得到，所以在算法中是已知的，不需要再采用GLRT法來估計參數(shù)。進一步，按照文獻[19—20]的思路，結(jié)合預(yù)測誤差的概率分布函數(shù)，推導(dǎo)了不同PFA對應(yīng)的檢測門限和PD的計算公式，并總結(jié)了提高PD的方法。根據(jù)本文的理論分析結(jié)論，只要給出PFA和|Ya|的值，即可計算得到PD，這為銫原子鐘的頻率異常檢測提供了理論支撐。后續(xù)還可以將本文算法的性能與大部分采用頻差作為檢驗統(tǒng)計量的算法[11,14-16]進行對比分析，進一步優(yōu)化算法設(shè)計。

1 基本原理

1.1 原子鐘信號

典型銫鐘的時差表示為[21-24]

X1(t)=x0+y0t+σ1W1(t)+

(1)

式中，等號右側(cè)前2項為一次多項式表示的確定性分量，后2項為隨機性分量，即原子鐘噪聲；X1代表原子鐘時差；x0和y0分別代表時差和頻差的初值；W1(t)和W2(t)分別代表2個獨立的維納過程，并且有W(t)～N(0,t)，即每個維納過程服從數(shù)學(xué)期望為0，方差為時間t的正態(tài)分布；σ1和σ2分別是這2個維納過程的擴散系數(shù)(diffusion coefficients)，用于表明噪聲的強度；第3項代表調(diào)頻白噪聲(white frequency modulation noise，WFM)；第4項代表調(diào)頻隨機游走噪聲(walk random frequency modulation noise，RWFM)，即W1(t)、W2(t)的積分在X1上分別表現(xiàn)為WFM和RWFM。

WFM和RWFM都是有色噪聲[25]，一般采用Allan方差來表征頻率穩(wěn)定度[26]。Allan方差表達式為[21-23]

(2)

式中，τ為平滑時間。式(2)等號右側(cè)第1項為WFM分量，第2項為RWFM分量，表明WFM和RWFM在對數(shù)Allan方差圖中的斜率分別為-1和1。

式(1)中的原子鐘噪聲服從正態(tài)分布，數(shù)學(xué)期望為0，方差為[19,21]

(3)

式中，等號右側(cè)第1項是WFM的方差；第2項是RWFM的方差，它們都隨時間的變大而變大。觀察式(2)和式(3)，得到

(4)

圖1展示了這100臺仿真銫鐘的時差(綠色曲線)，以及式(3)的1倍和2倍平方根計算得到的1σ和2σ噪聲標準差(黑色加粗曲線)。從圖1清楚地看出噪聲標準差隨時間的變大而變大。

1.2 鐘差預(yù)測基本原理

觀測鐘差可以認為是時差疊加上觀測噪聲，表示為[17-18]

Z(tk)=X1(tk)+σ·ξ(tk)

(5)

式中，Z(tk)是在tk時刻的觀測鐘差；ξ(tk)是在tk時刻的觀測噪聲，一般認為是調(diào)相白噪聲(white phase modulation noise，WPM)；σ代表觀測噪聲的強度。

(6)

預(yù)測誤差由[預(yù)測值-真實值]來表示。由式(1)和式(6)，預(yù)測誤差表示為

(7)

(8)

[W1(t0)-W1(t0-T)]+

(9)

維納過程具有以下性質(zhì)

cov[W(s),W(t)]=min(s,t)

當t≥s時

式中，cov代表協(xié)方差；min代表最小值。當t=s時，協(xié)方差實際上就是方差。

根據(jù)維納過程和白噪聲的性質(zhì)，推導(dǎo)得到

(10)

(11)

(12)

于是得到

(13)

同理，根據(jù)維納過程的性質(zhì)，推導(dǎo)得到

(14)

(15)

將式(13)—(15)代入式(8)，得到

(16)

預(yù)測不確定度u(tp)為式(16)的平方根。1σ預(yù)測不確定度和2σ預(yù)測不確定度從統(tǒng)計意義上分別代表了預(yù)測誤差的68%和95%的置信區(qū)間。

1.3 觀測間隔選取

本地的鐘差測量不確定度可做到優(yōu)于0.1 ns(σ2=1×10-20s2)，當預(yù)測時間tp不是太短時，其影響基本可以忽略。假設(shè)σ2=0，將式(2)代入式(16)，得到

(17)

式中，u2(tp)包含2項；對于給定的tp，第2項的值不變；當選擇觀測間隔T等于Allan方差曲線(一般為V字形)取最小值對應(yīng)的平滑時間值(記為Tmin)時，第1項取最小值，于是u2(tp)取最小值。對于與圖1相同參數(shù)的典型銫鐘，使第1項取最小值的T約等于100 d。

圖1 100臺仿真銫鐘的時差、1σ和2σ噪聲標準差Fig.1 Time differences and 1σ and 2σ noise standard deviation of the 100 simulated cesium clocks

圖2 鐘差預(yù)測Fig.2 Diagram of clock prediction

從圖3得出結(jié)論：預(yù)測時間較短即當tp?T時，u(tp)僅略大于σy(tp)tp；這時預(yù)測性能主要由原子鐘噪聲決定。

從直觀上理解：假設(shè)時差和頻差的估計誤差都一直為零，這相當于估計值與真實值完全相同(實際上是不可能的)，此時按照式(7)，如果忽略觀測噪聲，預(yù)測誤差完全就是原子鐘噪聲，噪聲方差隨時間變大，如圖1所示。

圖3 預(yù)測不確定度及其前3項和第4項平方根的值Fig.3 Prediction uncertainty and its first three terms and the fourth term values

表1 不同T和tp對應(yīng)的u(tp)

從表1可以看出，當tp?T且T≤Tmin時，u(tp)僅略大于σy(tp)tp。

2 頻率異常檢測算法

2.1 基本思路

通過對比分析預(yù)測誤差和預(yù)測不確定度，提供了一種檢測原子鐘頻率異常檢測的方法[19]，即當異常發(fā)生在未來時間段[t0,t0+tp]時，在某個t時刻(t>t0)及之后時間段的預(yù)測誤差一直大于某個門限值，認為在t時刻頻率發(fā)生了跳變。

本節(jié)首先選取門限值等于3σ預(yù)測不確定度的理論值為示例進行說明，再推廣到一般情況。預(yù)測誤差落在±1σ、±2σ和±3σ預(yù)測不確定度范圍內(nèi)的理論概率分別為68.26%、95.44%和99.74%[28]。從理論上講，預(yù)測誤差大于3σ預(yù)測不確定度的概率很小，僅為1-99.74%=0.26%。

2.2 頻率跳變檢測方法

本文方法等價于做如下二元假設(shè)檢驗。

H0——未發(fā)生頻率跳變。

H1——發(fā)生頻率跳變。

選取預(yù)測誤差ε(tp)作為檢驗統(tǒng)計量。首先需要確定在這兩種假設(shè)下，檢驗統(tǒng)計量的數(shù)學(xué)分布。

H0情況下，顯然預(yù)測誤差服從N(0,u2(tp))分布。

H1情況下，假設(shè)跳變發(fā)生在[t0,t0+tp]的中間時刻，記為t0+tp-Tjump時刻，其中0

(18)

所以，在[t0,t0+tp]時間段內(nèi)積累的時差為

xjump(T)=Yatp=Y0Tjump

(19)

這時預(yù)測誤差服從N(Yatp,u2(tp))分布。

綜上，在這兩種假設(shè)下統(tǒng)計量的數(shù)學(xué)分布為

(20)

將H0為真時判H1成立的概率稱為虛警概率，用PFA表示[25]。把H1為真時判H1成立的概率稱為檢測概率，用PD表示[25]。PFA的值也被稱為顯著性水平[28]。

本文通過約束PFA，求解不同Ya對應(yīng)的PD。

如前文所述，ε(tp)落在±u(tp)、±2u(tp)和±3u(tp)范圍內(nèi)的概率分別為68.26%、95.44%和99.74%，對應(yīng)的PFA分別為1-68.26%=31.74%、1-95.44%=4.56%和1-99.74%=0.26%[28]。

以約束PFA=0.26%為例。直觀上理解，假如|Ya|tp=3u(tp)，顯然有PD=50%；那么當|Ya|tp>3u(tp)時，PD>50%；當|Ya|tp<3u(tp)時，PD<50%。因此，當約束PFA時，PD是與|Ya|正相關(guān)的；|Ya|越大，PD越大。

具體而言，以約束PFA=0.26%為例，根據(jù)式(20)，推導(dǎo)出檢測概率的表達式為

(21)

由式(21)可以求解得到任意Ya對應(yīng)的PD。

2.3 示例分析和更一般的情況

從圖4中看出，PD隨著|Ya|的增大而增大；當|Ya|=3u(tp)/tp=7.26×10-14時，PD=50%(圖中圓圈)，即直觀認識與式(21)計算結(jié)果是符合的。

服從標準正態(tài)分布的隨機變量X的分布函數(shù)用φ(x)來表示[28]。幾個關(guān)鍵數(shù)值為φ(-3)=0.001 3，φ(-2)=0.022 8，φ(-1)=0.158 7，φ(0)=0.5，φ(1)=0.841 3，φ(2)=0.977 2，φ(3)=0.998 7[28]。

圖4 PD-|Ya|曲線Fig.4 PD-|Ya| curve

前文所述的“ε(tp)落在±u(tp)、±2u(tp)和±3u(tp)范圍內(nèi)的概率分別為68.26%、95.44%和99.74%”，對應(yīng)于φ(1)-φ(-1)=68.26%、φ(2)-φ(-2)=95.44%和φ(3)-φ(-3)=99.74%[28]。

顯然，在H1情況下，統(tǒng)計量ε(tp)的分布函數(shù)為

(22)

所以，對于該實例，在約束PFA=0.26%時，不同|Ya|對應(yīng)的PD如表2所示。

表2 不同|Ya|對應(yīng)的PD

由圖4和表2看出：當|Ya|>4u(tp)/tp時，可以比較有效地檢測出異常；當|Ya|>1.05×10-13時，檢測概率在90%以上；當|Ya|較小(例如<3u(tp)/tp)時，檢測效果并不顯著。

上述示例中約束PFA=0.26%，對應(yīng)的ε(tp)門限值為±3u(tp)。在更一般的情況下，可以約束PFA=5%、PFA=1%、PFA=0.1%等。這時，對應(yīng)的ε(tp)門限值±γ需要查詢分布函數(shù)φ(x)的函數(shù)表來確定，具體方式是使φ(γ)-φ(-γ)=1-PFA。例如，當約束PFA=5%時，查閱φ(x)的函數(shù)表，得到γ=1.96u(tp)。然后，需要將式(21)中積分下限(上式)或上限(下式)分別修改為γ和-γ，得到式(23)。根據(jù)式(23)，可以求解得到任意幅度|Ya|對應(yīng)的PD值。

(23)

3 仿真分析、實測數(shù)據(jù)分析

3.1 仿真分析

圖5和圖6分別給出了在Ya=3u(tp)/tp和4u(tp)/tp的情況下時，兩次各10 000臺銫鐘的預(yù)測誤差ε(tp)。從圖5和圖6中明顯看出在t0+tp/2時刻發(fā)生的頻率跳變。

分別統(tǒng)計在這兩種情況下ε(tp)落在±3u(tp)外的銫鐘數(shù)量，作為本文方法檢測出頻率跳變異常的數(shù)量。在這兩種情況下，能檢測出頻率跳變異常的數(shù)量分別為5067臺和8502臺。

表3列出了在這兩種情況下，理論PD值、實際檢測數(shù)量和相應(yīng)的試驗PD值。

圖5 Ya=3u(tp)/tp時10 000臺銫鐘的預(yù)測誤差Fig.5 Prediction errors of 10 000 cesium clocks when Ya=3u(tp)/tp

圖6 Ya=4u(tp)/tp時10 000臺銫鐘的預(yù)測誤差Fig.6 Prediction errors of 10 000 cesium clocks when Ya=4u(tp)/tp

表3說明，試驗PD值和理論PD值基本符合。

綜上，仿真試驗驗證了理論分析的結(jié)論。

表3 不同Ya對應(yīng)的理論和試驗PD值

3.2 實測數(shù)據(jù)分析

圖7 銫鐘鐘差和鐘差的一次多項式殘差Fig.7 Cesium clock difference and its one-order residual

人為在第39.5 d處引入幅度為8×10-14的頻率跳變。共進行20次預(yù)測，都選取T=20 d和tp=3 d。第1次預(yù)測從第20 d開始，之后每次預(yù)測開始時間往后移1 d，最后1次預(yù)測開始時間為第39 d。共獲得20組預(yù)測誤差曲線。選取γ=3u(tp)，即約束PFA=0.26%。

圖8給出了這20組預(yù)測誤差曲線以及3σ預(yù)測不確定度曲線，其中第20組為紅色(當tp等于1 d、2 d和3 d時標注圓圈)，第19組為藍色(當tp等于1 d、2 d和3 d時標注方塊)，其余18組為綠色。從圖8看出：對于第20組，相當于在tp=0.5 d時發(fā)生了異常，當tp=1 d時未檢測出異常，tp=2 d和tp=3 d時成功檢測出了異常；對于第19次預(yù)測，相當于在tp=1.5 d時發(fā)生了異常，當tp=2 d時未檢測出異常，當tp=3 d時成功檢測出了異常。

圖8 20組預(yù)測誤差曲線以及3σ預(yù)測不確定度曲線Fig.8 20 prediction error curves and 3σ prediction uncertainty curve

表4列出了對于第20組和第19組，當tp分別等于1 d、2 d和3 d時，按照式(16)計算得到的γ、按照式(18)計算得到的Ya、按照式(23)計算得到的理論PD值、實際檢測情況(圖8)。

表4 不同Ya對應(yīng)的理論PD值和檢測情況

綜上，實測數(shù)據(jù)試驗驗證了理論分析的結(jié)論。根據(jù)本文的理論分析和仿真試驗，還可以進一步得到以下結(jié)論。

(1)PD和PFA、|Ya|之間存在函數(shù)關(guān)系，約束PFA、|Ya|之中任意一個量的值，可以得到PD和另一個量的函數(shù)曲線。

圖4給出了PD-|Ya|曲線。實際上，當約束|Ya|時，還可以畫出PD-PFA曲線，這就是著名的接收機工作特性(receiver operating characteristic，ROC)曲線，分析在不同|Ya|情況下不同PFA值對應(yīng)的PD值，可以更合理地選取PFA的值。

(2) 根據(jù)ROC曲線，對于相同的|Ya|，假如PFA值更大，相應(yīng)的PD也更大。

例如，當T=20 d和tp=1 d時，根據(jù)表2，假如將PFA值由0.26%提高至4.56%，相當于將檢測門限γ由±3u(tp)降低至±2u(tp)，原來|Ya|=3u(tp)/tp=7.26×10-14時，才有PD=50%，現(xiàn)在只需要|Ya|=2u(tp)/tp=4.84×10-14時，就有PD=50%，而當|Ya|依然為3u(tp)/tp=7.26×10-14時，PD=84.13%。

(3) 當tp?T且T≤Tmin時，對于相同的PFA，假設(shè)|Ya|保持不變，增加tp，可以有效增大PD。

例如：當T=20 d，約束PFA=0.26%時，如果tp=1 d，當Ya=3u(1 d)/1 d=7.26×10-14時，有PD=50%；如果tp=2 d，當Ya依然為7.26×10-14時，根據(jù)式(21)計算，有PD=87.36%；如果tp=2 d，此時只要Ya=3u(2 d)/2 d=5.26×10-14，即有PD=50%。從直觀上理解：預(yù)測時間tp越長，由頻率跳變引起的積累時差越大，即式(20)中的|Ya|tp變大，所以PD增大。

(4) 進一步，當tp?T且T≤Tmin時，對于同樣的Y0，增加tp，使得|Ya|也相應(yīng)地變大，|Ya|和tp同時增大使PD增大的效果更顯著。

例如，依然采用上述情況，假設(shè)在t0+12 h時刻發(fā)生頻率跳變，幅度為Y0=6u(1 d)/1 d=1.45×10-13，如果tp=1 d，按照(18)計算得到Y(jié)a=Y0/2=7.26×10-14，這時根據(jù)式(21)計算，有PD=50%；如果tp=2 d，按照(18)計算得到Y(jié)a=3Y0/4=1.14×10-13，這時根據(jù)式(21)計算，有PD=99.93%，即基本可以檢測出異常。所以，對于同樣的Y0，增加tp，相當于使式(20)中的|Ya|和tp都增大，所以|Ya|tp將變得更大，從而PD更高。

4 結(jié)束語

本文提出一種基于鐘差預(yù)測的銫原子鐘頻率異常的檢測算法。本算法采用假設(shè)檢驗的方法，檢測概率(PD)和虛警概率(PFA)、平均頻率跳變幅度(|Ya|)之間存在函數(shù)關(guān)系。為了充分利用頻率異常體現(xiàn)在時差上的積累效應(yīng)，本算法將時差的預(yù)測誤差作為檢驗統(tǒng)計量，對檢驗統(tǒng)計量做二元假設(shè)檢驗。在兩種情況下，檢驗統(tǒng)計量都服從正態(tài)分布。通過約束PFA，查詢正態(tài)分布的分布函數(shù)，得到對應(yīng)的檢測門限值。根據(jù)在異常情況下的檢驗統(tǒng)計量分布函數(shù)和檢測門限值，推導(dǎo)給出了求解|Ya|對應(yīng)的PD的表達式。仿真試驗和實測數(shù)據(jù)試驗都驗證了理論分析的結(jié)論。

在實際應(yīng)用中，由于算法已經(jīng)約束了PFA，一旦算法給出報警，那么極有可能發(fā)生了頻率異常(發(fā)生異常的概率為1-PFA)，此時需要進行人工干預(yù)。|Ya|越大，PD也就越高。根據(jù)本文的分析結(jié)論，提高PD的方法有：①提高PFA(相當于降低γ，盡管會增加系統(tǒng)報警的頻度)；②增加tp(當tp?T且T≤Tmin時，相當于同時增大了|Ya|和tp)。本文成果還可以應(yīng)用于精密定軌與時間同步領(lǐng)域[29]。