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

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

引 言

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

目前，許多學(xué)者參與到近紅外光譜的特征變量篩選研究當(dāng)中，取得了一定成果。李倩倩等[3]應(yīng)用無信息變量消除法(uninformative variables elimination，UVE)剔除了光譜中不含有效信息的光譜點，然后構(gòu)建總氮和總糖的PLS定量模型。但由于過度依賴光譜和指標(biāo)關(guān)聯(lián)關(guān)系，并且需要單個變量逐一提取，容易忽略光譜變量全局貢獻(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ù)計算各光譜點的特征貢獻(xiàn)度，作為和聲搜索算法改進(jìn)的擾動權(quán)重。然后應(yīng)用和聲搜索算法進(jìn)行最優(yōu)光譜變量選擇，為了避免陷入局部最優(yōu)或者收斂速度過慢，引入變量特征貢獻(xiàn)度(變量擾動權(quán)重)對初始和聲記憶庫生成方法進(jìn)行優(yōu)化，并對和聲搜索算法的參數(shù)進(jìn)行動態(tài)調(diào)整，提高了算法的自適應(yīng)能力，為近紅外光譜的特征變量選擇提供了更好的解決方案，有效降低了建模復(fù)雜度并提高了模型性能。

1 理論與算法

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

為了計算光譜變量的特征貢獻(xiàn)度[5]，首先應(yīng)用PLS方法計算各光譜點的載荷(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)計算每個光譜變量的貢獻(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)中解向量的每一個分量都是用該決策變量的下界和上界之間均勻分布的隨機(jī)數(shù)(范圍[Lxi，Uxi]，1≤i≤N)來初始化[7]。第j個解向量的第i個分量通過式(5)給出

(5)

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

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

HM(j，1：N)=xj

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

(6)

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

(7)

本工作中，輸入變量個數(shù)為1 555個光譜點，各光譜點依次橫向在和聲向量中排列。如果該光譜點被選定，則其值設(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]之間生成一個隨機(jī)數(shù)r1，與和聲庫取值概率(harmony memory considering rate,HMCR)進(jìn)行比較，若r1

(8)

(2)若這個和聲變量是從HM中得到的，則對這個和聲變量進(jìn)行微調(diào)。在[0-1]之間生成一個隨機(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在整個迭代過程中是不變的，在某種程度上會影響搜索過程的遍歷性和種群的多樣性，容易陷入局部最優(yōu)[11]。因此，本文以光譜變量特征貢獻(xiàn)度gi作為懲罰項，引入α作為BW和gi的平衡因子，并根據(jù)式(11)使BW隨迭代次數(shù)而動態(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):通過多概率擾動生成新和聲；和聲庫取值概率(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ù)測模型

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

(12)

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

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

2 實驗部分

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

為驗證算法的有效性，共選取800個具有代表性的煙葉樣品，隨機(jī)選取其中的600個樣品作為訓(xùn)練集，剩余的200個作為測試集。進(jìn)一步，又將訓(xùn)練集隨機(jī)劃分成3個數(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目篩，放入干燥皿中。每個樣品稱重15 g，放置于干凈的5 cm樣品杯中，用壓樣器自然壓實后進(jìn)行近紅外光譜掃描。為保證光譜一致性，每個樣品均重復(fù)裝樣掃描三次，采用三次掃描的平均光譜作為該樣品的最終光譜，樣品三次掃描的平均誤差為：1.5×10-3A，標(biāo)準(zhǔn)差為：0.4×10-3。煙堿、總糖、總氮三個化學(xué)指標(biāo)含量均按照煙草行業(yè)規(guī)定的標(biāo)準(zhǔn)方法測定。

2.3 光譜預(yù)處理

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

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

3 結(jié)果與討論

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

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

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

在600個樣品的訓(xùn)練集上，采用近紅外光譜的全譜段和各指標(biāo)檢測結(jié)果構(gòu)建PLS模型，得到第一主成分下各光譜點的載荷矢量，然后求得各光譜變量的特征貢獻(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個訓(xùn)練子集上各運行100次。以煙堿指標(biāo)為例，圖6表示1 555個變量在總計300次訓(xùn)練中被選擇的累積頻次，頻次越高說明該光譜變量對預(yù)測指標(biāo)越重要，可以選做特征變量。將各變量被選擇的累積頻次從高到低排序，以25為下降梯度不斷增加變量數(shù)，分別計算訓(xùn)練集上PLS模型預(yù)測結(jié)果的RMSEC。圖7為隨變量數(shù)的增加RMSEC變化趨勢圖。由圖可知，開始時，隨著變量數(shù)增加，RMSEC逐漸降低，當(dāng)選擇變量頻次為200時達(dá)到最小值。如果變量繼續(xù)增加，會引入噪聲和冗余信息，RMSEC隨之變大，對模型效果產(chǎn)生不利影響。因此以累積頻次200作為臨界點，最終得到煙堿指標(biāo)光譜的特征變量個數(shù)為199。對于總糖和總氮指標(biāo)，臨界點分別200和150，光譜特征變量個數(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é)果對比

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

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

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

4 結(jié) 論

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