陳 勇，吳彩娥，熊智新*

（南京林業(yè)大學(xué)輕工與食品學(xué)院，江蘇 南京 210037）

小麥粉中蛋白質(zhì)質(zhì)量分?jǐn)?shù)約為7%～15%，決定著小麥粉的營養(yǎng)品質(zhì)和加工品質(zhì)，所以小麥粉蛋白含量的檢測具有重要意義。小麥蛋白含量測定有化學(xué)法和物理法，前者雖然具有較高的準(zhǔn)確度，但需要對實驗材料進(jìn)行繁雜的化學(xué)處理，比較耗時費力，同時實驗使用的化學(xué)藥品會造成一定的環(huán)境污染。與傳統(tǒng)的化學(xué)方法相比，近紅外光譜分析是一種綠色、無損快速檢測技術(shù)，已經(jīng)廣泛應(yīng)用于農(nóng)業(yè)、食品、醫(yī)藥、化工等行業(yè)。目前，國內(nèi)外學(xué)者利用近紅外光譜法對小麥粉中蛋白質(zhì)、灰分和水分等含量進(jìn)行了檢測分析。然而，近紅外光譜有吸收信號弱、譜峰重疊嚴(yán)重以及易受外界環(huán)境干擾等缺點，并且一條光譜往往包含著數(shù)量眾多的波長點，這給建立高質(zhì)量的預(yù)測模型帶來了很大的挑戰(zhàn)。通過合適的波長選擇方法可以在眾多近紅外光譜波長中篩選出具有特征信息的波長，并對數(shù)據(jù)進(jìn)行有效地降維，降低建立預(yù)測模型的復(fù)雜程度，在一定程度上可提高預(yù)測模型的準(zhǔn)確性。

常見的波長選擇方法有無信息變量消除法、競爭自適應(yīng)加權(quán)采樣法、連續(xù)投影算法以及群體智能優(yōu)化算法。群體智能優(yōu)化算法具有強大的全局搜索能力，使其在特征變量篩選方面潛力巨大，其中遺傳算法、粒子群算法、灰狼算法等在波長選擇上已有很多成功研究案例。其中，蜻蜓算法（dragonfly algorithm，DA）是Mirjalili等在2016年通過對自然界蜻蜓行為進(jìn)行觀察、總結(jié)和抽象后，提出的一種新智能群體優(yōu)化算法。之后，分別采用連續(xù)DA、二進(jìn)制DA（binary DA，BDA）及多目標(biāo)DA，通過對幾類典型函數(shù)進(jìn)行優(yōu)化驗證了算法的有效性。2019年，Chen Yuanyuan等首次把BDA應(yīng)用于近紅外特征波長篩選，并比較了單群BDA（single-BDA）、多群BDA（multi-BDA）、基于集成學(xué)習(xí)BDA（ensemble learning-based BDA，ELB-BDA），表明后兩種算法可以明顯提高波長選擇的穩(wěn)定性以及分析模型的泛化能力。然而，使用BDA進(jìn)行波長篩選最終得到的波長位置和波長數(shù)量具有強隨機性。盡管使用multi-BDA和ELB-BDA在一定程度上降低了波長數(shù)量并進(jìn)一步篩選出特征波長，提高了算法的穩(wěn)定性，但是這兩種策略操作復(fù)雜，需要多次進(jìn)行single-BDA運算，導(dǎo)致運行時間遠(yuǎn)多于single-BDA，不利于快速建模。并且由于multi-BDA和ELB-BDA的內(nèi)核還是single-BDA算法，選出的波長分布隨機性仍然較大。為了加快變量選擇速度，Yun Yonghuan等基于達(dá)爾文自然進(jìn)化論中簡單有效的“適者生存”原則，在變量組合集群分析的變量選擇策略中采用指數(shù)衰減函數(shù)不斷縮小搜尋空間以加快算法速度。借鑒這一思想，在single-BDA算法中引入指數(shù)衰減函數(shù)，使算法能在迭代過程中由快趨緩地剔除無用變量，有利于保留重要變量，但不利的是由于衰減曲線末端趨向平緩可能造成無效計算，或末端迭代刪除變量數(shù)大于1而漏掉某一組或多組變量組合中可能存在的最優(yōu)變量。為此，本研究嘗試用線性衰減函數(shù)代替末期的指數(shù)衰減函數(shù)，保證算法末期每迭代一次剔除一個變量。Single-BDA算法的改進(jìn)因此分為兩個階段：第1階段使用指數(shù)衰減函數(shù)對變量進(jìn)行快速挑選；第2階段采用線性衰減函數(shù)對變量進(jìn)行精細(xì)挑選。本研究把改進(jìn)算法命名為衰減消去BDA（attenuation elimination-BDA，AEBDA），以期能從小麥粉近紅外光譜中盡快挑選出數(shù)量少且穩(wěn)定的特征波長，并建立精度較高的蛋白質(zhì)近紅外分析模型。

1 材料與方法

1.1 材料

市購不同品牌以及不同批次的160 個小麥粉的樣品，置于保鮮袋中于室溫（20～23 ℃）保存。

1.2 儀器與設(shè)備

Micro NIR Pro 1700便攜式近紅外光譜儀 美國Viavi Solutions公司；D200杜馬斯定氮儀 海能未來科技有限公司。

1.3 方法

1.3.1 小麥粉蛋白質(zhì)含量的測定

參照GB 5009.5ü2016《食品中蛋白質(zhì)的測定》中的燃燒法進(jìn)行測定。

1.3.2 近紅外光譜采集

對小麥粉樣品在室溫（20～23 ℃）環(huán)境下不做前處理進(jìn)行光譜掃描。光譜儀機身采用金屬試管架夾持固定，探頭向下垂直對準(zhǔn)深1 cm的圓盤樣品池，樣品池頂部與探頭底部相距1 cm。采集光譜時，小麥粉樣品鋪平樣品池，按120°的角度間隔采集，得到3 條不同檢測點的光譜，取平均作為該樣品的最終采集光譜。儀器的波長范圍為908～1 650 nm，光譜分辨率為6.25 nm，有125 個波長通道，并用與儀器配套的測控軟件MicroNIRTM Pro v2.5采集和儲存信號。

1.3.3 預(yù)處理方法

常見的預(yù)處理方法有移動平滑法（moving average filter，MAF）、卷積平滑法（Savitzky-Golay filter，SGF）、一階導(dǎo)數(shù)（1st D）、標(biāo)準(zhǔn)正態(tài)變量變換（standard normal variate transformation，SNV）等。使用這4 種預(yù)處理方法的3 種組合（MAF與SNV、SGF與SNV、1st D與SNV，預(yù)處理窗口寬度為5）以求達(dá)到最好的建模預(yù)測效果。

1.3.4 建模與模型評估方法

建模方法采用偏最小二乘回歸（partial least square regression，PLSR）算法，交互驗證采用留一法；single-BDA和AE-BDA算法中采用的適應(yīng)度函數(shù)為PLSR回歸建模中交互驗證標(biāo)準(zhǔn)偏差（root mean square error of cross validation，RMSECV）。

使用RMSECV、預(yù)測標(biāo)準(zhǔn)偏差（root mean square errors of prediction，RMSEP）、決定系數(shù)（）對所建立的校正模型進(jìn)行評估，計算方法如下：

式（1）～（3）中：y為第個樣品化學(xué)值的測定值；y為預(yù)測過程中第個樣品的預(yù)測值；為所用樣品集樣品數(shù)；為化學(xué)值測定值的均值。

在建模預(yù)測過程中，越接近1，表示模型的回歸或者預(yù)測結(jié)果越好；如果為負(fù)值，表示擬合效果極差。RMSECV與RMSEP的值越小表示所建立的模型的穩(wěn)定性與精確度越好。

1.3.5 波長變量選擇方法

1.3.5.1 Single-BDA波長選擇

蜻蜓通常通過5 種主要策略改變它們的位置：避撞、結(jié)隊、聚集、覓食、躲避。應(yīng)用到波長選擇方法中，食物即最優(yōu)解，敵人即極差解，通過不斷改變位置即迭代更新運算直至到達(dá)食物的位置即求得最優(yōu)解。據(jù)此建立的位置更新計算方法如表1所示。

表1 蜻蜓5 種主要位置更新策略的數(shù)學(xué)建模Table 1 Mathematical modeling of five major position shift strategies of dragonfly

采用single-BDA算法進(jìn)行波長選擇時，位置向量X的每個元素值只能是0或1，因此蜻蜓位置更新并非如表1所示的直接在原有位置X后加上更新速度值Δ，而是只能在0和1之間切換。這就需要對Δ進(jìn)行連續(xù)域到離散域的轉(zhuǎn)換。最簡單、最有效的方法是采用傳遞函數(shù)（式（4）），以Δ作為函數(shù)輸入，返回一個[0,1]之間的數(shù)值，表示位置變化的概率。

計算所有個體位置的變化概率后，應(yīng)用式（5）更新蜻蜓在空間中的搜索位置：

式（5）中：為0到1之間的隨機數(shù)，負(fù)號表示邏輯取反運算。

1.3.5.2 AE-BDA波長選擇

AE-BDA算法是經(jīng)single-BDA算法結(jié)合指數(shù)和線性衰減函數(shù)改進(jìn)而來，其在波長篩選過程中主要分為前后兩個部分：快速挑選和精細(xì)挑選，兩部分流程如圖1所示。設(shè)定最優(yōu)波長組合比例，AE-BDA算法在迭代運行過程中，統(tǒng)計每次迭代single-BDA計算完畢產(chǎn)生的所有個體中前h個最優(yōu)RMSECV個體中波長出現(xiàn)的次數(shù)。出現(xiàn)的頻率越高，表示波長越重要。然后，按指數(shù)或線性衰減函數(shù)所計算的當(dāng)前AE-BDA迭代應(yīng)保留的變量數(shù)，確定需要剔除的無用波長。AE-BDA算法采用均等抽樣法（equal sampling method，ESM），產(chǎn)生算法迭代中single-BDA第一代初始二進(jìn)制矩陣，保證每個波長具有相同的被選中機會。

圖1 AE-BDA流程圖Fig.1 Flow chart of AE-BDA

波長快速挑選部分：算法開始，在全波段中無用波長數(shù)所占比例高，隨著迭代的進(jìn)行其所占比例逐步降低，根據(jù)指數(shù)衰減函數(shù)曲線特點，波長剔除速度由快到緩。在這一步中，設(shè)置種群數(shù)中個體數(shù)為，原始變量數(shù)為，快速篩選部分篩選變化率，算法運行結(jié)束保留的波長個數(shù)為p（1為快速挑選部分迭代次數(shù)），如果當(dāng)前迭代變量數(shù)hr＜p，則快速挑選部分結(jié)束；通過不斷的迭代更新，波長數(shù)量呈指數(shù)級減少，達(dá)到所設(shè)置保留的波長個數(shù)p為止。每次迭代所保留的波長數(shù)為：

式（6）中：為當(dāng)前迭代次數(shù)；為波長篩選率；p為第次迭代后所剩的理論波長數(shù)。

精細(xì)挑選部分：根據(jù)線性衰減函數(shù)曲線特點，保證在算法末期每代可剔除一個變量。設(shè)置最終波長數(shù)p（2為精細(xì)挑選部分迭代次數(shù)），篩選過程中每次迭代并剔除出現(xiàn)次數(shù)最少的1 個波長，利用當(dāng)次迭代得到的波長組合建立PLSR模型，記錄校正集的RMSECV。每次迭代所保留的波長數(shù)為：

式（7）中：為當(dāng)前迭代次數(shù)。

當(dāng)波長數(shù)減少到最小設(shè)定值p時，迭代終止。選取全部RMSECV記錄中最小的波長組合作為AE-BDA波長優(yōu)選結(jié)果。

圖2中，曲線1為AE-BDA算法波長數(shù)量的衰減趨勢圖（以初始變量數(shù)＝500 個，p＝30 個，最終波長數(shù)p＝10 個為例）。在快速挑選部分，受指數(shù)函數(shù)性質(zhì)影響，波長數(shù)量經(jīng)過1次迭代快速降低到p個，然后再由精細(xì)挑選進(jìn)一步壓縮波長數(shù)量至最小波長數(shù)p個；曲線2為指數(shù)衰減函數(shù)變量衰減曲線，曲線1相比較于曲線2在算法末端（p至p階段）逐個波長剔除，不會導(dǎo)致某一波長點數(shù)被跳過。圖2的柱狀圖為每一次迭代剔除的無用變量數(shù)量，在快速挑選部分，隨著迭代次數(shù)的增加每次迭代剔除的無用變量數(shù)也呈指數(shù)減少；在精細(xì)挑選部分，每一代消去1 個無用變量，保證無用變量被精細(xì)均勻剔除。

圖2 AE-BDA算法波長數(shù)衰減圖Fig.2 Decreasing number of variables in AE-BDA

1.4 數(shù)據(jù)處理

使用MATLAB 2016a軟件以及實驗室自主開發(fā)的NIRSA 5.9.4化學(xué)計量學(xué)軟件，采用主成分分析和馬氏距離相結(jié)合的方法檢測異常樣本，主成分分析結(jié)合Kennard-Stone的方法對校正集和預(yù)測集進(jìn)行劃分。

2 結(jié)果與分析

2.1 小麥粉近紅外光譜圖

圖3為小麥粉近紅外光譜部分原始光譜，共有160 條小麥粉近紅外光譜，每條光譜含125 個波長通道。

圖3 小麥粉近紅外光譜部分原始光譜Fig.3 Selected original near-infrared spectra of wheat flour

2.2 小麥粉蛋白質(zhì)含量

由表2可知，所選樣本的蛋白質(zhì)含量基本覆蓋小麥粉中蛋白質(zhì)質(zhì)量分?jǐn)?shù)（7%～15%）且分布較為均勻，表明實驗所選樣本具有代表性。

表2 小麥粉蛋白質(zhì)含量統(tǒng)計表Table 2 Statistics of protein content of wheat flour

2.3 樣本異常值剔除與樣本集劃分

樣品、采集環(huán)境和儀器在一定程度上會產(chǎn)生異常樣本數(shù)據(jù)，嚴(yán)重影響所建模型的穩(wěn)定性與預(yù)測能力，所以在建模前必須將異常樣本數(shù)據(jù)從集合中剔除。3 種預(yù)處理組合各自檢測出的異常樣本集的并集（共8 個樣本）作為總異常樣本集予以剔除。最后，劃分剩余的152 個樣本為校正集和預(yù)測集；共得到98 個校正集樣本，54 個預(yù)測集樣本。如表3所示，校正集與預(yù)測集的樣本分布較寬，具有良好的代表性。

表3 校正集與預(yù)測集小麥粉蛋白質(zhì)含量統(tǒng)計表Table 3 Statistics of protein content of wheat flour in calibration and prediction sets

2.4 預(yù)處理方法對全譜段PLSR建模及預(yù)測的影響

由表4可知，MAF與SNV預(yù)處理組合的PLSR模型除

表4 不同預(yù)處理方法全譜PLSR模型性能Table 4 Performance of full-band PLSR models developed using different preprocessing methods

值略低于SGF與SNV預(yù)處理組合的PLSR模型外，各指標(biāo)均優(yōu)于其他建模方法。這說明采用MAF和SNV組合的預(yù)處理方法建立PLSR全譜校正模型可以提高小麥粉近紅外全譜段模型穩(wěn)定性和預(yù)測精度。因此，后續(xù)波長選擇算法研究中光譜都采用MAF與SNV組合進(jìn)行預(yù)處理。

2.5 2 種算法對波長選擇建模及預(yù)測的影響

分別用single-BDA和AE-BDA 2 種方法進(jìn)行波長篩選后建模，16 次建模及預(yù)測實驗結(jié)果的比較列于表5。實驗S1～S8波長選擇方法為single-BDA，實驗A1～A8波長選擇方法為AE-BDA。其中，AE-BDA算法設(shè)置種群數(shù)為200 個，衰減率為0.9，p設(shè)置為20 個，p設(shè)置為10 個，每次迭代中內(nèi)核single-BDA迭代次數(shù)為10 次。經(jīng)過程序運行后AE-BDA共迭代28 次（快速挑選迭代18 次，精細(xì)挑選迭代10 次），其內(nèi)核single-BDA共迭代280 次。為了便于改進(jìn)前后方法對比，single-BDA設(shè)置種群數(shù)量為200 個，迭代次數(shù)為280 次，篩選率為20%。

由表4和表5可知，無論是single-BDA還是AE-BDA算法，挑選后的波長進(jìn)行建模及預(yù)測效果均明顯優(yōu)于全譜段建模。使用single-BDA進(jìn)行波長選擇時，實驗S1～S8所挑選出的平均波長數(shù)為28.38 個，占原始波長數(shù)（125 個）的22.7%；而使用AE-BDA方法進(jìn)行波長選擇時，實驗A1～A8所挑選出平均波長數(shù)為15.8 個，為原始波長數(shù)的12.6%，相較于single-BDA方法篩選出建模波長數(shù)量更少。同時可以看到，AE-BDA方法挑選出的波長建模和預(yù)測評價指標(biāo)相應(yīng)都較接近，且其預(yù)測效果總體上好于single-BDA方法。盡管single-BDA方法建立的模型評價指標(biāo)及RMSECV總體上比AE-BDA好，但預(yù)測效果卻相對略差，這說明single-BDA方法選擇的波長過多，建立的模型可能存在過擬合現(xiàn)象，導(dǎo)致預(yù)測結(jié)果變差。

表5 Single-BDA與AE-BDA建模及預(yù)測性能比較Table 5 Comparison of modeling and prediction performance between single-BDA and AE-BDA

由圖4和圖5可知，盡管兩種方法都是采用隨機算法更新決策變量，但AE-BDA與single-BDA相比所挑選的波長點少，分布范圍更加集中，且AE-BDA中有3 個特征波長（1 385.07、1 508.95 nm和1 589.48 nm）在8 次實驗中均被挑選出來，而single-BDA方法挑選出頻次最多的2 個波長（908.1、1 502.76 nm）只有7 次。這說明AE-BDA挑選波長的隨機性要弱于single-BDA方法，結(jié)果更加穩(wěn)定。

圖4 Single-BDA（A）和AE-BDA（B）實驗入選波長位置Fig.4 Selected wavelength locations in single-BDA (A) and AE-BDA (B) experiments

圖5 Single-BDA（A）和AE-BDA（B）實驗中波長出現(xiàn)頻次圖Fig.5 Frequency of appearance of wavelengths in single-BDA (A) and AE-BDA (B) experiments

AE-BDA相較于single-BDA方法計算效率更高。由圖6可知，種群數(shù)相同（200 個）的情況下，隨著迭代次數(shù)的增加，AE-BDA參與迭代計算的總波長數(shù)呈指數(shù)下降，意味計算量也呈指數(shù)級降低，計算過程逐步加速；而single-BDA算法在迭代過程中，進(jìn)入每一次迭代的總波長數(shù)是不變的，計算量不變，運算速度保持相對恒定，運算時間較長。

圖6 Single-BDA和AE-BDA在迭代過程中初始波長數(shù)變化曲線Fig.6 Change in initial number of wavelength in single-BDA and AE-BDA as a function of number of iterations

雖然近紅外光譜帶重疊嚴(yán)重且很難清晰地指出各譜峰對應(yīng)的具體基團信息，但根據(jù)近紅外主要譜帶歸屬仍可對所選波長的合理性做出必要分析。本研究僅討論8 次AE-BDA實驗中每次都被挑選出的3 個波長點。蛋白質(zhì)分子中含有氨基，NüH鍵伸縮振動的一級倍頻吸收帶位于1 510 nm波長附近，1 508.95nm波長處為-氨基酸中的氨基的吸收峰；甲基和亞甲基中CüH鍵的伸縮振動的一級倍頻和變形振動的基頻吸收帶位于1 360～1 395nm波長處，1 385.07 nm波長處為CüH吸收峰；波長1 589.48 nm處位于OüH鍵的一級倍頻吸收譜區(qū)，由于模型中其標(biāo)準(zhǔn)回歸系數(shù)為負(fù)且較大，可以判斷為主要由非蛋白質(zhì)分子中OüH鍵吸收引起（波長點位置見圖7）。這可以作為對小麥粉中淀粉、脂肪等強背景組分的修正。通過進(jìn)一步對照分析可以看出，所選其他波長點大部分也位于以上譜帶區(qū)域內(nèi)或者附近。

圖7 出現(xiàn)8 次的波長點在全譜段的位置Fig.7 Locations of wavelengths that appeared eight times in the full-band spectrum

3 結(jié) 論

針對常規(guī)的single-BDA在近紅外光譜波長篩選過程中出現(xiàn)的問題，引入指數(shù)衰減函數(shù)和線性衰減函數(shù)對其進(jìn)行改進(jìn)，形成AE-BDA，并結(jié)合PLSR建立了小麥粉蛋白質(zhì)近紅外分析模型。通過對AE-BDA所挑選出的波長進(jìn)行特征分析表明，這些波長較好地分布于蛋白質(zhì)分子的主要官能團吸收區(qū)域。蛋白質(zhì)分子中重要的NüH鍵的吸收波長（1 508.95 nm）每次都能被挑選出來；其他反映蛋白質(zhì)分子的甲基、亞甲基以及次甲基的CüH鍵吸收波長也因出現(xiàn)頻次高而作為挑選出的特征波長用于建立小麥粉蛋白質(zhì)近紅外分析模型，提高了模型的可解釋性。

改進(jìn)后的AE-BDA運算速度快，所選波長穩(wěn)定，且一定程度上克服了過擬合，所建模型精度高。這為小麥粉蛋白質(zhì)近紅外光譜提供了一種可靠、有效的波長挑選方法，對于其他分析對象的近紅外光譜的變量選擇也有重要的借鑒意義。在后續(xù)研究中可以優(yōu)化算法中各因子的取值范圍，探討其對AE-BDA的影響，進(jìn)一步提高算法的運行效率和自動化水平。