張 峰， 湯曉君*， 仝昂鑫， 王 斌， 湯春瑞， 王 杰

1. 西安交通大學(xué)電力設(shè)備電氣絕緣國(guó)家重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710049 2. 中煤科工集團(tuán)重慶設(shè)計(jì)研究院(集團(tuán))有限公司， 重慶 400042

引 言

紅外光譜分析技術(shù)由于其響應(yīng)速度快、 分析成分多、 靈敏度高、 無損檢測(cè)等優(yōu)點(diǎn)， 廣泛應(yīng)用于食品質(zhì)量檢測(cè)、 電力設(shè)備故障診斷、 石油化工、 煤礦災(zāi)害預(yù)警等領(lǐng)域[1-4]。 隨著光譜儀分辨率的逐步提高， 每個(gè)樣本獲得的譜線多達(dá)幾千條。 這些譜線無可避免的會(huì)包括無用變量甚至是干擾變量。 通常， 當(dāng)分析物為單組分物質(zhì)時(shí)， 選擇主吸收峰位置處的單個(gè)譜線利用線性或者非線性擬合方法就能得到滿意的結(jié)果， 但是當(dāng)分析物為混合物時(shí)， 尤其是混合物之間的吸收光譜重疊嚴(yán)重時(shí)， 隨機(jī)選擇幾個(gè)變量不僅會(huì)增加對(duì)其他組分的交叉靈敏度， 甚至?xí)档捅旧斫M分的預(yù)測(cè)精度。 因此， 光譜變量特征選擇對(duì)于降低模型運(yùn)行時(shí)間， 提高模型的預(yù)測(cè)精度及降低對(duì)其他組分的交叉靈敏度具有重要意義[5-8]。

偏最小二乘法(partial least squares, PLS)具有運(yùn)行時(shí)間短與優(yōu)化的參數(shù)少等特點(diǎn)， 廣泛應(yīng)用于線性分析模型中[10]。 在PLS模型中， 回歸系數(shù)β是一個(gè)重要參數(shù)， 基于該參數(shù)， 演變了多種變量選擇方法， 如蒙特卡洛無信息變量消除法(Monte Carlo non-information variable elimination， MCUVE)[10]、 穩(wěn)定性競(jìng)爭(zhēng)自適應(yīng)重加權(quán)采樣法(stability competitive adaptive reweighted sampling, SCARS)[11-12]、 變量子集迭代優(yōu)化(iteratively variable subset optimization, IVSO)[13]。 其中MCUVE與SCARS方法類似， 均是通過蒙特卡洛方法產(chǎn)生大量的樣本空間， 獲取變量的穩(wěn)定性， 根據(jù)穩(wěn)定性的值來進(jìn)行變量篩選， 不同的是， SCARS利用了指數(shù)衰減函數(shù)來逐步對(duì)變量進(jìn)行剔除， 而MCUVE則是通過閾值來控制變量的選擇個(gè)數(shù)， 閾值的設(shè)定采用遍歷方法來進(jìn)行優(yōu)化選擇。 然而， 這兩種方法均未考慮到變量之間的相互影響， 當(dāng)變量之間的共線性較高時(shí)， 兩種方法選擇的變量效果并不好。 IVSO算法利用蒙特卡洛方法在變量空間中生成大量的子集， 對(duì)變量空間分別建立PLS模型， 對(duì)每個(gè)模型中的回歸系數(shù)進(jìn)行歸一化處理， 隨后將這些模型的回歸系數(shù)進(jìn)行求和并歸一化來得到新的權(quán)重值。 利用加權(quán)二進(jìn)制采樣技術(shù)根據(jù)權(quán)重值的大小來分配在下次迭代中變量被選擇的概率， 這樣回歸系數(shù)值較大的變量更有機(jī)會(huì)在下次迭代中被選擇。 然而回歸系數(shù)會(huì)隨著樣本的不同而變化， 這樣導(dǎo)致了IVSO算法選擇的變量隨機(jī)性大。

針對(duì)上述變量選擇方法存在的問題， 提出了一種基于變量影響值與集群分析相結(jié)合(impact of variable and population analysis, IVPA)的變量選擇方法。 為了評(píng)價(jià)提出方法的性能， 將該方法應(yīng)用于烷烴類氣體的中紅外數(shù)據(jù)集中， 并與SCARS和IVSO方法的預(yù)測(cè)結(jié)果進(jìn)行了對(duì)比。 結(jié)果表明， IVPA方法對(duì)異丁烷的預(yù)測(cè)最準(zhǔn)確， 對(duì)其他組分的交叉靈敏度最低， 能夠?qū)庾V交疊嚴(yán)重的組分進(jìn)行變量選擇。

1 IVPA算法基本原理

1.1 樣本空間中變量影響值

將掃描獲得的光譜數(shù)據(jù)用矩陣Xn×p表示，n為樣本的個(gè)數(shù)，p為譜線數(shù)量。yn×m為n個(gè)樣本所對(duì)應(yīng)的分析物濃度信息。 當(dāng)分析物為單組分時(shí)，m取值為1。 當(dāng)分析物為多組分時(shí)，m大于1， 此時(shí)可以用PLS2模型或者利用PLS1模型經(jīng)過m次循環(huán)來建立模型。 采用PLS1建立紅外光譜定量分析模型可以表示如式(1)

y=Xβ+e

(1)

式(1)中，β為回歸系數(shù)向量， 表示為β=[β1,β2,…,βp]T，y為分析第i種組分時(shí)的濃度向量， 其中i≤m，e為隨機(jī)誤差向量。 將其中一個(gè)變量在其值的基礎(chǔ)上乘以小于1的系數(shù)α1， 得到新的變量V1， 然后乘以大于1的系數(shù)α2， 得到變量V2， 對(duì)獲取的新的變量分別建立PLS模型， 得到對(duì)應(yīng)的交互驗(yàn)證均方根誤差UEi與DEi。 該過程循環(huán)p次， 這樣每個(gè)樣變量的影響值可以由式(2)計(jì)算

IVi=UEi-DEi

(2)

1.2 IVPA變量選擇方法

IVPA算法融合了變量的影響值與變量空間中最優(yōu)模型出現(xiàn)的頻率兩個(gè)重要指標(biāo)。 在每次迭代過程中， 同時(shí)計(jì)算變量的影響值與頻率， 根據(jù)變量影響值采用加權(quán)自舉采樣技術(shù)將變量劃分為精英變量與普通變量， 在普通變量中剔除頻率較低的變量， 這樣保證了每次迭代中剔除了變量影響值較低與最優(yōu)模型中出現(xiàn)的頻率較低的變量。 當(dāng)?shù)_(dá)到設(shè)定的次數(shù)時(shí)， 算法結(jié)束。 IVPA算法詳細(xì)步驟如下：

Step1： 迭代開始， 初始變量個(gè)數(shù)為p。 依據(jù)式(2)， 計(jì)算每個(gè)變量的影響值IVi， 根據(jù)IVi的值采用加權(quán)自舉采樣技術(shù)將變量劃分為精英變量與普通變量。 此時(shí)， 精英變量的個(gè)數(shù)約為p的0.632倍[14]；

Step2： 根據(jù)集群分析理念， 采用蒙特卡洛算法從p個(gè)變量中隨機(jī)選擇p1個(gè)變量， 生成W個(gè)變量空間，p1值的可以選擇為變量個(gè)數(shù)的0.2倍， 降低變量之間的共線性帶來的影響。 對(duì)W個(gè)變量空間建立PLS模型， 獲取每個(gè)模型的交叉均方根誤差(root mean squared error of cross validation, RMSECV)值， 將W個(gè)RMSECV值從小到大進(jìn)行排序， 選擇預(yù)測(cè)結(jié)果較好的βW個(gè)模型， 根據(jù)式(3)計(jì)算每個(gè)變量出現(xiàn)的頻率fi

(3)

式(3)中，F(xiàn)ij為變量i是否出現(xiàn)在預(yù)測(cè)結(jié)果較好的第j個(gè)模型中， 若出現(xiàn)， 則為1， 否則為0，β為比例系數(shù)， 取值小于0.5。

Step3： 利用指數(shù)衰減函數(shù)確定每次迭代時(shí)保留的變量數(shù)目， 指數(shù)衰減函數(shù)可以表示如式(4)

ri=ae-ki

(4)

式(4)中，ri為第i迭代后剩余的變量數(shù)， 其中a=(p/2)(1/(N-1)),k=ln(p/2)/(N-1)，i為當(dāng)前的迭代次數(shù)，N為設(shè)置的總迭代次數(shù)，p為上次迭代保留的變量數(shù)量。 通過剩余變量數(shù)可以計(jì)算出每次迭代中需要剔除的變量。 當(dāng)普通變量的數(shù)目大于剔除的變量數(shù)時(shí)， 剔除普通變量中頻率低的變量； 當(dāng)普通變量數(shù)目小于剔除的變量時(shí)， 剔除全部的普通變量， 并且從精英變量中剔除平均影響值低的變量；

Step4： 記錄每次迭代過程中變量空間中W個(gè)模型的最小均方根誤差(root mean squared error, RMSE)， 同時(shí)更新變量數(shù)p的值，p=ri；

Step5： 若迭代次數(shù)小于設(shè)定值時(shí)， 返回step1， 否則， 執(zhí)行Step6；

Step6： 選擇最小的RMSE所對(duì)應(yīng)的變量子集作為最優(yōu)變量組合。

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

2.1 光譜數(shù)據(jù)的獲取

實(shí)驗(yàn)所采用的傅里葉變換紅外光譜儀型號(hào)為Spectrum Two， 其氣體池長(zhǎng)度為10 cm， 掃描范圍400～4 000 cm-1, 波數(shù)分辨率為1 cm-1， 掃描類型為吸光度。 為了降低隨機(jī)噪聲的干擾， 每個(gè)樣品掃描次數(shù)設(shè)置為8。 待分析物的目標(biāo)氣體為甲烷(CH4)、 乙烷(C2H6)、 丙烷(C3H8)、 異丁烷(i-C4H10)、 正丁烷(n-C4H10)共5種烷烴類氣體。 實(shí)驗(yàn)時(shí)， 采用標(biāo)氣進(jìn)行光譜掃描， 五種氣體的濃度范圍分別為0.01%， 0.02%， 0.05%， 0.1%， 0.2%和0.5%。 五種組分氣體在0.05%下的吸光度光譜如圖1所示。

圖1 五種氣體在濃度為0.05%時(shí)的吸光度光譜

從圖1中可以看出， 5種氣體在3 000 cm-1附近吸收峰最高， 為烷烴類氣體的主吸收峰， 在1 200～1 700 cm-1范圍內(nèi)， 為烷烴類氣體的次吸收峰。 可以看出， 無論是在主吸收峰還是在次吸收峰， 5種氣體之間的光譜交疊嚴(yán)重。 從圖中還可以看出， 5種氣體的光譜發(fā)生不同程度的基線漂移， 根據(jù)前期研究成果[15]， 采用一種懲罰最小二乘方法對(duì)光譜數(shù)據(jù)進(jìn)行了基線校正。

2.2 IVPA參數(shù)確定

IVPA算法中需要優(yōu)化的參數(shù)有： (1)算法的迭代次數(shù)N； (2)樣本空間中變量乘以小于1的系數(shù)α1； (3)樣本空間中變量乘以大于1的系數(shù)α2； (4)變量空間中子模型的個(gè)數(shù)W； (5)變量空間中預(yù)測(cè)效果較好的模型占全部模型W的比例系數(shù)β。 可以按照以下步驟來進(jìn)行最優(yōu)參數(shù)選取。

(1)系數(shù)α1與α2的優(yōu)化。 固定算法中其他參數(shù)的值，N設(shè)置為50，W設(shè)置為1 000，β設(shè)置為0.1。 同時(shí)α1在0.1到1之間以間隔為0.1變化，α2在1.1到3.1之間以間隔為0.2變化， 共經(jīng)歷99次IVPA運(yùn)算后， 獲得99個(gè)RMSE值， 選擇RMSE最小的值對(duì)應(yīng)的α1與α1作為最終優(yōu)化的值。 RMSE值隨α1與α2變化趨勢(shì)如圖2所示。 從圖2中可以看出， 當(dāng)α1取值為0.7，α2為1.9時(shí)， 獲得的RMSE值最小。

圖2 IVPA算法中參數(shù)α1與α2的優(yōu)化選擇

(2)迭代次數(shù)N的選取。 將第一步獲得的最優(yōu)α1與α2值作為IVPA作為常量， 同樣，W設(shè)置為1 000，β設(shè)置為0.1。N的取值為20到160， 間隔為20，N每次取值時(shí)， IVPA算法循環(huán)執(zhí)行30次。 經(jīng)歷240次IVPA循環(huán)后， 獲得8組N取不同值時(shí)的RMSE值， 每組RMSE個(gè)數(shù)為30個(gè)。 圖3為RMSE隨N的變化趨勢(shì)。 可知， 當(dāng)N為80時(shí)， RMSE獲得的均值最低， 并且上下波動(dòng)幅度不大。

圖3 IVPA算法中參數(shù)N的優(yōu)化選擇

(3)變量空間中子模型個(gè)數(shù)W與比例系數(shù)β的確定。 通常W的值設(shè)置較大， 保證了每個(gè)變量都會(huì)分配到對(duì)應(yīng)的子模型中， 體現(xiàn)了集群分析中群體的理念。β通常設(shè)置為0.05到0.3之間， 為集群分析中最優(yōu)選擇的概念。 將前2步選擇的α1，α2與N作為已知參數(shù)，β的取值范圍為0.05～0.3之間， 間隔為0.05，W以間隔500在1 000到5 000之間取值。 這樣， 經(jīng)過60次計(jì)算后， 以RMSE最小值所對(duì)應(yīng)的W與β作為最終選擇的參數(shù)。 圖4為RMSE隨W與β的變化趨勢(shì)圖。 由圖可知， 當(dāng)W為3 500，β等于0.05時(shí)， RMSE最小。 因此， 最終獲得的參數(shù)為： 算法迭代次數(shù)N為80、 系數(shù)α1為0.7， 系數(shù)α2為1.9， 子模型個(gè)數(shù)W為3 500， 預(yù)測(cè)效果較好的模型占全部模型比例系數(shù)為0.05。

圖4 IVPA算法中參數(shù)W與β的優(yōu)化選擇

3 結(jié)果與討論

為了評(píng)價(jià)算法的性能， 以光譜交疊最嚴(yán)重的異丁烷光譜為例， SCARS， IVSO及IVPA 3種變量選擇方法進(jìn)行變量選擇， 選擇的變量如圖5所示。

圖5 三種方法選擇的變量分布圖

從圖5中可以看出， SCARS選擇的變量最多， 并且分布比較散， 除此之外， 還選擇了無吸收峰范圍內(nèi)的變量， 這些變量被認(rèn)為是干擾變量與無用變量； IVSO選擇的變量次之， 集中在2 950 cm-1附近； IVPA選擇的變量最少， 選擇的變量覆蓋了異丁烷的吸收峰的主要區(qū)域。

將濃度為0.05%的五種烷烴氣體的光譜作為輸入變量， 根據(jù)以上述3種變量選擇方法建立的模型對(duì)五種氣體濃度進(jìn)行預(yù)測(cè)， 其預(yù)測(cè)結(jié)果如表1所示。 從表中可以看出3種方法對(duì)正丁烷的交叉靈敏度最高， 這是因?yàn)檎⊥榕c異丁烷的分子結(jié)構(gòu)近似， 其吸收光譜形狀接近； 對(duì)甲烷的交叉靈敏度較低， 這是由于甲烷在大于3 000 cm-1范圍的吸光度較低， 3種方法選擇的變量集中在大于3 000 cm-1范圍內(nèi)。 還可以看出IVPA對(duì)異丁烷的預(yù)測(cè)結(jié)果最準(zhǔn)確， 最大交叉靈敏度為1.02%， 最小交叉靈敏度為0.11%。 SCARS與IVSO對(duì)異丁烷預(yù)測(cè)的結(jié)果相當(dāng)， 除此之外， IVPA對(duì)其他4種氣體的交叉靈敏度最低， 表明所提出的方法能夠有效的對(duì)光譜交疊嚴(yán)重的變量進(jìn)行提取， 并能很好的對(duì)其進(jìn)行分析。

表1 三種方法對(duì)500 ppm烷烴氣體的預(yù)測(cè)結(jié)果

4 結(jié) 論

利用傅里葉變換紅外光譜儀對(duì)烷烴類氣體進(jìn)行定量分析時(shí)， 分析物的種類多且光譜交疊嚴(yán)重， 嚴(yán)重制約著定量分析結(jié)果的準(zhǔn)確性及運(yùn)算效率。 提出了一種基于變量影響值與集群分析的變量選擇方法， 將該方法與SCARS、 IVSO三種變量選擇方法來對(duì)烷烴類5種氣體進(jìn)行了變量提取， 以異丁烷為例， 對(duì)其進(jìn)行了預(yù)測(cè)， 并分析了對(duì)甲烷、 乙烷、 丙烷、 正丁烷的交叉靈敏度。 結(jié)果表明， IVPA方法對(duì)異丁烷的預(yù)測(cè)結(jié)果最好， 對(duì)其他氣體的交叉靈敏度最低， 能夠?qū)化B嚴(yán)重的光譜進(jìn)行變量提取， 具有實(shí)際的應(yīng)用價(jià)值。