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

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

關鍵詞 小波變換; 紅外光譜; 譜峰檢測

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

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

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

1 引 言

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

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

2.1 基線校準

在進行光譜分析前，有必要對光譜信號進行基線校準。傳統的峰值提取方法是將基線校準作為單獨步驟。采用小波變換提取峰值則避免了此步驟。假設F（t）為原始光譜譜圖，B（t）為基線函數，則譜圖中的有用信號函數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）

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

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

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

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）中后一部分結果為0。因此，含基線的原始光譜譜圖在進行小波變換后，其小波系數中已不包含基線信息。所以采用小波變換法進行譜峰檢測時，無需進行基線校準。圖1為小波變換對基線移除的效果圖。

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

2.2 信號濾波

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

2.3 譜峰對應的小波變換特征

Marr小波正比于高斯函數的二階導數。令θ（t）為高斯函數。Ψ（t）為Marr小波母函數。

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

從另一個角度解釋小波變換，即將原信號f（t）與伸縮小波卷積得到（此處暫不考慮時間平移）。以Ψ（t）為小波母函數，則對信號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小波對信號進行小波變換的系數矩陣。信號的峰位處，各尺度的小波系數均有模極值存在圖中的虛線部分

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

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

據小波變換后各尺度小波系數之間的相互關系檢測譜峰。具體步驟如下：

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

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

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

圖3 SF6紅外光譜

Fig.3 IR spectrum for SF6

判定此山脊不成立，即該處無譜峰；否則該處可能對應譜峰。重復以上步驟，搜索整個譜圖， 顯示出搜索到的所有譜峰。

4 結果與討論

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

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

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］中的基于小波變換的峰值提取法。

基于局部極大值的峰值提取算法是利用譜圖的一階導數提取峰值，該算法在控制噪聲造成的干擾信息時存在困難，需要依賴前期的譜圖去噪與基線校正；并且在平頂譜峰處易造成連續辨認，如圖4B的a處。基于文獻［7，14］的小波變換算法能較好處理平頂譜峰，處理干擾信息的效果也優于基于局部極大值的算法， 但該算法不能完全消除噪聲影響。由圖4B和4C中的b處可見，兩圖均將峰臂上的噪聲信號識別為譜峰，而圖4D沒有。本算法將基線校準、去噪和峰值提取融為一體，降低了峰值提取算法的復雜度。由于小波系數模極大值位于信號變化處產生，所以能較好地處理平頂譜峰，如圖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（光譜學與光譜分析）， 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（光譜學與光譜分析）， 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）