惠 恬，趙高長，蘇 軍，龔瑩瑩

(西安科技大學 理學院，陜西 西安 710600)

0 引言

精準的時間尺度需要根據(jù)參與計算的原子鐘模型及其頻率穩(wěn)定度來設計[1-2]，由于銫原子鐘自身存在的噪聲會對其輸出信號的頻率穩(wěn)定度和精度造成影響，因此在計算時間尺度之前，對銫原子鐘數(shù)據(jù)進行消噪處理以提高其頻率穩(wěn)定度，具有重要意義。

目前，可以通過濾波的方法對原子鐘的噪聲進行處理。例如，文獻[3]應用小波方法、Vondrak方法和卡爾曼(Kalman)濾波等7種方法，對原子鐘數(shù)據(jù)進行消噪處理，均達到了比較好的降噪效果，其中小波方法的效果最好。在實際應用中，Kalman濾波[4-5]和小波方法[6-7]都有一定的缺陷。Kalman濾波的過程過于復雜，小波分析則需要選擇一個合適的小波函數(shù)，分解層數(shù)還要選擇合適的閾值。經(jīng)驗模態(tài)分解(empirical mode decomposition，EMD)方法是一種新的自適應信號處理技術，采用該方法對原子鐘數(shù)據(jù)進行降噪處理，可以消除隨機噪聲的影響，提高原子鐘的短期穩(wěn)定度，獲得更加精確的時間尺度，同時避免了Kalman濾波和小波閾值方法的缺點。文獻[8]采用EMD方法，對高性能5071A原子鐘信號進行了去噪和頻率預測。文獻[9-10]在文獻[8]的基礎上，采用EMD方法對銫原子鐘的白色調(diào)頻噪聲進行去噪處理，并應用經(jīng)過處理后的原子鐘數(shù)據(jù)進行時間尺度計算，研究了EMD方法對合成頻率穩(wěn)定度的影響。然而，傳統(tǒng)的EMD方法將信號分解成多個固有模函數(shù)(intrinsic mode function，IMF)分量后，將含有噪聲的高頻IMF分量舍棄，再將剩余的低頻IMF分量與殘差合成作為平滑后的原子鐘數(shù)據(jù)。實際上，每一個IMF分量中既有噪聲也有有用信號，簡單地舍棄一個或多個IMF分量，會導致相應的有用信號也同時被舍棄，造成信號的嚴重失真，影響降噪效果。

因此，針對EMD方法分解后的噪聲分量篩選及處理的問題，本文在傳統(tǒng)EMD方法的基礎上，結(jié)合小波閾值降噪法對其進行改進。針對原子鐘頻率數(shù)據(jù)非線性、非平穩(wěn)的特點，利用窗口做出劃分，對每一個窗口進行降噪處理。分別從時域和頻域?qū)υ摲椒ǖ慕翟胄Ч龀龇治?，并與傳統(tǒng)的EMD方法及常見的小波閾值降噪方法做出比較，以期提升原子鐘數(shù)據(jù)降噪效果。

1 原子鐘數(shù)據(jù)的模型

銫原子鐘的核心部件是原子頻標，在實際運行的過程中，原子鐘受到內(nèi)部噪聲及測量噪聲的影響，導致原子鐘的輸出信號存在相位偏差和頻率偏差[11-12]，降低了原子鐘的頻率穩(wěn)定度。相位數(shù)據(jù)和頻率數(shù)據(jù)都可以用來描述原子鐘的頻率穩(wěn)定度，且相位數(shù)據(jù)轉(zhuǎn)化為頻率數(shù)據(jù)，可以通過計算相位數(shù)據(jù)的一次差分除以取樣時間獲得，即：

(1)

其中：yi為頻率數(shù)據(jù)；xi和xi+1分別為i時刻和i+1時刻對應的相位數(shù)據(jù)；τ0為取樣時間間隔。

通常情況下，原子鐘的相位相對變化量x(t)和相對瞬時頻率偏差y(t)，都是由確定性變化部分和隨機變化部分組成，表達式分別見式(2)和式(3)：

(2)

y(t)=y0+Dt+εy(t)，

(3)

其中：x0為初始相位偏差；y0為初始頻率偏差；D為線性頻漂；εx(t)和εy(t)為隨機變化量，即隨機噪聲部分。關于原子鐘的隨機變化部分，目前國際上普遍認為原子鐘的隨機噪聲模型，可以看作是5種獨立的噪聲線性疊加而得到的符合冪律譜模型的經(jīng)驗公式[13-14]。

銫原子鐘數(shù)據(jù)具有非線性、非平穩(wěn)的特點。設含噪聲的銫原子鐘頻率數(shù)據(jù)為X(t)，可以表達為：

X(t)=f(t)+e(t)，

(4)

其中：X(t)為含有噪聲的原子鐘數(shù)據(jù)；f(t)為有用數(shù)據(jù)；e(t)為噪聲數(shù)據(jù)。為了提高結(jié)果的置信度，充分利用銫原子鐘數(shù)據(jù)，采取重疊ALLAN方差對原子鐘的頻率穩(wěn)定度進行評估。

2 降噪方法

EMD方法是在Hilbert-Huang變換(Hilbert-Huang transform，HHT)的基礎上，將信號分解至不同頻段，是一種循環(huán)迭代的算法[15-16]。EMD方法將原始信號分為若干個IMF和一個余項之和，不同頻率的IMF分量按頻率由高到低的順序排列，每一個IMF可以看作是一個函數(shù)且滿足以下條件：

(Ⅰ)在所有的數(shù)值序列中，最大點數(shù)和過零點數(shù)必須相等或最多相差1；

(Ⅱ)在任意點上，最大極值點和最小極值點的包絡的平均值為0。

EMD方法分解與原始信號都在同一尺度的時域中，保持了信號的非平穩(wěn)性，具有簡單、高效、自適應的特點。在頻域上，經(jīng)EMD分解出的IMF分量的頻率幾乎按2的負冪次方的形式遞減[17]，有用信號比較平穩(wěn)，對應IMF的低頻分量，而含有噪聲的信號波動較大，對應IMF的高頻分量。因此，可以通過頻譜分析找出高頻IMF分量，利用小波閾值對高頻IMF分量進行二次處理后再重構(gòu)，得到降噪信號。利用改進EMD方法對原子鐘頻率數(shù)據(jù)降噪，具體步驟如下：

(Ⅱ)找出子窗口內(nèi)數(shù)據(jù)Xi所有的局部極大值點和極小值點，利用3次樣條插值分別連接所有的極大值點和極小值點，作為上下包絡線u1(t)和v1(t)，計算上下包絡線的包絡均值m1(t)曲線，

(5)

(Ⅲ)計算窗口內(nèi)數(shù)據(jù)Xi和包絡均值m1(t)的差作為第一分量h1(t)，

h1(t)=X(t)-m1(t)。

(6)

(Ⅳ)若該分量沒有滿足IMF分量的兩個條件，則將該分量看作是窗口內(nèi)的數(shù)據(jù)，重復第Ⅱ步，得到h11(t),h12(t),…，h1k(t)，直到h1k(t)=h1(k-1)(t)-m1k(t)滿足IMF分量的兩個條件，稱h1k(t)為一階IMF分量，是該窗口內(nèi)數(shù)據(jù)Xi中最高頻率分量，記為：

c1(t)=h1k(t)。

(7)

(Ⅴ)將c1(t)從初始信號中分離，得到第1個剩余分量r1(t)，

r1(t)=X(t)-c1(t)。

(8)

(Ⅵ)由于r1(t)中仍存在較長周期的原始信號分量，因此重復上述步驟，得到：

(9)

(Ⅶ)直到rn(t)變成單調(diào)序列或不再滿足IMF調(diào)節(jié)時結(jié)束[18]，得到經(jīng)過EMD分解后的窗口內(nèi)的數(shù)據(jù)，可以表示為：

(10)

(Ⅷ)對IMF分量做出頻譜分析，將噪聲IMF分量篩選出來。

(Ⅸ)利用軟閾值函數(shù)對噪聲IMF分量進行小波閾值降噪，軟閾值函數(shù)為：

(11)

(Ⅹ)將降噪處理后的IMF分量與信號分量進行信號重構(gòu)，得到降噪后的窗口內(nèi)數(shù)據(jù)為：

(12)

(Ⅺ)在每個窗口內(nèi)重復以上步驟，得到降噪后的數(shù)據(jù)，

(13)

為了驗證改進的EMD方法的有效性，分別從信噪比(signal-to-noise ratio,SNR)、均方根誤差(root mean square error,RMSE)兩方面進行評價[19-20]。信噪比的數(shù)值越大，說明降噪效果越好。均方根誤差的數(shù)值越小，說明降噪效果越好。

3 實證分析

在實際計算中，選取一個月的銫原子鐘數(shù)據(jù)，采樣間隔τ0=1 h，在應用銫原子鐘頻率數(shù)據(jù)前，已對原始數(shù)據(jù)做出預處理。

采用EMD方法對銫原子鐘頻率數(shù)據(jù)進行分解，如圖1所示，將銫原子鐘頻率數(shù)據(jù)分解為8個IMF分量和1個殘差。對IMF分量做出頻譜分析，如圖2所示。由圖2可以看出：通過EMD分解后，各個IMF分量的頻率由高到低依次遞減，IMF1中包含了大部分噪聲。根據(jù)信號與IMF分量的對應關系，可以看出IMF1～IMF5為高頻分量即噪聲分量，其余為低頻分量即信號分量。

圖1 EMD分解后的IMF分量

在小波閾值降噪中，閾值選取和閾值函數(shù)的選擇十分重要，不同的選取方式對去噪效果有很大影響。常見的閾值選取方式有基于無偏估計原理的自適應閾值(rigrsure)、固定閾值(sqtwolog)、啟發(fā)式閾值(heursure)和極大極小閾值(minimax)。分別根據(jù)4種不同的閾值選取方式對高頻IMF分量做處理，由于處理結(jié)果存在的差異不是非常明顯，需要利用信噪比及均方根誤差兩個評價指標來判定降噪結(jié)果，見表1。根據(jù)降噪的評價指標可知：選用極大極小閾值進行處理的結(jié)果最好。

表1 不同閾值選取方式的降噪結(jié)果

表2 不同分組的降噪效果對比

根據(jù)表2，將銫原子鐘頻率數(shù)據(jù)分為5組分別做處理，去噪前后對比如圖3所示，藍色線條對應原始數(shù)據(jù)，紅色線條對應降噪數(shù)據(jù)。由圖3可以看出：改進EMD降噪后，噪聲得到了一定的抑制且保留了原始數(shù)據(jù)的主要圖像特征。從頻域上分析，分別做出降噪前后的頻譜，如圖4所示，其中，圖4a為降噪前的原始信號頻譜，圖4b為降噪后的信號頻譜。由圖4可以看出：經(jīng)過降噪處理后，波動明顯減小，說明噪聲能夠得到一定壓制。同時從時域上分析，為了提高結(jié)果的置信度，充分利用銫原子鐘數(shù)據(jù)，采取重疊ALLAN方差對原子鐘的頻率穩(wěn)定度進行評估。計算了降噪前后銫原子鐘頻率數(shù)據(jù)對應不同時間間隔的重疊ALLAN方差，如圖5所示，藍色線條代表數(shù)據(jù)降噪前對應不同時間間隔的重疊ALLAN方差，紅色線條為數(shù)據(jù)降噪后對應不同時間間隔的重疊ALLAN方差。由圖5可以看出：相同的時間間隔，通過改進EMD方法降噪的數(shù)據(jù)頻率穩(wěn)定度小于原始數(shù)據(jù)的頻率穩(wěn)定度。綜上所述，改進后的方法能夠有效提高其頻率穩(wěn)定度。

圖3 原始數(shù)據(jù)與降噪數(shù)據(jù)

(a) 原始信號頻譜

分別應用改進后的EMD方法、傳統(tǒng)的EMD方法及小波閾值方法，對同一組銫原子鐘數(shù)據(jù)做降噪處理，結(jié)果如圖6所示。圖6a為通過傳統(tǒng)EMD方法降噪的對比結(jié)果圖，原子鐘頻差數(shù)據(jù)經(jīng)過EMD分解后，將高頻的IMF分量舍棄，剩余的低頻IMF分量重構(gòu)作為降噪后的信號。由圖6a可以看出：簡單地去除一個或多個IMF分量的傳統(tǒng)EMD方法，會導致去除的IMF分量中的有用信息也被去掉，該方法較為粗糙，只能保持原始數(shù)據(jù)的大概趨勢。圖6b和圖6c分別是小波閾值算法及改進EMD方法的降噪結(jié)果的對比圖。由圖6b和圖6c可以看出：兩種方法都可以對噪聲產(chǎn)生一定的抑制作用，但是效果的差異并不是很明顯。利用信噪比及均方根誤差兩個評價指標來判定降噪結(jié)果，見表3。由表3可以發(fā)現(xiàn)：改進后的EMD方法降噪的效果優(yōu)于傳統(tǒng)的EMD方法及小波閾值降噪，說明該降噪方法可以有效降低原子鐘頻率數(shù)據(jù)中的噪聲。

圖6 不同方法降噪結(jié)果對比

表3 不同方法降噪的結(jié)果對比

4 結(jié)論

(1)傳統(tǒng)EMD方法分解后的噪聲數(shù)據(jù)篩選會造成部分有用信號的浪費，進而導致原數(shù)據(jù)的失真，在傳統(tǒng)EMD方法的基礎上，結(jié)合小波閾值對噪聲數(shù)據(jù)進行二次降噪，可以極大程度地保留噪聲數(shù)據(jù)中含有的少數(shù)有用數(shù)據(jù)，避免了有用數(shù)據(jù)的浪費及原始數(shù)據(jù)失真的情況。

(2)利用改進EMD方法對銫原子鐘頻差數(shù)據(jù)進行降噪，既保留了傳統(tǒng)EMD方法的優(yōu)點，又避免了Kalman濾波計算的復雜性和小波基函數(shù)的選擇，具有較好的自適應性，適用于具有非線性、非平穩(wěn)特點的銫原子鐘數(shù)據(jù)，減少了噪聲對時間尺度計算的影響，進一步得到精確的時間尺度，具有很強的實用價值。

(3)對原始數(shù)據(jù)進行加窗處理，效果比直接對數(shù)據(jù)進行處理的效果好，窗口數(shù)越多，得到的效果越好，但當數(shù)據(jù)分組超過一定數(shù)量時，由于每一組的數(shù)據(jù)量過少得到的效果反而很差。對噪聲數(shù)據(jù)進行二次降噪時，選用極大極小值的信噪比為3.453 9，均方根誤差為7.582 4e-14，均優(yōu)于其他3種閾值選擇方式，說明選用極大極小值的效果更好。

(4)改進EMD方法與傳統(tǒng)EMD方法和小波閾值方法相比較，信噪比從0.676 1和3.321 8提高到3.652 3，均方根誤差從1.044 0e-13和7.698 6e-14降低到7.411 1e-14，說明該方法可以抑制原子鐘頻差數(shù)據(jù)中的噪聲，進一步提高了原子鐘的頻率穩(wěn)定度。