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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 蒙特卡洛仿真

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 案例驗(yàn)證

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

4.1 實(shí)驗(yàn)參數(shù)設(shè)置

如表1所示，設(shè)置菌落總數(shù)生長(zhǎng)模型實(shí)驗(yàn)參數(shù)。

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

4.2 結(jié)果與分析

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

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

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

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

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

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

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

5 結(jié) 論

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