張維蔚，陳坤才＊，張玉華，陳燕珊，黃德演

（1.廣州市疾病預(yù)防控制中心，廣東廣州 510440）（2.廣州大學(xué)經(jīng)濟(jì)與統(tǒng)計(jì)學(xué)院，廣東廣州 510006）（3.廣東毓秀科技有限公司，廣東廣州 510623）

食品安全問題是影響國家發(fā)展的全球性的重要健康問題[1]。隨著生產(chǎn)水平的提高，市場上的食品種類也越來越多，食品安全監(jiān)管難度也隨之增大。食品安全事件仍時(shí)有發(fā)生，例如2018 年桂林學(xué)術(shù)交流會(huì)上500 余人餐后出現(xiàn)嘔吐、腹瀉、發(fā)燒等癥狀[3]和2021 年中央電視臺(tái)曝光“養(yǎng)羊大縣”青縣“瘦肉精”羊肉事件，該縣每年羊出欄量高達(dá)70 萬只[4]等。由此可見，食品安全問題依舊是國家治理的一個(gè)重要方向。

在所有食品種類中，肉及肉制品是國民日常所需的重要食品種類之一，是食品安全的重要監(jiān)管方向。相關(guān)數(shù)據(jù)顯示，我國2015 至2020 年間食品日常監(jiān)督管理抽檢數(shù)據(jù)中，肉及肉制品不合格占比高達(dá)8.00%[5]。劉明等[6]還指出肉及肉制品在大型活動(dòng)中，風(fēng)險(xiǎn)級別為中度風(fēng)險(xiǎn)。此外，肉及肉制品在我們?nèi)粘ｏ嬍骋约按笮突顒?dòng)中被廣泛食用，出現(xiàn)頻次超高，是需要重點(diǎn)研究的食品種類。

在食品安全的風(fēng)險(xiǎn)評估方向上，研究者提出了多種研究方法。除了常規(guī)的描述性統(tǒng)計(jì)分析[7]外，比較經(jīng)典的方法有蒙特卡羅法[8]、貝葉斯網(wǎng)絡(luò)[9]和層次分析法[10]，其中蒙特卡羅法和貝葉斯網(wǎng)絡(luò)分別為定量和定性評估，而層次分析法可以綜合定性和定量評估。除此之外，Lin 等[10]嘗試使用灰色關(guān)聯(lián)分析（GRA）綜合解釋結(jié)構(gòu)模型（ISM）來構(gòu)建風(fēng)險(xiǎn)指數(shù)，該方法剖析了檢測指標(biāo)間的相關(guān)性并繪制了結(jié)構(gòu)圖，使得構(gòu)建的風(fēng)險(xiǎn)指數(shù)更具客觀性。上面方法實(shí)際上都是通過生成指標(biāo)權(quán)重，繼而得到風(fēng)險(xiǎn)指數(shù)，對食品安全進(jìn)行評估，但這不能對風(fēng)險(xiǎn)趨勢進(jìn)行預(yù)測。因此結(jié)合數(shù)據(jù)時(shí)間特征及馬爾可夫鏈假設(shè)，研究者們[11]嘗試引入隱馬爾可夫模型（HMM）對食品安全進(jìn)行評估和預(yù)測，結(jié)果表明HMM 具有高準(zhǔn)確性和強(qiáng)時(shí)效性。在隱馬爾可夫模型研究中，估計(jì)隱狀態(tài)數(shù)目一直是一個(gè)難點(diǎn)。在這方面，劉鶴飛等[14]利用可逆跳躍蒙特卡羅算法給了HMM 的隱狀態(tài)數(shù)目的估計(jì)。Rousseeuw 等[15]使用譜聚類算法估計(jì)了離散狀態(tài)HMM 的參數(shù)。在此背景下，Zheng等[16]綜合了之前研究者的方法，根據(jù)觀測狀態(tài)首中時(shí)期望構(gòu)建特征矩陣，再由奇異值分解以及K 均值聚類確定隱狀態(tài)數(shù)目。

上述所闡述的方法已被應(yīng)用到許多種類食品的研究，但鮮有應(yīng)用到對于肉及肉制品的風(fēng)險(xiǎn)評估研究。鑒于此，本文針對廣州食品安全監(jiān)測數(shù)據(jù)關(guān)于肉及肉制品的檢測數(shù)據(jù)，綜合已有方法（即GRAISM-HMM），利用GRA 和ISM 構(gòu)建風(fēng)險(xiǎn)指數(shù)，運(yùn)用首中時(shí)期望特征矩陣對HMM 隱狀態(tài)進(jìn)行估計(jì)，從而對廣州市肉及肉制品安全進(jìn)行風(fēng)險(xiǎn)評估。

1 數(shù)據(jù)與方法

1.1 實(shí)驗(yàn)數(shù)據(jù)

本文數(shù)據(jù)來自廣州食品安全監(jiān)測2015 年至2021 年抽檢的肉及肉制品數(shù)據(jù)，樣品共806 份，檢測目的是檢測肉及肉制品中的化學(xué)污染物，檢測方法按照原國標(biāo)方法測定。檢測項(xiàng)目受每年食品安全政策及民生食品安全關(guān)注焦點(diǎn)影響，每年都有所變動(dòng)，每年具體抽檢指標(biāo)見表1 所示。

表1 2015年至2021年各年份檢測指標(biāo)表Table 1 Chemical detecting indicators from 2015 to 2021

在這些數(shù)據(jù)中，一些肉及肉制品種類具有不同的檢測方案，需要分類處理，去除未超出檢出限的指標(biāo)，各年份肉及肉制品分別的檢測指標(biāo)見表2 所示。

表2 2015至2021年不同肉及肉制品種類的有效檢測指標(biāo)Table 2 Effective detection indicators for different meat and meat product types from 2015 to 2021

1.2 實(shí)驗(yàn)環(huán)境

本研究基于Python 3.7.6 環(huán)境構(gòu)建模型，實(shí)驗(yàn)設(shè)備為Windows 11 筆記本電腦，Intel(R) Core? i5-11400H 處理器，NVIDIA GeForce RTX 3050 顯示適配器，16 G（3 200 MHz）內(nèi)存。

1.3 實(shí)驗(yàn)方法

1.3.1 灰色關(guān)聯(lián)模型

灰色關(guān)聯(lián)模型（GRA 模型）是一種多元分析方法，其原理是根據(jù)不同序列的曲線幾何形狀的相似程度來判斷其聯(lián)系是否緊密，曲線越接近，相應(yīng)序列之間的關(guān)聯(lián)度就越大，反之則越小。

首先設(shè)定待確定關(guān)聯(lián)度的兩個(gè)序列分別為比較序列X=（x1，x2，…xp）’和參考序列Y=（y1，y2，…yp）’，其中x1，x2，…xp，y1，y2，…yp分別為各因子，其樣本量均為n。由于各因子量綱不同，通常先進(jìn)行初值化處理。

式中：

xi（j）′——第i 個(gè)比較序列因子xi的第j 個(gè)樣本初值化結(jié)果；

xi（j）——第i 個(gè)比較序列因子xi的第j 個(gè)樣本；

yk（j）′——第k 個(gè)比較序列因子yk的第j 個(gè)樣本初值化結(jié)果；

yk（j）——第k 個(gè)比較序列因子yk的第j 個(gè)樣本。

接下來求解序列差異矩陣δ以及矩陣最大值δikmax和最小值δikmin，具體公式如下：

最后，對第j 個(gè)樣本，求解xi與yk的關(guān)聯(lián)系數(shù)Yik以及X 和Y 的關(guān)聯(lián)度矩陣C=(Cik)p×q，公式如下：

式中：

Cik——xi與yk的關(guān)聯(lián)度，關(guān)聯(lián)度值越大，則xi對yk的影響程度越大。

1.3.2 解釋結(jié)構(gòu)模型

解釋結(jié)構(gòu)模型（ISM）是一種用于分析復(fù)雜系統(tǒng)層級結(jié)構(gòu)的模型方法，該方法常用于將影響目標(biāo)變量的各因素進(jìn)行分層處理。方法步驟如下：

（1）構(gòu)建影響因素的影響矩陣A，即鄰接矩陣。鄰接矩陣元素為“0”或“1”，其中“1”表示兩因素之間有直接連接關(guān)系。如果由相關(guān)性構(gòu)建，“1”通常表示為強(qiáng)相關(guān)，否則為“0”。

（2）構(gòu)建可達(dá)矩陣M。首先將鄰接矩陣A與單位矩陣E相加，然后進(jìn)行n 次布爾運(yùn)算，即當(dāng)滿足(A+E)n-1≠(A+E)n=(A+E)n+1停止運(yùn)算，由此得到可達(dá)矩陣M=(A+E)n。其中可達(dá)矩陣第i 行第j 列元素Mij為“1”表示影響因素i 與影響因素j 之間存在路徑，否則表示不存在路徑。

（3）構(gòu)建可達(dá)集R、先行集Q 與它們的交集B?？蛇_(dá)集表示為M 中某影響因素對應(yīng)行中，元素為1的因素集合。先行集表示為M 中某影響因素對應(yīng)列中，元素為1 的因素集合。

（4）構(gòu)建相關(guān)層次。首先確定第一層因素，即可達(dá)集只有本身的影響因素，然后刪除該列，再次搜索可達(dá)集只有本身的影響因素放入第二層，直至最后一個(gè)因素被剔除。

ISM 用分層結(jié)構(gòu)顯式地展現(xiàn)了各影響因素的關(guān)系，為對各因素賦予權(quán)重提供了有效的依據(jù)。

1.3.3 隱馬爾可夫模型

隱馬爾可夫模型（HMM）是經(jīng)典的概率統(tǒng)計(jì)模型，包含馬爾可夫鏈和一般隨機(jī)過程。它描述的是由隱藏的馬爾可夫鏈以一定概率生成狀態(tài)序列，再由狀態(tài)序列以一定概率隨機(jī)生成觀測序列的過程，圖像表示如下：

圖1 隱馬爾可夫模型圖Fig.1 Hidden Markov model

其中｛O1,O2,…,OT｝表示不同時(shí)刻的觀測序列，｛q1,q2,…,qT｝表示不同時(shí)刻的隱藏狀態(tài)序列。隱狀態(tài)之間能夠相互轉(zhuǎn)換，但其不能被觀測，隱狀態(tài)以一定概率生成觀測序列使得我們能觀測到，因此HMM 的一個(gè)重要任務(wù)就是去求解這些轉(zhuǎn)換概率。若觀測也是離散狀態(tài)，先假設(shè)共有N個(gè)隱藏狀態(tài)｛R1,R2,…,RN｝，共有M 類離散觀測｛u1,u2,…uM｝，則模型參數(shù)可以由λ=(π,Λ,B)表示，其中π=(π1,π2,…πN)表示初始時(shí)刻的狀態(tài)分布概率矩陣，Λ=(αij)N×N表示隱狀態(tài)之間的轉(zhuǎn)移概率矩陣，其中aij=P(qt+1=Rj|qt=Ri),1≤i≤N,1≤j≤N,表示任意t時(shí)刻隱藏狀態(tài)為Ri時(shí)，下一時(shí)刻轉(zhuǎn)換成Rj的概率。B 表示隱狀態(tài)對于觀測的發(fā)射矩陣，bij=P(ot+1=uj|qt=Ri),1≤j≤M,1≤i≤N表示任意t時(shí)刻隱狀態(tài)為Ri時(shí)，觀測狀態(tài)為uj的概率。

在估計(jì)相關(guān)概率前，需要先確定隱狀態(tài)的數(shù)目N。如果兩個(gè)觀測狀態(tài)ui和uj是來自同一個(gè)隱狀態(tài)，那么ui到某個(gè)固定點(diǎn)的首中時(shí)（即首次命中時(shí)間）與uj到這個(gè)固定點(diǎn)的首中時(shí)分布相同。基于以上事實(shí)，有這樣一個(gè)定理：假設(shè)HMM 的轉(zhuǎn)移概率矩陣是遍歷的，若首中時(shí)期望矩陣有s個(gè)非零特征值，則該HMM 模型有s個(gè)隱狀態(tài)。這是由于若有兩個(gè)觀測狀態(tài)屬于同一個(gè)隱狀態(tài)，那它們所在的首中時(shí)期望矩陣行之間存在相關(guān)性，因此會(huì)存在0 特征值，所以非零特征值對應(yīng)隱狀態(tài)個(gè)數(shù)。參照Zheng等[16]所提出的估計(jì)算法，得到隱狀態(tài)數(shù)估計(jì)算法如表3 所示。Zheng 等[16]還對首中時(shí)期望均值進(jìn)行聚類分析來獲得隱狀態(tài)數(shù)目，但本文觀測狀態(tài)并不多，無需進(jìn)行聚類分析，所以運(yùn)用表3 的算法。

表3 HMM隱狀態(tài)數(shù)估計(jì)算法Table 3 Estimation algorithm of hidden states’ number of HMM

估計(jì)隱狀態(tài)數(shù)后，便可以HMM 建模。HMM主要解決三類問題：

（1）評估問題：當(dāng)給定觀測序列O=｛O1,O2,…,OT｝和已知模型參數(shù)λ時(shí)，計(jì)算P(O|λ)，即計(jì)算參數(shù)λ下，該觀測序列O被觀測到的概率。

（2）解碼問題：當(dāng)觀測序列和模型參數(shù)已知情況，計(jì)算最大可能出現(xiàn)的隱藏狀態(tài)序列。

（3）學(xué)習(xí)問題：已知觀測序列，求解模型參數(shù)λ使得P(O|λ)最大。

1.3.4 GRA-ISM-HMM聯(lián)合模型

該節(jié)闡述對以上三種模型的連接操作。首先通過灰色關(guān)聯(lián)分析可以得到各影響因素的相互影響程度，在本文指代各檢測指標(biāo)相互之間的關(guān)聯(lián)度。設(shè)定閾值為0.9[17]，當(dāng)關(guān)聯(lián)度超過閾值，判定這兩個(gè)檢測指標(biāo)之間有連接關(guān)系，進(jìn)而可以得到ISM 的鄰接矩陣。但在實(shí)驗(yàn)中發(fā)現(xiàn)大多數(shù)檢測指標(biāo)關(guān)聯(lián)度都非常大，這導(dǎo)致根據(jù)原ISM 的分層困難，因此在參考經(jīng)典ISM 的同時(shí)，增加一個(gè)條件即隨著分層結(jié)構(gòu)箭頭指向，關(guān)聯(lián)度應(yīng)逐級遞減。增加該條件后，便可以得到分層結(jié)果，由此我們可以逐級賦予權(quán)重，權(quán)重逐級遞減，從而得到各年份的風(fēng)險(xiǎn)指數(shù)。采用分位數(shù)四分法對風(fēng)險(xiǎn)指數(shù)進(jìn)行分級，可以得到離散觀測狀態(tài)，接著再使用HMM 便可以得到觀測下各年份隱藏的風(fēng)險(xiǎn)等級情況，具體流程圖如圖2 所示。

圖2 GAR-ISM-HMM 流程圖Fig.2 GAR-ISM-HMM flow chart

2 結(jié)果與分析

由于數(shù)據(jù)量綱不同，對數(shù)據(jù)進(jìn)行歸一化處理后進(jìn)行灰色關(guān)聯(lián)分析。若某指標(biāo)所有檢測結(jié)果都低于檢出限，歸一化后便都為0，應(yīng)不進(jìn)行灰色關(guān)聯(lián)分析。以2015 年為例，有效檢測指標(biāo)為鉛、鎘、鉻、喹乙醇及其代謝物、呋喃唑酮、呋喃妥因六個(gè)，得到的灰色關(guān)聯(lián)矩陣結(jié)果如圖3 所示。

圖3 2015 年灰色關(guān)聯(lián)分析結(jié)果圖Fig.3 Results of grey correlation analysis in 2015

根據(jù)灰色關(guān)聯(lián)度矩陣，可以利用ISM 構(gòu)建分層結(jié)果，得到分層結(jié)構(gòu)圖如圖4 所示，以2015 年為例。

圖4 2015 年分層結(jié)構(gòu)圖Fig.4 Hierarchical structurein 2015

表4 2015年各指標(biāo)權(quán)重表Table 4 Weights of each indicator in 2015

根據(jù)賦予的權(quán)重便可得到各年份加權(quán)綜合風(fēng)險(xiǎn)指數(shù)。各年份風(fēng)險(xiǎn)指數(shù)結(jié)果如圖5 所示。

圖5 2015-2020 年份綜合風(fēng)險(xiǎn)指數(shù)圖Fig.5 Composite risk index from 2015 to 2020

圖5 中“all”是所有年份風(fēng)險(xiǎn)指數(shù)數(shù)據(jù)的匯總，用來觀察所有樣本的數(shù)據(jù)分布情況。由于使用歸一化處理，且每年各指標(biāo)權(quán)重和為1，因此如果一個(gè)樣本各指標(biāo)數(shù)值都很大的情況下，加權(quán)綜合風(fēng)險(xiǎn)指數(shù)會(huì)越接近1。也就是說，數(shù)值“1”可作為風(fēng)險(xiǎn)指數(shù)結(jié)果參考頂點(diǎn)。首先從圖5 可以看出，所有樣本的風(fēng)險(xiǎn)指數(shù)都未超過0.5，且大多數(shù)樣本風(fēng)險(xiǎn)指數(shù)不超過0.2，由此可得廣州每年的肉及肉制品安全處于較低的風(fēng)險(xiǎn)水平，但依舊有極少樣本風(fēng)險(xiǎn)較高，關(guān)乎健康安全，這些樣本也需要重點(diǎn)關(guān)注。其次，從圖5 可以明顯看出，2019 年加權(quán)綜合風(fēng)險(xiǎn)指數(shù)相對明顯較高。但由于廣州肉及肉制品生化檢測還在逐步完善的過程，每年檢測指標(biāo)變動(dòng)較大，且檢測樣本種類也有改變。所以，2019 年風(fēng)險(xiǎn)指數(shù)高的原因需要進(jìn)一步探討。探討前我們可以先根據(jù)HMM模型探索該觀測下各年份的隱藏狀態(tài)，觀察隱藏風(fēng)險(xiǎn)狀態(tài)與觀測是否一致。

在進(jìn)行離散HMM 建模前，需要先將得到的風(fēng)險(xiǎn)指數(shù)劃分風(fēng)險(xiǎn)等級，根據(jù)樣本特點(diǎn)，在這里采用分位數(shù)劃分法。因?yàn)榇嬖谙鄬Ξ惓８叩娘L(fēng)險(xiǎn)指數(shù)樣本，所以傳統(tǒng)的5-標(biāo)度法（5 等分）并不能很好的劃分風(fēng)險(xiǎn)等級，而且分位數(shù)劃分法能夠使得各樣本落入各等級的概率相對均等。對5-分位數(shù)劃分、4-分位數(shù)劃分和3-分位數(shù)劃分結(jié)果進(jìn)行比較，最終選擇4-分位數(shù)劃分，即劃分為四個(gè)等級。最終劃分結(jié)果如表5 所示。

表5 風(fēng)險(xiǎn)等級劃分表Table 5 Risk classification

由于HMM 的離散觀測從“0”開始，表5 將第1 個(gè)等級劃分為“0”。根據(jù)表5 的風(fēng)險(xiǎn)等級以及各年份平均風(fēng)險(xiǎn)指數(shù)，2015 年至2020 年六年的觀測狀態(tài)序列為 [122231]。

得到觀測序列后，開始估計(jì)HMM 隱狀態(tài)數(shù)目。計(jì)算各個(gè)觀測狀態(tài)的之間的首中時(shí)，得到矩陣如下：

矩陣（8）中“∞”表示兩種狀態(tài)相互轉(zhuǎn)換概率為0。顯然矩陣（8）含非0 特征值3 個(gè)，因此隱狀態(tài)個(gè)數(shù)應(yīng)取3 個(gè)為宜。

對HMM 模型進(jìn)行多次訓(xùn)練，取得分最高的參數(shù)估計(jì)結(jié)果，結(jié)果如下所示：

π、Λ、B的概率都為四舍五入后的概率，如矩陣中的“0”實(shí)際結(jié)果均＜0.000 1，“1”實(shí)際結(jié)果均＞0.999 9，為了方便記錄，矩陣?yán)镉洖椤?”和“1”，從矩陣B 可以得到觀測狀態(tài)與隱狀態(tài)的明顯對應(yīng)關(guān)系，即隱狀態(tài)“0”對應(yīng)觀測狀態(tài)“2”，隱狀態(tài)“1”對應(yīng)觀測狀態(tài)“1”，隱狀態(tài)“2”對應(yīng)觀測狀態(tài)“3”。對得到的最高概率隱狀態(tài)序列 [1 0 0 0 2 1]進(jìn)行重新排級，可以得到隱藏的風(fēng)險(xiǎn)等級為 [0 1 1 1 2 0]。從中可以看出2019 年的風(fēng)險(xiǎn)確實(shí)相對較高，但2020 隱藏風(fēng)險(xiǎn)等級為“0”，且結(jié)果顯示2021 的風(fēng)險(xiǎn)等級較大的概率往風(fēng)險(xiǎn)等級“1”轉(zhuǎn)換，風(fēng)險(xiǎn)等級“1”為中等風(fēng)險(xiǎn)，對應(yīng)觀測狀態(tài)“2”，即預(yù)測2021 年的風(fēng)險(xiǎn)指數(shù)均值大概率落在 [0.047 8,0.118 3]之間。因此可以認(rèn)為廣州肉及肉制品風(fēng)險(xiǎn)預(yù)測趨勢良好。

利用網(wǎng)絡(luò)爬蟲采集外籍人才的網(wǎng)絡(luò)招聘信息，進(jìn)行探索性調(diào)查。對大量招聘信息進(jìn)行文本挖掘，整理出有用的信息，為企業(yè)問卷設(shè)計(jì)和抽樣提供依據(jù)。

為了驗(yàn)證這一結(jié)果，我們同樣利用GRA-ISM 模型求解2021 年的綜合風(fēng)險(xiǎn)指數(shù)，得到2021 年綜合肉及肉制品綜合風(fēng)險(xiǎn)指數(shù)的箱線圖結(jié)果如圖6 所示。

圖6 2021 年綜合風(fēng)險(xiǎn)指數(shù)圖Fig.6 Composite risk index in 2021

從圖6 中可以2021 年肉及肉制品的所有樣本的綜合風(fēng)險(xiǎn)指數(shù)大多分布在 [0,0.1]之間，總體風(fēng)險(xiǎn)較低。另一方面，對于各年份觀測到的風(fēng)險(xiǎn)等級，我們以平均風(fēng)險(xiǎn)指數(shù)劃分，因此求解2021 年平均風(fēng)險(xiǎn)指數(shù)，得到結(jié)果約為0.090 0，在 [0.047 8,0.118 3)內(nèi)。根據(jù)表5 可得，2021 年風(fēng)險(xiǎn)等級為“2”即觀測狀態(tài)“2”，與HMM 預(yù)測結(jié)果一致。

上文也說到各年份檢測指標(biāo)變動(dòng)較大，檢測種類也有所變化，所以我們也應(yīng)該重點(diǎn)探討2019 年為何風(fēng)險(xiǎn)提升，是否當(dāng)前檢測策略需要提升。下面對2019 年風(fēng)險(xiǎn)指數(shù)較高進(jìn)行逐一分析。首先，2019 檢測指標(biāo)較多，較多的危害指標(biāo)可能導(dǎo)致風(fēng)險(xiǎn)指數(shù)增加，因此考慮采取相同檢測指標(biāo)對各年份進(jìn)行分析。根據(jù)分析數(shù)據(jù)可得，各年份都檢測的指標(biāo)有金屬污染物指標(biāo)，即鉛、鎘、鉻三個(gè)指標(biāo)。根據(jù)前面分析的ISM 結(jié)構(gòu)模型，設(shè)定三個(gè)指標(biāo)權(quán)重分別為 [0.36,0.31,0.33]。最終各年份的風(fēng)險(xiǎn)指數(shù)圖如圖7 所示。

圖7 2015～2020 年份重金屬污染風(fēng)險(xiǎn)指數(shù)圖Fig.7 Heavy metal pollution risk index from 2015 to 2020

從圖7 中可以看出，2019 年金屬污染物風(fēng)險(xiǎn)指數(shù)依舊較高，但相對綜合風(fēng)險(xiǎn)指數(shù)，各年份風(fēng)險(xiǎn)相對有所變化，2019 年與各年份風(fēng)險(xiǎn)指數(shù)差異變小，即可以認(rèn)為風(fēng)險(xiǎn)指數(shù)變動(dòng)與檢測指標(biāo)變動(dòng)有一定因果關(guān)系，但還有其他因素影響。然而，從另一方面講，這也說明這幾年對風(fēng)險(xiǎn)檢測指標(biāo)的探索是正確的，適當(dāng)增加檢測指標(biāo)能更好挖掘安全風(fēng)險(xiǎn)。此外，2018 年重金屬化學(xué)污染風(fēng)險(xiǎn)較小，但綜合風(fēng)險(xiǎn)較大。

前文也說到各年份的檢測種類也有所變化?，F(xiàn)對數(shù)據(jù)中的包裝形式、采樣地點(diǎn)、產(chǎn)地、肉及肉制品分類等進(jìn)行分析，發(fā)現(xiàn)2019 年肉及肉制品檢測分類與其他年份相比有較明顯的差異，2019 年檢測的肉及肉制品分類增加了臘肉、臘腸兩個(gè)分類。因此，我們有理由懷疑肉及肉制品種類是影響肉制品風(fēng)險(xiǎn)變動(dòng)的原因之一，對各分類的風(fēng)險(xiǎn)指數(shù)進(jìn)行均值和方差計(jì)算，得到結(jié)果如表6 所示。

表6 各肉制品分類風(fēng)險(xiǎn)指數(shù)均值與方差表Table 6 Sample mean and sample variance of risk index for each meat product classification

從表6 可以看出各肉及肉制品分類的風(fēng)險(xiǎn)指數(shù)之間確實(shí)存在差異，其中臘肉、臘腸確實(shí)風(fēng)險(xiǎn)指數(shù)較高，這說明了2019 年增加臘肉、臘腸兩個(gè)分類確實(shí)是風(fēng)險(xiǎn)指數(shù)增加的原因之一。而除了臘肉、臘腸，其中雞肉的風(fēng)險(xiǎn)指數(shù)也較高，但雞肉每年都是重點(diǎn)檢測項(xiàng)目，不是年度差異的原因。

在現(xiàn)有對肉及肉制品的風(fēng)險(xiǎn)研究當(dāng)中，國內(nèi)外大多數(shù)研究都主要集中在對菌落類的研究[18]，對肉及肉制品中的化學(xué)風(fēng)險(xiǎn)物質(zhì)研究較少，且多具體為對某一種類（例如雞肉[20]）或某一類化學(xué)物質(zhì)（例如亞硝胺[21]）的研究，并沒有與本文類似的研究包括重金屬、農(nóng)藥殘留等綜合化學(xué)風(fēng)險(xiǎn)物質(zhì)的報(bào)道，但依舊有一些對比價(jià)值，具體對比研究如下。以2011 至2014 年全國肉及肉制品數(shù)據(jù)為例，Xin等[22]使用矩陣風(fēng)險(xiǎn)模型對全國肉及肉制品進(jìn)行風(fēng)險(xiǎn)評估，最終得出：燉肉和肉干是風(fēng)險(xiǎn)較高的肉制品種類，而細(xì)菌數(shù)目（菌落總數(shù)、大腸桿菌組）和防腐劑（包括山梨酸、苯甲酸等）是主要風(fēng)險(xiǎn)項(xiàng)目。對比發(fā)現(xiàn)，本文研究得到臘肉、臘腸是主要風(fēng)險(xiǎn)較高的種類，與該文結(jié)果相對一致；而在風(fēng)險(xiǎn)檢測項(xiàng)目上，從表1 可以看出，廣州肉及肉制品監(jiān)測在2018 到2019 年增加了防腐劑（包括山梨酸、苯甲酸等）的檢測，但在2020 剔除了這類指標(biāo)，這是由于我們的數(shù)據(jù)中，2018 年有鹵肉、叉燒等種類，2019 年有臘肉、臘腸等種類，對這些腌制品種類，山梨酸、苯甲酸等為主要風(fēng)險(xiǎn)項(xiàng)目，而2020 年所收集數(shù)據(jù)為在雞鴨肉生肉種類，因此剔除這類指標(biāo)。這也說明了我們的數(shù)據(jù)盡管包含多種種類，但對這些種類的檢測指標(biāo)是具有針對性的。對于雞肉類，劉峰等[23]通過高效液相色譜-串聯(lián)質(zhì)譜法和食品安全指數(shù)法，分析了2016～2020 年寧夏市售雞肉和雞蛋中獸藥殘留狀況，結(jié)果發(fā)現(xiàn)寧夏雞肉中獸藥殘留主要風(fēng)險(xiǎn)物質(zhì)為強(qiáng)力霉素，總體來說，寧夏雞肉和雞蛋的風(fēng)險(xiǎn)較低，但存在一定的風(fēng)險(xiǎn)隱患，但該文未考慮雞肉中獸藥殘留狀況在不同年份的變動(dòng)情況。而在本文中，我們主要研究不同年份肉及肉制品的綜合風(fēng)險(xiǎn)情況。此外，廣州每年檢測指標(biāo)有所變動(dòng)，以強(qiáng)力霉素為例，見表1 可知廣州監(jiān)測數(shù)據(jù)只在2019 年雞肉種類對其檢測，樣本少且只在一個(gè)時(shí)間段，難以分析，且2020 年剔除了強(qiáng)力霉素這個(gè)指標(biāo)。究其原因，細(xì)化到具體數(shù)據(jù)可以發(fā)現(xiàn)2019 年只有一個(gè)樣本的強(qiáng)力霉素檢測結(jié)果異常，其數(shù)值為9.28，其他都為1，即廣州從2019 年抽檢樣本分析認(rèn)為強(qiáng)力霉素該指標(biāo)不是影響風(fēng)險(xiǎn)的重要指標(biāo)，這與寧夏數(shù)據(jù)分析結(jié)果有所差異，可能是地域差異導(dǎo)致。在臘肉分類上，通過收集我國某縣城集貿(mào)市場和其周邊5 個(gè)鄉(xiāng)鎮(zhèn)集市的臘肉樣本數(shù)據(jù)，朱聯(lián)旭等[24]分析了臘肉中鉻（Cr）、鉛（Pb）、汞（Hg）、砷（As）和鎘（Cd）5 種重金屬和苯并 [a]芘的含量在不同臘肉斷面的分布規(guī)律，并對其進(jìn)行風(fēng)險(xiǎn)評估，得出結(jié)論：食用臘肉不會(huì)因?yàn)橹亟饘俸捅讲?[a]芘導(dǎo)致嚴(yán)重健康風(fēng)險(xiǎn)。本文對各年份肉及肉制品數(shù)據(jù)進(jìn)行重金屬風(fēng)險(xiǎn)評估，發(fā)現(xiàn)2019 年較其他年份的重金屬風(fēng)險(xiǎn)差異小于綜合風(fēng)險(xiǎn)差異，而2019 年增加了臘肉、臘腸種類，這也側(cè)面驗(yàn)證了臘肉、臘腸的重金屬風(fēng)險(xiǎn)較小，與朱聯(lián)旭等[7]的結(jié)果不沖突。綜上，本文研究結(jié)果與現(xiàn)有研究沒有沖突。

3 結(jié)論

從2015 至2020 年，廣州肉及肉制品化學(xué)污染物綜合風(fēng)險(xiǎn)指數(shù)總體較低，得到的風(fēng)險(xiǎn)指數(shù)都在 [0,0.45]之間。其中，2019 年風(fēng)險(xiǎn)最高，但這與2019 年檢測指標(biāo)變化以及增加臘肉、臘腸兩類肉制品有關(guān)。

HMM 結(jié)果顯示廣州2015～2020 年觀測到的風(fēng)險(xiǎn)等級為與隱藏風(fēng)險(xiǎn)等級一致，同時(shí)也顯示2019年風(fēng)險(xiǎn)異常。此外，HMM 預(yù)測顯示廣州肉及肉制品化學(xué)污染風(fēng)險(xiǎn)在往良好的態(tài)勢發(fā)展。

從2015 至2020 年，與綜合風(fēng)險(xiǎn)指數(shù)相比，肉及肉制品重金屬化學(xué)污染物風(fēng)險(xiǎn)各年份間差異較小，且較綜合風(fēng)險(xiǎn)指數(shù)有一定變化，特別是2018 年，重金屬化學(xué)污染風(fēng)險(xiǎn)較小，但綜合風(fēng)險(xiǎn)較大。這說明在2018 年，除重金屬外其他檢測指標(biāo)對綜合風(fēng)險(xiǎn)指數(shù)有較高影響。

不同肉及肉制品種類綜合風(fēng)險(xiǎn)指數(shù)有一定差異，其中雞肉、臘肉、臘腸風(fēng)險(xiǎn)相對較高，是需要重點(diǎn)關(guān)注的品類。

數(shù)據(jù)上，檢測數(shù)據(jù)中含有較多未超過檢出限的數(shù)據(jù)，這也顯示了廣州市場上肉及肉制品某些檢測指標(biāo)控制得不錯(cuò)。

廣州肉及肉制品化學(xué)污染檢測依舊有需要改進(jìn)，例如，每年檢測指標(biāo)變動(dòng)給風(fēng)險(xiǎn)評估帶來難度，因而需要盡快從數(shù)據(jù)中獲取信息，確立好良好的檢測方案。再如每年檢測種類也有所變動(dòng)，則需要確定重點(diǎn)風(fēng)險(xiǎn)種類，利用有限的資源有策略地開展檢測。另外，有些文獻(xiàn)顯示一些食品在不同季度有不同風(fēng)險(xiǎn)等級的情況[25]，因此在有能力的情況下可以對肉及肉制品也采集不同季度的樣本進(jìn)行比較。

文章結(jié)論雖然基于廣州市的檢測數(shù)據(jù)，但同樣可以給其他省市肉及肉制品監(jiān)管和風(fēng)險(xiǎn)評估提供借鑒。人民日益增長的美好生活需要對食品安全也提出了更高的要求，需要更有力地監(jiān)管和更精準(zhǔn)的風(fēng)險(xiǎn)評估方法。