尹 佳，黃 茜，陳 翔，陳 晨，陳 鋰，張 濤，徐 成，黃亞平，郭鵬程，文 紅,*

（1.湖北省食品質量安全監(jiān)督檢驗研究院，國家市場監(jiān)管重點實驗室（動物源性食品中重點化學危害物檢測技術），湖北省食品質量安全檢測工程技術研究中心，湖北 武漢 430075；2.武漢理工大學計算機與人工智能學院，湖北 武漢 430070；3.武漢理工大學理學院，湖北 武漢 430070）

近30多年來，中國已經成為世界上食物生產量和消費需求量增長最快的國家之一[1]?！吨袊澄锱c營養(yǎng)發(fā)展綱要（2014—2020）》[2]指出我國近年來在食物總量上已經實現(xiàn)了“食物供需基本平衡”，即已基本解決糧食數(shù)量安全問題。但目前在有效滿足人們飲食需要的同時，部分企業(yè)還存在重產量而忽視質量的現(xiàn)象，暴露出一些食品安全問題。因此，如何在能夠“保障食物有效供給”的前提下，保障食品安全和營養(yǎng)健康、控制有毒有害物質殘留超標、快速對風險發(fā)出預警，便成為下一步急需解決的問題，同時，食品安全預測研究也具有重要的社會意義。

實踐表明，通過對還未發(fā)生的潛在風險進行分析，提前判斷，主動預防，被認為是保證食品安全的最有效方式[3]。目前，對于食品質量安全風險預警的研究，統(tǒng)計方法和機器學習預測相結合的方法已占主導地位[4]，樓皓等[5]使用差分自回歸移動平均（autoregressive integrated moving average，ARIMA）-支持向量機（support vector machine，SVM）組合模型對中國出口歐盟的食品風險進行預測，結果表明組合模型比單一模型具有更高的精度；Yan Shengyang等[6]提出了一種基于反向傳播（back propagation，BP）神經網絡和遺傳算法的食品安全綜合指數(shù)預測方法；白寶光等[7]構建了預警指標體系和BP神經網絡預警模型，提高了乳制品風險預警的精準性；Yonar等[8]利用Holt線性趨勢模型對多個國家的小麥產量進行了預測；Sivamani等[9]使用季節(jié)性ARIMA模型為畜禽類食品供應提供準確的預測。上述主流方法主要是點估計，而區(qū)間估計在食品安全風險預警領域應用較少，但在其他領域已有廣泛的應用。趙會茹等[10]利用核密度估計（kernel density estimation，KDE）方法對長短期記憶神經網絡預測模型預測結果進行區(qū)間估計，建立了未來短期電力負荷的區(qū)間預測方法；張明宇等[11]基于改進的小波神經網絡，對上證指數(shù)開盤數(shù)據建立了一種新型股指區(qū)間預測模型；丁藤等[12]對風速建立ARIMA-廣義自回歸條件異方差（generalized autoregressive conditional heteroskedast，GARCH）模型，實現(xiàn)了通過少量歷史數(shù)據對超短期內風速的區(qū)間預測；Guo Tianli等[13]建立了自激發(fā)門限自回歸（self-exciting threshold autoregressive，SETAR）-GARCH混合模型，對非平穩(wěn)非線性地下水深度進行了評價。以上關于食品安全預警方面的模型主要是點估計模型，此類模型雖然有一定的預測精度，但其通常只能適用于線性變化的趨勢，且只能提供確定性的預測結果，而不能提供有關不確定性的信息[14]。而區(qū)間估計是在點估計的基礎上，給出總體參數(shù)估計的一個區(qū)間范圍，不僅可以提供更多的信息，而且可以量化預測結果的不確定性[15-18]，但其在食品安全預測方面應用鮮有報道。

因此，本實驗提出了一種點估計和區(qū)間估計組合預測模型——小波包分解（wavelet packet decomposition，WPD）-ARIMA-GARCH模型，應用于醬鹵肉制品安全風險預警的區(qū)間預測。在點估計部分提出了一種基于數(shù)據分解WPD[19]的ARIMA方法，對于缺失數(shù)據采用高斯過程回歸（Gaussian process regression，GPR）進行插值，采用WPD進行分解，依據ARIMA模型自動定階，對分解后的分量進行點估計預測，將各預測分量重構后輸出點估計部分預測結果；在區(qū)間估計部分，選擇GARCH模型[20]對點估計結果殘差進行區(qū)間估計，提供每一個確定預測結果的上限和下限，以期通過建立的組合模型量化醬鹵肉制品質量安全風險預測結果的不確定性，為其風險防控提供一定的參考。最后，將建立的最優(yōu)模型與其他預測模型如支持向量回歸（support vector regression，SVR）和KDE等進行了對比，驗證所建立組合模型結果的有效性。

1 材料與方法

本實驗選取2014—2019年來自國家市場監(jiān)督管理總局公布的以及檢測機構內部自行檢測獲得的部分醬鹵肉制品數(shù)據作為數(shù)據源，其中，前5 年作為訓練集，2019年數(shù)據作為測試集。為了更加全面了解醬鹵肉制品中存在的風險，將所有檢測項目均納入分析。由于存在檢驗項目結果信息不完全、數(shù)據的結構不統(tǒng)一以及稠密性差異等問題[21-22]，在建立模型之前，有必要對原始檢測數(shù)據進行清洗、集成、變換等預處理。

1.1 數(shù)據預處理

將檢測項目的檢測值（微生物項目除外）結合國家標準（例如GB 2760—2014《食品安全國家標準 食品添加劑使用標準》等）采用公式（1）進行去量綱化處理。根據去量綱化的結果結合專家打分法將項目風險等級劃分為5 級，1級為安全無風險，2級為輕微風險，3級為輕度風險，4級為中度風險，5級為重度風險（即不符合國家標準要求）。其中較低的檢測值（1～2級）視為相對安全的情況，3級需要引起關注，而4～5級則需要進行早期預警，并采取相應的預防措施。具體劃分標準見表1。

表1 檢驗項目的風險等級劃分標準Table 1 Criteria for risk classification of inspection items

式中：Yi為預處理后的風險等級評價值；Xstandard為國家標準中規(guī)定的標準值；Xi為檢驗項目的檢測值。

微生物項目參考表1按照檢測值進行風險等級劃分。

本實驗采用改進的softmax函數(shù)計算單一產品的綜合風險等級（公式（2）），通過函數(shù)中指數(shù)權重的變化來調節(jié)風險等級的權重[23]，將自然周作為一個數(shù)據集，通過加權求和（公式（3））計算自然周的綜合風險指數(shù)。

式中：level(A)為單一產品A的綜合風險等級；i為食品A中檢測項目的風險等級數(shù)值；ωi為風險等級i在食品A中的占比。

式中：weekrisk為自然周的綜合風險指數(shù)；i為該時間段的風險等級數(shù)值；ωi為該時間段風險等級i的占比。

通過對原始周綜合風險時間序列的數(shù)據特征進行分析，結合專家打分法，將自然周綜合風險指數(shù)weekrisk≤8劃分為低風險，8～21劃分為中風險，＞21劃分為高風險。其中中風險需要給予一定關注，高風險需要引起重視，并采取早期預警和預防措施。

1.2 缺失數(shù)據插值

由于基于時間序列的預測模型需要以數(shù)據的連續(xù)性和相對完整性為基礎，而本實驗通過自然周分箱后的檢測數(shù)據存在部分缺失情況，故帶入模型訓練前需要進行插值。

1.2.1 數(shù)據缺失情況分析

以某地區(qū)的數(shù)據作為插值實驗的數(shù)據集，該地區(qū)周綜合風險的數(shù)據情況如圖1所示。

圖1 某地區(qū)醬鹵肉的周綜合風險Fig.1 Weekly risk of soy sauce and pot-roast meat products from a certain region

分析該地區(qū)整體數(shù)據分布情況，數(shù)據點分為兩種類型：①平穩(wěn)較多，偶爾波動；②波動較多，平穩(wěn)較少。本實驗選擇該地區(qū)270～300 周的數(shù)據集A和310～350 周的數(shù)據集B作為這兩類情況的代表數(shù)據集。同時，分析缺失點數(shù)據，存在兩種情況：①連續(xù)型缺失，即中間連續(xù)一段時間都沒有數(shù)據；②間斷型缺失，即中間偶爾會出現(xiàn)缺失點。故該地區(qū)的缺失數(shù)據有以下4 種情況：

·D1：A集上的連續(xù)缺失；

·D2：A集上的間隔缺失；

·D3：B集上的連續(xù)缺失；

·D4：B集上的間隔缺失。

1.2.2 插值方法

針對數(shù)據缺失問題，本實驗采用GPR進行插值。

GPR使用多維高斯過程條件分布的性質來進行回歸預測，其具體的表達如式（4）所示：

式中：Xi為特征向量；yi為對應的值。

對于一個新的觀測值X*，GPR給出的預測值y*如式（5）所示：

式中：K為協(xié)方差核矩陣，例如徑向基核函數(shù)、線性核函數(shù)等。根據高斯函數(shù)的性質，當y*＝k(X*,X)K(X,X)-1y，即預測值的平均值時，P(y*|y)達到峰值。

1.2.3 插值方法對比

線性插值、二次樣條插值是常用的插值方法，在本實驗中用作插值對比方法。為了合理地評價各個插值方法的插值效果，假設原始序列中某些數(shù)據是缺失，通過計算各插值方法給出的預測值和真實值之間的誤差來評價插值方法的優(yōu)劣。在選擇數(shù)據樣本時，制定了以下幾條準則：

a）數(shù)據必須是完整無缺失的；

b）數(shù)據量必須保證T≥L；其中，T表示樣本數(shù)據段的數(shù)據量，L表示該樣本假設缺失數(shù)據的數(shù)據量；

c）數(shù)據樣本應當具有代表性。

本實驗將缺失個數(shù)L設置為10，隨機設置缺失值，在每個數(shù)據集Di上重復5 次實驗，根據式（6）計算各個插值方法的均方誤差（mean-square error，MSE）和加權平均誤差mseweight。

式中：ωi為Di的權重。考慮到D3的情況較少，故將其權重設為0.1，其他3 種權重均設為0.3。

1.3 WPD-ARIMA-GARCH組合模型的建立

1.3.1 WPD

小波分解是一種信號時頻分析方法，在小波分析理論[24]的基礎之上，WPD引入了最優(yōu)基選擇的概念，將頻帶經過多層次的劃分之后，根據被分析信號的特征，自適應地選取最優(yōu)基函數(shù)，使之與信號相匹配，以提高信號的分析能力[25-26]。WPD既對低頻部分進行分解，也對高頻部分進行分解，可以對原始信號進行更加細致的劃分，分解示意圖如圖2所示。

圖2 三級WPD示意圖Fig.2 Schematic diagram of three-scale WPD

1.3.2 ARIMA模型

ARIMA模型是Box和Jenkins在20世紀70年代初提出的一種著名的時間序列預測方法，并改進衍生出諸多精度優(yōu)良的模型[27]。在ARIMA模型中，最常用的是求和自回歸移動平均模型ARIMA(p,d,q)[28]，其中，p為自回歸階數(shù)，d為成為平穩(wěn)時間序列時所做的差分次數(shù)，q為移動平均階數(shù)。其基本原理是將非平穩(wěn)時間序列經過差分運算轉化為平穩(wěn)時間序列，再將因變量僅對它的滯后值以及隨機誤差項的現(xiàn)值和滯后值進行回歸所建立的模型[29]。

1.3.3 GARCH模型

GARCH模型由Bollerslev在1986年提出[30]，相較于自回歸條件異方差（autoregressive conditional heteroskedasticity，ARCH）[31]，其允許自決定條件方差，可以避免必要的高滯后順序的困境，在波動性分析和預測方面應用廣泛[32-34]。同時，GARCH模型也是一種概率模型，本實驗將其用于置信區(qū)間的區(qū)間估計，預測區(qū)間表示為式（7）：

式中：μt為時間t的預測平均值；Z為Z分數(shù)；α為置信水平；σt為時間t的標準偏差；n為樣本量。

1.3.4 WPD-ARIMA-GARCH組合模型

本實驗提出了一種WPD-ARIMA-GARCH組合模型，采用點估計和區(qū)間估計結合的方式對醬鹵肉制品的質量安全風險進行預測，具體流程如圖3所示。在點估計階段，采用WPD-ARIMA模型，對于缺失數(shù)據采用GPR進行插值，采用WPD進行分解，依據ARIMA模型自動定階，對分解后的分量進行點估計預測，將各預測分量重構后輸出點估計部分預測結果；在區(qū)間估計階段，采用GARCH模型對點估計結果殘差進行區(qū)間估計。

圖3 WPD-ARIMA-GARCH組合模型流程圖Fig.3 Flow chart of WPD-ARIMA-GARCH model

ARIMA-GARCH模型的運行過程如表2所示。

表2 ARIMA-GARCH模型的運行過程Table 2 Operation process of ARIMA-GARCH model

1.4 預測模型的評價指標

1.4.1 點估計模型評價指標

本實驗選用以下3 個指標——MSE、平均絕對誤差（mean absolute error，MAE）和平均絕對百分比誤差（mean absolute percentage error，MAPE）作為點估計模型的評價指標，具體計算如式（8）～（10）所示。

式中：n為樣本數(shù)量；yi和分別為第i個實際檢測值和預測值。

1.4.2 區(qū)間估計模型評價指標

本實驗選用以下3 個指標——預測區(qū)間覆蓋率（prediction interval coverage probability，PICP）、預測區(qū)間平均寬度（prediction interval normalized average width，PINAW）、覆蓋寬度標準（coverage width-based criterion，CWC）作為區(qū)間估計模型的評價指標。具體計算如式（11）～（14）所示[35]。

式中：n為樣本數(shù)量；當?shù)趇個真實值落在[Li,Ui]中時ci＝1，否則ci＝0（[Li,Ui]表示第i個預測區(qū)間，Li和Ui分別是該區(qū)間的下界和上界）。

蚤狀幼體Ⅰ期已有口器，開始攝食，主要投喂蛋黃、酵母。每隔3小時投喂一次，蚤狀幼體后期可投喂輪蟲及少量的鹵蟲幼體。水溫控制在22～24℃，充氣量微充氣略顯沸騰狀。每天換水30～40cm，換水網箱網目為80目。

式中：R為目標值的取值范圍。

式中：η和μ為懲罰系數(shù)，本實驗中η的取值為0.7。

1.5 對比模型

1.5.1 點估計對比模型

本實驗構建了SVR、WPD-SVR、ARIMA模型作為點估計預測階段對比模型。SVR是一種常用的點估計模型，其核心是誤差函數(shù)和核函數(shù)，通過最大化間隔帶的寬度與最小化總損失優(yōu)化模型[36]。在點估計模型中，設置lookback＝6，即使用6 個連續(xù)的數(shù)據預測第7個數(shù)據。

1.5.2 區(qū)間估計對比模型

本實驗分別采用KDE和GPR模型作為區(qū)間估計對比模型。KDE屬于非參數(shù)檢驗方法之一，是在單變量KDE的基礎上，通過對KDE變異系數(shù)的加權處理，可以建立不同風險值的預測模型[37]。KDE不利用有關數(shù)據分布的先驗知識，對數(shù)據分布不附加任何假定，是一種從數(shù)據樣本本身出發(fā)研究數(shù)據分布特征的方法。GPR是一個貝葉斯方法的概率模型，不僅可以作為缺失數(shù)據插值方法之一，還可以獲得整個回歸函數(shù)的分布，計算出其預測區(qū)間。

2 結果與分析

2.1 插值方法結果比較

在缺失數(shù)據插值部分，對于每個數(shù)據集Di，GPR、線性插值和二次樣條插值的均方誤差MSE(Di)和加權平均誤差mseweight結果見表3。

表3 不同的插值方法結果比較Table 3 Comparison of interpolation results using different methods

插值后該地區(qū)的周風險時間序列見圖4，其中2014—2018年的檢測數(shù)據作為訓練集，2019年檢測數(shù)據作為測試集。

圖4 某地區(qū)醬鹵肉的周綜合風險時間序列Fig.4 Weekly risk time series of soy sauce and pot-roast meat products from a certain region

2.2 點估計

2.2.1 WPD-ARIMA模型

將采用GPR插值后的周風險數(shù)據作為原始序列，輸入建立的WPD-ARIMA模型，對其進行WPD。

WPD用光滑性較好的8 階Daubechies小波基，對原始序列進行3 層分解，將原始信號分解為[‘AAA’‘AAD’‘ADD’‘ADA’‘DDA’‘DDD’‘DAD’‘DAA’]8 個子序列S1～S8。圖5為某地區(qū)分解后各級分量示意圖。

圖5 某地區(qū)的WPD示意圖Fig.5 WPD results of soy sauce and pot-roast meat products from a certain region

由于S1序列數(shù)據的波動范圍較大，通過增強迪基-福勒（Augmented Dickey-Fuller，ADF）檢驗，該序列為非平穩(wěn)序列，通過差分處理后進行建模。其余子序列P值均小于檢驗水平的臨界值（α＝0.05），滿足平穩(wěn)性要求，可以直接進行建模。

對差分處理后的子序列，根據赤池信息準則（Akaike information criterion，AIC）和貝葉斯信息準則（Bayesian information criterion，BIC）值越低越好的原則，依據ARIMA模型的自動定階確定最優(yōu)參數(shù)后（表4），對WPD得到的各個分量進行預測，將各預測分量重構后輸出最終的預測結果，并使用測試集用來驗證該模型的精確度。

表4 某地區(qū)醬鹵肉WPD各分量的ARIMA模型最優(yōu)參數(shù)Table 4 Optimal parameters for the ARIMA model of WPD components for soy sauce and pot-roast meat products from a certain region

2.2.2 點估計模型比較

為評價所構建的WPD-ARIMA模型在醬鹵肉制品預測中的效果，本研究分別采用SVR、ARIMA、WPDSVR模型進行對比分析。圖6和表5為某地區(qū)4 種點估計方法的預測結果。

圖6 某地區(qū)醬鹵肉的點估計結果Fig.6 Point estimation results of soy sauce and pot-roast meat products from a certain region

表5 基于各種點估計方法的結果比較Table 5 Comparison of results based on various point estimation methods

通過上述結果分析發(fā)現(xiàn)，未經WPD分解的ARIMA模型和SVR模型擬合度明顯較差，誤差也明顯增大，原因可能是WPD能夠自適應地選擇基函數(shù)，對原始數(shù)據進行更細致地劃分，使分解后的各個分量有更好的光滑性，從而使得組合模型的誤差更小，擬合度更高；WPD-SVR模型的誤差略高于WPD-ARIMA模型，可能是因為SVR對參數(shù)和核函數(shù)選擇敏感，性能的優(yōu)劣主要取決于核函數(shù)的選取。因此，本研究選擇WPD-ARIMA模型作為點估計模型，其不僅有最好的擬合度，而且MSE、MAE和MAPE均明顯小于其他3 種模型。

2.3 區(qū)間估計

2.3.1 區(qū)間估計WPD-ARIMA-GARCH預測結果

在經過WPD-ARIMA點估計模型預測后，本研究采用GARCH模型對其殘差進行區(qū)間估計，WPD-ARIMAGARCH預測結果如圖7所示。

圖7 某地區(qū)測試集WPD-ARIMA-GARCH的預測結果Fig.7 Prediction results of WPD-ARIMA-GARCH for test set

2.3.2 區(qū)間估計模型比較

為評價GARCH模型區(qū)間估計的效果，本研究分別采用KDE和GPR模型作為區(qū)間估計對比模型進行比較。WPD-ARIMA-KDE和WPD-ARIMA-GPR的預測結果如圖8、9所示。各種區(qū)間預測組合模型的預測結果如表6所示。

圖8 某地區(qū)測試集WPD-ARIMA-KDE的預測結果Fig.8 Prediction results of WPD-ARIMA-KDE for test set from a certain region

圖9 某地區(qū)測試集WPD-ARIMA-GPR的預測結果Fig.9 Prediction results of WPD-ARIMA-GPR for test set from a certain region

表6 基于不同區(qū)間預測方法的結果比較Table 6 Comparison of results based on different interval prediction methods

預測結果中的PINAW和CWC值越小，置信區(qū)間越窄；PICP值越大，置信區(qū)間涵蓋的目標值越多，模型的預測效果越好。實驗結果表明，GARCH模型區(qū)間估計預測效果明顯優(yōu)于KDE模型區(qū)間估計，同時，GARCH模型的90%和95%置信區(qū)間均可以覆蓋所有真實值，且90%置信區(qū)間的PINAW和CWC更小。雖然GPR的PINAW和CWC較低，但其PICP僅為0.58左右，這意味著其置信區(qū)間僅涵蓋近58%的真實值。綜合考慮準確性和置信區(qū)間，GARCH模型90%的置信區(qū)間效果最佳。

采用建立的WPD-ARIMA-GARCH模型，對測試集進行預測，結果如表7所示。結果表明，預測值均與實際真實值相近，且90%置信區(qū)間可以涵蓋所有真實值和預測值的結果。根據測試集的預測結果，結合公式（3）的計算結果以及專家打分法對風險的劃分，測試集中序號7和24的預測結果為高風險項，考慮本實驗中設置的lookback＝6，即需要對2019年第13周和第31周重點關注，也就是對2019年的3月底和7月底進行提前預警并采取預防措施；同時，2019年第8、16、18、21、27、34周的預測值在8和21之間，也應當給予適當?shù)年P注。該結果與實際食品抽樣檢測中風險一致，說明該預測結果準確且有效。

表7 WPD-ARIMA-GARCH的測試集預測結果Table 7 Prediction results of WPD-ARIMA-GARCH for test set

為了測試WPD-ARIMA-GARCH模型的通用性，本實驗采用所建立的模型，對其他10 個不同地區(qū)的醬鹵肉制品質量安全風險進行預測，預測結果的評價指標如表8所示。

表8 其他10 個地區(qū)醬鹵肉制品的預測結果Table 8 Prediction results of soy sauce and pot-roast meat products from 10 other different regions

結果表明，10 個地區(qū)的點估計平均誤差MSE、MAE和MAPE分別為1.626、0.806和20.824，90%置信區(qū)間下，區(qū)間估計預測PICP、PINAW和CWC分別為0.992、0.024和0.024，表明本研究構建的WPD-ARIMA-GARCH模型對醬鹵肉制品風險有較好的預測效果，該模型通過WPD和ARIMA的耦合，使其對高頻數(shù)據處理更加容易，提高預測精度，減小誤差，準確地模擬時間序列變量的波動性變化，不僅考慮了待分析數(shù)據在時間序列上的依存性，而且考慮了隨機波動的干擾，GARCH模型對殘差進行區(qū)間估計，進一步驗證結果的可靠程度，且提供每個確定預測結果的上限和下限，可以量化預測結果的不確定性。

3 結論

通過對原始數(shù)據分析和劃分，將點估計和區(qū)間估計結合起來，構建了WPD-ARIMA-GARCH組合模型，并創(chuàng)新性地應用于醬鹵肉制品的質量安全風險預測。在點估計方面，WPD-ARIMA模型與對照模型（WPDSVR、ARIMA、SVR）相比，不僅擬合度最好，而且MSE、MAE、MAPE均最小。在區(qū)間估計方面，與其他區(qū)間估計模型（如KDE、GPR）相比，基于GARCH模型的區(qū)間估計可以產生更窄的PINAW，同時可以保證更高的PICP，效果最佳。根據實驗建立的WPD-ARIMAGARCH組合模型，對某地區(qū)醬鹵肉制品進行風險預測，預測結果表明，2019年的3月底和7月底醬鹵肉制品安全風險較高，需要進行提前預警并采取必要的預防措施，同時，也應當給予2019年第8、16、18、21、27、34周適當?shù)年P注，該結果與實際抽樣檢測結果一致，因此，本實驗建立的模型準確且有效。同時該模型在10 個不同地區(qū)的平均MSE、MAE和MAPE分別為1.626、0.806和20.824，其90%置信區(qū)間的PINAW和CWC值均為0.024，該置信區(qū)間可以覆蓋所有真實值，具有高預測精度和較低的誤差，可以對醬鹵肉制品質量安全中潛在的風險起到防控和監(jiān)督的作用，并能夠在日常檢測的過程中提供相應的技術支持。