肖 鋒，翟陽陽，武寶明，王雅圣，魏 高，孟國營，翟 紅，李 光

(1.浙江大學(xué)工業(yè)控制技術(shù)國家重點(diǎn)實(shí)驗(yàn)室智能系統(tǒng)與控制研究所，浙江杭州310027;2.首都醫(yī)科大學(xué)，北京100069;

3

.汾西礦業(yè)(集團(tuán))有限責(zé)任公司，山西介休032000)

1 引言

煤的揮發(fā)分是煤樣在規(guī)定條件下隔絕空氣加熱，并進(jìn)行水分校正后的質(zhì)量損失［1］。去掉揮發(fā)分的殘?jiān)薪乖]發(fā)分不是指煤中原有的揮發(fā)性物質(zhì)，而是指煤在嚴(yán)格規(guī)定條件下，加熱分解的產(chǎn)物。煤的揮發(fā)分主要是由水分、碳?xì)浠衔锖吞細(xì)溲趸衔锝M成，物理吸附水和礦物質(zhì)分解生成的二氧化碳不屬于揮發(fā)分。揮發(fā)分是對煤進(jìn)行分類的主要指標(biāo)，根據(jù)煤的揮發(fā)分可以大致判斷其變質(zhì)程度、加工利用性質(zhì)和熱值的高低。實(shí)驗(yàn)室化學(xué)檢測揮發(fā)分的方法規(guī)范性很強(qiáng)，實(shí)驗(yàn)條件需要嚴(yán)格控制，因此對檢測人員的要求非常高，同時(shí)檢測時(shí)間比較長。

將一束不同波長的近紅外光照射到樣本表面，有一部分光射入樣本，在樣本顆粒和內(nèi)部經(jīng)過多次反射、折射、衍射和吸收后返回樣本表面，這樣的反射稱為漫反射。被吸收的光轉(zhuǎn)化為分子能量，使近紅外光譜能表征物質(zhì)性質(zhì)。近紅外光譜的吸收區(qū)域主要為含氫基團(tuán)的合頻以及倍頻吸收，而煤的揮發(fā)分均與含氫基團(tuán)相關(guān)，近紅外光譜檢測煤的揮發(fā)分理論上是可行的。

近紅外光譜技術(shù)最早20世紀(jì)60年代被用來檢測谷物的水分，其快速無損方便的檢測特點(diǎn)使它得到快速的發(fā)展，除了農(nóng)作物，紡織、化工、制藥和煙草等行業(yè)，都研究近紅外光譜的檢測方法。M.P.Fuller和 P.R.Griffiths［3］利用傅里葉近紅外漫反射光譜檢測粉末物質(zhì)的性質(zhì)，F(xiàn)ysh和Swinkels［4］定性分析了煤的近紅外漫反射光譜特性。日本人Mikio Kaihara［5］研究了煤的近紅外光譜與煤的化學(xué)性質(zhì)的相關(guān)性，利用偏最小二乘方法對光譜進(jìn)行建模分析，模型精度較差。J.M.Andres和M.T.Bona分析了煤的近紅外［6－7］和中紅外［8］的光譜特性，證實(shí)近紅外光譜方法檢測煤有機(jī)物成分能取得較好的效果。Hongfei Cheng［9］等人分析了高嶺土、煤系高嶺土的中紅外和近紅外光譜差異，對建模過程中光譜的波段選取有一定的指導(dǎo)意義。

本文利用近紅外漫反射光譜對煤的揮發(fā)分進(jìn)行定量分析?；貧w建模前對光譜進(jìn)行波段篩選和數(shù)學(xué)預(yù)處理，提高光譜信噪比。最后對肥煤單獨(dú)建立回歸模型，進(jìn)行分析。

2 方法

2.1 樣本

本實(shí)驗(yàn)中共挑選139個(gè)煤樣，覆蓋瘦煤、肥煤、焦煤、氣煤和貧瘦煤等多個(gè)不同的煤種。煤樣提供方根據(jù)國家標(biāo)準(zhǔn)——煤的工業(yè)分析方法［2］測定了全部樣本的揮發(fā)分。樣本揮發(fā)分的詳細(xì)參數(shù)如表1所示。同時(shí)，煤樣送到實(shí)驗(yàn)室進(jìn)行近紅外光譜采集之前，均被粉碎到3 mm以下。

2.2 近紅外漫反射光譜采集

實(shí)驗(yàn)采用德國Bruker公司MATRIX－I型工業(yè)現(xiàn)場級Fourier變換近紅外光譜儀采集煤樣光譜。使用Bruker公司提供的OPUS 6.0軟件進(jìn)行光譜儀的控制和近紅外光譜的采集。處理光譜數(shù)據(jù)的預(yù)處理操作和回歸模型建立均使用MATLAB R2007b。

采集光譜時(shí)，為了使裝入樣品杯的煤更好地代表該煤樣的性質(zhì)，需要對煤樣縮分?？s分是指煤樣制備時(shí)將試樣分成具有代表性的幾個(gè)部分，其中一份或多份用于保留。在此，利用堆錐四分法對煤樣縮分［1］，即將煤樣均勻堆成圓錐體，然后壓成厚度均勻的圓餅狀，等分成四個(gè)扇形，取相對的兩個(gè)扇形部分裝入樣品杯作為采集光譜的煤樣。對每一個(gè)煤樣重復(fù)縮分和裝樣15次取平均光譜作為該煤樣的最終光譜。光譜采集過程中，為了進(jìn)一步減小環(huán)境造成的誤差，實(shí)驗(yàn)室溫度保持在(21±1)℃，相對濕度保持在55%±5%。

煤樣光譜采集波段選擇10000～4000 cm－1，分辨率8cm－1，64次疊加平均。光譜獲取模式選為吸收模式(absorbance)。

2.3 數(shù)據(jù)預(yù)處理

實(shí)驗(yàn)中采集的近紅外漫反射光譜包含大量環(huán)境噪聲，也有物理性質(zhì)的影響。在利用光譜對煤的揮發(fā)分進(jìn)行回歸建模時(shí)，利用光譜預(yù)處理方法可以降低噪聲，以及煤樣顆粒度不勻帶來的影響，提高光譜的信噪比。在本研究中，采用了波段選取方法、多元散射校正平滑、一階求導(dǎo)、二階求導(dǎo)和無關(guān)變量消除的光譜預(yù)處理方法。

波段選擇:實(shí)驗(yàn)中光譜采集選擇10000～4000 cm－1波段。如果找出光譜中與揮發(fā)分顯著相關(guān)的波段進(jìn)行建模，可以提高建模精度。實(shí)驗(yàn)中首先利用全波段進(jìn)行建模，然后根據(jù)光譜特征定性選擇存在光譜吸收峰的波段，最后根據(jù)光譜各波段與揮發(fā)分的相關(guān)性進(jìn)行波段選擇。

多元散射校正:近紅外漫反射吸收光譜容易受到待測樣品的顆粒度等影響，造成基線偏移。多元散射校正［12］可以有效消除散射影響導(dǎo)致的基線漂移，提高光譜的信噪比。多元散射校正的方法通過建立一個(gè)待測樣品的“理想光譜”，以該光譜為標(biāo)準(zhǔn)對其他樣品的光譜進(jìn)行修正。由于“理想光譜”無法得到，實(shí)際運(yùn)用中選取所有光譜的平均光譜作為“理想光譜”。

表1 煤樣的揮發(fā)分參數(shù)Tab.1 volatile matter parameters of coal

平滑和求導(dǎo)算法:平滑和求導(dǎo)均是常用的預(yù)處理方法，平滑可以用來降低光譜噪聲，求導(dǎo)能提高光譜靈敏度。實(shí)驗(yàn)中選取Savitzky-Golay平滑和Savitzky-Golay求導(dǎo)的方法對光譜進(jìn)行預(yù)處理。

無關(guān)變量消除:也是一種常用波段選擇方法［11］。無關(guān)變量消除基于偏最小二乘回歸系數(shù)，通過光譜中混入噪聲與原光譜一起作為因變量建模，比較噪聲的系數(shù)和光譜的系數(shù)，認(rèn)為光譜系數(shù)小于噪聲系數(shù)的波段是于回歸建模無關(guān)的變量，予以消除。

3 實(shí)驗(yàn)結(jié)果

3.1 煤近紅外漫反射光譜

不同煤種采集的光譜會有一些細(xì)微的差別，如圖1所示，選取樣本量較多的肥煤、焦煤、瘦煤和貧瘦煤。瘦煤的近紅外光譜與其他煤種的斜率相差較大，且各自譜峰的幅度和位置均有差別。從圖1也可以發(fā)現(xiàn)光譜1350～2400 nm(7406～4166 cm－1)波段譜峰變化顯著。

圖1 4個(gè)不同煤樣光譜圖Fig.1 spectra of 4 different kinds of coal

3.2 模型建立

研究中利用的是偏最小二乘算法［10］建立煤樣光譜和揮發(fā)分的回歸模型。模型好壞的標(biāo)準(zhǔn)選用內(nèi)部交叉檢驗(yàn)參數(shù)，交叉檢驗(yàn)均方差(root mean square error of cross validation，RMSECV)和外部檢驗(yàn)參數(shù)，預(yù)測均方差(root mean square error of prediction，RMSEP)。

首先利用全部波段(即10000～4000 cm－1波段)進(jìn)行建模分析，然后選取譜峰特征明顯的7406～4166 cm－1波段建模，最后利用相關(guān)性定量進(jìn)行選擇波段，求取光譜各個(gè)波段與煤揮發(fā)分的相關(guān)性，如圖2所示。同時(shí)對3種選取波段分別進(jìn)行多元散射校正。表2所示為不同預(yù)處理得到的結(jié)果。

圖2 光譜與煤揮發(fā)分相關(guān)性Fig.2 correlation between the spectra and volatile matter of coal

表2 選取波段和多元散射校正預(yù)處理的建模結(jié)果Tab.2 model results of different spectral ranges and MSC pre-treatment

從表2的結(jié)果可以看出選取波段之后RMSECV和RMSEP均明顯變小，模型精度和穩(wěn)定性得到改善。選擇譜峰明顯波段進(jìn)行建模可以得到較好的結(jié)果。各個(gè)波段光譜與定標(biāo)值的相關(guān)性都很低。利用相關(guān)性進(jìn)行波段選擇時(shí)不同的相關(guān)性閾值可以得到不同的結(jié)果，選擇0.11時(shí)結(jié)果最好，此時(shí)對應(yīng)7251～4000 cm－1波段，與根據(jù)光譜定性選擇波段的方法得到結(jié)果類似。

多元散射校正的預(yù)處理方法沒有使結(jié)果變好，選取波段之后進(jìn)行多元散射校正的結(jié)果更差。Bruker近紅外光譜儀的近紅外光從樣品杯底部照射，取其漫反射光譜。樣品表面比較平整，顆粒度不勻造成的散射影響小，因此使用多元散射校正的方法沒有使結(jié)果變好。

選取7406～4166 cm－1波段，然后 Savitzky-Golay平滑和求導(dǎo)，得到結(jié)果如表3所示。

表3 Savitzky-Golay平滑和求導(dǎo)預(yù)處理建模結(jié)果Tab.3 model results of Savitzky-Golay smoothing and Savitzky-Golay derivative

平滑的方法可以用來消除噪聲，有效信息也會被削弱。對比光譜經(jīng)過平滑和未平滑的結(jié)果，可以發(fā)現(xiàn)結(jié)果接近。平滑沒有明顯優(yōu)化建模效果。求導(dǎo)之后模型內(nèi)部交叉驗(yàn)證結(jié)果變差，預(yù)測能力降低，尤其二階求導(dǎo)之后。求導(dǎo)算法放大了噪聲，減低信噪比，模型的穩(wěn)定性和預(yù)測精度沒有提高。

光譜采集和定標(biāo)值測量過程中不可避免出現(xiàn)一些誤差，誤差的存在會打破光譜和定標(biāo)值的對應(yīng)關(guān)系，建立模型時(shí)會對模型造成比較大的影響，這樣的樣本點(diǎn)稱為異常點(diǎn)，在建模過程中應(yīng)該剔除。實(shí)驗(yàn)中選用學(xué)生氏殘差(Studentized Residual)剔除異常點(diǎn)。學(xué)生氏殘差是多元統(tǒng)計(jì)分析中常用的異常點(diǎn)剔除算法。計(jì)算公式如下所示:其中，yi表示樣本i的定標(biāo)值;表示樣本i的預(yù)測值;Se為樣本預(yù)測偏差的標(biāo)準(zhǔn)差。

圖3所示是針對7406～4166 cm－1波段和7406～4166 cm－1波段Savitzky-Golay平滑運(yùn)用學(xué)生氏殘差剔除異常點(diǎn)的結(jié)果，學(xué)生氏殘差閾值取2.5。

圖3 學(xué)生氏殘差剔除異常點(diǎn)結(jié)果(從左到右依次是7406～4166 cm－1波段剔除異常點(diǎn)前，7406～4166cm－1波段剔除異常點(diǎn)后，7406～4166 cm－1波段Savitzky-Golay平滑剔除異常點(diǎn)前，7406～4166 cm－1波段Savitzky-Golay平滑剔除異常點(diǎn)后)Fig.3 results of deleting outliers with semi-studentized residual(From left to right:7406 ～4166 cm－1with outliers，7406 ～4166 cm－1 deleting outliers，7406 ～4166 cm －1and Savitzky-Golay smoothing with outliers，7406 ～4166 cm －1and Savitzky-Golay smoothing deleting outliers)

可以看出，雖然模型預(yù)測結(jié)果(RMSEP)只是稍微改善，但是RMSECV明顯變小，學(xué)生氏殘差剔除異常點(diǎn)可以提高模型的穩(wěn)定性和預(yù)測精度。

對于選取波段7406～4166 cm－1進(jìn)行Savitzky-Golay平滑剔除異常點(diǎn)的情況進(jìn)行無關(guān)變量消除，RMSECV=0.57，RMSEP=0.78。其預(yù)測集的預(yù)測效果如圖4所示。

3.3 分類建模探索

全體煤樣一起建模時(shí)，不同煤種中粘結(jié)度，膠質(zhì)層指數(shù)等物理性質(zhì)上的不同會反應(yīng)在光譜中，這些信息在建模時(shí)會作為噪聲影響模型的結(jié)果，而一般的預(yù)處理方法無法有效的消除這些因素的影響。因此考慮如果對不同種類的煤樣分類建模。相同種類的煤樣物理化學(xué)性質(zhì)接近，建模時(shí)能降低噪聲的干擾。

本次樣品中包含6種煤樣，如表1所示。焦煤和瘦煤揮發(fā)分的標(biāo)準(zhǔn)差較小，覆蓋范圍小，氣煤和1/3焦煤的樣本數(shù)量又太少。綜合考慮揮發(fā)分的分布范圍和樣本個(gè)數(shù)，選用肥煤進(jìn)行分類建模探索。

進(jìn)行回歸建模時(shí)同樣將肥煤分成預(yù)測集和校正集，均勻選取6個(gè)樣本作為預(yù)測集，25個(gè)作為校正集。建模結(jié)果如表4所示。

表4 肥煤回歸建模結(jié)果Tab.4 regression model results of fat coal

由于求導(dǎo)和多元散射校正的方法對全部煤樣的建模效果不好，因此沒有采用。對肥煤分別進(jìn)行全波段，7406 ～4166 cm－1和7406～4166 cm－1平滑的分析，同時(shí)利用學(xué)生氏殘差進(jìn)行了異常點(diǎn)剔除。平滑采用移動窗口寬度7，多項(xiàng)式次數(shù)2的Savitzky-Golay平滑。其中，選取7406～4166 cm－1波段然后運(yùn)用無關(guān)變量消除后，RMSECV和RMSEP為0.40，0.67。對比建模的結(jié)果可以發(fā)現(xiàn)對肥煤的建模精度和穩(wěn)定性明顯優(yōu)于整體煤樣。證明對煤樣先進(jìn)行分類然后回歸建模是可行的。

4 討論和結(jié)論

實(shí)驗(yàn)中，為了最大程度降低光譜采集過程中的噪聲，嚴(yán)格控制了實(shí)驗(yàn)室的溫度和濕度范圍，盡量減小環(huán)境變化對光譜的影響。同時(shí)對煤樣進(jìn)行15次重復(fù)裝樣，降低由于煤樣表面不均勻或者環(huán)境原因產(chǎn)生的隨機(jī)誤差影響。

數(shù)據(jù)分析過程中，選取偏最小二乘算法對光譜進(jìn)行回歸分析。同時(shí)，嘗試不同的預(yù)處理方法，如平滑、求導(dǎo)和波段選取，以及針對顆粒物品的多元散射校正對光譜進(jìn)行處理。分析證實(shí)在對煤樣光譜與揮發(fā)分含量進(jìn)行回歸建模時(shí)，波段選取能改善模型結(jié)果，通過選擇含有效信息更多的波段可以提高光譜的信噪比。進(jìn)行學(xué)生氏殘差剔除異常點(diǎn)之后建模結(jié)果明顯改善。利用無關(guān)變量消除再次挑選與揮發(fā)分相關(guān)性大的波段進(jìn)行建模，進(jìn)一步改善了模型結(jié)果。最后對分類建模的方法進(jìn)行了嘗試，選取肥煤進(jìn)行類似的回歸建模分析，結(jié)果顯示單獨(dú)對肥煤建模的效果優(yōu)于全體煤樣。可以先對煤樣分類然后回歸建模的方法。

揮發(fā)分主要由水分、碳?xì)溲趸锖吞細(xì)浠衔锝M成，而近紅外光譜對含氫基團(tuán)物質(zhì)檢測能力較好。理論上證實(shí)了煤的近紅外光譜能包含揮發(fā)分信息，而近紅外光譜方法用于檢測煤的灰分和硫分時(shí)，因?yàn)榛曳种饕蔁o機(jī)物構(gòu)成，而硫分則包含了有機(jī)硫和無機(jī)硫，導(dǎo)致近紅外光譜對煤灰分和硫分的檢測精度不如對揮發(fā)分的檢測精度。

近紅外漫反射光譜方法檢測煤揮發(fā)分含量，無損快速的檢測特征讓其具有良好的運(yùn)用前景。需要增加樣本數(shù)量以提高模型的穩(wěn)定性和精度，同時(shí)運(yùn)用分類建模的方法進(jìn)一步提高檢測精度。

［1］GB/T 3715－2007．煤質(zhì)及煤分析有關(guān)術(shù)語［S］．北京:中國標(biāo)準(zhǔn)出版社，2007．

［2］GB/T 212－2008，煤的工業(yè)分析方法［S］．北京:中國標(biāo)準(zhǔn)出版社，2008．

［3］M P Fuller，P R Griffiths．Diffuse reflectance measurements by infrared Fourier transform spectrometry［J］．Analytical Chemistry，1978，50(13):1906 －1910．

［4］A Fysh，A J Swinkels．Near-infrared diffuse reflectance spectroscopy of coal［J］．Applied Spectroscopy，1985，39(2):354－357．

［5］Mikio Kaihara，Toshibiro Sato，et al．Application of near infrared spectroscopy to rapid analysis of coals［J］．Spectroscopy Letters:An International Journal for Rapid Communication，2002，35(3):369 －376．

［6］M T Bona，J M Andres．ASTM clustering for improving coal analysis by near－ infrared spectroscopy［J］．Talanta，2006，70:711 －719．

［7］M T Bona，J M Andres．Coal analysis by diffuse reflectance near－infrared spectroscopy:Hierarchical cluster and linear discriminant analysis［J］．Talanta，2007，72(4):1423－1431．

［8］M T Bona，J M Andres．Application of chemometric tools for coal classification and multivariate calibration by transmission and drift mid － infrared spectroscopy［J］．Analytica Chimica Acta，2008，624:68 －78．

［9］Cheng Hongfei，Zhang Jinshan，et al．A spectroscopic comparison of selected Chinese kaolinite，coal bearing kaolinite and halloysite:a mid－infrared and near－infrared study［J］．Spectrochimica Acta Part A:Molecular and Biomolecular Spectroscopy，2010，77(4):856 －861．

［10］ S Wold．PLS-regression:a basic tool of chemometrics［J］．Chemometrics and Intelligent Laboratory Systems，2001，58:109－130．

［11］ Shao Xueguang，Wang Fang，Cheng Da，et al．A method for near-infrared spectral calibration of complex plant wavelet transform and elimination of uninformative variables［J］．Analytical and Bioanalytical Chemistry，2004，378(5):1382－1387．

［12］ Inge S．Helland，Tormod Nas，Tomas Isaksson．Related versions of the multiplicative scatter correction method for preprocessing spectroscopic data［J］．Chemometrics and Intelligent Laboratory Systems，1995，25，233 －241．