本技術(shù)方案屬于鮮肉檢測技術(shù)領域，具體是一種鮮肉新鮮度k值的thz光譜快速無損檢測方法。

背景技術(shù)：

肉是人們主要的蛋白質(zhì)、脂肪等的攝取源。鮮肉的新鮮程度在食品安全中尤其重要。以豬肉為例，它是我國肉制品消費的主要種類，我國豬肉消費量占肉類消費總量的77％，如何控制生產(chǎn)銷售環(huán)節(jié)中的豬肉品質(zhì)、保證食品安全，受到了生產(chǎn)企業(yè)、消費者和國家質(zhì)檢部門的高度關(guān)注。與包括豬、牛、羊、魚肉等鮮肉的品質(zhì)包括營養(yǎng)成分、風味、嫩度、保水性和新鮮度等要素，其中，新鮮度是一個重要的方面。

新鮮度測定通常有兩種方法，感官評定法和理化分析法。感官評定法主觀因素大，重復性差；理化分析法精確可信，重復性高。理化分析所測定的新鮮度指標包括肉色、微生物指標、揮發(fā)性鹽基總氮指標、k值、ph值和聚胺類化合物等，其中，基于atp分解過程的k值被證實是可行的且受到越來越多的重視。

鮮肉中的一些物質(zhì)隨保存時間的增加而變質(zhì)，在宰后的肉品中，三磷酸腺苷(atp)會自動分解，過程為：三磷酸腺苷(atp)→二磷酸腺苷(adp)→磷酸腺苷(amp)→次黃嘌呤核苷酸(imp)→肌苷酸(hxr)→次黃嘌呤(hx)。

據(jù)atp分解過程，有如下測定k值的指標方程：

式中atp、adp、amp、imp、hxr、hx——各物質(zhì)在樣品中的含量

肉品越是新鮮k值越低，反之k值越高。k值通?？梢圆捎酶咝б合嗌V分析法(highperformanceliquidchromatography,hplc)獲得。雖然此種理化分析法準確且可信，但是過程耗時、對被測樣品具有破壞性(有損檢測理化參數(shù)的方法是把肉剁成肉泥，再進一步加工化驗)，不能滿足生產(chǎn)和流通環(huán)節(jié)中的快速和廣泛性測量要求，用快速無損檢測方法測定豬肉的k值正成為科研工作的一個重要研究熱點。

目前，鮮肉新鮮度指標的無損檢測方法有近紅外光譜分析、可見-近紅外高光譜圖像分析、電子鼻技術(shù)等。光譜分析法基于與鮮肉新鮮度相關(guān)的化學物質(zhì)在可見和近紅外波段有特征表達的原理，建立起無損檢測豬肉新鮮度的數(shù)學模型，實現(xiàn)了新鮮度的快速和比較準確的測定。但是近紅外光譜儀實驗過程中，需要對樣品倉抽真空處理，無法實現(xiàn)在線快速檢測，光纖式光譜儀雖然不需要抽真空處理，但是需要用光纖采集樣本的光譜信息，耦合光纖細小，檢測區(qū)域有限，無法實現(xiàn)對樣本表面的一次性全面采集?？梢?近紅外高光譜圖像分析技術(shù)雖然能夠?qū)崿F(xiàn)全面采集，但是每次試驗開始前以及長時間試驗過程中，都需要用標準反射率標定板對光譜采集系統(tǒng)進行校準，以防止相機感光器件的溫度漂移以及光源的亮度變化對檢測系統(tǒng)的影響，不利于在線連續(xù)檢測。電子鼻技術(shù)使用氣體傳感器采集肉的腐敗氣味，基于腐敗程度越深、氣味釋放程度越強的原理測定肉的新鮮度，但是此方法需要將樣品放置在密閉的檢測倉中，并需要穩(wěn)定放置一定時間，以使得檢測倉中氣味分布均勻，氣味濃度達到傳感器的檢測下限，檢測的等待時間長。

技術(shù)實現(xiàn)要素：

太赫茲波(terahertz,thz)位于毫米波和紅外線之間，屬于遠紅外波段，頻率在0.1thz～10thz范圍內(nèi)(本案選擇頻率在0.1thz～2thz范圍的太赫茲波)。從能量上看，thz波處于電子和光子之間。thz波的多元化特性使得很多化學分子在thz波段下表現(xiàn)出比在其它波段下所不具備的分子運動特性。核苷酸類(如adp、atp)及其相關(guān)物質(zhì)(如imp)屬于生物小分子，在太赫茲波段的吸收主要是由于其分子自身的轉(zhuǎn)動、振動或分子集團的整體振動，有太赫茲波譜特征結(jié)構(gòu)，在thz波段存在多個吸收峰。

因此，本發(fā)明以k值為研究參數(shù)，通過探索鮮肉新鮮度的thz光譜特性，建立thz光譜數(shù)據(jù)與鮮肉新鮮度k值之間的關(guān)系模型，實現(xiàn)對鮮肉k值的快速無損檢測，提出一種鮮肉新鮮度快速無損檢測新方法，具體技術(shù)方案如下：

一種鮮肉新鮮度k值的thz光譜快速檢測方法，其特征是包括步驟：

一)肉品采集和處理：將鮮肉均勻切割成肉片后冷藏；取樣切割時，避開肉的脂肪和結(jié)締組織；

二)測定肉品的thz光譜數(shù)據(jù)：采用thz衰減全反射attenuatedtotalreflectance,atr檢測模式快速檢測樣品；光譜采集時，保持環(huán)境的溫、濕度穩(wěn)定；

三)采用thz光譜數(shù)據(jù)與鮮肉k值之間的關(guān)系預測模型快速檢測鮮肉的新鮮度k值：所述關(guān)系模型是主成分回歸pcr預測模型、或者是反向傳播神經(jīng)網(wǎng)絡回歸bpann預測模型；

將步驟二)得到的thz光譜數(shù)據(jù)代入關(guān)系預測模型，則快速計算出未知樣本的k值；

所述步驟二)中：

采用thz衰減全反射atr檢測模式對鮮肉樣本的thz光譜進行取樣；對鮮肉樣本的thz光譜進行取樣得到樣本光譜數(shù)據(jù)進行處理；

所有建模樣本被隨機分為校正集和預測集；(校正集中樣本數(shù)和預測集中樣本數(shù)的比約為2:1，)校正集的光譜數(shù)據(jù)用于建立預測模型；預測集的光譜數(shù)據(jù)用來檢驗預測模型的準確性。xn×m是特指校正集中的光譜數(shù)據(jù)。

樣本光譜數(shù)據(jù)矩陣為式中，n為校正集樣本數(shù)，m為光譜波點數(shù)，x是樣本的atr反射率；

先對光譜信息進行壓縮：把矩陣xn×m分解為兩個矩陣u和p，表示為un×a＝xn×mpm×a；式中，矩陣u表示數(shù)據(jù)xij(i＝1,2,…,n；j＝1,2,…,m)在新的坐標體系中的矢量位置；矩陣p的列向量表示新的坐標系與原坐標系的線性變換系數(shù)；u的維數(shù)a小于x的維數(shù)m；取矩陣u的各行數(shù)據(jù)ii(i＝1,2,…,a)作為步驟三)的預測模型的輸入；

所述步驟三)中：

a、所述pcr預測模型：

對矩陣u進行多元線性回歸，得到多元線性方程：y＝ub+e，即為pcr預測模型，用來預測k值；式中，y即為樣本集的實測k值，e是服從正態(tài)分布的隨機誤差變量矩陣；

所述參數(shù)b的獲得方法為：在pcr預測模型的建模過程中，把多個樣本隨機分為多個校正集和多個預測集；校正集的光譜數(shù)據(jù)用于建立預測模型；預測集的光譜數(shù)據(jù)用來檢驗預測模型的準確性；將校正集樣品的實測k值組成矩陣y，得到方程的解：b＝(u′u)^-1u′y，式中，u′為u的轉(zhuǎn)置矩陣。得到參數(shù)矩陣b后就可以用預測集樣本的thz光譜數(shù)據(jù)代入上述多元線性方程求出預測集中樣本的預測k值。

b、所述bpann預測模型的建模步驟包括：

1)構(gòu)建多層反向傳播神經(jīng)網(wǎng)絡：

把多層反向傳播神經(jīng)網(wǎng)絡bpann視為非線性函數(shù)，bpann輸入和預測值分別為該非線性函數(shù)的自變量和因變量，bpann表達了從自變量到因變量的函數(shù)映射關(guān)系；

bpann的拓撲結(jié)構(gòu)由三層網(wǎng)絡組成，分別是輸入層、隱含層和輸出層，每層包括多個神經(jīng)元，相鄰層的神經(jīng)元之間是帶權(quán)值的單向連接；

ii是輸入層的輸入值，i＝1,2,…,a；hj是隱含層的輸出值，j＝1,2,…,b；輸出層是一個節(jié)點的值為bpann的預測值，即鮮肉的預測k值；hwij為輸入層的神經(jīng)元與隱含層的神經(jīng)元之間的連接權(quán)值，owj為隱含層的神經(jīng)元與輸出層的神經(jīng)元之間的連接權(quán)值；

隱含層的節(jié)點數(shù)式中：a為輸入層節(jié)點數(shù)，c為1～10范圍內(nèi)的常數(shù),1為輸出層節(jié)點數(shù)；

2)對步驟1)的模型進行訓練，步驟包括：

2.1)網(wǎng)絡初始化：根據(jù)輸入維數(shù)和輸出維數(shù)確定a，b；初始化hwij、owj、隱含層閾值thaj、輸出層閾值thb，權(quán)值和閾值是范圍在-1～1的隨機數(shù)；給定學習速率η，η＝0.01；

2.2)隱含層輸出計算：隱含層的輸出式中，函數(shù)f(x)為隱含層激勵函數(shù)，該函數(shù)的表達式為：

2.3)輸出層輸出計算：bpann的預測輸出

2.4)誤差計算：bpann的預測誤差式中，y是采用hplc方法實測的k值；

2.5)權(quán)值更新：hwij＝hwij+ηhj(1-hj)iiowje，i＝1,2,…,a；j＝1,2,…,b；owj＝owj+ηhje，j＝1,2,…,b；

2.6)閾值更新：thaj＝thaj+ηhj(1-hj)owje，j＝1,2,…,b；thb＝thb+e；

2.7)根據(jù)迭代次數(shù)和最終迭代誤差判斷算法迭代是否結(jié)束，若沒有結(jié)束，返回步驟2.2)。

所述bpann模型的建模步驟還包括步驟3)對bpann模型進行優(yōu)化，得到bpann改進模型，方法是把對bpann預測模型作為弱預測器進行優(yōu)化，得到強預測器，步驟包括：

3.1)初始化：n為校正集的樣本數(shù)，初始化校正集中訓練數(shù)據(jù)的分布權(quán)重為dt(j)＝1/n(j＝1,2,…,n)，確定弱預測器預測誤差的閾值φ，根據(jù)樣本的輸入輸出維數(shù)確定bpann弱預測器結(jié)構(gòu)，初始化網(wǎng)絡參數(shù)；

3.2)訓練弱預測器：訓練第t個弱預測器時，用訓練數(shù)據(jù)訓練弱預測器，得到弱預測器輸出的預測k值

然后計算預測誤差t是弱預測器的個數(shù)；

3.3)計算弱預測器的權(quán)重

3.4)調(diào)整訓練數(shù)據(jù)分布權(quán)重并歸一化：式中，

3.5)將3.2～3.4)步驟循環(huán)t次后，得到t個弱預測器，預測模型為ft。按更新的權(quán)重疊加，得到強預測器的預測模型

所述φ值是采用留一交叉驗證法選取模型的性能最好時的φ值，即在一定范圍內(nèi)，以固定的步長，考察φ取不同的值時，建模步驟在留一交叉驗證法下的rmsecv的值，取rmsecv的值為最小時的φ作為最優(yōu)閾值；

留一交叉驗證法的運算過程是：在校正集的n個樣本中，每個樣本都單獨作為預測集，其余n-1個樣本作為校正集；以校正集為依據(jù)，通過建模步驟得到這個預測集樣本的預測值；經(jīng)過n次循環(huán)，得到所有樣本的預測值；

獲得留一交叉驗證法的性能指標rmsecv，

所述步驟二)中，光譜頻段為0.1thz～2thz。

所述步驟二)中，先對樣本光譜數(shù)據(jù)進行預處理，預處理得到的數(shù)據(jù)作為矩陣x的數(shù)據(jù)；樣本光譜數(shù)據(jù)進行預處理的方法為：

先采用一階導數(shù)fd預處理光譜數(shù)據(jù)：式中，x是樣本在波點數(shù)為i處的atr反射率；當i＝1時，其一階微分值為……；當i＝m時，其一階微分值為

再采用sg多項式平滑光譜：選擇奇數(shù)個光譜數(shù)據(jù)序列中相鄰的數(shù)據(jù)點，由這些數(shù)據(jù)點構(gòu)成一個窗口，對窗口的中心點進行平滑；平滑過程中，對窗口中的數(shù)據(jù)點擬合一個多項式，并根據(jù)所得該多項式對平滑的點進行計算；得到平滑后的數(shù)據(jù)點表達式為：式中，窗口的寬度為2d+1，s為常數(shù)，λj為sg平滑系數(shù)序列；

該多項式的擬合運用最小二乘方法。

本發(fā)明的樣品前處理的方法，只需要切割樣本以形成肉片放入儀器中測量，比現(xiàn)有的方法簡單且快速了許多。

附圖說明

圖1是atr附件光路示意圖；

圖2-1是hplc測得的k值變化趨勢圖；

圖2-2是一個肉樣連續(xù)測量6次所得的6條thz原始光譜圖；

圖2-3是經(jīng)svsrs處理后的光譜圖；

圖3-1是80個樣本的thz光譜原始譜線圖；

圖3-2是經(jīng)過15點一階求導和sg平滑后的thz光譜譜線圖；

圖4-1是不同的差分寬度下pcr模型的預測性能變化圖；

圖4-2是pcr模型的k值預測散點圖；

圖5-1是bpann模型的結(jié)構(gòu)圖；

圖5-2是bpann模型的算法流程圖；

圖5-3是bpann模型的rmsecv與差分寬度和隱含層節(jié)點數(shù)的關(guān)系圖(rmsecv的最小值用紅圈標注)；

圖5-4是bpann模型的k值預測散點圖；

圖6-1是bpann改進模型模型的算法流程圖；

圖6-2是bpann改進模型模型中rmsecv與閾值φ關(guān)系圖(閾值φ＝0.05～0.25)；

圖6-3是bpann改進模型算法中rmsecv與閾值φ關(guān)系圖(閾值φ＝0.13～0.14)；

圖6-4是bpann改進模型模型的k值預測散點圖。

具體實施方式

下面對本技術(shù)方案的研發(fā)、實驗過程進行說明：

1、thz光譜本身特性：

太赫茲技術(shù)的一個顯著特點是安全性。相比于x射線具有千電子伏的光子能量，太赫茲輻射的光子能量低于各種化學鍵的鍵能，因此它不會引起有害的電離反應。另外，由于水對太赫茲波有強烈的吸收，太赫茲輻射不能穿透人體的皮膚。因此即使強烈的太赫茲輻射，對人體的影響也只停留在皮膚表層，而不像微波可以穿透到人體的內(nèi)部。thz安全性有利于設計出安全便捷的無損檢測設備，并保護測試人員免受輻射危害。

太赫茲技術(shù)的另一個顯著特點是光譜可分辨性。與近紅外和中紅外輻射相比，處于遠紅外波段的thz輻射的光子能量相對較低，但這一波段仍然包含了豐富的光譜信息。有機分子，由于其轉(zhuǎn)動和振動(包括集體振動)的躍遷，在這一頻段表現(xiàn)出強烈的吸收和色散特性。

2、鮮肉新鮮度k值相關(guān)的分子物質(zhì)在thz波段的體現(xiàn)：

鮮肉新鮮度——k值測量的是核苷酸類(如adp、atp)及其相關(guān)物質(zhì)(如imp)的含量比，它們屬于生物小分子，在太赫茲波段的吸收主要是由于其分子自身的轉(zhuǎn)動、振動或分子集團的整體振動，其太赫茲波譜特征結(jié)構(gòu)較為明晰。嘌呤多晶體在0.2～2.5thz波段存在多個吸收峰。

3、thz檢測方式的選擇：

由于生物分子中包含較多的原子，并且具有高密度的集合振動模式；分子中的各個原子之間的作用導致生物分子振動具有較大的非簡諧性，因此生物分子的太赫茲光譜往往被非均勻展寬并重疊在一起。

可以采用兩種方式解決這一問題。一種，使用低溫冷卻方式(如液氮冷卻)壓窄吸收峰的展寬，分離出各吸收峰，觀察其吸收特征；另一種，在常溫測量方式下獲得thz譜線，對譜線使用多種預處理方式，比如標準正態(tài)變換(standardnormalvariation,snv)、多元散射校正(multi-scattercalibration,msc)、導數(shù)(包括一階導數(shù)和二階導數(shù))、savitzky-golay(sg)多項式平滑、小波變換等。盡量減少采集到的原始光譜中的高頻隨機噪聲、基線漂移、光散射等噪聲信息，使用線性或非線性數(shù)學模型分析譜線和被測物質(zhì)的復雜關(guān)系，建立可靠和穩(wěn)定的thz光譜分析模型。

考慮到本方案的目的是對被測鮮肉樣品實現(xiàn)無損快速檢測，應在常溫方式下測量，所以選擇后一種方式。

4、thz檢測樣本的選擇

鮮肉(例如豬肉)樣品中的表皮，肌肉和脂肪在thz波段下具有不同的折射率和透射率。將一塊肉的切片中不同的組織成分——肥肉和瘦肉，放置于thz光譜分析系統(tǒng)中，在atr反射模式的檢測范圍內(nèi)(0.1至4thz)，瘦肉的吸收系數(shù)大于肥肉的吸收系數(shù)，而且隨著thz波的頻率增加，兩者之間的差異越來越大?？梢?，肉的不同組織成分，thz波的光譜差異較大。

一般來說，消費者在購買肉制品時更喜歡購買較瘦的肉，所以瘦肉的品質(zhì)會決定著整個肉的品質(zhì)。本技術(shù)方案選擇樣本中的瘦肉部分作為thz波譜的檢測對象。

下面結(jié)合鮮豬肉的實施例對本方法進行進一步說明：

1、材料和方法

1.1試驗材料

試驗材料為冷鮮豬肉的通脊肉，每天上午從當?shù)爻匈徺I，用冷藏箱運回實驗室后，將豬肉均勻切割成2.5cm×2.5cm×0.5cm的肉片。取樣時避開豬肉的脂肪和結(jié)締組織，以防止這些成分對thz檢測結(jié)果的干擾。將豬肉樣品用保鮮袋包好并編號后置于4℃冰柜中貯藏待測。連續(xù)采集8天的肉樣，以日期先后為編號，依次為7d～0d。在0d的下午，完成8個肉樣的thz光譜采集以及k值測定。重復試驗10次，得到80個肉樣的光譜數(shù)據(jù)和理化值。

1.2thz光譜采集

thz檢測設備的型號是tas7500sp(advantest公司，日本)，在室溫環(huán)境下工作，頻率分辨率為7.6ghz，檢測頻率范圍是0.1thz～4thz，共498個采樣頻率點(波點)，樣品的譜線經(jīng)過2048次自動掃描并平均后得到。

thz與水有強烈的相互作用，對于富含水的肉樣透射深度只有幾百微米，所以透射模式不適合肉制品新鮮度無損檢測，反射模式也會因為反射波被水吸收而無法檢測。用thz衰減全反射(attenuatedtotalreflectance,atr)檢測模式，可以克服樣品中富含的游離水以及結(jié)合水對thz波的強烈吸收，使得樣品表面微米級厚度化學物質(zhì)的thz特性能夠反映在thz全反射波的光譜里。atr檢測附件的光路如圖1所示。

每個樣品從冰箱取出后，去除包裝，將樣本平整放入atr檢測窗表面，如圖1所示。樣品的上下兩個表面分別采集3次thz光譜數(shù)據(jù)，每個樣本能夠獲得6份thz光譜數(shù)據(jù)，將6份光譜數(shù)據(jù)算術(shù)平均，作為該樣本的最終thz光譜數(shù)據(jù)。因為thz光譜儀對溫度和濕度比較敏感，所以光譜采集時，保持實驗室內(nèi)溫度、濕度基本一致。

1.3k值測定

采集完樣本的thz光譜數(shù)據(jù)后立即對同一樣本進行k值測定。采用hplc方法檢測肉制品中的k值。利用流動相中atp分解的關(guān)聯(lián)產(chǎn)物(atp、adp、amp、imp、hxr和hx)在固定相中的流速不同，將這六種化學組分分離，分別測得這些組分的含量，根據(jù)公式1計算出被測樣品的k值。

用于k值檢測的樣品前處理過程如下：將被測樣品剁碎成肉泥，從中取(2.00±0.05)g放入50ml離心管內(nèi),加入冷卻后的10％的高氯酸溶液20ml,渦旋振蕩1min,以8000r/min速度離心10min,取出上清液。再用5％的高氯酸溶液20ml重提沉淀物中的待測物,以8000r/min速度離心10min,合并上清液。用10mol/l的naoh溶液調(diào)提取液ph值近6.0,然后再用1.0mol/l的naoh溶液繼續(xù)調(diào)節(jié)ph至6.0～6.4,再用超純水定容至50ml。用0.45μm的微孔濾膜過濾,濾液于4℃下保存,待測。

hplc條件如下：finnigansurveyor液相色譜儀(賽默飛公司，美國)，aq-c18色譜柱(賽默飛公司，美國)，流動相為0.05mol/l的k3po4緩沖液(ph＝6.5)，緩沖液用超純水配制，樣品的進樣量1μl,流速200μl/min,檢測波長254nm。atp分解的關(guān)聯(lián)產(chǎn)物由外標法定量，測定范圍在0～0.5mmol/l。

試驗中使用的試劑詳情如下：atp關(guān)聯(lián)物(atp、adp、amp、imp、hxr和hx)共6種標準品(純度≥99％,sigma-aldrich公司，美國)；磷酸鉀、氫氧化鈉、高氯酸(分析純,國藥集團化學試劑有限公司，中國)；試驗用水為milliporeacademic制備的超純水。

2、結(jié)果與分析

2.1k值分析

采用hplc法獲得了80個樣本的k值，并做了統(tǒng)計分析，結(jié)果如圖2-1所示。不同試驗批次的肉樣，存儲時長一樣，k值也會因為個體差異而不同。通過樣本的平均值可以看出，隨著儲存天數(shù)的增加，新鮮度不斷降低，樣本的k值在逐步增加，但增加量不是固定的。其中，5d～6d、6d～7d的增加量相對較大，說明在上述存儲周期內(nèi)，肉樣的品質(zhì)變化較大；感官評價也證實肉樣在6d～7d存儲周期內(nèi)，觸摸較黏，有明顯的酸味，因此可推測肉樣在此周期內(nèi)新鮮度顯著下降。本試驗所有樣本的k值覆蓋了豬肉的不同新鮮度，覆蓋范圍廣，使得thz光譜預測k值的數(shù)學模型有較好的魯棒性。

80個樣本被隨機分為54個校正集和26個預測集，個數(shù)比大致為2:1。其中校正集用于建立thz光譜預測k值的數(shù)學模型；預測集用來檢驗所建立的模型預測未知樣本k值的準確性。通過表1可以看出，校正集、預測集和樣本總集合的k值范圍基本相同，平均值和標準差也沒有明顯區(qū)別，因此校正集和預測集的樣本分割是合適的。

表1校正集與預測集的k值統(tǒng)計信息

2.2光譜預處理

光譜質(zhì)量評價方法：肉樣的同一表面在相同的測試條件下連續(xù)測量6次，獲得6條thz光譜，理論上這6條光譜應該完全重合，但由于儀器的噪聲及測試誤差的影響，連續(xù)重復測得的6條光譜不可能完全重合，為了評價光譜的質(zhì)量，可以計算同一樣品表面連續(xù)多次重復測試的ns條光譜的標準方差光譜(standardvariancespectrumofrepeatspectral,svsrs)，svsrs越小，說明光譜質(zhì)量越好。

式中xij——樣本第i次測量在波點j處的atr反射率,——波點j處ns次測量atr反射率的均值。

圖2-2是某個肉樣表面在0.1thz～4thz的6次thz光譜，圖2-3是該光譜對應的svsrs，從圖中可以看出在2thz～4thz的svsrs值明顯增高，重復性差。故選擇0.1thz～2thz作為建模的光譜頻段。

圖3-1是80個樣本在0.1thz～2thz區(qū)域的thz光譜譜線，可以看出，不同豬肉樣本的原始光譜強度有很大的差異。光譜的差異不僅包含了樣本成分的差異，還包括測量誤差，基線漂移和背景噪聲。為了消除干擾信息，除了盡可能保持試驗環(huán)境因素一致外，還須對光譜數(shù)據(jù)進行預處理，以減弱或去除各種干擾因素。

本發(fā)明采用一階導數(shù)(firstorderderivative，fd)預處理光譜數(shù)據(jù)，一階導數(shù)能夠減少基線偏移、漂移和背景干擾造成的數(shù)據(jù)偏差，使得與新鮮度密切相關(guān)的光譜特性變得更為顯著。

本發(fā)明使用的一階導數(shù)fd預處理光譜數(shù)據(jù)：式中，x是樣本在波點數(shù)為i處的atr反射率；當i＝1時，其一階微分值為……；當i＝m時，其一階微分值為

由于導數(shù)計算會增加噪聲，故導數(shù)預處理之后采用savitzky–golay(sg)多項式平滑光譜。在求導過程中，差分寬度選擇十分重要：如果差分寬度太小，噪聲會很大，影響所建模型的質(zhì)量；如果差分寬度太大，平滑過度，會失去大量的細節(jié)信息。sg多項式平滑選擇奇數(shù)個數(shù)據(jù)序列中相鄰的數(shù)據(jù)點，由這些數(shù)據(jù)點構(gòu)成一個窗口，對窗口的中心點進行平滑。平滑中，對窗口中的數(shù)據(jù)點擬合一多項式，并根據(jù)所得多項式對平滑的點進行計算。多項式的擬合運用最小二乘方法。根據(jù)計算可以得到平滑后的數(shù)據(jù)點表達式為：

式中，窗口的寬度為2d+1；s為常數(shù)，其數(shù)值與窗口寬度有關(guān)；λj為sg平滑系數(shù)序列。

圖3-2為經(jīng)過15點一階求導和sg平滑后得到的一階導數(shù)光譜圖。

2.3豬肉k值預測模型

由豬肉腐敗過程中會使k值相關(guān)的atp關(guān)聯(lián)產(chǎn)物的含量發(fā)生變化，而這些生物分子對thz波具有靈敏的光譜響應，thz光譜能夠反映出這些生物分子含量的變化。因此thz光譜數(shù)據(jù)和豬肉的新鮮度之間存在一種間接的關(guān)聯(lián)性。本方法使用主成分回歸(principalcomponentsregression，pcr)、非線性算法-誤差反向傳播神經(jīng)網(wǎng)絡回歸(backpropagationartificialneuralnetwork,bpann)以及bpann改進模型分別驗證這種關(guān)聯(lián)性。

檢測模型評價方法：

評價模型質(zhì)量的指標是采用校正集相關(guān)系數(shù)(rc)、校正集均方根誤差(rootmeansquarederrorofcalibrationset，rmsec)、預測集相關(guān)系數(shù)(rp)和預測集均方根誤差(rootmeansquarederrorofpredictionset，rmsep)。這些參數(shù)定義如下：

式中——樣本集(包括校正集和預測集)中第i個樣品的預測值；

yi——樣本集(包括校正集和預測集)中第i個樣品的實測值；

nc——校正集的樣品數(shù)；

np——預測集的樣品數(shù)；

ym——校正集或者預測集樣品的平均值。

相關(guān)系數(shù)rc和rp越大，和/或rmsec和rmsep越小，則模型的預測能力越好。

本例使用matlabr2009b(mathworks公司，美國)軟件對光譜數(shù)據(jù)進行計算與建模。

2.3.1主成分回歸pcr預測模型

可以用矩陣形式表示80個樣本的thz光譜數(shù)據(jù)：

式中，n為樣品數(shù)，m為光譜波點數(shù)，x是樣本的atr反射率；xn×m光譜數(shù)據(jù)是指校正集中的數(shù)據(jù)。

本試驗中的每一條光譜都含有250個波點，光譜的波點數(shù)遠大于樣本數(shù)，如果直接用于回歸分析，會出現(xiàn)過擬合，降低模型的預測精度和穩(wěn)定性。同時光譜波點的數(shù)據(jù)間存在較高的共線性和相關(guān)性，也會使得回歸模型的結(jié)果產(chǎn)生失真。可以通過主成分分析法(principalcomponentanalysis，pca)對光譜信息進行壓縮，通過少數(shù)幾個主成分的得分來近似反映原光譜數(shù)據(jù)，消除光譜數(shù)據(jù)點的信息冗余和相關(guān)性。

pca算法是一種常用并有效的信息壓縮方法，它從多個樣本構(gòu)成的變量協(xié)方差矩陣出發(fā)，采用特征分解的方法(即坐標系的線性變換)獲取方差最大的主成分來代替原有變量，得到了新變量——主成分得分。主成分之間相互獨立、互不相關(guān)，可消除原始數(shù)據(jù)中存在的相關(guān)性和信息冗余。因此，pca在多變量統(tǒng)計分析和光譜信息壓縮、提取等領域得到了廣泛的應用。

運用主成分分析，原變量矩陣x可以分解為兩個矩陣u和p，表示為un×a＝xn×mpm×a；

式中，式中，矩陣u是得分矩陣(scorematrix)，表示數(shù)據(jù)xij(i＝1,2,…,n；j＝1,2,…,m)在新的坐標體系中(主成分坐標體系)的矢量位置(即坐標值)；矩陣p是載荷矩陣(loadingmatrix)，其列向量表示新的坐標系與原坐標系的線性變換系數(shù)；u的維數(shù)a小于x的維數(shù)m；因此，高維的光譜數(shù)據(jù)空間被降為低維的線性無關(guān)的主成分特征空間。取矩陣u的各行數(shù)據(jù)ii(i＝1,2,…,a)作為預測模型的輸入；

pcr預測模型：對矩陣u進行多元線性回歸，得到多元線性方程：y＝ub+e，即為pcr預測模型，用來預測k值；式中，y即為樣本集的實測k值，e是服從正態(tài)分布的隨機誤差變量矩陣；

所述參數(shù)b的獲得方法為：在pcr預測模型的建模過程中，把多個樣本隨機分為多個校正集和多個預測集；校正集的光譜數(shù)據(jù)用于建立預測模型；預測集的光譜數(shù)據(jù)用來檢驗預測模型的準確性；將校正集樣品的實測k值組成矩陣y，得到方程的解：b＝(u′u)^-1u′y，式中，u′為u的轉(zhuǎn)置矩陣。得到參數(shù)矩陣b后就可以用預測集樣本的thz光譜數(shù)據(jù)代入上述多元線性方程求出預測集中樣本的預測k值。

綜上所述，pcr可以分為兩步：1測定主成分數(shù)，即用主成分分析法將光譜矩陣x降維；2對于降維的x矩陣再進行線性回歸分析。

本發(fā)明用pca法將預處理后的光譜矩陣進行主成分分解，取累積貢獻率達到了95％的主成分集合作為預測模型的輸入。光譜預處理的差分寬度根據(jù)預測模型達到最高相關(guān)系數(shù)rp和最低rmsep時選擇。模型計算表明當一階差分寬度為15時，累積貢獻率95％的主成分集合數(shù)為26，pcr的預測性能最佳(rp＝0.63,rmsep＝16.78％)，如圖4-1所示。圖4-2是pcr模型在此參數(shù)條件下的k值預測散點圖。

2.3.2誤差反向傳播神經(jīng)網(wǎng)絡回歸bpann預測模型

bpann是一個能揭示輸入與輸出信號間復雜聯(lián)系的有力分析工具。bpann是一種多層前饋神經(jīng)網(wǎng)絡，這個網(wǎng)絡可以被理解為一個非線性函數(shù)，網(wǎng)絡輸入和預測值分別為該函數(shù)的自變量和因變量。bpann表達了從自變量到因變量的函數(shù)映射關(guān)系。本發(fā)明所用的bpann的拓撲結(jié)構(gòu)如圖5-1所示。

bpann的拓撲結(jié)構(gòu)由三層網(wǎng)絡組成，分別是輸入層、隱含層和輸出層，每層包括多個神經(jīng)元，相鄰層的神經(jīng)元之間是帶權(quán)值的單向連接；

ii是輸入層的輸入值，i＝1,2,…,a；hj是隱含層的輸出值，j＝1,2,…,b；輸出層是一個節(jié)點的值為bpann的預測值，即鮮肉的預測k值；hwij為輸入層的神經(jīng)元與隱含層的神經(jīng)元之間的連接權(quán)值，owj為隱含層的神經(jīng)元與輸出層的神經(jīng)元之間的連接權(quán)值；

輸出層是一個節(jié)點，yl(l＝1)是bpann的預測值即預測k值；wij和wjk為bpann的神經(jīng)元連接權(quán)值。

為了降低神經(jīng)網(wǎng)絡的復雜度，采用上述pca方法將預處理后的光譜矩陣進行主成分分解，取累積貢獻率達到95％的主成分集合作為神經(jīng)網(wǎng)絡的輸入，單隱含層的節(jié)點數(shù)由經(jīng)驗公式得出：

式中：隱含層的節(jié)點數(shù)式中：a為輸入層節(jié)點數(shù)，c為1～10范圍內(nèi)的常數(shù)，1為輸出層節(jié)點數(shù)。據(jù)此隱含層的節(jié)點數(shù)范圍選擇為4～17。

bp神經(jīng)網(wǎng)絡預測前首先要訓練網(wǎng)絡，通過訓練使網(wǎng)絡具有聯(lián)想記憶和預測能力。bpann的訓練過程包括以下幾個步驟：

1)網(wǎng)絡初始化：根據(jù)輸入維數(shù)和輸出維數(shù)確定a，b；初始化hwij、owj、隱含層閾值thaj、輸出層閾值thb，權(quán)值和閾值是范圍在-1～1的隨機數(shù)；給定學習速率η，η＝0.01；

2)隱含層輸出計算：隱含層的輸出式中，函數(shù)f(x)為隱含層激勵函數(shù)，該函數(shù)的表達式為：

3)輸出層輸出計算：bpann的預測輸出

4)誤差計算：bpann的預測誤差式中，y是采用hplc方法實測的k值；

5)權(quán)值更新：hwij＝hwij+ηhj(1-hj)iiowje，i＝1,2,…,a；j＝1,2,…,b；

owj＝owj+ηhje，j＝1,2,…,b；

6)閾值更新：thaj＝thaj+ηhj(1-hj)owje，j＝1,2,…,b；thb＝thb+e；

7)根據(jù)迭代次數(shù)和最終迭代誤差判斷算法迭代是否結(jié)束，若沒有結(jié)束，返回步驟2)。

bpann的主要特點是信號前向傳播，誤差反向傳播。在前向傳遞中，輸入信號從輸入層經(jīng)隱含層逐層處理，直至輸出層。每一層的神經(jīng)元狀態(tài)只影響下一層神經(jīng)元狀態(tài)。如果輸出層得不到期望輸出，則轉(zhuǎn)入反向傳播，根據(jù)預測誤差調(diào)整網(wǎng)絡權(quán)值和閾值，從而使bp神經(jīng)網(wǎng)絡預測輸出不斷逼近期望輸出。當訓練結(jié)束后，就可以用訓練好后的網(wǎng)絡預測k值。

可見，基于bpann的k值預測建模包括bpann的構(gòu)建，訓練和預測三步。該算法流程用圖5-2表示。

訓練神經(jīng)網(wǎng)絡前，輸入的光譜矩陣被歸一化至-1～1的范圍內(nèi)，實測k值被歸一化至0～1的范圍內(nèi)。因為隱含層的節(jié)點數(shù)會影響神經(jīng)網(wǎng)絡預測的性能，所以隱含層的節(jié)點數(shù)采用留一交叉驗證法選取，即選取rmsecv值最小時的節(jié)點數(shù)。根據(jù)這一原則，考察預測模型在預處理差分寬度3～47范圍內(nèi)，隱含層的節(jié)點數(shù)在4～17的范圍內(nèi)，rmsecv的數(shù)值變化情況，如圖5-3所示。可以看出，預處理差分寬度為13且隱含層節(jié)點數(shù)是9時，rmsecv數(shù)值最小為18.10％，bpann模型的預測性能最佳。圖5-4是bpann模型在此參數(shù)條件下的k值預測散點圖。2.3.3bpann改進模型

本改進算法的核心思想是針對同一個校正集訓練不同的預測器(弱預測器)，然后把這些弱預測器集合起來，構(gòu)成一個更強的最終預測器(強預測器)。該算法的特點是弱預測器的性能差，但強預測器的性能卻很優(yōu)。其原因在于它能夠合理的劃分訓練集合，以及合理的合并弱預測器以形成強預測器。其合理之處具體體現(xiàn)在以下兩方面：使用加權(quán)后選取的訓練數(shù)據(jù)而非隨機選取的訓練樣本，這樣將訓練的焦點集中在比較難預測的訓練數(shù)據(jù)樣本上；使用加權(quán)的投票機制將弱預測器聯(lián)合起來，讓預測效果好的弱預測器具有較大的權(quán)重，而預測效果差的預測器具有較小的權(quán)重。本改進模型把bpann作為弱預測器，反復訓練bpann預測樣本輸出，通過改進算法將所得的多個bpann弱預測器組成強預測器。

具體改進算法如下：

1.初始化：n為校正集的樣本數(shù)，初始化校正集中訓練數(shù)據(jù)的分布權(quán)重為dt(j)＝1/n(j＝1,2,…,n)，確定弱預測器預測誤差的閾值φ，根據(jù)樣本的輸入輸出維數(shù)確定bpann弱預測器結(jié)構(gòu)，初始化網(wǎng)絡參數(shù)；

2.訓練弱預測器：訓練弱預測器：訓練第t個弱預測器時，用訓練數(shù)據(jù)訓練弱預測器，得到弱預測器輸出的預測k值

然后計算預測誤差t是弱預測器的個數(shù)；

3.計算弱預測器的權(quán)重wt：

4.測試數(shù)據(jù)權(quán)重調(diào)整，

式中，bt是歸一化因子

5.將2,3,4步驟循環(huán)t次后，得到t個弱預測器，預測模型為ft。按更新的權(quán)重疊加，得到強預測器的預測模型

改進算法流程如圖6-1所示，可見它是一種整合算法，通過上述加權(quán)策略整合集合內(nèi)各弱預測器的預測能力而獲得總體預測性能的提高。本發(fā)明采用bpann改進算法，以前面建立的bpann模型參數(shù)為基礎，使用改進算法提升其預測性能，構(gòu)建k值的bpann改進預測模型。

預測誤差的閾值φ對bpann改進模型的性能有顯著的影響，采用留一交叉驗證法選取rmsecv值最小時的φ：首先，以較寬的步長0.01在較大的范圍選擇φ為0.05～0.25，rmsecv值的情況如圖6-2所示?？梢钥闯?，在0.13～0.14的范圍內(nèi)，模型的性能比較好。接著，在此范圍內(nèi)，以較小的0.001步長考察0.13～0.14，可以看出圖6-3中，φ取0.136時，rmsecv值最小，bpann改進模型能夠獲得最優(yōu)的預測性能。將此參數(shù)代入bpann改進算法中迭代，獲得了10個弱bpann預測器，這些預測器的性能如表2所示?？梢钥闯觯喝躅A測器的性能有差異，其中弱預測器2和8的rmsec較小，預測性能較好，其預測誤差和εt較小，迭代所生成的權(quán)重wt較大，代表著在整合后的強預測器中這兩個弱預測器有較大比重；相反，弱預測器6的rmsec較大，預測性能最差，其預測誤差和εt最大，迭代所生成的權(quán)重wt最小，代表著在整合強預測器的過程中這個弱預測器貢獻很?。簧鲜?0個弱預測器的預測性能都不理想，但是當用權(quán)重wt加權(quán)整合弱預測器后，所獲得的強預測器卻有較好的性能，其預測性能優(yōu)于各弱預測器。通過這個強預測器獲得bpann改進模型的k值預測散點圖，如圖6-4所示。

表2bpann改進模型的預測器性能

2.4thz光譜預測模型的比較

上述三種預測模型采用預測集數(shù)據(jù)驗證模型的k值預測準確度，方法為：對于預測集數(shù)據(jù)，取經(jīng)過一階微分預處理后的預測集樣本光譜數(shù)據(jù)：式中，n′為預測集樣本數(shù)，m為光譜波點數(shù)，x是預測集樣本的atr反射率；

使用校正集pca算法中的線性變換系數(shù)矩陣p對預測集光譜信息進行壓縮：表示為un′×a＝xn′×mpm×a，取矩陣u的各行數(shù)據(jù)ii(i＝1,2,…,a)作為已建立好的預測模型的輸入，經(jīng)過預測模型的計算得到了模型輸出，即預測k值。通過公式(5)和公式(6)的計算，能夠得到模型的預測能力，作為評價模型性能的指標。

表3系統(tǒng)的將bpann改進模型和傳統(tǒng)的bpann以及主成分回歸(principalcomponentsregression，pcr)模型在thz光譜預測k值的性能方面進行比較。，可以看出光譜數(shù)據(jù)預處理的差分寬度選擇13(bpann、bpann改進模型)或者15(pcr模型)是合適的,差分寬度過小或者過大都會降低模型的預測性能。還可以看出非線性模型(bpann、bpann改進模型)的預測性能明顯優(yōu)于線性模型(pcr模型)。而且，bpann改進模型比bpann的k值預測性能有了提高。上述結(jié)果可以通過以下幾個方面來解釋：

表3三種k值預測模型的回歸結(jié)果比較

1、從thz光譜和豬肉k值的非線性趨勢上看，k值由atp的6種關(guān)聯(lián)物含量的比值決定，這些分子物質(zhì)的含量與thz光譜中的譜線數(shù)據(jù)是非線性關(guān)系。因此非線性模型的擬合效果會好于線性模型。

2、從豬肉變質(zhì)的化學機制上看，豬肉的變質(zhì)是一個復雜的化學過程，在若干種致腐細菌的作用下，肌肉中的蛋白質(zhì)先水解為多肽，再水解成氨基酸，并進一步分解為各種有機物質(zhì)?，F(xiàn)有技術(shù)的研究表明，蛋白質(zhì)、多肽和氨基酸等生物分子在thz波段都具有各自的特征吸收，其特征吸收主要來自于分子的集體振動模，所以thz光譜能夠反映這些變化。但是豬肉中的生物分子種類繁多，在常溫下各特征譜具有重疊性，導致thz光譜并不是僅僅表達了k值的變化信息，而是豬肉的化學成分復雜變化的綜合信息。因此，k值和thz光譜的關(guān)系應該是復雜非線性關(guān)系，而線性模型不能夠解釋這種復雜的回歸關(guān)系。

3、從建模算法的原理結(jié)構(gòu)上看，非線性模型比線性模型更善于自學習和自調(diào)整，bpann的網(wǎng)絡拓撲結(jié)構(gòu)更適合分析復雜的化學成分。當遇到復雜的回歸關(guān)系時，所建預測模型會更優(yōu)秀。bpann改進算法通過在迭代過程中合理整合弱預測器，使bpann弱預測模型逐漸提升為強預測模型。因此，bpann改進模型比bpann模型的預測性能有了提高。

2.5模型的實際意義

表3中的rp與rmsep指標表示模型對未知樣本k值預測的準確性。未知樣本，即未知其k值的樣本，不同于校正集和預測集的樣本，其k值是通過理化分析hplc方法獲得的。通過一階微分預處理以及pca光譜信息壓縮方法，將處理后的未知樣本的thz光譜數(shù)據(jù)代入模型，能夠快速計算出未知樣本的k值。本方案把校正集的樣本用于建模，預測集的樣本用于判斷模型優(yōu)劣，未知樣本的k值預測就是模型的實用化。模型的實際意義就是通過thz光譜快速檢測未知樣本的k值，而不是僅僅用已知樣本(包括校正集和預測集)建模與訓練。

3、結(jié)論

本文研究結(jié)果表明，采用atr全反射模式，在0.1thz～2thz范圍內(nèi)獲取冷鮮豬瘦肉表面的thz光譜數(shù)據(jù)，經(jīng)一階微分和sg平滑濾波處理后構(gòu)建的數(shù)學模型，能快速無損檢測豬肉k值，以評價豬肉新鮮度。通過對三種預測模型pcr、bpann、bpann改進模型的比較研究表明，非線性bpann和bpann改進模型比線性pcr模型能夠更好地擬合thz光譜數(shù)據(jù)和新鮮度k值之間的復雜非線性關(guān)系。而bpann改進模型在處理復雜關(guān)系模型上具有一定的優(yōu)勢，將bpann的預測性能提升到了rp為0.78，rmsep為13.45％。

基于相同的thz檢測原理，thz光譜分析法還能夠無損檢測其它肉類(如雞肉、牛肉、魚肉等)的k值。本發(fā)明為進一步基于此方法開發(fā)設計便攜式快速無損檢測設備提供了理論基礎。