牟士杭，楊 暉，范彥平，李 然，韓 韌，華云松

（上海理工大學(xué) 光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093）

引 言

散斑能見度光譜(speckle visibility spectroscopy,SVS)法是建立在擴(kuò)散波光譜法(DWS)理論的基礎(chǔ)上發(fā)展出來的一種用于研究時(shí)變動(dòng)力學(xué)(timevarying dynamics)的新方法。該方法有效解決了DWS法不適用于非平穩(wěn)隨機(jī)過程的問題。在測量的被動(dòng)層中會(huì)產(chǎn)生許多隨機(jī)噪聲，被動(dòng)層顆粒溫度運(yùn)動(dòng)產(chǎn)生的信號(hào)較小，加上產(chǎn)生隨機(jī)噪聲，產(chǎn)生信噪比較小的信號(hào)，因此要想分析原始信號(hào)產(chǎn)生的顆粒溫度，提高信噪比并最大限度地提高濾除噪聲信號(hào)所用的方法非常重要。

近十年來，小波變換在國際上發(fā)展迅速并成為熱門的研究領(lǐng)域。小波變換能同時(shí)進(jìn)行時(shí)間和空間頻率的局部化的分析，在時(shí)域和頻域方面有良好的分辨性，基于多尺度分析理論，在不同尺度下的信號(hào)具有不同的局部特征。小波基函數(shù)對(duì)被分析信號(hào)的低頻段有著很好的頻率分辨率，對(duì)高頻段有著高時(shí)間分辨率；在小波多分辨率理論基礎(chǔ)上，把信號(hào)攜帶的信息分解到細(xì)節(jié)部分，然后得到的細(xì)節(jié)部分再對(duì)其進(jìn)行分析；連續(xù)信號(hào)的細(xì)節(jié)信息和噪聲信號(hào)的細(xì)節(jié)信息具有不同的特征[1-4]。

本文利用小波分析的特有屬性建立一種低信噪比被動(dòng)層顆粒溫度信號(hào)分析的方法，實(shí)現(xiàn)顆粒溫度信號(hào)中噪聲的抑制，提高顆粒溫度檢測精度，以便更好地分析顆粒運(yùn)動(dòng)的特征和規(guī)律。

1 散斑能見度光譜法原理

式中：m為粒子的總數(shù)；vi為第i個(gè)粒子的速度，為所有粒子運(yùn)動(dòng)的平均速度。

激光散斑法的原理如圖1所示。激光器發(fā)出的光通過凹透鏡擴(kuò)散到滾筒中被動(dòng)層粒子上，然后經(jīng)粒子間多次散射后，在空間產(chǎn)生干涉現(xiàn)象，發(fā)生如圖1(b)所示的“散斑”[7]現(xiàn)象。當(dāng)滾筒轉(zhuǎn)動(dòng)時(shí)，粒子也發(fā)生移動(dòng)，散斑圖像產(chǎn)生隨機(jī)波動(dòng)，被稱為“動(dòng)態(tài)散斑”。計(jì)算散斑圖像的對(duì)比度就可給出測量區(qū)域中顆粒位移變化，就能測出顆粒速度波動(dòng)，即顆粒脈動(dòng)速度[7]。對(duì)于不規(guī)則運(yùn)動(dòng)的粒子，散斑的波動(dòng)率和粒子的脈動(dòng)速度相關(guān)，由于線陣CCD照相機(jī)上的散斑圖像波動(dòng)隨時(shí)間變化，就可得到相對(duì)應(yīng)的顆粒速度波動(dòng)隨著時(shí)間變化而相應(yīng)地改變，如圖1(c)所示[7-8]。

2 顆粒溫度信號(hào)小波去噪原理

2.1 小波函數(shù)的基本特點(diǎn)

函數(shù)x(t)的連續(xù)小波變換在所選小波基函數(shù)的條件下，其擴(kuò)展信號(hào)為

式中：a為尺度因子；b為平移因子；

為小波基函數(shù)。

由式(2)可知，a和b這兩個(gè)參數(shù)對(duì)小波變換有很大的影響。當(dāng)a很小時(shí)，小波更集中于時(shí)域，并獲得信號(hào)細(xì)節(jié)；當(dāng)a很大時(shí)，小波形狀被擴(kuò)展并且可以獲得信號(hào)的宏觀信息，就有效濾除噪聲和得到高不失真信號(hào)[9]。

小波變換的反演公式為

在實(shí)際應(yīng)用的過程中，由于處理的信號(hào)絕大部分是經(jīng)測量信號(hào)采集系統(tǒng)得到的離散的數(shù)據(jù)信號(hào)，因此本文將采取小波變換離散形式。

圖1 激光散斑法的原理示意圖Fig. 1 Schematic diagram of laser speckle method

2.2 小波基函數(shù)的選取

目前采用信號(hào)檢測來確定小波函數(shù)的常用方法是最大限度地讓小波的波形接近待測瞬態(tài)的信號(hào)波形[10]。通過實(shí)驗(yàn)可以得到，基于小波變換的整個(gè)時(shí)頻域的分辨率分析，DB(Daubechies)小波基函數(shù)是一個(gè)較好的處理低信噪比顆粒溫度信號(hào)的方法[10-11]。如圖2所示左邊為DB12的尺度函數(shù)與小波函數(shù)，右邊是粒子溫度信號(hào)圖，兩個(gè)圖之間具有很強(qiáng)的相似性，都是由小的正弦波組成，具有正交性[10]。

圖2 DB12 小波函數(shù)與顆粒溫度信號(hào)對(duì)比Fig. 2 Comparison between DB12 wavelet function and particle temperature signal

實(shí)際過程中通常需要結(jié)合具體情況來具體分析，采用正演與反演的實(shí)驗(yàn)方法選擇合適的小波函數(shù)。本文采用DB N(N為分解層數(shù))作為小波基函數(shù)。DB系列小波函數(shù)是正交和雙正交的，長度為2N?1，濾波長度為2N，消失矩為N。經(jīng)過數(shù)據(jù)處理和實(shí)驗(yàn)對(duì)比，我們可以知道，小波分解層數(shù)并不是越多越好。當(dāng)分解層數(shù)超過某數(shù)時(shí)，許多真實(shí)的信號(hào)會(huì)損失；在分解層數(shù)達(dá)不到某值時(shí)，不能有效地濾除噪聲。在實(shí)際應(yīng)用中，分解層的數(shù)量應(yīng)該通過比較峰值信噪比等指標(biāo)來確定。

3 被動(dòng)層顆粒溫度信號(hào)去噪算法的設(shè)計(jì)

令x0(t)表示在理想條件下被動(dòng)層粒子溫度的瞬態(tài)信號(hào)，在有噪聲n(t)的條件下，粒子溫度信號(hào)為

在低信噪比的條件下，對(duì)粒子溫度信號(hào)測量，需從混合信號(hào)x(t)中踢除噪聲信號(hào)，使原始信號(hào)被最大限度地提取。

根據(jù)采樣定理，將x(t)、n(t)和x0(t)分別進(jìn)行 N點(diǎn)離散化，可得到一個(gè)N點(diǎn)序列x(k)、n(k)和x0(k) (k =1,2,···,N)。這樣式(4)將變?yōu)橐粋€(gè)N點(diǎn)序列

根據(jù)小波變換分解原理，采用為Mallat算法分解圖。Mallat算法是在基于小波信號(hào)分解塔式算法的條件下得到的小波快速算法，其算法在小波變換里的作用與快速傅里葉變換在傅里葉變換中起到的作用是相同的[12-13]。

根據(jù)顆粒溫度信號(hào)與噪聲在小波變換下的不同特性，運(yùn)用的Mallat算法流程圖如圖3所示[14]，圖中G(w)表示低頻信號(hào)，H(w)表示高頻信號(hào)。

圖3 信號(hào) x(t)的 Mallat分解結(jié)構(gòu)Fig. 3 Mallat decomposition structure of signal x(t)

基于小波變換分解原理，采用Mallat小波算法進(jìn)行分析，基于信號(hào)和噪聲在小波變換下的不同特征，算法的設(shè)計(jì)過程如下[15]。

步驟1采用一維小波分析對(duì)顆粒溫度數(shù)據(jù)信號(hào)進(jìn)行分解。顆粒溫度有噪的信號(hào)表示為x(t)=x0(t)+n(t)，x(t)為含噪的顆粒溫度信號(hào)，x0(t)為真實(shí)的被動(dòng)層顆粒溫度信號(hào)，n(t)為噪聲信號(hào)。n(t)為高斯白噪聲，正常的情況下是高頻信號(hào)，x0(t)是低頻顆粒溫度信號(hào)，然后確定一個(gè)合適的小波函數(shù)并確定分解層數(shù)，最后一維分解計(jì)算。

步驟2通過小波分析對(duì)實(shí)測到的粒子溫度信號(hào)進(jìn)行分解，得到低頻信號(hào)與許多高頻信號(hào)。

步驟3小波分解后對(duì)所述顆粒溫度數(shù)據(jù)圖的高頻信號(hào)不變，在這個(gè)基礎(chǔ)上對(duì)分解好的低頻信號(hào)進(jìn)行又一次的分解處理，重新獲得高頻與低頻，然后做如下處理[16]：

對(duì)所述的顆粒溫度信號(hào)圖小波分解之后的高頻部分保留恒定；再線性變換較低分辨率顆粒溫度信號(hào)，在低頻信號(hào)的小波系數(shù)和數(shù)值關(guān)系自適應(yīng)幅度增強(qiáng)條件下，重新得到低頻信號(hào)。

步驟4小波分解高頻系數(shù)的閾值量化，高頻系數(shù)為不同尺度下軟閾值所決定的[17]。軟閾值去噪法：dm(x)=式中：x為實(shí)際測量顆粒溫度系數(shù)；dm(x)為閾值化分析計(jì)算出一維小波系數(shù)；sgn為一符號(hào)函數(shù)；T為閾值，表示為一個(gè)軟閾值過程。在通常的情況下噪聲的特性是按照高斯分布N(0, 1)，在這里閾值T的計(jì)算依據(jù)為最小均方差的極值方法的恒定閾值，對(duì)于噪聲未知和非白噪聲該方法很適合[10]。此法把不知道的有噪聲信號(hào)當(dāng)作和未知回歸函數(shù)的估計(jì)式類似，并且可以使給定函數(shù)內(nèi)最大均方誤差最小化，然后對(duì)染噪顆粒溫度信號(hào)降噪處理，并且在處理期間對(duì)每一層選擇合適的閾值降噪，最后重構(gòu)信號(hào)[17]。

步驟5進(jìn)行小波重構(gòu)，小波分解后細(xì)節(jié)部分低頻系數(shù)與高頻系數(shù)再經(jīng)過小波重構(gòu)得出處理后的顆粒溫度信號(hào)[18]。

步驟6根據(jù)Mallat重構(gòu)算法恢復(fù)信號(hào)并對(duì)信號(hào)小波逆變換。在噪聲抑制后凈化顆粒溫度信號(hào)，對(duì)新的低頻部分和顆粒溫度信號(hào)圖小波分解之后的高頻部分進(jìn)行小波逆變換以獲得處理后的顆粒溫度信號(hào)圖。

4 計(jì)算機(jī)仿真與結(jié)果評(píng)估

用計(jì)算機(jī)仿真時(shí)，選擇原始信號(hào)x0(t)為矩形波，使用不同的小波函數(shù)進(jìn)行模擬比較。選擇具有良好正交性與雙正交性效應(yīng)的DB5的小波來分析模擬信號(hào)，含噪信號(hào)為x(t)=x0(t)+n(t)，信號(hào)仿真點(diǎn)數(shù)為N=2000，仿真結(jié)果如圖4所示。圖4(a)是一個(gè)無噪聲的矩形波信號(hào)，(b)是在(a)的基礎(chǔ)上疊加有噪聲的矩形波信號(hào)，對(duì)(a)和(b)中信號(hào)進(jìn)行傅里葉變換以獲得原始方波信號(hào)的頻譜和含噪聲信號(hào)的頻譜，然后作小波變換。采用DB5小波，作3層分解，用stein自適應(yīng)軟閾值對(duì)系數(shù)作閾化處理，得到的消噪效果如圖4(c)所示。另外，在作傅里葉變換時(shí)采取低通濾波器去噪，圖4(d)是傅里葉變換去噪的結(jié)果圖，我們可從圖中發(fā)現(xiàn)，低通濾波器在除去噪聲信號(hào)過程中把一些真實(shí)的信號(hào)給丟失了，濾波器的寬度對(duì)去噪效果影響很大。比較圖4(c)和(d)可以看出，對(duì)于非平穩(wěn)信號(hào)去噪，小波變換的消噪效果優(yōu)于傅里葉變換效果[19]。

圖4 算法在加噪的方波信號(hào)中的特性分析Fig. 4 Analysis of the characteristics of the algorithm in the noise-induced square wave signal

5 數(shù)據(jù)測量與處理

5.1 測量系統(tǒng)的搭建

實(shí)驗(yàn)裝置如圖5所示，滾筒內(nèi)直徑是150 mm其運(yùn)動(dòng)是由4個(gè)底座固定在滾軸上的驅(qū)動(dòng)輪帶動(dòng)的，用的是德國Dunker公司生產(chǎn)的直流電動(dòng)機(jī)和減速器，通過閉環(huán)控制電機(jī)帶動(dòng)滾筒以一定的速度轉(zhuǎn)動(dòng)[7]。在實(shí)驗(yàn)過程中用的激光發(fā)射器Nova ProDPSS(波長為532 μm，功率為300 mW)是由RGB公司設(shè)計(jì)的，入射光被一個(gè)凹透鏡擴(kuò)展，被一個(gè)平面鏡反射，并最終入射到圓柱體中被動(dòng)層中的粒子上，在圖中z方向向下的位置，如圖5所示[7]。采用的CCD照相機(jī)是由DALSA公司設(shè)計(jì)的，其像素為1024，單個(gè)像素的長度是14 nm，最大線速度為68000 Hz，CCD表面配有波長為0.532 nm的濾光片，其作用是去掉周圍所帶來的無用光[7]。在實(shí)驗(yàn)中采用的顆粒是一種規(guī)則粒徑的0.15 mm玻璃砂[20]。

圖5 滾筒和 SVS 測量系統(tǒng)示意圖Fig. 5 Roller and SVS measurement system structure diagram

滾筒被動(dòng)層顆粒溫度在測量時(shí)，會(huì)伴隨很多其的噪聲，從而產(chǎn)生很多無用的噪聲加在被動(dòng)層顆溫度數(shù)據(jù)信號(hào)上，尤其是在高噪聲環(huán)境中，有用顆粒溫度信號(hào)被噪聲信號(hào)嚴(yán)重污染[15]，結(jié)果使粒子溫度測量信號(hào)的信噪比非常低，如果利用一般信號(hào)處理方法對(duì)被動(dòng)層顆粒溫度信號(hào)進(jìn)行檢測是非常困難的，因此要引入小波分析來提高系統(tǒng)測精度，以利顆粒溫度信號(hào)更好的分析與觀察[21]。

5.2 實(shí)驗(yàn)數(shù)據(jù)處理與分析

將小波算法應(yīng)用于被動(dòng)層粒子信號(hào)處理，選取電機(jī)電壓為0.9 V，圖7為經(jīng)過測量處理得出的含噪聲的顆粒溫度信號(hào)圖，即是運(yùn)用小波變換分解的顆粒溫度信號(hào)圖。選取DB12小波函數(shù)族將干擾后的顆粒溫度信號(hào)進(jìn)行多層分解，然后再重構(gòu)顆粒溫度信號(hào)，經(jīng)過12層小波分解后，信號(hào)波形如圖 7所示，圖中的 d1，d2，d3，d4，d5，d6，d7是干擾引起的噪聲成分，是高頻信號(hào)。

對(duì)圖6的顆粒溫度信號(hào)進(jìn)行小波變換，再使用軟閾值處理法對(duì)小波變換后的系數(shù)閾值處理，處理后零元素的個(gè)數(shù)占總元素個(gè)數(shù)之比為88.04%，去除噪聲并重構(gòu)顆粒溫度信號(hào)，可得到濾波后的顆粒溫度信號(hào)。與測量得出的實(shí)際顆粒溫度信號(hào)相比，降噪重構(gòu)后的顆粒溫度信號(hào)能量保留成分是實(shí)際測量顆粒溫度信號(hào)的98.4%，而且降噪后得到的顆粒溫度信號(hào)具有較好的規(guī)律性與光滑度，顆粒溫度信號(hào)分析更方便。

圖6 測量得出的被動(dòng)層顆粒溫度信號(hào)圖Fig. 6 The temperature signal of the passive layer is measured

從圖7中d9與d10可以看出，信號(hào)有局部突變的變化，此時(shí)可以認(rèn)為滾筒中的顆粒物質(zhì)運(yùn)動(dòng)產(chǎn)生大的轉(zhuǎn)變，即產(chǎn)生崩塌，據(jù)此我們就可以準(zhǔn)確定位或追蹤崩塌的位置。顆粒溫度信號(hào)的細(xì)節(jié)部分都集中在a12上，從圖形中可以更好地分析被動(dòng)層中的顆粒溫度。如圖8所示，降噪后的顆粒溫度信號(hào)中所包含的顆粒溫度瞬時(shí)峰值明顯降低，由此可以看出有用的信號(hào)集中于信號(hào)的底部，說明了小波變換可以有效地去除顆粒溫度被動(dòng)層的噪聲，使得顆粒溫度信號(hào)靈敏度得到較大的提高，具有12階消失矩的DB12小波能夠抑制多項(xiàng)式信號(hào)成分并且保留大部分測量得到顆粒溫度信號(hào)有用成分[22]。

在測量實(shí)驗(yàn)中可得出，小波基的不同選取對(duì)粒子溫度信號(hào)的去噪有很大的影響。選擇不合適時(shí)，即使是簡單的含噪聲信號(hào)都不會(huì)把噪聲信號(hào)很好且有效地濾除掉。除此以外，小波分解層數(shù)的數(shù)量和小波系數(shù)的閾值化方法也對(duì)去噪效果有很大影響?？傊?，利用小波變換對(duì)被動(dòng)層顆粒溫度信號(hào)去噪處理，可以極大地提高信號(hào)的信噪比和可檢測性，并且在被動(dòng)層顆粒溫度的測量信號(hào)處理方面有著非常大的實(shí)際應(yīng)用價(jià)值。

6 結(jié)束語

準(zhǔn)確測量微弱的顆粒運(yùn)動(dòng)信號(hào)對(duì)研究滾筒中顆粒運(yùn)動(dòng)具有非常重要的意義。利用小波變換對(duì)顆粒溫度信號(hào)多尺度分解，然后進(jìn)行細(xì)化處理，最后對(duì)信號(hào)小波重構(gòu)，得到的結(jié)果可最大限度地消除電機(jī)所產(chǎn)生的振動(dòng)和白光所造成的隨機(jī)噪聲。去噪后的效果圖能使我們更直觀地分析信號(hào)，能從細(xì)節(jié)信號(hào)中分辨出被動(dòng)層崩塌信號(hào)的間斷點(diǎn)，便于分析與觀察顆粒的運(yùn)動(dòng)特性，并提高顆粒溫度測量系統(tǒng)準(zhǔn)確性和精確度。

圖7 小波分解的層數(shù)示意圖Fig. 7 The number of layers of wavelet decomposition