激光拉曼光譜解譜和處理的方法研究

薛曉康　李曉宇　丁 卯

(1 上海化學品公共安全工程技術研究中心，上海 200062;2 上海化工研究院檢測中心，上海 200062)

引言

拉曼光譜可以被看作是一項“指紋”技術，因為它可以提供非常豐富的結構信息。因此拉曼光譜可以被用作物質的定性識別。并且拉曼光譜具有制樣簡單，不破壞樣品，在幾乎所有的環境下都可以采集。

由于拉曼光譜具有上述的優點，故在化學品成分分析中被廣泛應用。但是拉曼光譜激光源通常是可見光，所以有易產生噪聲，熒光干擾嚴重的缺點。這些缺點會影響對樣品的定性定量分析，然而現在的硬件技術無法避免這些缺點，所以這時就需要使用數學算法對拉曼光譜圖進行后期的處理以達到過濾噪聲和熒光的目的。

中國專利(CN103217409B)公開了一種拉曼光譜的預處理方法[1]。其使用基于小波變換的自適應閾值去噪聲，采用非對稱最小二乘的基線校正算法去除熒光背景。本文使用基于自適應迭代重加權懲罰最小二乘法的算法進行基線校正，使用基于懲罰最小二乘法的算法進行平滑以及使用連續小波變換進行峰檢測。從而改善了基于非對稱最小二乘法的傳統基線校正方法的兩個缺陷：首先，平滑參數需要優化以便得到最優結果；其次，非對稱參數對于所有的基線數據點都是一成不變的。因這樣基線可能會出現負值部分[2]。

1　實驗部分

1.1　儀器試劑

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

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

1.2　實驗方法

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

2　結果與討論

以化學計量學為基礎，信號處理技術為工具，配合計算機算法的數據處理方法。具體步驟如下：1)對拉曼光譜原始信號進行基于自適應迭代重加權懲罰最小二乘法的基線校正。2)對進行完第一步的拉曼光譜信號進行基于懲罰最小二乘法的平滑。3)對進行完第一步和第二步的信號進行基于連續小波變換的峰檢測。

2.1　基線校正[2]

對拉曼光譜原始信號進行基于自適應迭代重加權懲罰最小二乘法算法的具體步驟如圖1所示。

圖1　自適應迭代重加權懲罰最小二乘法基線校正結構圖Figure 1　Structure of baseline correction by adaptive iteratively reweighted penalized least squares (airPLS).

自適應迭代重加權懲罰最小二乘法的表達式(1)為：

(1)

式中，Q為原始基線與擬合后的基線保真度與粗糙度間的平衡。t為迭代次數。w為權重向量，通過自適應迭代方法得到。x為原始信號向量，z為擬合向量，x與z的長度記為m，λ為粗糙度系數。

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

(2)

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

迭代會在達到最大迭代次數或滿足條件(3)式時結束(圖2-3)：

|dt|<0.001×|x|

(3)

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

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

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

2.2　基線平滑[3]

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

1)將公式(1)中的加權系數去除，即得到峰平滑的數學表達式(4)：

(I+λD′D)z=y

(4)

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

(5)

(6)

(7)

H矩陣的列可以通過平滑其所對應的單位矩陣找到。

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

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

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

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

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

2.3　峰檢測[4]

對進行基線平滑的拉曼光譜信號進行基于連續小波變換峰檢測算法的具體步驟如下：

1)進行峰檢測的條件有很多，比如信噪比、峰強度閾值、峰形、脊線、極大值、峰寬等。本實驗是使用信噪比和脊線作為峰檢測條件，用連續小波變換作為算法。連續小波變換是對信號函數與經過縮放與平移的小波母函數乘積在整個時間域的積分。其公式如式(8)：

(8)

S(t)是信號，a是縮放系數，b是位移系數。Ψ(t)是小波母函數，Ψa,b(t)即為經過縮放和平移的小波函數。結果C(a,b)是一個小波系數的二維矩陣(2D)。

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

(9)

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

3)當將此方法進行峰探測時，連續小波變換系數在任意縮放系數下都會在峰中心周圍有一個極大值。極大值會在和峰寬匹配時達到最大。當把對連續小波變換系數的縮放倍數作為第三維度放到連續小波變換二維系數圖中時，在峰位置就會出現一條清晰的脊線。所以峰檢測算法此時就包含三個步驟：通過連接極大值來識別脊線；識別出代表峰的脊線以及優化峰參數[5]。

現存的峰檢測方法都無法在不影響假陽性率的情況下同時檢測出強峰和弱峰。本論文中提供的方法可以在背景中通過峰形進行不同尺度間的峰檢測，同時假陽性的頻率并沒有提高。

2.4　峰寬估算

對進行完基線平滑的拉曼光譜信號進行基于連續小波變換峰檢測的算法還可估算出拉曼光譜圖中峰的寬度。

1)此處使用的算法是基于哈爾小波函數的微分運算。根據哈爾小波函數的特點，一組信號的n次導數可以通過應用n次連續小波變換來實現。哈爾小波函數的數學表達式為式(10)[6]：

(10)

2)峰寬評估步驟如下：

①使用在峰檢測中同樣的縮放系數對此哈爾小波進行連續小波變換。二維連續小波變換系數以M×N的矩陣表示。

②然后對此矩陣中所有值取絕對值。

③對于在峰檢測階段檢測到的每一個峰都有兩個參數：峰指數和峰尺寸。二維連續小波變換矩陣中對應峰尺寸的行被用來從峰指數中尋找每個區域的極小值。

④如果極小值不存在，那么峰的起點或終點就是三倍于其峰尺寸的最小值或下一個峰指數。如果極小值存在，那么峰起點或終點就是最近的那個極小值。

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

當需要進行峰檢測光譜的背景非常復雜時，峰寬估算就可以幫助峰檢測算法進行背景估算進而準確地進行峰檢測。當背景不是很復雜，并且峰都很好地分開的時候，峰寬估算就用來進行傳統意義上的峰寬計算了。

3　結論

基于懲罰最小二乘法的光譜平滑具有快速，可以連續控制平滑度并且可以進行交叉驗證得到最客觀的平滑值。改善了基于非對稱最小二乘法的傳統基線校正方法的兩個缺陷：首先，平滑參數需要優化以便得到最優結果；其次，非對稱參數對于所有的基線數據點都是一成不變的，這樣基線可能會出現負值部分。同時，基于連續小波變換的峰檢測算法可以自動地并且同時考慮峰形和峰高對峰進行檢測，最大地降低了峰檢測假陽性的概率。

[1] 張煒,何石軒,杜春雷,等. 一種拉曼光譜預處理方法:中國,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.