薛曉康　李曉宇　丁 卯

(1 上?；瘜W(xué)品公共安全工程技術(shù)研究中心，上海 200062;2 上海化工研究院檢測(cè)中心，上海 200062)

引言

拉曼光譜可以被看作是一項(xiàng)“指紋”技術(shù)，因?yàn)樗梢蕴峁┓浅ＸS富的結(jié)構(gòu)信息。因此拉曼光譜可以被用作物質(zhì)的定性識(shí)別。并且拉曼光譜具有制樣簡(jiǎn)單，不破壞樣品，在幾乎所有的環(huán)境下都可以采集。

由于拉曼光譜具有上述的優(yōu)點(diǎn)，故在化學(xué)品成分分析中被廣泛應(yīng)用。但是拉曼光譜激光源通常是可見(jiàn)光，所以有易產(chǎn)生噪聲，熒光干擾嚴(yán)重的缺點(diǎn)。這些缺點(diǎn)會(huì)影響對(duì)樣品的定性定量分析，然而現(xiàn)在的硬件技術(shù)無(wú)法避免這些缺點(diǎn)，所以這時(shí)就需要使用數(shù)學(xué)算法對(duì)拉曼光譜圖進(jìn)行后期的處理以達(dá)到過(guò)濾噪聲和熒光的目的。

中國(guó)專(zhuān)利(CN103217409B)公開(kāi)了一種拉曼光譜的預(yù)處理方法[1]。其使用基于小波變換的自適應(yīng)閾值去噪聲，采用非對(duì)稱(chēng)最小二乘的基線校正算法去除熒光背景。本文使用基于自適應(yīng)迭代重加權(quán)懲罰最小二乘法的算法進(jìn)行基線校正，使用基于懲罰最小二乘法的算法進(jìn)行平滑以及使用連續(xù)小波變換進(jìn)行峰檢測(cè)。從而改善了基于非對(duì)稱(chēng)最小二乘法的傳統(tǒng)基線校正方法的兩個(gè)缺陷：首先，平滑參數(shù)需要優(yōu)化以便得到最優(yōu)結(jié)果；其次，非對(duì)稱(chēng)參數(shù)對(duì)于所有的基線數(shù)據(jù)點(diǎn)都是一成不變的。因這樣基線可能會(huì)出現(xiàn)負(fù)值部分[2]。

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

1.1　儀器試劑

激光拉曼光譜儀(美國(guó)必達(dá)泰克公司)；數(shù)據(jù)采集軟件：BWspec3.27；激發(fā)波長(zhǎng)785 nm，光譜掃描范圍175～3 200 cm-1，激發(fā)功率0～315 mW，分辨率5 cm-1，4 mL石英比色皿。

化學(xué)試劑和樣品1-苯基-3-甲基-5-吡唑啉酮(CAS 89-25-8)均為分析純。

1.2　實(shí)驗(yàn)方法

利用數(shù)據(jù)采集軟件BWspec3.27，設(shè)置積分時(shí)間36 000 ms，采集3次取平均值，激光功率90%，采集樣品1-苯基-3-甲基-5-吡唑啉酮(CAS 89-25-8)原始拉曼光譜數(shù)據(jù)。

2　結(jié)果與討論

以化學(xué)計(jì)量學(xué)為基礎(chǔ)，信號(hào)處理技術(shù)為工具，配合計(jì)算機(jī)算法的數(shù)據(jù)處理方法。具體步驟如下：1)對(duì)拉曼光譜原始信號(hào)進(jìn)行基于自適應(yīng)迭代重加權(quán)懲罰最小二乘法的基線校正。2)對(duì)進(jìn)行完第一步的拉曼光譜信號(hào)進(jìn)行基于懲罰最小二乘法的平滑。3)對(duì)進(jìn)行完第一步和第二步的信號(hào)進(jìn)行基于連續(xù)小波變換的峰檢測(cè)。

2.1　基線校正[2]

對(duì)拉曼光譜原始信號(hào)進(jìn)行基于自適應(yīng)迭代重加權(quán)懲罰最小二乘法算法的具體步驟如圖1所示。

圖1　自適應(yīng)迭代重加權(quán)懲罰最小二乘法基線校正結(jié)構(gòu)圖Figure 1　Structure of baseline correction by adaptive iteratively reweighted penalized least squares (airPLS).

自適應(yīng)迭代重加權(quán)懲罰最小二乘法的表達(dá)式(1)為：

(1)

式中，Q為原始基線與擬合后的基線保真度與粗糙度間的平衡。t為迭代次數(shù)。w為權(quán)重向量，通過(guò)自適應(yīng)迭代方法得到。x為原始信號(hào)向量，z為擬合向量，x與z的長(zhǎng)度記為m，λ為粗糙度系數(shù)。

在迭代開(kāi)始，我們給定w一個(gè)初始值即：w0=1。迭代開(kāi)始之后，在每一個(gè)迭代步驟t，w均可由表達(dá)式(2)得到：

(2)

向量dt包含有在t迭代步，x和zt-1的負(fù)差值。當(dāng)在t-1步迭代時(shí)，如果第i個(gè)數(shù)據(jù)點(diǎn)比zt-1大時(shí)，這個(gè)數(shù)據(jù)點(diǎn)可以被看作是峰上的一點(diǎn)，所以此時(shí)的權(quán)重可以設(shè)置為零以便使其不進(jìn)入下一步迭代。在本發(fā)明中，這種方法可以在權(quán)重向量w中自動(dòng)地逐步排除峰上的點(diǎn)并保留基線上的點(diǎn)。

迭代會(huì)在達(dá)到最大迭代次數(shù)或滿足條件(3)式時(shí)結(jié)束(圖2-3)：

|dt|<0.001×|x|

(3)

圖2　原始拉曼光譜圖Figure 2　Raw Raman spectrum.

圖3　僅通過(guò)airPLS校正的光譜圖Figure 3　Spectrum corrected by airPLS only.

通過(guò)圖2和圖3可以看出，airPLS算法不管對(duì)直線的基線(175～1 682 cm-1)還是彎曲的基線(1 682～3 699 cm-1)，都可以很好地進(jìn)行校正，說(shuō)明airPLS算法的靈活性很高。同時(shí)也可以發(fā)現(xiàn)airPLS算法在校正基線時(shí)完整地保留了那些很小的峰[2](比如：426～677 cm-1)。這也是airPLS算法的強(qiáng)大之處。

2.2　基線平滑[3]

對(duì)進(jìn)行完基線校正的拉曼光譜信號(hào)進(jìn)行基于懲罰最小二乘法平滑算法的具體步驟如圖4-6：

1)將公式(1)中的加權(quán)系數(shù)去除，即得到峰平滑的數(shù)學(xué)表達(dá)式(4)：

(I+λD′D)z=y

(4)

式中I為單位矩陣；D為微分矩陣；z為平滑后光譜的向量；Δz=Dz；y為原始曲線向量；λ為平滑度。

(5)

(6)

(7)

H矩陣的列可以通過(guò)平滑其所對(duì)應(yīng)的單位矩陣找到。

圖4　僅通過(guò)懲罰最小二乘法平滑的光譜圖Figure 4　Spectrum smoothed by penalized least squares only.

圖5　既通過(guò)懲罰最小二乘法平滑又通過(guò)airPLS校正的光譜圖Figure 5　Spectrum corrected by airPLS and smoothed by penalized least squares.

圖6　通過(guò)Savitzky-Golay平滑的光譜圖(多項(xiàng)式級(jí)數(shù)：3；SG窗口尺寸：15)Figure 6　Spectrum smoothed by Savitzky-Golay(polynomial order:3; SG window size:15).

峰的基線校正和基線平滑順序可以互換，互換處理順序不會(huì)影響處理結(jié)果。

雖然Savitzky-Golay平滑更加知名，但是基于最小二乘法的平滑更加快速和靈活。將此平滑整合到現(xiàn)代軟件中后，將會(huì)在速度、靈活性和交叉驗(yàn)證方面得到極大的提升并且此基本算法在Matlab中很容易編輯。

2.3　峰檢測(cè)[4]

對(duì)進(jìn)行基線平滑的拉曼光譜信號(hào)進(jìn)行基于連續(xù)小波變換峰檢測(cè)算法的具體步驟如下：

1)進(jìn)行峰檢測(cè)的條件有很多，比如信噪比、峰強(qiáng)度閾值、峰形、脊線、極大值、峰寬等。本實(shí)驗(yàn)是使用信噪比和脊線作為峰檢測(cè)條件，用連續(xù)小波變換作為算法。連續(xù)小波變換是對(duì)信號(hào)函數(shù)與經(jīng)過(guò)縮放與平移的小波母函數(shù)乘積在整個(gè)時(shí)間域的積分。其公式如式(8)：

(8)

S(t)是信號(hào)，a是縮放系數(shù)，b是位移系數(shù)。Ψ(t)是小波母函數(shù)，Ψa,b(t)即為經(jīng)過(guò)縮放和平移的小波函數(shù)。結(jié)果C(a,b)是一個(gè)小波系數(shù)的二維矩陣(2D)。

2)由于小波系數(shù)反映了信號(hào)s和Ψa,b(t)間的相似程度，所以小波母函數(shù)的需要具有拉曼光譜峰最基本的特點(diǎn)。本文選擇了“墨西哥帽”函數(shù)作為小波母函數(shù)(圖7)。數(shù)學(xué)表達(dá)為式(9)：

(9)

圖7　“墨西哥帽”函數(shù)示意圖Figure 7　Schema of “Mexican hat” function.

3)當(dāng)將此方法進(jìn)行峰探測(cè)時(shí)，連續(xù)小波變換系數(shù)在任意縮放系數(shù)下都會(huì)在峰中心周?chē)幸粋€(gè)極大值。極大值會(huì)在和峰寬匹配時(shí)達(dá)到最大。當(dāng)把對(duì)連續(xù)小波變換系數(shù)的縮放倍數(shù)作為第三維度放到連續(xù)小波變換二維系數(shù)圖中時(shí)，在峰位置就會(huì)出現(xiàn)一條清晰的脊線。所以峰檢測(cè)算法此時(shí)就包含三個(gè)步驟：通過(guò)連接極大值來(lái)識(shí)別脊線；識(shí)別出代表峰的脊線以及優(yōu)化峰參數(shù)[5]。

現(xiàn)存的峰檢測(cè)方法都無(wú)法在不影響假陽(yáng)性率的情況下同時(shí)檢測(cè)出強(qiáng)峰和弱峰。本論文中提供的方法可以在背景中通過(guò)峰形進(jìn)行不同尺度間的峰檢測(cè)，同時(shí)假陽(yáng)性的頻率并沒(méi)有提高。

2.4　峰寬估算

對(duì)進(jìn)行完基線平滑的拉曼光譜信號(hào)進(jìn)行基于連續(xù)小波變換峰檢測(cè)的算法還可估算出拉曼光譜圖中峰的寬度。

1)此處使用的算法是基于哈爾小波函數(shù)的微分運(yùn)算。根據(jù)哈爾小波函數(shù)的特點(diǎn)，一組信號(hào)的n次導(dǎo)數(shù)可以通過(guò)應(yīng)用n次連續(xù)小波變換來(lái)實(shí)現(xiàn)。哈爾小波函數(shù)的數(shù)學(xué)表達(dá)式為式(10)[6]：

(10)

2)峰寬評(píng)估步驟如下：

①使用在峰檢測(cè)中同樣的縮放系數(shù)對(duì)此哈爾小波進(jìn)行連續(xù)小波變換。二維連續(xù)小波變換系數(shù)以M×N的矩陣表示。

②然后對(duì)此矩陣中所有值取絕對(duì)值。

③對(duì)于在峰檢測(cè)階段檢測(cè)到的每一個(gè)峰都有兩個(gè)參數(shù)：峰指數(shù)和峰尺寸。二維連續(xù)小波變換矩陣中對(duì)應(yīng)峰尺寸的行被用來(lái)從峰指數(shù)中尋找每個(gè)區(qū)域的極小值。

④如果極小值不存在，那么峰的起點(diǎn)或終點(diǎn)就是三倍于其峰尺寸的最小值或下一個(gè)峰指數(shù)。如果極小值存在，那么峰起點(diǎn)或終點(diǎn)就是最近的那個(gè)極小值。

⑤重復(fù)步驟③-④，直到得到所有的峰寬。

當(dāng)需要進(jìn)行峰檢測(cè)光譜的背景非常復(fù)雜時(shí)，峰寬估算就可以幫助峰檢測(cè)算法進(jìn)行背景估算進(jìn)而準(zhǔn)確地進(jìn)行峰檢測(cè)。當(dāng)背景不是很復(fù)雜，并且峰都很好地分開(kāi)的時(shí)候，峰寬估算就用來(lái)進(jìn)行傳統(tǒng)意義上的峰寬計(jì)算了。

3　結(jié)論

基于懲罰最小二乘法的光譜平滑具有快速，可以連續(xù)控制平滑度并且可以進(jìn)行交叉驗(yàn)證得到最客觀的平滑值。改善了基于非對(duì)稱(chēng)最小二乘法的傳統(tǒng)基線校正方法的兩個(gè)缺陷：首先，平滑參數(shù)需要優(yōu)化以便得到最優(yōu)結(jié)果；其次，非對(duì)稱(chēng)參數(shù)對(duì)于所有的基線數(shù)據(jù)點(diǎn)都是一成不變的，這樣基線可能會(huì)出現(xiàn)負(fù)值部分。同時(shí)，基于連續(xù)小波變換的峰檢測(cè)算法可以自動(dòng)地并且同時(shí)考慮峰形和峰高對(duì)峰進(jìn)行檢測(cè)，最大地降低了峰檢測(cè)假陽(yáng)性的概率。

[1] 張煒,何石軒,杜春雷,等. 一種拉曼光譜預(yù)處理方法:中國(guó),ZL201310094703.0[P].2013-03-22.

[2] ZHANG Z M, CHEN S, LIANG Y Z. Baseline correction using adaptive iteratively reweighted penalized least squares[J].Analyst, 2010, 135: 1138-1146.

[3] EILERS P H C. A Perfect Smoother[J].AnalyticalChemistry, 2003, 75: 3631-3636.

[4] DU P, KIBBE W A, LIN S M. Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching[J].Bioinformatics, 2006, 22: 2059-2065.

[5] ZHANG Z M, CHEN X Q, LU H M, et al. Mixture analysis using reverse searching and non-negative least squares[J].ChemometricsandIntelligentLaboratorySystems, 2014, 137: 10-20.

[6] ZHANG Z M, CHEN S, LIANG Y Z, et al. An intelligent background-correction algorithm for highly flourescent samples in Raman spectroscopy[J].JournalofRamanSpectroscopy, 2010, 41: 659-669.