新閾值函數(shù)去噪方法在可降水量預(yù)報(bào)中的應(yīng)用

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

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

1 引言

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

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

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

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

2 小波閾值去噪方法

2.1 基本原理

假設(shè)一個(gè)一維含噪信號(hào)為

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

2.2 小波閾值去噪的步驟

小波閾值去噪過(guò)程一般分為3個(gè)步驟［7］：

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

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

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

2.3 閾值函數(shù)

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

硬閾值函數(shù)為

圖1 小波分解樹(shù)圖

軟閾值函數(shù)為

其中dj，k為小波系數(shù)，sng（·）為符號(hào)函數(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)時(shí)會(huì)損失一些有用的高頻信息，造成重構(gòu)信號(hào)的信噪比較低，均方誤差較大；而硬閾值方法，雖然能夠保留原信號(hào)中的一些突變成分，但同時(shí)也可能產(chǎn)生新的不連續(xù)點(diǎn)，且對(duì)數(shù)據(jù)變化反應(yīng)也很敏感，在重構(gòu)信號(hào)時(shí)會(huì)出現(xiàn)一定的振蕩，導(dǎo)致大的方差［9］。

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

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

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

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

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

其中，Px為原始信號(hào)的能量，Px＾-x為噪聲的能量，N為原始信號(hào)的長(zhǎng)度。信噪比越高說(shuō)明信號(hào)中的噪聲越??；最小均方誤差越小說(shuō)明去噪后信號(hào)與原始純凈信號(hào)越接近，去噪效果越好。

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

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

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

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

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

4.1 基本原理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 小結(jié)

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

表3 泉州等8個(gè)地方的Kolmogorov熵

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

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

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

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

［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］孫斌，周云龍等.氣液兩相流壓差波動(dòng)信號(hào)小波去噪中閾值規(guī)則的確定［J］.信號(hào)處理，2006，22（1）：96-99

［7］汪同慶，郭子義，劉家兵，葉俊勇.一種基于改進(jìn)閾值函數(shù)的小波域超聲信號(hào)去噪方法［J］.無(wú)損檢測(cè)，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］.噪聲于振動(dòng)控制，2008，（1）：79-81

［10］朱育峰，嚴(yán)紹瑾，彭永清.從一維氣候時(shí)間序列直接提取K熵及其可預(yù)報(bào)時(shí)間的估計(jì)［J］.南京氣象學(xué)院學(xué)報(bào)，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］趙貴兵，石炎福，段文峰，等.從混沌序列同時(shí)計(jì)算關(guān)聯(lián)維和Kolmogorov熵［J］.計(jì)算物理，199：16（5）：309～315

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

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