金學波，張佳帥，郭天洋，王小藝，蘇婷立，賴燕群，孔建磊，白玉廷

（1.北京工商大學人工智能學院，北京 100048；2.北京工商大學 中國輕工業(yè)工業(yè)互聯(lián)網(wǎng)與大數(shù)據(jù)重點實驗室，北京 100048；3.北京工商大學食品與健康學院，北京 100048；4.北京服裝學院，北京 100105；5.中糧糧油工業(yè)（荊州）有限公司，湖北 荊州 434300）

小麥是一種在世界各地廣泛種植的谷類作物，在我國小麥的生長范圍廣泛，是人們喜愛的主食之一，也是牲畜的主要飼料用料之一。據(jù)《國家統(tǒng)計局關于2021年糧食產(chǎn)量數(shù)據(jù)的公告》，2021年全國小麥播種面積23 568.4千 hm，總產(chǎn)量為12 226萬 t。

在小麥的生產(chǎn)中，重金屬污染問題引起了人們的廣泛關注。若食用被重金屬污染的小麥，當污染物生物累積性達到一定濃度后會引起慢性中毒甚至癌癥、生理畸形等病癥。如重金屬鎘在攝入后可在人體的肝、腎等器官組織中蓄積，造成各器官組織的損傷，當達到一定濃度后會破壞骨骼中鈣的正常補給，嚴重者會造成骨軟化癥等惡性疾病。因此，通過預測方法掌握食品加工鏈重金屬污染發(fā)展規(guī)律并及時加以預防控制具有極大的實際應用意義。

目前用于預測建模的方法主要是多元線性回歸方法。例如，牛耀星等采用一級動力學模型結合Arrhenius方程建立了基于金針菇子實體品質指標的貨架期預測模型，并對模型的預測精確度進行了驗證及評價。李婕使用Huang模型、Baranyi模型、修正的Gompertz模型對雙孢蘑菇表面背景菌群和假單胞菌的生長情況進行建模，模型適合描述溫度對雙孢蘑菇樣品表面假單胞菌生長速率的影響。宋波等采用逐步多元回歸分析建立不同鎘濃度等級下水稻籽粒重金屬鎘含量預測模型，模型能夠預測稻米中鎘的累積量，為廣西鎘異常區(qū)內其他水稻產(chǎn)地的安全生產(chǎn)提供了參考。王燁鋒采用多元逐步線性回歸分析方法對土壤-水稻體系中鎘的積累情況建立預測模型，決定系數(shù)（）均達0.7以上。通過實測值反推模型，3 組模型擬合散點值均勻分布在＝線兩側，并通過采集田間樣品對模型進行了驗證。也有一些學者進行了神經(jīng)網(wǎng)絡建模的嘗試，取得了良好的效果。例如，于灝等針對小麥中鎘的含量分別建立了多元回歸模型與淺層神經(jīng)網(wǎng)絡模型，結果表明基于神經(jīng)網(wǎng)絡的小麥鎘含量預測模型的精度總體優(yōu)于多元線性回歸模型。

綜上可知，目前對于原糧產(chǎn)地小麥中鎘含量分析的研究較多，卻鮮有研究涉及小麥加工過程中鎘含量的變化。而在加工過程中，鎘含量的變化與小麥粉的品質直接相關。本研究基于小麥原糧（即原麥）的鎘含量數(shù)據(jù)預測小麥加工成品——小麥粉中的鎘含量，以此為依據(jù)指導加工環(huán)節(jié)中的重金屬危害管控，從而避免小麥粉的重金屬含量超標，保證食品安全。

另一方面，加工過程中的危害物含量具有強非線性、強隨機性的特征，因此線性回歸、淺層神經(jīng)網(wǎng)絡等模型建模精度不高，預測性能有待提升。將深度學習網(wǎng)絡應用于強非線性、強隨機性數(shù)據(jù)預測領域的研究有很多。Haider等建立了一種基于長短期記憶（long and short term memory，LSTM）網(wǎng)絡的小麥產(chǎn)量預測模型。牛哲文等建立了基于門控循環(huán)單元（gated recurrent unit，GRU）的預測模型，使用具有波動性和不確定性的風速、風向等歷史天氣數(shù)據(jù)對風力發(fā)電的功率進行了預測。但是，由于這些深度網(wǎng)絡模型對噪聲非常敏感，并不適用于小麥加工過程中鎘含量預測模型的建立。

為了提升小麥加工過程中鎘含量預測的準確度，本研究基于小麥加工鏈鎘含量數(shù)據(jù)進行深度學習網(wǎng)絡建模，并使用正則化方法優(yōu)化其損失函數(shù)，通過加入噪聲懲罰項防止訓練時模型對噪聲的過度擬合，并使用貝葉斯超參數(shù)優(yōu)化方法對模型超參數(shù)進行優(yōu)化。

1 材料與方法

1.1 材料

小麥為某企業(yè)收購原糧。

1.2 設備

Y7000P便攜式電腦（配備16 GB運行內存和RTX3050Ti顯卡） 聯(lián)想控股股份有限公司。

1.3 方法

1.3.1 小麥加工鏈及關鍵節(jié)點

小麥加工鏈是指從原麥到小麥粉的所有生產(chǎn)環(huán)節(jié)，簡化示意圖如圖1所示。該加工鏈共有清理、潤麥、皮磨等10 個加工環(huán)節(jié)。原糧委托某公司按圖1進行加工，并對各加工環(huán)節(jié)進行抽檢。抽檢方法為：每隔2 h對全加工鏈抽檢1 次，共抽取16 次，每次在各環(huán)節(jié)分別抽取3～4 份（每份1 kg）平行樣品。為便于表示各環(huán)節(jié)的樣本，設定原糧、清理、潤麥、皮磨一（1M芯）、皮磨二（2M芯）、皮磨三（3M芯）、皮磨四（4M芯）、皮磨五（5M芯）、皮磨六（6M芯）、包裝、存放環(huán)節(jié)分別為～。

圖1 小麥加工鏈簡化示意圖Fig. 1 Flow chart of wheat processing chain

1.3.2 鎘含量的測定

取1.3.1節(jié)樣品，委托某公司參照GB/T 20380.4—2006《淀粉及其制品 重金屬含量 第4部分：電熱原子吸收光譜法測定鎘含量》測定鎘含量（以鎘元素計），單位為mg/kg。

1.4 數(shù)據(jù)處理與分析

1.4.1 鎘含量數(shù)據(jù)集構建

取樣時設置3～4 個平行，各環(huán)節(jié)每次取樣的鎘含量結果以平行樣品鎘含量的平均值表示，共獲得16 條全加工鏈的鎘含量數(shù)據(jù)。由于這些數(shù)據(jù)數(shù)量不足以訓練神經(jīng)網(wǎng)絡，利用等差插值的方法對數(shù)據(jù)進行增廣處理。

1.4.2 基于正則化損失函數(shù)的GRU深度網(wǎng)絡模型分析

首先用門控循環(huán)單元建立深度學習模型，其次，設計基于正則化方法優(yōu)化的損失函數(shù)，該損失函數(shù)在模型原損失函數(shù)的基礎上添加訓練數(shù)據(jù)中噪聲的懲罰項，達到在訓練過程中減小對噪聲擬合程度的目的。因此訓練過程中損失函數(shù)會在一個穩(wěn)定值附近波動，該穩(wěn)定值用來量化數(shù)據(jù)中的噪聲。

GRU在時序預測任務中有很好的效果，此前已經(jīng)有大量的研究中都使用GRU作為預測模型并且證明了其性能。GRU模型中只有兩個門，分別是更新門和重置門。更新門控制前一時刻的狀態(tài)信息帶入到當前狀態(tài)的程度，該門的數(shù)值越大就說明前一時刻的信息帶入越多。重置門用于控制前一時刻狀態(tài)的遺忘程度，數(shù)據(jù)越大說明忽略越少。

GRU模型結構如圖2所示，根據(jù)該結構進行前向傳播時每個單元中參數(shù)計算，如式（1）～（4）所示。

式中：x為時刻的輸入；z為更新門衰減系數(shù)；r為重置門衰減系數(shù)；h是-1時刻的狀態(tài)向量；h為時刻狀態(tài)向量；W、U為更新門權重；W、U為重置門權重；W、U為待選a的權重；b、b、b為偏置向量；為Sigmoid激活函數(shù)；*表示逐元素乘法；tanh為tanh激活函數(shù)。

圖2 GRU模型結構圖Fig. 2 Structure diagram of GRU model

但GRU模型在訓練時易對鎘含量數(shù)據(jù)中含有的非線性噪聲過度擬合，降低模型的泛化能力，從而影響模型的預測性能，因此需要采取措施防止模型過度擬合。正則化的作用在于限制模型中的參數(shù)，讓模型的參數(shù)不會過大，從而減少模型過度擬合的可能。在此步驟中，為減少模型因過度學習加工鏈數(shù)據(jù)中的隨機噪聲而導致模型過度擬合，設計噪聲平滑損失函數(shù)，如式（5）所示。

噪聲平滑損失函數(shù)的計算公式中包括兩個部分，第一部分衡量數(shù)據(jù)擬合程度，表示將實際序列與擬合序列之間的殘差平方和最小化的目標；第二部分衡量數(shù)據(jù)平滑程度，表示訓練過程中對序列平滑度的需求。其中，擬合程度用平均絕對誤差表示，平滑程度通過計算P的范數(shù)來實現(xiàn)。是正則化項，將其視為在擬合與平滑兩個目標之間的權衡取舍。

P的計算是實現(xiàn)噪聲平滑損失函數(shù)的關鍵，計算時首先需要定義矩陣，該矩陣對數(shù)據(jù)進行平滑度懲罰。矩陣如式（6）所示。

矩陣的維度由輸入數(shù)據(jù)決定，假設輸入數(shù)據(jù)長度為，則的維度為（，），然后通過與輸入數(shù)據(jù)組成矩陣的乘法計算得到矩陣P，最后通過矩陣范數(shù)的計算得到可以表示數(shù)據(jù)每3 個點間平滑程度的范數(shù)。

在訓練過程中，本實驗所設計的損失函數(shù)可以有效地考慮噪聲對模型的影響從而達到減少過度擬合的作用。當公式（5）中為0時，網(wǎng)絡學習原始數(shù)據(jù)，隨著不斷增大，擬合線趨向平滑。所以，本實驗所設計的損失函數(shù)可以使模型對噪聲具有更強的魯棒性，有效減小噪聲在訓練過程中的影響，從而提高模型的預測準確度。

1.4.3 應用網(wǎng)絡的貝葉斯優(yōu)化方法分析

深度學習模型的超參數(shù)選擇直接決定了模型的性能。本研究基于序列模型的優(yōu)化（sequence model based optimization，SMBO）方法進行模型超參數(shù)優(yōu)化，優(yōu)化的超參數(shù)主要包括GRU中的神經(jīng)元數(shù)、Dropout率、訓練次數(shù)、批處理大小、優(yōu)化器。

貝葉斯優(yōu)化方法在確定模型參數(shù)時，使用代理模型來擬合真實的目標函數(shù)，并根據(jù)擬合結果主動選擇最“潛在”的評估點。所以需要定義一個目標函數(shù)（）和優(yōu)化的超參數(shù)空間。目標函數(shù)（式（7））表示通過貝葉斯優(yōu)化需要達到的最小化的目標，通過設定模型的目標函數(shù)，找到此度量上產(chǎn)生最佳得分的模型超參數(shù)。

貝葉斯優(yōu)化的代理函數(shù)如式（8）所示。

式中：*為貝葉斯優(yōu)化確定的最優(yōu)參數(shù)；為輸入的一組超參數(shù)；為多維超參數(shù)的參數(shù)空間；（）為目標函數(shù)。

貝葉斯優(yōu)化主要由兩個步驟組成：首先通過第+1步估計和更新高斯過程，然后通過最大化代理函數(shù)來指導超參數(shù)的采樣。在高斯過程中，設置目標函數(shù)（）服從高斯分布（式（9））。

式中：（）為目標函數(shù)；（）為目標函數(shù)（）的均值；（，’）為目標函數(shù)（）的協(xié)方差矩陣；表示高斯分布。初始的（，’）如式（10）所示。

式中：為超參數(shù)的數(shù)量；表示協(xié)方差矩陣；（，）、（，w）、（w，）、（w，w）為矩陣內的不同超參數(shù)的組合。

在進行貝葉斯優(yōu)化時，高斯過程的協(xié)方差矩陣會隨著迭代過程改變，假設在第+1步輸入的一組參數(shù)為w，則此時協(xié)方差矩陣如式（11）所示。

式中：’為+1時刻的協(xié)方差矩陣；為初始的協(xié)方差矩陣；w為+1時刻的超參數(shù)值；為超參數(shù)的組合矩陣，為的轉置，的組合公式為＝[（w，），（w，），…，（w，w）]，此時可以得到目標函數(shù)的后驗概率（式（12））。

得到后驗概率后，通過超參數(shù)搜索方法來尋找最優(yōu)超參數(shù)，本實驗使用上置信界（upper confidence bound，UCB）采集函數(shù)（式（13））來完成超參數(shù)搜索。

1.4.4 預測效果評價標準與模型對比驗證

以80%的鎘含量數(shù)據(jù)用作訓練集，對模型的性能進行訓練優(yōu)化；以20%的鎘含量數(shù)據(jù)數(shù)據(jù)作為驗證集，驗證模型的性能。采用均方根誤差（root mean square error，RMSE）、平均絕對誤差（mean absolute error，MAE）、皮爾遜相關系數(shù)評估預測的性能，計算分別如式（14）～（16）所示。三者均反映模型預測值與真實值之間的關系。RMSE、MAE越小，說明模型預測結果與真實情況越接近。越大，說明模型預測結果與真實情況越接近。

采用線性回歸、循環(huán)神經(jīng)網(wǎng)絡（recurrent neural network，RNN）、LSTM和GRU模型與本研究所得模型進行對比。通過比較各模型基于小麥加工鏈鎘含量數(shù)據(jù)的預測效果評價指標（RMSE、MAE、）驗證模型的預測性能，驗證本研究所得模型在實際小麥加工鏈鎘含量預測中的可靠性。

2 結果與分析

2.1 小麥加工鏈中小麥的鎘含量數(shù)據(jù)集構建

小麥加工鏈各環(huán)節(jié)樣品的鎘含量數(shù)據(jù)經(jīng)增廣處理后最終獲得2 057 條全加工鏈的鎘含量數(shù)據(jù)。以其中3 條全加工鏈的鎘含量變化為例，如圖3所示，在清理環(huán)節(jié)、潤麥環(huán)節(jié)中鎘含量略有增加，之后各環(huán)節(jié)鎘含量的變化相對平穩(wěn)。

圖3 各加工環(huán)節(jié)不同樣本中鎘含量Fig. 3 Cadmium contents in different samples measured during processing

2.2 不同模型預測的性能對比

基于2.1節(jié)小麥加工鏈10 個環(huán)節(jié)中2 057 條全加工鏈的鎘含量數(shù)據(jù)進行實驗。將本研究所得模型與線性回歸、RNN、LSTM、GRU模型進行對比，采用RMSE、MAE、評估模型的性能，結果如表1所示。本研究所得模型的RMSE（0.259 5）最小、MAE（0.161 6）最小且最大，比GRU模型的預測性能更好（RMSE降低了46.37%）。

表1 不同預測方法基于鎘數(shù)據(jù)的性能對比Table 1 Performance comparison of different prediction methods based on cadmium data

通過驗證集的驗證與對比分析可知，本研究所得模型在預測任務中表現(xiàn)更為出色，擁有更強的含噪時序數(shù)據(jù)分析能力。由此可見，本研究所得模型可以有效地分析含噪時間序列數(shù)據(jù)，對于小麥加工鏈危害物等含噪數(shù)據(jù)有良好的預測性能。

2.3 模型預測結果

基于2.1節(jié)中20%的原糧鎘含量數(shù)據(jù)，對后續(xù)清理、潤麥、皮磨等10 個環(huán)節(jié)鎘含量進行預測，其中6 個樣本預測結果如圖4所示。其中，清理環(huán)節(jié)和潤麥環(huán)節(jié)樣品的鎘含量略有上升，尤其是清理環(huán)節(jié)，由于鎘元素在小麥表皮的富集，這兩個環(huán)節(jié)的浸潤對鎘元素有一定的析出作用。皮磨一及皮磨二環(huán)節(jié)后鎘含量持續(xù)下降，說明皮磨對危害物鎘有較好的去除作用。預測結果還表明，～環(huán)節(jié)鎘含量基本趨于平穩(wěn)。

GB 2762—2017《食品安全國家標準 食品中污染物限量》規(guī)定小麥原糧及成品糧的鎘含量不得高于0.1 mg/kg。由本研究所得模型可知，當小麥原糧中鎘含量低于GB 2762—2017規(guī)定的0.1 mg/kg時，加工鏈中清理環(huán)節(jié)和潤麥環(huán)節(jié)的鎘含量可能會使其暫時超標，但最終成品小麥粉符合GB 2762—2017要求；而當原麥中鎘含量高于0.1 mg/kg時，成品小麥粉中鎘含量將高于0.1 mg/kg。因此考慮檢測數(shù)據(jù)偏差的情況下，為確保小麥粉的食用安全，原糧中危害物鎘含量高于0.1 mg/kg時不能進行后續(xù)加工流程。

圖4 模型預測小麥加工鏈中樣品鎘含量的變化規(guī)律Fig. 4 Trend of cadmium content in wheat processing chain predicted by the model

3 結 論

本實驗基于小麥加工鏈各個環(huán)節(jié)小麥中的鎘含量。采用深度學習的建模方法對加工鏈各環(huán)節(jié)樣品中的鎘含量進行建模預測。所建模型首先使用正則化方法修改深度學習模型的損失函數(shù)，通過加入噪聲懲罰項來減小噪聲項對模型預測性能的影響。同時使用貝葉斯優(yōu)化算法優(yōu)化模型超參數(shù)。經(jīng)過對比驗證，本研究所建模型對于實際的小麥加工鏈鎘含量數(shù)據(jù)預測準確度較高。

應用該預測模型可以準確地預測出加工鏈中從原麥到小麥粉的各個關鍵環(huán)節(jié)鎘含量的變化，能夠對小麥加工鏈的鎘危害物防控實踐進行指導。預測結果表明，通過加工，原麥中的鎘含量呈整體下降趨勢。但是，在清理和潤麥的環(huán)節(jié)，鎘含量會略有上升。由預測模型可知，如果原麥的鎘含量滿足GB 2762—2017要求，即小麥中的鎘含量不高于0.1 mg/kg，則經(jīng)過加工的成品小麥粉基本滿足GB 2762—2017的要求。