范賢光、王秀芬、王 昕*、許英杰、闕 靖、王小東、何 浩、李 韋、左 勇

1. 廈門大學(xué)航空航天學(xué)院、福建 廈門 361005 2. 北京長城計量測試技術(shù)研究所國防科技工業(yè)第一計量測試研究中心、北京 100095

基于特征提取的低信噪比拉曼光譜去噪方法研究

范賢光1、王秀芬1、王 昕1*、許英杰1、闕 靖1、王小東1、何 浩1、李 韋1、左 勇2

1. 廈門大學(xué)航空航天學(xué)院、福建 廈門 361005 2. 北京長城計量測試技術(shù)研究所國防科技工業(yè)第一計量測試研究中心、北京 100095

為了提高拉曼光譜檢測系統(tǒng)的時間分辨率、常常需要采用較短的采樣積分時間、此時帶有分子結(jié)構(gòu)振動譜的有用拉曼信號可能完全淹沒在噪聲中、嚴(yán)重影響信號的進(jìn)一步分析、因此有必要對測量所得的光譜信號進(jìn)行噪聲消除處理。傳統(tǒng)的消噪方法是基于信號與噪聲在頻域或統(tǒng)計特性之間的差異、通過平滑濾波或取平均值的方法來消除噪聲、一般適用于噪聲強(qiáng)度不高的情況、對于信噪比較低的情況處理效果并不理想。針對傳統(tǒng)去噪方法的不足、從信號重構(gòu)的角度、利用基于小波變換的譜峰識別、半峰寬檢測提取光譜特征參數(shù)、再利用最小二乘擬合的方法、能夠有效地提取淹沒于強(qiáng)噪聲背景下的有用拉曼信號。在仿真中、運(yùn)用該算法得到的光譜曲線光滑、峰位置準(zhǔn)確、信噪比改善明顯。在實(shí)驗(yàn)中、分別利用該方法處理頭孢呋辛酯片和羅紅霉素拉曼光譜數(shù)據(jù)、得到了清晰的譜峰位置、幅值及半峰寬信息、實(shí)現(xiàn)了對短積分時間、強(qiáng)噪聲背景的拉曼信號的有效還原、提高了檢測系統(tǒng)的時間分辨率。仿真和實(shí)驗(yàn)結(jié)果表明、該方法需要調(diào)整參數(shù)少、易于實(shí)現(xiàn)、在信噪比比較低的情況下依然能夠得到良好的去噪效果、為進(jìn)一步分析光譜數(shù)據(jù)提供準(zhǔn)確可靠的信息。

拉曼光譜； 去噪； 譜峰識別； 半峰寬

引 言

拉曼光譜(Raman spectroscopy)是一種基于拉曼散射效應(yīng)的散射光譜[1]、由印度物理學(xué)家Raman發(fā)現(xiàn)。拉曼光譜技術(shù)具有非侵入、無損傷和無輻射等優(yōu)點(diǎn)、是一種強(qiáng)有力的分析工具、正廣泛地應(yīng)用于生物領(lǐng)域的活體分析[2-3]、食品安全檢測[4-5]、材料生產(chǎn)[6]和毒品快速檢測[7]等領(lǐng)域。拉曼檢測系統(tǒng)中的噪聲主要由CCD散粒噪聲、暗電流噪聲、讀出噪聲等組成。在拉曼光譜的實(shí)際應(yīng)用過程中、為了觀察快速變化的動態(tài)過程、如觀察腫瘤細(xì)胞的分裂過程等、需要采用較短的積分時間(如0.5 s)以提高光譜測量的時間分辨率。此時、獲得的拉曼信號中、帶有被測物質(zhì)的指紋信息譜很容易淹沒于強(qiáng)噪聲背景下。因此、需要通過一定的方法對強(qiáng)噪聲背景的拉曼光譜進(jìn)行去噪處理、同時、也提高了檢測系統(tǒng)的時間分辨率。

傳統(tǒng)的拉曼光譜信號噪聲消除方法有很多、如移動窗口平均平滑法、經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition、EMD)等。這些算法都能在一定程度上實(shí)現(xiàn)拉曼光譜信號的噪聲消除、然而也同樣存在各自的問題。移動窗口平均平滑法、基于信號和噪聲統(tǒng)計特性之間的差異、其基本假設(shè)是噪聲為零均值噪聲、通過對原始信號取平均值達(dá)到提高信噪比的目的、是消除噪聲最常用方法。通過選擇一個寬度為奇數(shù)2ω+1的平滑窗口、以中心波長點(diǎn)k為參照點(diǎn)從左到右移動窗口、將窗體覆蓋區(qū)域內(nèi)所有測量的平均值代替中心波長點(diǎn)對應(yīng)的測量值、直至完成所有點(diǎn)的平滑、該方法平滑窗口寬度ω影響平滑結(jié)果、寬度太小、平滑效果不佳、寬度太大、則平滑掉特征峰信息、造成光譜失真； 且存在邊界問題。對于寬度為數(shù)2ω+1的平滑窗口、光譜左右兩端各有ω個點(diǎn)不能被處理。

EMD將信號分解成有限個由高到低、時間尺度由小到大的本征模態(tài)函數(shù)(instrinsic mode function,IMF)分量及單一趨勢項(xiàng)、缺少嚴(yán)格的數(shù)學(xué)根據(jù)、更多的是一種基于經(jīng)驗(yàn)的判斷、且停止準(zhǔn)則難以確定、同時存在端點(diǎn)效應(yīng)、模態(tài)混疊效應(yīng)、算法穩(wěn)定性并不高。

傳統(tǒng)消噪方法適用于噪聲強(qiáng)度不高的情況、對淹沒于強(qiáng)噪聲背景下的拉曼信號效果不佳。因此、本文提出一種基于特征提取的光譜去噪方法。首先、通過對原始拉曼信號進(jìn)行小波變換、搜索脊線進(jìn)行譜峰位置和強(qiáng)度檢測、然后基于信號與噪聲在不同分解尺度的不同形態(tài)獲取半峰寬、最后根據(jù)獲得的峰值和半峰寬進(jìn)行最小二乘擬合、從而實(shí)現(xiàn)強(qiáng)噪聲背景下的拉曼信號信噪分離。

1 理 論

1.1 基于脊線提取的譜峰識別

采用連續(xù)Mexian Hat小波變換進(jìn)行譜峰識別、該小波定義如式(1)所示。

(1)

對拉曼信號f(x)做連續(xù)Mexian Hat小波變換有

=f(x)*φa(b)

(2)

其中、a為伸縮因子、b為平移因子、φ*表示φ的共軛。Wf(a，b)正比于f(x)與標(biāo)準(zhǔn)差為a的高斯概率密度函數(shù)卷積的二階導(dǎo)。由于拉曼光譜近似為高斯函數(shù)形式、其二階導(dǎo)數(shù)在譜峰位置必然表現(xiàn)為局部極大值。同一譜峰在不同尺度上的局部極大值位置相近、該局部最大值隨小波尺度變大而變大、且當(dāng)尺度與峰的大小接近的時候達(dá)到最大、如圖1所示、小波系數(shù)的三維圖呈現(xiàn)山脊的形狀、搜索出小波系數(shù)矩陣中的脊線即可識別出譜峰[8-9]。

圖1 CWT系數(shù)的脊線圖

小波脊線的提取有兩種最基本的方法、即基于小波系數(shù)的模極大值或相位信息的提取方法[10]。理論上講、以小波變換的相位變換涉及到相位求導(dǎo)、噪聲非常敏感、只有在比較高信噪比下才能準(zhǔn)確提取信號小波變換的脊線。因此本文根據(jù)模值信息提取小波脊線Carmona方法、通過有小尺寸的極大值向大尺寸蔓延的方式進(jìn)行脊線搜索。在獲得脊線中選取順利連接到最大尺度的脊線、運(yùn)用閾值法剔除由于噪聲引起的強(qiáng)度較弱的部分、留下的部分的最大尺度點(diǎn)即為識別的峰值點(diǎn)。由于噪聲強(qiáng)度大、在提取的譜峰中可能存在由噪聲引起的“假峰”、由于噪聲的零均值隨機(jī)性、故采用參照平均信號的方式剔除部分假峰。

另外、對信號進(jìn)行小波變換時、要求信號是無限長的、對于有限長信號、存在一個邊界延拓問題[11]。基線產(chǎn)生兩端干擾峰、因此采用邊界元素復(fù)制法、延拓兩端的信號長度。

1.2 半峰寬檢測

半峰寬是指峰高一半處的寬度、峰寬則代表譜峰兩側(cè)觀測拐點(diǎn)處所做的兩條切線與基線的兩交點(diǎn)之間的距離。目前、主要采用以Harr小波為母函數(shù)的連續(xù)小波變換求導(dǎo)的方法進(jìn)行峰寬檢測、但是理論上講、該方法得到的既不是峰寬也不是半峰寬、而是一個與譜峰相近的寬度。另外、對積分時間短、信噪低的拉曼信號進(jìn)行哈爾小波變換、得到的波形平滑度很差、因此峰寬檢測誤差也很大。

圖2 不同尺度下、單峰(a)和重疊峰(b)的小波系數(shù)

本文采用文獻(xiàn)[12]所述的方法求取半峰寬。圖2(a)和(b)分別為單峰處和重疊峰處小波變換局部極大值隨尺度變化的曲線。對于單峰、當(dāng)尺度較小時、局部極大值隨著尺度的增加而較快的增大、直到尺度與譜峰寬度相匹配、由于受到基線的影響、局部極大值隨尺度的變化趨于平緩、利用這點(diǎn)我們可以估計單峰的譜峰寬度[12]。對于重疊峰、當(dāng)尺度增大到一定程度時、由于受到相鄰譜峰的影響、小波系數(shù)會呈現(xiàn)下降趨勢、為了避開基線及相鄰峰對峰位的干擾、我們找到峰位置處的全部尺度的極大值、選取符合要求的第一個或第二個極大值作為半峰寬。

1.3 基于最小二乘的信號重構(gòu)

最小二乘法是一種由觀測數(shù)據(jù)估算線性模型中未知參數(shù)的方法,其基本思想是選擇估算量使得模型輸出與實(shí)際測量輸出之差的平方和達(dá)到最小、能有效避免正負(fù)誤差相抵、且數(shù)學(xué)處理方便。

基于譜峰的位置、強(qiáng)度和半峰寬信息、利用最小二乘法對原始拉曼信號進(jìn)行重構(gòu)。首先、單個拉曼峰可近似為高斯峰[13-14]、即

(3)

式中、γG為高斯峰的半寬、w為峰位置、v為拉曼位移。

考慮拉曼光譜本身的熒光背景、拉曼光譜可以表示為N個高斯峰G(v，[γGi，wi])(i=1,2,…,N)、基線B(v)及擬合殘差的累積形式r(v)、即

(4)

其中、ai為拉曼譜峰擬合系數(shù)、φB為基線擬合系數(shù)。

具體擬合步驟如下：

(1) 根據(jù)峰值檢測和峰寬檢測步驟確定的峰值點(diǎn)個數(shù)m、確定拉曼峰高斯擬合函數(shù)個數(shù)N=m；

(2) 根據(jù)譜峰識別獲取的峰位置Peak[i](i=1,2,…,N)和半峰寬檢測獲得的各峰的半寬Halfwidth[i](i=1,2,…,N)、按照式(4)創(chuàng)建光譜曲線擬合表達(dá)式。其中、采用二次多項(xiàng)式B(v)=b0+b1v+b2v2擬合拉曼光譜的基線。

(3) 構(gòu)造正規(guī)方程組、求解各項(xiàng)系數(shù)a1、a2、…、aN、b0、b1、b2、即得重構(gòu)的拉曼光譜表達(dá)式。

2 仿 真

根據(jù)光譜硬建模理論及拉曼譜峰近似高斯峰的假設(shè)[15]、采用高斯函數(shù)簇的線性組合逼近實(shí)際拉曼光譜曲線、構(gòu)建的仿真拉曼光譜曲線如圖3(a)所示、峰位置、峰幅值及半峰寬參數(shù)見表1所示。對仿真信號加上強(qiáng)度為60 dB噪聲后的信號如圖3(b)所示。采用本文算法對上述加噪信號進(jìn)行處理、峰值檢測小波分解尺度為scales=1∶32、捕捉到的峰位置為Peak=[300,802,1 499]、半峰寬檢測中獲取的半峰寬Halfwidth=[13,5,13]。圖3(c)中實(shí)線部分為算法處理后的信號、虛線為仿真理想信號、作為評價算法性能的參考信號。

表1 仿真信號參數(shù)表

如圖3(c)所示、采用本算法處理的信號與理想仿真信號基本上重合、譜線平滑、識別的譜峰位置準(zhǔn)確、基本保留了信號的譜峰強(qiáng)度、且峰形合理。

圖3 仿真信號(a)、噪聲信號(b)和仿真信號與去噪信號比較圖(c)

Fig.3 Ideal simulation signal (a),noisy signal (b)and Comparison between ideal signal and signal by algorithm proposed (c)

采用信噪比(SNR)和均方誤差(MSE)評價算法性能、信噪比和均方誤差定義如式(5)和式(6)所示。

(6)

表2 處理前后信號的SNR和MSE

3 實(shí)驗(yàn)部分

3.1 材料及儀器

本文采用nanophoton公司開發(fā)的第三代拉曼成像系統(tǒng)拉曼-11分別測量頭孢呋辛酯片(積分時間0.2 s)和羅紅霉素(積分時間1 s)。

3.2 結(jié)果分析

利用本文算法分別對頭孢呋辛酯片(cefuroxime axetil tablets,CA)和羅紅霉素(roxithromycin、RC)的拉曼光譜進(jìn)行處理。

圖4(a)為頭孢呋辛酯片(積分時間為0.2 s)的原始測量信號、圖4(b)為采用本文算法重構(gòu)出的去噪信號。其中、峰值檢測算法中、小波變換尺度為scales=1∶32、得到的峰值點(diǎn)為PeakCA=[71,335,487,604]、HalfwidthCA=[7,6,14,15]。

圖4 頭孢辛酯片原始拉曼信號(積分時間：0.2 s、a)和本文處理后的信號(b)

Fig.4 Original Raman signal of Cefuroxime Axetil Tablets with 0.2 s integration time(a)and signal by algorithm proposed(b)

圖5(a)為羅紅霉素(積分時間為1 s)的原始測量信號、圖5(b)為采用本文算法重構(gòu)出的去噪信號信號。其中、峰值檢測算法中、小波變換尺度為scales=1∶16、得到的峰值點(diǎn)分別為: PeakRC=[52,104,181,529]、HalfwidthRC=[14,13,8,8]。

由圖4(a)、圖5(a)可以看出、帶有被測物質(zhì)頭孢呋辛酯片和羅紅霉素指紋譜峰信息已經(jīng)完全淹沒在噪聲背景中、我們無法從原始測量拉曼光譜中提取出可供進(jìn)一步分析的有用信息。由圖4(b)、圖5(b)可知、采用本文算法重構(gòu)的光譜信號峰形合理、譜線光滑、峰位置、幅值及半峰寬信息清晰、實(shí)現(xiàn)了對短積分時間、強(qiáng)背景噪聲拉曼信號的還原、提高了檢測系統(tǒng)的時間分辨率。

圖5 羅紅霉素原始拉曼信號(積分時間：1 s、a)和本文處理后的信號(b)

Fig.5 Original Raman signal of Roxithromycin with 1 s integration time (a) and signal by algorithm proposed (b)

為了比較本文方法與傳統(tǒng)去噪方法、利用移動窗體平均平滑法處理圖4(a)所示的頭孢呋辛酯片信號和圖5(a)所示的羅紅霉素信號、結(jié)果如圖6所示。移動窗體平均平滑法的窗體寬度對平滑效果影響很大、且沒有相應(yīng)的準(zhǔn)則指導(dǎo)選取。經(jīng)過多次試驗(yàn)、移動窗體平均平滑法采用的窗體寬度為15時平滑效果相對較好。由圖6可以看出經(jīng)過移動窗體平均平滑法處理后的信號仍存在明顯的尖峰噪聲、信噪比低、帶有被測物質(zhì)分子振動的譜峰信息仍舊無法直觀獲取、而采用本文算法處理的信號譜線平滑、信噪比改善明顯、譜峰信息直觀清晰。

綜上所述、通過對原始信號取平均值的方法提高信噪比的傳統(tǒng)去噪方法對于背景噪聲強(qiáng)、信噪比低的拉曼信號去噪效果并不理想、處理得到的光譜曲線仍舊存在明顯的噪聲，

圖6 算法處理結(jié)果對比(本文算法和移動窗口濾波)、頭孢辛酯片(a)和羅紅霉素(b)

信噪比改善不大、仍舊無法提取出被測物質(zhì)的有效信息。而本文提出的基于特征提取的信號重構(gòu)方法、能夠?qū)Ρ粡?qiáng)噪聲背景淹沒的信號實(shí)現(xiàn)有效的噪聲消除、獲取的拉曼光譜曲線平滑、信噪比改善明顯、峰位置定位準(zhǔn)確、譜峰強(qiáng)度也得到了很好的保留。

4 結(jié) 論

提出了一種基于特征提取和信號重構(gòu)的強(qiáng)噪聲背景下的拉曼光譜噪聲去除方法。該算法先利用連續(xù)Mexian Hat小波變換進(jìn)行譜峰識別、然后根據(jù)峰位置附近小波變換系數(shù)隨尺度變化規(guī)律、提取拉曼譜峰半峰寬、最后根據(jù)獲取的峰位置和半峰寬、采用高斯函數(shù)作為擬合函數(shù)、利用最小二乘擬合進(jìn)行光譜重構(gòu)、實(shí)現(xiàn)了強(qiáng)噪聲背景下的拉曼光譜信噪分離、能夠有效地提取淹沒于強(qiáng)噪聲背景下的拉曼信號。與傳統(tǒng)的基于噪聲去除方法相比、本文提出的算法能夠在信噪比較低情況下、有效地重構(gòu)出被強(qiáng)背景噪聲淹沒的拉曼信號、大大提高了檢測系統(tǒng)的時間分辨率。此外、本文提出的算法調(diào)整參數(shù)少、易于實(shí)現(xiàn)、能夠很好地與現(xiàn)有的拉曼儀器軟件相結(jié)合。

[1] WEI Na,FENG Xu-qiao,ZHANG Xiao-fang,et al(韋 娜,馮敘橋,張孝芳、等). Spectroscopy and Spectral Analysis(光譜學(xué)與光譜分析),2013,33(3): 694.

[2] Keng H C,Hiro-o H、Shinsuke S. The Royal Society of Chemistry,2011,47(33): 9423.

[3] NIU Li-yuan,LIN Man-man,LI Xue,et al(牛麗媛,林漫漫,李 雪、等). Laser & Optoelectronics Progress(激光與光電子學(xué)進(jìn)展),2012,49(6): 063001.

[4] Jagtiani A V,Sawant R,Carletta J. Measurement Science and Technology,2008,19(6): 1.

[5] Ismail H B,Havva T T,Hüseyin E G,et al. Royal Society of Chemsitry,2015,5: 56606.

[6] WU Juan-xia,XU Hua,ZHANG Jin(吳娟霞,徐 華,張 錦). Acta Chimica Sinica(化學(xué)學(xué)報),2014,72(3): 301.

[7] WANG Lei,GUO Shu-xia,DAI Yin-zhen,et al(王 磊,郭淑霞,戴吟臻、等). Chinese Journal of Analytical Chemistry(分析化學(xué)研究學(xué)報),2015,43(1): 33.

[8] Du P,Kibbe W A,Lin S M. Bioinformatics,2006,22(17): 2059.

[9] JIANG Yong-hua,MO Xiao-qiang,YOU Jia-yi,et al(蔣永華,莫小強(qiáng),尤佳伊). China Measurement &Test(中國測試),2013,39(3).

[10] CHEN Yun-gu,SU Ben-yue(陳蘊(yùn)谷,蘇本躍). Chinese Journal of Quantum Electronics(量子電子學(xué)報),2012,29(6): 665.

[11] FENG Lin,JING Yi-chen,MEI Lin-li. Journal of Propulsion and Power(推進(jìn)與動力雜志),2004,20(2): 319.

[12] DENG Xian-lai,JI Guo-yi(鄧先來,紀(jì)國宜). Noise and Vibration Control(噪聲與振動控制),2012,3: 72.

[13] JIANG Cheng-zhi,SUN Qiang,LIU Ying,et al(姜承志,孫 強(qiáng),劉 英、等). Acta Optica Sinca(光學(xué)學(xué)報),2014,34(6): 299.

[14] LI Jin-rong,DAI Lian-kui,RUAN Hua(李津蓉,戴連奎,阮 華). CIESC Journal(化工學(xué)報),2012,63(7): 2128.

[15] LI Jin-rong,DAI Lian-kui,WU Xiao-li(李津蓉,戴連奎,武曉莉). Chinese Journal of Analytical Chemistry(分析化學(xué)),2014,42(10): 1518.

*Corresponding author

Research of the Raman Signal De-Noising Method Based on Feature Extraction

FAN Xian-guang1、WANG Xiu-fen1、WANG Xin1*、XU Ying-jie1、QUE Jing1、WANG Xiao-dong1、HE Hao1、LI Wei1，ZUO Yong2

1. School of Aerospace Engineering,Xiamen University,Xiamen 361005,China 2. Changcheng Institute of Metrology & Measurement,The 1st Metrology &Measurement Research Centre of National Defense Science Industry of China,Beijing 100095,China

To improve time resolution of the Raman measurement system,we need to adopt short scanning time. In this case,the weak Raman signal with vibrational spectrum of the molecular structure is easily to be buried by the high background noise,which influences the further analysis seriously. So it is necessary to de-noise the raw Raman signals. Conventional methods manage to de-noise signal by means of smoothing or averaging based on the difference between signal and noise in frequency characteristic or statistical features. They are commonly applied in the situation where the background noise is not so strong,and cannot give satisfactory results to the Raman signals with low signal-to-noise ratio. In this paper,the algorithm proposed detects peak positions and get peak half-width based on wavelet transform,and then reconstructs the Raman signals by least square fitting algorithm with characteristic parameters obtained,which extracts the useful signal from high background noise efficiently. In the simulation,the Raman curve fitted by the proposed algorithm was smooth,and the peak positions obtained were accurate,so the signal-to-noise ratio improved significantly. In the experiment,we adopted this algorithm to de-noise the tested Raman signal of Cefuroxime Axetil Tablets and Roxithromycin,respectively. The peak positions,peak half-width and amplitude were obtained and proved to be accurate. Therefore,the useful pure Raman signal could be recovered from the high background noise efficiently by the proposed algorithm,which improved the time resolution of Raman system. Both the simulation and the experiment showed that the proposed method could be easily performed with only a few parameters. Comparing with conventional methods,it could achieve satisfactory results under high background noise and provide accurate and reliable information for further analysis.

Raman spectroscopy; De-noise; Peak detection; Peak half-width

Aug. 31,2015; accepted Dec. 11,2015)

2015-08-31、

2015-12-11

國家自然科學(xué)基金項(xiàng)目(21503171)、中央高?；究蒲袠I(yè)務(wù)費(fèi)項(xiàng)目(20720150091、20720150094)、福建省高端裝備制造協(xié)同創(chuàng)新中心資金項(xiàng)目資助

范賢光、1980年生、廈門大學(xué)航天航空學(xué)院儀器與電器系副教授 e-mail: fanxg@xmu.edu.cn *通訊聯(lián)系人 e-mail: xinwang@xmu.edu.cn

TH744.1

A

10.3964/j.issn.1000-0593(2016)12-4082-06