改進(jìn)和聲搜索算法的近紅外光譜特征變量選擇

張 磊，丁香乾，宮會麗，吳麗君，白曉莉，羅 林

1.中國海洋大學(xué)信息科學(xué)與工程學(xué)院，山東 青島 266100 2.云南中煙工業(yè)有限責(zé)任公司技術(shù)中心，云南 昆明 650024

引 言

近紅外光譜分析作為光譜測量和化學(xué)計(jì)量學(xué)結(jié)合的分析技術(shù)，是近年來分析化學(xué)領(lǐng)域發(fā)展最迅速的高新分析技術(shù)，以其簡便、快速、低成本、綠色環(huán)保、涵蓋信息量大和多組分同時(shí)測定等優(yōu)點(diǎn)，廣泛應(yīng)用于煙草、石油、紡織、食品等領(lǐng)域的定量檢測和定性分析，是現(xiàn)階段開展廣域范圍、大規(guī)模樣品檢測最佳的技術(shù)手段[1]。然而，近紅外光譜由于變量維度高、包含眾多與測定量無關(guān)的冗余信息和高頻噪聲等，對其直接建模不但會增加模型的復(fù)雜度，同時(shí)也會影響模型的預(yù)測性能和泛化能力[2]。因此，如何從高維的光譜變量中篩選與預(yù)測指標(biāo)密切相關(guān)的特征變量，從而構(gòu)建準(zhǔn)確、穩(wěn)定、高效的預(yù)測模型就顯得尤為重要。

目前，許多學(xué)者參與到近紅外光譜的特征變量篩選研究當(dāng)中，取得了一定成果。李倩倩等[3]應(yīng)用無信息變量消除法(uninformative variables elimination，UVE)剔除了光譜中不含有效信息的光譜點(diǎn)，然后構(gòu)建總氮和總糖的PLS定量模型。但由于過度依賴光譜和指標(biāo)關(guān)聯(lián)關(guān)系，并且需要單個(gè)變量逐一提取，容易忽略光譜變量全局貢獻(xiàn)度。徐寶鼎等[4]應(yīng)用粒子群算法(particle swarm optimization，PSO)進(jìn)行近紅外光譜特征選擇，采用特征分層的方式劃分初始粒子群。但優(yōu)化過程中只針對固定的初始粒子群內(nèi)部進(jìn)行優(yōu)選，忽略了不同粒子群光譜之間的比較。

鑒于此，本文提出了基于改進(jìn)和聲搜索(harmony search,HS)算法的光譜特征變量選擇方法。首先通過偏最小二乘(partial least square,PLS)載荷系數(shù)計(jì)算各光譜點(diǎn)的特征貢獻(xiàn)度，作為和聲搜索算法改進(jìn)的擾動(dòng)權(quán)重。然后應(yīng)用和聲搜索算法進(jìn)行最優(yōu)光譜變量選擇，為了避免陷入局部最優(yōu)或者收斂速度過慢，引入變量特征貢獻(xiàn)度(變量擾動(dòng)權(quán)重)對初始和聲記憶庫生成方法進(jìn)行優(yōu)化，并對和聲搜索算法的參數(shù)進(jìn)行動(dòng)態(tài)調(diào)整，提高了算法的自適應(yīng)能力，為近紅外光譜的特征變量選擇提供了更好的解決方案，有效降低了建模復(fù)雜度并提高了模型性能。

1 理論與算法

1.1 光譜變量特征貢獻(xiàn)度

為了計(jì)算光譜變量的特征貢獻(xiàn)度[5]，首先應(yīng)用PLS方法計(jì)算各光譜點(diǎn)的載荷(X-Loading)系數(shù)pi。通過定義光譜矩陣提取因子ti和預(yù)測指標(biāo)矩陣提取因子ui之間最大的協(xié)方差，求解如式(1)的最優(yōu)化問題

max{Cov(t1，u1)}=max〈E0w1，F(xiàn)0c1〉

(1)

利用拉格朗日乘數(shù)法求出w1和c1滿足式(2)

(2)

(3)

最后應(yīng)用式(4)計(jì)算每個(gè)光譜變量的貢獻(xiàn)度：其中，n為光譜維數(shù)。

(4)

1.2 和聲搜索算法的改進(jìn)和光譜特征變量選擇過程

和聲搜索算法是一種相對較為新穎的啟發(fā)式優(yōu)化算法，在2001年由Geem等提出，主要用于解決組合優(yōu)化和特征變量選擇問題，在許多研究領(lǐng)域展現(xiàn)了良好的性能。與演奏家們的目標(biāo)是合奏出最優(yōu)美的和聲一樣，HS的標(biāo)準(zhǔn)就是優(yōu)化問題的目標(biāo)函數(shù)[6]。

和聲搜索算法的改進(jìn)和特征變量選擇步驟如下。

1.2.1 初始化和聲記憶庫

和聲記憶庫(harmony memory，HM)中解向量的每一個(gè)分量都是用該決策變量的下界和上界之間均勻分布的隨機(jī)數(shù)(范圍[Lxi，Uxi]，1≤i≤N)來初始化[7]。第j個(gè)解向量的第i個(gè)分量通過式(5)給出

(5)

其中，j=1，2，3，…，HMS，rand[0,1]是0-1之間均勻分布的隨機(jī)數(shù)。

第j個(gè)解向量的目標(biāo)函數(shù)值用f(xj)來表示，置為解向量的最后一個(gè)分量。對應(yīng)HM的矩陣結(jié)構(gòu)由式(6)生成

HM(j，1：N)=xj

HM(j，N+1)=f(xj)

(6)

由于上述和聲記憶庫的初始解向量完全是隨機(jī)生成的，其在整個(gè)解空間Xi中的分布可能存在極端情況，經(jīng)過不斷迭代很容易陷入局部最優(yōu)的狀態(tài)，在一定程度上影響算法的搜索性能[8]。為了增強(qiáng)全局搜索能力，在各分量生成時(shí)引入該分量的特征貢獻(xiàn)度gi作為激勵(lì)因子，采用隨機(jī)遍歷和激勵(lì)因子共同作用的方式生成初始解向量，如式(7)所示

(7)

本工作中，輸入變量個(gè)數(shù)為1 555個(gè)光譜點(diǎn)，各光譜點(diǎn)依次橫向在和聲向量中排列。如果該光譜點(diǎn)被選定，則其值設(shè)為1，否則設(shè)為0。每一行的最后一位表示該解向量目標(biāo)函數(shù)。HM矩陣結(jié)構(gòu)如圖1所示。

圖1 和聲記憶庫(HM)矩陣結(jié)構(gòu)Fig.1 Harmony memory matrix structure

1.2.2 從和聲記憶庫生成新的光譜特征選擇向量

(1)在[0-1]之間生成一個(gè)隨機(jī)數(shù)r1，與和聲庫取值概率(harmony memory considering rate,HMCR)進(jìn)行比較，若r1

(8)

(2)若這個(gè)和聲變量是從HM中得到的，則對這個(gè)和聲變量進(jìn)行微調(diào)。在[0-1]之間生成一個(gè)隨機(jī)數(shù)r2，與音調(diào)微調(diào)概率(pitch adjusting rate,PAR)進(jìn)行比較，若r2

(9)

其中，r3是[0-1]之間的隨機(jī)數(shù)，若r3>0.5，取“+”號；若r3≤0.5，取“-”號。

(10)

傳統(tǒng)HS算法中，BW在整個(gè)迭代過程中是不變的，在某種程度上會影響搜索過程的遍歷性和種群的多樣性，容易陷入局部最優(yōu)[11]。因此，本文以光譜變量特征貢獻(xiàn)度gi作為懲罰項(xiàng)，引入α作為BW和gi的平衡因子，并根據(jù)式(11)使BW隨迭代次數(shù)而動(dòng)態(tài)調(diào)整，從而適應(yīng)光譜選擇問題特定階段的搜索。

BW=(1-α)BW-α(BW-gi)

(11)

其中，t為當(dāng)前迭代次數(shù)，T為最大迭代數(shù)。產(chǎn)生新和聲向量的過程如圖2所示。

圖2 新和聲向量生成過程 P(E1),P(E2),P(E3):通過多概率擾動(dòng)生成新和聲；和聲庫取值概率(E1)，音調(diào)微調(diào)概率(E2)，S型函數(shù)概率(E3)Fig.2 New harmony vectors generation process

P(E1),P(E2),P(E3):Probabilities to improvise a new harmony vector by memory considerations HMCR (E1),pitch adjustment rate PAR (E2),and activation function Sigmoid (E3)

1.2.3 應(yīng)用光譜特征選擇向量構(gòu)建PLS預(yù)測模型

光譜特征向量中各分量表示光譜點(diǎn)選擇狀態(tài)，非零列表示該光譜點(diǎn)被選中。最后一個(gè)分量應(yīng)用各訓(xùn)練子集均方根誤差(root mean square error,RMSE)的均值[12]作為光譜特征向量的目標(biāo)函數(shù)，如式(12)所示

(12)

其中，n為樣本個(gè)數(shù)，ei為第i個(gè)樣本的預(yù)測誤差，等于第i個(gè)樣本的預(yù)測值pi和真實(shí)值ai之間的差值。

針對各訓(xùn)練子集數(shù)據(jù)，基于HM中每個(gè)光譜特征向量對應(yīng)選擇的特征光譜點(diǎn)，分別構(gòu)建PLS預(yù)測模型，最終計(jì)算各訓(xùn)練子集模型RMSE的均值作為目標(biāo)函數(shù)值。如果該目標(biāo)函數(shù)值比HM中最差的光譜特征向量好，即f(Xnew)

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

2.1 樣本數(shù)據(jù)

為驗(yàn)證算法的有效性，共選取800個(gè)具有代表性的煙葉樣品，隨機(jī)選取其中的600個(gè)樣品作為訓(xùn)練集，剩余的200個(gè)作為測試集。進(jìn)一步，又將訓(xùn)練集隨機(jī)劃分成3個(gè)數(shù)量相等的訓(xùn)練子集，用于增強(qiáng)變量選擇算法的泛化能力。樣本數(shù)據(jù)分布情況如表1所示。

表1 樣本數(shù)據(jù)分布情況Table 1 Sample data distribution

2.2 儀器設(shè)備與采集方法

選用布魯克公司的MATRIX-I型傅里葉變換近紅外光譜儀采集近紅外光譜，掃描范圍為4 000～10 000 cm-1，分辨率為8 cm-1，掃描次數(shù)為64。將煙葉樣品放置在60 ℃烘箱中烘干2 h，用旋風(fēng)磨磨成粉末，過40目篩，放入干燥皿中。每個(gè)樣品稱重15 g，放置于干凈的5 cm樣品杯中，用壓樣器自然壓實(shí)后進(jìn)行近紅外光譜掃描。為保證光譜一致性，每個(gè)樣品均重復(fù)裝樣掃描三次，采用三次掃描的平均光譜作為該樣品的最終光譜，樣品三次掃描的平均誤差為：1.5×10-3A，標(biāo)準(zhǔn)差為：0.4×10-3。煙堿、總糖、總氮三個(gè)化學(xué)指標(biāo)含量均按照煙草行業(yè)規(guī)定的標(biāo)準(zhǔn)方法測定。

2.3 光譜預(yù)處理

采用一階導(dǎo)數(shù)加Savitzky-Golay平滑的方法進(jìn)行光譜預(yù)處理[14]，移動(dòng)窗口寬度為9，多項(xiàng)式數(shù)為3。圖3為采用2.2中方法采集的原始光譜圖，圖4為預(yù)處理后光譜圖。

圖3 原始光譜圖Fig.3 Original spectra

3 結(jié)果與討論

為充分驗(yàn)證算法的有效性，選取煙草化學(xué)分析中比較重要且常用的煙堿、總糖、總氮三個(gè)指標(biāo)，分別開展對應(yīng)的光譜特征變量選擇和預(yù)測建模實(shí)驗(yàn)。

圖4 預(yù)處理后光譜圖Fig.4 Pre-processed spectra

3.1 光譜變量特征貢獻(xiàn)度計(jì)算

在600個(gè)樣品的訓(xùn)練集上，采用近紅外光譜的全譜段和各指標(biāo)檢測結(jié)果構(gòu)建PLS模型，得到第一主成分下各光譜點(diǎn)的載荷矢量，然后求得各光譜變量的特征貢獻(xiàn)度。不同預(yù)測指標(biāo)對應(yīng)光譜變量的特征貢獻(xiàn)度均不相同，各指標(biāo)光譜變量貢獻(xiàn)度如圖5所示。

3.2 光譜特征變量選擇結(jié)果

改進(jìn)和聲搜索算法根據(jù)各訓(xùn)練子集RMSE均值最小化作為約束標(biāo)準(zhǔn)進(jìn)行光譜特征變量選擇，以向量中取值為1的分量作為篩選出的變量。算法在3個(gè)訓(xùn)練子集上各運(yùn)行100次。以煙堿指標(biāo)為例，圖6表示1 555個(gè)變量在總計(jì)300次訓(xùn)練中被選擇的累積頻次，頻次越高說明該光譜變量對預(yù)測指標(biāo)越重要，可以選做特征變量。將各變量被選擇的累積頻次從高到低排序，以25為下降梯度不斷增加變量數(shù)，分別計(jì)算訓(xùn)練集上PLS模型預(yù)測結(jié)果的RMSEC。圖7為隨變量數(shù)的增加RMSEC變化趨勢圖。由圖可知，開始時(shí)，隨著變量數(shù)增加，RMSEC逐漸降低，當(dāng)選擇變量頻次為200時(shí)達(dá)到最小值。如果變量繼續(xù)增加，會引入噪聲和冗余信息，RMSEC隨之變大，對模型效果產(chǎn)生不利影響。因此以累積頻次200作為臨界點(diǎn)，最終得到煙堿指標(biāo)光譜的特征變量個(gè)數(shù)為199。對于總糖和總氮指標(biāo)，臨界點(diǎn)分別200和150，光譜特征變量個(gè)數(shù)分別為213和276。

圖5 不同指標(biāo)對應(yīng)的光譜變量貢獻(xiàn)度(a):煙堿;(b):總糖;(c):總氮Fig.5 Contributions of each spectral variable(a):Nicotine;(b):Total sugar;(c):Total nitrogen

圖6 光譜變量被選擇的累積頻次Fig.6 Cumulative frequency of spectral variables selected

3.3 預(yù)測結(jié)果對比

為了驗(yàn)證算法的有效性，將改進(jìn)和聲搜索算法與全光譜以及UVE和PSO等光譜特征變量選擇方法進(jìn)行對比，分別構(gòu)建煙堿、總糖、總氮三個(gè)指標(biāo)的PLS預(yù)測模型，以選擇的光譜特征變量數(shù)、訓(xùn)練集和測試集的決定系數(shù)及均方根誤差作為評價(jià)標(biāo)準(zhǔn)，結(jié)果如表2所示。選擇多次實(shí)驗(yàn)中最小的PRESS值對應(yīng)的主因子數(shù)作為最終PLS模型的主因子數(shù)。

圖7 RMSEC隨變量增加變化趨勢Fig.7 RMSEC trends with variables increasing

表2 各類型特征光譜在三個(gè)指標(biāo)訓(xùn)練集和預(yù)測集上的性能Table 2 Performance of training and prediction sets using various types of characteristic spectra of three indicators

決定系數(shù)R2越大，預(yù)測均方根誤差RMSEP越小，代表模型的預(yù)測性能越好。由表2可見，各種特征變量選擇方法相對全光譜來說，均減少了波長點(diǎn)數(shù)，同時(shí)也提高了預(yù)測精度。相較而言，本文提出的改進(jìn)和聲搜索算法對三個(gè)預(yù)測指標(biāo)都提取了最少的光譜特征變量(提取的光譜特征變量只占全光譜的13%～18%)，具有最高的R2和最小的RMSEP值。同時(shí)，經(jīng)過專家對篩選譜段對應(yīng)的官能團(tuán)分析，所選譜段均能夠較好的反應(yīng)其指標(biāo)項(xiàng)特征信息，充分表明了本算法提取的光譜特征變量能有效降低冗余信息和噪聲、消減變量間的多重共線性，使得模型更加穩(wěn)健、泛化能力更強(qiáng)。

4 結(jié) 論

基于改進(jìn)的和聲搜索算法提出了一種近紅外光譜特征變量選擇方法。首先利用PLS載荷系數(shù)計(jì)算光譜變量對預(yù)測指標(biāo)的貢獻(xiàn)度，作為變量權(quán)重。然后，利用和聲搜索算法進(jìn)行特征變量篩選，過程中引入變量權(quán)重對和聲搜索算法的初始化和參數(shù)動(dòng)態(tài)調(diào)整進(jìn)行改進(jìn)。最后針對篩選后的變量通過PLS建模在煙堿、總糖、總氮三個(gè)指標(biāo)的訓(xùn)練集和測試集進(jìn)行驗(yàn)證，并與全光譜和幾種常用的特征變量選擇方法進(jìn)行對比。實(shí)驗(yàn)結(jié)果表明，采用本算法進(jìn)行光譜特征變量選擇對模型性能優(yōu)于全光譜和其他光譜變量選擇方法，由此說明對近紅外光譜進(jìn)行特征變量選擇的必要性以及本算法的有效性，既保證了模型預(yù)測性能又降低了建模的復(fù)雜度，為近紅外定量預(yù)測模型的構(gòu)建和優(yōu)化提供了參考。