集成模糊層級劃分的LightGBM食品安全風(fēng)險預(yù)警模型：以肉制品為例

高亞男，王文倩，王建新

（北京林業(yè)大學(xué)信息學(xué)院，北京 100083）

食品安全問題是關(guān)系公眾生命健康、影響經(jīng)濟發(fā)展和人類社會穩(wěn)定的全球議題。近年來，國內(nèi)外由食品安全問題導(dǎo)致的公共健康事故屢見不鮮，如南韓的“鉛螃蟹”和“寄生蟲泡菜”[1]、歐洲的“馬牛肉”[2]、國內(nèi)的“塑化劑保健品”和“硫磺生姜”[3]等。為提高食品行業(yè)的安全監(jiān)督監(jiān)控水平，世界各國均頒布了一系列法律文件用于支持食品行業(yè)風(fēng)險控制。歐盟通過面向市場和面向食品安全兩個階段完成了《歐盟食品法》的頒布，進一步建立了歐洲食品安全局（European Food Safety Authority，EFSA）來全方位覆蓋食品供應(yīng)鏈的風(fēng)險評估[4]。在我國，2009年頒布的《食品安全法》取代了《食品衛(wèi)生法》，更加突出強調(diào)食品風(fēng)險評估的必要性[5]。

在國家市場監(jiān)督管理總局的指導(dǎo)下，各地食品藥品監(jiān)督管理機構(gòu)積極進行監(jiān)管技術(shù)創(chuàng)新和技術(shù)難點攻克，不斷加強食品監(jiān)管信息化建設(shè)。然而，相較其他國家，悠久的飲食文化和特有的地理位置造就了我國現(xiàn)有的復(fù)雜食品體系，在食品種類、加工方法、包裝貯存方法、添加劑與添加方式等方面的多樣性遠超國外[6]。對此，本研究在綜合分析數(shù)據(jù)特點后，使用界限值劃分層級結(jié)構(gòu)和貝葉斯先驗概率結(jié)合的方法計算模型期望輸出，并應(yīng)用訓(xùn)練成本低、準確率高的LightGBM進行風(fēng)險值預(yù)測。在實際應(yīng)用中，專家可以對預(yù)測結(jié)果進行校正，幫助提高產(chǎn)出規(guī)則的準確性，直至預(yù)測模型逐漸逼近最優(yōu)決策方案。

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

1.1 食品安全數(shù)據(jù)來源與特點

食品安全數(shù)據(jù)的來源大致可以概括為3 個方面：1）靜態(tài)數(shù)據(jù)；2）動態(tài)數(shù)據(jù)；3）專家經(jīng)驗數(shù)據(jù)。其中，靜態(tài)數(shù)據(jù)指一經(jīng)定義在一段時間內(nèi)不改變的數(shù)據(jù)，如法令法規(guī)中的標準數(shù)據(jù)（合格標準線、最低最高檢出限、檢驗依據(jù)、判定依據(jù)等）、銷售生產(chǎn)企業(yè)的信息數(shù)據(jù)（企業(yè)規(guī)模、成立年限、主要食品類別、生產(chǎn)銷售區(qū)域、采購地點等）等。該類數(shù)據(jù)大多存在于各地工商管理局的企業(yè)注冊信息記錄文件中。而每年由企業(yè)食品安全體系檢測（生成、加工、流通、消費等環(huán)節(jié)的例行檢查記錄和違法違規(guī)記錄等）和常規(guī)食品抽檢產(chǎn)生的檢查數(shù)據(jù)和抽檢數(shù)據(jù)（食品類別、檢驗項目、檢出含量、生產(chǎn)時間、抽檢時間等）受時間、地域等因素的影響，檢查結(jié)果和抽檢記錄信息保持動態(tài)積累和更新。數(shù)據(jù)通過各地食品安全監(jiān)督管理系統(tǒng)采集錄入并存儲在相應(yīng)數(shù)據(jù)庫中以供檢索和查詢。相較于前兩種，專家依據(jù)經(jīng)驗或文獻調(diào)研給出的數(shù)據(jù)，如某類食品的污染物出現(xiàn)概率、污染指標度量值、項目風(fēng)險等級值等，既可以在一段時間內(nèi)穩(wěn)定不變，也可以隨不同應(yīng)用背景而動態(tài)調(diào)整。

食品安全相關(guān)數(shù)據(jù)涉及食品信息、生產(chǎn)企業(yè)信息、銷售企業(yè)信息、檢驗機構(gòu)信息、國家法令法規(guī)、專家經(jīng)驗指標、各污染物檢驗標準等方面[7]。其屬性類型大體可以概括為離散字符型屬性（企業(yè)規(guī)模、生產(chǎn)省份、檢驗項目等）、離散數(shù)值型屬性（生產(chǎn)日期、抽樣日期、樣品狀態(tài)等）、連續(xù)數(shù)值型屬性（企業(yè)營業(yè)額、檢出含量等）。屬性類別多、屬性值格式雜亂使數(shù)據(jù)具有高維度和高復(fù)雜度的特點，因而數(shù)據(jù)多呈現(xiàn)線性不可分狀態(tài)，且分布規(guī)律隱藏在不完整數(shù)據(jù)、缺失數(shù)據(jù)、錯誤錄入數(shù)據(jù)等噪聲干擾中，這更增加了風(fēng)險分析的難度。

1.2 食品安全風(fēng)險評估的國內(nèi)外研究進展

食品安全事件造成的惡劣影響以及食品安全數(shù)據(jù)的特點促使相應(yīng)的國家法律法規(guī)頒布，對食品安全風(fēng)險預(yù)警的需求日益提高。有效的風(fēng)險預(yù)警模型能通過對先驗知識的挖掘得出既定規(guī)律，并分析出風(fēng)險要素、風(fēng)險程度或預(yù)測出風(fēng)險值，對政府合理分配有限資源、正確定位風(fēng)險點、從源頭解決食品安全問題具有重要意義[8]。國內(nèi)外學(xué)者為此展開了大量研究和應(yīng)用實踐。Delphi法是一種用于融合來自不同地區(qū)不同領(lǐng)域眾多專家意見的方法，包括重復(fù)多輪意見反饋和根據(jù)中間反饋結(jié)果修改后續(xù)調(diào)查問卷兩個步驟[9]。依據(jù)安全事項重要性劃分而進行三輪專家經(jīng)驗反饋的Delphi法被用于巴西食品交易風(fēng)險評估工具中。最終的評估工具由39 個風(fēng)險事項組成，其中包含專家在原有因素列表上的增加事項和細化事項[10]。但從輸入屬性指標到結(jié)果輸出指標，該方法中的全部評價指標均由專家人為確定，評估效率較低且人工成本高，不能滿足風(fēng)險預(yù)警系統(tǒng)時效性的需求。層次分析法（analytic hierarchy process，AHP）通過建立目標層、準則層、方案層的三層結(jié)構(gòu)來處理復(fù)雜的多目標決策問題，因此可用來計算屬性指標的權(quán)重和劃分評價等級，進而具有屬性約減和指標去噪的功能[11]。Rathore等[12]使用單一的AHP方法建立了印度食品供應(yīng)鏈的定量風(fēng)險評估模型，通過輸入供應(yīng)鏈專家的指標偏好初始化比較矩陣，其結(jié)果是輸出方案指標權(quán)重，以此來查找薄弱環(huán)節(jié)。Delphi法和AHP法是目前比較成熟的風(fēng)險評估方法。除此之外，風(fēng)險矩陣[13]和灰色理論[14]等方法也被融合應(yīng)用到食品安全風(fēng)險評估領(lǐng)域中并取得一定進展。然而，上述方法存在屬性覆蓋范圍小、指標精度低、人為成本高、評估過程長、無法動態(tài)更新、適應(yīng)性弱等不足，造成評估結(jié)果的準確率低且缺失精準定位風(fēng)險點的能力。

計算機軟硬件的高速發(fā)展正在推動食品行業(yè)的信息化建設(shè)進程，食品安全相關(guān)數(shù)據(jù)的積累給智能計算方法在食品安全風(fēng)險評估領(lǐng)域的應(yīng)用創(chuàng)造了條件。章德賓等[15]將167 個食品檢驗項目含量作為輸入，5 個污染大類的預(yù)測值作為輸出，構(gòu)建具有兩層隱含層的BP神經(jīng)網(wǎng)絡(luò)，按照預(yù)測值大小預(yù)測出下一周的污染主要來源為“致病菌超標”和“獸藥殘留”。同樣運用BP神經(jīng)網(wǎng)絡(luò)模型，王星云等[16]選定重金屬“鉛”抽檢記錄中的“生產(chǎn)企業(yè)省份”、“抽樣地點”等13 個屬性作為輸入，預(yù)測記錄的“檢驗結(jié)果”。在這里，“檢驗結(jié)果”屬性的值是合格與不合格，是抽檢記錄的既有屬性，不需要人為制定標簽，這提高了預(yù)測結(jié)果的準確性。除此之外，Wang Jing等[17]將某乳品企業(yè)的產(chǎn)品作為分析對象，提取運輸時間、溫度、季節(jié)、包裝方式等7 個影響因素作為規(guī)則挖掘?qū)傩皂?，利用關(guān)聯(lián)規(guī)則Apriori算法產(chǎn)生規(guī)則庫，并使用支持度和置信度過濾規(guī)則保留出現(xiàn)最頻繁的規(guī)則組合，該結(jié)果被作為供應(yīng)鏈環(huán)節(jié)制定中應(yīng)盡量避免產(chǎn)生的流程組合。陳夢音等[18]先運用主成分分析篩選牦牛奶渣質(zhì)量的評價指標，隨后在保留外觀色澤和營養(yǎng)品質(zhì)兩個主要指標的基礎(chǔ)上使用聚類做數(shù)據(jù)分簇，最后用P值評價簇差異性，得出簇內(nèi)規(guī)律。Han Yongming等[19]使用灰色關(guān)聯(lián)分析法確定抽檢數(shù)據(jù)的指標權(quán)重后，制定標簽風(fēng)險值，再運用隱馬爾可夫法做模型訓(xùn)練和風(fēng)險值預(yù)測，最后從時間維度分析得出的準確率變化與食品質(zhì)量變化趨勢相符。AHP法作為復(fù)雜度低的指標權(quán)重計算和風(fēng)險分級的重要方法，通常與其他理論方法結(jié)合來提高風(fēng)險評估的準確性和可解釋性。Su Xiaoyan等[20]首先用AHP計算由專家預(yù)定義的風(fēng)險指標的權(quán)重，再結(jié)合Dempster-Shafer理論合成最終風(fēng)險值。相似地，在通過AHP計算各風(fēng)險指標的權(quán)重制定風(fēng)險值標簽后，Geng Zhiqiang等分別應(yīng)用極限學(xué)習(xí)機[21]和深度徑向基神經(jīng)網(wǎng)絡(luò)[22]，將各指標值作為輸入來預(yù)測風(fēng)險值大小。除此之外，支持向量機[23]也被用于少量抽檢數(shù)據(jù)的風(fēng)險評估工作中。圖1概括了上述提到的風(fēng)險預(yù)警方法及其大體算法流程。

先前的研究大多利用抽檢數(shù)據(jù)中的數(shù)值型特征作為輸入進行風(fēng)險值預(yù)測，如各檢驗項目含量、生產(chǎn)抽樣季度、溫度、產(chǎn)量等。特征數(shù)量少會使風(fēng)險預(yù)測值與實際風(fēng)險值的誤差增大，造成與事實相反的“假風(fēng)險”，甚至誤導(dǎo)決策。由于不能全面考慮風(fēng)險點屬性，使得有價值的數(shù)據(jù)蘊含的關(guān)系和規(guī)律被忽略。除此之外，在訓(xùn)練數(shù)據(jù)標簽的制定上，只考慮了專家因素對指標權(quán)重的影響而忽略了大量歷史數(shù)據(jù)中已存在的顯著統(tǒng)計量信息，沒有校驗二者誤差的方法而直接制定標簽，不可避免地存在錯誤標簽。為了保證實驗數(shù)據(jù)的可信賴性和降低錯誤預(yù)警的風(fēng)險，用歷史數(shù)據(jù)的先驗風(fēng)險概率結(jié)合模糊層級劃分法計算各屬性的具體特征值權(quán)重，將離散型和連續(xù)型屬性特征值加權(quán)求和再進行歸一化的結(jié)果作為風(fēng)險值。為了學(xué)習(xí)到上述既定規(guī)則，使用獨熱編碼處理離散型屬性，用LightGBM模型用于結(jié)果預(yù)測，并隨之融合了專家干預(yù)策略進行準確性驗證和結(jié)果校正。

圖1 相關(guān)算法研究分類Fig.1 Classification of the previous algorithm studies

1.3 提出的風(fēng)險預(yù)警模型

圖2 算法總體流程Fig.2 Flow chart of the algorithm proposed in this paper

如圖2所示，原始數(shù)據(jù)被分成兩步處理：一方面用于模型期望輸出值計算；另一方面用于模型輸出特征處理。兩步的處理結(jié)果特征被劃分為訓(xùn)練集和測試集，分別用于模型訓(xùn)練和測試驗證結(jié)果。

1.3.1 模糊綜合風(fēng)險值計算

1.3.1.1 先驗風(fēng)險統(tǒng)計量

包含有限屬性值（字符型、數(shù)值型等）個數(shù)的屬性稱為離散型屬性，如食品類別、生產(chǎn)省份、經(jīng)濟地帶等。假設(shè)屬性D包含nD個屬性值，表示為D＝{ai}（i＝1, 2, 3,...,nD），則采用歷史數(shù)據(jù)先驗概率歸一化后的結(jié)果值定義離散屬性值權(quán)重，如公式（1）所示。

式中：B表示抽樣批次數(shù)；R表示記錄數(shù)。B（a,b）、R（a,b）分別表示滿足a和b條件下的抽樣批次數(shù)和記錄數(shù)，如R（False,ai）表示在所有不合格記錄中，屬性值ai所占記錄數(shù)。B（False,:）、B（Total,:）分別表示不合格記錄數(shù)和總記錄數(shù)。這樣可以得到每一個離散屬性的值和權(quán)重的鍵值對應(yīng)關(guān)系。

1.3.1.2 模糊層級劃分

可由兩種及兩種以上其他類型屬性構(gòu)造生成的屬性為間接屬性。如，由“生產(chǎn)日期”、“抽樣日期”和“保質(zhì)期”可構(gòu)造得到的“保質(zhì)周期剩余率”屬性表示食品保質(zhì)期內(nèi)剩余時間占保質(zhì)期的比例、由“生產(chǎn)地址”和“抽樣地址”可構(gòu)造得到“運輸距離”屬性表示食品的近似運輸距離等。除此之外，部分連續(xù)數(shù)值型屬性具有相關(guān)法規(guī)文件規(guī)定的限值作為層級劃分參考。而剩余的連續(xù)數(shù)值型屬性只有大量數(shù)值型數(shù)據(jù)特征。

圖3 一般連續(xù)數(shù)值型屬性值分布規(guī)律（以“保質(zhì)周期剩余率”為例）Fig.3 Distribution of generally continuous attributes (taking the quality-cycle residual rate as an example)

連續(xù)數(shù)值型屬性的權(quán)重計算過程由層級劃分和模糊隸屬度計算兩步構(gòu)成。層級劃分中，針對有限值參考的屬性，根據(jù)限值制定劃分界限。如各檢驗項目含量具有“最小檢出限”和“最大允許限”兩個限值參考。其他的連續(xù)數(shù)值型屬性，首先利用屬性數(shù)據(jù)構(gòu)造擬合近似正態(tài)分布，得到分布參數(shù)，然后根據(jù)正態(tài)分布的3σ原則劃分層級界限。以“保質(zhì)周期剩余率”屬性為例，圖3A表示其原始數(shù)據(jù)分布規(guī)律，可以看出整體符合保質(zhì)周期剩余率與分布頻率成正比的“長尾”分布。即大部分抽檢食品處于生產(chǎn)完成初期、食品質(zhì)量優(yōu)質(zhì)階段，極少部分抽檢食品處于即將過期階段。通過圖3B可以看出，保質(zhì)周期剩余率通過分布反轉(zhuǎn)、平移、取對數(shù)構(gòu)造后的頻率分布可近似擬合為正態(tài)分布。雖然正態(tài)擬合存在誤差，但該方法通過均值μ和標準差σ劃定的分布界限仍可消除極端異常點的影響，使得屬性值層級分布符合雙向厚長尾的特點。

圖4 連續(xù)數(shù)值型屬性模糊層級劃分Fig.4 Hierarchy partition of continuous attributes

有界限值的連續(xù)數(shù)值型屬性采用如圖4所示的劃分方法。其中，min表示檢驗項目含量屬性的“最小檢出限”，max表示檢驗項目含量屬性的“最大允許限”，屬性值被劃分為4 個層次等級。在此基礎(chǔ)上，公式（2）被用于計算污染物屬性p在第i條抽檢記錄中的含量（xpi）屬于第j等級的程度（Cpj），來確定xpi的層級。其中，lpj表示污染物屬性p在層次等級j上的分界值，由此得到含量xpi對應(yīng)的層級Jpi的隸屬度CpJ＝max（Cpj）（j＝1, 2, 3,...,C）。依據(jù)界限值和分布擬合參數(shù)制定的層級劃分界限，可以將連續(xù)數(shù)值型屬性離散化為層級分布，完成連續(xù)屬性向離散屬性的轉(zhuǎn)化，進而根據(jù)公式（1）計算各層級先驗概率權(quán)重。

1.3.1.3 綜合風(fēng)險值計算

設(shè)有n條食品相關(guān)數(shù)據(jù)記錄，每條記錄共包含m個屬性（包含m1個離散屬性和m2個連續(xù)數(shù)值型屬性），數(shù)據(jù)集表示為則由上述權(quán)重計算過程可得到離散屬性的權(quán)重矩陣為其被作為離散屬性值對風(fēng)險值的貢獻程度，即類似的連續(xù)數(shù)值型屬性層級對應(yīng)的權(quán)重矩陣為設(shè)第j個連續(xù)型屬性值序列為＝{xj1，xj2，xj3，...，xjm2}T，1≤j≤m2，則其對風(fēng)險值的貢獻程度由公式（3）定義。

自然對數(shù)化處理后，第i條抽檢記錄的模糊綜合風(fēng)險值通過公式（6）計算得到。至此模型的期望輸出標簽計算完成。

1.3.2 基于One-Hot編碼的特征重構(gòu)

One-Hot編碼是對有限個狀態(tài)值進行處理的編碼方式，每個狀態(tài)值對應(yīng)獨立的寄存器位，并且只有一位有效，被廣泛應(yīng)用于離散特征處理中[24]。與通常定義的處理方式不同，對含有M個狀態(tài)值的離散屬性特征采用（M－1）個寄存器位編碼，所有的寄存器位對第M個狀態(tài)值無效。因此，每個狀態(tài)值映射到（M－1）維的特征集合為S＝{s1,s2,...,sM－1,sM}，如圖5所示。其中第j個特征值的第i個特征位取值表示為δ（si,sj），i≥1，j≤M－1，δ為Kronecker符號函數(shù)。并且第M個值的狀態(tài)位都為無效。如抽樣地區(qū)的所屬“經(jīng)濟地帶”屬性包含“東北”、“東部”、“西部”、“中部”4 個屬性值，按照所定義的One-Hot編碼方式可將其映射到三維特征空間，編碼為“100”、“010”、“001”和“000”。

圖5 離散屬性的One-Hot編碼Fig.5 One-Hot coding of discrete attributes

1.3.3 預(yù)測模型

1.3.3.1 梯度提升樹

在分類問題中，決策樹以Gini系數(shù)或信息增益作為指標衡量屬性特征對類別的區(qū)分程度，從而用來確定分割節(jié)點。處理回歸時，Cart回歸樹采用均方誤差或指數(shù)誤差等作為指標選擇最佳分割點。在有監(jiān)督學(xué)習(xí)中，期望得到一個準確率高、泛化力強且穩(wěn)定的模型，但由于數(shù)據(jù)量和方法自身缺陷限制，往往只能得到多個有偏模型，即弱監(jiān)督模型，也可稱為子模型。集成學(xué)習(xí)通過投票組合這些弱監(jiān)督模型，以平滑噪聲糾正誤差，得到表現(xiàn)更優(yōu)的強監(jiān)督模型，即“bagging”的思想[25]。以決策樹為子模型，隨機森林通過兩次隨機選取完成模型生成：隨機選取訓(xùn)練集和隨機選取子模型分裂特征。其中，有放回的隨機抽樣保證了子模型訓(xùn)練集的獨立同分布，選擇部分特征實現(xiàn)分裂有助于防止過擬合?？梢姡跇?gòu)造隨機森林過程中，各弱監(jiān)督模型生成過程中，對數(shù)據(jù)有相同的采樣優(yōu)先級。

梯度提升樹（gradient boosting decision tree，GBDT）同樣以Cart回歸樹作為弱學(xué)習(xí)器，但與RF的子模型的并行構(gòu)造不同，GBDT采用子模型關(guān)聯(lián)遞進構(gòu)造的“boosting”思想，各子模型之間通過擬合負梯度，解決了損失函數(shù)為一般函數(shù)的問題，以期達到損失最速下降的目的[26]。設(shè)有N條訓(xùn)練數(shù)據(jù)。設(shè)迭代次數(shù)為M，即子模型生成個數(shù)為M。第m次迭代的損失函數(shù)為L（y,fm－1（x）+Tm（x））），其中，fm－1（x）表示第m－1次迭代生成的強學(xué)習(xí)器，Tm（x）表示第m個子模型的輸出，x表示輸入特征，y表示對應(yīng)的期望輸出。其算法步驟如下：

初始化弱學(xué)習(xí)器為常數(shù)，如公式（7）所示。

對第m（m＝1, 2,...,M）次迭代有：

計算每個樣本i＝1, 2,...,N在當(dāng)前第m－1個強學(xué)習(xí)器上的負梯度為對第i個樣本，將gim代替yi作為期望輸出，組成新的訓(xùn)練樣本集{（xi,gim）}（i＝1, 2,...,N），學(xué)習(xí)得到第m棵回歸樹Tm（x）。設(shè)Tm（x）有J個葉子節(jié)點，即所有訓(xùn)練數(shù)據(jù)被分割成J塊區(qū)域，第j塊表示為Rjm（j＝1, 2,...,J）。

對所有葉子節(jié)點區(qū)域計算最優(yōu)擬合值，如式（8）所示：

更新強學(xué)習(xí)器，如式（9）所示：

得到最終GBDT模型，如式（10）所示：

可以看出，第m個強學(xué)習(xí)器擬合的是損失對f（x）在前（m－1）個弱學(xué)習(xí)器累加處的一階泰勒展開。因此，這樣的函數(shù)估計存在誤差?；诖?，Chen Tianqi等[27]提出了XGBoost模型，將其擴展成二階泰勒展開，并加入正則項修正，大大提高了模型的精確度。

1.3.3.2 LightGBM原理與優(yōu)勢分析

上述提到的GBDT類模型，在函數(shù)空間中使用梯度下降方法調(diào)優(yōu)來得到最佳模型。然而在構(gòu)造子模型時，對每個特征都要遍歷所有樣本點來選取最優(yōu)分割點，該操作非常耗時。對此，Ke Guolin等[28]提出的LightGBM模型分別在訓(xùn)練數(shù)據(jù)上采用基于梯度的單邊采樣和在特征上采用互斥特征捆綁來提高學(xué)習(xí)器的訓(xùn)練速度和泛化能力。每次更新弱學(xué)習(xí)器時，單邊采樣在不改變特征值分布和不損失精度的同時，壓縮訓(xùn)練數(shù)據(jù)集，減少計算量。此外，采樣增加了弱學(xué)習(xí)器的多樣性，進而提高了泛化能力。然而，高維數(shù)據(jù)的特征往往是稀疏的。這就使得數(shù)據(jù)中存在大量的互斥特征，即對一條記錄來說，通常不會同時出現(xiàn)非0值，如經(jīng)過One-Hot編碼處理后的特征。因此，LightGBM對互斥特征進行捆綁來提高運行效率。其中，涉及兩個問題：哪些特征應(yīng)該被綁在一起和綁定的方式。通過構(gòu)建加權(quán)無向圖，LightGBM將構(gòu)建特征集合建模成圖著色問題，使用類似貪心算法求取結(jié)果，復(fù)雜度為O（#feature2）。雖然在處理高維特征時復(fù)雜度較高，但只需要處理一次。關(guān)于捆綁方式，算法使用劃分直方圖的方法捆綁特征集合中的互斥特征。通過記錄每個特征塊的個數(shù)，平移界定范圍來合并互斥特征。通過這種方式，弱學(xué)習(xí)器查找分割節(jié)點時只需遍歷所有塊值來確定分割點即可，遍歷復(fù)雜度為O（#塊數(shù)*feature2）。雖然特征值離散化后找到的分割點不如原始特征值精確，但決策樹本身為弱學(xué)習(xí)器，bin分割帶來的誤差剛好起到正則項的作用。

1.3.4 模型應(yīng)用的具體流程

圖6 專家干預(yù)流程Fig.6 Flow chart of experts' modification

圖7 風(fēng)險預(yù)警系統(tǒng)展示Fig.7 Display of food safety risk prewarning system

在實際應(yīng)用中，初始風(fēng)險值數(shù)據(jù)需要根據(jù)1.3.1節(jié)的計算規(guī)則進行預(yù)先計算，并用于最初的模型生成。但由于計算規(guī)則的不易更改和表示能力有限的特殊性，在產(chǎn)生初始模型的基礎(chǔ)上，引入專家干預(yù)策略來輔助重新訓(xùn)練更加可靠的模型。如圖6所示，首先對現(xiàn)有風(fēng)險值數(shù)據(jù)進行風(fēng)險分析，得出簡單的統(tǒng)計規(guī)律。若分析結(jié)果與專家經(jīng)驗存在明顯差異，即產(chǎn)生了異常風(fēng)險，專家可通過修改異常風(fēng)險數(shù)據(jù)對應(yīng)的屬性值權(quán)重來重新計算風(fēng)險值。修改完成后，加載參與風(fēng)險分析的現(xiàn)有數(shù)據(jù)和修改后的數(shù)據(jù)，重新訓(xùn)練下一代模型。當(dāng)有新數(shù)據(jù)錄入時，使用最新版本模型進行預(yù)測，隨后進行風(fēng)險分析。循環(huán)往復(fù)，直至沒有異常風(fēng)險產(chǎn)生。此時，相較初始規(guī)則，新版本模型表示的規(guī)則雖不能被公式量化，但風(fēng)險預(yù)測更加精確，且數(shù)據(jù)魯棒性更高。針對上述的方法流程，開發(fā)了相應(yīng)的食品安全風(fēng)險分級預(yù)警系統(tǒng)，用于功能可視化操作。部分系統(tǒng)界面如圖7所示。

2 結(jié)果與分析

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

以2015—2017年上半年的肉制品相關(guān)數(shù)據(jù)為例，從中抽取部分屬性用于模型優(yōu)化和結(jié)果評估。對數(shù)據(jù)的預(yù)處理過程可概括為如下步驟。

2.1.1 屬性選取

選取特征值數(shù)量介于“食品大類”和“食品細類”之間的“食品次亞類”作為代表食品類別信息的屬性。同理，選取剩余屬性中具有代表性的屬性：“包裝分類”、“區(qū)域類型”、“抽樣地點”、“抽樣環(huán)節(jié)”、“抽樣省份”、“生產(chǎn)省份”、“經(jīng)濟地帶”、“樣品狀態(tài)”、“樣品類型”、“生產(chǎn)企業(yè)年銷售額”。除此之外，選取原始數(shù)據(jù)中界限值明顯、含量值有效記錄數(shù)量大的檢驗項目作為輸入，包括微生物類（“大腸菌群”、“菌落總數(shù)”）、防腐劑類（“山梨酸”）、重金屬類（“鉛”、“鉻”、“鎘”）和其他類（“過氧化值”、“亞硝酸鹽”）。

2.1.2 屬性特征值規(guī)范

對檢驗項目名稱做規(guī)范合并，如“山梨酸及其鉀鹽”統(tǒng)一為“山梨酸”、“亞硝酸鹽（以亞硝酸鈉計）”和“亞硝酸鹽殘留量”統(tǒng)一為“亞硝酸鹽”以及“鉛（以Pb計）”統(tǒng)一為“鉛”等。由于高風(fēng)險的檢出必要性，在檢驗項目含量中做如下規(guī)范：將“未檢出”規(guī)范為當(dāng)前記錄的“標準最低檢出限”；將＜0.005等類似含量規(guī)范為＜符號右邊數(shù)值。部分其他原始數(shù)據(jù)中離散型屬性的特征值規(guī)范結(jié)果分布如表1所示。

表1 離散型屬性數(shù)據(jù)特征值分布Table 1 Feature distribution of discrete attributes

2.1.3 屬性轉(zhuǎn)換和生成

將“抽樣日期”、“生產(chǎn)日期”兩個時間戳屬性分別轉(zhuǎn)換為從{1, 2, 3, 4}中取值的“抽樣季度”和“生產(chǎn)季度”的季度信息；將檢驗項目的檢出含量轉(zhuǎn)換為同一單位度量；保質(zhì)期屬性值中，如“2 年”、“常溫下9 個月”、“2～8 ℃下6 個月”、“25 ℃下3 個月”，根據(jù)抽檢季度特征轉(zhuǎn)換為對應(yīng)時間/d。除此之外，利用公式（11）可計算得到“保質(zhì)周期剩余率”屬性。

表2 處理后的數(shù)據(jù)分布Table 2 Distribution of processed data

通過上述處理，結(jié)合1.3.1節(jié)中風(fēng)險值標簽的計算方法，可得到如表2所示的實驗數(shù)據(jù)。數(shù)據(jù)包含10 個連續(xù)數(shù)值型屬性和11 個離散型屬性，其中離散型屬性又分為9 個離散字符型屬性和3 個離散數(shù)值型屬性。所有數(shù)據(jù)共計9 757 批次。

2.2 模型對比實驗

將所有數(shù)據(jù)按0.2訓(xùn)練測試比分為包含7 806 批次數(shù)據(jù)的訓(xùn)練集和1 951 批次數(shù)據(jù)的測試集。訓(xùn)練集仍按0.2驗證比劃分為6 245 批次數(shù)據(jù)用于訓(xùn)練模型和1 561 批次驗證數(shù)據(jù)用于模型超參調(diào)優(yōu)。在得到具有較優(yōu)參數(shù)模型的基礎(chǔ)上，測試集用于測試不同模型的預(yù)測誤差，并從中選取最適方法模型。本實驗環(huán)境是i7-8550U CPU，四核16G RAM的Windows 10操作系統(tǒng)。實驗調(diào)用了LightGBM官網(wǎng)提供的LGBMRegressor方法[29]做回歸預(yù)測。

2.2.1 LightGBM參數(shù)調(diào)優(yōu)結(jié)果

圖8 LightGBM參數(shù)調(diào)整Fig.8 Parameter adjustments of LightGBM

由1.3.3.2節(jié)可見，LightGBM模型有大量與訓(xùn)練速率、收斂速率、準確率和擬合程度相關(guān)的參數(shù)待調(diào)整。暫且只選取直接影響準確率和模型泛化能力的部分參數(shù)進行調(diào)整，包括：“l(fā)earning_rate”、“max_bin”、“max_depth”、“num_leaves”和“n_estimators”。實驗選擇均方誤差（mean square error，MSE）作為結(jié)果評價指標。首先，根據(jù)數(shù)據(jù)集大小，由經(jīng)驗給出初始參數(shù)的大體范圍，然后固定其他參數(shù)只調(diào)整“num_leaves”，選取最優(yōu)模型對應(yīng)的參數(shù)，固定當(dāng)前參數(shù)進而調(diào)整下一個參數(shù)，循環(huán)運行直至損失符合要求。圖8展示了模型參數(shù)的調(diào)整流程及產(chǎn)出結(jié)果?？梢钥闯?，隨著參數(shù)的優(yōu)化，模型的驗證損失由0.052 7下降到0.052 4，因此選擇learning_rate=0.01、max_bin=200、max_depth=12、num_leaves=21、n_estimators=900對應(yīng)的LightGBM模型作為最優(yōu)模型。

2.2.2 優(yōu)化后的模型結(jié)果

圖9是模型以選定的最優(yōu)參數(shù)訓(xùn)練生成的第900個弱學(xué)習(xí)器的一部分分支。從圖中第一層的最左側(cè)分支可以看出，“年銷售額”屬性通過歸一化后的特征閾值－0.095 3和分割的信息增益值0.259 5可以將6 552 條數(shù)據(jù)記錄分割成包含預(yù)測標簽值≈sigmod（－0.006 1）＝0.498 5的4 954 條記錄的“大腸菌群”節(jié)點和預(yù)測標簽值≈sigmod（0.008 5）＝0.502 1的1 598 條記錄的“菌落總數(shù)”節(jié)點，分割效果明顯。可見，“年銷售額”屬性對預(yù)測結(jié)果的分割能力較其他屬性更強，即在確定風(fēng)險值的過程中，“年銷售額”屬性起到了重要作用。同理，可以得出“菌落總數(shù)”、“大腸菌群”等屬性對風(fēng)險的貢獻程度也較高。

圖9 第900個弱學(xué)習(xí)器Fig.9 The 900th weak learning classifier

圖10 訓(xùn)練過程中的屬性重要程度排序Fig.10 Feature importance rank during training process

上述結(jié)果與圖10所示的訓(xùn)練過程中模型輸出的屬性重要度排序基本吻合。不能完全吻合的原因是，圖9只是900 個弱學(xué)習(xí)器中的其中一個，而圖10則反映了所有弱學(xué)習(xí)器中屬性作為分割節(jié)點的總次數(shù)排序，即屬性的重要程度排序?？梢钥闯?，生產(chǎn)企業(yè)的“年銷售額”在訓(xùn)練過程中作為弱學(xué)習(xí)器的分割點次數(shù)為4 662 次，超過其他屬性的分割節(jié)點次數(shù)。這說明其相對其他屬性對預(yù)測風(fēng)險值作出的貢獻最多，對當(dāng)前標簽數(shù)據(jù)具有最好的劃分能力。除此之外，通過計算連續(xù)型屬性值與風(fēng)險值的Pearson系數(shù)和比較離散型屬性O(shè)ne_Hot編碼對應(yīng)的平均風(fēng)險值大小，可以得出正負相關(guān)性。由于年銷售額與企業(yè)規(guī)模和企業(yè)成熟度有非常強的正相關(guān)關(guān)系，因此這個結(jié)果與領(lǐng)域?qū)＜业姆治龊皖A(yù)期相符，即規(guī)模較大的企業(yè)風(fēng)險程度較低。其次，對肉制品風(fēng)險影響較大的因素還有微生物類檢驗項目（“菌落總數(shù)”、“大腸菌群”）、重金屬類檢驗項目（“鉛”、“鉻”）、防腐劑類檢驗項目（“山梨酸”）、其他類（“亞硝酸鹽”）和保質(zhì)期相關(guān)屬性（“保質(zhì)周期剩余率”）。還可以得出，在食品種類中，“熟肉干制品”類食品的風(fēng)險程度較高，在“貴州省”抽樣和在第一季度生產(chǎn)的肉制品的風(fēng)險程度較低，“西部”地區(qū)相較其他地區(qū)的肉制品風(fēng)險程度較高等。

2.2.3 各模型預(yù)測結(jié)果對比

為了驗證LightGBM的有效性，用相同的訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)進行不同方法之間的結(jié)果對比。同樣地，仍將訓(xùn)練數(shù)據(jù)按照20%分割得到驗證集調(diào)整其他方法的超參。對比方法包括BP神經(jīng)網(wǎng)絡(luò)、RBF網(wǎng)絡(luò)、隨機森林、通常的GBDT、XGBoost和LightGBM模型，表3列出了在當(dāng)前數(shù)據(jù)集規(guī)模的前提下，各模型調(diào)整的參數(shù)結(jié)果、調(diào)參過程中的最小驗證集損失和在同一測試集上的預(yù)測誤差。除此之外的其他參數(shù)，均使用默認值。

表3 對比方法參數(shù)設(shè)置Table 3 Parameters of comparative methods

2.3 水產(chǎn)品的風(fēng)險分析應(yīng)用

2.3.1 數(shù)據(jù)介紹

以2015—2017年上半年的水產(chǎn)品相關(guān)數(shù)據(jù)為測試用例。由于該數(shù)據(jù)在“年銷售額”、“保質(zhì)期”屬性上的缺失值和錯誤值過多，如果按其進行數(shù)據(jù)過濾，數(shù)據(jù)量將大大降低，這必然會引起實驗結(jié)果的偏差。因此，為了驗證模型的通用性，得到更精確的規(guī)律，在水產(chǎn)品數(shù)據(jù)中增加了“恩諾沙星”、“孔雀石綠”、“存儲狀態(tài)”等屬性。最終，數(shù)據(jù)被擴充為包含39 個連續(xù)型屬性和13 個離散型屬性，共52 維屬性字段，287 430 條記錄，計7 370 批次。

2.3.2 風(fēng)險分析

假設(shè)不計算風(fēng)險值而直接依據(jù)“是否合格”來判定水產(chǎn)品污染物中的風(fēng)險分布，如圖11所示（比例代表當(dāng)前污染物導(dǎo)致的不合格產(chǎn)品記錄在總不合格產(chǎn)品記錄中的占比）?？梢钥闯觯亟饘俸汪~藥類污染物是導(dǎo)致水產(chǎn)品不合格的主要原因。其中，重金屬中的“鎘”問題最為突出，占比24.6%。魚藥類中問題嚴重的污染物主要包括“孔雀石綠”、“恩諾沙星”、呋喃類（“呋喃西林”、“呋喃唑酮”）和抗生素類（“氯霉素”）。雖然通過這種方式可以大致看出污染物中的風(fēng)險分布情況，但對合格程度和不合格程度的度量仍存在欠缺。這種風(fēng)險評估的不完備性會導(dǎo)致大量潛在風(fēng)險被忽略，而不能及時作出響應(yīng)決策來預(yù)防食品安全事件發(fā)生。

圖11 污染物導(dǎo)致的不合格記錄占比分布Fig.11 Proportion distribution of unqualified records caused by specific pollutants

因此，進一步依據(jù)1.3.1節(jié)部分闡述的規(guī)則進行計算得到的風(fēng)險值來重新度量水產(chǎn)品污染物中的風(fēng)險分布，如圖12所示（比例代表當(dāng)前污染物導(dǎo)致的不合格產(chǎn)品批次平均風(fēng)險值在總不合格產(chǎn)品批次平均風(fēng)險值中的占比）。可以看出，“氯霉素”導(dǎo)致的不合格產(chǎn)品的不合格程度最大，即其不合格批次的檢測含量遠遠超過了“合格限值”。除此之外，“鎘”和呋喃類污染物的不合格程度也較大。由此，可以推斷，雖然“氯霉素”導(dǎo)致的不合格批次數(shù)量相對較少，但其“不合格程度”較大，存在與“鎘”同樣嚴重的污染問題。

圖12 污染物導(dǎo)致的不合格批次平均風(fēng)險值占比分布Fig.12 Proportion distribution of mean risk values of unqualified samples caused by specific pollutants

圖13 水產(chǎn)品中污染物不合格程度度量對比Fig.13 Comparison of unqualified degree of specific pollutants in aquatic products

為了驗證風(fēng)險值標簽的可靠性，擬用公式（12）來量化污染物的不合格程度。圖13展示了所得數(shù)據(jù)對比情況。可以看出，不考慮規(guī)則，僅從不合格含量和限值角度分析，“氯霉素”的不合格程度仍為最高，且遠遠超出排名第二位的污染物呋喃西林唑酮代謝物，與圖12依據(jù)風(fēng)險值繪制結(jié)果基本吻合。

圖14 水產(chǎn)品中的屬性重要程度排序Fig.14 Ranking of attributes in terms of importance for aquatic products

圖14 展示了水產(chǎn)品風(fēng)險預(yù)警模型訓(xùn)練過程中產(chǎn)生的屬性值重要程度排序分布?？梢钥闯觯廴疚飳︼L(fēng)險的貢獻程度綜合了圖11、12的分布特點。雖然“氯霉素”的不合格程度高，但其不合格批次數(shù)量占比小，在貢獻程度上排名第五位。相比較而言，權(quán)衡了不合格批次數(shù)量和不合格程度兩種度量的“恩諾沙星”和“孔雀石綠”對風(fēng)險的貢獻程度則較為接近。呋喃類污染物的貢獻程度也較符合上述分析結(jié)果。

除此之外，仍可看出離散型屬性值對風(fēng)險的貢獻程度。如“食品類別”中的“海水蟹”相較“海水魚蝦類”產(chǎn)生風(fēng)險的可能性更高，且“海水魚蝦類”對風(fēng)險起到了“抑制”作用；“抽樣省_新疆”、“抽樣省_四川”相較“抽樣省_北京”、“抽樣省_山東”等省份產(chǎn)生風(fēng)險的可能性較大；“存儲狀態(tài)”中“常溫”產(chǎn)生的風(fēng)險遠遠高于未在圖中出現(xiàn)的“避光”、“冷藏”；在“超市”里的“海水魚蝦類”產(chǎn)品存在風(fēng)險的可能性較低等。

3 結(jié) 論

綜上，本研究方法的產(chǎn)出成果主要包含風(fēng)險值預(yù)測、風(fēng)險分析結(jié)論和屬性值重要程度分布3 個部分。風(fēng)險值預(yù)測能對新錄入的數(shù)據(jù)快速計算更加精確的風(fēng)險值，并融合了專家動態(tài)干預(yù)策略，打破了原有方法中固定規(guī)則的局限性。在得出風(fēng)險值的基礎(chǔ)上，統(tǒng)計意義上的風(fēng)險分析結(jié)果給專家提供了干預(yù)參考。最后，屬性值重要程度分布全面闡述了連續(xù)型屬性和離散型屬性對風(fēng)險的貢獻程度，決策制定人員可依據(jù)分布規(guī)律進行屬性值任意組合，制定更加精準的風(fēng)險防控策略，如定點抽檢、高頻抽檢、污染物溯源等。

食品安全問題是關(guān)系民生健康的重要議題。本文首先對食品安全相關(guān)數(shù)據(jù)和以往應(yīng)用的智能化方法進行歸納分析。其次，根據(jù)數(shù)據(jù)特點和已存在方法的不足提出了先驗風(fēng)險概率與模糊層級劃分相結(jié)合的風(fēng)險值計算規(guī)則，并應(yīng)用LightGBM模型結(jié)合專家經(jīng)驗干預(yù)策略進行風(fēng)險值校正和預(yù)測。依托肉制品和水產(chǎn)品數(shù)據(jù)的實驗詳細闡述了模型的使用和參數(shù)調(diào)整方法，并從已有數(shù)據(jù)屬性出發(fā)，驗證了模型的優(yōu)越性和結(jié)論規(guī)律的總體合理性。

但方法仍然存在很多不足。在數(shù)據(jù)應(yīng)用方面，如何將時間序列考慮進來得到更加精細的時序相關(guān)風(fēng)險規(guī)律是一項至關(guān)重要的工作。在模型應(yīng)用方面，如何綜合運用食品各方數(shù)據(jù)，解決食品安全“數(shù)據(jù)孤島”問題，協(xié)同訓(xùn)練得到更加符合實際需求的、更通用的模型，也是后續(xù)工作的重點內(nèi)容。