優(yōu)化Savitzky-Golay濾波器的參數(shù)及其在傅里葉變換紅外氣體光譜數(shù)據(jù)平滑預處理中的應用

趙安新，湯曉君，張鐘華, 3，劉君華

1.西安科技大學，陜西 西安 710054 2.西安交通大學電力設備電氣絕緣國家重點實驗室，陜西 西安 710049 3.中國計量科學研究院，北京 100013

優(yōu)化Savitzky-Golay濾波器的參數(shù)及其在傅里葉變換紅外氣體光譜數(shù)據(jù)平滑預處理中的應用

趙安新1, 2，湯曉君2，張鐘華2, 3，劉君華2

1.西安科技大學，陜西 西安 710054 2.西安交通大學電力設備電氣絕緣國家重點實驗室，陜西 西安 710049 3.中國計量科學研究院，北京 100013

對傅里葉變換中紅外光譜數(shù)據(jù)的平滑預處理中，通常采用Savitzky-Golay濾波器的方法進行光譜數(shù)據(jù)的平滑預處理，然而Savitzky-Golay濾波器的多項式擬合階次和窗寬等參數(shù)的合理選則始終是一個難題，并無統(tǒng)一的選擇依據(jù)，通常在一定數(shù)值范圍內(nèi)，采用多組數(shù)據(jù)進行遍歷嘗試，最終選擇一組相對較優(yōu)的數(shù)據(jù)作為Savitzky-Golay濾波器的多項式擬合階次和窗寬參數(shù)。文中探索了Savitzky-Golay濾波器的多項式擬合階次和窗寬等參數(shù)的優(yōu)化選取這一問題，并對其主要頻率指標參數(shù)與多項式擬合階次和窗寬等參數(shù)進行了定性定量分析，得出了截至頻率、阻帶起始點頻率、第一旁瓣峰值頻率及第一旁瓣峰值幅度與窗寬和階次之間具體的計算方程表達式。隨后，根據(jù)采集的中紅外氣體組分的光譜數(shù)據(jù)特征，依據(jù)上述計算方程式，優(yōu)化計算選取多項式擬合階次和窗寬分別為8和11時，其Savitzky-Golay濾波器的中紅外氣體組分的光譜數(shù)據(jù)平滑效果最優(yōu)。最后通過實際采集的0.1%, 0.2%, 0.5%, 1%, 2%, 5%濃度的CH4光譜數(shù)據(jù)進行平滑預處理，在次吸收峰區(qū)域，原始光譜的折算吸光率相對最大誤差和最小誤差分別為17.230 5%和0.243 0%，平滑處理后的光譜的折算吸光率相對最大誤差和最小誤差分別為0.088 0%和0.020 6%。可見經(jīng)過Savitzky-Golay濾波器進行所探索的光譜數(shù)據(jù)預處理之后其相對誤差基本穩(wěn)定，并且相對較低，為后期光譜數(shù)據(jù)的準確定性和定量分析奠定了基礎。

Savitzky-Golay濾波器； 光譜數(shù)據(jù)預處理； 折算吸光率； 傅里葉變換紅外光譜

引 言

在光譜分析數(shù)據(jù)預處理的方法中通常采用Savitzky-Golay濾波器(S-G)方法來進行原始數(shù)據(jù)的平滑與去噪[1-3]。Savitzky-Golay濾波器(S-G)方法[4]最早由Savitzky和Golay于1964年提出，其基本原理是基于最小二乘擬合，用擬合值代替原始數(shù)值，以達到去除高頻噪聲點，平滑原數(shù)據(jù)序列的作用。在S-G平滑濾波過程中，濾波器窗寬及多項式擬合階次[5]決定了S-G濾波器平滑和去噪效果。多項式擬合階次越高，能保持的中心矩越高，即保持的細節(jié)信息越多。而濾波器窗寬越大，平滑效果越好，對噪聲的衰減越大。對多項式擬合階次與濾波器窗寬的選擇，即是對保持更多信號有效細節(jié)信息與濾除更多噪聲信號的綜合考慮，因此構造一個合適的量化指標來評價各參數(shù)下濾波器平滑的效果顯得至關重要。Browne等[6]提出了一種多尺度多項式濾波器，Larivee等[7]選擇將平滑后信號的信息熵作為評價標準，Phillip Barak[8]將擬合信號的殘差平方和作為平滑效果的評判標準，Jakubowska等[9]將Barak提出的自適應階次S-G濾波算法用于平滑復雜的仿真Feldberg典型伏安信號。Vivo-Truyols等[10]認為當擬合殘差的自相關與儀器噪聲的自相關最接近時，平滑效果最佳，此時濾波器所選窗寬為最優(yōu)值。因此發(fā)展一種簡單易用、適合大眾的方法很重要。之前雖已有人[11-13]對S-G濾波器頻率特性進行了較深入的研究，然而并未發(fā)現(xiàn)有人得出濾波器窗寬、多項式擬合階次等參數(shù)與頻率特性的具體量化關系，因此無法指導用戶對參數(shù)進行選擇，這是S-G平滑濾波器一直以來存在的缺陷。本研究從頻域出發(fā)，探索了不同參數(shù)下的S-G平滑濾波器的頻率特性，據(jù)此擬合出濾波器的主要頻率特性指標(截止頻率、阻帶起始點頻率、最大旁瓣峰值頻率、最大旁瓣峰值幅值等指標)與階次及窗寬的關系。以此關系式，可對實際光譜的平滑濾波需求定量選擇出S-G平滑濾波器的階次與窗寬。

1 實驗部分

1.1 材料及儀器

實驗采用CH4氣體傅里葉紅外光譜作為分析處理目標，設定樣本氣CH4的濃度分別為0.1%，0.2%，0.5%，1%，2%和5%。

采用Bruke公司的alpha型傅里葉變換紅外光譜儀，探測器為氘代硫酸三肽DTGS檢測器，掃描范圍為500～4 000 cm-1，光譜分辨率為1 cm-1，譜線值為吸光度光譜。

1.2 折算吸光率

由朗伯-比耳定律(Lambert-Beer)[14-15]可知，氣體在任一波數(shù)ν處的吸光度A(ν)為

(1)

式(1)中：T(ν)為波數(shù)ν處的透射率光譜值；I0為紅外光透過背景(通常是空光路)的光強，可近似等于入射光強；I為紅外光透過氣體樣品的光強；a(ν)為氣體在波數(shù)ν處的吸光度系數(shù)，在不同波數(shù)處的數(shù)值不相同，但對于特定波數(shù)，可認為不變；b為光程(樣品厚度)；c為氣體樣品濃度。

由朗伯-比耳(Lambert-Beer)定律可知任一波長光的吸收強度(吸光度)與樣品中各組分的濃度(氣體濃度)成正比，與光程差(氣體池長度)成正比。

對于給定的光譜儀，若氣體池不變，則認為光程b是恒定不變的，因此可用折算吸光率δ(ν)[16]來等價吸光率a(ν)與光程差b，即

δ(ν)=a(ν)×b

(2)

對于單組份樣品，式(1)可轉(zhuǎn)化為

(3)

進而可知單組份氣體折算吸光率δ(ν)的計算式(4)。

(4)

理論上在特定波數(shù)ν處，氣體濃度在一定范圍內(nèi)時，吸光率a(ν)恒定不變，對應的折算吸光率δ(ν)也不變。然而由于干涉信號在采樣過程中存在信號截斷，導致光譜存在能量泄露。另外，氣體對紅外光也會存在一定程度的散射，再加上光譜噪聲的影響，會導致折算吸光率δ(ν)在光譜存在吸收的區(qū)域隨著氣體濃度的增加而呈現(xiàn)下降趨勢。

2 結(jié)果與討論

2.1 主要Savitzky-Golay濾波器頻率參數(shù)分析

為探索Savitzky-Golay濾波器多項式擬合階次與窗寬之間的時域和頻域特性關系，使用計算機模擬軟件，分別模擬Savitzky-Golay濾波器在多項式擬合階次d=2, 4, 6, 8, 10, 12, 14，窗寬范圍N為(d+3, 121)區(qū)間內(nèi)奇數(shù)(odd)的特性指標，計算其主要頻率特性指標(截止頻率fc、阻帶起始點頻率fb、第一旁瓣峰值頻率fp及第一旁瓣峰值幅度Ap)。

圖1為Savitzky-Golay濾波器主要頻率特性指標截止頻率fc、阻帶起始點頻率fb、第一旁瓣峰值頻率fp及第一旁瓣峰值幅度Ap與窗寬之間的關系。四種頻率特性指標參數(shù)與窗寬基本都呈現(xiàn)指數(shù)分布形式，在一定擬合階次下，各項指標參數(shù)隨著窗寬的增大按照指數(shù)分布的形式逐漸減小，在窗寬10～40區(qū)間為其轉(zhuǎn)折點，并且隨著擬合階次的增加，其轉(zhuǎn)折點逐漸右移且變化越來越平緩； 在窗寬大于40的區(qū)域，各項頻率指標參數(shù)變化越來越平緩； 在窗寬大于100的區(qū)域，各項頻率指標參數(shù)基本不再隨著窗寬改變而變化。此項特征對于根據(jù)特定的平滑和濾波對象，選擇多項式擬合階次和窗寬提供了依據(jù)，如果已知平滑對象的各項頻率指標參數(shù)，可以根據(jù)此特征選定Savitzky-Golay濾波器的參數(shù)。同時，如果知道Savitzky-Golay濾波器的參數(shù)，可以確定平滑對象的各項頻率指標參數(shù)。根據(jù)式(5)—式(8)計算的數(shù)據(jù)，分別擬合了Savitzky-Golay濾波器主要頻率特性指標截止頻率fc、阻帶起始點頻率fb、第一旁瓣峰值頻率fp及第一旁瓣峰值幅度Ap與窗寬和階次之間的關系方程式。

圖1 Savitzky-Golay濾波器主要特性指標與窗寬N在不同多項式擬合階次下的變換關系

(5)

(6)

(7)

(8)

2.2 傅里葉變換紅外光譜數(shù)據(jù)預處理

圖2為濃度分別為0.1%，0.2%，0.5%，1%，2%，5%濃度CH4氣體的傅里葉變化紅外吸收光譜，該光譜采用分段比基線校正方法對其進行了數(shù)據(jù)基線校正的預處理[16]，光譜經(jīng)過基線校正預處理之后，其基線的漂移和部分畸變得到了一定程度的抑制。其在中紅外波段主要存在兩個主要的吸收區(qū)域，即是主吸收峰區(qū)域(2 850～3 200 cm-1)和次吸收峰區(qū)域(1 200～1 400 cm-1)。從主、次吸收峰處峰高隨濃度的變化趨勢可以看出，在此濃度范圍內(nèi)吸光度和濃度基本呈正比狀態(tài)，初步驗證了1.2節(jié)朗伯-比耳定律(Lambert-Beer)在推導折算吸光率時得出的結(jié)論。在次級吸收峰區(qū)域其正比例關系優(yōu)于主吸收峰，次級吸收峰的諧波較少，主吸收峰區(qū)域的諧波較多，相互之間存在一定的干擾，因此在2.3節(jié)的計算和推導過程中采用次級吸收峰區(qū)域進行計算。

圖2 經(jīng)過基線校正后濃度0.1%，0.2%，0.5%，1%，2%，5%的甲烷CH4紅外吸收光譜

Fig.2 Normal absorbance infrared spectra of concentration 0.1%, 0.2%, 0.5%, 1%, 2% and 5% for CH4in the region of 800～3 850 cm-1after baseline correction

2.3 光譜平滑與折算吸光度分析

以式(4)計算吸光度，0.1%, 0.2%, 0.5%, 1%, 2%和5%濃度的CH4氣體在原始光譜、Savitzky-Golay平滑和去噪后光譜次吸收峰區(qū)域的折算吸光率如表1所示，表1列出了不同濃度(0.1%, 0.2%, 0.5%, 1%, 2%和5%)甲烷(CH4)目標氣體在次吸收峰區(qū)域內(nèi)的折算吸光率及其與相應標準濃度之間的相對誤差。從標準光譜的折算吸光度計算結(jié)果來看，初步驗證了研究中所的得出的結(jié)論，即不同濃度下折算吸光率相差不大，但隨著濃度的增加折算吸光率逐漸減小。原始光譜濃度越小與標準光譜之間的誤差越大。利用Savitzky-Golay濾波器進行平滑和去噪后，其折算吸光度與標準光譜之間的誤差相對較小，基本與濃度沒有特定的關系，誤差維持在一定范圍內(nèi)，初步驗證了本方法的有效性。圖3為濃度為0.1%的CH4氣體，經(jīng)過Savitzky-Golay平滑和去噪后的傅里葉變化紅外吸收光譜，從圖3觀察，其平滑之后在中紅外波段波形基本無變化，在邊界具有水汽吸收區(qū)域基本抑制了水汽的影響，同時光譜變得較平滑。Savitzky-Golay平滑和去噪后的吸收光譜和標準吸收光譜基本重合，偏差很小。同時其在次級吸收峰區(qū)域的折算吸光度分別為：25.901 7和26.248 9，標準光譜在次吸收峰區(qū)域的折算吸光率為26.272 0，其相對誤差分別為1.409 5%和0.088 0%。說明平滑和濾波之后，其處理效果較好，同時在次級吸收峰其折算吸光率更能反映處理效果。而且次級吸收峰由于在低濃度范圍，其吸光度和濃度的正比例關系優(yōu)于主吸收峰區(qū)域，計算的結(jié)果也驗證了此結(jié)論。為后期的定性和定量計算奠定了基礎。

圖3 經(jīng)過Savitzky-Golay平滑和去噪后濃度0.1%的甲烷CH4紅外吸收光譜

Fig.3 Normal absorbance infrared spectra of concentration 0.1% for CH4in the region of 800～3 850 cm-1after Savitzky-Golay smoothing and denoising

表1 不同濃度甲烷(CH4)氣體在主吸收峰區(qū)域折算吸光率及其相對誤差(標準光譜、原始光譜和平滑處理后光譜)

Table 1 The converted absorbance and relative error in the second peak absorption region for the different concentrations of methane CH4(standard spectra, raw sqectra and smoothing and denoising spectra)

標準光譜原始光譜處理后光譜測量值測量值相對誤差/%測量值測量值0.1%27.690132.461217.230527.68440.02060.2%27.020830.762213.846427.01410.02490.5%27.068427.13420.243027.05530.04831.0%26.272025.90171.409526.24890.08802.0%25.423025.94462.051625.36930.21115.0%18.938419.25091.650218.83170.5634

3 結(jié) 論

針對傅里葉變換紅外光譜數(shù)據(jù)平滑預處理的常用方法Savitzky-Golay濾波器參數(shù)的選擇問題，探索了其主要頻率指標參數(shù)與參數(shù)擬合階次和窗寬之間的變化關系，提出其定量的計算方程式，為其在后續(xù)的平滑數(shù)據(jù)預處理中的應用提供依據(jù)。其次，利用實測的不同濃度的甲烷CH4光譜，采用上述結(jié)論對其進行平滑預處理。最后，利用折算吸光度對原始光譜數(shù)據(jù)和平滑處理后的光譜數(shù)據(jù)計算結(jié)果對比可知，原始光譜數(shù)據(jù)的折算吸光率相對最大誤差和最小誤差分別為17.230 5%和0.243 0%，平滑處理后的光譜的折算吸光率相對最大誤差和最小誤差分別為0.088 0%和0.020 6%。大大降低了光譜折算吸光率的計算誤差，初步驗證了本方法的有效性，同時為后續(xù)的光譜數(shù)據(jù)定性和定量分析奠定了基礎。

[1] Zimmermann B, Kohler A.Appl.Spectrosc., 2013, 67: 892.

[2] Serafinczuk J, Pietrucha J, Schroeder G, et al.Opt.Appl., 2011, 41: 315.

[3] Staggs J.Fire Safety J, 2005, 40: 493.

[4] Gorry P A.Anal.Chem., 1990, 62: 570.

[5] Luo J W, Ying K, Bai J.Signal Process, 2005, 85: 1429.

[6] Browne M, Mayer N, Cutmore T R H.Digit Signal Process, 2007, 17: 69.

[7] Larivee R J, Brown S D.Anal.Chem., 1992, 64: 2057.

[8] Barak P.Anal.Chem., 1995, 67: 2758.

[9] Jakubowska M, Kubiak W W.Anal.Chim.Acta, 2004, 512: 241.

[10] Vivó-Truyols G, Schoenmakers P J.Anal.Chem., 2006, 78: 4598.

[11] Luo J W, Ying K, He P, et al.Digit Signal Process, 2005, 15: 122.

[12] Zuo C, Chen Q, Yu Y, et al.Opt.Express, 2013, 21: 5346.

[13] Browne M, Mayer N, Cutmore T R H.Digit Signal Process, 2007, 17: 69.

[14] Bunaciu A A, Aboul-Enein H Y, Fleschin S.Appl.Spectrosc.Rev., 2011, 46: 251.

[15] Tang X, Zhang L, Wang E, et al.Spectrosc.Spect.Anal., 2012, 32: 2730.

[16] Tang X, Wang J, Zhang L, et al.Spectrosc.Spect.Anal., 2013, 33: 334.

(Received Jan.20, 2015; accepted Apr.18, 2015)

Optimizing Savitzky-Golay Parameters and Its Smoothing Pretreatment for FTIR Gas Spectra

ZHAO An-xin1, 2，TANG Xiao-jun2，ZHANG Zhong-hua2, 3，LIU Jun-hua2

1.Xi’an University of Science and Technology, Xi’an 710054, 2.State Key Laboratory of Electrical Insulation and Power Equipment, Xi’an Jiaotong University，Xi’an 710049, 3.National Institute of Metrology, China, Beijing 100013, China

In the smoothing pretreatment for the quantitative analysis of hydrocarbon mixed gases by Fourier transform infrared analysis (FTIR), the Savitzky-Golay filter is usually used as one of the smoothing preprocessing methods in the Fourier transform infrared spectrum data smoothing pretreatment.However, the parameters of the Savitzky-Golay filter such as the polynomial order and frame size are not easy to decide.There is no one unified choice basis.Users usually adopt multiple sets in the special data set to try, and then select a set of relatively optimal data as the optimizing parameters of the Savitzky-Golay filter.The optimal selection method of the Savitzky-Golay filter parameters was explored, and the concrete calculation equations were deduced according to the relation among the normalized cut-off frequency, the normalized beginning frequency of the stopband, the normalized first side lobe peak frequency of the stopband, the normalized first side lobe peak amplitude with the polynomial order and frame size of the Savitzky-Golay filter parameters.Then when the polynomial order and frame size are set as 8 and 11 respectively according the above conclusion and the characteristics of the actual spectral data, the Savitzky - Golay filter smoothing effect is optimum.Through the acquisition the concentration of 0.1%, 0.2%, 0.5%, 1%, 2%, 5% for the actual CH4spectra, the relative maximum and minimum error of the raw spectra converted absorbance were 17.230 5% and 0.243 0% respectively, and the relative maximum and minimum error of the smooth spectra converted absorbance were 0.088 0% and 0.088 0% respectively in the second absorption peak.The relative error of converted absorbance was basically stable through the Savitzky-Golay filter after the spectral data preprocessing and it was relatively low, so, it laid a foundation for the late spectral data accurate qualitative and quantitative analysis.

Savitzky-Golay filter; Spectra pretreatment; Converted absorbance; Fourier transform infrared spectroscopy

2015-01-20，

2015-04-18

國家重大科學儀器設備開發(fā)專項(2012YQ240127)，國家自然科學基金項目(51277144)，電力設備電氣絕緣國家重點實驗室基金項目(EIPE11307)

趙安新，1981年生，西安科技大學博士研究生 e-mail: zhaoanxin@126.com

O433.4

A

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