變量影響值與集群分析相結合的中紅外波長選擇方法

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

1. 西安交通大學電力設備電氣絕緣國家重點實驗室， 陜西 西安 710049 2. 中煤科工集團重慶設計研究院(集團)有限公司， 重慶 400042

引 言

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

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

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

1 IVPA算法基本原理

1.1 樣本空間中變量影響值

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

y=Xβ+e

(1)

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

IVi=UEi-DEi

(2)

1.2 IVPA變量選擇方法

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

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

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

(3)

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

Step3： 利用指數(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為當前的迭代次數(shù)，N為設置的總迭代次數(shù)，p為上次迭代保留的變量數(shù)量。 通過剩余變量數(shù)可以計算出每次迭代中需要剔除的變量。 當普通變量的數(shù)目大于剔除的變量數(shù)時， 剔除普通變量中頻率低的變量； 當普通變量數(shù)目小于剔除的變量時， 剔除全部的普通變量， 并且從精英變量中剔除平均影響值低的變量；

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

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

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

2 實驗部分

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

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

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

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

2.2 IVPA參數(shù)確定

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

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

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

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

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

(3)變量空間中子模型個數(shù)W與比例系數(shù)β的確定。 通常W的值設置較大， 保證了每個變量都會分配到對應的子模型中， 體現(xiàn)了集群分析中群體的理念。β通常設置為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次計算后， 以RMSE最小值所對應的W與β作為最終選擇的參數(shù)。 圖4為RMSE隨W與β的變化趨勢圖。 由圖可知， 當W為3 500，β等于0.05時， RMSE最小。 因此， 最終獲得的參數(shù)為： 算法迭代次數(shù)N為80、 系數(shù)α1為0.7， 系數(shù)α2為1.9， 子模型個數(shù)W為3 500， 預測效果較好的模型占全部模型比例系數(shù)為0.05。

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

3 結果與討論

為了評價算法的性能， 以光譜交疊最嚴重的異丁烷光譜為例， SCARS， IVSO及IVPA 3種變量選擇方法進行變量選擇， 選擇的變量如圖5所示。

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

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

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

表1 三種方法對500 ppm烷烴氣體的預測結果

4 結 論

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