范賢光， 王海濤， 王 昕， 許英杰， 王秀芬， 闕 靖

廈門大學(xué)物理與機(jī)電工程學(xué)院， 福建 廈門 361005

基于非均勻B樣條的拉曼光譜基線校正算法

范賢光， 王海濤， 王 昕*， 許英杰， 王秀芬， 闕 靖

廈門大學(xué)物理與機(jī)電工程學(xué)院， 福建 廈門 361005

基線校正是一種常用的消除光譜熒光干擾的方法， 是拉曼光譜數(shù)據(jù)處理的必要步驟之一。 傳統(tǒng)的多項(xiàng)式擬合基線校正算法， 簡單且易于實(shí)現(xiàn)， 但是擬合階次難以確定， 靈活性較差。 使用非均勻B樣條代替多項(xiàng)式進(jìn)行擬合， 在保留原有算法優(yōu)點(diǎn)的基礎(chǔ)上， 利用原始拉曼譜圖的峰位置信息自適應(yīng)地確定非均勻B樣條的節(jié)點(diǎn)向量， 然后以固定階次擬合光譜基線。 B樣條自身具有分段光滑的特性， 而計(jì)算樣條節(jié)點(diǎn)的節(jié)點(diǎn)向量自適應(yīng)選取算法中的峰位置信息通過使用兩次具有不同母函數(shù)的連續(xù)小波變換(continuous wavelet transform, CWT)來獲取， 既加強(qiáng)了原始光譜數(shù)據(jù)與B樣條算法本身的聯(lián)系， 也克服了傳統(tǒng)多項(xiàng)式擬合的不足。 為了驗(yàn)證本文算法的有效性， 選取了甲基對(duì)硫磷和某品牌菜籽油兩種被測物進(jìn)行實(shí)驗(yàn)， 并使用該算法進(jìn)行了基線校正， 并與兩種其他的基線校正算法與進(jìn)行了對(duì)比。 實(shí)驗(yàn)結(jié)果表明， 該方法利用固定的擬合階次就能達(dá)到較好的校正效果， 所需要的參數(shù)較少， 校正結(jié)果不會(huì)出現(xiàn)過擬合或欠擬合的現(xiàn)象， 是一種有效的拉曼光譜基線校正算法。

拉曼光譜； 基線校正； 非均勻B樣條； 節(jié)點(diǎn)向量； 峰位置

引 言

拉曼光譜(Raman spectroscopy)作為一種鑒定物質(zhì)結(jié)構(gòu)的分析測試手段， 廣泛應(yīng)用于材料、 化工、 石油、 高分子、 生物、 環(huán)保、 地質(zhì)等領(lǐng)域[1]。 由于熒光干擾等因素的存在， 光譜往往會(huì)發(fā)生較大漂移[2]， 因此， 需要通過一定的方法對(duì)拉曼光譜進(jìn)行基線校正。

目前， 主要的基線校正的方法有多項(xiàng)式擬合[3]、 小波變換[4]、 求導(dǎo)[5]等方法。 多項(xiàng)式擬合是最常用的基線校正算法， 通常使用最小二乘法， 通過循環(huán)迭代獲得拉曼光譜的基線。 此方法易于實(shí)現(xiàn)， 但精度不高， 當(dāng)譜峰較多時(shí)， 擬合階次難以確定， 容易出現(xiàn)過擬合， 并且由于沒有考慮局部性， 當(dāng)某些數(shù)據(jù)點(diǎn)漂移時(shí)， 會(huì)影響整體擬合效果。

小波變換法利用小波變換對(duì)原始光譜進(jìn)行分解， 對(duì)于分解到小波域的信號(hào)， 采取低頻置零， 高頻閾值過濾的操作， 截取到的有用信號(hào)再進(jìn)行重構(gòu)。 此方法的前提條件是， 基線相對(duì)于原信號(hào)是變化緩慢的低頻信號(hào)， 基線和光譜特征峰的頻率有明顯界限。 此外， 小波函數(shù)的選取和變換尺度都較難確定， 且對(duì)校正效果的好壞有直接影響。 求導(dǎo)的方法遵循小波變換法“信號(hào)分解-信號(hào)處理-信號(hào)重構(gòu)”這樣的處理原則， 但是校正基線后， 會(huì)改變?cè)脊庾V的形狀。

針對(duì)以上問題， 本文提出一種基于非均勻B樣條的基線校正方法， 利用小波變換尋峰算法自適應(yīng)地確定非均勻B樣條的節(jié)點(diǎn)向量， 并使用固定階次的非均勻B樣條， 通過循環(huán)迭代擬合基線， 從而達(dá)到較好的基線校正效果。

1 理 論

1.1 非均勻B樣條擬合

B樣條方法以其低階光滑的特性， 被廣泛的應(yīng)用于自由曲線曲面的造型中[6-7]。 B樣條的曲線方程為

(1)

其中，di(i=0, 1, …,n)表示第i個(gè)控制系數(shù)， B樣條曲線的控制多邊形就是由這些控制點(diǎn)依次連接而成。Ni, k(u)(i=0, 1, …,n)表示第i個(gè)k次(k+1階)B樣條基函數(shù)， 定義為

(2)

其中，ui稱為節(jié)點(diǎn)， 它是節(jié)點(diǎn)向量U={u0,u1, …,un}這個(gè)單調(diào)非減集合的一個(gè)元素。 定義域被節(jié)點(diǎn)細(xì)分為一個(gè)個(gè)的節(jié)點(diǎn)區(qū)間。 如果這些節(jié)點(diǎn)區(qū)間是均勻的， 則稱使用這些節(jié)點(diǎn)的B樣條曲線為均勻B樣條曲線， 否則， 稱之為非均勻B樣條曲線。 計(jì)算一條k次B樣條曲線， 由方程(1)可知， 共需n+1個(gè)控制系數(shù)di(i=0, 1, …,n)以及n+1個(gè)k次B樣條基函數(shù)。 控制系數(shù)一般通過最小二乘法獲得[8]， B樣條基函數(shù)則利用節(jié)點(diǎn)向量通過式(2)獲得。 因此， 節(jié)點(diǎn)向量對(duì)于計(jì)算整個(gè)B樣條曲線有著非常重要的作用， 節(jié)點(diǎn)的分布直接影響B(tài)樣條算法的擬合效果。

1.2 節(jié)點(diǎn)向量自適應(yīng)選取方法

拉曼光譜基線變化較大的位置一般發(fā)生在拉曼光譜的譜峰位置附近， 因此本文初步選取峰位置點(diǎn)作為B樣條算法的節(jié)點(diǎn)。 文獻(xiàn)[9]提出了一種利用連續(xù)小波變換實(shí)現(xiàn)拉曼譜峰識(shí)別的方法， 本文利用該方法獲得拉曼光譜的譜峰位置， 進(jìn)而獲得非均勻B樣條節(jié)點(diǎn)向量。 小波變換是一種利用小波把信號(hào)從時(shí)域轉(zhuǎn)換到小波域的手段, 素有“數(shù)學(xué)顯微鏡”的美譽(yù)[10-11]。 連續(xù)小波變換的數(shù)學(xué)表達(dá)式可以定義為

(3)

其中，s(t)是原信號(hào)，a為尺度因子，b為時(shí)間平移因子，Ψa, b(t)是小波母函數(shù)，C(a,b)為小波變換系數(shù)， 表示不同尺度不同時(shí)間平移下小波函數(shù)的線性組合。 小波系數(shù)表示了選定小波母函數(shù)與處于特定時(shí)移b和特定尺度a下與原信號(hào)的相似度。 譜峰識(shí)別的算法中， 選取墨西哥帽小波作為小波母函數(shù)

(4)

使用墨西哥帽小波作為小波母函數(shù)， 對(duì)拉曼信號(hào)進(jìn)行連續(xù)小波變換后， 在每個(gè)尺度a上的小波系數(shù)， 靠近峰中心位置附近時(shí)， 會(huì)出現(xiàn)模極大值。 并且這個(gè)模極大值隨尺度的變化而變化， 當(dāng)小波變換尺度與原譜峰峰形最匹配時(shí)， 模極大值達(dá)到最大。 將小波系數(shù)中的模極大值連接起來， 則形成如圖1所示的脊線。 該圖的橫坐標(biāo)為拉曼位移， 縱坐標(biāo)為連續(xù)小波變換尺度。 那些又長并且由較大的模極大值連接成的脊線就是較大的譜峰的位置。 由于噪聲的存在， 最初檢測到的譜峰會(huì)混雜有假峰， 因此， 需要對(duì)假峰進(jìn)行剔除處理。 文獻(xiàn)[12]介紹了譜峰位置檢測的大量準(zhǔn)則， 如信噪比、 脊線、 模極大值和峰寬等。 本文中的峰位置檢測按照基于信噪比和脊線的方法， 其中信噪比用來剔除假峰， 脊線長度閾值用來提取較大的峰位置。

獲得所有譜峰位置后， 其峰位置點(diǎn)P=[p0,p1, …,pn]還需進(jìn)行進(jìn)一步的處理， 才能作為非均勻B樣條的節(jié)點(diǎn)。 這是由于為保證B樣條擬合的效果， 其節(jié)點(diǎn)的分布不宜過密， 也不宜過稀。 過密的節(jié)點(diǎn)分布， 一方面增加擬合的計(jì)算量， 另一方面容易導(dǎo)致過擬合； 而過稀疏的節(jié)點(diǎn)分布， 容易導(dǎo)致欠擬合。 為保證節(jié)點(diǎn)的合理分布， 設(shè)定兩個(gè)參數(shù)d1和d2， 對(duì)峰位置P=[p0,p1, …,pn]進(jìn)行處理， 從而確定B樣條節(jié)點(diǎn)向量。

Fig.1 Identified ridge lines based on CWT coefficient image

節(jié)點(diǎn)向量的獲取方法可分為以下幾個(gè)步驟：

(1) 使用墨西哥帽小波作為母函數(shù)對(duì)原拉曼信號(hào)進(jìn)行連續(xù)小波變換；

(2) 連接不同尺度下的模極大值序列， 構(gòu)建脊線， 具體參考文獻(xiàn)[9]；

(3) 根據(jù)參考文獻(xiàn)[13]的方法剔除假峰， 得到最終的峰位置向量P=[p0,p1, …,pn]；

(5) 處理后的峰值點(diǎn)向量P賦值給U， 作為非均勻B樣條的節(jié)點(diǎn)向量。

1.3 非均勻B樣條基線校正的實(shí)現(xiàn)

本文利用改進(jìn)的尋峰算法確定非均勻B樣條的節(jié)點(diǎn)向量， 然后利用非均勻B樣條循環(huán)逼近曲線基線， 原始光譜扣除基線后， 即得到基線校正后的曲線。 非均勻B樣條基線校正算法的原理可由圖2表示， 整個(gè)算法的基本步驟如下：

(1) 采樣獲得m個(gè)點(diǎn)構(gòu)成的原始光譜數(shù)據(jù)S0， 令y0=S0；

(2) 使用1.2節(jié)節(jié)點(diǎn)確定算法計(jì)算非均勻B樣條節(jié)點(diǎn)向量U；

(3) 對(duì)y0使用4階(3次)非均勻B樣條擬合， 得到y(tǒng)n作為初次擬合的基線；

(4)yn與y0逐點(diǎn)比較， 取較小的點(diǎn)賦值給yn；

(5)yn與y0的相對(duì)誤差若小于迭代閾值ζ， 則判定yn為光譜曲線的基線， 否則， 令y0=yn， 返回步驟3繼續(xù)執(zhí)行， 直到滿足條件為止；

(6) 校正后的光譜數(shù)據(jù)Scorrect=S0-yn。

Fig.2 Process of baseline correction by Non-Uniform B-spline fitting

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

2.1 材料及儀器

本文所選取的被測物質(zhì)是濃度為10 ppm的甲基對(duì)硫磷和濃度為10 ppm某品牌菜籽油。 所使用的實(shí)驗(yàn)儀器為Ocean Optics公司的QE65Pro。

2.2 結(jié)果分析

利用本文算法對(duì)甲基對(duì)硫磷(parathion-methyl， PM)和菜籽油(colza oil， CO)進(jìn)行處理， 兩種物質(zhì)的原始拉曼譜圖和擬合出的基線如圖3所示。 其中， 峰值檢測算法參數(shù)為： 小波變換尺度scales=1∶70, 剔除脊線的參數(shù)SNR=5和ridgeLength=10， 得到的峰值點(diǎn)分別為

PPM=[481, 614, 730, 828, 913, 1 006, 1 076, 1 196, 1 275, 1 357, 1 433, 1 505, 1 646, 2 370]

PCO=[390, 840, 972, 1 093, 1 255, 1 292, 1 436, 1 655, 1 736, 1 869, 1 889, 1 941, 1 975, 2 025, 2 087, 2 146, 2 183, 2 235, 2 342, 2 403]

選取的循環(huán)迭代閾值ζ=0.5， 非均勻B樣條擬合階數(shù)為4。 按照上述給定參數(shù)， 根據(jù)本文算法擬合出的甲基對(duì)硫磷和菜籽油的拉曼譜圖的基線如圖3所示。 由圖3(a)可以看出， 本文所測的甲基對(duì)硫磷的拉曼譜圖在1 000～1 500 cm-1附近峰的分布較密， 拉曼位移從1 700～2 500 cm-1沒有明顯的拉曼峰。 根據(jù)峰值檢測算法計(jì)算出的甲基對(duì)硫磷的峰值點(diǎn)PPM， 峰值較密附近存在多達(dá)7個(gè)峰值點(diǎn)， 1 700～2 500 cm-1僅有1個(gè)峰值點(diǎn)， 因此PPM不能直接作為B樣條節(jié)點(diǎn)向量， 需要對(duì)PPM進(jìn)行進(jìn)一步處理。

根據(jù)節(jié)點(diǎn)向量自適應(yīng)選取方法得到參數(shù)d1=100，d2=250， 則處理后的非均勻B樣條節(jié)點(diǎn)為

UPM=[481, 614, 730, 913, 1 196, 1 357, 1 505, 1 646, 1 846, 2 180, 2 370]

UCO=[390, 615, 840, 972, 1 093, 1 255, 1 436, 1 655, 1 869, 1 975, 2 087, 2 161, 2 235, 2 342]

處理后的B樣條節(jié)點(diǎn)向量UPM， 節(jié)點(diǎn)分布合理， 得到的校正基線也較為平滑， 沒有在峰位置過密的1 000～1 500 cm-1區(qū)域出現(xiàn)過擬合， 也沒有在峰值點(diǎn)較少的1 700～2 500 cm-1上出現(xiàn)欠擬合。 校正后的甲基對(duì)硫磷的拉曼光譜圖如圖4(a)所示。 本文算法在進(jìn)行基線校正時(shí)， 并沒有對(duì)原始拉曼數(shù)據(jù)進(jìn)行平滑去噪等預(yù)處理操作， 圖3(b)為實(shí)驗(yàn)所測的菜籽油的拉曼譜圖， 箭頭標(biāo)示處存在高頻噪聲， 對(duì)比圖4(b)菜籽油基線校正后的拉曼譜圖， 可以發(fā)現(xiàn)菜籽油的譜圖的高頻噪聲處基線擬合依然平滑， 校正后的譜圖也沒有發(fā)生過擬合的現(xiàn)象， 證明了本算法的有效性。

Fig.3 Raman spectroscopy and its baseline by Non-Uniform B-spline fitting

選取傳統(tǒng)多項(xiàng)式擬合算法與改進(jìn)的懲罰最小二乘算法[13]與本文算法進(jìn)行比較。 圖5給出了使用5階多項(xiàng)式擬合的基線校正算法所得出的菜籽油拉曼譜圖及其校正基線。 所

Fig.4 Raman spectroscopy after baseline correction by Non-Uniform B-spline fitting

Fig.5 Raman spectroscopy and its baseline by polynomial fitting of colza oil

測得的菜籽油拉曼譜圖漂移較大， 且存在高頻噪聲， 利用5階多項(xiàng)式擬合菜籽油基線時(shí)， 800 cm-1處開始的基線出現(xiàn)明顯的過擬合現(xiàn)象， 1 800 cm-1處又出現(xiàn)欠擬合現(xiàn)象。 圖6為利用改進(jìn)的懲罰最小二乘法獲得的甲基對(duì)硫磷的校正基線及其校正后的拉曼譜圖。 文獻(xiàn)[13]將該方法做成了R語言開源包baselineWavelet。 該方法先利用兩個(gè)小波變換獲取峰的信息， 然后在峰位置處使用直線連接峰首尾作為基線， 非峰部分使用懲罰最小二乘法擬合基線， 拼接后得到最終基線。 此方法計(jì)算出的基線與原始譜圖結(jié)合緊密， 不易出現(xiàn)過擬合的現(xiàn)象。 但是， 由于基線是兩種方法拼接而成， 該方法的基線整體不夠平滑， 如圖6(a)所示。 比較圖4(a)和圖6(b)， 本文算法與改進(jìn)的懲罰最小二乘算法校正效果差別不大， 但是， 此方法相較于本文算法， 默認(rèn)參數(shù)眾多， 兩次小波變換計(jì)算量也較大。

Fig.6 Raman spectroscopy of parathion-methyl

(a): Raman spectroscopy and its baseline; (b): Raman spectroscopy after baseline correction by baseline wavelet

綜上所述， 采用基于非均勻B樣條的基線校正方法， 充分利用了B樣條的低階光滑特性， 克服了傳統(tǒng)多項(xiàng)式基線校正算法擬合階數(shù)難確定的缺陷， 在基線偏移嚴(yán)重和存在高頻噪聲時(shí)， 仍能夠擬合出光滑的基線， 并沒有出現(xiàn)欠擬合和過擬合的現(xiàn)象， 對(duì)于存在基線漂移的拉曼光譜信號(hào)， 達(dá)到了較好的基線校正效果。

3 結(jié) 論

提出了一種基于非均勻B樣條的拉曼光譜基線校正方法， 先利用基于小波的譜峰識(shí)別算法獲得拉曼譜圖的峰位置， 然后對(duì)這些峰位置點(diǎn)進(jìn)行處理， 從而確定出非均勻B樣條的節(jié)點(diǎn)向量， 最后利用B樣條以固定的階次循環(huán)迭代擬合拉曼譜圖的基線。 與傳統(tǒng)的多項(xiàng)式擬合算法相比， 本文的自適應(yīng)節(jié)點(diǎn)確定算法充分利用了不同拉曼譜圖的信息， 讓擬合出的基線與原始拉曼譜圖結(jié)合更緊密； 同時(shí)， 本文的方法能夠以固定的階次擬合光譜基線， 并且能夠有效地避免過擬合和欠擬合。 因此， 本文所提出的方法， 能夠較好地實(shí)現(xiàn)拉曼光譜的基線校正， 并且調(diào)整參數(shù)少、 易于實(shí)現(xiàn)， 可以作為一種有效的拉曼光譜基線校正算法。

[1] CHU Xiao-li(褚小立). Molecular Spectroscopy Analytical Technology Combined with Chemometrics and Its Applications(化學(xué)計(jì)量學(xué)方法與分子光譜分析技術(shù)). Beijing: Chemical Industry Press(北京： 化學(xué)工業(yè)出版社), 2011. 311.

[2] McCreer R L. Raman Spectroscopy for Chemical Analysis. Wiley-Interscience, 2000.

[3] FENG Xin-wei, ZHU Zhong-liang, SHEN Meng-jie, et al. Computer and Applied Chemsity, 2009, 26(6): 759.

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

[5] O’Grady A, Dennis AC, Denvir D. Analytical Chemistry, 2001, 73: 2058.

[6] Piegl L, Tiller W. The NURBS Book(非均勻有理B樣條). Translated by ZHAO Gang, MU Guo-wang, WANG La-zhu(趙 罡， 穆國旺， 王拉柱， 譯). Beijing: Tsinghua University Press(北京： 清華大學(xué)出版社), 2010. 60.

[7] GUO Jian-feng, ZHU Chang-qing(郭建峰， 朱長青). Engineering of Surveying and Mapping(測繪工程)， 2000, 9(4): 25.

[8] SHI Fa-zhong(施法中). CAGD & NURBS(計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理B樣條). Beijing: Higher Education Press(北京： 高等教育出版社), 2013. 303.

[9] Du P, Kibbe W A, Lin S M. Bioinformatics, 2006, 22: 2059.

[10] Daubechies I. Ten Lectures on Wavelets(小波十講). Translated by LI Jian-ping(李建平， 譯). Beijing: National Defense Industry Press(北京： 國防工業(yè)出版社)， 2011.

[11] SUN Yan-kui(孫延奎). Wavelet Analysis and Applications. Beijing: China Machine Press， 2005.

[12] Yang C, He Z Y, Yu W C. BMC Bioinformatics， 2009. 10.

[13] Zhang Z M, Chen S, Liang Y Z, et al. Journal of Raman Spectroscopy, 2010, 41: 659.

*Corresponding author

Baseline Correction Algorithm for Raman Spectroscopy Based on Non-Uniform B-Spline

FAN Xian-guang, WANG Hai-tao, WANG Xin*, XU Ying-jie, WANG Xiu-fen, QUE Jing

School of Physics and Mechanical & Electrical Engineering, Xiamen University, Xiamen 361005, China

As one of the necessary steps for data processing of Raman spectroscopy, baseline correction is commonly used to eliminate the interference of fluorescence spectra. The traditional baseline correction algorithm based on polynomial fitting is simple and easy to implement, but its flexibility is poor due to the uncertain fitting order. In this paper, instead of using polynomial fitting, non-uniform B-spline is proposed to overcome the shortcomings of the traditional method. Based on the advantages of the traditional algorithm, the node vector of non-uniform B-spline is fixed adaptively using the peak position of the original Raman spectrum, and then the baseline is fitted with the fixed order. In order to verify this algorithm, the Raman spectra of parathion-methyl and colza oil are detected and their baselines are corrected using this algorithm, the result is made comparison with two other baseline correction algorithms. The experimental results show that the effect of baseline correction is improved by using this algorithm with a fixed fitting order and less parameters, and there is no over or under fitting phenomenon. Therefore, non-uniform B-spline is proved to be an effective baseline correction algorithm of Raman spectroscopy.

Raman spectroscopy; Baseline correction; Non-uniform B-spline; Node vector; Peak position

Dec. 23, 2014; accepted Apr. 18, 2015)

2014-12-23，

2015-04-18

國家重大科學(xué)儀器設(shè)備開發(fā)專項(xiàng)(2011YQ03012417)資助

范賢光， 1980年生， 廈門大學(xué)物理與機(jī)電工程學(xué)院機(jī)電工程系副教授 e-mail: fanxg@xmu.edu.cn *通訊聯(lián)系人 e-mail: xinwang@xmu.edu.cn

O657.3

A

10.3964/j.issn.1000-0593(2016)03-0724-05