伍貽威

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

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

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

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

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

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

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

1 基本原理

1.1 原子鐘信號(hào)

典型銫鐘的時(shí)差表示為[21-24]

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

(1)

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

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

(2)

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

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

(3)

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

(4)

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

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

觀測(cè)鐘差可以認(rèn)為是時(shí)差疊加上觀測(cè)噪聲，表示為[17-18]

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

(5)

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

(6)

預(yù)測(cè)誤差由[預(yù)測(cè)值-真實(shí)值]來(lái)表示。由式(1)和式(6)，預(yù)測(cè)誤差表示為

(7)

(8)

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

(9)

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

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

當(dāng)t≥s時(shí)

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

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

(10)

(11)

(12)

于是得到

(13)

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

(14)

(15)

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

(16)

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

1.3 觀測(cè)間隔選取

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

(17)

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

圖1 100臺(tái)仿真銫鐘的時(shí)差、1σ和2σ噪聲標(biāo)準(zhǔn)差Fig.1 Time differences and 1σ and 2σ noise standard deviation of the 100 simulated cesium clocks

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

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

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

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

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

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

2 頻率異常檢測(cè)算法

2.1 基本思路

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

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

2.2 頻率跳變檢測(cè)方法

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

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

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

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

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

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

(18)

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

xjump(T)=Yatp=Y0Tjump

(19)

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

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

(20)

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

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

如前文所述，ε(tp)落在±u(tp)、±2u(tp)和±3u(tp)范圍內(nèi)的概率分別為68.26%、95.44%和99.74%，對(duì)應(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%；那么當(dāng)|Ya|tp>3u(tp)時(shí)，PD>50%；當(dāng)|Ya|tp<3u(tp)時(shí)，PD<50%。因此，當(dāng)約束PFA時(shí)，PD是與|Ya|正相關(guān)的；|Ya|越大，PD越大。

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

(21)

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

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

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

服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量X的分布函數(shù)用φ(x)來(lái)表示[28]。幾個(gè)關(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%”，對(duì)應(yīng)于φ(1)-φ(-1)=68.26%、φ(2)-φ(-2)=95.44%和φ(3)-φ(-3)=99.74%[28]。

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

(22)

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

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

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

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

(23)

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

3.1 仿真分析

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

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

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

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

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

表3說(shuō)明，試驗(yàn)PD值和理論P(yáng)D值基本符合。

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

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

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

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

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

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

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

表4列出了對(duì)于第20組和第19組，當(dāng)tp分別等于1 d、2 d和3 d時(shí)，按照式(16)計(jì)算得到的γ、按照式(18)計(jì)算得到的Ya、按照式(23)計(jì)算得到的理論P(yáng)D值、實(shí)際檢測(cè)情況(圖8)。

表4 不同Ya對(duì)應(yīng)的理論P(yáng)D值和檢測(cè)情況

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

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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

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