基于HMM的食品安全風險預測方法研究

馬永軍，劉宇姍，侯陽陽，劉 洋

（天津科技大學計算機科學與信息工程學院，天津 300222）

基于HMM的食品安全風險預測方法研究

馬永軍，劉宇姍，侯陽陽，劉 洋

（天津科技大學計算機科學與信息工程學院，天津 300222）

以乳品為例，提出了一種基于隱馬爾可夫模型（HMM）的食品安全風險預測方法.依據HACCP對乳品供應鏈進行分析，找出供應鏈各個環節的關鍵控制點，并將其危害影響因素作為HMM的量化指標.分析時將供應鏈分為4個環節，每個環節最終的安全等級概率作為下一環節的初始狀態概率分布，利用HMM預測供應鏈風險等級和計算風險值，從而為供應鏈風險評價提供依據.

食品安全；HMM模型；乳品；供應鏈；風險預測

食品安全是關系國計民生的重大問題，近年來食品安全問題的預警和治理已受到高度關注.其中，風險預測作為風險預警和治理的基礎受到格外重視.風險預測是利用適當的風險預測工具和方法來確定安全風險等級和優先控制順序的過程，從而為實施有效的風險管理措施提供決策支持.

目前常用的食品安全風險預測方法主要分為定性和定量兩大類，其中定性的方法主要有灰色系統法、德爾菲法、層次分析法等.灰色系統法［1］是對不確定系統的研究方法；德爾菲法［2］能夠充分發揮各位專家的不同意見，評定過程公正；層次分析法［3］是食品安全風險分析與預測中應用最為廣泛的方法.但是它們的共性缺點就是在運算過程中受主觀因素影響較大，使得結果的精確度不夠，并且不能根據系統狀態的變化實時改變風險預測結果［4］.定量的分析方法主要有貝葉斯網絡法、人工神經網絡和支持向量機（SVM）等.貝葉斯網絡法［5］是一種基于概率推理的圖形化概率網絡，能在有限的不確定信息條件下進行學習和推理，但是貝葉斯網絡構造繁瑣，實際應用時還需反復交叉不斷完善，易用性不好.人工神經網絡［6］具有高度的自適應能力，但是學習速度較慢，而且算法容易陷入局部極值.SVM［7］本質上是將問題轉換為凸優化問題，可以保證找到全局極值，但在將低維空間向高維空間映射時，又存在計算耗時問題，因此不適合對于大樣本數據的分析處理.

食品安全分析和食品供應鏈密不可分，食品供應鏈主要由原料生產、生產加工、儲運和消費4個環節組成，食品安全風險存在于食品供應鏈的各個環節中，而且存在食品安全風險遷移問題，即一個環節的風險會向下一個環節傳遞.在目前的方法中，對風險的判斷一般是對4個環節的風險進行簡單平均或加權平均，未充分考慮食品安全風險的形成機理，割裂了4個環節之間的風險傳遞關系.

HMM模型是序列數據處理和統計學習的一種重要概率模型，具有建模簡單、數據計算量小、運行速度快等特點［8］，其在風險預測分析領域的應用較廣泛［9］.HMM模型被用來對網絡系統安全進行實時風險預測，表現出了很好的適用性和擴展性［10］.此外，HMM模型中的初始狀態分布矢量參數，對應初始狀態時各風險等級出現的概率，該參數可以有效表示食品供應鏈各環節對下一環節的影響程度.因此本文以乳品為例，利用HMM模型進行食品安全風險分析與預測.

1 HMM模型

HMM模型是在馬爾可夫鏈的基礎上發展起來的，在HMM中觀察到的事件是狀態的隨機函數，該模型是一個雙重隨機過程，即一個觀察狀態，一個隱藏狀態.

HMM的數學表達式為

其中：N為HMM中狀態的個數；M為HMM中對應的觀測值個數；π是初始狀態的分布矢量；A為狀態轉移概率矩陣；B是觀察向量的概率矩陣.

初始狀態的分布矢量可表示為

也可以簡單表示為

狀態轉移概率矩陣可表示為

元素aij表示從狀態i轉移到狀態j的概率，即

觀察向量的概率矩陣可表示為

元素bjk表示在j狀態的情況下，觀察狀態k出現的概率，即

HMM需要解決3個問題：

（1）評估問題：給定五元組模型λ=（N， M，π，A， B）和觀察序列O=（O1， O2，…，OT），計算這個觀察序列出現的概率.

（2）解碼問題：給定五元組模型λ=（N， M，π，A， B）和觀察序列O=（O1， O2，…，OT），求可能性最大的隱藏狀態序列.

（3）學習問題：也叫訓練問題.給定觀察值序列O=（O1， O2，…，OT），據此確定HMM模型λ=（N， M，π，A， B）.即如何調整λ=（N， M，π，A， B），使得P（ O| λ）最大.

上述的3個問題分別對應3個經典的解決算法［11-12］，分別是前向-后向算法、Viterbi算法、Baum-Welch算法.

2 食品安全風險預測方法

食品供應鏈主要包括原料生產、生產加工、儲藏運輸和消費4個環節，每個環節都會受很多因素的影響.假設影響因素為

其中XN（1≤N≤4）為風險預測的觀測值，XNi為第N個環節的第i個影響因素，通過一定的方法反映出風險預測值和其影響因素之間的關系，數學模型表達式為

HMM模型有2個隨機過程：可見的狀態序列，本文指供應鏈各個環節檢測到的數據；系統的真實狀態，本文指供應鏈各環節的真實風險值.供應鏈每個環節的HMM三元組參數模型包括狀態轉移概率矩陣A，觀測概率矩陣B和最初狀態分布π，字母表示為λ=（π， A， B）.除了第一個環節外，其他3個環節的參數模型都受上一環節的影響，本環節最終的風險概率直接影響下一環節的初始狀態概率分布.例如，原料生產環節通過HMM模型計算，最終得出各個安全狀態出現的概率分布為δ=｛δ1， δ2，…，δi｝（i為真實狀態的個數），那么下一環節，即生產加工環節的初始狀態分布概率π=δ.

設有A1—A5共5個狀態，表示供應鏈風險等級，其中A1表示正常安全狀態，供應鏈受威脅的概率為P=0；A2、A3、A4表示危險等級逐級加深，對應供應鏈受威脅的概率分別為0＜P≤0.2（低風險）、0.2＜P≤0.5（中級風險）、0.5＜P≤0.8（中高級風險）；A5表示重大危險狀況，供應鏈受威脅的概率為0.8＜P≤1.

這里的HMM為離散HMM，在建立離散的HMM模型系統時，可以假定：系統狀態包含系統的所有信息，并且在當前狀態下的觀察是獨立的.則HMM模型的狀態轉移模型可以表示為圖1.從一個節點移動到另一個節點，表示系統在源節點顯示的狀態，會轉變成目的節點的狀態，該模型圖是一個完全連接圖，表明任意安全狀態都有轉變為其他任意安全狀態的可能.

圖1 HMM狀態轉移模型Fig. 1 Transfer model of HMM state

據以上分析，提出圖2所示的風險預測流程.

圖2 風險預測流程Fig. 2 Risk assessment process

在用HMM模型進行風險預測時，首先對數據進行預處理，采用的歸一化映射為

其中x， y∈Rn.歸一化后的原始數據在［0，1］范圍內，即yi∈［0，1］，i=1，2，…，n .然后，將歸一化的數據通過DBSCAN聚類去除噪音點，得到初始數據.DBSCAN算法不需要事先知道要形成的簇類的數量，可以發現任意形狀的簇類，并且能夠識別出噪聲點，適用于類別數未知的聚類問題［13］.把聚類之后的數據通過Baum-Welch算法訓練得到對應的HMM三元組模型后，即可通過Viterbi算法對供應鏈進行相應的風險預測.

3 實驗分析

在食品安全方面，乳品一直受高度關注.所以本文以乳品為例進行分析.乳品在從農場到餐桌的過程經歷了4個環節：原奶生產、乳品加工、儲藏運輸和消費.HACCP［14］（hazard analysis critical control point）表示危害分析和關鍵控制點，它是對食品加工過程中可能發生的食源性危害進行識別和評估，進而采取一定控制的食品安全控制體系.利用HACCP確定關鍵控制點的判斷流程如圖3所示.

圖3 確定關鍵控制點的判斷樹Fig. 3 Critical control point decision tree

根據HACCP體系并結合乳品供應鏈的4個環節，對各個環節中的生物性、物理性和化學性的潛在危害進行分析，確定了關鍵控制點.以第一環節為例，原奶生產環節的關鍵點選擇見表1.

由于乳品一般都難以長久保持并容易受到微生物污染，供應鏈上任何環節的潛在危險因素都會影響食品質量安全.乳品加工環節是整個供應鏈最重要的部分，負責將收集到的原奶加工完成，儲運和消費環節的風險因素也基本來源于供應鏈的前幾個階段.因此，乳品生產供應鏈的各個環節緊密聯系，對每一個環節的檢查監管都要格外重視.最后一個環節的風險值即為整條供應鏈的風險大小.

表1 乳品原奶生產環節關鍵控制點分析Tab. 1 Critical control point analysis of raw material production

通過對天津某乳品企業的調研，采集大量數據，并建立了乳品追溯系統，利用HMM模型對乳品安全風險進行預測.

乳品供應鏈各環節的隱藏狀態有5個：A1—A5，分別代表5個安全等級.對歸一化的數據通過DBSCAN聚類進行聚類處理后，共得到4類，也就是系統的4個觀察狀態：V1—V4；利用HMM模型依次對4個環節進行參數評估和風險概率計算，每個環節最終的風險概率直接影響下一環節的初始狀態概率分布.

在訓練HMM模型時，利用Baum-Welch算法，通過讓概率P（O|λ）達到局部最大值得到HMM的參數模型.

定義供應鏈所有特征觀察值序列前向變量［11］為

由式（8）可知，供應鏈所有特征觀察值序列概率為

通過式（9）求取HMM的三元組參數模型λ=（π， A， B），再由重估公式［15］可得

輸出概率P（O|λ）會隨著重估次數的增加越來越大，直至參數πi、aij、bij收斂或算法達到停止條件為止.

計算得到第一環節的初始狀態分布概率為π=｛0.4612，0.1132，0.2139，0.0547，0.1570｝，狀態轉移矩陣為

觀測向量概率矩陣為

有了HMM的三元組模型，就可以對相關數據進行解碼分析.在此給定一組觀察序列｛O2， O1， O3， O4，O4， O3， O2， O2， O1， O1， O4， O3｝，觀察序列的狀態分布如圖4所示.

圖4 觀察狀態序列Fig. 4 Sequence of observed state

繼而利用Viterbi算法求解觀察序列的最優路徑，得到供應鏈最有可能的真實狀態.首先將πi與所有特征觀察值轉移概率bi相乘，得到初始化路徑δt（ i）.

將初始化的路徑δt（ i）與狀態轉移概率矩陣的元素aij相乘，取最大值與所有特征觀察值轉移概率bj相乘，得到當前路徑最大值δt（j）.

根據式（14）、（15）結果，可求出對應的食品安全各環節真實狀態的最大概率P?和最大概率對應狀態序列qt?分別為

因此，可以評估出該環節對應的真實狀態｛T1， T2， T1， T2， T4， T4， T1， T2， T4， T2， T5， T2｝.

真實狀態序列的分布如圖5所示.

圖5 真實狀態序列Fig. 5 Sequence of true state

由結果可知，系統在t=5，6，9，11時的風險級別較大，在t=1，3，7時刻，風險級別較小.通過Viterbi算法可以計算出原奶生產環節在最后時刻處在各狀態的概率為δ=｛0.1672，0.5342，0.1745，0.0723，0.0518｝，則在乳品生產環節中，初始狀態分布概率就為π=｛0.1672，0.5342，0.1745，0.0723，0.0518｝.

參照對第一環節的計算過程，依次利用Baum-Welch算法對各環節進行參數計算，然后利用Viterbi算法計算系統的真實狀態.

最后，計算得到最終環節——消費環節各安全等級的概率為π=｛0.1547，0.6073，0.0645，0.1043，0.0692｝.

假設在A1-A5不同安全等級的情況下，系統對應的開銷，即風險對系統的影響分值為C=｛0，5，10，15，20｝.參考風險評價模型的典型代表OCTAVE［10］給出一個計算風險值的公式

其中：Rt表示在t時刻系統所處的總體風險值；αt（ i）為t時刻系統處在安全狀態Ai的概率；N是安全狀態的數目；c（ i）是與狀態Ai關聯的開銷.

根據式（19），可以計算出乳品的最終風險值R= 6.049，5.

因此，可以利用上述方法對乳品各環節的不同時刻進行風險預測與開銷分析，有助于生產者和決策者更全面了解乳品供應鏈各環節的風險態勢，及時查出問題發生的緣由，根據已有流程或者規定采取相應措施，并不斷更新和完善風險規避的方法，最大程度減少乳品安全對群眾的影響.

4 結 語

本文利用HMM模型對乳品質量安全風險進行分析與預測.該方法根據食品安全風險的形成機理，考慮了4個環節之間的風險傳遞關系以及各個環節風險的動態性，結合定性與定量評價，可以更精確地描述乳品風險的大小.該風險預測模型除了可用于對乳品進行評估外，還可以擴展應用到食品安全的其他領域.今后還需通過更多實驗對各環節風險傳遞的影響因子進行優化，改進HMM模型的參數訓練方法，以提高參數模型的訓練速度和評估準確度.

［1］ 黃秋婷，丁怡，邱佩麗. 基于灰色殘差修正模型的食品安全監測預警分析方法研究［J］. 現代農業科技，2013（5）：304-305.

［2］ Soon J M，Davies W P，Chadd S A，et al. A Delphi-based approach to developing and validating a farm food safety risk assessment tool by experts［J］. Expert Systems with Applications，2012，39（9）：8325-8336.

［3］ 珺安. 基于層次分析法的乳品質量安全預警系統研究［D］. 哈爾濱：東北農業大學，2012.

［4］ 孔大維. 基于供應鏈管理的食品可追溯系統的構建研究［D］. 成都：成都理工大學，2013.

［5］ 張麗，滕飛，王鵬. 基于貝葉斯網絡的食品供應鏈風險評價研究［J］. 食品研究與開發，2014（18）：179-182.［6］ Xiao J，Ma Z S，Zhang D J. The research on the BP neural network in the food security risk early warning applied under supply chain environment［C］// Proceedings of Second International Symposium on Knowledge Acquisition and Modeling. Piscataway：IEEE，2009：281-284.

［7］ 曹秋勤. 基于支持向量機的蔬菜質量安全預測及溯源模型的研究與應用［D］. 廣州：華南理工大學，2014.

［8］ 王志堂，蔡淋波. HMM模型（HMM）及其應用［J］. 湖南科技學院學報，2009，30（4）：42-44.

［9］ 張建華. 基于深度學習的語音識別應用研究［D］. 北京：北京郵電大學，2015.

［10］ 龍百元. 基于層次分析法和HMM模型的風險預測方法研究［D］. 長沙：湖南大學，2007.

［11］ Ieee L R R F. A tutorial on hidden Markov models and selected applications in speech recognition［J］. Proceedings of the IEEE，1989，77（2）：257-286.

［12］ 謝錦輝. 隱Markov模型（HMM）及其在語音處理中的應用［M］. 武漢：華中理工大學出版社，1995.

［13］ 慕升弟. 基于協議HMM模型的網絡流量識別技術研究［D］. 重慶：重慶大學，2013.

［14］ CAC. Draft Principle and application of the hazard analysis and critical control point system：Report of twenty-fifth session of the codex committee on food hygiene［R］. Washington DC：CAC，1991.

［15］ Schulz L L，Tonsor G T. Cow-calf producer perceptions regarding individual animal traceability［J］. Journal of Agricultural and Applied Economics，2010，42（4）：659-677.

責任編輯：常濤

Food Safety Risk Prediction Based on HMM

MA Yongjun，LIU Yushan，HOU Yangyang，LIU Yang
（College of Computer Science and Information Engineering，Tianjin University of Science & Technology，Tianjin 300222，China）

Taking dairy products as the research subjects，a food safety risk assessment method based on the Hidden Markov Model（HMM） was developed.The critical control point for each sector of the supply chain was found out with the method of hazard analysis critical control point.Those critical points can be used as quantitative indicators of HMM.According to the characteristics HMM model，the supply chain was divided into four sectors.The ultimate probability of each sector is the initial state probability distribution of the next part.Finally the level of risk in the supply chain and risk level can be calculated with HMM，which can provide a reference for supply chain risk assessment.

food safety；HMM；dairy products；supply chain；risk prediction

TP311

A

1672-6510（2016）05-0069-05

10.13364/j.issn.1672-6510.20150209

2015-11-16；

2016-04-18

教育部規劃項目（12YJAZH091）

馬永軍（1970-），男，吉林長春人，教授，yjma@tust.edu.cn.