基于蒙特卡洛仿真的小麥倉儲環(huán)節(jié)菌落總數(shù)風(fēng)險評估方法

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

（1.北京工商大學(xué)計算機與信息工程學(xué)院，北京 100048；2.北京工商大學(xué)，中國輕工業(yè)工業(yè)互聯(lián)網(wǎng)與大數(shù)據(jù)重點實驗室，北京 100048；3.北京工商大學(xué)，北京市食品添加劑工程技術(shù)研究中心，北京 100048；4.北京市農(nóng)業(yè)信息技術(shù)研究中心，北京 100097）

小麥?zhǔn)鞘澜绲诙蠹Z食作物，也是中國繼玉米和大米之后的第三大作物。20世紀(jì)80年代初以來，全世界的小麥產(chǎn)量顯著增加，中國每年的小麥產(chǎn)量超過1.2億 t，已成為世界上最大的小麥生產(chǎn)國。然而，近年來小麥?zhǔn)称钒踩珕栴}層出不窮并引起了社會各界的廣泛關(guān)注[1-4]。風(fēng)險評估是一種結(jié)構(gòu)化的科學(xué)過程，用于估計風(fēng)險的概率和嚴(yán)重程度以及隨之而來的不確定性，以確定因接觸食物中的生物、化學(xué)或物理危害物而對健康生活造成的潛在危害和相關(guān)風(fēng)險[5-6]。國際食品安全組織正大力推動將風(fēng)險評估技術(shù)應(yīng)用于評估食品中危害物安全的問題[7]。目前的中國食品安全法要求建立國家食品安全風(fēng)險評估制度，以評估中國食品中的危害物風(fēng)險[8]。

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

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

1 菌落動態(tài)生長模型構(gòu)建

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

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

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

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

結(jié)合公式（1）、（2），建立菌落動態(tài)生長模型，如公式（3）所示。

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

定義菌落動態(tài)生長模型的菌落生長指標(biāo)x根據(jù)公式（4）計算。

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

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

需要注意的是，不同菌落生長特性不同，同種菌落的生長環(huán)境也會有差別，因此，不同菌落動態(tài)生長模型參數(shù)也不同。為使模型更為精確、合理，模型參數(shù)θ的取值需結(jié)合實測數(shù)據(jù)和相關(guān)參數(shù)估計方法確定。

2 小麥倉儲環(huán)節(jié)中菌落總數(shù)的概率分布預(yù)測

2.1 蒙特卡洛仿真

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

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

2.2 基于蒙特卡洛仿真的小麥倉儲環(huán)節(jié)中菌落總數(shù)的概率分布預(yù)測

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

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

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

基于以上設(shè)置，可利用蒙特卡洛仿真對菌落總數(shù)演化過程進行模擬。針對蒙特卡洛仿真中的每一次模擬，同一指標(biāo)在同一預(yù)測時間點的數(shù)值均不相同，這些數(shù)值同時具有規(guī)律性和隨機性。規(guī)律性體現(xiàn)在預(yù)測值是基于確定的菌落動態(tài)生長模型、確定的菌落生長指標(biāo)和模型參數(shù)初值產(chǎn)生的；隨機性體現(xiàn)在每一個菌落生長指標(biāo)在每一個預(yù)測時間點均受噪聲影響，導(dǎo)致得到的預(yù)測值不盡相同。當(dāng)模擬次數(shù)足夠大時，結(jié)合概率統(tǒng)計規(guī)律，可以獲得各菌落生長指標(biāo)在不同時間點上取值的概率分布?；诿商乜宸抡娴男←渹}儲環(huán)節(jié)中菌落總數(shù)概率分布預(yù)測的具體步驟為：

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

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

3）對于時刻k∈[0,K]，基于k－h(huán)時刻的菌落生長指標(biāo)和模型參數(shù)，根據(jù)公式（6）進行單步預(yù)測，分別得到；

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

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

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

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

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

3 基于“危害度”的小麥倉儲環(huán)節(jié)中菌落總數(shù)風(fēng)險評估

在系統(tǒng)工程領(lǐng)域，相關(guān)學(xué)者基于可靠性理論中性能可靠度的定義和安全性分析中安全概率的定義[32-33]，提出了“健康度”的概念，并以此作為評估系統(tǒng)工作性能及表現(xiàn)的度量指標(biāo)[34-35]。本實驗參考“健康度”的思想提出了新的度量指標(biāo)“危害度”，并基于該指標(biāo)對小麥倉儲環(huán)節(jié)中菌落總數(shù)風(fēng)險進行定量評價。

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

式中：Rk為k時刻該動態(tài)系統(tǒng)處于有風(fēng)險狀態(tài)的概率，也稱危害度；對于小麥倉儲環(huán)節(jié)的菌落總數(shù)而言，小麥倉儲環(huán)節(jié)中菌落生長指標(biāo)x的演化過程可以看成一個動態(tài)系統(tǒng)。

結(jié)合公式（10），小麥倉儲環(huán)節(jié)中菌落總數(shù)“危害度”可通過公式（12）計算。

式中：SF＝{xk|Nk＞Θ,Nk∈xk}，Θ為菌落總數(shù)超標(biāo)閾值/（CFU/g）；Rc,k∈[0,1]，Rc,k＝0表示k時刻小麥倉儲環(huán)節(jié)中菌落總數(shù)無超標(biāo)風(fēng)險；Rc,k＝1表示k時刻小麥倉儲環(huán)節(jié)中菌落總數(shù)存在超標(biāo)風(fēng)險。

4 案例驗證

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

4.1 實驗參數(shù)設(shè)置

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

表1 參數(shù)設(shè)置Table 1 Parameter settings

4.2 結(jié)果與分析

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

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

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

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

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

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

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

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

圖4 時區(qū)[15,30]下小麥倉儲環(huán)節(jié)中菌落總數(shù)危害度變化曲線Fig.4 Variation in the hazard degree of total bacterial count during wheat storage

5 結(jié) 論

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