基于多尺度小波變換的紅外光譜譜峰識別算法

摘 要 傳統(tǒng)的譜峰檢測方法一般分為3個步驟：譜線平滑、基線校正和譜峰識別?，F(xiàn)有的基于小波變換的峰值檢測方法能較好地將基線校正和譜峰識別兩個步驟融為一步。在此基礎(chǔ)之上，本研究將譜線平滑也很好地融入到小波變換的峰值檢測算法中，使整個峰值檢測算法成為一個整體。在峰值提取時，原始譜圖直接處理，不再是處理加工過的譜圖，減小了譜峰檢測結(jié)果出錯的可能性。另外，對基于小波變換的譜峰檢測算法中模極大值檢測算法存在不確定閾值的問題進行了修改，使得基于小波變換的譜峰檢測算法更為完善。

關(guān)鍵詞 小波變換; 紅外光譜; 譜峰檢測

2010-08-19收稿；2011-01-21接受

本文系國家自然科學(xué)基金項目（No.50677047）資助

* E-mail:cai-tao@hotmail.com

1 引 言

傳統(tǒng)的譜峰檢測算法一般分為3個步驟：譜線平滑、基線校正與譜峰識別。圍繞這3個步驟，已提出大量算法。如采用移動平均法進行譜線平滑［1］；采用Kaiser 濾波器進行譜線平滑［2，3］。這兩種方法均是通過將原始信號與平滑函數(shù)直接相乘達到平滑譜線的目的，方法簡單，但信號失真較大。采用小波變換法平滑譜線［4，5］，由于小波變換系數(shù)存在較大的冗余度，所以在實現(xiàn)譜線平滑的同時，能較好地減少信息失真。在基線校準(zhǔn)方面，文獻［6］求出原始譜圖信號的一階導(dǎo)數(shù)，利用信號的單調(diào)性實現(xiàn)校正。文獻［3］采用線性插值法實現(xiàn)校正。文獻［7］則采用小波變換實現(xiàn)基線校正，方法簡單有效。在譜峰檢測方面，基于信噪比［2，3，6，7］、局部極大值［3，7，8］、斜率［9］、峰面積［4］、峰寬［10］、光譜變量降維［11］和譜峰模型［12］的方法都有報道。但這些算法的成效均在很大程度上依賴于前期的譜線平滑與基線校準(zhǔn)。Baggerly等［13］通過實驗證實了不同的譜線平滑與基線校準(zhǔn)算法在很大程度上決定傳統(tǒng)譜峰檢測法的結(jié)果。文獻［7，14］提出了基于小波變換的峰值檢測方法，并在該算法中很好地融合了基于小波變換的基線校準(zhǔn)。但該算法在光譜平滑方面融合的還不夠，算法中模極大值的搜索算法存在重復(fù)搜索的問題。文獻［15］對比了現(xiàn)有諸多譜峰檢測法，得出小波變換法的綜合性能是最佳的。本研究參考文獻［7，14］，并加以改進。將基于小波變換的譜線平滑也融合到譜峰檢測過程中，并在兩個不同階段分別實現(xiàn)一次譜線平滑，使譜峰檢測的精確度提高。改進了模極大值搜索算法，消除了原算法需要設(shè)定搜索閥值而帶來不確定因素的隱患，并且減小了原算法的復(fù)雜度。通過對SF6氣體紅外光譜譜圖的分析處理，證實了算法的優(yōu)越性

2 譜圖的小波變換原理和特征

2.1 基線校準(zhǔn)

在進行光譜分析前，有必要對光譜信號進行基線校準(zhǔn)。傳統(tǒng)的峰值提取方法是將基線校準(zhǔn)作為單獨步驟。采用小波變換提取峰值則避免了此步驟。假設(shè)F（t）為原始光譜譜圖，B（t）為基線函數(shù)，則譜圖中的有用信號函數(shù)F′（t）可表示為：

F′（t）＝F（t）－B（t）（1）

如果對原始譜圖進行小波變換，即：

Wf（a，b）＝∫RF（t）#8226;Ψ（a，b）dt＋∫RB（t）#8226;Ψ（a，b）dt（2）

一般情況下，基線是單調(diào)的，而且變換緩慢，所以基線函數(shù)B（t）可以表示為：

B（t）＝B′（t）＋C（3）

其中，B′（t）為過零點的斜坡函數(shù)，C為常數(shù)。而小波母函數(shù)在支撐域的積分為0。由此可知， 圖1 小波變換對基線的移除：f（t）為不含基線的原始信號，Wf為尺度為6時，其對應(yīng)的小波系數(shù)；g（t）為含有基線的原始信號；Wg為尺度為6時，其對應(yīng)的小波系數(shù)

Fig.1 Baseline remove: f（t） is original signal， Wf is wavelet coefficients of f（t） when the scale is 6， g（t） is the original contained baseline signal， Wg is wavelet coefficients of g（t） when scale is 6式（2）中后一部分結(jié)果為0。因此，含基線的原始光譜譜圖在進行小波變換后，其小波系數(shù)中已不包含基線信息。所以采用小波變換法進行譜峰檢測時，無需進行基線校準(zhǔn)。圖1為小波變換對基線移除的效果圖。

第6期蔡 濤等： 基于多尺度小波變換的紅外光譜譜峰識別算法

2.2 信號濾波

李氏指數(shù)是數(shù)學(xué)上用來表征函數(shù)局部特征的一種度量，它反映了曲線的光滑程度。即函數(shù)在某一點的李氏指數(shù)表征了該點的奇異性大小，指數(shù)越大，該點的光滑度越高；指數(shù)越小，該點光滑度越低。如果將該理論引入到小波變換中［16］，就會發(fā)現(xiàn)噪聲對應(yīng)的李氏指數(shù)一般為負(fù)數(shù)，而有用信號的李氏指數(shù)則為正數(shù)。對應(yīng)到小波系數(shù)中就是噪聲產(chǎn)生的模極大值會隨著尺度的變大而減小，有用信號產(chǎn)生的模極大值會隨著尺度的變大而變大。利用該規(guī)律，即可較好地識別真實信號和具有較高似真度的噪聲信號。

2.3 譜峰對應(yīng)的小波變換特征

Marr小波正比于高斯函數(shù)的二階導(dǎo)數(shù)。令θ（t）為高斯函數(shù)。Ψ（t）為Marr小波母函數(shù)。

Ψ（t）∝d2θ（t）dt2（4）

從另一個角度解釋小波變換，即將原信號f（t）與伸縮小波卷積得到（此處暫不考慮時間平移）。以Ψ（t）為小波母函數(shù)，則對信號f（t）在尺度a的小波變換可寫成：

Wf（t）＝f（t）#8226;Ψ（1a）∝f#8226;a2d2θadt2（t）＝a2d2（f#8226;θa）（t）dt2（5）

圖2 采用Marr小波對信號進行小波變換的系數(shù)矩陣。信號的峰位處，各尺度的小波系數(shù)均有模極值存在圖中的虛線部分

Fig.2 Coefficients matrix of original signal with Marr wavelet. There is a ridge along wave peak which remarked with dotted line

由于f#8226;θa（t）可看作是高斯函數(shù)在尺度a下對信號f（t）進行平滑的結(jié)果。而小波變換則wf（t）可看作信號f（t）在尺度a下由高斯函數(shù)平滑后再取二階導(dǎo)數(shù)的結(jié)果。高斯函數(shù)平滑能夠去除部分噪聲，而信號f（t）的二階導(dǎo)數(shù)則能凸顯出信號的位置與形態(tài)的顯著變化。圖2為采用Marr小波作為母函數(shù)對信號進行小波運算的結(jié)果。在信號發(fā)生突變處，各尺度的小波系數(shù)也有相應(yīng)的模極值。而且相鄰尺度的模極值連起來如同一條山脊，并收斂于對應(yīng)的譜峰位置。此處提到的模極值包括一個局部模極大值和緊鄰區(qū)域內(nèi)左右各一個局部模極小值。利用該性質(zhì)，本研究提出了新的譜峰檢測法。文獻［14］做了類似工作，但所用的小波母函數(shù)為Ridger小波，Ridger小波是奇函數(shù)，無法通過小波變換實現(xiàn)基線校準(zhǔn)。

3 多尺度小波變換的計算過程

據(jù)小波變換后各尺度小波系數(shù)之間的相互關(guān)系檢測譜峰。具體步驟如下：

對光譜信號進行小波變換，小波母函數(shù)選擇Marr小波，分解尺度為32。但每隔一個尺度計算一次系數(shù)。

搜索系數(shù)矩陣中的局部模極值。設(shè)小波系數(shù)矩陣為m×n階矩陣，其中m對應(yīng)分解尺度，n對應(yīng)光譜分辨率。coef（i，j）表示系數(shù)矩陣中i行j列的系數(shù)。則小波系數(shù)矩陣中的模極大值應(yīng)滿足以下準(zhǔn)則：

式（6）定義系數(shù)局部模極大值為大于或等于該極大值前2個系數(shù)，且大于或等于該極大值后2個系數(shù)。在尋找小波系數(shù)模極大值的同時，實現(xiàn)去噪。將式（6）中的大于等于準(zhǔn)則改為小于等于準(zhǔn)則，可以尋找出系數(shù)矩陣中的局部模極小值。尋找由各尺度局部模極大值構(gòu)成的山脊。采用文獻［7］算法構(gòu)建模極值序列，包括一組模極大值序列和兩組模極小值序列。每組模極值序列在每個尺度只取一個值。且相鄰兩個尺度的模極值都應(yīng)在一個允許的區(qū)間［nk－2k－1， nk＋2k－1］ ［17］內(nèi)尋找，其中nk是尺度k的小波系數(shù)模極大位置。具體過程如下：將最大尺度上的第一個模極大值作為山脊的起始值，并記錄該極大值對應(yīng)的時間。為防止下面幾個尺度連續(xù)在相應(yīng)的區(qū)域內(nèi)未出現(xiàn)對應(yīng)模極大值，定義了空缺參數(shù)g，初始值為0。當(dāng)下個尺度無對應(yīng)極大值時，g值加1； 若g＞3， g則值歸0，

圖3 SF6紅外光譜

Fig.3 IR spectrum for SF6

判定此山脊不成立，即該處無譜峰；否則該處可能對應(yīng)譜峰?！C重復(fù)以上步驟，搜索整個譜圖， 顯示出搜索到的所有譜峰。

4 結(jié)果與討論

選擇兩種已有的譜峰檢測算法與本算法進行比較。譜圖選取實驗室在進行SF6氣體絕緣性研究時獲取的SF6氣體紅外光譜譜圖，如圖3所示， 驗證了算法的有效性。

圖4為各種算法峰值提取的效果圖。其中A為原始信號，D為采用本算法得到的峰信息。B為文獻［6］中基于局部極大值提取算法，圖4 A： 原始譜圖； B： 基于局部極大值算法的提取結(jié)果； C： 基于小波變換算法的提取結(jié)果； D： 本算法的提取結(jié)果

Fig.4 A： Original spectrum， B： Result with Local Maximum method， C： Result with wavelet transform method， D： Result with this method

C為文獻［7，14］中的基于小波變換的峰值提取法。

基于局部極大值的峰值提取算法是利用譜圖的一階導(dǎo)數(shù)提取峰值，該算法在控制噪聲造成的干擾信息時存在困難，需要依賴前期的譜圖去噪與基線校正；并且在平頂譜峰處易造成連續(xù)辨認(rèn)，如圖4B的a處?；谖墨I［7，14］的小波變換算法能較好處理平頂譜峰，處理干擾信息的效果也優(yōu)于基于局部極大值的算法， 但該算法不能完全消除噪聲影響。由圖4B和4C中的b處可見，兩圖均將峰臂上的噪聲信號識別為譜峰，而圖4D沒有。本算法將基線校準(zhǔn)、去噪和峰值提取融為一體，降低了峰值提取算法的復(fù)雜度。由于小波系數(shù)模極大值位于信號變化處產(chǎn)生，所以能較好地處理平頂譜峰，如圖4D中的a處。

References

1 Li X， Gentleman R， Shi Q. Bioinformatics and Computational Biology Solutions Using R and Bioconductor. NewYork: New York Springer. 2005:91～109

2 Mantini D， Petrucci F， Pieragostino D， Delboccio P. BMC Bioinformatics，2007， 8（3）: 101

3 Yasui Y， Pepe M， Thompson M L. Biostatistics， 2003， 4（3）: 449～463

4 Bellew M， Coram M， Fitzgibbon M. Bioinformatics， 2006， 22（15）: 1902～1909

5 Karpievitch Y V，Hill E G， Smolka A J. Bioinformatics， 2007， 23（2）: 264～265

6 Coombes K R， Tsavachidis S， Morris J S. Proteomics， 2005， 5（16）: 4107～4117

7 Du P， Kibbe W A， Lin S M. Bioinformatics， 2006， 22（17）: 2059～2065

8 Smith C A， Want E J， Maille G O. Anal． Chem．，2006， 78（3）: 779～787

9 Coombes K R， Fritsche H A， Clarke C. Clinical Chemistry， 2003， 49（10）: 1615～1623

10 Katajamaa M， Miettinen J， Oresic M. Bioinformatics， 2006， 22（5）: 634～636

11 JIN Xiang-Jun， ZHANG Yong， XIE Yun-Fei， CONG Qian， ZHAO Bing（金向軍， 張 勇， 謝云飛， 叢 茜， 趙 冰）. Spectroscopy and Spectral Analysis（光譜學(xué)與光譜分析）， 2009，29（3）： 656～660

12 Lange E， Gropl C， Reinert K. In Pac Symp Biocomput Maui．， Hawaii， USA， 1996: 243～254

13 Baggerly K A， Morris J S， Coombers K R. Bioinformatics， 2004， 20（5）: 777～785

14 Wee A， Grayden D B， Zhu Y G. Electrophoresis， 2008， 29（20）: 4214～4225

15 Yang C， He Z Y， Yu W C. BMC Bioinformatics， 2009， 10: 4

16 Mallat S， Zhong S. IEEE Transactions on Pattern Analysis and Machine Intelligence， 1992， 14（7）: 710～732

17 Tao Wei-Liang， Wang Xian-Pei， LIU Yan， YUAN Lei（陶維亮， 王先培， 劉 艷， 袁 磊）. Spectroscopy and Spectral Analysis（光譜學(xué)與光譜分析）， 2009， 29（5）: 1241～1245

Peak Detection for Infrared Spectrum Based on

Continuous Wavelet Transform

CAI Tao*， WANG Xian-Pei， DU Shuang-Yu， YANG Jie

（Laboratory of System Integrated and Faults Diagnosis， Wuhan University， Wuhan 430079）

Abstract In spectrum analysis， peak detection is an essential step for subsequent analysis. Traditionally， peak detection procedure is divided into three consequent parts: smoothing， baseline correction and peak finding. The existing peak detection method based on continuous wavelet transform can combine the baseline correction and peak finding into one part. It simplified the traditional peak detection procedure， but this method still has two parts. A method based on continuous wavelet transform which finishes the three parts at a time was proposed in this study. The baseline′s function of original signal is monotone and linear， so after wavelet transform， there is no information of baseline in the wavelet coefficients. What we need do is deal with the coefficients. First， remove the noise in the coefficients based on Liapunov Exponent. Then， find the ridge mentioned in this study. The position of ridge is the peak′s position. The proposed method further simplifies the peak detection procedure.

Keywords Continuous wavelet transform; Infrared spectrum; Peak detection

（Received 19 August 2010; accepted 21 January 2011）