基于蒙特卡洛仿真的小麥倉儲環節菌落總數風險評估方法

王小藝，陳 謙，趙峙堯,*，熊 科，史 策，裴鵬鋼

（1.北京工商大學計算機與信息工程學院，北京 100048；2.北京工商大學，中國輕工業工業互聯網與大數據重點實驗室，北京 100048；3.北京工商大學，北京市食品添加劑工程技術研究中心，北京 100048；4.北京市農業信息技術研究中心，北京 100097）

小麥是世界第二大糧食作物，也是中國繼玉米和大米之后的第三大作物。20世紀80年代初以來，全世界的小麥產量顯著增加，中國每年的小麥產量超過1.2億 t，已成為世界上最大的小麥生產國。然而，近年來小麥食品安全問題層出不窮并引起了社會各界的廣泛關注[1-4]。風險評估是一種結構化的科學過程，用于估計風險的概率和嚴重程度以及隨之而來的不確定性，以確定因接觸食物中的生物、化學或物理危害物而對健康生活造成的潛在危害和相關風險[5-6]。國際食品安全組織正大力推動將風險評估技術應用于評估食品中危害物安全的問題[7]。目前的中國食品安全法要求建立國家食品安全風險評估制度，以評估中國食品中的危害物風險[8]。

現有的食品安全風險評估方法主要分為兩類：定性風險評估方法和定量風險評估方法[9]。定性風險評估是一種基于專家經驗的風險評估方法，主要是根據評估專家的經驗和知識對風險進行分析判斷，并據此來描述評估結果[10]。目前已有的定性風險評估方法包括德爾菲法、決策樹、危害分析臨界控制點（hazard analysis critical control point，HACCP）等[11]。這些方法多被應用于食品供應鏈中物流、經濟、信息等因素的風險管理優化，通過建立相應的風險指標體系評價食品安全風險水平。定量評估方法是指依據理化實驗抽檢或數據驅動模型獲取相關定量數據，并通過分析計算出風險指標的量化值，以描述食品安全風險水平[10-13]。其中，理化實驗抽檢方法包括色譜法、聚合酶鏈式反應法、光學傳感器法等[12]；常用的數據驅動方法包括灰色關聯度分析、支持向量機、蒙特卡洛仿真等[13-14]。這些方法多用于危害物的風險定量評估。結合已有研究分析[15-16]，定性風險評估方法大多依賴風險評估者的主觀經驗，受專家自身素質的影響大。然而，由于現代食品安全監控數據具有維度高、隨機性強等特性，該類方法難以滿足風險評估的精度要求。對于定量風險評估方法，理化實驗抽檢方法一般對實驗設備及實驗環境要求較高，實驗周期較長，相關學者致力于快速無損檢測技術的研究，并取得了一定成效，但單點檢測無法滿足食品安全預警要求，數據驅動模型的評估預測能力強弱更多取決于數據的準確性高低以及體量大小。同時，已有的食品安全定量風險評估缺少直觀明確的度量指標，大多數已有的食品安全風險評估方法都是通過危害物含量直接量化風險，模型演化過程的隨機性及環境的不確定性會導致危害物含量的測量精度較低。

考慮到定性和定量風險評估方法各自的優缺點和風險度量指標研究的不足，本實驗以小麥倉儲環節中菌落總數為評估對象，研究一種針對生物性危害物的風險評估方法。首先，建立可以適應不同生長環境的菌落動態生長模型；在此基礎上，通過蒙特卡洛仿真對小麥倉儲環節中菌落總數概率分布進行預測；最后，提出風險度量指標“危害度”，并基于該指標對小麥倉儲環節中菌落總數風險進行定量評價。該風險評估方法可以為風險管理和決策部門提供實時、直觀的風險評估結果，為風險防控提供理論支持。

1 菌落動態生長模型構建

在微生物生長允許的環境條件下，微生物生長絕大部分需經歷遲滯期、指數期和穩定期[17-18]。預測微生物學中常使用數學機理模型來描述微生物的動態變化，這些機理模型一般可分為初級模型、二級模型和三級模型[19-22]：初級模型用來描述在恒定環境下，微生物隨時間的生長狀況，其動力學參數決定了特定微生物在特定恒定環境下的生長特性；二級模型用來模擬環境條件對微生物動力學參數的影響，例如溫度和相對濕度對初級模型中生長速率和遲滯期的影響；三級模型是通過多種一級模型與二級模型組合成的系統軟件，來描述動態環境下不同微生物的生長狀況[23]。本實驗中菌落屬于一類微生物，其含量用菌落總數表示，可用來評估食品被有害菌落污染的程度及其衛生狀況。

對于菌落的動態生長，引入采樣時刻k＝{h,2h,…,K}，其中h為采樣周期/d，K為采樣時長/d，菌落動態生長模型按公式（1）[20-21]計算。

式中：Nk和Nk－h分別為k和k－h時刻的菌落總數/（CFU/g）；f[k－h,k]為在[k－h,k]時間段內菌落生長數量增量/（CFU/g）。

本實驗中f[k－h,k]通過菌落生長Logistic初級模型[21]獲得，如公式（2）所示。

結合公式（1）、（2），建立菌落動態生長模型，如公式（3）所示。

式中：λk受溫度T/℃和濕度αw變化的影響，可通過二階多項式函數（二級模型）描述[24]，C0～C5為函數參數；μk受溫度T和相對濕度αw變化的影響，可通過多因子基數模型（二級模型）描述[25]，τk和ρk分別為多因子基數模型中溫度因子和濕度因子，μopt為菌落生長比生長速率最適值；Tmin、Tmax分別為菌落生長溫度區間最小值和最大值，Topt和αw,opt分別為菌落生長最適溫度和最適相對濕度，αw,min為菌落生長濕度區間最小值，Tk和αw,k分別為k時刻的溫度和相對濕度，它們的數值變化可通過環境演化模型fT和fαw描述。

定義菌落動態生長模型的菌落生長指標x根據公式（4）計算。

假設以上模型參數θ為慢變量（近似恒定），公式（3）可簡化為公式（6）。

式中：f為菌落動態生長模型函數；ωx、ωθ分別為菌落生長指標和模型參數變化中受到的高斯噪聲干擾[26]，分別滿足ωx,k∶N（0,Σx）、ωθ,k∶N（0,Σθ），k∈[0,K]，∑x、∑θ為對應的噪聲協方差陣。

需要注意的是，不同菌落生長特性不同，同種菌落的生長環境也會有差別，因此，不同菌落動態生長模型參數也不同。為使模型更為精確、合理，模型參數θ的取值需結合實測數據和相關參數估計方法確定。

2 小麥倉儲環節中菌落總數的概率分布預測

2.1 蒙特卡洛仿真

蒙特卡洛仿真是以中心極限定理和大數定律為基礎的概率統計方法，通過對模型演化進行規律性模擬，以估計參數的不確定性。根據中心極限定理，可以從有限的隨機樣本群中隨機抽取服從特定分布規律的樣本，根據大數定律，可以通過統計反復實驗過程中某事件發生的頻率來近似地獲取隨機事件的概率[27]。基于以上兩點定理，蒙特卡洛仿真通常適用于優化、積分、隨機概率分布等數學問題中的解析無法被精確計算的情景[28-29]。

一般的蒙特卡洛仿真過程為：1）確定模型的參數變量及基本特征；2）根據問題和實測數據構造隨機模型；3）根據模型中各個隨機變量的分布產生隨機數并進行多次隨機采樣；4）通過多次模擬獲取實驗結果并統計分析獲取問題的概率解以及解的精度估計[30]。

2.2 基于蒙特卡洛仿真的小麥倉儲環節中菌落總數的概率分布預測

考慮到蒙特卡洛仿真作為定量風險評估方法在處理微生物生長演化過程的隨機性及環境不確定性問題上的優勢[31]，本實驗基于蒙特卡洛仿真對菌落總數概率分布進行預測。

首先，假設小麥倉儲環境為恒溫恒濕，則公式（3）中菌落動態生長環境演化模型可用公式（7）、（8）表示。

式中：Tset為設置的恒定溫度/℃；αw,set為設置的恒定相對濕度；ωT和ωαw分別為環境溫度和濕度演化模型的高斯噪聲干擾，分別滿足ωT,∶kN（0,QT）、ωαw∶,k N（0,Qαw），k∈[0,K]，QT、Qαw為對應的噪聲協方差陣。

基于以上設置，可利用蒙特卡洛仿真對菌落總數演化過程進行模擬。針對蒙特卡洛仿真中的每一次模擬，同一指標在同一預測時間點的數值均不相同，這些數值同時具有規律性和隨機性。規律性體現在預測值是基于確定的菌落動態生長模型、確定的菌落生長指標和模型參數初值產生的；隨機性體現在每一個菌落生長指標在每一個預測時間點均受噪聲影響，導致得到的預測值不盡相同。當模擬次數足夠大時，結合概率統計規律，可以獲得各菌落生長指標在不同時間點上取值的概率分布。基于蒙特卡洛仿真的小麥倉儲環節中菌落總數概率分布預測的具體步驟為：

1）設置菌落生長指標初值x0，模型參數初值θ0，蒙特卡洛仿真所用粒子數D/個；

2）對第d∈[1,D]個粒子，在k＝0時刻，按照公式（9）賦予初值；

3）對于時刻k∈[0,K]，基于k－h時刻的菌落生長指標和模型參數，根據公式（6）進行單步預測，分別得到；

4）若k≤K，k←k＋h，返回步驟3，否則進行步驟5；

5）若d≤D，d←d＋1，返回步驟2，否則進行步驟6；

6）對k∈[0,K]，按照公式（10）計算菌落總數的概率密度函數P（xk|（x0,θ0））；

式中：P（xk|（x0,θ0））為基于初始值（x0,θ0）的條件概率密度函數；δ（·）表示狄拉克函數。

大數據已經發展成社會和時代發展的主要屬性，伴隨國家“互聯網+”策略的落實，在全新的階段，需要將大數據置于重點位置，這也是必然的變革趨勢[1]。劉延東副總理在首屆國際教育信息化大會上明確強調了“互聯網+”戰略的重要地位，給出的倡議為:需要更加注重教育領域的信息化發展，利用創新技術來推動教學工作的發展，確保受教育者能夠公平地享用信息技術，并為不同文明的交流提供有利條件。

3 基于“危害度”的小麥倉儲環節中菌落總數風險評估

在系統工程領域，相關學者基于可靠性理論中性能可靠度的定義和安全性分析中安全概率的定義[32-33]，提出了“健康度”的概念，并以此作為評估系統工作性能及表現的度量指標[34-35]。本實驗參考“健康度”的思想提出了新的度量指標“危害度”，并基于該指標對小麥倉儲環節中菌落總數風險進行定量評價。

對于一般動態系統，其風險指標所在的n維空間Rn可劃分成一個無風險空間SH和一個有風險空間SF（SHUSF＝Rn）。對于某一時刻k，系統的危害度Rk可按公式（11）計算。

式中：Rk為k時刻該動態系統處于有風險狀態的概率，也稱危害度；對于小麥倉儲環節的菌落總數而言，小麥倉儲環節中菌落生長指標x的演化過程可以看成一個動態系統。

結合公式（10），小麥倉儲環節中菌落總數“危害度”可通過公式（12）計算。

式中：SF＝{xk|Nk＞Θ,Nk∈xk}，Θ為菌落總數超標閾值/（CFU/g）；Rc,k∈[0,1]，Rc,k＝0表示k時刻小麥倉儲環節中菌落總數無超標風險；Rc,k＝1表示k時刻小麥倉儲環節中菌落總數存在超標風險。

4 案例驗證

以定量菌落總數為例，基于MATLAB軟件并通過實測數據演化曲線對比來驗證小麥倉儲環節中菌落總數概率分布預測及風險評估方法的有效性。其中，實測數據演化曲線來源于文獻[36]，其研究對象為小麥在儲藏過程中的霉菌活動特性，包括以蠕孢霉和鏈格孢霉為主的混合霉菌菌落總數演化特性。

4.1 實驗參數設置

如表1所示，設置菌落總數生長模型實驗參數。

表1 參數設置Table 1 Parameter settings

4.2 結果與分析

圖1為混合霉菌動態生長模型中比生長速率μ和生長遲滯期λ的溫度-相對濕度響應曲面圖。由圖1A可知，μ隨溫度和相對濕度的增大呈先增大后減小的趨勢，在最適溫度和相對濕度時，μ最大；由圖1B可知，λ隨相對濕度增大而減小；隨溫度升高呈先減小后稍微增大的趨勢，在最適溫度附近時λ最小。這兩個菌落生長二級模型的響應曲面圖符合菌落生長特性，進一步證明了菌落動態生長模型的有效性。

圖1 比生長速率（A）與生長遲滯期（B）的溫度-相對濕度響應曲面圖Fig.1 Temperature-humidity response surface plots of specific growth rate (A) and growth lag (B)

根據2.2節中基于蒙特卡洛仿真的小麥倉儲環節菌落總數演化算法對其進行概率分布預測，結果見圖2。為更直觀地表現菌落總數的分布特征，通過蒙特卡洛仿真中大量重復實驗的粒子群統計并計算特定時刻的菌落總數概率分布函數（圖3）。

隨著時間的延長，菌落總數的生長演化規律呈現發散狀態，這是由于小麥倉儲環節中菌落的生長受到演化過程的隨機性和環境不確定性的影響。并且隨著時間的延長，菌落總數的演化曲線整體呈現先平坦后陡增的趨勢，這說明菌落生長是一個從慢變到爆發的過程，這符合菌落在生長遲滯期和指數期的生長特性（圖2）。菌落總數的概率函數隨時間的延長呈現由高窄向矮寬的變化過程，體現出菌落生長的隨機性隨時間延長而增加的特點，曲線越高則表示菌落總數出現概率越大（圖3）。結合圖2、3可知，菌落生長趨勢符合菌落生長規律并穩定在一定范圍內，其具有確定性；同時，小麥的倉儲環境雖然穩定，但是由于環境不確定性的影響，菌落總數概率分布預測具有一定的隨機性。

本實驗通過25 ℃、相對濕度0.6條件下‘鄭麥004’在倉儲環節中霉菌菌落總數的35 d實測數據演化曲線（數據間隔為7 d）驗證模型的有效性。實測數據演化曲線為圖2、3中的深藍色曲線，該曲線符合菌落生長規律，并基本被相同環境條件下的菌落總數蒙特卡洛仿真演化曲線覆蓋，這說明了該模型的有效性。各時刻菌落總數的實測數據接近蒙特卡洛仿真中相同時刻出現概率較高的菌落總數數值，這進一步驗證了模型的有效性。

圖2 基于蒙特卡洛仿真的菌落總數演化曲線Fig.2 Evolution curves of total bacterial count based on Monte Carlo simulation

圖3 菌落總數的概率分布函數Fig.3 Probability distribution function of total bacterial count

在不同時間菌落總數的概率分布結果的基礎上，計算不同時刻菌落總數的危害度，定量評價菌落總數風險的變化趨勢，結果見圖4。結果表明，危害度呈現先平坦后陡增、并逐漸趨近于1的趨勢。其符合菌落生長規律，能夠直觀明確地對菌落總數風險進行量化評價。

圖4 時區[15,30]下小麥倉儲環節中菌落總數危害度變化曲線Fig.4 Variation in the hazard degree of total bacterial count during wheat storage

5 結 論

本實驗考慮了蒙特卡洛仿真在微生物風險評估中的優勢，彌補了模型演化過程的隨機性及環境的不確定性對風險定量評估的不利影響，并結合度量指標“危害度”進行風險定量評價，降低風險評估的偶然性、提高其準確度和直觀性。未來研究主要分為兩個部分：1）本實驗研究的是單一環節中危害物的風險演化過程，而食品供應鏈是一個面向市場的組織網絡，針對其復雜的環節切換及環境變化，需進一步研究完整供應鏈上危害物的風險評估方法；2）基于蒙特卡洛仿真的風險評估方法需建立一定規模的粒子群，并對每個粒子進行完整的演化模擬實驗，這必將增大模型計算復雜度從而導致運算效率較低，針對這一缺點，需進一步從數學計算機理出發，通過改進數學解析算法降低風險評估方法的計算量，以提高其評估速度。