韓仲志 劉 杰

（青島農(nóng)業(yè)大學(xué)理學(xué)與信息科學(xué)學(xué)院 山東青島266109）

黃曲霉毒素（Aflatoxin）是一種劇毒、強(qiáng)致癌物，其毒性為砒霜的68 倍，是目前發(fā)現(xiàn)的最強(qiáng)的I類化學(xué)致癌物質(zhì)。其廣泛存在于花生、玉米籽粒表面及其制品中。中美兩國(guó)都對(duì)其進(jìn)行了強(qiáng)制限量：食品級(jí)和飼料級(jí)的限量標(biāo)準(zhǔn)分別為20 μg/kg 和100 μg/kg[1-2]。目前對(duì)黃曲霉毒素的檢測(cè)主要是生化方法，包括薄層層析法、高效液相色譜法、微柱法、酶聯(lián)免疫吸附法等[3]。雖然檢測(cè)精度高，但檢測(cè)手段和檢測(cè)儀器復(fù)雜，檢測(cè)速度慢，價(jià)格昂貴，不能在線檢測(cè)。

近年來，光譜成像技術(shù)作為一項(xiàng)新的化學(xué)計(jì)量學(xué)手段，廣泛應(yīng)用于農(nóng)產(chǎn)品及食品檢測(cè)中[4-5]。黃曲霉毒素具有紫外熒光特性和表面淺表面分布特性。高光譜技術(shù)為一種圖譜合一的新興手段，在黃曲霉毒素檢測(cè)中受到普遍關(guān)注。M.Atas 等[6]研究了市場(chǎng)上購(gòu)買的辣椒粉黃曲霉毒素污染的高光譜成像分析技術(shù)，提出基于神經(jīng)網(wǎng)絡(luò)權(quán)值優(yōu)化的特征選擇方法，并找到對(duì)檢測(cè)起關(guān)鍵作用的特征波長(zhǎng)，在420 nm 處留一法識(shí)別率達(dá)到85%；H.Yao[7]團(tuán)隊(duì)主要集中在玉米的黃曲霉毒素污染的高光譜研究，他們通過人工種植的方法發(fā)現(xiàn)熒光強(qiáng)度隨著黃曲霉素含量的增加而減少，并且存在峰移現(xiàn)象，還基于此兩點(diǎn)發(fā)現(xiàn)申請(qǐng)了專利[8]，然而這種設(shè)備的檢測(cè)還是手工進(jìn)行，并不能做到自動(dòng)檢測(cè)。王偉[9]研究組使用美國(guó)農(nóng)業(yè)部（USDA）的高光譜數(shù)據(jù)研究黃曲霉毒素的檢測(cè)問題，指出根據(jù)普通CCD 成像，黃曲霉毒素自然污染籽粒的檢出率為87.5%，他們采用Fisher 優(yōu)選了700～1 100 nm的5 個(gè)高光譜波長(zhǎng)，對(duì)人工污染黃曲霉毒素濃度的預(yù)測(cè)正確率達(dá)到88.3%[10]。

上述研究忽略了一個(gè)關(guān)鍵因素，黃曲霉毒素在籽粒表面分布并不均勻，這是由于黃曲霉毒素是黃曲霉菌的代謝物，呈顆粒狀不均勻分布。通過光學(xué)探測(cè)器得到的像元往往是特定分辨率下各種物質(zhì)平均像元，得到的光譜是各種物質(zhì)（包括黃曲霉毒素）的混合光譜。單個(gè)像素內(nèi)黃曲霉毒素豐度是定量計(jì)算整個(gè)籽粒含量的重要基礎(chǔ)，有必要對(duì)單像素內(nèi)的黃曲霉毒素豐度進(jìn)行解析。本研究擬借用遙感領(lǐng)域像元分解的方法，探討花生籽粒表面黃曲霉毒素的豐度，進(jìn)而對(duì)黃曲霉毒素進(jìn)行定量反演。

1 試驗(yàn)材料

1.1 樣品制備

試驗(yàn)所用花生樣品為市場(chǎng)上購(gòu)買4 粒紅小花生，經(jīng)青島市海潤(rùn)農(nóng)大檢測(cè)中心檢測(cè)未檢出黃曲霉毒素，挑選1 g 左右表面光滑250 籽粒備用。在青島農(nóng)業(yè)大學(xué)食品學(xué)院實(shí)驗(yàn)室對(duì)購(gòu)買的黃曲霉毒素B1標(biāo)準(zhǔn)樣品，通過乙氰配比20，50，100，500 和1 000 μg/kg 共5 種黃曲霉毒素溶液，然后使用移液器將1 μL 的黃曲霉毒素溶液滴到花生表面，自然風(fēng)干。如表1，共制備4 粒紅花生5 個(gè)濃度，每個(gè)濃度50 粒，共250 ?；ㄉ瑢?duì)這250 個(gè)籽粒以此編號(hào)為1～250，正面朝上，放置在無(wú)熒光紙板上，冷藏備用。

表1 試驗(yàn)材料濃度與編號(hào)Table 1 Concentration and number of experimental materials

1.2 圖像采集

圖像采集在中國(guó)科學(xué)院青島光電所光譜實(shí)驗(yàn)室的暗室環(huán)境下進(jìn)行。圖1（a）列出了試驗(yàn)中使用的采集設(shè)備，試驗(yàn)儀器為便攜式成像儀，該成像儀主要由液晶可調(diào)式濾波器（Liquid Crystal Tunable Filter，LCTF）和一臺(tái)普通CCD 相機(jī)組成，譜段范圍為400～720 nm，光譜帶寬為10 nm，使用一臺(tái)UV 365 nm 大功率LED 紫外燈 （美國(guó)陸陽(yáng)LUYOR-3404 臺(tái)式紫外燈，樣品處照度7 000 流明）作為光源。高光譜圖像的分析采用ENVI4.7（ITT Visual Information Solutions，Boulder，Colo.），Matlab2012b（Math Works，USA）完成。計(jì)算機(jī)配置為：CPU 為Intel-E4600，主頻2.39 GHz，內(nèi)存為DDR3-3.25 GB。高光譜成像儀配套采集軟件可根據(jù)圖像整體亮度情況自動(dòng)調(diào)整每個(gè)波長(zhǎng)下的曝光時(shí)間從而避免過飽和和欠飽和現(xiàn)象。

紫外燈下，將花生籽粒污染面朝上每9 個(gè)一組，每次采集9 顆花生籽粒，依此采完250 個(gè)籽粒的圖像，每次獲得的高光譜圖像數(shù)據(jù)塊大小為1 392×1 040×33。圖1（a）是基于RGB 圖像合成原理，通過33，5，16 三個(gè)波段（720，440，550 nm）合成的假彩色圖像。

圖1 圖像采集設(shè)備圖Fig.1 Image acquisition equipment

圖2 假彩色合成圖像Fig.2 Composited color image

2 研究方法

2.1 算法

高光譜圖像像素的混合過程通常通過線性混合模型 （linear spectral random mixture model，LSRMM）[11]來描述，高光譜某個(gè)混合像元的輻射值：

其中，ei，j為第j個(gè)純像元的光譜曲線，j=1，2，…m，m 為純像元個(gè)數(shù)；其中i=1，2，…L，L 為波段數(shù)；aj為第j 個(gè)像元的豐度，aj滿足非負(fù)約束和為一約束。

根據(jù)線性混合模型，需要確定混合像元的端元波譜，N-FINDR 端元提取算法[12]是一種從已獲得的高光譜圖像自動(dòng)化的獲取端元的方法。高光譜的所有像素在高維空間形成一個(gè)凸體，每個(gè)端元對(duì)應(yīng)于凸體的頂點(diǎn)，非端元?jiǎng)t分布在凸體的內(nèi)部、棱或面上，該算法的迭代思想是，對(duì)經(jīng)PCA 或MNF 降維后的數(shù)據(jù)，隨機(jī)選擇一組像素作為初始端元，然后將圖像中每個(gè)像素替換現(xiàn)有的端元向量中的端元計(jì)算替換后的體積，若體積增大則用新的端元向量替換原端元向量。重復(fù)遍歷所有像元，最終的端元向量即是該圖像的端元。

著名的科學(xué)雜志《Nature》刊登了兩位科學(xué)家D.D.Lee 和H.S.Seung[13]對(duì)數(shù)學(xué)中非負(fù)矩陣研究的突出成果。該文提出了一種新的矩陣分解思想——非負(fù)矩陣分解（Non-negative Matrix Factorization，NMF）算法，即NMF 是在矩陣中所有元素均為非負(fù)數(shù)約束條件之下的矩陣分解方法。

豐度圖像的灰度分布直方圖可能提供一些有價(jià)值的信息。可用于黃曲霉素含量的預(yù)測(cè)中，對(duì)一副灰度圖像的灰度量化特征表示為：

式中：n=1，2，…，N 表示灰度值量化為灰度階的數(shù)量。I（x，y）——坐標(biāo)（x，y）處灰度值。將這些直方圖量化特征作為支持向量機(jī)回歸的輸入特征。

支持向量機(jī)（Support Vector Machine，SVM）是Corinna Cortes 和Vapnik 等[14]首先提出的，它在解決小樣本、非線性及高維模式識(shí)別中表現(xiàn)出許多特有的優(yōu)勢(shì)，并能夠推廣應(yīng)用到函數(shù)擬合等其它機(jī)器學(xué)習(xí)問題中。支持向量機(jī)算法可描述為：

式中：αk——拉格朗日乘子；K（x，xk）——核函數(shù)，本研究選自RBF 核函數(shù)；b--偏置系數(shù)。

2.2 研究框架

黃曲霉毒素在花生表面分布呈微小顆粒狀不均勻分布，通過光學(xué)探測(cè)器得到的像元是一定區(qū)域（一個(gè)像素）內(nèi)各種物質(zhì)平均像元。將單個(gè)像元進(jìn)行亞像元分解，分解為花生種皮光譜、黃曲霉素光譜和背景光譜各個(gè)組分的豐度，通過黃曲霉素組分豐度圖像可更為精確地反演黃曲霉含量。圖3是本文提出的方法流程圖。

圖3 方法流程圖Fig.3 Method flow chart

對(duì)通過高光譜儀獲得的高光譜圖像，首先進(jìn)行圖像預(yù)處理和ROI 提取，將單個(gè)花生籽粒高光譜圖像塊提取出來，然后通過N-FINDR 算法，提取花生種皮、黃曲霉毒素及背景的端元光譜，反復(fù)比較提取較為純凈的3 個(gè)端元光譜，并觀察3 個(gè)端元光譜波形，確定黃曲霉素的光譜，然后使用非負(fù)矩陣分解（NMF）方法對(duì)籽粒高光譜圖像進(jìn)行非負(fù)矩陣分解，得到黃曲霉毒素豐度圖像，對(duì)豐度圖像求取直方圖量化特征，使用該特征進(jìn)行回歸預(yù)測(cè)輸入?yún)?shù)，將該豐度與實(shí)際數(shù)據(jù)進(jìn)行支持向量機(jī)回歸，為比較預(yù)測(cè)效果，與偏最小二乘回歸預(yù)測(cè)效果進(jìn)行比較，進(jìn)行精度評(píng)價(jià)和檢測(cè)方法的有效性。

3 結(jié)果與分析

3.1 圖像預(yù)處理

背景區(qū)域?qū)D像分析造成干擾，需要去除，采用RGB 圖像分量運(yùn)算（R-1.2G）去除背景，并將背景置零；根據(jù)籽粒區(qū)域平均亮度一致原則將籽粒區(qū)域進(jìn)行光照補(bǔ)償，然后標(biāo)記籽粒區(qū)域，將籽粒區(qū)域定義為感興趣區(qū)域（ROI），根據(jù)籽粒左上角坐標(biāo)與籽粒區(qū)域長(zhǎng)寬，在33 個(gè)波段上將籽粒提取出來，順序存儲(chǔ)，將單籽粒高光譜數(shù)據(jù)存為.mat 文件，對(duì)應(yīng)的偽彩色圖像存為.jpg 文件，共得到250個(gè)單個(gè)籽粒（ROI）高光譜圖像塊，預(yù)處理過程通過Matlab 編程自動(dòng)實(shí)現(xiàn)。圖4（a）是分割得到其中一個(gè)單籽粒高光譜圖像塊ROI 對(duì)應(yīng)的假彩色圖像，可見圖像上黃曲霉素區(qū)域的熒光現(xiàn)象比較明顯。

3.2 N-FINDR 純端元提取

端元提取就是在高光譜圖像上找到純像元（端元），圖4（b，c）是基于N-FINDR 對(duì)圖4（a）提取的端元所在位置圖和端元曲線。圖4（b）綠色圈點(diǎn)表示尋找到的端元的位置。圖4（c）從端元曲線來看端元1 在400～420 波段內(nèi)有一明顯的波峰，這于黃曲霉素的熒光曲線相符合，可見端元1 即為黃曲霉素端元，端元3 為純黑色很顯然是背景的端元，因?yàn)轭A(yù)處理是已將背景置為0，那么端元2 則為花生種皮光譜端元?？梢钥闯?，基本將3 種端元尋找出來。根據(jù)經(jīng)驗(yàn)，多個(gè)籽粒反復(fù)比較得到黃曲霉毒素、花生種皮和背景的純像元光譜。

圖4 端元及其光譜曲線Fig.4 Endmember and spectrum

由于高光譜數(shù)據(jù)量巨大，在端元提取之前首先要進(jìn)行最大噪聲分離（MNF）[15]以減少數(shù)據(jù)量，這樣可以用前幾個(gè)波段來代替整個(gè)高光譜影像，可以有效加快算法的執(zhí)行效率。圖5是最大噪聲分離后的信噪比分布，可以看出，前面MNF 波段信噪比高，5 個(gè)以后信噪比已經(jīng)低于1，前2 個(gè)MNF 分量的分布相當(dāng)集中（圖5b）。

圖5 MNF 前兩個(gè)端元分布與等高線圖Fig.5 First two end-members distribution and contour map

3.3 NMF 豐度提取

使用非負(fù)矩陣分解（NMF）進(jìn)行亞像元分解，就是得到各個(gè)端元的豐度，圖6是各個(gè)端元的豐度，可以明顯看出，圖6（a）是黃曲霉毒素的豐度，圖6（b）是花生種皮的豐度，圖6（c）是背景的豐度，通過黃曲霉素豐度圖（圖6a）像素值求和即可計(jì)算出每個(gè)籽粒黃曲霉素的豐度和。

圖6 各端元的豐度圖Fig.6 Endmember abundance image

對(duì)黃曲霉毒素豐度圖像進(jìn)行求和。圖7（a）是基于MNF 進(jìn)行豐度分解后得到的250 個(gè)籽粒黃曲霉毒素的豐度和分布，單從圖像上看并無(wú)規(guī)律。可見僅通過豐度求和較難準(zhǔn)確預(yù)測(cè)黃曲霉毒素含量。豐度概率分布可能蘊(yùn)含著大量的信息，對(duì)豐度圖像進(jìn)行直方圖量化是體現(xiàn)豐度概率分布的有效手段。

3.4 直方圖量化

直方圖量化就是將豐度圖像（0～1）轉(zhuǎn)化為灰度圖像（0～255），然后將灰度值直方圖分成不同的區(qū)間（如0～50，50～100，100～150 等），統(tǒng)計(jì)灰度值落入每個(gè)區(qū)間的像素個(gè)數(shù)。圖7（b）是對(duì)圖6（a）黃曲霉毒素豐度圖的直方圖量化為12 個(gè)區(qū)間得到的直方圖量化結(jié)果，通過量化可以發(fā)現(xiàn)，灰度主要集中在低灰度區(qū)域，這與整副圖像偏黑有關(guān)。各個(gè)區(qū)域直方圖分布代表了不同豐度黃曲霉毒素的分布，這些分布信息與黃曲霉毒素含量直接相關(guān)，可將該量化特征（12 個(gè)）作為回歸模型的輸入特征。

圖7 黃曲霉毒素豐度統(tǒng)計(jì)及直方圖量化Fig.7 Aflatoxin abundance and quantification histograms

3.5 支持向量機(jī)回歸

將直方圖量化特征作為黃曲霉毒素回歸模型的輸入?yún)?shù)，進(jìn)行含量預(yù)測(cè)，與黃曲霉毒素的實(shí)際含量比較得到預(yù)測(cè)誤差。圖8是比較了支持向量機(jī)回歸（SVR）與傳統(tǒng)的偏最小二乘回歸（PLS）得到的誤差圖，前200 個(gè)樣本用來訓(xùn)練，后5 個(gè)樣本用來測(cè)試，可以看出SVR 效果較PLS 效果好，且訓(xùn)練集誤差小于測(cè)試集誤差。

表2是采用五折交叉驗(yàn)證法得到的兩種模型的預(yù)測(cè)結(jié)果，可以看出SVR 訓(xùn)練集誤差只有0.89%。平均相對(duì)誤差為12.16%，由于黃曲霉素限量很低，只有20～100 μg/kg，1 μg/kg 相當(dāng)于1 t 里的1 g[15]，能夠達(dá)到這樣一個(gè)精度，可以在很大程度上將高黃曲霉毒素含量的籽粒鑒別出來，這是實(shí)際生產(chǎn)應(yīng)用可以接受的。

圖8 SVR 和PLS 兩種方法預(yù)測(cè)誤差Fig.8 Predicted error of SVR and PLS methods

表2 兩種模型的預(yù)測(cè)結(jié)果Table 2 Predicted result of the two models

4 結(jié)論

環(huán)境中的灰塵、背景板在一定程度上都具有熒光現(xiàn)象，在高光譜偽彩色合成圖像上表現(xiàn)為高亮黃綠熒光區(qū)域，這些因素的存在對(duì)黃曲霉素的準(zhǔn)確預(yù)測(cè)造成了干擾，實(shí)用的技術(shù)中需要對(duì)這些因素進(jìn)行預(yù)先的去除，比如加裝必要的除塵設(shè)備和配置沒有熒光現(xiàn)象的背景板，特氟龍材質(zhì)的背景板是一種很好的選擇。

基于亞像元分解的方法被廣泛應(yīng)用遙感影像的定量解析中，本研究基于MNF 的亞像元分解，提出了一種能在食品安全檢測(cè)領(lǐng)域應(yīng)用的黃曲霉豐度檢測(cè)方法，并構(gòu)建了基于豐度圖像的直方圖量化特征進(jìn)行非線性回歸，進(jìn)而預(yù)測(cè)黃曲霉素精度，該方法最優(yōu)總體預(yù)測(cè)誤差以降低到12.16%。由于黃曲霉毒素含量在mg/kg 數(shù)量級(jí)，也就是百萬(wàn)分之一數(shù)量級(jí)，這一精度，已能夠較為準(zhǔn)確地識(shí)別花生是否被黃曲霉毒素污染，進(jìn)而將污染籽粒剔除。相比較生化的檢測(cè)方法，雖然預(yù)測(cè)精度較低，但該方法預(yù)測(cè)速度快，可以在瞬間完成，對(duì)黃曲霉的在線快速檢測(cè)及其相關(guān)便攜式裝備的研發(fā)具有積極意義。