郭培源，徐 盼，董小棟，許晶晶

（北京工商大學(xué)計(jì)算機(jī)與信息工程學(xué)院，食品安全大數(shù)據(jù)技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100048）

隨著人們生活水平的提高以及愈來愈多食品問題的發(fā)生，食品安全逐漸成為人們重點(diǎn)關(guān)注的問題，因而食品檢測(cè)的方法及其性能近年來成了研究熱點(diǎn)之一。香腸作為種類豐富且深受大家喜愛的食品，其檢測(cè)研究具有重要意義。受加工過程中不達(dá)標(biāo)的衛(wèi)生條件以及運(yùn)輸貯藏過程中環(huán)境因素的影響，香腸品質(zhì)會(huì)有所下降，而帶來很多安全問題和隱患[1]。評(píng)價(jià)香腸品質(zhì)常用的指標(biāo)是揮發(fā)性鹽基氮含量、蛋白質(zhì)分解產(chǎn)生、過氧化值和酸價(jià)含量[1-2]，以及亞硝酸鹽含量，少量的添加可給香腸上色且提供獨(dú)特的風(fēng)味，過量使用，則會(huì)對(duì)身體造成危害。除此之外，香腸的菌落總數(shù)也是一項(xiàng)重要的評(píng)價(jià)指標(biāo)，因?yàn)榧?xì)菌會(huì)加速香腸的腐敗，從而給人體的健康帶來不利影響[3-4]。

香腸中菌落總數(shù)的測(cè)定，傳統(tǒng)是采用理化實(shí)驗(yàn)的方法，即通過培養(yǎng)皿計(jì)數(shù)獲得，但是用理化實(shí)驗(yàn)獲取香腸菌落總數(shù)周期長(zhǎng)、耗試劑、操作繁瑣，且對(duì)樣品具有破壞性[5]。而近年來興起的高光譜成像技術(shù)是一種無損檢測(cè)技術(shù)，與理化實(shí)驗(yàn)相比，它具有在線實(shí)時(shí)、對(duì)樣品無破壞性、準(zhǔn)確便捷等優(yōu)點(diǎn)，現(xiàn)已廣泛應(yīng)用于食品檢測(cè)領(lǐng)域[6-7]。王莉等[8]采用波長(zhǎng)范圍400～1 000 nm可見近紅外高光譜對(duì)冷鮮羊肉的菌落總數(shù)和揮發(fā)性鹽基氮含量進(jìn)行新鮮度的檢測(cè)研究，其中，采用偏最小二乘回歸得到了最佳預(yù)測(cè)結(jié)果，建立的菌落總數(shù)和揮發(fā)性鹽基氮含量預(yù)測(cè)模型的校正集相關(guān)系數(shù)分別為0.906 7和0.914 7，預(yù)測(cè)集相關(guān)系數(shù)分別為0.874 3和0.880 2。劉善梅等[9]采用高光譜成像技術(shù)對(duì)生鮮豬肉的含水率進(jìn)行無損檢測(cè)，建立偏最小二乘回歸預(yù)測(cè)模型，交叉驗(yàn)證和預(yù)測(cè)相關(guān)系數(shù)分別為0.926和0.924，均方根誤差（root mean square error，RMSE）分別為0.467%和0.438%。張雷蕾等[10]在470～1 000 nm波長(zhǎng)范圍內(nèi)，從高光譜圖像中提取反射光譜，對(duì)光譜進(jìn)行多元散射校正（multiplicative scatter correction，MSC）處理，并采用偏最小二乘建模分析，實(shí)現(xiàn)對(duì)豬肉的新鮮度評(píng)價(jià)。Jin Huali等[11]利用偏最小二乘方法分別在400～1 000 nm全波段上和1 000～2 500 nm中選取的6 個(gè)特征波長(zhǎng)上進(jìn)行建模預(yù)測(cè)花生油中的成分含量，兩種方法的效果都很好，但是后者的效果優(yōu)于前者。Xiong Zhenjie等[12]采用偏最小二乘-連續(xù)投影算法的方法實(shí)現(xiàn)了紅肉中色素含量的定量檢測(cè)，并采用圖像處理的方法將色素在紅肉中的分布進(jìn)行可視化研究。

雖然高光譜成像技術(shù)已廣泛應(yīng)用于食品檢測(cè)領(lǐng)域，但利用高光譜技術(shù)檢測(cè)香腸內(nèi)化學(xué)物質(zhì)的含量以及對(duì)香腸進(jìn)行分級(jí)的相關(guān)研究與應(yīng)用非常少。本實(shí)驗(yàn)采用400～1 000 nm高光譜儀采集香腸的高光譜數(shù)據(jù)，并分別采用迭代決策樹（gradient boosting decision tree，GBDT）和支持向量回歸（support vector regression，SVR）方法建立香腸菌落總數(shù)的預(yù)測(cè)模型，以期為香腸菌落總數(shù)的快速定量和品質(zhì)控制提供參考。

1 材料與方法

1.1 材料與儀器

廣式香腸，購于北京永輝超市，將香腸切塊，每塊香腸厚2 cm，獲取50 份樣品，其中每份樣品有200 g，每天取一份樣本，將樣品平放于電移臺(tái)上，采用“推掃式”成像的方法獲取香腸的光譜值，然后進(jìn)行菌落數(shù)測(cè)定。

肉制品光譜檢測(cè)儀購自北京卓立漢光公司，波段范圍400～1 000 nm，采樣分辨率5 nm，共有128 個(gè)波段。高光譜成像系統(tǒng)硬件由裝有圖像采集卡的計(jì)算機(jī)、CCD相機(jī)、成像光譜儀、光源等系統(tǒng)組成。

1.2 高光譜檢測(cè)原理

高光譜成像技術(shù)既可以獲取含有物質(zhì)內(nèi)外理化信息的光譜值，同時(shí)也能通過成像設(shè)備獲取樣品各個(gè)波段的圖像數(shù)據(jù)，這種圖譜合一的三維數(shù)據(jù)稱為“數(shù)據(jù)立方體”[13]，如圖1所示。其中，圖像代表兩維的空間維度，而波長(zhǎng)代表一維的光譜維度?！皵?shù)據(jù)立方體”中每個(gè)波段可獲取一幅二維圖像，樣品的每個(gè)像素可以獲取一條光譜曲線[14-15]。圖像信息能夠全面反映物體的外在特征，而光譜信息則能夠反映物體的內(nèi)在物理結(jié)構(gòu)和化學(xué)成分等信息[16]。

圖1 高光譜成像技術(shù)檢測(cè)原理Fig. 1 Detection principle of hyperspectral imaging technology

1.3 方法

1.3.1 高光譜數(shù)據(jù)采集

由于高光譜圖像采集系統(tǒng)獲得的原始高光譜圖像在各個(gè)波段范圍內(nèi)的光源強(qiáng)度、光源亮度分布不均勻，并且暗電流和噪聲等因素會(huì)對(duì)光譜信息產(chǎn)生影響[17]，因而需要對(duì)采集到的高光譜圖像進(jìn)行黑白板校正處理[18]，得到樣品的光譜反射值，具體如式（1）所示：

式中：R為校正后圖像；IR為原始圖像；ID為黑板校正圖像；Iw為白板校正圖像。

使用高光譜分析處理軟件ENVI5.1，在每個(gè)樣本的高光譜圖像上選取感興趣區(qū)域（region of interest，ROI），對(duì)ROI采用N維可視化工具獲取平均光譜曲線[19]，如圖2所示。對(duì)50 個(gè)樣本中每個(gè)樣本選取10 個(gè)ROI，共得到500 個(gè)光譜數(shù)據(jù)。

圖2 香腸樣本的ROI（a）及其平均光譜曲線（b）Fig. 2 Region of interest in sausage samples (a) and its average spectral curve (b)

1.3.2 光譜預(yù)處理

在采樣過程中，由于樣品的不均勻性、高頻隨機(jī)噪聲、基線漂移、光散射等因素會(huì)對(duì)建模效果產(chǎn)生負(fù)面影響[20]，所以為了減少或消除此類因素的影響，需對(duì)采集的原始高光譜數(shù)據(jù)進(jìn)行不同的預(yù)處理，本研究采用MSC的預(yù)處理方法。MSC是高光譜建模最常用的預(yù)處理方法，分析結(jié)果較好[21]。它可以有效地消除樣品顆粒分布不均勻或者樣品大小不同等情況造成的散射誤差。

首先計(jì)算樣品得到的所有高光譜的平均光譜，將得到的平均光譜作為基準(zhǔn)光譜。每個(gè)光譜與基準(zhǔn)光譜進(jìn)行一元線性回歸運(yùn)算，求得各光譜相對(duì)于基準(zhǔn)光譜的線性平移量（回歸常數(shù)）和傾斜偏移量（回歸系數(shù)）。在每個(gè)原始光譜中減去回歸常數(shù)且除以回歸系數(shù)后，每個(gè)光譜的基線平移和偏移都得到了修正，而樣品成分含量對(duì)應(yīng)的光譜信息在數(shù)據(jù)處理的過程中沒有受到影響，進(jìn)而提高原始光譜的信噪比。平均光譜、回歸方程和MSC運(yùn)算的算法過程如（2）～（4）：

式中：Ai,j為香腸樣品的平均光譜曲線；Ai（MSC）為經(jīng)多元散射校正后的光譜；A為n×p維定標(biāo)光譜數(shù)據(jù)矩陣；n為光譜數(shù)量；p為波長(zhǎng)點(diǎn)數(shù)；A為原始光譜的平均矢量；Ai為1×p維矩陣，表示每個(gè)光譜矢量；mi為Ai與A線性回歸得到的相對(duì)偏移系數(shù)，Bi為Ai與A線性回歸獲得的平移量。

1.3.3 主成分分析

高光譜數(shù)據(jù)的數(shù)據(jù)量較為龐大，且相鄰波段的圖像相互重疊，具有很大的關(guān)聯(lián)性，因此高光譜數(shù)據(jù)降維處理的效果將影響后續(xù)的建模效果[22]，而主成分分析（principal component analysis，PCA）在數(shù)據(jù)降維方面具有獨(dú)特的優(yōu)勢(shì)，所得的主成分分量之間相互獨(dú)立，可以有效地消除高光譜數(shù)據(jù)中的冗余信息[23]。一般情況下，PC1包含波段中80%的方差信息，前3 個(gè)主成分包含了所有波段中90%以上的信息量[24]。本實(shí)驗(yàn)對(duì)光譜的128 個(gè)波段進(jìn)行PCA，選取前5 個(gè)主成分建立模型，這5 個(gè)主成分的方差累計(jì)貢獻(xiàn)率達(dá)到95%以上。

1.3.4 SVR

SVR是將復(fù)雜實(shí)際問題通過非線性變換轉(zhuǎn)換到高維特征空間，在高維空間中構(gòu)造線性決策函數(shù)來實(shí)現(xiàn)原空間中的非線性決策[24-27]。它是找到一個(gè)超平面，使到超平面最遠(yuǎn)的樣本點(diǎn)的“距離”最小。本研究選用SVR方法實(shí)現(xiàn)香腸菌落總數(shù)的預(yù)測(cè)，所用的高斯核函數(shù)[28]如式（5）所示：

式中：σ為核函數(shù)的寬度參數(shù)；x、xi分別為超平面最遠(yuǎn)的樣本點(diǎn)及中心點(diǎn)。

高斯核函數(shù)對(duì)數(shù)據(jù)中存在的噪聲有著較好的抗干擾能力，由于其具有很強(qiáng)的局部性，其參數(shù)決定了函數(shù)作用范圍，隨著參數(shù)σ的增大而減弱。

SVR高斯核函數(shù)的c和g參數(shù)，常用粒子群（particle swarm optimization，PSO）算法、網(wǎng)格搜索、遺傳算法（genetic algorithm，GA）3 種方法進(jìn)行尋優(yōu)。其中，c為懲罰系數(shù)，即對(duì)誤差的寬容度。c越高，說明越不能容忍出現(xiàn)誤差，容易出現(xiàn)過擬合；c越小，容易欠擬合。c過大或過小，泛化能力都差。g為選擇高斯核函數(shù)后自帶的一個(gè)參數(shù)，隱含地決定了數(shù)據(jù)映射到新的特征空間后的分布。g越大，支持向量越少；g越小，支持向量越多，支持向量的個(gè)數(shù)影響訓(xùn)練與預(yù)測(cè)的速度。

1.3.5 迭代決策樹

迭代決策樹（gradient boosting decision tree，GBDT）作為回歸樹的一種，相對(duì)于一般決策樹算法具有防止過擬合、泛化能力較強(qiáng)等優(yōu)點(diǎn)。模型訓(xùn)練的時(shí)候，對(duì)于輸入的一個(gè)樣本，首先會(huì)賦予一個(gè)初值，然后會(huì)遍歷每一棵決策樹，每棵樹都會(huì)對(duì)預(yù)測(cè)值進(jìn)行調(diào)整修正，每一次迭代是為了改進(jìn)上一次結(jié)果，減少上一次模型的殘差，并且在殘差減少的梯度方向上建立新的組合模型[29]。其基本思想是通過構(gòu)建M個(gè)弱分類器，經(jīng)過多次迭代最終組合而成一個(gè)強(qiáng)分類器。

GBDT又被稱為提升樹，其可以表示為決策樹的加法模型：

式中：T（xi,θm）為決策樹；θm為決策樹的參數(shù)；M為樹的個(gè)數(shù)。

針對(duì)不同問題的GBDT，其主要區(qū)別在于使用的損失函數(shù)不同，包括用平方誤差損失函數(shù)的回歸問題，用指數(shù)損失函數(shù)的分類問題，以及用一般損失函數(shù)的一般決策問題。本研究使用平方誤差損失函數(shù)實(shí)現(xiàn)回歸。

提升樹的訓(xùn)練流程如下：

輸入：訓(xùn)練數(shù)據(jù)集T＝{（x1,y2）, （x2,y2）, …, （xN,輸出：提升樹fM（x）。

1）初始化f0（x）＝0；2）對(duì)m＝1, 2, …,M；a）按式（6）計(jì)算殘差rmi＝y(tǒng)i-fm-1（xi），i＝1, 2, …,N；b）擬合殘差rmi學(xué)習(xí)一個(gè)回歸樹，得到T（x,θm）；c）更新fm（x）＝fm-1（x）+T（x,θm）；3）得到回歸問題提升

1.3.6 模型評(píng)價(jià)指標(biāo)

1.3.6.1 RMSE

RMSE是觀測(cè)值與真值偏差的平方和與觀測(cè)次數(shù)n比值的平方根，它能很好地反映測(cè)量的精密度，RMSE越小，模型的預(yù)測(cè)效果越好。具體如式（7）所示：

式中：Xobs,i為觀測(cè)值；Xmodei,i為真實(shí)值。

1.3.6.2 決定系數(shù)R2

決定系數(shù)用于判斷回歸方程的擬合程度，R2越接近1，模型的預(yù)測(cè)效果越好。具體如式（8）所示：

式中：yi為真實(shí)值；y為均值；為估計(jì)值。

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

實(shí)驗(yàn)采用10折交叉驗(yàn)證的方法對(duì)原始數(shù)據(jù)進(jìn)行處理，即每次用9 個(gè)子集的并集作為訓(xùn)練集，余下的1 個(gè)子集作為測(cè)試集，這樣總共獲得10 組訓(xùn)練/測(cè)試集，從而進(jìn)行10 組訓(xùn)練和測(cè)試，最終訓(xùn)練與測(cè)試的結(jié)果返回的是10 組的均值。因而本實(shí)驗(yàn)中每次有450 個(gè)香腸光譜樣本作為訓(xùn)練集，50 個(gè)作為測(cè)試集。對(duì)經(jīng)MSC預(yù)處理和PCA降維處理后的光譜分別采用SVR和GBDT方法建立香腸菌落總數(shù)的預(yù)測(cè)模型，并驗(yàn)證模型的預(yù)測(cè)效果。采用RMSE和R2作為評(píng)價(jià)模型預(yù)測(cè)效果的指標(biāo)，獲得了較好的實(shí)驗(yàn)結(jié)果，其訓(xùn)練集和測(cè)試集均方根誤差分別為0.001和0.003，決定系數(shù)R2分別為0.998和0.996。

2 結(jié)果與分析

2.1 光譜預(yù)處理的結(jié)果

原始光譜以及經(jīng)MSC后的光譜圖如圖3所示。

圖3 原始（A）和經(jīng)MSC預(yù)處理后（B）的光譜圖Fig. 3 Original spectra (A) and spectra preprocessed by MSC (B)

2.2 SVR建模結(jié)果

表1 不同參數(shù)尋優(yōu)算法對(duì)應(yīng)的SVR建模結(jié)果Table 1 Modeling results of SVR with different parameter optimization algorithms

3 種尋優(yōu)算法得到的最優(yōu)c和g參數(shù)及其對(duì)應(yīng)模型的預(yù)測(cè)結(jié)果如表1所示。SVR方法采用PSO算法進(jìn)行c和g參數(shù)尋優(yōu)得到的預(yù)測(cè)模型最優(yōu)。PSO算法的收斂速度快，受問題維數(shù)的變化影響較小，使得求解過程更容易計(jì)算[30]。采用PSO算法進(jìn)行c和g參數(shù)尋優(yōu)的SVR方法對(duì)香腸光譜訓(xùn)練集和測(cè)試集樣本的菌落總數(shù)預(yù)測(cè)結(jié)果如圖4所示。

圖4 SVR模型的訓(xùn)練集（a）與測(cè)試集（b）的預(yù)測(cè)結(jié)果Fig. 4 Prediction results of training (a) and test (b) sets with SVR model

2.3 GBDT建模結(jié)果

迭代1 000、1 500 次和2 000 次的GBDT建模結(jié)果如表2所示。迭代2 000 次得到的建模結(jié)果最好，并且迭代過程很快。迭代2 000 次的GBDT方法對(duì)香腸光譜樣本的訓(xùn)練和測(cè)試結(jié)果如圖5所示。

表2 不同迭代次數(shù)的GBDT建模結(jié)果Table 2 Modeling results of GBDT with different iterations

圖5 GBDT模型的訓(xùn)練集（a）與測(cè)試集（b）的預(yù)測(cè)結(jié)果Fig. 5 Prediction results of training (a) and test (b) sets with GBDT model

2.4 兩種方法建模結(jié)果的比較

由SVR和GBDT的建模結(jié)果，比較采用PSO算法進(jìn)行參數(shù)尋優(yōu)的SVR建模結(jié)果與迭代2 000 次的GBDT建模結(jié)果可知，GBDT建模結(jié)果要遠(yuǎn)優(yōu)于SVR的，GBDT所得的RMSE非常小，比SVR所得的要小一個(gè)數(shù)量級(jí)，并且GBDT所得的R2幾乎為1。除此之外，SVR建模所需的訓(xùn)練時(shí)間很長(zhǎng)，GBDT訓(xùn)練時(shí)間則很短。因而基于高光譜成像技術(shù)利用GBDT方法預(yù)測(cè)香腸菌落總數(shù)的方法可行且有效。

3 結(jié) 論

本實(shí)驗(yàn)通過高光譜成像系統(tǒng)采集50 個(gè)香腸樣本的高光譜數(shù)據(jù)，并利用高光譜分析處理軟件ENVI5.1，在每個(gè)香腸樣本的高光譜圖像中選擇10 個(gè)ROI，從而獲得500 個(gè)香腸樣本的平均光譜數(shù)據(jù)。實(shí)驗(yàn)采用MSC方法對(duì)光譜預(yù)處理，并采用PCA方法從128 個(gè)光譜波段中選擇5 個(gè)特征波段，從而提高了模型的預(yù)測(cè)精度。以處理過的光譜數(shù)據(jù)作為輸入，理化實(shí)驗(yàn)所得的香腸菌落總數(shù)值作為輸出，分別采用SVR和GBDT方法建立香腸菌落總數(shù)的預(yù)測(cè)模型。實(shí)驗(yàn)結(jié)果可知，迭代2 000 次的GBDT建模結(jié)果最優(yōu)。本實(shí)驗(yàn)中，GBDT模型迭代2 000 次時(shí)，訓(xùn)練集和測(cè)試集的RMSE都很小，R2也都接近1，但是當(dāng)?shù)螖?shù)多于2 000 次時(shí)，是否會(huì)產(chǎn)生過擬合、建模效果需要進(jìn)一步論證。除此之外，進(jìn)一步需要探索研究地是，將每個(gè)像素點(diǎn)下預(yù)測(cè)出的菌落總數(shù)定量反演到香腸樣本表面圖像上，生成可視化分布圖，使香腸新鮮度的動(dòng)態(tài)變化趨勢(shì)更加直觀、形象地呈現(xiàn)。