基于WPD-ARIMA-GARCH組合模型的醬鹵肉制品安全風(fēng)險(xiǎn)區(qū)間預(yù)測(cè)

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

（1.湖北省食品質(zhì)量安全監(jiān)督檢驗(yàn)研究院，國(guó)家市場(chǎng)監(jiān)管重點(diǎn)實(shí)驗(yàn)室（動(dòng)物源性食品中重點(diǎn)化學(xué)危害物檢測(cè)技術(shù)），湖北省食品質(zhì)量安全檢測(cè)工程技術(shù)研究中心，湖北 武漢 430075；2.武漢理工大學(xué)計(jì)算機(jī)與人工智能學(xué)院，湖北 武漢 430070；3.武漢理工大學(xué)理學(xué)院，湖北 武漢 430070）

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

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

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

1 材料與方法

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

1.1 數(shù)據(jù)預(yù)處理

將檢測(cè)項(xiàng)目的檢測(cè)值（微生物項(xiàng)目除外）結(jié)合國(guó)家標(biāo)準(zhǔn)（例如GB 2760—2014《食品安全國(guó)家標(biāo)準(zhǔn) 食品添加劑使用標(biāo)準(zhǔn)》等）采用公式（1）進(jìn)行去量綱化處理。根據(jù)去量綱化的結(jié)果結(jié)合專家打分法將項(xiàng)目風(fēng)險(xiǎn)等級(jí)劃分為5 級(jí)，1級(jí)為安全無(wú)風(fēng)險(xiǎn)，2級(jí)為輕微風(fēng)險(xiǎn)，3級(jí)為輕度風(fēng)險(xiǎn)，4級(jí)為中度風(fēng)險(xiǎn)，5級(jí)為重度風(fēng)險(xiǎn)（即不符合國(guó)家標(biāo)準(zhǔn)要求）。其中較低的檢測(cè)值（1～2級(jí)）視為相對(duì)安全的情況，3級(jí)需要引起關(guān)注，而4～5級(jí)則需要進(jìn)行早期預(yù)警，并采取相應(yīng)的預(yù)防措施。具體劃分標(biāo)準(zhǔn)見(jiàn)表1。

表1 檢驗(yàn)項(xiàng)目的風(fēng)險(xiǎn)等級(jí)劃分標(biāo)準(zhǔn)Table 1 Criteria for risk classification of inspection items

式中：Yi為預(yù)處理后的風(fēng)險(xiǎn)等級(jí)評(píng)價(jià)值；Xstandard為國(guó)家標(biāo)準(zhǔn)中規(guī)定的標(biāo)準(zhǔn)值；Xi為檢驗(yàn)項(xiàng)目的檢測(cè)值。

微生物項(xiàng)目參考表1按照檢測(cè)值進(jìn)行風(fēng)險(xiǎn)等級(jí)劃分。

本實(shí)驗(yàn)采用改進(jìn)的softmax函數(shù)計(jì)算單一產(chǎn)品的綜合風(fēng)險(xiǎn)等級(jí)（公式（2）），通過(guò)函數(shù)中指數(shù)權(quán)重的變化來(lái)調(diào)節(jié)風(fēng)險(xiǎn)等級(jí)的權(quán)重[23]，將自然周作為一個(gè)數(shù)據(jù)集，通過(guò)加權(quán)求和（公式（3））計(jì)算自然周的綜合風(fēng)險(xiǎn)指數(shù)。

式中：level(A)為單一產(chǎn)品A的綜合風(fēng)險(xiǎn)等級(jí)；i為食品A中檢測(cè)項(xiàng)目的風(fēng)險(xiǎn)等級(jí)數(shù)值；ωi為風(fēng)險(xiǎn)等級(jí)i在食品A中的占比。

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

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

1.2 缺失數(shù)據(jù)插值

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

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

以某地區(qū)的數(shù)據(jù)作為插值實(shí)驗(yàn)的數(shù)據(jù)集，該地區(qū)周綜合風(fēng)險(xiǎn)的數(shù)據(jù)情況如圖1所示。

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

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

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

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

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

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

1.2.2 插值方法

針對(duì)數(shù)據(jù)缺失問(wèn)題，本實(shí)驗(yàn)采用GPR進(jìn)行插值。

GPR使用多維高斯過(guò)程條件分布的性質(zhì)來(lái)進(jìn)行回歸預(yù)測(cè)，其具體的表達(dá)如式（4）所示：

式中：Xi為特征向量；yi為對(duì)應(yīng)的值。

對(duì)于一個(gè)新的觀測(cè)值X*，GPR給出的預(yù)測(cè)值y*如式（5）所示：

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

1.2.3 插值方法對(duì)比

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

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

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

c）數(shù)據(jù)樣本應(yīng)當(dāng)具有代表性。

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

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

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

1.3.1 WPD

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

圖2 三級(jí)WPD示意圖Fig.2 Schematic diagram of three-scale WPD

1.3.2 ARIMA模型

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

1.3.3 GARCH模型

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

式中：μt為時(shí)間t的預(yù)測(cè)平均值；Z為Z分?jǐn)?shù)；α為置信水平；σt為時(shí)間t的標(biāo)準(zhǔn)偏差；n為樣本量。

1.3.4 WPD-ARIMA-GARCH組合模型

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

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

ARIMA-GARCH模型的運(yùn)行過(guò)程如表2所示。

表2 ARIMA-GARCH模型的運(yùn)行過(guò)程Table 2 Operation process of ARIMA-GARCH model

1.4 預(yù)測(cè)模型的評(píng)價(jià)指標(biāo)

1.4.1 點(diǎn)估計(jì)模型評(píng)價(jià)指標(biāo)

本實(shí)驗(yàn)選用以下3 個(gè)指標(biāo)——MSE、平均絕對(duì)誤差（mean absolute error，MAE）和平均絕對(duì)百分比誤差（mean absolute percentage error，MAPE）作為點(diǎn)估計(jì)模型的評(píng)價(jià)指標(biāo)，具體計(jì)算如式（8）～（10）所示。

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

1.4.2 區(qū)間估計(jì)模型評(píng)價(jià)指標(biāo)

本實(shí)驗(yàn)選用以下3 個(gè)指標(biāo)——預(yù)測(cè)區(qū)間覆蓋率（prediction interval coverage probability，PICP）、預(yù)測(cè)區(qū)間平均寬度（prediction interval normalized average width，PINAW）、覆蓋寬度標(biāo)準(zhǔn)（coverage width-based criterion，CWC）作為區(qū)間估計(jì)模型的評(píng)價(jià)指標(biāo)。具體計(jì)算如式（11）～（14）所示[35]。

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

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

式中：R為目標(biāo)值的取值范圍。

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

1.5 對(duì)比模型

1.5.1 點(diǎn)估計(jì)對(duì)比模型

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

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

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

2 結(jié)果與分析

2.1 插值方法結(jié)果比較

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

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

插值后該地區(qū)的周風(fēng)險(xiǎn)時(shí)間序列見(jiàn)圖4，其中2014—2018年的檢測(cè)數(shù)據(jù)作為訓(xùn)練集，2019年檢測(cè)數(shù)據(jù)作為測(cè)試集。

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

2.2 點(diǎn)估計(jì)

2.2.1 WPD-ARIMA模型

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

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

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

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

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

表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 點(diǎn)估計(jì)模型比較

為評(píng)價(jià)所構(gòu)建的WPD-ARIMA模型在醬鹵肉制品預(yù)測(cè)中的效果，本研究分別采用SVR、ARIMA、WPDSVR模型進(jìn)行對(duì)比分析。圖6和表5為某地區(qū)4 種點(diǎn)估計(jì)方法的預(yù)測(cè)結(jié)果。

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

表5 基于各種點(diǎn)估計(jì)方法的結(jié)果比較Table 5 Comparison of results based on various point estimation methods

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

2.3 區(qū)間估計(jì)

2.3.1 區(qū)間估計(jì)WPD-ARIMA-GARCH預(yù)測(cè)結(jié)果

在經(jīng)過(guò)WPD-ARIMA點(diǎn)估計(jì)模型預(yù)測(cè)后，本研究采用GARCH模型對(duì)其殘差進(jìn)行區(qū)間估計(jì)，WPD-ARIMAGARCH預(yù)測(cè)結(jié)果如圖7所示。

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

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

為評(píng)價(jià)GARCH模型區(qū)間估計(jì)的效果，本研究分別采用KDE和GPR模型作為區(qū)間估計(jì)對(duì)比模型進(jìn)行比較。WPD-ARIMA-KDE和WPD-ARIMA-GPR的預(yù)測(cè)結(jié)果如圖8、9所示。各種區(qū)間預(yù)測(cè)組合模型的預(yù)測(cè)結(jié)果如表6所示。

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

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

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

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

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

表7 WPD-ARIMA-GARCH的測(cè)試集預(yù)測(cè)結(jié)果Table 7 Prediction results of WPD-ARIMA-GARCH for test set

為了測(cè)試WPD-ARIMA-GARCH模型的通用性，本實(shí)驗(yàn)采用所建立的模型，對(duì)其他10 個(gè)不同地區(qū)的醬鹵肉制品質(zhì)量安全風(fēng)險(xiǎn)進(jìn)行預(yù)測(cè)，預(yù)測(cè)結(jié)果的評(píng)價(jià)指標(biāo)如表8所示。

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

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

3 結(jié)論

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