基于優化的灰色關聯分析-極限學習機食用油污染物風險評價模型研究

于家斌，范依云，王小藝,3, *，趙峙堯，金學波，白玉廷，王 立，陳慧敏

（1.北京工商大學人工智能學院，北京 100048；2.北京工商大學 中國輕工業工業互聯網與大數據重點實驗室，北京 100048；3.北京服裝學院文理學院，北京 100029）

糧油食品安全問題不僅影響著公眾健康，還關系到社會穩定和經濟發展等問題。近些年來，糧油食品安全事件屢見不鮮，如“鎘大米”[1]、“地溝油”[2]等，這些事件嚴重危害了公眾健康，影響了社會穩定。我國在2015年新修訂實施了《中華人民共和國食品安全法》，其中第二章食品安全風險監測與評價的第十七條規定了國家應構建食品安全風險評價體系[3]。而食用油是糧油食品的重要組成部分，是人們日常飲食的重要一環。因此，對食用油中各類污染物進行風險評價具有重要意義。

目前，國內外風險評價方法主要分為定性評價、定量評價以及定性-定量結合評價3 種。定性評價是依據專家經驗、調查問卷等方式對風險進行分析判斷得出結果，代表性的方法包括德爾菲法、決策樹、危害分析臨界控制點法[4-6]等，主要應用于食品監管、參與者對食品風險態度以及參與者食品安全意識分析等，定性評價過于依賴評估者的主觀判斷，無法基于大量高維的檢測數據進行風險評價。定量評價方法主要包括兩類，一類是通過檢測方法獲取風險因子數據直接描述風險水平，如氣相色譜法、光學傳感器法等[7]。另一類是通過數據驅動算法挖掘大量檢測數據中的風險規律，從而得到風險程度量化值，如熵權系數法、蒙特卡洛仿真法、支持向量機法、人工神經網絡、深度學習等[8-11]，但這類定量評價過于依賴數據，數據中存在噪聲會影響評價結果。定性-定量結合方法中具有代表性的方法有層次分析法，其主要是通過建立層次結構模型并結合專家意見構建風險判斷矩陣[12]，此外還有如云模型[13]、模糊綜合評價法[14]等方法依據隸屬度理論將定量風險值以及區間轉化為定性語言集，定性-定量結合評價法常用于屬性約簡[15]、風險分級[16]。綜上，現有風險評價方法在實際應用時仍存在一定局限性，如食用油檢測數據具有高維性、非線性、離散性的特征，導致傳統數據驅動方法評價效率和精度不高；其次，對于數據中的噪聲，現有方法并沒有針對性的處理，一般只是對數據進行簡單地清洗；另外，在利用數據驅動算法建模時，需設置多個參數，例如高亞男等[17]在采用LightGBM模型預測風險值時，需調整葉子節點個數、最大深度和學習率等參數，現有方法一般是利用人工試湊法來進行參數調優，這樣憑借經驗進行參數的試錯，是主觀且繁瑣的。

鑒于現有風險評價方法存在的問題，國內外專家學者進行了大量研究。針對數據中的噪聲問題，有學者利用基于小波理論的小波閾值法進行了濾波處理研究，例如Chen Jian等[18]利用小波軟閾值法對電力系統中低頻振蕩信號產生的噪聲進行濾波處理并取得良好效果。在針對食用油檢測數據這類以高維性、非線性、離散性為基礎的數據進行評價模型構建時，灰色關聯分析（grey relational analysis，GRA）能夠從大量復雜和非線性數據中分析不同指標之間的灰色關聯系數從而制定指標的權重，如Han Yongming等[19]通過對乳制品9 類指標的GRA確定了指標權重。極限學習機（extreme learning machine，ELM）作為一種單隱層前饋神經網絡，其隨機產生輸入層權值和隱層偏置的特性[20]，使其對數據中的噪聲有一定抑制能力，并且相比于傳統基于梯度下降的數據驅動評價算法具有更好的泛化能力和更快的速度，例如Zaher等[21]采用ELM預測泡沫混凝土的抗壓強度，效果優于支持向量機回歸等傳統數據驅動模型。針對模型參數設置主觀性強的問題，貝葉斯優化算法常用于超參數調整，是解決模型擬合問題的理想選擇，例如Acerbi等[22]提出一種混合貝葉斯優化算法并進行了模型擬合實驗，結果表明該算法要優于其他常見的優化算法。

綜上，針對食用油檢測數據的特點以及傳統風險評價方法中存在的問題，本實驗首先在風險評價模型前端加入基于小波閾值法的數據濾波模塊；接著結合GRA和ELM的優勢，對濾波后數據進行GRA得到每種風險指標權重并融合為多指標綜合風險值作為風險值標簽，利用ELM網絡訓練并輸出綜合風險值；在上述過程中利用實用貝葉斯優化（practical Bayesian optimization，PBO）算法分別優化濾波模塊和ELM網絡參數；最后對綜合風險值進行模糊綜合分析得到風險等級劃分結果。

1 材料與方法

1.1 數據特點以及預處理

本實驗收集整理了國內某食用油主要產出省2017—2019年共11 345 行、150 組花生油檢測數據，數據主要呈現以下特點[23]：1）高維性：每個食用油樣品的抽檢信息由不同屬性組成，如抽樣編號、產地信息、產品信息、生產日期、檢驗項目、檢驗結果、關鍵限值等。這些屬性相互獨立，并且很多屬性如抽樣編號等對于風險程度的判斷是冗余的，需要濾除掉。2）離散性：一個食用油樣品大約有幾十種檢測指標，一個指標維度與另一個指標維度不同，而各項指標的國家標準也不同。這意味著這些指標具有不同的離散域，檢測結果在離散域中也是無序的。因此，需要從這些指標中篩選出影響風險程度的關鍵指標。3）非線性：每個檢查信息包括數值信息和描述信息。數值信息包含大量空值，描述信息也包含了一些離散值。在完整的檢驗信息中，不同屬性之間的信息分布是不對稱和不平衡的，增加了風險分析難度。因此，需要進行指標篩選以及數據預處理工作。

首先進行指標篩選，食用油檢測指標種類繁多，包括酸值超標檢測、過氧化值檢測、溶劑殘留量檢測、金屬污染檢測、化學污染檢測和生物毒素檢測等。但并不是上述所有的檢測項都能作為食用油的風險評價指標，需要從其中篩選出對食用油風險程度影響較大的指標。本實驗分析各項指標的檢出情況，參考文獻[24]，選取酸值、過氧化值、砷含量、鉛含量、黃曲霉毒素B1含量和苯并[a]芘含量這6 項指標作為關鍵的風險評價指標。食用油的酸值和過氧化值是食用油的基礎類風險指標，砷、鉛是食用油中常見的重金屬類污染物，黃曲霉毒素B1則是食用油中常見的真菌毒素類污染物，苯并[a]芘是食用油中主要的化學類污染物，這6 項指標能夠代表相應類別的污染物，同時相對于其他風險指標更容易被檢出，因此通過這6 項指標能夠基本判斷食用油的風險程度。

其次進行數據預處理，步驟如下：1）刪除無用信息。每個樣本都有許多描述性字符，有必要簡化這些信息以強調影響風險程度的特性，剔除與食品安全關系不大的食品感官類指標，如口感、顏色、形狀等。2）刪除多余符號。例如某個食用油樣品總砷的檢測結果為“＜0.01 mg/kg”，總砷的限值不高于0.25 mg/kg，該樣品總砷的檢測結果并未超過限值，則去掉結果中的“＜”符號，將檢驗值記錄為“0.01 mg/kg”。3）使用極小值替代未檢出指標數據。本實驗對未檢出指標數據的處理參考文獻[15]，具體方法為針對檢測結果為未檢出的指標數據，利用極小的數字“0.001”代替“未檢出”，而不是用“0”代替。這是為了保證得到的輸入矩陣能夠被評價模型識別并且是正定的。至此指標的篩選以及數據預處理工作完成。

由于篇幅限值，僅列舉了部分原始檢測數據，如表1所示，對其進行預處理之后的數據如表2所示。

表1 食用油部分原始檢測數據Table 1 Selected raw chemical data of edible oil samples

表2 預處理后的食用油數據Table 2 Preprocessed chemical data of edible oil samples

1.2 風險評價模型

如圖1所示，本實驗提出的風險評價模型主要分為基于小波閾值法的數據濾波模塊和基于GRA-ELM的評價模塊，PBO算法則負責優化濾波模塊和ELM網絡中的重要參數，ELM網絡輸出的綜合風險值通過模糊綜合分析得到食用油樣品所屬風險等級。

圖1 風險評價模型整體算法流程圖Fig.1 Overall algorithm flow chart of risk assessment models

1.2.1 基于小波閾值的濾波模塊

食用油檢測數據可能會含有噪聲，這些噪聲會影響最終評價結果。Huber[25]對模型魯棒性進行了3 個層面的解釋：1）模型具有較高的精度或有效性；2）對于模型假設出現的較小偏差（如噪聲），只能對模型產生較小的影響；3）對于模型假設出現的較大偏差（如突變點），不會對模型性能產生較大的影響。實際進行風險因子含量檢測時，傳感器可能會產生散粒噪聲、熱噪聲、低頻噪聲等，其中散粒噪聲和熱噪聲影響較大，這兩種噪聲都是白噪聲。因此，為了充分驗證模型的魯棒性并盡可能模擬實際情況會產生的噪聲干擾，本研究對原始食用油檢測數據的噪聲方差進行了修改，如式（1）所示。

將上述得到的數據輸入到濾波模塊進行濾波處理，濾波模塊主體是小波閾值法，小波閾值法主要原理為設置一個臨界閾值，經小波變換后，所得小于閾值的小波系數主要是由噪聲引起，這部分系數需要被濾除掉，反之大于閾值的系數則是由信號引起，該部分系數需要被保留，最后通過小波重構得到濾波處理后的信號。本研究小波閾值法去噪分為3 個步驟：1）使用離散化小波變換將信號轉變到小波域；2）對各尺度上的小波系數做閾值量化處理；3）小波重構得到濾波后信號。

具體地，小波閾值去噪過程首先需進行離散小波變換（discrete wavelet transform，DWT）過程，DWT實際上是一個分解過程，分解則需選取小波分解層數和母小波函數類型，這個過程由貝葉斯優化算法自適應選取。在分解完成后會得到對應層數的高頻和低頻分量，這時需要選擇合適的閾值σ進行噪聲的平滑處理，經小波變換轉到小波域之后，白噪聲仍然呈現出較強的隨機性，因此小波域中更容易區分噪聲與信號，小波變換后各個分量中原始數據的有效信號對應較大的系數，而此時原始數據中的噪聲對應較小的系數。假設在小波域中噪聲的方差為σ，依據高斯分布特性，99.9%的噪聲系數都在[－υ，υ]范圍內，接著設定合適的閾值，通過比較信號絕對值和閾值的大小，將小于閾值的點重置為零，大于或等于閾值的點重置為該點信號與閾值的差，即實現了噪聲的平滑處理，上述過程如式（2）所示。

式中：Di,t表示第i層分解結果的高頻子序列索引為t的數據；λi是間接調節閾值的比例系數，其取值范圍為λi∈（0,1），i＝1,...,N；σ為噪聲估計方差；median表示被平滑序列的中位數；Di為分解后索引為i的高頻子序列；0.674 5為高斯噪聲標準方差調整系數；T為被去噪序列長度；γ為根據估計的噪聲方差所計算的每個子序列噪聲的估計閾值；而最終的閾值υ是通過比例系數λ的值調整估計閾值γ得到的；由于各個子序列中的噪聲含量不同，因此在平滑處理每個分量時對應的λ均不同，而每個λ則是通過PBO算法自適應選取的，這樣可以對不同數據中的噪聲進行平滑處理。具體PBO算法會在1.2.3節詳細介紹。優化最后對每層高頻分量濾波處理后，通過小波重構得到估計真值，如式（3）所示。

1.2.2 基于GRA-ELM的風險評價模塊

1.2.2.1 灰色關聯分析

GRA是一種多因素分析方法，它通過比較統計集合之間的幾何關系來劃分復雜系統中多因素之間的關系[26]。GRA主要包括以下步驟：1）利用GRA對指標進行賦權；2）計算綜合風險值。

首先需要獲取參考向量和比較向量，設參考向量為x1＝{x1(1),x1(2),...,x1(n)}，n為食用油樣本個數。比較向量為xj＝{xj(1),xj(2),...,xj(n)}，j＝1,2,...,m，m為風險指標個數。對m組指標進行歸一化處理，消除量綱的影響，如公式（4）所示。

然后計算關聯系數，k時刻的fj（k）和f1（k）的灰色關聯系數如公式（5）所示。

式中：ξj（k）為灰色關聯系數；ρ為調節參數，可以使各個系數的差異性增強，ρ∈（0，1）。序列f1和序列fj間的關聯系數見公式（6）。

為了保證結果的準確性，用每個風險指標各充當一次參考序列，則可得到所有風險指標的相關系數矩陣，見公式（7）。

本研究指標間關聯系數矩陣計算結果如表3所示，例如酸值對于自身灰色關聯系數為1，與過氧化值、砷含量、鉛含量、黃曲霉毒素B1含量和苯并[a]芘含量的灰色關聯系數分別為0.688 1、0.820 6、0.770 0、0.966 1和0.722 3。

表3 風險指標灰色關聯系數矩陣Table 3 Gray correlation coefficient matrix of risk indicators

根據表2可得權重向量，如公式（8）所示。

式中：m為風險指標個數；γij為表3中第i行指標與第j列指標之間的灰色關聯系數。

則權重向量W＝[0.173 980 203，0.150 610 512，0.172 256 898，0.168 421 495，0.174 607 178，0.160 123 714]。

接下來需要計算綜合風險評價值，首先需要進行無量綱化處理，具體為用風險因子的實際檢測值與指標關鍵限值α的比值表示相對風險值，如公式（9）所示。

式中：Pij為第j類風險因子的第i條檢測含量結果xij經無量綱化處理后得到的相對風險值，其中i＝1,2,...,n，j＝1,2,...,m；aj為第j類風險因子的指標關鍵限值。將公式（8）所得權重向量W與公式（9）中得到的各指標相對風險矩陣相乘，即得到公式（10）。

式中：Y＝[y1,y2,...,ym]T為綜合風險評價序列；P為相對風險矩陣；W＝[w1,w2,...,wm]T為權重向量。至此模型的期望輸出標簽計算完成。

1.2.2.2 極限學習機網絡

由Huang Guangbin等[27]提出的ELM是一種新型神經網絡學習算法，ELM的網絡結構與單隱層前饋神經網絡一樣，分為輸入層、隱含層和輸出層，但ELM的訓練階段不再是傳統神經網絡中基于梯度的算法（后向傳播），而是隨機生成輸入層權值和隱藏層偏置，輸出層權重則是通過最小化由訓練誤差項和輸出層權重范數的正則項構成的損失函數，并依據Moore-Penrose的廣義逆矩陣理論計算求出，得到所有節點的權值和偏差后即完成了ELM的訓練過程，測試數據利用求得的輸出層權重即可計算出網絡的預測輸出。

由于隨機生成輸入層權值和隱藏層偏置，ELM網絡不需要像傳統神經網絡那樣通過一次次迭代得到最終的解，因此ELM網絡的計算復雜度低，可調參數非常少，一般為隱層節點數，并且ELM網絡在保證精度的同時，學習速度更快、泛化能力更強。圖2為ELM的結構。

圖2 ELM網絡結構圖Fig.2 Network structure diagram of ELM

假設有N個任意樣本（xi,ti），其中Xi＝[xi1,xi2,...,xin]T∈Rn，ti＝[ti1,ti2,...,tin]T∈Rm，對于一個有L個隱層節點的單隱層神經網絡可以按公式（11）表示。

式中：g(x)為激活函數；Wi＝[wi1,wi,2,...,wi,n]T為輸入權重；βi為輸出權重；bi是第i個隱層節點偏置；Wi×Xj表示Wi和Xj的內積。

ELM學習的目標是使得輸出的誤差損失最小，可按式（12）表示。

式中：oj表示網絡輸出；tj表示目標輸出。

矩陣表示為式（13）。

式中：H是隱層節點的輸出；β是輸出權重；T是目標輸出。

這等價于最小化式（15）所示的損失函數。

在ELM中，輸入權重Wi和隱層偏置bi隨機生成后，隱層的輸出矩陣H就被唯一確定，ELM的訓練過程可以轉化為求解Hβ＝T，輸出權重β可以被確定，如式（16）所示。

式中：H?是矩陣H的Moore-Penrose廣義逆。在本研究中，原始的食用油安全檢測數據作為ELM的輸入，由公式（8）～（10）所得到的綜合風險評價數據作為期望輸出，對ELM進行訓練，對ELM隱層節點數的參數優化是由PBO算法完成的。

1.2.3 實用貝葉斯優化算法

由1.2.1節可知，小波閾值法中母小波函數、小波分解層數和每層小波高頻分量閾值是濾波模塊中的重要參數，由1.2.2.2節可知，ELM的隱層節點數目也是重要的參數[28]，這些參數會對模型最終的預測性能產生極大的影響，由于人工選取效率低且存在主觀性過強的問題，本研究利用PBO算法，對建模過程中的參數進行優化，該方法將網格自適應直接搜索（mesh adaptive direct search，MADS）與通過局部高斯過程（Gaussian process，GP）執行的BO搜索相結合，主要分為搜索和輪詢兩個階段。簡而言之，PBO是在一系列快速的、局部的BO步驟（MADS的搜索階段）和系統的網格探索（輪詢階段）之間交替進行，當搜索階段失敗時，意味著GP模型沒有成功優化參數，這時會切換到輪詢階段，輪詢階段執行的是無模型優化，在這個階段會收集優化目標的信息，以便下次在搜索階段構建更好的GP模型，直到達到優化目標。

在使用貝葉斯優化時，首先需定義目標函數，但由于ELM的訓練過程實際上是一個黑盒過程，因此本實驗中使用網絡的均方根誤差（root mean square error，RMSE）作為PBO算法的目標函數，具體如式（17）所示。

PBO算法偽代碼如表4所示。

表4 PBO算法偽代碼Table 4 PBO algorithm pseudocode

1.2.4 模糊綜合分析

由1.2.2節中評價模型所得的綜合風險評價值是對風險的定量評價，但為了給相關監管部門提供更直觀的參考依據，定性分析同樣重要[29]。模糊綜合分析來源于模糊數學中的隸屬度理論[30]，它能將待評價對象進行定性和定量分析間的相互轉換，在風險評價領域有極大的實用價值。對于本研究，模糊綜合分析主要分為兩步：1）構建模糊評語集；2）構建隸屬度函數。

食品風險因子的關鍵限值是由國家制定的用來判定食品質量是否滿足上市要求的統一標準，但單一的限量標準無法對食品質量的潛在風險進行評判，也無法進行多指標的綜合評價，例如，某個花生油樣品所選取的6 類風險因子均未超標，這并不能說明該樣品完全沒有風險。因此本研究依據5標度模式[31]，以國家標準為基礎設定5 個等級安全風險評語集，即低風險、中風險、警報風險、高風險和極高風險，每個風險等級對應的量化值分別為0.25a、0.5a、0.75a、a和2a，其中a代表國家限量標準，同時也是高風險對應的量化值，對于前4 個風險等級，以等距法設定對應量化值，為突出極高風險的危害性，將其對應量化值設為2 倍的國家限量標準。設立警報風險是由于模糊綜合分析的對象是多指標的綜合風險評價值，當達到警報風險時說明樣品的多個風險因子都十分接近國家限量標準，應該予以重視。

本實驗所檢測食用油的6 類風險因子均屬于負效應類指標，即指標檢測值需小于關鍵限值，適用于降半梯形隸屬度函數，對于1.3節中所得第i條綜合風險值yi，其隸屬度hiq如式（19）～（21）所示。

當q＝1時，

當q＝2,3,4時，

當q＝5時，

式中：Zq表示5標度評價等級所對應的量化值除以國家限量標準后的無量綱值，即{Z1,Z2,Z3,Z4,Z5}＝{0.25,0.5,0.75,1,2}，分別與模糊評語級低風險、中風險、警報風險、高風險和極高風險相對應，計算出第i條綜合風險值yi所對應5 個風險等級的隸屬度hiq后，根據最大隸屬度原則即可判斷出食用油所屬的風險等級。

2 結果與分析

實驗所用計算機為Windows10 64位操作系統，處理器為Intel(R) Core(TM) i5-10400 CPU @ 2.90 GHz，運行內存為32 GB，顯卡為NVIDIA RTX 2060（80 W），實驗基于Matlab R2015b軟件。為了驗證本實驗提出的基于優化的GRA-ELM風險評價模型的優越性和有效性，將經修改噪聲方差后的150 個花生油樣本風險指標數據作為模型輸入數據，將由式（8）～（10）得到的綜合風險評價值作為模型期望輸出數據，其中70%作為訓練集，30%作為測試集，首先將ELM網絡與傳統神經網絡模型（誤差反向傳播（error back propagation，BP）、徑向基函數（radial basis function，RBF））進行對比以證明ELM網絡的泛化能力和魯棒性，再分別測試濾波模塊和PBO算法的有效性。

2.1 不同情況性能對比

2.1.1 不同神經網絡性能對比

為了驗證ELM網絡的泛化能力和魯棒性，將修改噪聲方差后的原始檢測數據輸入RBF和BP網絡模型進行對比實驗，統一去掉濾波模塊，ELM網絡的輸入層、隱含層和輸出層的節點數分別設為6、20和1，單層BP網絡學習率設為0.1，動量因子為0.9，迭代次數為500，激勵函數為Sigmoid函數，隱含層節點設為20。RBF網絡的學習率設為0.1，迭代次數為500。評判指標為模型預測值與期望值的決定系數（coefficient of determination，R2）和RMSE，分別如式（22）、（23）所示。

GRA-ELM、GRA-RBF和GRA-BP模型回歸曲線如圖3所示，模型絕對誤差曲線如圖4所示，GRA-ELM、GRA-RBF和GRA-BP模型性能對比結果如表5所示，圖3中期望值曲線是利用式（8）～（10）對測試集樣本計算得到的。

表5 GRA-ELM、GRA-RBF和GRA-BP性能對比結果Table 5 Comparative performance of GRA-ELM,GRA-RBF and GRA-BP models

圖3 GRA-ELM、GRA-RBF和GRA-BP泛化情況對比Fig.3 Comparison of the generalization results of GRA-ELM,GRARBF and GRA-BP models

圖4 GRA-ELM、GRA-RBF和GRA-BP模型絕對誤差曲線Fig.4 Absolute error curves of GRA-ELM,GRA-RBF and GRA-BP models

從圖3可直觀觀察到，3 種模型都受到了噪聲的干擾，其中BP網絡受到的干擾最嚴重，無法準確地預測風險值。而表5的結果中RBF網絡的效果雖略優于BP網絡，但從圖3中可以看出，RBF的綜合風險值中由于噪聲的干擾出現了多個負值，是明顯的錯誤值。相比于BP和RBF網絡，沒有濾波模塊的ELM網絡雖在一定程度上會受到噪聲影響，但其預測值依然能與期望值保持較小的差距，由表5可知，ELM的RMSE以及R2分別為0.132 0、0.796 4，而BP、RBF的RMSE以及R2分別為0.193 8、0.373 8和0.172 0、0.576 2，圖4中不同模型絕對誤差曲線的對比更進一步說明ELM網絡在魯棒性和泛化能力上的優勢。分析其原因，傳統神經網絡如BP和RBF等會在一次次迭代訓練中對數據中存在的噪聲進行學習，得到的模型會因為過分學習訓練數據中的噪聲引起過擬合問題，而ELM網絡的輸入權值和隱層偏置是隨機生成的，與訓練數據無關，噪聲的影響會低很多。

2.1.2 濾波模塊和實用貝葉斯優化算法的有效性

為了驗證濾波模塊和PBO算法的有效性，對不加濾波模塊的模型（GRA-ELM）、加上濾波模塊的模型（WT-GRA-ELM）以及再經過PBO后的模型（WT-PBO-GRA-ELM）進行對比實驗，以經驗法確定的參數與PBO后參數對比如表6所示，模型對比結果如表7所示，模型回歸曲線如圖5所示，圖5中期望值曲線是利用式（8）～（10）對測試集樣本計算得到的。

表6 PBO前后的參數對比Table 6 Comparison of parameters before and after PBO

表7 WT-PBO-GRA-ELM、WT-GRA-ELM和GRA-ELM性能對比結果Table 7 Comparative performance of WT-PBO-GRA-ELM,WT-GRA-ELM and GRA-ELM models

圖5 WT-PBO-GRA-ELM、WT-GRA-ELM和GRA-ELM泛化情況對比Fig.5 Comparison of the generalization results of WT-PBO-GRAELM,WT-GRA-ELM and GRA-ELM models

由表7可知，加上濾波模塊后RMSE達到了0.094 5，相比無濾波時降低了28.41%，R2提升了12.04%，效果明顯。從圖5中也可看出，加上濾波模塊后ELM的擬合度更高，這充分說明濾波模塊對于最終風險評價結果準確率的重要性。而經過PBO后，總體性能得到進一步提升，RMSE和R2分別達到了0.056 3和0.946 1。由表6可知，經過PBO算法的參數優化，小波分解后每層高頻分量的閾值是不同的，這能夠更精確地濾除數據中的噪聲。從圖5可看出，WT-PBO-GRA-ELM模型對期望值曲線的擬合效果最好，從而說明了PBO算法的有效性。

2.2 風險等級劃分結果

針對表1中的食用油原始檢測數據，采用式（19）～（21）對本實驗所建立WT-PBO-GRA-ELM評價模型輸出的綜合風險值進行模糊綜合分析后，得到的風險隸屬度矩陣如表8所示，根據最大隸屬度原則可以得出各個樣本的風險等級：樣本1、3、4、5和6為低風險，樣本7為中風險，樣本2為極高風險。

表8 風險隸屬度矩陣Table 8 Risk membership matrix

以花生油樣本2為例，分別對GRA-ELM、WT-GRAELM模型輸出的綜合風險值進行模糊綜合分析，所得風險隸屬度如圖6所示，GRA-ELM、WT-GRA-ELM模型的風險分級結果都為高風險。而表1中樣本2的黃曲霉毒素B1檢測結果為20.60 μg/kg，是黃曲霉毒素B1標準限值的2 倍以上，由1.2.4節中國家限量標準倍數對應的風險語言集可知，2 倍標準限值對應的風險語言集為極高風險，同時，樣本2其他風險指標檢測結果也較為接近標準限值，因此樣本2的真實風險等級為極高風險，只有WT-PBOGRA-ELM模型判定結果與真實情況一致，說明本實驗所構建的風險評價模能準確識別出食用油風險程度以及優先次序，可以為監管部門制定有針對性的防護策略，并為確立優先監管領域和合理分配風險管理措施資源提供科學依據。

圖6 GRA-ELM、WT-GRA-ELM和WT-PBO-GRA-ELM隸屬度結果對比Fig.6 Comparison of the membership results of GRA-ELM,WT-GRA-ELM and WT-PBO-GRA-ELM models

3 結論

食用油檢測數據具有高維性、復雜性、非線性和離散性的特征，同時，在測定食用油各風險指標含量時可能會出現噪聲。本研究在分析國內某食用油主要產出省2017—2019年11 345 行花生油安全日常檢測數據及其相關信息的基礎上，建立了基于優化的GRA-ELM風險評價模型。首先對指標進行篩選，得到6 類風險指標，再對數據進行預處理，得到150 組樣本的花生油檢測數據。將預處理后數據輸入到小波閾值濾波模塊中，對濾波后數據進行GRA，得到指標的權重，并與指標相對風險值結合得到綜合風險評價期望值，接著輸入到ELM網絡中訓練，再利用PBO算法進行參數優化。最后對評價模型輸出的綜合風險值進行模糊綜合分析，實現對花生油樣品風險程度分級。在不同模型對比實驗中，首先通過與BP、RBF網絡模型的對比，ELM網絡的R2和RMSE分別為0.132 0和0.796 4，性能優于BP網絡和RBF網絡，證明ELM網絡擁有更好的泛化能力和魯棒性；其次通過對比有無濾波模塊以及有無PBO參數優化情況下的結果，最終得出加上濾波模塊且進行PBO參數優化后的WT-PBOGRA-ELM模型R2和RMSE分別為0.056 3和0.946 1，證明了濾波模塊和PBO算法的有效性；最后通過對比模型評價花生油樣本2的結果，得出只有WT-PBO-GRA-ELM模型得到的風險等級結果為極高風險，與該樣本真實風險等級一致，證明了本研究風險評價模型的有效性。

未來本研究還可以在以下幾個方面更加深入，首先在數據方面，可以將地理因素、時序因素加入模型從而提取出相關風險規律。其次在指標賦權方面，可以將專家打分等主觀賦權與危害物本身的毒理學特性和基于數據的客觀賦權有機地結合在一起形成更加綜合合理的風險因子權重體系。此外，還可以進一步研究對于整條食用油加工鏈的風險評價。