原子鐘數據的改進經驗模態分解降噪

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

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

0 引言

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

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

因此，針對EMD方法分解后的噪聲分量篩選及處理的問題，本文在傳統EMD方法的基礎上，結合小波閾值降噪法對其進行改進。針對原子鐘頻率數據非線性、非平穩的特點，利用窗口做出劃分，對每一個窗口進行降噪處理。分別從時域和頻域對該方法的降噪效果做出分析，并與傳統的EMD方法及常見的小波閾值降噪方法做出比較，以期提升原子鐘數據降噪效果。

1 原子鐘數據的模型

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

(1)

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

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

(2)

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

(3)

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

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

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

(4)

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

2 降噪方法

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

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

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

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

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

(5)

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

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

(6)

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

c1(t)=h1k(t)。

(7)

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

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

(8)

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

(9)

(Ⅶ)直到rn(t)變成單調序列或不再滿足IMF調節時結束[18]，得到經過EMD分解后的窗口內的數據，可以表示為：

(10)

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

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

(11)

(Ⅹ)將降噪處理后的IMF分量與信號分量進行信號重構，得到降噪后的窗口內數據為：

(12)

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

(13)

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

3 實證分析

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

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

圖1 EMD分解后的IMF分量

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

表1 不同閾值選取方式的降噪結果

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

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

圖3 原始數據與降噪數據

(a) 原始信號頻譜

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

圖6 不同方法降噪結果對比

表3 不同方法降噪的結果對比

4 結論

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

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

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

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