基于光譜成像的豬肉新鮮度空間分布預測評價方法

趙茂程 吳澤本 汪希偉,2 邢曉陽 陳加新 唐于維一

(1.南京林業大學機械電子工程學院， 南京 210037；2.南京林業大學機電產品包裝生物質材料國家地方聯合工程研究中心， 南京 210037)

0 引言

豬肉新鮮度的評價指標包括感官、揮發性鹽基氮(Total volatile basic nitrogen，TVB-N)含量和微生物含量等指標[1-2]，本文選擇TVB-N含量客觀評價豬肉新鮮度。

高光譜成像技術在農林產品質量安全檢測等方面獲得了較多應用[3-7]。光譜圖像中每個像素點均包含特定信息，因此可以實現指標空間分布預測的可視化[8-17]，即基于興趣區域的均值光譜建立化學計量學預測模型，并將其應用于像素光譜中，從而得到指標空間分布預測圖像。這一方法已被應用于農林產品品質檢測工作中[18]。但像素位置缺乏微觀指標參考值，無法對可視化圖像進行直接評價，因此目前的研究僅停留在建立可視化圖像層面，未對指標空間分布預測的效果進行評價。

然而，將眾多準度相近的化學計量學模型應用到像素光譜進行指標空間分布預測時，像素預測結果存在明顯差異，甚至出現預測結果不具備統計意義的情況[19-20]。究其原因，其差異主要表現在兩方面：一是準度，即各像素位置微觀預測值的統計均值與理化檢測值的偏差程度；二是精度，即根據指標理論允許范圍，異常點在興趣區域內所占比值。因此，本文以豬肉新鮮度指標空間分布預測為例，建立基于不同光譜預處理及化學計量學模型所得新鮮度指標空間分布預測圖像，從準度和精度兩方面，對指標空間分布預測評價的方法進行研究。

1 材料與方法

1.1 實驗材料

實驗所用豬肉樣本購買于江蘇省南京市麥德龍下關商場，分割出5條豬肉背長肌作為研究對象，用內置冰袋的保溫箱運送至南京林業大學逸夫科技實驗樓。將每條豬肉背長肌分割成18塊厚度約為10 mm的肉塊，共獲取90個實驗樣本。分割完成后，將實驗樣本放置于4℃的冰箱中冷藏保存，以維持豬肉內部的自然組織結構狀態。

1.2 高光譜圖像采集及處理

1.2.1高光譜成像系統

高光譜成像系統由光譜成像單元、照明系統、計算機和輔助支架組成，其中光譜成像單元包括工作在可見-近紅外波段(550～1 000 nm)的CVA-200型聲光可調諧濾波器(BRIMEROSE公司，美國)、ORCA-R2型可見-近紅外相機(HAMAMATSU公司，日本)和可變焦鏡頭(18～200 mm，尼康公司，日本)；照明系統包括C3K型UPS穩壓電源(山特電子(深圳)有限公司)輸出的12個50 W鹵素燈和1個半球形漫反射穹頂。由計算機(Yangtian A4600t， Windows XP系統)通過基于Visual Studio 2008平臺編寫的采集程序實現高光譜圖像采集控制。

1.2.2高光譜圖像采集

每隔24 h從冰箱中取出10個實驗樣本，用保溫箱運送至高光譜成像系統處采集光譜數據。在550～1 000 nm波段范圍內，區域灰度均值范圍設定為1 600～2 400，采用變曝光時間采集方式，相機曝光時間如圖1所示。采集步長為3 nm，曝光物距為860 mm，圖像分辨率為516像素×672像素，單次采集實驗樣本的高光譜圖像。

圖1 曝光時間Fig.1 Exposure time

由于970～1 000 nm波段內圖像信噪比過低，故舍棄該部分圖像，僅保留550～970 nm共141個波段的光譜圖像。

在豬肉樣本高光譜圖像采集結束后，先采集標稱值75%的標準反射率標定板(SRT-75-100型，Labsphere公司，美國)的光譜圖像；后蓋上鏡頭蓋，采集此時的光譜圖像，即暗噪圖像。將以上兩種光譜圖像作為當日采集豬肉樣本光譜圖像相對反射率校正的基準。

1.2.3高光譜圖像預處理

(1)光譜均值濾波

本文采用光譜均值濾波提高像素光譜質量。設置固定寬度的濾波器窗口，沿原始光譜滑動，不同窗口的跨度相互重疊，取窗口內k個波段的平均值代替第1個波段的數據。濾波后第i個波段的光譜數據計算方法為

(1)

式中Rout(i)——光譜均值濾波后第i個波段的光譜圖像

Rin(h)——第h個波段的原始光譜圖像

m0——原始光譜圖像的波段數量

本文采用的聲光可調諧濾波器每個波段的帶寬最大為20 nm。為使光譜均值濾波帶寬與光學濾鏡帶寬相當，本文設置濾波器寬度為6、18、30、42、54 nm，對應波段數量分別為2、6、10、14、18。

(2)相對反射率校正

為進一步提高信噪比，消除圖像采集過程中暗電流、背景光強度及光源分布不均勻等產生的噪聲影響，需要對采集的高光譜圖像進行相對反射率校正[21]。利用當日系統暗噪圖像和標稱值為75%的標準反射率標定板的光譜圖像，對采集的樣本光譜進行相對反射率校正，計算公式為

(2)

式中RT——相對反射率校正后的光譜圖像

Ro——相對反射率校正前的光譜圖像

Rb——當日系統暗噪圖像

Rw——標稱值為75%的標準反射率標定板光譜圖像

w——根據標準反射率，用于Rw的標準反射率標定板各波段實際反射率

(3)肌肉有效興趣區域創建

TVB-N含量指標主要反映肌肉中蛋白質的分解程度，因此將興趣區域定義為肌肉有效興趣區域(Eligible muscle region of interest，EMROI)，即僅保留肌肉區域，排除脂肪、筋膜及肉皮等組織成分的干擾。豬肉中肥肉和肌肉部分與背景在圖像特定通道上有不同的反射率，因此，首先采用固定閾值法在727 nm處人工設置反射率閾值為0.3，對高光譜數據集內所有圖像進行分割，實現豬肉與背景的分離；其次，由于部分豬肉樣本的側面也被采集至圖像中，用半徑為8像素的圓形模板對所有分離后豬肉圖像進行邊緣腐蝕操作，以消除對預測模型的影響，僅保留待檢的上表面區域；最后在805 nm處人工設置反射率閾值0.45對所有腐蝕后豬肉圖像進行分割，以剔除肥肉部分，實現肌肉的精確分離。

(4)特征光譜提取

通過求取不同波長下EMROI內所有空間像素位置處的光譜數據均值，提取得到豬肉樣本不同預處理方法下的特征光譜{si}，即均值光譜。計算公式為

(3)

式中si——特征光譜(均值光譜)在第i個波段的反射率

aij——EMROI內第j個像素點在第i個波段的反射率

n——EMROI內像素點數量

m——光譜圖像的波段數量

(5)特征波段選取

本文采集的高光譜圖像數據量較大，包含550～970 nm范圍內的141個波段信息，其中包含很多冗余信息，降低了模型運算效率。因此，需要采用合適的方法選取相關性較高的特征波段。采用連續投影算法(Successive projection algorithm，SPA)[22]對基于EMROI的全波段光譜進行特征波段提取，該算法是一種前向循環的使矢量空間共線性最小化特征變量篩選方法，可以篩選出有效信息，降低數據之間共線性影響，增加模型的魯棒性和泛化性。

為便于工業現場應用，實現樣本的快速檢測，常需要組建專用的多光譜成像系統。基于掃描式光譜成像，光譜波段數量一般不超過6個；基于快照式光譜成像，相機在一次曝光周期內同時采集圖像的波段數量通常為20個左右。因此，本文分別優選出6個和20個特征波段，用于建立豬肉新鮮度預測模型。

1.3 揮發性鹽基氮含量測定

參照對鮮凍肉的評價標準[1]，TVB-N含量(質量比)不超過15 mg/(100 g)可認定為新鮮肉。將采集完高光譜圖像的樣本用保溫箱運送至理化檢測實驗室，以測定TVB-N含量。參照對揮發性鹽基氮測定的國標方法[23]，取樣本表面5 mm的純肌肉部分，剁成肉粉，采用半微量定氮法進行測定。在同等條件下，對獲得的兩次獨立測定結果計算平均值。若兩次結果的絕對差值不超過平均值的10%，以該平均值作為該實驗樣本的TVB-N含量理化檢測值。

1.4 基于均值光譜的豬肉新鮮度預測模型建立與評價

偏最小二乘回歸(Partial least squares regression，PLSR)[24]是一種廣泛應用于光譜分析的線性回歸建模方法，是典型的相關分析和主成分分析的集成和發展，可以同時實現提取變量特征、分析變量間相關性和回歸建模，是光譜分析領域應用最廣泛的方法之一[25]。其思路是從自變量集合中選取主成分，然后建立主成分與自變量的回歸方程，公式為

(4)

式中Y——預測值

β0——增益常數

βi——第i個波段的增益系數

Xi——第i個波段的平均反射率

本文基于Matlab工具箱中偏最小二乘回歸函數plsregress建立預測模型，得到回歸參數向量β，作為化學計量學模型的波段增益。其中，βi值越大，對預測結果的貢獻度越高，相應波段的放大倍率也越高。

1.5 基于像素光譜的豬肉新鮮度空間分布預測評價方法

1.5.1豬肉新鮮度空間分布預測

將基于均值光譜的化學計量學模型{β0，β1,…，βm}應用到像素光譜進行新鮮度指標空間分布預測，計算出相應像素點的TVB-N含量預測值，得到灰度圖像。各像素位置微觀預測值的計算公式為

(5)

式中Pj——EMROI內第j個像素點的微觀預測值

根據數值大小結合偽色彩處理技術，繪制得到豬肉新鮮度可視化圖像。肉樣表面預測值由高(腐敗)至低(新鮮)分別以紅色至藍色顯示，背景設置成黑色，肥肉部分設置成白色，樣本原始尺寸(即包含樣本側面的未腐蝕圖像尺寸)在肉樣外側用白色線條表示，色帶范圍設置為0～25 mg/(100g)。

1.5.2空間分布預測評價方法

由于無法對各像素位置的微觀新鮮度TVB-N含量指標進行直接評價，因此從準度及精度兩方面，對基于不同光譜預處理及化學計量學模型所得肉樣檢測區域新鮮度指標空間分布預測結果進行評價。

(1)準度評價

(2)精度評價

求取不同波長下EMROI內所有空間像素位置處光譜數據均值得到均值光譜的過程，相當于經過一個尺寸為EMROI大小的空間均值濾波器降噪，可有效提高信噪比。在將基于均值光譜所建化學計量學模型應用到像素光譜時，由于未經過該空間均值濾波，像素光譜質量明顯低于區域均值光譜，導致模型的預測精度下降，從而造成豬肉新鮮度TVB-N含量預測值出現負值。根據新鮮度理化指標的理論允許范圍，以TVB-N含量微觀預測值小于零的像素點在興趣區域內所占比值作為精度評價指標，即異常預測點占比越小，空間預測精度越高。

1.5.3空間分布預測精度影響因素

(1)像素光譜質量

通過理想均值濾波對光譜圖像進行6～54 nm的5種不同帶寬的光譜濾波降噪，可以有效提高像素光譜質量。將暗噪圖像各波段EMROI內所有像素點亮度值的標準差，記為暗噪標準差，以傳感器輸出值(Digital number，DN)作為計量單位，從而量化圖像隨機噪聲抑制效果，以評價像素光譜質量。通常來講，暗噪標準差越小，圖像的噪聲水平越低，像素光譜質量越高。

(2)模型波段增益

對于一般的化學計量學模型，波段增益是光譜圖像每一個波段的放大倍率。本文應用PLSR建立預測模型，m個波段的增益系數分別為PLSR回歸系數{β1，β2，…，βm}。

1.6 數據處理軟件

光譜圖像預處理、PLSR預測模型建立、空間分布預測、數據分析及圖形繪制均采用Matlab R2017b軟件；數據匯總整理由Excel 2016軟件完成。

2 結果與分析

2.1 樣本集劃分

本文共采集90個樣本。圖像采集時，第21號樣本20%表面肌肉區域被不慎掉落的脂肪組織覆蓋，造成圖像光譜數據與理化檢測值不對應，被判為無效樣本排除，剩余有效樣本共89個。選擇合適的定標集方法劃分樣本可以有效提高模型的預測能力。本文按照TVB-N含量排序，采用留出法[26]確定樣本集，其中67個樣本為訓練集，22個樣本為預測集。如表1所示，訓練集、預測集和總樣本的TVB-N含量變化范圍為3.49～26.19 mg/(100 g)，樣本標簽值跨度范圍在豬肉TVB-N含量正常范圍內，并且數據集樣本之間差異明顯，具有廣泛性；平均值和標準差相近，保證訓練集與預測集處于相同的分布狀況。數據集具有較強的泛化性和一致性，為建立可靠的豬肉新鮮度預測模型提供了數據基礎。

表1 豬肉樣本訓練集與預測集TVB-N含量統計結果Tab.1 Statistical results of TVB-N value in training set and prediction set of pork samplesmg/(100 g)

2.2 光譜曲線分析

對550～970 nm內141個波段求取EMROI內所有空間像素位置處的光譜數據均值，可獲得豬肉樣本不同預處理方法下的光譜反射率曲線。圖2分別表示豬肉樣本原始光譜反射率曲線及濾波器寬度為6、18、30、42、54 nm的濾波后光譜反射率曲線。由圖2a可以看出，不同豬肉樣本之間的光譜曲線具有相似的趨勢。575 nm附近的吸收峰是脫氧肌紅蛋白與氧合肌紅蛋白的合頻，肌紅蛋白的存在狀態決定了鮮肉的肉色[27]。

圖2 89個豬肉樣本EMROI內光譜反射率曲線Fig.2 Spectral reflectance curves of 89 pork samples in EMROI

對比圖2a～2f發現，隨著濾波器寬度的增加，光譜均值濾波有效降低了圖像噪聲，抑制了原始光譜曲線的細小波動，使光譜曲線變得更加光滑，預測模型的穩健性得以提升；但曲線峰谷的數量及波動幅度出現下降，特別是575 nm附近的吸收峰，當濾波器寬度為54 nm時該吸收峰已經消失，說明噪聲抑制能力增強的同時可能會導致部分細節信息的丟失。

2.3 豬肉新鮮度預測模型的建立

采用偏最小二乘回歸法，分別基于原始光譜和濾波后5種不同帶寬的光譜對全波段、精選出的20、6個特征波段建立新鮮度預測模型。為了得到可靠穩定的模型，采用交叉驗證法對訓練集進行驗證，以減小訓練過程中的過擬合。設置交叉驗證次數為5，分別建立PLSR豬肉新鮮度預測模型，模型性能如表2所示。

表2 PLSR豬肉新鮮度預測模型準度Tab.2 Accuracy of PLSR pork freshness model

2.4 豬肉新鮮度空間分布預測評價

2.4.1豬肉新鮮度空間分布預測準度評價

圖3 均值光譜預測值與像素光譜預測值對照結果Fig.3 Comparison results of mean spectral prediction and pixel spectral prediction

從原理上對樣本基于均值光譜的預測值和各像素位置微觀預測值的統計均值相等的原因進行探討。偏最小二乘回歸建模過程中，利用線性化學計量學算法，樣本基于均值光譜的預測值公式為

(6)

將基于均值光譜所建模型應用到像素光譜進行新鮮度指標空間分布預測時，利用線性化學計量學算法，各像素位置微觀預測值的統計均值公式為

(7)

式中P2——各像素位置微觀預測值的統計均值

將P2化簡，可得

(8)

由式(8)可知，對于線性化學計量學算法，樣本基于均值光譜的預測值恒等于各像素位置微觀預測值的統計均值。

文獻[19]利用類似光譜成像系統進行過新鮮度分布預測。與本文不同，其所得樣本表面新鮮度像素級定量可視化預測均值與有效興趣區域特征光譜預測值之間存在明顯差異，而該差異是由在光譜預處理環節采用的標準正態化處理(Standard normal variate，SNV)算法所引入的非線性所致。

2.4.2豬肉新鮮度空間分布預測精度評價

(1)各種新鮮度分布可視化圖像感官對比

根據各像素點預測值結合偽色彩處理技術，繪制得到豬肉新鮮度可視化圖像。同一肉樣經不同光譜預處理及化學計量學預測模型所得的新鮮度指標空間分布可視化圖像存在明顯差異。整體來看，隨濾波器寬度的增加，大部分化學計量學模型的可視化效果呈現出逐漸提升的趨勢；但存在個別模型游離在趨勢之外，可視化效果不升反降。選取一塊肥瘦分明的典型豬肉樣本，利用不同的化學計量學模型對新鮮度空間分布預測，在同一維度上，即特征波段數量相等或同為全波段建模時，各選擇3個典型模型建立可視化圖像，如圖4所示。圖4a～4c中，圖4c的可視化效果最好；圖4d～4f中，圖4f的可視化效果最好；圖4g～4i中，圖4h的可視化效果最好。肉樣表面預測值由高(腐敗)至低(新鮮)分別以紅色至藍色顯示，背景設置成黑色，肥肉部分設置成白色，樣本原始尺寸(即包含樣本側面的未腐蝕圖像尺寸)在肉樣外側用白色線條表示，色帶范圍設置為0～25 mg/(100g)。顏色越接近紅色，預測值越高，豬肉越腐?。活伾浇咏{色，預測值越低，豬肉越新鮮。可視化圖像右側上方附以該模型的各波段暗噪標準差和化學計量學模型增益，暗噪標準差用藍色線條表示，各波段增益值用紅色線條表示；右側下方附以肉樣表面TVB-N含量預測值分布直方圖，預測值等于零用紅色線條標出，預測值小于零像素占比用紅色數字標出。

圖4 豬肉新鮮度預測模型空間可視化圖像Fig.4 Spatial visualization images of pork freshness prediction model

(2)新鮮度空間分布預測精度對比

將準度相當的化學計量學模型應用到像素光譜，指標空間分布預測精度會存在明顯差異，圖5為豬肉新鮮度18種空間分布預測的精度與像素光譜質量及預測模型波段增益的關系圖。圖5a、5c、5e、5g、5i、5k、5m、5p、5r分別與圖4a～4i所示模型對應。從中可以看出，大部分模型隨濾波器寬度的增加，預測精度逐漸提升；但當濾波器寬度增加到一定程度時，預測精度開始下降。圖5a～5f、圖5g～5l、圖5m～5r分別表示基于6個特征波段、20個特征波段以及全波段建模時，原始光譜及濾波器寬度依次為6、18、30、42、54 nm濾波后光譜的空間分布預測精度與像素光譜質量及預測模型波段增益的關系。橢圓的橫軸半徑表示該模型各波段增益標準差，縱軸半徑表示暗噪標準差；箱線圖上下邊緣分別表示預測值小于零像素占比的最大值和最小值，箱子上下兩端邊位置分別表示預測值小于零像素占比的上下四分位數，箱子中間紅色線條表示預測值小于零像素占比的中值，極端異常值用紅色“+”標出。可以發現，在同一維度上，即特征波段數量相等或同為全波段建模時，橢圓的橫軸半徑和縱軸半徑越小，箱子面積就越小，箱子位置越低，模型的預測精度越高。

圖5 豬肉新鮮度18種空間分布預測的精度與像素光譜質量及預測模型波段增益的關系Fig.5 Relationships between precision of 18 spatial distribution prediction of pork freshness and spectral quality of pixels and band gain of prediction model

基于6個特征波段建模時，隨濾波器寬度的增加，橢圓縱軸半徑逐漸減小，暗噪標準差逐漸降低，由20.58 DN下降至7.91 DN，像素光譜質量逐漸提升；圖5a～5e模型橢圓橫軸半徑無明顯差異，各波段增益標準差均處于107.17～149.45的范圍內。此時箱子面積逐步減小，箱子位置逐漸降低，預測精度逐步提升。而圖5f模型相較于圖5e模型，橢圓橫軸半徑有所增大，各波段增益標準差達到205.59，增長幅度為37.56%；在暗噪標準差減小幅度僅為4.24%的情況下，箱子面積增大，箱子位置升高，預測精度略有下降。

基于20個特征波段建模時，相同的預測精度變化趨勢再次出現。隨濾波器寬度的增加，橢圓縱軸半徑逐漸減小，暗噪標準差逐漸降低，由20.51 DN下降至7.86 DN，像素光譜質量逐漸提升；圖5g～5k模型橢圓橫軸半徑雖略有起伏，但各波段增益標準差均處于45.66～90.97的范圍內。此時箱子面積逐步減小，箱子位置逐漸降低，預測精度逐步提升。而圖5l模型相較于圖5k模型，橢圓橫軸半徑有所增大，各波段增益標準差達到99.55，增長幅度為31.51%；在暗噪標準差減小幅度僅為4.97%的情況下，箱子面積增大，箱子位置升高，預測精度略有下降。

基于全波段建模時，隨濾波器寬度的增加，橢圓縱軸半徑逐漸減小，暗噪標準差逐漸降低，由20.41 DN下降至7.83 DN，像素光譜質量逐漸提升；圖5m～5p模型橢圓橫軸半徑相差無幾，各波段增益標準差最大僅為43.85。此時，箱子面積逐步減小，箱子位置逐漸降低，預測精度逐步提升。而圖5q模型相較于圖5p模型，橢圓橫軸半徑有所增大，各波段增益標準差達到114.85，增長幅度為1 036.00%，遠超過7.23%的暗噪標準差減小幅度；圖5r模型橢圓橫軸半徑進一步增大，各波段增益標準差達到287.54，增長幅度為150.36%，而暗噪標準差減小幅度僅為4.63%；與此同時，圖5q和圖5r模型的箱子面積明顯增大，箱子位置升高，預測精度顯著下降。

總之，預測精度與像素光譜質量的變化趨勢非常接近，即像素光譜質量越高，預測精度越高；而異常模型的出現則是由于各波段增益標準差明顯增大，增長幅度遠大于暗噪標準差減小幅度所造成的。

像素光譜質量與空間分布預測精度的相關度散點圖如圖6所示。暗噪標準差與預測值小于零像素占比中值的擬合方程為Y1=1.295X-6.672，大部分藍色圓點分布趨勢貼近擬合直線，其相關系數為0.72，表明兩者具有良好的一致性，即像素光譜質量與空間分布預測精度相關度較高。而左上方界外點的出現，則是受到各波段增益的影響，說明單純考慮化學計量學模型的準度指標，忽視其波段增益數值偏大，可能導致其應用到像素光譜時精度極差。暗噪標準差與像素預測值標準差中值的擬合方程為Y2=0.591X+0.456，橙色三角形分布非常貼近擬合直線，且沒有界外點出現，其相關系數為0.76，說明兩者同樣具有良好的一致性，即隨像素光譜質量的提升，暗噪標準差減小，此時像素預測值標準差隨之減小，興趣區域內各像素點預測值分布范圍逐漸減小，分布更加精確集中，并不是整體向預測值增大的方向移動。

圖6 像素光譜質量與空間分布預測精度相關度Fig.6 Pearson correlation coefficient between pixel spectral quality and spatial distribution prediction precision

綜上所述，新鮮度指標空間分布預測的精度明顯受到像素光譜質量及化學計量學模型波段增益值的共同影響，其中像素光譜質量占主導作用。因此，在應用線性化學計量學算法，將基于興趣區域均值光譜所建化學計量學模型應用到像素光譜進行新鮮度指標空間分布預測時，應充分提高像素光譜信噪比，同時注意限制化學計量學的波段增益值，以提高空間分布預測的精度。

3 結論

(1)雖然受限于缺乏像素位置微觀指標參考值無法進行直接評價，但通過對像素位置微觀預測值的統計均值進行準度評價及根據指標理論正確值的有效范圍進行精度評價的方法，能夠對基于光譜成像的化學計量學指標空間分布預測質量進行綜合評價。

(2)在應用線性化學計量學算法，將基于興趣區域均值光譜所建化學計量學模型應用到像素光譜進行新鮮度指標空間分布預測時，其準度不會下降。

(3)新鮮度指標空間分布預測的精度明顯受到像素光譜質量及化學計量學模型波段增益值的共同影響，其中像素光譜質量占主導作用(R=0.72)。

(4)在應用線性化學計量學算法，將基于興趣區域均值光譜所建化學計量學模型應用到像素光譜進行新鮮度指標空間分布預測時，高準度化學計量學模型不一定適用于像素光譜進行空間分布預測，區域預測效果與空間分布預測效果并無必然聯系；可以通過提高像素光譜信噪比和限制模型波段增益提高預測的精度。