基于改進(jìn)HMM的食品安全風(fēng)險(xiǎn)評(píng)估方法

熊 慧 唐宏亮 - 丁 永

(1. 武漢交通職業(yè)學(xué)院，湖北 武漢 430065；2. 湖北第二師范學(xué)院，湖北 武漢 430205；3. 江西理工大學(xué)，江西 贛州 341000)

目前的食品安全風(fēng)險(xiǎn)評(píng)估是食品檢驗(yàn)檢測(cè)部門對(duì)食品質(zhì)量進(jìn)行的抽樣檢測(cè)，僅評(píng)價(jià)食品是否合格，不能深入挖掘食品安全信息下更深層次的風(fēng)險(xiǎn)隱患。風(fēng)險(xiǎn)評(píng)估是用于估計(jì)風(fēng)險(xiǎn)概率和相關(guān)不確定性的過程，挖掘潛在危險(xiǎn)和相關(guān)風(fēng)險(xiǎn)。由于食品安全問題不僅會(huì)對(duì)人體造成傷害，還會(huì)影響社會(huì)的穩(wěn)定發(fā)展。因此有必要對(duì)食品安全風(fēng)險(xiǎn)評(píng)估方法進(jìn)行研究。

王小藝等[1]將蒙特卡羅仿真用于食品風(fēng)險(xiǎn)評(píng)估方法，構(gòu)建菌落動(dòng)態(tài)生長模型，并對(duì)風(fēng)險(xiǎn)進(jìn)行量化分析。結(jié)果表明，該風(fēng)險(xiǎn)評(píng)估方法可以為風(fēng)險(xiǎn)管理和決策部門提供較為直觀的風(fēng)險(xiǎn)評(píng)估結(jié)果。高亞男等[2]提出了一種先驗(yàn)風(fēng)險(xiǎn)概率與模糊層次劃分相結(jié)合的方法，計(jì)算不同類型屬性的模糊綜合風(fēng)險(xiǎn)值作為預(yù)測(cè)模型的標(biāo)簽值。結(jié)果表明，該方法可以較為準(zhǔn)確地預(yù)測(cè)風(fēng)險(xiǎn)值，為決策者提供有價(jià)值的決策依據(jù)。張利鋒等[3]根據(jù)國家相關(guān)標(biāo)準(zhǔn)程序，結(jié)合食物攝入量，對(duì)17種抗生素的殘留量進(jìn)行采集，并采用點(diǎn)評(píng)估法對(duì)雞肉和雞蛋飲食中接觸抗生素的風(fēng)險(xiǎn)進(jìn)行評(píng)估。結(jié)果顯示，抗生素總體檢出率為14.4%，不合格率為6.76%。雖然雞肉和雞蛋含有一定的抗生素殘留，但暴漏的風(fēng)險(xiǎn)非常低，對(duì)人類健康的危害幾乎為0。然而，上述研究中的食品安全風(fēng)險(xiǎn)評(píng)估方法精度低且不準(zhǔn)確，需要人的主觀判斷，缺乏客觀性，無法結(jié)合時(shí)間特征進(jìn)行動(dòng)態(tài)評(píng)估，適應(yīng)能力有待進(jìn)一步提高。

布谷鳥搜索算法具有良好的全局搜索功能和跳出局部?jī)?yōu)化的能力；隱馬爾可夫模型應(yīng)用在評(píng)估中，結(jié)合數(shù)據(jù)時(shí)間特征，使評(píng)估結(jié)果具有更好的實(shí)時(shí)性和針對(duì)性。文章擬提出一種通過布谷鳥搜索(Cuckoo Search，CS)算法優(yōu)化隱馬爾可夫模型(Hidden Markov Model，HMM)的食品風(fēng)險(xiǎn)評(píng)估方法，以期為食品評(píng)估技術(shù)的發(fā)展提供參考。

1 食品安全風(fēng)險(xiǎn)指標(biāo)構(gòu)建

由于檢測(cè)指標(biāo)多且食品數(shù)據(jù)復(fù)雜度高，采用灰色關(guān)聯(lián)分析算法優(yōu)化的解釋結(jié)構(gòu)模型進(jìn)行多層次劃分[4]。以滅菌乳檢測(cè)指標(biāo)為例。其中，所使用的滅菌乳數(shù)據(jù)有9個(gè)檢測(cè)指標(biāo)，分為4個(gè)層次。分層的結(jié)果如圖1所示。第1層權(quán)重最高，第4層權(quán)重最低。表1各指標(biāo)的權(quán)重值，然后將其與檢測(cè)指標(biāo)對(duì)應(yīng)的相關(guān)系數(shù)平均值相乘，得到檢測(cè)指標(biāo)的權(quán)重結(jié)果。

文中通過與國家標(biāo)準(zhǔn)的檢測(cè)差值評(píng)估食品安全風(fēng)險(xiǎn)，將對(duì)身體有益物質(zhì)(蛋白質(zhì)、脂肪和非脂乳固體)定為正指標(biāo)，含量必須高于國家標(biāo)準(zhǔn)。將對(duì)身體有害物質(zhì)(汞、鉛、砷、鎘和黃曲霉毒素)定為負(fù)指標(biāo)，含量必須低于國家標(biāo)準(zhǔn)[5]。將酸度定為中性指標(biāo)，介于國家標(biāo)準(zhǔn)范圍內(nèi)。根據(jù)5-標(biāo)度法，安全風(fēng)險(xiǎn)數(shù)據(jù)分為5類(1～5)。值越高，質(zhì)量安全風(fēng)險(xiǎn)越高。文中第1類計(jì)算結(jié)果為[-0.077 6，0.060 96]，第2類計(jì)算結(jié)果為(0.060 96，0.199 52]，第3類計(jì)算結(jié)果(0.199 52，0.338 08]，第4類計(jì)算結(jié)果(0.338 08，0.476 64]，第5類計(jì)算結(jié)果(0.476 64，0.615 20]。

圖1 分層結(jié)果Figure 1 System structure

表1 權(quán)重結(jié)果Table 1 Optimal parameters

2 風(fēng)險(xiǎn)評(píng)估方法

2.1 隱馬爾可夫模型

隱馬爾可夫模型由Baum等提出，主要由兩個(gè)過程組成，馬爾可夫鏈和觀測(cè)過程，分別描述了狀態(tài)轉(zhuǎn)移過程和觀測(cè)概率矩陣[6]。隱馬爾可夫模型可以被描述為一個(gè)五元素表達(dá)式λ=(Q,O,π,A,B)，各元素如式(1)～式(5)所示[7]。

Q={q1,q2,…,qN}，

(1)

O={o1,o2,…,oM}，

(2)

π={[π(i)N]π(iq)=P(qi)}，

(3)

A=π[aij]N*N,aij=P(qi|qj)}，

(4)

B=π[bj(k)]N*M,bj(k)=P(ok|qj)}，

(5)

式中：

Q、O——無模型的隱含和觀測(cè)狀態(tài)；

π——初始狀態(tài)概率矩陣；

A、B——無狀態(tài)轉(zhuǎn)移和觀測(cè)狀態(tài)的轉(zhuǎn)移概率矩陣；

N——模型中隱藏層狀態(tài)的數(shù)量；

M——觀察特征的數(shù)量；

q——隱含序列狀態(tài)；

o——觀測(cè)序列狀態(tài)。

CHH主要解決基本問題：

學(xué)習(xí)問題：利用Baum-Welch算法對(duì)觀測(cè)序列O進(jìn)行訓(xùn)練，得到條件概率P(O,λ)最大的模型λ[8]。

評(píng)估問題：利用正反向算法計(jì)算觀測(cè)序列O在模型λ中出現(xiàn)的概率P(O,λ)[9]。

預(yù)測(cè)問題：已知模型λ和觀測(cè)序列O，使用Viterbi算法找到序列中發(fā)生概率最大的狀態(tài)。

(1) 前后向算法：為了解決評(píng)估問題，采用前后向算法求解P(O,λ)的值，這是解決隨機(jī)推理問題最常用的方法。由兩部分組成：前向和后向。

前向概率：時(shí)刻t的隱含狀態(tài)為qi，觀測(cè)狀態(tài)為o1,o2,…,ot的概率，如式(6)所示[8]。

at(i)=P(o1,o2,…,ot,it=qi|λ)。

(6)

計(jì)算最終結(jié)果如式(7)所示。

(7)

后向概率：時(shí)刻t的隱含狀態(tài)qi，時(shí)刻t+1到時(shí)刻T觀測(cè)狀態(tài)為ot+1,ot+2,…,oT的概率，如式(8)所示[10]。

βt(i)=P(ot+1,ot+2,…,oT,it=qi|λ)。

(8)

計(jì)算最終結(jié)果如式(9)所示。

(9)

給定模型λ和觀測(cè)序列O，在時(shí)刻t處于狀態(tài)qi的概率表示如式(10)所示。

(10)

模型λ和觀測(cè)序列O已知，在時(shí)刻t的狀態(tài)qi，且時(shí)刻t+1處于狀態(tài)qj的概率如式(11)所示[11]。

(11)

(2) Baum-Welch算法：Baum-Welch算法在解決學(xué)習(xí)問題中應(yīng)用廣泛，步驟為初始化、參數(shù)更新和終止[12]。隨機(jī)初始化模型參數(shù)πi、aij、bj(k)，并通過正向和反向算法計(jì)算γt(i)和ξt(i,j)。參數(shù)計(jì)算如式(12)～式(14)所示。

πi=γ1(i)，

(12)

(13)

(14)

直到參數(shù)收斂，輸出最終結(jié)果。

(3) Viterbi算法：在Viterbi算法中，求解最優(yōu)狀態(tài)序列的主要步驟類似于前后向算法[13]。解析過程是狀態(tài)序列的初始化、遞歸、終止和獲取狀態(tài)序列。在初始化過程中，可以根據(jù)變量定義獲得初始化條件。當(dāng)時(shí)間達(dá)到最大值時(shí)，Viterbi算法終止，同時(shí)解得P*的最大值。

2.2 改進(jìn)隱馬爾可夫模型

HMM模型中Baum-Welch算法對(duì)初值有很強(qiáng)的依賴性，因?yàn)槲闹型ㄟ^CS算法優(yōu)化的HMM初始參數(shù)，在通過Baum-Welch算法進(jìn)行局部校正，使其快速收斂到全局最優(yōu)解[14]。CS算法是一種基于杜鵑寄生育雛行的自然啟發(fā)式算法，具有良好的全局搜索功能和跳出局部?jī)?yōu)化的能力，滿足全局收斂的要求。

布谷鳥搜索算法優(yōu)化隱馬爾可夫模型的過程：

步驟1：數(shù)據(jù)預(yù)處理。采用灰色關(guān)聯(lián)分析算法和解釋結(jié)構(gòu)模型對(duì)數(shù)據(jù)進(jìn)行預(yù)處理。

步驟2：通過CS算法優(yōu)化的HMM初始參數(shù)。

步驟3：采用Baum-Welch算法進(jìn)行局部校正，以提高風(fēng)險(xiǎn)評(píng)估算法的準(zhǔn)確性[15]。

步驟4：使用Viterbi算法完成風(fēng)險(xiǎn)分類。

圖2為CS算法優(yōu)化HMM的流程圖。

3 試驗(yàn)結(jié)果與分析

3.1 試驗(yàn)參數(shù)

為了對(duì)CS-HMM方法的優(yōu)越性進(jìn)行驗(yàn)證，將改進(jìn)的HMM方法與未改進(jìn)前的HMM方法進(jìn)行了比較。利用2018年11月—2021年2月951條滅菌乳的檢測(cè)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)驗(yàn)證。該設(shè)備是聯(lián)想PC，操作系統(tǒng)是windows 10 64位旗艦，Intel i52450m CPU，頻率2.5 GHz，8 GB內(nèi)存[16]。

圖2 改進(jìn)隱馬爾科夫模型流程圖Figure 2 Flow chart of improved hidden markov mode

風(fēng)險(xiǎn)指數(shù)作為訓(xùn)練數(shù)據(jù)，布谷鳥搜索算法用于全局搜索，找出HMM模型的最優(yōu)初始值。計(jì)算HMM風(fēng)險(xiǎn)評(píng)估模型的參數(shù)A、B、π，如式(15)～式(17)所示。

(15)

(16)

(17)

采用布谷鳥搜索算法計(jì)算的全局最優(yōu)解作為Baum-Welch算法的初值進(jìn)行精細(xì)搜索。以300個(gè)訓(xùn)練集數(shù)據(jù)為例，訓(xùn)練結(jié)果如式(18)～式(20)所示。

(18)

(19)

(20)

3.2 試驗(yàn)分析

基于食品安全風(fēng)險(xiǎn)評(píng)估模型的評(píng)估標(biāo)準(zhǔn)，通過誤差矩陣分析模型的準(zhǔn)確性，將樣本分為正樣本TP和TN。當(dāng)判斷準(zhǔn)確度時(shí)，各類別單獨(dú)為正，其余均為負(fù)，TPR為真正率，TNR為真負(fù)率。初始風(fēng)險(xiǎn)閾值設(shè)置為0，并在0.2與1.0之間逐步增加。圖3為2個(gè)不同大小的訓(xùn)練集評(píng)估準(zhǔn)確度結(jié)果。

從圖3可以看出，不同訓(xùn)練集模型準(zhǔn)確度隨閾值增加逐漸降低，當(dāng)閾值達(dá)到0.6時(shí)，下降率顯著增加。為了保證評(píng)估的準(zhǔn)確性，選擇閾值為0.6時(shí)的TPR和TNR。

分了分析數(shù)據(jù)集對(duì)評(píng)估結(jié)果影響，在閾值0.6時(shí)分析不同數(shù)據(jù)集對(duì)測(cè)試準(zhǔn)確度的影響，訓(xùn)練集從200～375，次增加25，測(cè)試集為全部數(shù)據(jù)。圖4為準(zhǔn)確度隨訓(xùn)練集變化曲線。

從圖4可以看出，模型準(zhǔn)確度最高出現(xiàn)在訓(xùn)練集為300時(shí)，此時(shí)TPR和TNR分別達(dá)到97.79%和99.47%，所以300個(gè)數(shù)據(jù)足以滿足模型訓(xùn)練。

為了比較該模型與傳統(tǒng)HMM模型的性能，在閾值為0.6時(shí)分析不同模型隨訓(xùn)練集的變化趨勢(shì)，結(jié)果如圖5所示。

由圖5可知，文中模型TPR1和TNR1在訓(xùn)練集大小為300時(shí)達(dá)到峰值，TPR1和TNR1值分別為97.79%和99.47%。傳統(tǒng)HMM模型在350時(shí)達(dá)到峰值，TPR2和TNR2值分別為81.49%和96.84%。文中模型的準(zhǔn)確度峰值優(yōu)于未改進(jìn)前，需要的訓(xùn)練集也較少，證明了文中模型的優(yōu)越性。

通過vitabi算法對(duì)訓(xùn)練后的模型進(jìn)行風(fēng)險(xiǎn)評(píng)估，如圖6所示模型改進(jìn)前后評(píng)估結(jié)果比較。

TPR(250)和TNR(250)表示訓(xùn)練集大小為250時(shí)模型的真正率和真負(fù)率，TPR(300)和TNR(300)未訓(xùn)練集大小為300時(shí)模型的真正率和真負(fù)率

圖4 不同訓(xùn)練集模型準(zhǔn)確度Figure 4 Model accuracy of different training sets

TPR1和TNR1、TPR2和TNR2分別為文中模型和HMM訓(xùn)練時(shí)結(jié)果的真正率和真負(fù)率

圖6 改進(jìn)前后評(píng)估結(jié)果比較Figure 6 Comparison of evaluation results beforeand after improvement

從圖6可以看出，與傳統(tǒng)HMM模型相比，文中模型的評(píng)估結(jié)果更接近真實(shí)值。這是因?yàn)槲闹心Ｐ筒捎肅S優(yōu)化HMM模型，使參數(shù)快速收斂到全局最優(yōu)解，且結(jié)合時(shí)間特征進(jìn)行預(yù)測(cè)，可以有效提高企業(yè)的食品安全質(zhì)量。

4 結(jié)論

文中提出了一種將隱馬爾可夫模型和布谷鳥搜索算法結(jié)合用于食品風(fēng)險(xiǎn)評(píng)估方法。通過布谷鳥搜索算法搜索全局最優(yōu)解作為隱馬爾可夫模型初值，在通過Baum-Welch算法進(jìn)行局部校正，使其快速收斂到全局最優(yōu)解。結(jié)果表明，該風(fēng)險(xiǎn)評(píng)估方法比未改進(jìn)前更準(zhǔn)確、有效。在訓(xùn)練集為300時(shí)，測(cè)試數(shù)據(jù)真正值和真負(fù)值的準(zhǔn)確度分別為97.79%和99.47%。該模型的食品安全風(fēng)險(xiǎn)評(píng)估尚處于試驗(yàn)階段，布谷鳥搜索算法在搜索全局最優(yōu)解過程中存在收斂速度慢的問題，后續(xù)工作中應(yīng)該改進(jìn)布谷鳥搜索算法的性能，如自適應(yīng)布谷鳥搜索算法、合作協(xié)同進(jìn)化布谷鳥搜索算法等，爭(zhēng)取盡快實(shí)際應(yīng)用，在實(shí)際應(yīng)用中不斷完善模型。