陳晞，杜言霞，吳勇凱，溫繼昌

（泉州市氣象局，福建泉州，362000）

1 引言

大氣系統(tǒng)是一個非線性系統(tǒng)，大氣觀測數(shù)據(jù)必然是非線性時間序列。實(shí)際獲得的大氣觀測數(shù)據(jù)總是不可避免地存在噪聲，噪聲的存在會淹沒信號中的部分細(xì)節(jié)，給分析和研究帶來困難。因此，降噪是時間序列分析的一個重要內(nèi)容。由于大氣觀測數(shù)據(jù)是非線性時間序列，線性濾波方法不僅不能將噪聲去除，而且可能使真實(shí)信號發(fā)生扭曲，從而引入新的噪聲。

基于小波變換的非線性濾波卻能很好的對這樣的非線性序列進(jìn)行去噪。小波濾波方法與線性濾波方法相比，能夠在去除噪聲的同時，很好地保留信號的突變部分，而且還不用考慮信號的頻帶是否重疊［1，2］。由于小波去噪的這些優(yōu)點(diǎn)，小波去噪方法得到了越來越廣泛的應(yīng)用。

小波去噪的機(jī)理是基于信號與噪聲的小波系數(shù)在尺度上的不同性質(zhì)，采用相應(yīng)的規(guī)則，對含噪信號的小波系數(shù)進(jìn)行取舍、抽取或切削等非線性處理，以達(dá)到去除噪聲的目的。Donoho和Johnstone在1994年提出的小波閾值法（WaveShrink）是小波域去噪的主要方法之一。該方法具有非常明顯的漸進(jìn)近似最優(yōu)性質(zhì)，可以在均方差意義下取得最優(yōu)的去噪效果［3，4］。

小波閾值去噪算法中，選取合適的閾值函數(shù)是最基本的問題之一，對去噪效果有非常重要的影響［5］。本文首先討論經(jīng)典的硬閾值函數(shù)（hard shrinkage function）、軟閾值函數(shù)（soft shrinkage function）和本文提出的新閾值函數(shù)，并運(yùn)用信噪比（SNR）和最小均方誤差（MMSE）對這幾種閾值函數(shù)去噪效果進(jìn)行衡量。然后將新閾值函數(shù)應(yīng)用于大氣可降水量時間序列，運(yùn)用Kolmogorov熵和神經(jīng)網(wǎng)絡(luò)預(yù)測對去噪效果進(jìn)行衡量。

2 小波閾值去噪方法

2.1 基本原理

假設(shè)一個一維含噪信號為

其中，fi為真實(shí)信號，σ為常數(shù)，ei為高斯白噪聲，xi為含噪信號。一般情況下，有用信號通常表現(xiàn)為低頻信號或是一些比較平穩(wěn)的信號，而噪聲信號則表現(xiàn)為高頻信號。小波分析可以將信號多層分解，如圖1所示，噪聲部分通常包含在d1、d2、d3等高頻系數(shù)中，利用閾值函數(shù)對小波系數(shù)進(jìn)行取舍、抽取或切削等非線性處理，然后對處理后的小波系數(shù)進(jìn)行重構(gòu)，可以達(dá)到去噪的目的［6］。

2.2 小波閾值去噪的步驟

小波閾值去噪過程一般分為3個步驟［7］：

①對含噪信號進(jìn)行小波變換，得到各層小波系數(shù)dj，k；

②利用閾值函數(shù)對小波系數(shù)dj，k進(jìn)行處理，在盡可能多的保留小波系數(shù)中有用信號部分的同時，盡可能多的除去噪聲部分。

在小波去噪的三個步驟中,如何對小波系數(shù)進(jìn)行非線性處理是小波去噪的關(guān)鍵步驟，本文在經(jīng)典軟、硬閾值函數(shù)的基礎(chǔ)上提出了新閾值函數(shù)。

2.3 閾值函數(shù)

經(jīng)典軟、硬閾值函數(shù)定義如下：

硬閾值函數(shù)為

圖1 小波分解樹圖

軟閾值函數(shù)為

其中dj，k為小波系數(shù)，sng（·）為符號函數(shù)，d＾j，k為處理后的小波系數(shù)，tj為第j層的閾值，可根據(jù)如下經(jīng)驗(yàn)公式計(jì)算［8］

軟、硬閾值函數(shù)去噪雖然在實(shí)際中得到了廣泛的應(yīng)用，也取得了較好的效果，但它們本身存在著缺點(diǎn)。軟閾值方法中，系數(shù)的估計(jì)和原系數(shù)存在著恒定的差值，發(fā)生了系數(shù)收縮，因而重構(gòu)時會損失一些有用的高頻信息，造成重構(gòu)信號的信噪比較低，均方誤差較大；而硬閾值方法，雖然能夠保留原信號中的一些突變成分，但同時也可能產(chǎn)生新的不連續(xù)點(diǎn)，且對數(shù)據(jù)變化反應(yīng)也很敏感，在重構(gòu)信號時會出現(xiàn)一定的振蕩，導(dǎo)致大的方差［9］。

因此，本文在經(jīng)典軟、硬閾值函數(shù)的基礎(chǔ)上提出了新閾值函數(shù)，定義如下：

式中M為任意正常數(shù)，一般取為正整數(shù)，本文取M＝5，當(dāng)M趨于無窮時，新閾值函數(shù)逼近于硬閾值函數(shù)。

從定義可以看出，新域值函數(shù)克服可硬閾值函數(shù)不連續(xù)性的缺點(diǎn)和軟閾值函數(shù)對大于閾值的小波系數(shù)進(jìn)行收縮而產(chǎn)生較大偏差的缺點(diǎn)。

3 仿真實(shí)驗(yàn)

為了比較各種閾值函數(shù)去噪的效果，本文選用MATLAB中wnoise函數(shù)產(chǎn)生的含噪信號bumps來進(jìn)行測試。運(yùn)用信噪比（SNR）和最小均方誤差（MMSE）來檢驗(yàn)去噪效果的好壞，信噪比和最小均方誤差定義為：

其中，Px為原始信號的能量，Px＾-x為噪聲的能量，N為原始信號的長度。信噪比越高說明信號中的噪聲越?。蛔钚【秸`差越小說明去噪后信號與原始純凈信號越接近，去噪效果越好。

3.1 確定最優(yōu)分解層數(shù)

為了確定對原始信號分解為多少層時去噪效果最好，本文選用db4小波對bumps信號進(jìn)行分解，并用新閾值函數(shù)對小波系數(shù)進(jìn)行處理處理，表1列出了不同分解層數(shù)的去噪結(jié)果。從表1可以看出，當(dāng)分解層數(shù)N為5時，信噪比最大，最小均方誤差最小，這時的去噪效果最好，因此本文選用的分解層數(shù)為5。

3.2 各種閾值函數(shù)去噪效果的比較

為了比較上述三種閾值函數(shù)的去噪效果，本文用db4小波對bumps信號進(jìn)行5層分解，并分別用上述三種閾值函數(shù)進(jìn)行處理，重構(gòu)獲得去噪信號。三種閾值函數(shù)去噪效果如表2所示。從表2可以看出，軟閾值函數(shù)去噪的信噪比最小，最小均方誤差最大，去噪效果最差；新域值函數(shù)去噪的信噪比最大，最小均方誤差最小，去噪效果最好。

4 新閾值函數(shù)去噪在大氣可降水量時間序列中的應(yīng)用

4.1 基本原理

大氣系統(tǒng)是非線性系統(tǒng)，對大氣可降水量時間序列的去噪也需要采用非線性去噪方法。前文介紹的新閾值函數(shù)在bumps信號去噪中取得了較好的效果，這說明了新域值函數(shù)在非線性時間序列去噪中的優(yōu)越性。本節(jié)將把新閾值函數(shù)應(yīng)用于大氣可降水量時間序列，檢驗(yàn)新閾值函數(shù)是否適用于大氣可降水量時間序列。

在實(shí)際中干凈的信號是不存在的，不能像前文一樣運(yùn)用信噪比和最小均方誤差來衡量去噪效果的好壞。因此，需要重新尋找一種能夠衡量信號去噪效果的物理量。

文獻(xiàn)10、11指出Kolmogorov熵表征了系統(tǒng)在n維相空間中的總體膨脹性質(zhì)，從而給出系統(tǒng)演化過程中可能具有的不可預(yù)報性。利用Kolmogorov熵值可以直接分辨系統(tǒng)是規(guī)則或是不規(guī)則運(yùn)動。當(dāng)Kolmogorov熵為0時，運(yùn)動有序；當(dāng)Kolmogorov熵趨于無窮時，運(yùn)動為隨即狀態(tài)，當(dāng)Kolmogorov熵介于0與無窮大之間時，對應(yīng)的系統(tǒng)便為混沌運(yùn)動。Kolmogorov熵值的倒數(shù)描寫了在一定精度內(nèi)，使系統(tǒng)原信息仍然有“效”的平均時間尺度，即系統(tǒng)的平均可預(yù)報時間［10，11］。

很容易知道噪聲會使時間序列的可預(yù)報性降低，可預(yù)報時間縮短。如果時間序列中含有的噪聲較大，則可預(yù)報時間較短，Kolmogorov熵較大；如果時間序列中所含噪聲較小，則可預(yù)報時間較長，Kolmogorov熵較小。因此，本節(jié)選用Kolmogorov熵來衡量新閾值函數(shù)對大氣可降水量時間序列去噪效果的好壞。

表1 不同分解層數(shù)時新域值函數(shù)的去噪效果

表2 硬閾值函數(shù)、軟閾值函數(shù)和新域值函數(shù)的去噪效果

同時，本文還將運(yùn)用RBF神經(jīng)網(wǎng)絡(luò)對去噪前后的可降水量時間序列進(jìn)行一步預(yù)測，觀察去噪信號的可預(yù)測性是否增強(qiáng)。如果去噪時間序列的預(yù)測效果比原始時間序列的去噪效果好，則說明了新閾值函數(shù)去噪減小了原始時間序列中的噪聲，從而適用于可降水量時間序列的去噪。

Grassberger和Procaccia提出了從時間序列計(jì)算吸引子關(guān)聯(lián)維D2的G-P算法［12］，并提出通過計(jì)算K2熵來逼近Kolmogorov熵的思想［13］，但該方法比較非常復(fù)雜。趙貴兵，石炎福等在G-P算法的基礎(chǔ)上提出用最小二乘法同時從時間序列計(jì)算出關(guān)聯(lián)維和Kolmogorov熵的方法，該方法得到了Kolmogorov熵的穩(wěn)定估計(jì)［14］。該方法計(jì)算簡單，且結(jié)果較好，因此，本文直接使用該方法來計(jì)算可降水量時間序列的Kolmogorov熵。

4.2 Kolmogorov熵驗(yàn)證可預(yù)測性

采用泉州、龍巖等8個地方2012年4月30日00時到2013年4月30日18時每日4次的分辨率為1°×1°的大氣可降水量數(shù)據(jù)，每個地方總共1464個數(shù)據(jù)點(diǎn)，數(shù)據(jù)來源為NCEP再分析格點(diǎn)資料。首先，采用新閾值函數(shù)對8個地方的原始含噪信號進(jìn)行濾波，獲得去噪信號。然后，采用最小二乘法分別計(jì)算這8個地方的原始信號和去噪信號的Kolmogorov熵，結(jié)果如表3所示。從表中可以看出，去噪信號與原始信號相比，Kolmogorov熵明顯減小，說明去噪信號的平均可預(yù)報時間增長，可預(yù)報性增強(qiáng)，去噪信號比原始信號相比噪聲明顯減小。因此，新閾值函數(shù)同樣適用于大氣可降水量時間序列。

4.3 RBF神經(jīng)網(wǎng)絡(luò)對可降水量時間序列進(jìn)行預(yù)測

近年來，人工神經(jīng)網(wǎng)絡(luò)在氣象預(yù)報領(lǐng)域獲得了較快的發(fā)展，并取得了較好的結(jié)果，部分成果甚至高于業(yè)務(wù)預(yù)報水平［15］。本文運(yùn)用RBF神經(jīng)網(wǎng)絡(luò)對大氣可降水量時間序列進(jìn)行預(yù)測，根據(jù)預(yù)測結(jié)果的好壞來判斷去噪時間序列的可預(yù)測性是否增強(qiáng)，噪聲是否減小。馬莎、肖明等在基于混沌RBF網(wǎng)絡(luò)的地下廠房圍巖變形預(yù)報一文中詳細(xì)介紹了RBF神經(jīng)網(wǎng)絡(luò)對非線性時間序列進(jìn)行預(yù)報的算法算法［16］。

首先將泉州地區(qū)的大氣可降水量時間序列運(yùn)用新閾值函數(shù)進(jìn)行去噪，并將原始時間序列和去噪時間序列歸一化到均值為0，幅度為1。然后分別對這兩個歸一化序列進(jìn)行處理，前1300個數(shù)據(jù)點(diǎn)作為訓(xùn)練樣本，之后100個數(shù)據(jù)作為測試序列。為了衡量測試結(jié)果的好壞，我們引入了預(yù)測絕對誤差ERR和最小均方誤差MMSE，其定義如下：

測試結(jié)果如圖2（a）、圖2（b）、圖2（c）所示。從圖可以看出，去噪時間序列的預(yù)測誤差明顯小于原始時間序列的預(yù)測誤差。計(jì)算可得，原始時間序列預(yù)測值的最小均方誤差為0.00321，去噪時間序列預(yù)測值的最小均方誤差為0.00097，所以去噪時間序列的可預(yù)測性明顯增強(qiáng)了。

同樣，用RBF網(wǎng)絡(luò)對其它7個地方的可降水量時間序列及其去噪序列進(jìn)行預(yù)測，其最小均方誤差如表4所示。表4顯示，這8個地方大氣的可將水量時間序列，經(jīng)新閾值函數(shù)去噪后，時間序列的可預(yù)測性明顯增強(qiáng)，與Kolmogorov熵所得出的結(jié)果一致，因此新閾值函數(shù)適用于可降水量時間序列的去噪。

5 小結(jié)

小波閾值去噪是目前應(yīng)用最為廣泛的一種小波去噪方法，它在去除噪聲的同時能較好的保留信號的局部特性。本文在經(jīng)典軟、硬閾值去噪的基礎(chǔ)上提出了新閾值函數(shù)去噪。測試結(jié)果表明，與經(jīng)典軟、閾值函數(shù)去噪相比，新閾值函數(shù)去噪提高了信號的信噪比，同時降低了最小均方誤差。將新閾值函數(shù)去噪算法應(yīng)用于大氣可降量時間序列，并運(yùn)用Kolmogorov熵和RBF神經(jīng)網(wǎng)絡(luò)預(yù)測對去噪效果進(jìn)行衡量，結(jié)果表明，新閾值函數(shù)去噪提高了時間序列的可預(yù)測性。因此，新閾值函數(shù)更適用于可降水量時間序列等大氣觀測數(shù)據(jù)的去噪過程。

表3 泉州等8個地方的Kolmogorov熵

圖2 （a）原時間序列中測試序列及其預(yù)測結(jié)果

圖2 （b）去噪時間序列中測試序列及其預(yù)測結(jié)果

圖2 （c）原時間序列和去噪時間序列的預(yù)測絕對誤差

表4 泉州等地區(qū)原始序列和去噪序列預(yù)測的最小均方差

［1］潘泉，張磊等.小波濾波方法及應(yīng)用［M］.北京：清華大學(xué)出版社，2005：5-9.

［2］Vidakovic B，lozoya C B.On time-dependent wavelet denoising［J］.IEEE Trans.Signal Processing，1998，46（9）：2549～2551

［3］Donoho D L.De-noising by soft—thresholding［J］.IEEE Trans Inform Theory，1995，41（3）：613—627

［4］Donoho D L，Johnstone I M.Ideal spatial adaptation by wavelet shrinkage［J］.Biometrika，1994，81（3）：425-455

［5］Bruce A G，Gao Hongye.Understanding wave shrink：cariance and bias estimation［J］.Biometrika，1996，83（4）：727—745

［6］孫斌，周云龍等.氣液兩相流壓差波動信號小波去噪中閾值規(guī)則的確定［J］.信號處理，2006，22（1）：96-99

［7］汪同慶，郭子義，劉家兵，葉俊勇.一種基于改進(jìn)閾值函數(shù)的小波域超聲信號去噪方法［J］.無損檢測，2006,29（11）：641-643

［8］Ibarra-Castanedo C，Galmiche F,Darabi A，et al.Thermographic nondestructive evaluation：overview of recent progress［J］.Proceedings of SPIE，2003

［9］林穎，常永貴等.基于一種新閾值函數(shù)的小波閾值去噪研究［J］.噪聲于振動控制，2008，（1）：79-81

［10］朱育峰，嚴(yán)紹瑾，彭永清.從一維氣候時間序列直接提取K熵及其可預(yù)報時間的估計(jì)［J］.南京氣象學(xué)院學(xué)報，1994，（1）.

［11］林振山.非線性力學(xué)與大氣科學(xué)［M］.南京：南京大學(xué)出版社，1993：305-310

［12］Grassberger P，Procaccia I.Characterization of strange attractors.Phys Rev Lett，1983，50（5）：346-349.

［13］Grassberger P，Procaccia I.Estimation of the Kolmogorov entropy from a chaotic signal.Phys Rev A，1983，28（4）：2591-2593.

［14］趙貴兵，石炎福，段文峰，等.從混沌序列同時計(jì)算關(guān)聯(lián)維和Kolmogorov熵［J］.計(jì)算物理，199：16（5）：309～315

［15］洪展.一次臺風(fēng)暴雨過程的水汽特征分析.氣象研究與應(yīng)用，2014，35（4）.

［16］馬莎、肖明等.基于混沌RBF網(wǎng)絡(luò)的地下廠房圍巖變形預(yù)報.武漢大學(xué)學(xué)報，2009.