王朝輝,趙層,趙倩,王艷輝,賴翰卿,王曉東,王靖會(huì)
(1.吉林農(nóng)業(yè)大學(xué)食品科學(xué)與工程學(xué)院,吉林長(zhǎng)春130118;2.廣東地球土壤研究院,廣東廣州510385;3.吉林省長(zhǎng)春市凈月開(kāi)發(fā)區(qū)福祉街道辦事處,吉林長(zhǎng)春130122;4.吉林農(nóng)業(yè)大學(xué)信息技術(shù)學(xué)院,吉林長(zhǎng)春130118)
中國(guó)稻米的產(chǎn)量占世界稻米產(chǎn)量的30%以上[1],而吉林省“梅河大米”是中國(guó)粳米地理標(biāo)志產(chǎn)品,其產(chǎn)地位于世界黃金糧食生產(chǎn)帶(北緯45°)。在實(shí)際生活中梅河大米的種類繁多,通常采用凱氏定氮法和分光光度法等化學(xué)類方法測(cè)定不同品種中大米的蛋白質(zhì)含量,但是這些傳統(tǒng)的化學(xué)方法不但對(duì)樣品本身具有破壞性,而且步驟繁瑣,檢測(cè)周期過(guò)長(zhǎng)。紅外光譜技術(shù)作為一種快速無(wú)損的檢測(cè)手段,已經(jīng)廣泛的應(yīng)用于大米主要成分(蛋白質(zhì)[2]、脂肪[3]、淀粉[4]、水分[5])的檢測(cè),但是其只能根據(jù)光譜信息得到組分的含量,無(wú)法實(shí)現(xiàn)更加直觀的表達(dá)即含量的可視化。
高光譜是包括圖像信息和光譜信息[6]的三維立方體數(shù)據(jù),得到的高光譜圖像既包含了大米的內(nèi)部信息(內(nèi)部物理結(jié)構(gòu)、化學(xué)成分的信息)也包含了大米的外部信息(粒型、缺陷等),能彌補(bǔ)近紅外不能快速識(shí)別某一物質(zhì)的空間分布情況[7-8]即圖像的缺失。目前,高光譜成像技術(shù)在農(nóng)產(chǎn)品方面的主要應(yīng)用有瘀傷識(shí)別分類[9]、籽粒的識(shí)別[10]、品種的檢測(cè)[11]、含水量預(yù)測(cè)[12]等,而高光譜成像技術(shù)在大米中的主要應(yīng)用有堊白檢測(cè)[13]、水分含量檢測(cè)[14]、摻假大米的檢測(cè)[15]等。雖然高光譜檢測(cè)技術(shù)在組分的檢測(cè)方面已有諸多應(yīng)用,但關(guān)于高光譜成像技術(shù)應(yīng)用于大米中蛋白質(zhì)含量分布可視化的研究還未見(jiàn)相關(guān)報(bào)道。
本文選取了吉林省梅河市4 個(gè)產(chǎn)區(qū)的3 個(gè)品種(稻花香、秋田小町、吉粳60)的大米為研究對(duì)象,利用高光譜成像技術(shù)對(duì)采集的大米進(jìn)行檢測(cè),求取大米感興趣區(qū)域的平均光譜,為了降低光譜的信噪比和得到相對(duì)穩(wěn)健的模型,對(duì)光譜進(jìn)行了卷積平滑(savitzkygolay,S-G)、均值中心化(mean centering)、多元散射校正 (multiplicative scatter correction,MSC)3 種算法的預(yù)處理,建立了大米蛋白質(zhì)含量的偏最小二乘回歸(partial least-squares regression,PLSR)、主成分回歸(principal component regression,PCR)、誤差反向傳播神經(jīng)網(wǎng)絡(luò)(back propagation neural networks,BP-NN)3 種預(yù)測(cè)模型。采用SPA 選取特征波長(zhǎng),建立特征波長(zhǎng)模型,并將大米高光譜圖像轉(zhuǎn)變?yōu)榈鞍踪|(zhì)含量分布圖,實(shí)現(xiàn)不同品種產(chǎn)地大米蛋白質(zhì)含量的可視化。
以吉林省梅河市種植的稻米品種(稻花香、秋田小町、吉粳60)為研究對(duì)象,通過(guò)空間分層隨機(jī)布點(diǎn)的原則采集了包括稻花香、秋田小町、吉粳60 3 個(gè)品種在內(nèi)的120 個(gè)稻米樣品。試驗(yàn)前先將采集的稻米進(jìn)行脫殼、去糙,之后剔除霉變、未成熟和粒型不完整的籽粒,選擇大小均勻、飽滿的籽粒進(jìn)行高光譜掃描和進(jìn)行大米蛋白質(zhì)含量的化學(xué)測(cè)定。采集的120 個(gè)大米樣品中86 個(gè)用于校正集,34 個(gè)用于檢驗(yàn)集。
ImSpector V10E-QE 型高光譜圖像系統(tǒng):芬蘭Spectral Imaging 有限責(zé)任公司;9589(EKE-ER)全光譜鹵素?zé)簦汉商mPhilips 公司;PSA200-11-XZolix 電控位移臺(tái):北京Zolix 公司;C8484-05G 相機(jī):日本Hamam atsu Photonics;GZ02DS20 可升降樣品臺(tái):北京廣正儀器有限公司;JLGJ4.5 型檢驗(yàn)礱谷機(jī)、JNMJ3 型檢驗(yàn)?zāi)朊讬C(jī):臺(tái)州市糧儀廠;JXFM110 錘式旋風(fēng)磨:上海嘉定糧油儀器有限公司;K1100 全自動(dòng)凱氏定氮儀、SH220 石墨消解儀:濟(jì)南海能儀器股份有限公司;PHS-3C 型酸度計(jì):上海儀電科學(xué)儀器股份有限公司。
測(cè)量的光譜范圍是408.360 0 nm~1 007.220 0nm,包含477 個(gè)波段,系統(tǒng)采集大米樣品前通過(guò)對(duì)比多次預(yù)實(shí)驗(yàn)結(jié)果來(lái)調(diào)整曝光時(shí)間、物距等各項(xiàng)參數(shù)。采集圖像時(shí)物距為13.5 cm,曝光時(shí)間為15 ms,位移臺(tái)移動(dòng)速度為1.62 mm/s。設(shè)定好儀器參數(shù)后將放有大米樣品的黑色底板平鋪在載物臺(tái)上(黑色底板的反射率接近于0,大米樣本按照網(wǎng)格單粒隨機(jī)的方式擺放至載物臺(tái)黑板上),進(jìn)行大米高光譜數(shù)據(jù)的采集。高光譜成像系統(tǒng)的結(jié)構(gòu)示意圖如圖1 所示。
圖1 高光譜成像系統(tǒng)示意圖Fig.1 Schematic diagram of hyper-spectral imaging system
在高光譜圖像采集過(guò)程中由于會(huì)受到外界(自然光光源分布不均)和系統(tǒng)(暗電流噪聲)等因素的干擾[16],導(dǎo)致得到的高光譜圖像包含大量的無(wú)效信息,為了減少由于上述因素導(dǎo)致的噪聲增強(qiáng)現(xiàn)象,因此在采集大米高光譜圖像之前,必須先進(jìn)行原始大米圖像的黑白校正。校正公式如下
式中:Ra為原始大米的圖像;Rd為黑板圖像;Rb為白板圖像;Rn校正后大米圖像。
將采集完光譜數(shù)據(jù)的大米樣品進(jìn)行磨粉過(guò)篩,采用GB5009.5-2016《食品安全國(guó)家標(biāo)準(zhǔn)食品中蛋白質(zhì)的測(cè)定》測(cè)量大米中的蛋白質(zhì)含量,之后把稱量好的大米樣品和催化劑放入準(zhǔn)備好的消化管中(大米樣品重量 0.5g 左右,CuSO4和 K2SO4的重量分別是 0.2 g 和3 g),準(zhǔn)備完成后在消化管中倒入10 mL 的濃硫酸進(jìn)行消化,并做空白管試驗(yàn),當(dāng)消化管內(nèi)溶液的顏色變?yōu)榫G色透明時(shí)證明其消化完全,則可通過(guò)全自動(dòng)凱氏定氮儀對(duì)其進(jìn)行蛋白質(zhì)含量的測(cè)定,每個(gè)樣品測(cè)量3 次取其平均值作為該大米樣品蛋白質(zhì)含量的化學(xué)值。
1.4.1 感興趣區(qū)域與平均光譜的提取
在ENVI5.0 軟件中對(duì)獲得的高光譜圖像進(jìn)行感興趣區(qū)域的選取,選取3×15 網(wǎng)格內(nèi)所有單粒米樣中間部位11×11 像素的正方形區(qū)域作為1 個(gè)感興趣區(qū)域(region of interest,ROI),ROI 全部像素點(diǎn)的平均光譜作為建立模型的光譜數(shù)據(jù)。
1.4.2 光譜預(yù)處理
為了獲得相對(duì)穩(wěn)定和準(zhǔn)確的模型,從原始光譜中摒棄噪音等過(guò)多無(wú)用信息的干擾,從而放大差異,提高分辨率。采用卷積平滑(savitzky-golay,S-G)、均值中心化(mean centering)、多元散射校正(multiplicative scatter correction,MSC)3 種光譜預(yù)處理方法。通過(guò)3種預(yù)處理算法的對(duì)比,以期找到適合本次試驗(yàn)的光譜預(yù)處理方法。
1.4.3 模型的建立
采用偏最小二乘回歸(partial least-squares regression,PLSR)、主成分回歸(principal component regression,PCR)、誤差反向傳播神經(jīng)網(wǎng)絡(luò)(back propagation neural networks,BP-NN)3 種常用的算法建立蛋白質(zhì)的預(yù)測(cè)模型,采用不同建模算法,最終呈現(xiàn)的建模效果也各不相同。本文以大米為試驗(yàn)對(duì)象,以校正集和檢驗(yàn)集樣品的相關(guān)系數(shù)(Rc2、Rp2)和均方根誤差(RMSEC、RMSEP)為評(píng)價(jià)模型性能的指標(biāo),通過(guò)對(duì)比最終得到穩(wěn)健性好的大米蛋白質(zhì)預(yù)測(cè)模型。
1.4.4 偏最小二乘回歸模型
偏最小二乘回歸 (partial least square regression,PLSR)是多因變量對(duì)多自變量的回歸模型,是常用的一種光譜建模分析方法[17]。PLSR 通過(guò)對(duì)多自變量和多因變量進(jìn)行壓縮分解,分解時(shí)得到兩者最強(qiáng)對(duì)應(yīng)關(guān)系的主成分,加強(qiáng)模型的預(yù)測(cè)能力。
1.4.5 主成分回歸模型
主成分回歸(principal component regression,PCR)是一種多元線性回歸分析方法,利用主成分分析將樣品光譜矩陣進(jìn)行分解,把得到的變量組合成相互無(wú)關(guān)的新變量[18],用較少的變量代表原有的變量,能有效解決預(yù)測(cè)變量高度相關(guān)和共線的問(wèn)題。
1.4.6 誤差反向傳播神經(jīng)網(wǎng)絡(luò)模型
誤差反向傳播神經(jīng)網(wǎng)絡(luò)(back propagation neural networks,BP-NN)是一種按誤差逆?zhèn)鞑ニ惴ㄓ?xùn)練的多層前饋網(wǎng)絡(luò)[19],是目前應(yīng)用最廣泛的神經(jīng)網(wǎng)絡(luò)模型之一。本文選用有導(dǎo)師學(xué)習(xí)算法來(lái)調(diào)整神經(jīng)元間的聯(lián)接權(quán),使得網(wǎng)絡(luò)輸出更符合實(shí)際。
1.4.7 連續(xù)投影算法(successive projection algorithm,SPA)特征波長(zhǎng)的選取
連續(xù)投影算法(successive projection algorithm,SPA)是一種使矢量空間共線性最小化的前向變量選擇算法[20],采用SPA 對(duì)原始高光譜數(shù)據(jù)進(jìn)行特征波長(zhǎng)的選取,解決了高光譜數(shù)據(jù)自身維度高,變量多等問(wèn)題,由于后期模型的建立使用消除了大量冗余信息的高光譜數(shù)據(jù),因此可以根據(jù)挑選的特征波長(zhǎng)(保留光譜信息中最優(yōu)的變量組)建立簡(jiǎn)化高效的模型。
1.4.8 大米中蛋白質(zhì)含量的可視化
通過(guò)對(duì)比得到最適合本研究的大米蛋白質(zhì)含量的預(yù)測(cè)模型,將像素點(diǎn)的光譜數(shù)據(jù)帶入預(yù)測(cè)模型,可得到各像素點(diǎn)位置的蛋白質(zhì)含量。通過(guò)各像素點(diǎn)位置蛋白質(zhì)含量的高低,使用MATLAB 對(duì)灰度圖像進(jìn)行各像素點(diǎn)顏色的設(shè)定。光譜預(yù)處理、預(yù)測(cè)模型的建立和大米蛋白質(zhì)可視化的實(shí)現(xiàn)均在MATLAB2016a 軟件中完成。
圖2 是選取ROI 全部像素點(diǎn)的平均值得到的反射率光譜(408.360 0 nm~1 007.220 0 nm),120 個(gè)大米全波長(zhǎng)光譜的趨勢(shì)大致相同,說(shuō)明大米在全波長(zhǎng)的情況下具有相同的反射特性,但在反射強(qiáng)度上存在差異。
圖2 120 個(gè)大米樣本光譜圖Fig.2 Spectral map of 120 rice samples
對(duì)凱氏定氮法測(cè)量的大米蛋白質(zhì)的化學(xué)值進(jìn)行統(tǒng)計(jì)分析,得到校正集、檢驗(yàn)集和全部樣品的最小值、最大值、平均值和標(biāo)準(zhǔn)偏差,結(jié)果如表1 所示。
表1 校正集和檢驗(yàn)集樣品蛋白質(zhì)含量Table 1 Calibration set and test set sample protein content
從表1 中可以看出校正集蛋白質(zhì)含量范圍大于檢驗(yàn)集的蛋白質(zhì)含量范圍,說(shuō)明校正集和檢驗(yàn)集的劃分合理。
對(duì)120 個(gè)大米樣品的蛋白質(zhì)含量進(jìn)行統(tǒng)計(jì)分析得到如圖3 所示的柱形圖。
圖3 不同品種和產(chǎn)地的蛋白含量柱狀圖Fig.3 Bar graph of protein content of different varieties and regions
如圖所示3 個(gè)品種的大米在吉樂(lè)鄉(xiāng)產(chǎn)區(qū)的蛋白質(zhì)含量明顯高于其他產(chǎn)區(qū),不同品種大米的蛋白質(zhì)含量在產(chǎn)區(qū)間的分布規(guī)律也不相同,稻花香從曙光到黑山頭蛋白質(zhì)含量表現(xiàn)為依次遞增,而秋田小町和吉粳60從曙光、灣龍?jiān)俚胶谏筋^蛋白質(zhì)含量表現(xiàn)為先增加后減少。對(duì)品種和產(chǎn)地進(jìn)行雙因素方差分析,得到關(guān)于品種因素、產(chǎn)地因素和品種產(chǎn)地的交互對(duì)應(yīng)的P 值,因?yàn)樗萌叩腜 值均遠(yuǎn)遠(yuǎn)小于0.01,因此品種和產(chǎn)地對(duì)大米蛋白質(zhì)含量均有顯著影響。
通過(guò)對(duì)比不同預(yù)處理算法建立的PLSR 模型的相關(guān)系數(shù)和均方根誤差,得到最優(yōu)的光譜預(yù)處理算法。表2 是光譜預(yù)處理的結(jié)果,通過(guò)對(duì)比可以發(fā)現(xiàn)MC 預(yù)處理后的光譜所建立模型最優(yōu),其校正集和檢驗(yàn)集的相關(guān)系數(shù)為0.927 5 和0.907 9,SG 和MSC 預(yù)處理后所得模型檢驗(yàn)集相關(guān)系數(shù)下降為0.807 0 和0.897 1。
表2 基于不同預(yù)處理方法的PLSR 模型結(jié)果Table 2 Results of PLSR models based on different preprocessing methods
本文以477 個(gè)波段的光譜數(shù)據(jù)為自變量,以測(cè)量的蛋白質(zhì)含量為因變量,建立PLSR、PCR、BP 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型,基于全波段不同建模算法的結(jié)果見(jiàn)表3。
表3 基于全波段不同建模算法的結(jié)果Table 3 Results based on different modeling algorithms for the full band
如表3 所示,建立PLSR 模型時(shí)確定的最優(yōu)潛在變量個(gè)數(shù)是13,其校正集的相關(guān)系數(shù)和均方根誤差RMSEC 分別為0.927 5 和0.034 6,其檢驗(yàn)集的相關(guān)系數(shù)和均方根誤差RMSEP 分別為0.907 9 和0.050 8;建立PCR 模型時(shí)確定的最優(yōu)主成分個(gè)數(shù)是17,其校正集的相關(guān)系數(shù)和均方根誤差RMSEC 分別為0.920 0和0.037 9,其檢驗(yàn)集的相關(guān)系數(shù)和均方根誤差RMSEP 分別為0.897 1 和0.055 7;建立BP 神經(jīng)網(wǎng)絡(luò)模型時(shí),神經(jīng)網(wǎng)絡(luò)隱含層神經(jīng)元數(shù)20,其校正集的相關(guān)系數(shù)和均方根誤差RMSEC 分別為0.968 5 和0.015 7,其檢驗(yàn)集的相關(guān)系數(shù)和均方根誤差RMSEP 分別為0.862 5 和0.074 5。通過(guò)對(duì)3 種建模算法的結(jié)果對(duì)比得到最優(yōu)的建模算法為PLSR。
蛋白質(zhì)PLSR 模型的實(shí)測(cè)值和預(yù)測(cè)值的散點(diǎn)圖見(jiàn)圖4。
圖4 蛋白質(zhì)PLSR 模型的實(shí)測(cè)值和預(yù)測(cè)值的散點(diǎn)圖Fig.4 Scatter plot of measured and predicted values of protein PLSR model
特征波長(zhǎng)的計(jì)算過(guò)程在MATLAB 中實(shí)現(xiàn),設(shè)定最小和最大的選定波長(zhǎng)數(shù)為5 到50,最終得到的結(jié)果如圖5 所示。
圖5 SPA 篩選特征波長(zhǎng)Fig.5 SPA screening characteristic wavelength
通過(guò)SPA 篩選出27 個(gè)特征波長(zhǎng),分別為:416.690 0、429.840 0、426.250 0、459.910 0、469.590 0、519.590 0、861.890 0、881.170 0、892.750 0、910.760 0、927.480 0、939.060 0、959.650 0、966.080 0、971.230 0、975.080 0、978.940 0、982.800 0、990.510 0、993.090 0、995.660 0、1 000.800 0、1 003.370 0、948.070 0、747.900 0、487.810 0、1 005.940 0 nm。
PLSR 通過(guò)計(jì)算分別得到合適的X 和Y 的潛變量P 和Q,通過(guò)迭代法調(diào)整得分矩陣T 和U,同時(shí)保持殘差矩陣E 和F 的值接近0,求得T 和U 的內(nèi)在聯(lián)系,間接得到X 和Y 的映射關(guān)系[21]。
基于SPA 篩選的27 個(gè)特征波長(zhǎng)的反射值作為自變量X,全自動(dòng)凱氏定氮儀測(cè)定的大米蛋白質(zhì)含量為因變量Y 建立PLSR 模型,得到27 個(gè)特征波長(zhǎng)的反射值和蛋白質(zhì)含量的內(nèi)在關(guān)系,并用此關(guān)系在已知波譜反射值的情況下對(duì)蛋白質(zhì)含量值進(jìn)行預(yù)測(cè)。圖6 所示SPA-PLSR 模型的校正集和檢驗(yàn)集的相關(guān)系數(shù)分別為0.923 2 和 0.903 9。
圖6 蛋白質(zhì)SPA-PLSR 模型的實(shí)測(cè)值和預(yù)測(cè)值的散點(diǎn)圖Fig.6 Scatter plot of measured and predicted values of the protein SPA-PLSR model
選取稻花香、秋田小町和吉粳60 的大米高光譜圖像進(jìn)行成像分析,通過(guò)ENVI5.0 從大米樣本的高光譜圖像中提取27 個(gè)特征波長(zhǎng)下的高光譜圖像,然后提取特征圖像中每個(gè)像素點(diǎn)的光譜反射值,將提取的所有像素點(diǎn)的光譜反射值代入2.5 已建立的PLSR 模型,最終預(yù)測(cè)到每個(gè)像素點(diǎn)的蛋白質(zhì)含量。將高光譜灰度圖像轉(zhuǎn)入MATLAB2016a 中,高光譜灰度圖像中所有像素點(diǎn)的灰度值用得到的對(duì)應(yīng)像素點(diǎn)的蛋白質(zhì)含量值代替,生成如圖7 所示的大米蛋白質(zhì)含量分布圖。
圖7 不同品種產(chǎn)地蛋白質(zhì)含量分布可視化圖Fig.7 Visualization of protein content distribution in different varieties of production areas
如圖7 所示3 個(gè)品種在吉樂(lè)鄉(xiāng)產(chǎn)區(qū)的蛋白質(zhì)含量值最高,在曙光、灣龍和黑山頭產(chǎn)區(qū)蛋白質(zhì)含量無(wú)明顯差異,因?yàn)榧獦?lè)鄉(xiāng)采集的大米地處山區(qū)與其他3 個(gè)產(chǎn)區(qū)相比距離最遠(yuǎn),地貌差異最大,因此形成的大米蛋白質(zhì)含量差異明顯,表明產(chǎn)地對(duì)大米蛋白質(zhì)的含量是有影響的,與圖3 所示的規(guī)律基本相同。
從全部大米蛋白質(zhì)含量分布圖可以看出米粒內(nèi)部大米蛋白質(zhì)的分布是不均勻的,但是總分布趨勢(shì)基本相同,胚內(nèi)多于胚乳,大米蛋白質(zhì)含量的高低主要表現(xiàn)為胚乳內(nèi)蛋白質(zhì)含量分布情況,這與Little and Dawson(1960)公布的大米蛋白質(zhì)染色圖得到的結(jié)論基本相同。
采用高光譜成像技術(shù)對(duì)大米中蛋白質(zhì)含量分布的可視化進(jìn)行了可行性研究,通過(guò)MC 光譜預(yù)處理方法和SPA 特征波段的選取,得到了簡(jiǎn)化高效的PLSR蛋白質(zhì)含量預(yù)測(cè)模型,基于該定量模型的基礎(chǔ)上對(duì)不同品種和不同產(chǎn)地大米中蛋白質(zhì)含量的分布進(jìn)行可視化研究。由于同一品種間大米形狀相近,難以通過(guò)普通RGB 圖像區(qū)分,通過(guò)對(duì)蛋白質(zhì)含量分布成像可以為大米產(chǎn)地的識(shí)別提供思路,而對(duì)比不同品種間大米的蛋白質(zhì)含量分布圖,可以為后期選育大米品種提供依據(jù)。