沈霞芬，蔡強*，吳利芳，肖革新

1.浙江清華長三角研究院（嘉興 314006）；2.國家食品安全風(fēng)險評估中心（北京 100022）

食品安全是個綜合概念，涉及衛(wèi)生、質(zhì)量、營養(yǎng)等多方面內(nèi)容，以及種植、養(yǎng)殖、生產(chǎn)、包裝、儲存、運輸、流通、消費等多個環(huán)節(jié)，單一的食品衛(wèi)生、食品質(zhì)量、食品營養(yǎng)等概念無法涵蓋以上全部內(nèi)容[1]。因此，關(guān)于食品安全，除了食品生產(chǎn)過程中所有指標(biāo)需要滿足特定的標(biāo)準(zhǔn)外，許多食品安全風(fēng)險的存在和發(fā)展由食品供應(yīng)鏈其他因素驅(qū)動，如空氣、土壤、水等農(nóng)業(yè)資源環(huán)境的影響，種養(yǎng)殖過程中化肥、農(nóng)獸藥等施用帶來的影響，人體對食物的攝入量等[2]。這些因素與供應(yīng)鏈之間的相互作用復(fù)雜，需要使用一種系統(tǒng)或整體的分析方法揭示其中的關(guān)系。目前，國內(nèi)在食品安全風(fēng)險領(lǐng)域的分析，主要以監(jiān)管部門監(jiān)督抽檢和日常檢查數(shù)據(jù)的風(fēng)險分析為主，食品加工流通環(huán)節(jié)外其他因素的分析尚未形成深入挖掘。黃湘鷺等[3]對2016—2017年國家食品安全監(jiān)督抽檢的不合格項目進行歸類分析，結(jié)果發(fā)現(xiàn)抽檢合格率呈現(xiàn)穩(wěn)中有升的趨勢。李宗亮[4]以食品的污染物和營養(yǎng)物質(zhì)為評價對象，以不合格率和不合格程度為評價指標(biāo)構(gòu)建食品風(fēng)險系數(shù)評價體系，并結(jié)合BP神經(jīng)網(wǎng)絡(luò)模型，將食品風(fēng)險系數(shù)應(yīng)用于各類食品及抽檢時待測參數(shù)的安全趨勢預(yù)測預(yù)警等領(lǐng)域，為監(jiān)管部門提供決策依據(jù)。相比之下，歐美發(fā)達國家的政府信息公開工作具有起步早、數(shù)據(jù)多、深度開放等特點，在食品安全分析方面也有著更進一步的研究[5]。Yamine等[6]使用貝葉斯網(wǎng)絡(luò)模型來識別和量化RASFF中關(guān)于源自印度、土耳其和荷蘭的水果和蔬菜的所有食品安全通告數(shù)據(jù)與產(chǎn)品原產(chǎn)國氣候因素、農(nóng)業(yè)因素和經(jīng)濟因素之間的關(guān)系，研究結(jié)果表明不能確定氣候因素對蔬果化學(xué)食品安全危害產(chǎn)生影響。Alberto等[7]根據(jù)聯(lián)合國糧食及農(nóng)業(yè)組織提供的關(guān)于進出口和生產(chǎn)的公開數(shù)據(jù)，構(gòu)建數(shù)學(xué)模型模擬糧食產(chǎn)品的分配情況，首次從食品安全風(fēng)險的角度論證網(wǎng)格結(jié)構(gòu)對食品分布影響。

從食品安全多個影響因素的綜合分析與關(guān)鍵風(fēng)險因子確定，以及單一指標(biāo)數(shù)據(jù)的預(yù)測2個角度探討食品安全風(fēng)險分析的數(shù)據(jù)模型，并給出相應(yīng)的應(yīng)用實例。

1 材料與方法

1.1 材料

1.1.1 數(shù)據(jù)來源

食品安全涉及從農(nóng)田到餐桌的整個過程，理論上需要收集覆蓋食品供應(yīng)鏈上的所有數(shù)據(jù)，同時融合時空、地理、環(huán)境、氣象及醫(yī)療健康等信息。但國內(nèi)目前尚未建立統(tǒng)一的數(shù)據(jù)標(biāo)準(zhǔn)，也尚未解決跨層級、跨部門、跨業(yè)務(wù)、全鏈條多源異構(gòu)數(shù)據(jù)的融合，所以數(shù)據(jù)收集十分困難。為既兼顧各指標(biāo)的合理性、時間上的統(tǒng)一又不失數(shù)據(jù)的準(zhǔn)確性，從國家統(tǒng)計局官網(wǎng)上收集1999—2018年關(guān)于中國環(huán)境、農(nóng)業(yè)因素及消費量的6個食品安全評價指標(biāo)數(shù)據(jù)，數(shù)據(jù)維度為6×20。同時，以“食物中毒事件”為關(guān)鍵詞，通過國家衛(wèi)生健康委員會官網(wǎng)（http：//www.nhc.gov.cn）檢索的“衛(wèi)生部辦公廳關(guān)于XXXX年全國食物中毒事件情況的通報”文件中，收集2008—2015年公開發(fā)布的中國食物中毒事件起數(shù)的月度數(shù)據(jù)，數(shù)據(jù)長度為96條。

1.1.2 分析軟件

采用R對食品安全6個影響因素進行主成分分析，并對食物中毒事件起數(shù)季度數(shù)據(jù)進行時間序列分析。以α=0.05為檢驗水準(zhǔn)，p＜0.05說明具有顯著性差異。

1.2 方法

1.2.1 多影響因素的綜合分析

采用主成分分析法從多個影響因素入手，提取主要成分，并基于成分的方差貢獻率計算綜合評分，以此作為各影響因素對年度食品安全的影響程度。采用灰色關(guān)聯(lián)度分析，計算食品安全影響因素與綜合評分之間的關(guān)聯(lián)度，確定最為關(guān)鍵的要素。

1) 構(gòu)建一個m×n的矩陣X，行表示年份，列表示食品安全的影響因素。

2) 采用Z-score對以上n個指標(biāo)進行標(biāo)準(zhǔn)化使其量綱統(tǒng)一，得到1組新變量（z1,z2,…,zn）[8]，并計算這組新變量的相關(guān)系數(shù)矩陣R。

3) 由相關(guān)系數(shù)矩陣R計算得到其特征值λ，并從大到小排序，λ1≥λ2≥…≥λn≥0，且計算各特征值對應(yīng)的特征向量εi=（εi1,εi2,…,εin），i=1,2,…,n，確定X的第i個主成分。

4) 采用灰色關(guān)聯(lián)度計算各評價指標(biāo)Xi與綜合評分Y之間關(guān)聯(lián)度[9]，關(guān)聯(lián)系數(shù)定義為：

式中：i為第i個評價指標(biāo)；k為各指標(biāo)的第k個數(shù)值，故Y(k)為綜合評分Y的第k個數(shù)值，Xi(k)為第i個評價指標(biāo)的第k個數(shù)值；ρ為分辨系數(shù)，取值為（0,1），ρ越小，區(qū)分度越低，一般取值為0.5。

由于關(guān)聯(lián)系數(shù)是比較數(shù)列與參考數(shù)列在各個時刻的關(guān)聯(lián)程度值，所以它是一個序列，信息過于分散不便于進行整體比較。因此需要將各時刻的關(guān)聯(lián)系數(shù)求其平均值，作為比較數(shù)列xi與參考數(shù)列Y之間關(guān)聯(lián)度的數(shù)量表示，關(guān)聯(lián)度ri定義為：

ri越大，表明指標(biāo)i對綜合評分的影響越大。

1.2.2 單一指標(biāo)數(shù)據(jù)的預(yù)測分析

對單位時間內(nèi)某個指標(biāo)數(shù)據(jù)進行統(tǒng)計，如在農(nóng)業(yè)生產(chǎn)中，農(nóng)獸藥和化肥的使用量可按季度、按月記錄；加工運輸環(huán)節(jié)中，食品生產(chǎn)件數(shù)、檢驗檢測信息等可按周、按天記錄，這些數(shù)據(jù)便為時間序列。對于數(shù)據(jù)都為正值的時間序列，可通過5個步驟完成模型構(gòu)建及分析。

1) 對于非線性的變量進行對數(shù)變換，消除數(shù)據(jù)的非線性，且不改變數(shù)據(jù)之間的統(tǒng)計性質(zhì)。

2) 對于非平穩(wěn)的變量進行差分變換，弱化隨機性使其平穩(wěn)化。利用ADF單位根檢驗，檢驗變換后的序列是否平穩(wěn)。具體方法：原假設(shè)H0，序列非平穩(wěn)；備擇假設(shè)H1，序列平穩(wěn)。檢驗統(tǒng)計量小于5%臨界值時拒絕原假設(shè)，即認為該序列平穩(wěn)。

3) 對d階差分后的平穩(wěn)時間序列構(gòu)建ARIMA模型，若模型殘差的Ljung-Box統(tǒng)計量的p值大于0.05，則模型通過檢驗[10]，再利用所構(gòu)建的模型進行預(yù)測。

2 結(jié)果與分析

2.1 食品安全綜合評價及關(guān)鍵影響因素確定

2.1.1 食品安全風(fēng)險的影響因素

選取的食品安全風(fēng)險指標(biāo)包含1999—2018年環(huán)境、農(nóng)業(yè)因素和消費量等數(shù)據(jù)。食品安全風(fēng)險與各影響因素值的大小呈正比。

2.1.2 食品安全風(fēng)險的綜合評分

從標(biāo)準(zhǔn)化后各影響因素的散點圖矩陣來看，除了廢氣排放量與其他幾個指標(biāo)的相關(guān)性較弱之外（各點越接近矩形的對角線表明相關(guān)性越強），其余指標(biāo)之間都具有一定的相關(guān)性，存在一定的信息重疊（圖1）。主成分分析選取的前2個主成分的特征值分別為3.34和2.27，前2個主成分的累計貢獻率達93%，基本包括原來那些變量所包含信息。

圖1 食品安全影響因素的散點圖矩陣

由2個主成分加權(quán)所得的綜合得分，代表著環(huán)境、農(nóng)業(yè)和消費因素對食品安全的影響程度（表2）。

表2 1999—2018年食品安全影響因素的綜合評分

結(jié)果表明，環(huán)境、農(nóng)業(yè)及消費量對中國食品安全風(fēng)險的影響具有階段性：1999—2002年中國食品安全風(fēng)險逐年降低；2002—2012年食品安全風(fēng)險基本上呈現(xiàn)逐年升高趨勢；直到2012年之后又基本上恢復(fù)到逐年降低趨勢。產(chǎn)生該結(jié)果的主要原因在于工業(yè)與經(jīng)濟的發(fā)展加速了環(huán)境污染及農(nóng)業(yè)生產(chǎn)過程中農(nóng)獸藥與化肥的使用量，導(dǎo)致2002—2012年食品安全的風(fēng)險逐漸增高，但以犧牲環(huán)境與國民健康的工業(yè)和經(jīng)濟發(fā)展打破中國可持續(xù)發(fā)展的基本理念?！笆舜蟆敝笳_始出臺各種環(huán)保與民生健康的政策[11-12]，使得工業(yè)廢棄物及農(nóng)業(yè)生產(chǎn)過程中化學(xué)品的投入量有所控制，因此食品安全風(fēng)險在2012年后呈現(xiàn)明顯降低趨勢。

表1 食品安全風(fēng)險因素

2.1.3 食品安全關(guān)鍵影響要素確定

灰色關(guān)聯(lián)度分析得到各影響因素與食品安全風(fēng)險綜合評分之間的關(guān)聯(lián)度分別為0.637，0.663，0.682，0.741，0.714和0.598。因此1999—2018年影響中國食品安全的因素排名依次為單位面積農(nóng)獸藥使用量、單位面積農(nóng)用化肥施用折純量、廢氣（二氧化硫+煙粉塵）排放量、廢水排放量、工業(yè)固體廢物傾倒丟棄量、居民主要食品人均消費量。該結(jié)果表明，單從環(huán)境、農(nóng)業(yè)和消費層面來看，農(nóng)業(yè)生產(chǎn)過程中投入的化學(xué)品量對于中國食品安全來說是最為顯著的影響因素，其次為環(huán)境因素。為減輕病蟲草害對產(chǎn)量的威脅并獲得高產(chǎn)，農(nóng)民對于農(nóng)獸藥和化肥不合理使用的行為比較普遍。陳曉明等[13]表明1991—2013年中國農(nóng)藥使用量增長135.5%，不合理的使用行為導(dǎo)致其殘留物在食物鏈中不斷傳遞和遷移[13]。史常亮等[14]指出中國每公頃耕地化肥施用量由1980年的127.8 kg增加到2013年的437.4 kg，是國際公認化肥施用安全上限225 kg/hm2的1.94倍。過度使用化肥產(chǎn)生的污染土壤的重金屬和放射性物質(zhì)，通過食物鏈不斷在植物和生物體內(nèi)富集，最終在人體內(nèi)積累危害人體健康。

2.2 時間序列分析預(yù)測食品安全風(fēng)險

2.2.1 季節(jié)性分布結(jié)果

2008—2015年全國食物中毒事件起數(shù)總體上呈現(xiàn)先降低后增高趨勢（圖2），并且存在明顯季度模式。第3季度的食物中毒事件呈現(xiàn)年度最高，第2季度次之，第1季度最低。因此，以季度作為時間單位，對食物中毒事件做時間序列預(yù)測將更具準(zhǔn)確性。

圖2 2008—2015年全國食物中毒事件起數(shù)的月度數(shù)據(jù)時序圖

2.2.2 平穩(wěn)性檢驗結(jié)果

消除非線性后的季度食物中毒事件序列，經(jīng)過1階差分后還是存在季節(jié)與非季節(jié)相關(guān)性（圖3a2）；經(jīng)過1階季節(jié)差分后的序列消除相關(guān)性（圖3a3）；經(jīng)過1階差分及1階季節(jié)差分后的序列在滯后1階和5階上還是存在一定程度相關(guān)性（圖3a4）。

圖3 食物中毒事件對數(shù)序列及其差分后序列的自相關(guān)（ACF，a）和偏自相關(guān)（PACF，b）函數(shù)圖

平穩(wěn)性檢驗結(jié)果表明，除了食物中毒事件起數(shù)的對數(shù)序列外，其他經(jīng)過差分處理后序列都平穩(wěn)（p＜0.05），符合ARIMA模型建立的條件（表3）。

表3 所有序列的平穩(wěn)性檢驗結(jié)果

2.2.3 ARIMA模型及預(yù)測結(jié)果

根據(jù)差分后序列的自相關(guān)（ACF）和偏自相關(guān)（PACF）函數(shù)圖（圖3）建立3個ARIMA模型（表4），模型殘差的Ljung-Box統(tǒng)計量的p值及殘差檢驗圖均表明3個模型通過檢驗。

表4 根據(jù)自相關(guān)（ACF）和偏自相關(guān)（PACF）圖確定ARIMA模型的擬合優(yōu)度指標(biāo)

表5為由R自動選取的ARIMA模型。其中，ARIMA (3,1,0)等同于表4中第一個模型，但由于調(diào)用的是R中的不同函數(shù)，所以在結(jié)果上略顯差異；考慮季節(jié)性的模型，R自動選擇的是ARIMA(0,1,1)×(0,1,0)4，雖然模型殘差的Ljung-Box統(tǒng)計量的p值大于0.05，但相比于其他幾個模型的p值小很多，且殘差檢驗圖顯示，殘差在滯后5階及其后面的階數(shù)上存在較小相關(guān)性。

表5 R自動選擇ARIMA模型的擬合優(yōu)度指標(biāo)

以AIC最小原則，選擇最好解釋數(shù)據(jù)但含最少自由參數(shù)的模型ARIMA(5,0,0)×(0,1,0)4，即：

作為全國食物中毒事件起數(shù)的時間序列模型。模型的擬合值與真實值的變動具有較高的一致性（圖4a）。該模型對2015年度中毒事件起數(shù)的預(yù)測顯示出較強的季節(jié)性，預(yù)測2015年第3季度的中毒事件數(shù)仍然是全年最高，達到66件；第1季度和第4季度的中毒事件起數(shù)相對較低，分別為16和27件；第2季度的中毒事件起數(shù)全年第二，為35件。產(chǎn)生該結(jié)果的主要原因是春夏秋季節(jié)天氣溫?zé)?、氣候濕潤，適宜細菌和有毒植物的生長繁殖，成為食物中毒事件的高發(fā)期。從預(yù)測效果上來看，該模型對2015年第2季度的預(yù)測完全吻合，除了第1季度稍有偏差外，其余3個季度預(yù)測效果較好，真實值都位于95%的預(yù)測區(qū)間內(nèi)[7]（圖4b）。

圖4 模型ARIMA(5,0,0)×(0,1,0)4的擬合結(jié)果（a）與預(yù)測結(jié)果（b）

3 結(jié)論

食品安全風(fēng)險分析是一個針對復(fù)雜系統(tǒng)的綜合性評價問題。主成分分析與灰色關(guān)聯(lián)度結(jié)合能夠客觀、全面地反映各項影響因素的信息，同時通過估計綜合評分與各影響因素之間的一致性程度判斷出最為關(guān)鍵的影響因素。在食品領(lǐng)域，由于氣候變化、飲食習(xí)慣等因素的影響，很多食品安全相關(guān)的數(shù)據(jù)具有周期性。時間序列ARIMA模型能夠根據(jù)某個指標(biāo)的歷史數(shù)據(jù)，通過計量的手段反映時序中動態(tài)依存關(guān)系，不僅從數(shù)量上揭示該指標(biāo)的發(fā)展變化規(guī)律，而且能夠準(zhǔn)確預(yù)測未來走勢。

然而，在食品領(lǐng)域，食品安全的影響因素不僅是環(huán)境、農(nóng)業(yè)和消費因素，還關(guān)系到食品生產(chǎn)和流通過程中食品質(zhì)量等問題，理論上需要實現(xiàn)對食品從農(nóng)田到餐桌的全面評價。但由于目前國內(nèi)食品鏈上各環(huán)節(jié)信息化水平不一，一些環(huán)節(jié)數(shù)據(jù)采集缺失，數(shù)據(jù)存儲沒有統(tǒng)一的標(biāo)準(zhǔn)，使得食品安全的綜合評價缺乏全面性。但從分析方法上來說，主成分分析與時間序列分析對于食品安全領(lǐng)域多因素的綜合評價和單一指標(biāo)的預(yù)測分析都具有較好的指示作用。