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

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

0 引言

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

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

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

(4)

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

2 降噪方法

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

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

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

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

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

(5)

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

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

(6)

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

c1(t)=h1k(t)。

(7)

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

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

(8)

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

(9)

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

(10)

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

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

(11)

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

(12)

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

(13)

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

3 實(shí)證分析

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

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

圖1 EMD分解后的IMF分量

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

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

表2 不同分組的降噪效果對(duì)比

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

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

(a) 原始信號(hào)頻譜

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

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

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

4 結(jié)論

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

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

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

(4)改進(jìn)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ù)中的噪聲，進(jìn)一步提高了原子鐘的頻率穩(wěn)定度。