王昕 康哲銘 劉龍 范賢光 ?
1) (廈門大學(xué)儀器與電氣系, 廈門 361005)
2) (福建省高等院校傳感器技術(shù)重點(diǎn)實(shí)驗(yàn)室, 廈門市光電子傳感器技術(shù)重點(diǎn)實(shí)驗(yàn)室, 廈門 361005)
拉曼光譜由印度物理學(xué)家Raman于1919年從水分子的散射現(xiàn)象中首次發(fā)現(xiàn), 并將其發(fā)表于雜志《Nature》上[1]. 拉曼光譜作為一種用于分析分子化學(xué)成分、結(jié)構(gòu)等信息的檢測技術(shù), 具有無侵入、特異性高、無標(biāo)記、無電離輻射、不受水的干擾等優(yōu)點(diǎn), 被廣泛應(yīng)用于生物醫(yī)學(xué)、材料生產(chǎn)、化學(xué)化工等領(lǐng)域[2,3]. 然而由于熒光干擾等因素的存在,基線漂移現(xiàn)象通常存在于光譜信號中, 而基線校正算法是解決該現(xiàn)象的必要手段[4?7].
目前, 常用的基線校正算法有多項(xiàng)式擬合算法(polynomial fitting)、B樣條擬合算法(B-spline fitting)、小波變換 (wavelet transform)和自適應(yīng)迭代重加權(quán)懲罰最小二乘(adaptive iterative reweighted penalty least square, airPLS)等算法[8?13].多項(xiàng)式擬合算法簡單易實(shí)現(xiàn), 但是在處理不同的光譜信號時(shí)難以確定擬合階數(shù). 小波變換則是通過用一系列的正弦波進(jìn)行疊加來代替每條拉曼光譜, 能夠有效地將拉曼光譜的低頻和高頻進(jìn)行區(qū)分. 但是, 由于不同的光譜具有不同的頻率和噪聲, 很難找出一個(gè)適用于所有拉曼光譜的分解辦法, 且其計(jì)算量相比其他算法也相對較高. 均勻B樣條擬合算法在面對基線漂移較為平緩的拉曼光譜信號時(shí),其基線校正效果良好. 然而由于均勻B樣條內(nèi)節(jié)點(diǎn)的選擇缺乏目的性, 隨著基線漂移程度越來越劇烈, 通常導(dǎo)致最終基線校正效果一般. 而若非均勻選擇內(nèi)節(jié)點(diǎn), 則需要針對每條光譜信號進(jìn)行手動(dòng)選擇, 工作量也會隨之增加, 且普適性和靈活性較差.范賢光等[14]通過小波變換找到光譜數(shù)據(jù)的峰值并作為內(nèi)結(jié)點(diǎn), 該方法可以解決內(nèi)結(jié)點(diǎn)選擇的問題,但其計(jì)算復(fù)雜度和計(jì)算量較高, 且過擬合和欠擬合仍不能完全消除.
因此, 本文提出一種基于中值濾波和非均勻B 樣條 (median filtering and un-uniform B-spline,MF-UUB)的拉曼光譜基線校正算法. 該算法不僅可以根據(jù)不同的拉曼光譜信號自適應(yīng)地選擇內(nèi)結(jié)點(diǎn), 且通過使用中值濾波和非均勻B樣條算法能獲得更好的基線校正效果.
B樣條算法因其具有低階平滑性而廣泛應(yīng)用于曲線和曲面擬合中[15]. 在本文中, 將該算法用于拉曼光譜的基線擬合.
假設(shè)輸入拉曼光譜數(shù)據(jù)r, 用節(jié)點(diǎn)符號t將數(shù)據(jù)點(diǎn)列在x軸上進(jìn)行區(qū)間劃分, 點(diǎn)序列如下所示:
其中t1到tN為內(nèi)節(jié)點(diǎn), 其余的為外節(jié)點(diǎn). 如果內(nèi)節(jié)點(diǎn)在區(qū)間內(nèi)均勻劃分則為均勻B樣條, 反之, 則為非均勻B樣條[16]; 而外節(jié)點(diǎn)對算法的影響較小, 其區(qū)間大小一般取內(nèi)節(jié)點(diǎn)區(qū)間的平均值.N為區(qū)域內(nèi)內(nèi)節(jié)點(diǎn)數(shù),k為B 樣條的階數(shù), 且一般取4.x的k階B樣條曲線如(2)式所示:
式中ci為控制系數(shù)列向量, B樣條曲線的控制多邊形就是由這些控制點(diǎn)依次連接而成;為B樣條基函數(shù), 定義為[17]:
局部加權(quán)線性回歸 (locally weighted linear regression, LWLR)[18]是線性回歸方法的一個(gè)改進(jìn), 可有效應(yīng)用于光譜信號的平滑. 因此, 本文通過使用局部加權(quán)線性回歸對原始光譜數(shù)據(jù)進(jìn)行平滑預(yù)處理, 以便后續(xù)有效地進(jìn)行內(nèi)節(jié)點(diǎn)的選取. 它通過預(yù)測數(shù)據(jù)與輸入數(shù)據(jù)的距離賦予不同權(quán)重, 距離越近則權(quán)重越大, 然后再通過最小二乘法進(jìn)行計(jì)算求解. 該方法的本質(zhì)是基于加權(quán)的均方誤差最小化進(jìn)行求解:
式中r為輸入的光譜數(shù)據(jù);y為預(yù)測數(shù)據(jù);θ為回歸系數(shù);wi為權(quán)重, 該值可通過權(quán)重函數(shù)計(jì)算得到,一般選用高斯函數(shù)作為權(quán)重函數(shù), 權(quán)重函數(shù)如下所示:
式中δ為帶寬參數(shù).
通過非均勻B樣條算法進(jìn)行基線擬合, 其擬合效果的好壞主要取決于原始拉曼光譜的基線能否盡可能地被內(nèi)節(jié)點(diǎn)劃分為幾個(gè)平緩的波段. 因此, 如何根據(jù)不同的拉曼光譜自適應(yīng)地選取內(nèi)節(jié)點(diǎn), 則為非均勻B樣條的核心問題. 為保證扣除基線后仍能最大程度的保留原始光譜的特征峰, 且盡可能地在基線變化較為劇烈的位置處選擇內(nèi)節(jié)點(diǎn),同時(shí)有效地將基線劃分為幾個(gè)平緩的波段. 在本文中, 內(nèi)節(jié)點(diǎn)的選取則根據(jù)原始光譜信號的波谷點(diǎn)進(jìn)行確定, 因此, 內(nèi)節(jié)點(diǎn)的選取步驟如下.
1)使用局部加權(quán)線性回歸對原始光譜數(shù)據(jù)進(jìn)行平滑預(yù)處理. 由于光譜數(shù)據(jù)中伴隨著隨機(jī)噪聲,若直接對光譜數(shù)據(jù)進(jìn)行篩選波谷點(diǎn), 則會導(dǎo)致篩選后的數(shù)據(jù)點(diǎn)中存在無效數(shù)據(jù)點(diǎn). 局部加權(quán)線性回歸的窗口長度可通過設(shè)置窗口數(shù)量計(jì)算得到, 窗口數(shù)量一般取102個(gè). 原始光譜數(shù)據(jù)與平滑后的數(shù)據(jù)如圖1所示.
2)篩選光譜數(shù)據(jù)的波谷點(diǎn). 對平滑后的光譜信號的每個(gè)數(shù)據(jù)點(diǎn)逐個(gè)進(jìn)行差分計(jì)算, 并設(shè)置兩波谷點(diǎn)間的最小距離閾值ξL, 從而初步篩選出波谷點(diǎn)(圖1中紅色數(shù)據(jù)點(diǎn))并設(shè)為K=[k0,k1,···,kn].最小距離閾值ξL通過內(nèi)節(jié)點(diǎn)預(yù)設(shè)值計(jì)算得到, 如(7)式:
式中sik為內(nèi)節(jié)點(diǎn)數(shù)量預(yù)設(shè)值, 一般設(shè)為30;P為原始光譜數(shù)據(jù)點(diǎn)數(shù); R ound[·]為四舍五入取整. 將從平滑后的光譜數(shù)據(jù)中篩選的波谷點(diǎn)K對應(yīng)至原始光譜數(shù)據(jù)上. 由于對應(yīng)點(diǎn)并非都準(zhǔn)確地位于原始數(shù)據(jù)的波谷上, 因此對每個(gè)數(shù)據(jù)點(diǎn)進(jìn)行加窗, 一般窗口距離取0.2倍的最小距離閾值ξL, 并求得每個(gè)窗口內(nèi)的最小值. 求得的數(shù)據(jù)點(diǎn)則為原始光譜數(shù)據(jù)的波谷點(diǎn)(圖1中藍(lán)色數(shù)據(jù)點(diǎn)), 并設(shè)為P=[p0,p1,···,pn].
3)扣除無效的數(shù)據(jù)點(diǎn). 計(jì)算兩波谷點(diǎn)區(qū)間內(nèi)最大值與最小值之差, 設(shè)為hi, 并設(shè)置最小高度閾值ξH, 其值一般取0.05倍的原始光譜數(shù)據(jù)中最大值與最小值之差:
式中α值一般取 0.05,r為原始拉曼光譜數(shù)據(jù). 若波谷點(diǎn)pi和pi+1 之間的高度差hi小于閾值ξH , 則扣除數(shù)據(jù)點(diǎn)pi+1. 無效的數(shù)據(jù)點(diǎn)(圖2中綠色數(shù)據(jù)點(diǎn))如圖2所示, 剩余的有效數(shù)據(jù)點(diǎn)設(shè)為d=[d0,d1,···,dn].
圖 1 波谷點(diǎn)選擇策略Fig. 1. Strategy of selection of trough points.
圖 2 內(nèi)節(jié)點(diǎn)選擇策略Fig. 2. Strategy of selection of internal knots.
4) 插入額外數(shù)據(jù)點(diǎn). 計(jì)算兩數(shù)據(jù)點(diǎn)di和di+1間的距離, 并設(shè)為li; 若li大于內(nèi)節(jié)點(diǎn)距離閾值 (3 倍的最小距離閾值ξL), 則在數(shù)據(jù)點(diǎn)di和di+1間插入另一數(shù)據(jù)點(diǎn)qi. 循環(huán)執(zhí)行該步驟, 直至兩數(shù)據(jù)點(diǎn)間的距離均小于內(nèi)節(jié)點(diǎn)距離閾值. 插入的數(shù)據(jù)點(diǎn)(圖2中紅色數(shù)據(jù)點(diǎn))設(shè)為q=[q0,q1,···,qn], 如圖2所示.
5)使用確定的數(shù)據(jù)點(diǎn)作為非均勻B樣條的內(nèi)節(jié)點(diǎn). 數(shù)據(jù)點(diǎn)d和q之和得到的數(shù)據(jù)點(diǎn)T=[t0,t1,···,tn]則作為內(nèi)節(jié)點(diǎn)(圖2中黑色與紅色數(shù)據(jù)點(diǎn)之和).
中值濾波(median filtering)[19]是一種非線性的平滑濾波技術(shù), 它通過把一串?dāng)?shù)字序列或者數(shù)字圖像中的某一個(gè)點(diǎn)的值用其相鄰區(qū)域內(nèi)所有值的中值進(jìn)行代替. 在本文中, 利用中值濾波算法對原始拉曼光譜數(shù)據(jù)進(jìn)行處理, 使光譜數(shù)據(jù)中的波峰趨于平緩, 從而優(yōu)化B樣條曲線的控制系數(shù)列向量ci,以便后續(xù)的非均勻B樣條能在波峰過渡到平滑波段的位置處更好地進(jìn)行基線擬合. 使其擬合的基線在該位置處變化更加平緩, 能有效地減少欠擬合現(xiàn)象的發(fā)生. 中值濾波的公式為:
式中 M ed[·]為取輸入數(shù)據(jù)點(diǎn)列的中值. 此外, 中值濾波的窗口數(shù)量越少則擬合的基線越趨于平緩, 但窗口數(shù)量過少則會影響擬合基線的準(zhǔn)確性. 因此,窗口數(shù)量一般取與內(nèi)節(jié)點(diǎn)數(shù)量一致.
將內(nèi)節(jié)點(diǎn)選取方法、中值濾波算法和非均勻B樣條擬合等算法應(yīng)用于拉曼光譜的熒光背景扣除, 該算法步驟如下:
1)使用局部加權(quán)線性回歸對原始光譜數(shù)據(jù)r進(jìn)行平滑預(yù)處理, 得到數(shù)據(jù)z;
2)對數(shù)據(jù)z進(jìn)行逐點(diǎn)差分計(jì)算和設(shè)置最小距離閾值ξL得到數(shù)據(jù)z的波谷點(diǎn), 并通過加窗處理得到數(shù)據(jù)r的波谷點(diǎn)P=[p0,p1,···,pn];
3)利用最小高度閾值ξH扣除無效波谷點(diǎn), 并設(shè)置內(nèi)節(jié)點(diǎn)距離閾值進(jìn)行插值, 得到內(nèi)節(jié)點(diǎn)T=[t0,t1,···,tn];
4)將內(nèi)節(jié)點(diǎn)數(shù)量作為中值濾波窗口數(shù)量, 對數(shù)據(jù)r進(jìn)行中值濾波處理, 得到數(shù)據(jù)u;
5)對數(shù)據(jù)u進(jìn)行非均勻B樣條擬合, 得到擬合基線v;
6)計(jì)算u和v的誤差絕對值, 若誤差值大于所設(shè)定的閾值ξ(閾值ξ一般取0.05), 則將v賦值給u并回到步驟4); 誤差絕對值公式為
7)從原始拉曼光譜r中扣除基線v.
由上述流程可知, 該算法需要調(diào)節(jié)的參數(shù)僅包含最小距離閾值ξL和局部加權(quán)線性回歸窗口數(shù)量.其中最小距離閾值ξL通過設(shè)置內(nèi)節(jié)點(diǎn)數(shù)量預(yù)設(shè)值sik計(jì)算得到. 且內(nèi)節(jié)點(diǎn)數(shù)量預(yù)設(shè)值sik一般取30個(gè), 局部加權(quán)線性回歸窗口數(shù)量一般取102個(gè)即可獲得良好的基線校正效果.
本文選用聚甲基丙烯酸甲酯(polymethyl methacrylate, PMMA)和 正 辛 烷 (N-OCTANE)作為待測樣品. N-OCTANE拉曼光譜在光譜儀積分時(shí)間為1 s下獲得, PMMA拉曼光譜在光譜儀積分時(shí)間為10 s下獲得.
圖3(a),(b)分別是N-OCTANE和PMMA的原始拉曼光譜, 其中N-OCTANE拉曼光譜的基線漂移程度較小, 而PMMA拉曼光譜的基線漂移程度則較為劇烈, 尤其在拉曼位移為 1500 cm–1處,基線的存在會淹沒原本就微弱的譜峰, 這些現(xiàn)象使得光譜信號難以被直觀地辨識. 通過本文所提的MF-UUB基線校正算法對兩種樣品的拉曼光譜分別進(jìn)行基線擬合, 得到的擬合基線同樣如圖3所示, 基線校正后的拉曼光譜如圖 4所示. 其中, 局部加權(quán)線性回歸的窗口數(shù)量設(shè)為102個(gè), 內(nèi)節(jié)點(diǎn)數(shù)量預(yù)設(shè)值設(shè)為30個(gè), 因此對于兩條拉曼光譜分別得到的內(nèi)節(jié)點(diǎn)數(shù)量為14個(gè)和17個(gè); 中值濾波的窗口數(shù)量分別與各自光譜數(shù)據(jù)的內(nèi)節(jié)點(diǎn)數(shù)量一致. 由圖3可知, 該算法能有效地隨波峰的變化進(jìn)行基線擬合, 且在線性波段和曲線波段的擬合中均有較好的擬合效果, 這證明了MF-UUB算法的靈活性.由圖4可知, 基線校正后的拉曼光譜能有效地保留所有特征峰信息, 且沒有因?yàn)闊晒獗尘暗目鄢霈F(xiàn)多余的無效特征峰或凸包. 對于大部分的拉曼光譜, 使用固定的參數(shù)設(shè)置, 便能獲得較好的基線擬合效果. 在后續(xù)對拉曼信號的進(jìn)一步分析中, 如拉曼組分識別的應(yīng)用, 識別的過程主要根據(jù)譜峰的位置等特點(diǎn)進(jìn)行識別, 而若無法有效且完整地扣除熒光背景, 將使得某些譜峰的位置發(fā)生偏移, 或在一定程度上改變特征峰形狀及大小, 這將影響對物質(zhì)的有效識別. 而本文所提算法則能有效地避免上述問題, 能為后續(xù)的進(jìn)一步分析提供更可靠有效的信息.
圖 3 原始拉曼光譜和基于 MF-UUB 算法擬合的基線(a) N-OCTANE 拉曼光譜; (b) PMMA 拉曼光譜Fig. 3. Original Raman spectra and fitted baseline by MFUUB: (a) N-OCTANE Raman spectrum; (b) PMMA Raman spectrum.
圖 4 基于 MF-UUB 算法基線校正后的光譜 (a) N-OCTANE拉曼光譜; (b)PMMA拉曼光譜Fig. 4. Corrected Raman spectra by MF-UUB: (a) N-OCTANE Raman spectrum; (b) PMMA Raman spectrum.
為驗(yàn)證中值濾波算法是否能有效地提升擬合基線的準(zhǔn)確性, 通過使用不加中值濾波處理的非均勻B樣條擬合算法進(jìn)行基線擬合, 擬合的基線如圖5所示. 由圖5(a)可知, 擬合的基線與本文所提算法擬合的基線在整體上基本一致. 但是在波峰過渡到平滑波段的位置, 如拉曼位移為1400—2000 cm–1等處, 擬合的基線存在一定的欠擬合現(xiàn)象. 而通過使用中值濾波算法處理后, 則可以有效地改善這一缺陷. 這是因?yàn)橥ㄟ^中值濾波處理后的拉曼光譜的波峰會趨于平緩, 使得擬合的基線變化幅度變小, 從而能更好地在波峰過渡到平滑波段的位置處進(jìn)行擬合.
圖 5 原始拉曼光譜和基于非均勻B樣條算法擬合的基線 (a) N-OCTANE 拉曼光譜; (b)PMMA 拉曼光譜Fig. 5. Original Raman spectra and fitted baseline by ununiform B-spline algorithm: (a) N-OCTANE Raman spectrum; (b) PMMA Raman spectrum.
此外, 如果通過非均勻B樣條擬合算法便能獲得較好的擬合基線, 則中值濾波算法不會對擬合的基線有太大的改變. 由圖5(b)可知, 該光譜的基線擬合效果好于圖5(a), 且不存在欠擬合現(xiàn)象. 因此, 通過中值濾波進(jìn)行處理后再進(jìn)行基線擬合的效果則提升相對有限, 僅在拉曼位移為1000—1300 cm–1處存在一定程度的改善.
為進(jìn)一步驗(yàn)證該算法能有效地進(jìn)行基線擬合且優(yōu)于傳統(tǒng)算法, 分別使用多項(xiàng)式擬合算法和均勻B樣條擬合算法對N-OCTANE和PMMA原始拉曼光譜進(jìn)行基線擬合. 此外, 由于多項(xiàng)式擬合算法的擬合階數(shù)難以確定, 因此同時(shí)選用5階、7階和9階多項(xiàng)式擬合算法對原始光譜進(jìn)行基線擬合. 而均勻B樣條擬合算法則采用4階3次B樣條, 且為保證驗(yàn)證的有效性, 均勻B樣條的內(nèi)節(jié)點(diǎn)數(shù)與本文所提算法保持一致, 分別設(shè)為14個(gè)和17個(gè).基于均勻B樣條算法和多項(xiàng)式擬合算法擬合的基線分別如圖6和圖7所示.
由圖6可知, 均勻B樣條擬合算法在整體上能有效地進(jìn)行基線擬合. 但由于內(nèi)節(jié)點(diǎn)的選擇缺乏目的性, 無法有效地在正確位置處將光譜數(shù)據(jù)的基線劃分為幾個(gè)平滑波段, 這導(dǎo)致擬合的基線存在欠擬合現(xiàn)象, 如圖 6(a)中拉曼位移為 300—700 cm–1以及圖6(b)中拉曼位移為1000—1200 cm–1等多處波段中. 由于無法有效且完整地扣除熒光背景,將導(dǎo)致扣除熒光背景后的拉曼光譜在某些波段處出現(xiàn)無效的凸包, 且在一定程度上改變特征峰的形狀.
由圖7可知, 多項(xiàng)式擬合算法在N-OCTANE原始拉曼光譜的基線擬合中表現(xiàn)效果好于PMMA原始拉曼光譜, 能相對有效地進(jìn)行基線擬合, 但整體上仍存在著一定的擬合缺陷. 隨著多項(xiàng)式擬合階數(shù)的增加, 擬合的基線在某些波段處會得到改善,如圖7(a)中拉曼位移為700—1200 cm–1以及圖7(b)中拉曼位移為 1100—1800 cm–1處. 但同時(shí)隨著擬合階數(shù)的增加, 擬合的基線在某些波段處的擬合效果也會變差, 如圖7(a)中拉曼位移為1200—1600 cm–1以及圖 7(b)中拉曼位移為 800—1100 cm–1處. 因此可以得出, 多項(xiàng)式擬合算法的階數(shù)難以確定, 且欠擬合和過擬合現(xiàn)象難以控制和預(yù)料. 針對基線漂移較為劇烈的原始拉曼光譜, 通過多項(xiàng)式擬合算法擬合的基線效果較為一般, 無法有效且完整地扣除熒光背景.
圖 6 原始拉曼光譜和基于均勻B樣條算法擬合的基線(a) N-OCTANE 原始拉曼光譜; (b) PMMA 原始拉曼光譜Fig. 6. Original Raman spectra and fitted baseline by uniform B-spline algorithm: (a) N-OCTANE Raman spectrum;(b) PMMA Raman spectrum.
圖 7 原始拉曼光譜和基于多項(xiàng)式擬合算法擬合的基線(a) N-OCTANE 拉曼光譜; (b)PMMA 拉曼光譜Fig. 7. Original Raman spectra and fitted baseline by polynomial fitting algorithm: (a) N-OCTANE Raman spectrum;(b) PMMA Raman spectrum.
此外, airPLS是一種應(yīng)用廣泛且效果良好的基線校正算法[20], 因此, 通過該算法擬合的基線如圖8所示. 由圖8可知, 該算法整體基線擬合效果良好, 在大部分波段處均能準(zhǔn)確有效地進(jìn)行基線擬合, 且擬合的基線與通過本文所提算法擬合的基線基本一致. 但是, 仍然有一些擬合缺陷產(chǎn)生, 如擬合的基線通過波谷點(diǎn)時(shí)存在著一定的過擬合現(xiàn)象.
圖 8 原始拉曼光譜和基于 airPLS 擬合的基線 (a) NOCTANE拉曼光譜; (b)PMMA拉曼光譜Fig. 8. Original Raman spectra and fitted baseline by air-PLS algorithm: (a) N-OCTANE Raman spectrum; (b) PMMA Raman spectrum.
綜上所述, 本文提出MF-UUB算法能有效地進(jìn)行基線校正. 與多項(xiàng)式擬合算法相比, 本文所提算法能以固定的階數(shù)進(jìn)行擬合; 與均勻B樣條算法相比, 本文所提算法能根據(jù)不同的拉曼光譜自適應(yīng)選擇內(nèi)節(jié)點(diǎn); 與非均勻B樣條算法和airPLS相比, 本文所提算法具有更好的基線擬合效果.
本文提出了一種簡單靈活且有效的拉曼光譜基線校正算法. 該方法能自適應(yīng)地根據(jù)不同光譜數(shù)據(jù)選取內(nèi)節(jié)點(diǎn), 且通過中值濾波算法優(yōu)化B樣條曲線控制系數(shù), 使非均勻B樣條算法對拉曼光譜數(shù)據(jù)進(jìn)行擬合時(shí)能獲得更好的基線擬合效果. 與多項(xiàng)式擬合, 均勻B樣條和airPLS算法相比, 該算法具有更好的基線校正效果, 且不存在欠擬合和過擬合等缺陷. 此外, 該算法具有廣泛的適用性和靈活性, 針對不同的基線漂移情況, 可以有效地進(jìn)行基線擬合. 因此, 本文所提算法能為拉曼光譜信號的進(jìn)一步分析提供更準(zhǔn)確可靠的信息.